TENDER becomes Tender
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TestAOD / HighLevelQA / NumberOfEventsMCvsDATA.C
1 /////////////////////////////////////////////////////////
2 // NumberOfEventsMCvsDATA.C                            //
3 //                                                     //
4 // plots the ratio (# of MC events)/(# of DATA events) //
5 //    against the run number                           //
6 //                                                     //
7 // Written by John Groh                                //
8 /////////////////////////////////////////////////////////
9
10 // for cout
11 #include <iostream>
12 using namespace std;
13
14 // possible cuts to use -   0    1    2    3
15 Double_t CentCutMin[4] =   {0,  30,  30,  30};
16 Double_t CentCutMax[4] =   {5,  40,  40,  40};
17 Double_t QvecCutMin[4] =   {0,   0,   0, 1.5};
18 Double_t QvecCutMax[4] = {100, 100, 0.4, 100};
19 Double_t EtaMin[4] =    {-0.8,-0.8,-0.8,-0.8};
20 Double_t EtaMax[4] =     {0.8, 0.8, 0.8, 0.8};
21 Double_t Nsigmapid = 3.;
22 UInt_t trkbit = 1024;
23
24 //--------------------
25 // Start of function -
26 //--------------------
27 // icut selects which cut to use (2 and 3 have too low statistics and cause divide-by-zero errors, so don't use)
28 // nSigmaCut is the size of the deviation with which to determine outliers
29 void NumberOfEventsMCvsDATA(Int_t icut = 1, const Float_t nSigmaCut = 3)
30 {  
31   // load libraries (I'm not sure I need all of these...)
32   gSystem->Load("libCore.so");
33   gSystem->Load("libPhysics.so");
34   gSystem->Load("libTree");
35   gSystem->Load("libMatrix");
36   gSystem->Load("libSTEERBase");
37   gSystem->Load("libESD");
38   gSystem->Load("libAOD");
39   gSystem->Load("libANALYSIS");
40   gSystem->Load("libOADB");
41   gSystem->Load("libANALYSISalice");
42   gSystem->Load("libTender");
43   gSystem->Load("libCORRFW");
44   gSystem->Load("libPWGTools");
45   gSystem->Load("libPWGLFspectra");
46   gSystem->Load("libProof.so");
47   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
48
49   // get number of runs used
50   Int_t nRunsTemp = 0;
51   Int_t dummy;
52   ifstream runList("output/AODQAChecks/RunList.txt");
53   while (!runList.eof())
54     {
55       runList >> dummy;
56       nRunsTemp++;
57     }
58   runList.close();
59
60   // I know this is silly, but you need a const for array sizes.
61   // Also, it seems to want to read one more line than there is in the file
62   const Int_t nRuns = nRunsTemp - 1;
63
64   // fill an array with the run numbers and another two with the file names
65   runList.open("output/AODQAChecks/RunList.txt");
66   Int_t runs[nRuns];
67   TString dataFilesDATA[nRuns];
68   TString dataFilesMC[nRuns];
69
70   for (Int_t irun=0; irun<nRuns; irun++)
71     {
72       runList >> runs[irun];
73       dataFilesDATA[irun] = Form("output/AODQAChecks/DATA/AnalysisResults%i.root",runs[irun]);
74       dataFilesMC[irun] = Form("output/AODQAChecks/MC/AnalysisResults%i.root",runs[irun]);
75     }
76   runList.close();
77
78   // choose which TDirectories in the .root files to use
79   TString snameDATA;
80   TString snameMC;
81   snameDATA = Form("OutputAODSpectraTask_Data_Cent%.0fto%.0f_QVec%.1fto%.1f_Eta%.1fto%.1f_%.1fSigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],QvecCutMin[icut],QvecCutMax[icut],EtaMin[icut],EtaMax[icut],Nsigmapid,trkbit);
82   snameMC = Form("OutputAODSpectraTask_MC_Cent%.0fto%.0f_QVec%.1fto%.1f_Eta%.1fto%.1f_%.1fSigmaPID_TrBit%d",CentCutMin[icut],CentCutMax[icut],QvecCutMin[icut],QvecCutMax[icut],EtaMin[icut],EtaMax[icut],Nsigmapid,trkbit);
83
84   TFile * fileDATA;
85   TDirectoryFile * DirDATA;
86   AliSpectraAODEventCuts * ecutsDATA;
87   TFile * fileMC;
88   TDirectoryFile * DirMC;
89   AliSpectraAODEventCuts * ecutsMC;
90
91   // histogram to be filled w/ the ratio (# of MC events)/(# of DATA events)
92   TH1F * hNEventsRatio = new TH1F("hNEventsRatio","",nRuns,0,nRuns);
93
94   // loop over all runs
95   for (Int_t irun=0; irun<nRuns; irun++)
96     {
97       Printf("\n--- Processing run %i ---",runs[irun]);
98       
99       // open the file
100       fileDATA = TFile::Open(dataFilesDATA[irun].Data());
101       if (!fileDATA)
102         {
103           Printf("\n!!! ERROR: Could not open DATA file %s !!!\n",dataFiles[irun].Data());
104           continue;
105         }
106       fileMC = TFile::Open(dataFilesMC[irun].Data());
107       if (!fileMC)
108         {
109           Printf("\n!!! ERROR: Could not open MC file %s !!!\n",dataFiles[irun].Data());
110           continue;
111         }
112       
113       // choose the right directory and event cut objects
114       DirDATA = (TDirectoryFile*)fileDATA->Get(Form("%s",snameDATA.Data()));
115       if (!DirDATA)
116         {
117           Printf("\n!!! ERROR: DirDATA is a null pointer. Skipping to next file !!!\n");
118           continue;
119         }
120       DirMC = (TDirectoryFile*)fileMC->Get(Form("%s",snameMC.Data()));
121       if (!DirMC)
122         {
123           Printf("\n!!! ERROR: DirMC is a null pointer. Skipping to next file !!!\n");
124           continue;
125         }
126       ecutsDATA = (AliSpectraAODEventCuts*)DirDATA->Get("Event Cuts");
127       if (!ecutsDATA)
128         {
129           Printf("\n!!! ERROR: ecutsDATA is a null pointer. Skipping to next file !!!\n");
130           continue;
131         }      
132       ecutsMC = (AliSpectraAODEventCuts*)DirMC->Get("Event Cuts");
133       if (!ecutsMC)
134         {
135           Printf("\n!!! ERROR: ecutsMC is a null pointer. Skipping to next file !!!\n");
136           continue;
137         }      
138
139       // calculate the ratio and fill the histogram
140       Int_t nEventsDATA = ecutsDATA->NumberOfEvents();
141       Int_t nEventsMC = ecutsMC->NumberOfEvents();
142       Float_t ratio = (Float_t)nEventsMC / (Float_t)nEventsDATA;
143       hNEventsRatio->SetBinContent(irun+1,ratio);
144       Printf("# of MC events: %i\n# of DATA events: %i\nRatio: %.2f",nEventsMC,nEventsDATA,ratio);
145
146     } // end loop over runs
147
148   // draw it
149   TCanvas * cNEventsRatio = new TCanvas("cNEventsRatio","cNEventsRatio");
150   for (Int_t irun=0; irun<nRuns; irun++)
151     hNEventsRatio->GetXaxis()->SetBinLabel(irun+1,Form("%i",runs[irun]));
152   hNEventsRatio->GetYaxis()->SetTitle("(# of MC events)/(# of DATA events)");
153   hNEventsRatio->GetYaxis()->SetTitleOffset(1.2);
154   hNEventsRatio->SetMarkerStyle(21);
155   hNEventsRatio->SetMarkerColor(kBlue);
156   hNEventsRatio->SetStats(kFALSE);
157   hNEventsRatio->DrawCopy("P");
158 }
159
160
161
162
163