Analyzing CORSIKA output

Hi. I am unable to edit CorsikaPlotter to analyze Corsika Output using the instruction given.

Edit CorsikPlotter.cc, with editor of your choice, e.g.
gedit CorsikaPlotter.cc
 Add ROOT 1D histogram to study shower zenith angle distribution

20 bins from 0 to 90 [degree], name: zenith
 Fill shower zenith into histogram (take care of units)

Note: check COAST Documentation link on ISAPP 2018 indico page
 Create TCanvas, draw histogram, save pdf format
 Type make in directory to compile.

I am just a beginner kindly help.

#include <crsRead/MCorsikaReader.h>

#include <crs/TSubBlock.h>
#include <crs/MRunHeader.h>
#include <crs/MEventHeader.h>
#include <crs/MEventEnd.h>
#include <crs/MParticleBlock.h>
#include <crs/MLongitudinalBlock.h>
#include <crs/MParticle.h>

#include <TCanvas.h>
#include <TH1D.h>

#include
using namespace std;

int
main (int argc, char **argv)
{
if (argc<2) {
cout << " Name of the file " << endl;
return 1;
}

// CREATE HISTOGRAMS HERE …

// …

string fname(argv[1]);
crsRead::MCorsikaReader cr(fname, 3);

int ShowerCounter = 0;
crs::MRunHeader Run;
while (cr.GetRun (Run)) {

crs::MEventHeader Shower;
while (cr.GetShower(Shower)) {

  
  // FILL SHOWER-LEVEL HISTOGRAMS HERE ....
  
  // ........................................
  

  ++ShowerCounter;           
  crs::TSubBlock Data;
  while (cr.GetData (Data)) {
    
    switch (Data.GetBlockType ()) {
      
        case crs::TSubBlock::ePARTDATA:
        {
          const crs::MParticleBlock& ParticleData = Data;
          crs::MParticleBlock::ParticleListConstIterator iEntry;
          for (iEntry = ParticleData.FirstParticle();
               iEntry != ParticleData.LastParticle();
               ++iEntry) {
            
            if (iEntry->IsParticle()) {
              
              crs::MParticle iPart(*iEntry);
              
              const int id    = iPart.GetParticleID();
              const double w  = iPart.GetWeight();
              const double e  = iPart.GetKinEnergy();
              const double x  = iPart.GetX();
              const double y  = iPart.GetY();

	  
	  // FILL PARTICLE-LEVEL HISTOGRAMS HERE ....

	  // ........................................
	  
	  
            }
            
          } // end particle loop
          
          break;
        }
        
        case crs::TSubBlock::eLONG:
          break;
          
        default:
          break;
    } // end data block
    
  } // loop data
  
  crs::MEventEnd ShowerSummary;
  cr.GetShowerSummary(ShowerSummary);
  const double Xmax = ShowerSummary.GetXmax();

} // loop shower

} // loop runs (usually just 1)

// PERFORM FINAL COMPUTATION ON HISTOGRAMS …
// AND PRINT/SAVE

// TCanvas* canvas = new TCanvas(“printout”);

// canvas->SaveAs(“canvas.pdf”);

// …

cout << " Read " << ShowerCounter << " showers from file " << endl;

return 0;
}

My Output is in this format.

/home/sam/Chege1/DAT000002

1 Like

You need to contact the “Corsika” support team and/or your colleagues and/or your supervisor (hopefully, they know “Corsika”).

Thank you. will consult them

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.