FANDOM


 
/***********************************************************
  *
  *    Read simulated electrons and positrons from a
  *    collection of .slcio files and plot some distributions
  *
  ***********************************************************/
 
 #include "stdlib.h"
 
 // LCIO includes
 #include "lcio.h"
 #include <IOIMPL/LCFactory.h>
 #include <IMPL/LCCollectionVec.h>
 #include <EVENT/MCParticle.h>
 #include <IMPL/MCParticleImpl.h>
 #include <IMPL/ReconstructedParticleImpl.h>
 #include <UTIL/LCTOOLS.h>
 #include <Exceptions.h>
 
 // ROOT includes
 #include "TROOT.h"
 #include "TLorentzVector.h"
 
 #include <sstream>
 #include <iomanip>
 
 using namespace std;
 
 void electronenergy(UInt_t nFirstJob, UInt_t nLastJob, const char * fileName)
 {
     // Obtain an LCReader
     IO::LCReader* lcReader = IOIMPL::LCFactory::getInstance()->createLCReader() ;
 
     TString fName = fileName;
     //if(!(fName.EndsWith("_"))) fName += "";
 
     // Declare an LCEvent pointer
     EVENT::LCEvent* evt;
 
     // Read each file (first job -to- last job)
     stringstream fNameStream, cfNameStream;
     for(UInt_t iJob = nFirstJob; iJob <= nLastJob; iJob++)
     {
         // Construct the actual file name ("BaseName.JobNumber.slcio")
         fNameStream.str("");
         fNameStream << fName.Data() << "." << iJob << ".slcio";
 
         try
         {
             const char* fileNameString = fNameStream.str().c_str();
             lcReader->open(fileNameString);
         }
         catch(lcio::IOException &ex)
         {
             // Could not open that particular file. Just ignore and try the next one...
             continue;
         }
 
         // Ok, file opened.
 
         // Repeat for each event, as long as there are events left in the reader
         while( (evt = lcReader->readNextEvent()) != 0 )
         {
             // Get the collection of Monte Carlo particles from the event
             EVENT::LCCollection* pMCParticles = evt->getCollection("MCParticle");
 
             // This type will help us out by doing some conversions for us
             TLorentzVector electron;
             electron *= 0.0;   // zero-out memory in the variable
 
             // Go over each particle in the collection, and look for particles of specific type
             for (Int_t i = 0; i < pMCParticles->getNumberOfElements(); i++)
             {
                 // Get the particle at index i
                 EVENT::MCParticle* pMCParticle = (EVENT::MCParticle*)pMCParticles->getElementAt(i);  // need to cast, as getElementAt returns an LCObject*
 
                 // If it's an undecayed particle, stable in the generator, just skip
                 if (pMCParticle->getGeneratorStatus() != 1)
                     continue;
 
                 // Check PDG code to see what kind of particle it is;
                 // only count if it's an electron (11) or a positron (-11)
                 if(pMCParticle->getPDG() == 11 || pMCParticle->getPDG() == -11)
                 {
                     const Double_t* p1 = pMCParticle->getMomentum();   // returns an array of doubles, length 3
                     Double_t e1 = pMCParticle->getEnergy();
 
                     // Set the data to the electron helper object
                     electron.SetPxPyPzE(p1[0], p1[1], p1[2], e1);
 
                     // Extract interesting quantities
                     Double_t energy = electron.Energy();
                     Double_t angle = electron.Theta() * (180.0 / M_PI);
 
                     // Process the data here...
                     // (... omitted ...)
 
                 }
             }
         }
 
         // Done reading this file, open the next one (on the next loop iteration).
         // Make sure to close the LCReader.
         lcReader->close();
     }
 
     // Do some more stuff if required...
     // (... omitted ...)
 }
 
 // Main application entry point - the program starts here
 Int_t main(int argc, char* argv[])
 {
     // Start from the first argument to the program
     Int_t iarg = 1;
 
     // Set default values for jobs, to be used if parameters not supplied
     UInt_t nFirstJob = 1;
     UInt_t nLastJob = 10;
 
     // Read the first and the last job from the two successive parameters (e.g. Prog.exe 100 200)
     if(argc>iarg) nFirstJob = atoi(argv[iarg]); iarg++;
     if(argc>iarg) nLastJob   = atoi(argv[iarg]); iarg++;
 
     // Set default base file name to use if the user doesn't supply one
     TString fName = "mokka_pairs";
 
     // Read the next parameter into the fName variable (e.g. Prog.exe 100 200 mokka_pairs_alternate)
     if(argc>iarg) fName = argv[iarg]; iarg++;
 
     // Call the processing function defined above
     electronenergy(nFirstJob, nLastJob, fName.Data());
 
     return 0;
 }

Ad blocker interference detected!


Wikia is a free-to-use site that makes money from advertising. We have a modified experience for viewers using ad blockers

Wikia is not accessible if you’ve made further modifications. Remove the custom ad blocker rule(s) and the page will load as expected.