eic-smear/scripts/read.cxx

0001 // read.cxx
0002 //
0003 // Created by TB on 6/13/11.
0004 // Copyright 2011 BNL. All rights reserved.
0005 //
0006 // Example of how to read a file produced by BuildTree for a simple analysis.
0007 // To run, in ROOT do:
0008 // root [0] .L /path/to/read.cxx
0009 // root [1] read("myInputFile.root", 10000 )
0010
0011 void read(TString inFileNames, int nEvents ) {
0012
0013    // If the analysis solely uses TTree::Draw statements,
0014    // you don't need to load
0015    // the shared library. You will receive warnings such as
0016    // Warning in <TClass::TClass>: no dictionary for class Particle
0017    // is available
0018    // but these can be ignored. However if you wish to work with the event
0019    // objects themselves, the shared library must be loaded:
0020    // Load the shared library, if not done automaticlly:
0021    //   gSystem->Load("/path/to/libeicsmear.so" );
0022
0023    // The TTrees are named EICTree.
0024    // Create a TChain for trees with this name.
0025    TChain tree("EICTree");
0026
0027    // Add the file(s) we want to analyse to the chain.
0028    // We could add multiple files if we wanted.
0029    tree.Add(inFileNames); // Wild cards are allowed e.g. tree.Add("*.root" );
0030 // tree.Add(/path/to/otherFileNames ); // etc...
0031
0032    // Create an object to store the current event from the tree.
0033    // This is how we access the values in the tree.
0034    // If you want to use generator-specific values, then
0035    // the event type should match the type in the TTree. Valid types are
0036    // EventPythia, EventPepsi, EventRapgap, EventDjangoh, EventMilou.
0037    // If you only need shared quantities like x, Q2 and the particle list
0038    // you can use EventBase and the macro will be general for any Monte Carlo.
0039    erhic::EventPythia* event(NULL);// = new EventPythia;
0040 // EventBase* event(NULL);
0041
0042    // Now associate the contents of the branch with the buffer.
0043    // The events are stored in a branch named event:
0044    tree.SetBranchAddress("event", &event ); // Note &event, not event.
0045
0046    // Now we can do some analysis...
0047
0048    // We record the largest particle pT we find here:
0049    double highestPt(-1. );
0050
0051    // Histograms for our analysis.
0052    TH2D Q2VsX("Q2VsX",
0053               "Q^{2} vs. Bjorken x;log_{10}(x);log_{10}(Q^{2})",
0054                100, -5., 0., 100, -2., 3. );
0055    TH1D ptHist("ptHist",
0056                "pT of charged pions",
0057                50, 0.0, 2.0 );
0058    TH1D deltaPhi("deltaPhi",
0059                  "Delta-phi of hadrons",
0060                  40, 0.0, 3.1415 );
0061
0062    // Loop over events:
0063    for(int i(0); i < nEvents; ++i ) {
0064
0065       // Read the next entry from the tree.
0066       tree.GetEntry(i);
0067
0068       // Fill the Q2 vs. x histogram:
0069       Q2VsX.Fill(TMath::Log10(event->GetX()),
0070                  TMath::Log10(event->GetQ2()));
0071
0072       // The event contains a vector (array) of particles.
0073       int nParticles = event->GetNTracks();
0074
0075       // We now know the number of particles in the event, so loop over
0076       // the particles:
0077       for(int j(0); j < nParticles; ++j ) {
0078          const Particle* particle = event->GetTrack(j);
0079          // Let's just select charged pions for this example:
0080          int pdg = particle->GetPdgCode();
0081          if(abs(pdg) != 211 ) continue;
0082
0083          ptHist.Fill(particle->GetPt());
0084
0085          // Update the highest pT:
0086          if(particle->GetPt() > highestPt ) {
0087             highestPt = particle->GetPt();
0088          } // if
0089       } // for
0090    } // for
0091
0092    std::cout << "The highest pT was " << highestPt << " GeV/c" << std::endl;
0093
0094    TCanvas canvas;
0095    Q2VsX.Draw("colz" );
0096    canvas.Print("Q2VsX.png" );
0097    ptHist.Draw();
0098    canvas.Print("pt.png" );
0099 }