EICMCAnalysis/analysis_eicsmear/eic_unsmeared_eta_p.C

0001 // A quick look at the momentum vs pseudorapidity of the unsmeared particles
0002 // Currently loads both smeared and unsmeared input files for later use
0003 // Only unsmeared input file used right now
0004
0005 // --------------------------
0006
0007 /* STL includes */
0008 #include <iostream>
0009
0010 /* ROOT includes */
0011 #include <TROOT.h>
0012 //#include <gSystem.h>
0013 #include <TFile.h>
0014 #include <TTree.h>
0015 #include <TH2F.h>
0016 #include <THnSparse.h>
0017
0018 /* eicsmear includes */
0019 //#include <eicsmear/erhic/EventPythia.h>
0020 //#include <eicsmear/erhic/ParticleMC.h>
0021
0022 using namespace std;
0023
0024 int
0025 eic_unsmeared_eta_p( TString filename_output,
0026               TString filename_mc,
0027               TString filename_mc_smeared = "",
0028               bool debug = false )
0029 {
0030   /* Uncomment this line when running without compilation */
0031   gSystem->Load("libeicsmear");
0032
0033   /* Open input files. */
0034   TFile *file_mc = new TFile(filename_mc, "OPEN");
0035   TFile *file_mc_smeared = new TFile(filename_mc_smeared, "OPEN");
0036
0037   /* Get trees from files. */
0038   TTree *tree = (TTree*)file_mc->Get("EICTree");
0039   TTree *tree_smeared = (TTree*)file_mc_smeared->Get("Smeared");
0040
0041   /* Output file. */
0042   TFile *file_out = new TFile(filename_output, "RECREATE");
0043
0044   /* Add friend to match branches in trees. */
0045   tree->AddFriend(tree_smeared);
0046
0047   erhic::EventPythia *event  = NULL;
0048   Smear::Event       *eventS = NULL;
0049
0050   tree->SetBranchAddress("event", &event);
0051   tree->SetBranchAddress("eventS", &eventS);
0052
0053   /* Create histogram for eta versus momentum. */
0054   TH2F* h_eta = new TH2F("h_eta","Momentum vs. Pseudorapidity",80,-4,4,45,-5,40 );
0055   TH2F* h_eta_accept = (TH2F*)h_eta->Clone("h_eta_accept");
0056   h_eta->GetXaxis()->SetTitle("#eta");
0057   h_eta->GetYaxis()->SetTitle("Momentum (GeV)");
0058
0059   /* Create canvas */
0060   TCanvas *c1 = new TCanvas;
0061   h_eta->Draw("AXIS");
0062
0063   /* Loop over all events in tree. */
0064   unsigned max_event = tree->GetEntries();
0065
0066   for ( unsigned ievent = 0; ievent < max_event; ievent++ )
0067     {
0068       if ( ievent%1000 == 0 )
0069         cout << "Processing event " << ievent << endl;
0070
0071       /* load event */
0072       tree->GetEntry(ievent);
0073
0074       /* Cut on EVENT kinematics */
0075       float y = event->GetTrueY();
0076       if ( y > 0.99 || y < 0.01 )
0077         continue;
0078
0079       float eta = event->ScatteredLepton()->GetEta();
0080       float p = event->ScatteredLepton()->GetP();
0081
0082       h_eta->Fill(eta,p);
0083
0084     } // end loop over events
0085
0086   /* Write histogram. */
0087
0088   h_eta->Write();
0089   h_eta_accept->Write();
0090
0091   /*draw histogram */
0092
0093   h_eta->Draw("COLZ");
0094   h_eta->Draw("sameaxis");
0095   return c1;
0096
0097   /* Close output file. */
0098   file_out->Close();
0099
0100   return 0;
0101 }