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");
  33   gSystem->Load("libPhysics");
  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");
  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