EICMCAnalysis/analysis_eicsmear/eic_sphenix_test_smearing.C

0001 /* STL includes */
0002 #include <iostream>
0003
0004 /* ROOT includes */
0005 #include <TROOT.h>
0006 //#include <gSystem.h>
0007 #include <TFile.h>
0008 #include <TTree.h>
0009 #include <TH2F.h>
0010 #include <THnSparse.h>
0011
0012 /* eicsmear includes */
0013 //#include <eicsmear/erhic/EventPythia.h>
0014 //#include <eicsmear/erhic/ParticleMC.h>
0015
0016 using namespace std;
0017
0018 int
0019 eic_sphenix_test_smearing( TString filename_output,
0020               TString filename_mc,
0021               TString filename_mc_smeared = "",
0022               bool debug = false )
0023 {
0024   /* Uncomment this line when running without compilation */
0025   gSystem->Load("libeicsmear");
0026
0027   /* Open input files. */
0028   TFile *file_mc = new TFile(filename_mc, "OPEN");
0029   TFile *file_mc_smeared = new TFile(filename_mc_smeared, "OPEN");
0030
0031   /* Get trees from files. */
0032   TTree *tree = (TTree*)file_mc->Get("EICTree");
0033   TTree *tree_smeared = (TTree*)file_mc_smeared->Get("Smeared");
0034
0035   /* Output file. */
0036   //TFile *file_out = new TFile(filename_output, "RECREATE");
0037
0038   /* Add friend to match branches in trees. */
0039   tree->AddFriend(tree_smeared);
0040
0041   erhic::EventPythia *event  = NULL;
0042   Smear::Event       *eventS = NULL;
0043
0044   tree->SetBranchAddress("event", &event);
0045   tree->SetBranchAddress("eventS", &eventS);
0046
0047   /* Create histogram */
0048   TH2F* h_e_smeared_vs_eta = new TH2F("h_e_smeared_vs_eta","Energy Smeared vs True Pseudorapidity",100,-5,5,70,0,35);
0049   h_e_smeared_vs_eta->GetXaxis()->SetTitle("#eta_{true}");
0050   h_e_smeared_vs_eta->GetYaxis()->SetTitle("E_{smeared} (GeV)");
0051
0052   TH2F* h_e_smeared_vs_true = new TH2F("h_e_smeared_vs_true","Energy Smeared vs True",60,0,30,70,0,35);
0053   h_e_smeared_vs_true->GetXaxis()->SetTitle("E_{true} (GeV)");
0054   h_e_smeared_vs_true->GetYaxis()->SetTitle("E_{smeared} (GeV)");
0055
0056   TH1F* h_eta = new TH1F("h_eta", ";#eta;dN/d#eta", 100, -5, 5);
0057   TH1F* h_eta_accept = (TH1F*)h_eta->Clone("h_eta_accept");
0058   h_eta_accept->SetLineColor(kGreen+4);
0059
0060   TH1F* h_e_eref_true = new TH1F("h_e_eref_true","True reference energy",300,0,30);
0061   h_e_eref_true->GetXaxis()->SetTitle("E_{true} (GeV)");
0062   h_e_eref_true->GetYaxis()->SetTitle("# entries");
0063   h_e_eref_true->SetLineColor(kRed);
0064
0065   TH1F* h_e_eref_smeared = new TH1F("h_e_eref_smeared","Smeared reference energy",300,0,30);
0066   h_e_eref_smeared->GetXaxis()->SetTitle("E_{smeared} (GeV)");
0067   h_e_eref_smeared->GetYaxis()->SetTitle("# entries");
0068   h_e_eref_smeared->SetLineColor(kBlue);
0069
0070   /* Loop over all events in tree. */
0071   unsigned max_event = tree->GetEntries();
0072
0073   for ( unsigned ievent = 0; ievent < max_event; ievent++ )
0074     {
0075       if ( ievent%1000 == 0 )
0076         cout << "Processing event " << ievent << endl;
0077
0078       /* load event */
0079       tree->GetEntry(ievent);
0080
0081       /* Cut on EVENT kinematics */
0082       float y = event->GetTrueY();
0083       if ( y > 0.99 || y < 0.01 )
0084        continue;
0085
0086       float energy = event->ScatteredLepton()->GetE();
0087       float energy_smeared = eventS->ScatteredLepton()->GetE();
0088
0089       float eta = event->ScatteredLepton()->GetEta();
0090
0091       /* Fill histograms */
0092       h_e_smeared_vs_eta->Fill(eta,energy_smeared);
0093       h_e_smeared_vs_true->Fill(energy,energy_smeared);
0094
0095       h_eta->Fill(eta);
0096       if ( energy_smeared > 0 )
0097     h_eta_accept->Fill( eta );
0098
0099       /* Fill histograms if truth energy within range around reference energy */
0100       float eref = 19.05;
0101       float erange = 0.1;
0102       if ( energy > (eref-erange/2.) && energy < (eref+erange/2.) )
0103     {
0104       h_e_eref_true->Fill(energy);
0105       h_e_eref_smeared->Fill(energy_smeared);
0106     }
0107
0108     } // end loop over events
0109
0110   float underflow = tree->GetEntry(0);
0111   float overflow = tree->GetEntry(max_event + 1);
0112   cout << "underflow: " << underflow << endl;
0113   cout << "overflow: " << overflow << endl;
0114
0115   /* Write histograms. */
0116   //h_e_smeared_vs_true->Write();
0117   //h_e_smeared_vs_eta->Write();
0118
0119   /*draw histograms */
0120   TCanvas *c1 = new TCanvas();
0121   h_e_smeared_vs_true->DrawClone("COLZ");
0122   gPad->RedrawAxis();
0123
0124   TCanvas *c2 = new TCanvas();
0125   h_e_smeared_vs_eta->DrawClone("COLZ");
0126   gPad->RedrawAxis();
0127
0128   TCanvas *c3 = new TCanvas();
0129   h_e_eref_smeared->Draw();
0130   h_e_eref_smeared->Fit("gaus");
0131   h_e_eref_true->Draw("sames");
0132   gPad->RedrawAxis();
0133
0134   TCanvas *c4 = new TCanvas();
0135   h_eta_accept->Divide(h_eta);
0136   h_eta_accept->Draw();
0137   gPad->RedrawAxis();
0138
0139   /* Close output file. */
0140   //file_out->Close();
0141
0142   return 0;
0143 }