PWGLF/SPECTRA/PiKaPr/TestAOD/HighLevelQA/AODQAChecks.C

   1 /////////////////////////////////////////////////////////////////////////////////////////
   2 // AODQAChecks.C
   3 //
   4 // Performs some run-by-run high-level QA checks on AOD data
   5 //    for runs marked with global quality "good" in the RCT
   6 //    for the LHC10h period
   7 //
   8 // Note that the following macros must be in the current
   9 //    directory: CheckEfficiencies.C, CheckIntegratedRawYield.C,
  10 //    CheckMatchingEfficiency.C, CheckNSigmaStability.C,
  11 //    FindOutliers.C, and PlotAndSave.C
  12 //
  13 // Written by John Groh, summer student (Summer 2012)
  14 // Advisors: Michele Floris and Alexander Kalweit
  15 /////////////////////////////////////////////////////////////////////////////////////////
  16
  17 // for cout
  18 #include <iostream>
  19
  20 // for processing .txt file with run numbers
  21 #include <fstream>
  22 using namespace std;
  23
  24 // global constants
  25 const Int_t nPart = 3;
  26 TString Particle[nPart] = {"Pion", "Kaon", "Proton"};
  27 const Int_t nCharge = 2;
  28 TString Charge[nCharge] = {"Pos", "Neg"};
  29 TString Sign[nCharge] = {"Plus", "Minus"};
  30 Int_t Color[nPart] = {1,2,4};
  31 Int_t Marker[nPart*nCharge] = {20,21,22,24,25,26};
  32 TString Names[nPart*nCharge] = {"#pi^{+}",
  33                                 "K^{+}",
  34                                 "p",
  35                                 "#pi^{-}",
  36                                 "K^{-}",
  37                                 "#bar{p}"};
  38
  39 // possible cuts to use -   0    1    2    3
  40 Double_t CentCutMin[4] =   {0,  30,  30,  30};
  41 Double_t CentCutMax[4] =   {5,  40,  40,  40};
  42 Double_t QvecCutMin[4] =   {0,   0,   0, 1.5};
  43 Double_t QvecCutMax[4] = {100, 100, 0.4, 100};
  44 Double_t EtaMin[4] =    {-0.8,-0.8,-0.8,-0.8};
  45 Double_t EtaMax[4] =     {0.8, 0.8, 0.8, 0.8};
  46 Double_t Nsigmapid = 3.;
  47 UInt_t trkbit = 1024;
  48
  49 // fixed Pt values chosen for efficiency and matching efficiency plots (GeV)
  50 const Float_t FixedPtEff = 0.4;
  51 const Float_t FixedPtMatchEff = 0.9;
  52
  53 // inclusions for functions called
  54 #include "CheckIntegratedRawYield.C"
  55 #include "CheckNSigmaStability.C"
  56 #include "CheckEfficiencies.C"
  57 #include "CheckMatchingEfficiency.C"
  58 #include "PlotAndSave.C"
  59 #include "FindOutliers.C"
  60
  61 //--------------------
  62 // Start of function -
  63 //--------------------
  64 // useMC selects whether to use Monte Carlo or Data - must be run seperately for each
  65 // icut selects which cut to use (2 and 3 have low statistics and cause divide-by-zero errors, so don't use)
  66 // nSigmaCut is the size of the deviation with which to determine outliers
  67 void AODQAChecks(Bool_t useMC = 1, Int_t icut = 1, const Float_t nSigmaCut = 3)
  68 {
  69   Printf("\n\n--- Running AODQAChecks() with useMC = %i ---\n\n",useMC);
  70
  71   gSystem->Load("libCore.so");
  72   gSystem->Load("libPhysics.so");
  73   gSystem->Load("libTree");
  74   gSystem->Load("libMatrix");
  75   gSystem->Load("libSTEERBase");
  76   gSystem->Load("libESD");
  77   gSystem->Load("libAOD");
  78   gSystem->Load("libANALYSIS");
  79   gSystem->Load("libOADB");
  80   gSystem->Load("libANALYSISalice");
  81   gSystem->Load("libTENDER");
  82   gSystem->Load("libCORRFW");
  83   gSystem->Load("libPWGTools");
  84   gSystem->Load("libPWGLFspectra");
  85   gSystem->Load("libProof.so");
  86   gSystem->AddIncludePath("-I$ALICE_ROOT/include");
  87
  88   TString fold = "AODQAChecks";
  89
  90   // get number of runs used
  91   Int_t nRunsTemp = 0;
  92   Int_t dummy;
  93   ifstream runList("output/AODQAChecks/RunList.txt");
  94   while (!runList.eof())
  95     {
  96       runList >> dummy;
  97       nRunsTemp++;
  98     }
  99   runList.close();
 100
 101   // I know this is silly, but you need a const for array sizes.
 102   // Also, it seems to want to read one more line than there is in the file
 103   const Int_t nRuns = nRunsTemp - 1;
 104
 105   // fill an array with the run numbers and another with the file names
 106   runList.open("output/AODQAChecks/RunList.txt");
 107   Int_t runs[nRuns];
 108   TString dataFiles[nRuns];
 109
 110   for (Int_t irun=0; irun<nRuns; irun++)
 111     {
 112       runList >> runs[irun];
 113       if (useMC) dataFiles[irun] = Form("output/%s/MC/AnalysisResults%i.root",fold.Data(),runs[irun]);
 114       else dataFiles[irun] = Form("output/%s/DATA/AnalysisResults%i.root",fold.Data(),runs[irun]);
 115     }
 116   runList.close();
 117
 118   // choose which TDirectory in the .root file to use
 119   TString sname;
 120   if (useMC)
 121     sname = 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);
 122   else
 123     sname = 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);
 124
 125   // create a .root file for output
 126   TFile * fout = new TFile(Form("results/%s/Res_%s_AODQAChecks.root",fold.Data(),sname.Data()),"RECREATE");
 127   TFile * file;
 128   TDirectoryFile * Dir;
 129   AliSpectraAODTrackCuts * tcuts;
 130   AliSpectraAODEventCuts * ecuts;
 131   AliSpectraAODHistoManager * hman;
 132
 133   // histograms for stability of nsigma distribution plots
 134   TH1F * TPCnsigMeanTrendPion = new TH1F("TPCnsigMeanTrendPion","",nRuns,0,nRuns);
 135   TH1F * TPCnsigMeanTrendKaon = new TH1F("TPCnsigMeanTrendKaon","",nRuns,0,nRuns);
 136   TH1F * TPCnsigMeanTrendProton = new TH1F("TPCnsigMeanTrendProton","",nRuns,0,nRuns);
 137   TH1F * TPCnsigSigmaTrendPion = new TH1F("TPCnsigSigmaTrendPion","",nRuns,0,nRuns);
 138   TH1F * TPCnsigSigmaTrendKaon = new TH1F("TPCnsigSigmaTrendKaon","",nRuns,0,nRuns);
 139   TH1F * TPCnsigSigmaTrendProton = new TH1F("TPCnsigSigmaTrendProton","",nRuns,0,nRuns);
 140   TH1F * TOFnsigMeanTrendPion = new TH1F("TOFnsigMeanTrendPion","",nRuns,0,nRuns);
 141   TH1F * TOFnsigMeanTrendKaon = new TH1F("TOFnsigMeanTrendKaon","",nRuns,0,nRuns);
 142   TH1F * TOFnsigMeanTrendProton = new TH1F("TOFnsigMeanTrendProton","",nRuns,0,nRuns);
 143   TH1F * TOFnsigSigmaTrendPion = new TH1F("TOFnsigSigmaTrendPion","",nRuns,0,nRuns);
 144   TH1F * TOFnsigSigmaTrendKaon = new TH1F("TOFnsigSigmaTrendKaon","",nRuns,0,nRuns);
 145   TH1F * TOFnsigSigmaTrendProton = new TH1F("TOFnsigSigmaTrendProton","",nRuns,0,nRuns);
 146
 147   // histograms for integrated raw yield stability plots
 148   TH1F * IntegRawYieldAll = new TH1F("IntegRawYieldAll","",nRuns,0,nRuns);
 149   TH1F * IntegRawYieldPiPlus = new TH1F("IntegRawYieldPiPlus","",nRuns,0,nRuns);
 150   TH1F * IntegRawYieldPiMinus = new TH1F("IntegRawYieldPiMinus","",nRuns,0,nRuns);
 151   TH1F * IntegRawYieldKPlus = new TH1F("IntegRawYieldKPlus","",nRuns,0,nRuns);
 152   TH1F * IntegRawYieldKMinus = new TH1F("IntegRawYieldKMinus","",nRuns,0,nRuns);
 153   TH1F * IntegRawYieldProton = new TH1F("IntegRawYieldProton","",nRuns,0,nRuns);
 154   TH1F * IntegRawYieldAntiproton = new TH1F("IntegRawYieldAntiproton","",nRuns,0,nRuns);
 155
 156   // objects for efficiency plots
 157   if (useMC)
 158     {
 159       TCanvas * cEfficienciesAllRuns = new TCanvas("cEfficienciesAllRuns","cEfficienciesAllRuns",50,25,700,500);
 160       cEfficienciesAllRuns->Divide(3,2);
 161       TH1F * EfficiencyPiPlus = new TH1F("EfficiencyPiPlus","",nRuns,0,nRuns);
 162       TH1F * EfficiencyKPlus = new TH1F("EfficiencyKPlus","",nRuns,0,nRuns);
 163       TH1F * EfficiencyProton = new TH1F("EfficiencyProton","",nRuns,0,nRuns);
 164       TH1F * EfficiencyPiMinus = new TH1F("EfficiencyPiMinus","",nRuns,0,nRuns);
 165       TH1F * EfficiencyKMinus = new TH1F("EfficiencyKMinus","",nRuns,0,nRuns);
 166       TH1F * EfficiencyAntiproton = new TH1F("EfficiencyAntiproton","",nRuns,0,nRuns);
 167     }
 168   else
 169     {
 170       TCanvas * cEfficienciesAllRuns = 0x0;
 171       TH1F * EfficiencyPiPlus = 0x0;
 172       TH1F * EfficiencyKPlus = 0x0;
 173       TH1F * EfficiencyProton = 0x0;
 174       TH1F * EfficiencyPiMinus = 0x0;
 175       TH1F * EfficiencyKMinus = 0x0;
 176       TH1F * EfficiencyAntiproton = 0x0;
 177     }
 178
 179
 180   // objects for matching efficiency plots
 181   TH1F * MatchEffPos = new TH1F("MatchEffPos","",nRuns,0,nRuns);
 182   MatchEffPos->Sumw2();
 183   TH1F * MatchEffNeg = new TH1F("MatchEffNeg","",nRuns,0,nRuns);
 184   MatchEffNeg->Sumw2();
 185
 186   // overall loop over data files (1 per run)
 187   for (Int_t irun=0; irun<nRuns; irun++)
 188     {
 189       Printf("\n--- Processing data from run %i ---\n", runs[irun]);
 190
 191       // open file and print all objects in file
 192       file = TFile::Open(dataFiles[irun].Data());
 193       if (!file)
 194         {
 195           Printf("\n\n!!! ERROR: Could not open file %s !!!\n\n",dataFiles[irun].Data());
 196           continue;
 197         }
 198       file->Print();
 199       Printf("sname: %s", sname.Data());
 200
 201       // choose the right directory, event, and histo managers
 202       Dir = (TDirectoryFile*)file->Get(Form("%s",sname.Data()));
 203       if (!Dir)
 204         {
 205           Printf("\n\n!!! ERROR: Dir is a null pointer. Skipping to next file !!!\n\n");
 206           continue;
 207         }
 208       ecuts = (AliSpectraAODEventCuts*)Dir->Get("Event Cuts");
 209       if (!ecuts)
 210         {
 211           Printf("\n\n!!! ERROR: ecuts is a null pointer. Skipping to next file !!!\n\n");
 212           continue;
 213         }
 214       tcuts = (AliSpectraAODTrackCuts*)Dir->Get("Track Cuts");
 215       if (!tcuts)
 216         {
 217           Printf("\n\n!!! ERROR: tcuts is a null pointer. Skipping to next file !!!\n\n");
 218           continue;
 219         }
 220       hman = (AliSpectraAODHistoManager*)Dir->Get("SpectraHistos");
 221       if (!hman)
 222         {
 223           Printf("\n\n!!! ERROR: hman is a null pointer. Skipping to next file !!!\n\n");
 224           continue;
 225         }
 226
 227       //---------------------------------------------------------------------------------------
 228       // Find the trends in the means and sigmas of fits to the peaks of fixed-Pt projections -
 229       // of the nsigma distributions of the TPC and TOF as a function of the run number       -
 230       //---------------------------------------------------------------------------------------
 231       CheckNSigmaStability(hman,
 232                            TPCnsigMeanTrendPion,
 233                            TPCnsigMeanTrendKaon,
 234                            TPCnsigMeanTrendProton,
 235                            TPCnsigSigmaTrendPion,
 236                            TPCnsigSigmaTrendKaon,
 237                            TPCnsigSigmaTrendProton,
 238                            TOFnsigMeanTrendPion,
 239                            TOFnsigMeanTrendKaon,
 240                            TOFnsigMeanTrendProton,
 241                            TOFnsigSigmaTrendPion,
 242                            TOFnsigSigmaTrendKaon,
 243                            TOFnsigSigmaTrendProton,
 244                            runs,
 245                            nRuns,
 246                            irun,
 247                            useMC);
 248
 249       //---------------------------------------------------------------------------------
 250       // plot the stability of the integrated raw yield as a function of the run number -
 251       //---------------------------------------------------------------------------------
 252       CheckIntegratedRawYield(ecuts,
 253                               hman,
 254                               IntegRawYieldAll,
 255                               IntegRawYieldPiPlus,
 256                               IntegRawYieldPiMinus,
 257                               IntegRawYieldKPlus,
 258                               IntegRawYieldKMinus,
 259                               IntegRawYieldProton,
 260                               IntegRawYieldAntiproton,
 261                               runs,
 262                               irun,
 263                               nRuns,
 264                               Names,
 265                               useMC);
 266
 267       //-----------------------------------------------------------------------------
 268       // find the trend in the MC correction factor as a function of the run number -
 269       //-----------------------------------------------------------------------------
 270       if (useMC) CheckEfficiencies(hman,
 271                                    cEfficienciesAllRuns,
 272                                    EfficiencyPiPlus,
 273                                    EfficiencyKPlus,
 274                                    EfficiencyProton,
 275                                    EfficiencyPiMinus,
 276                                    EfficiencyKMinus,
 277                                    EfficiencyAntiproton,
 278                                    FixedPtEff,
 279                                    runs,
 280                                    nRuns,
 281                                    irun);
 282
 283       //-------------------------------------------------------------------------------
 284       // for a fixed Pt, plot the matching efficiency as a function of the run number -
 285       //-------------------------------------------------------------------------------
 286       CheckMatchingEfficiency(tcuts,
 287                               MatchEffPos,
 288                               MatchEffNeg,
 289                               FixedPtMatchEff,
 290                               runs,
 291                               irun,
 292                               nRuns,
 293                               useMC);
 294
 295       //---------------------------------------------------------------------------
 296       // save eta-phi distributions for each run on a seperate page of a pdf file -
 297       //---------------------------------------------------------------------------
 298       // canvas for temporarily drawing
 299       TCanvas * cEtaPhi = new TCanvas("cEtaPhi","cEtaPhi");
 300       TH2F * hEtaPhi = (TH2F*)((TH2F*)tcuts->GetHistoEtaPhiHighPt()->Clone("hEtaPhi"));
 301       hEtaPhi->Rebin2D(4);
 302       hEtaPhi->Scale(1./hEtaPhi->GetEntries());
 303       gPad->SetGridy();
 304       gPad->SetGridx();
 305       hEtaPhi->DrawClone("COLZ");
 306
 307       // legend
 308       TLegend * lEtaPhi = new TLegend(0.4,0.85,.6,.95);
 309       lEtaPhi->SetFillColor(0);
 310       lEtaPhi->AddEntry(hEtaPhi,Form("Run %i",runs[irun]),"");
 311       if (useMC) lEtaPhi->AddEntry(hEtaPhi,"MC","");
 312       else lEtaPhi->AddEntry(hEtaPhi,"DATA","");
 313       lEtaPhi->DrawClone();
 314
 315       // save to a pdf and close the temporary canvas
 316       if (useMC)
 317         {
 318           if (irun == 0) cEtaPhi->SaveAs("Plots/MC/EtaPhi.pdf(","pdf");
 319           else if (irun < nRuns-1) cEtaPhi->SaveAs("Plots/MC/EtaPhi.pdf","pdf");
 320           else if (irun == nRuns-1) cEtaPhi->SaveAs("Plots/MC/EtaPhi.pdf)","pdf");
 321           cEtaPhi->Close();
 322         }
 323       else
 324         {
 325           if (irun == 0) cEtaPhi->SaveAs("Plots/DATA/EtaPhi.pdf(","pdf");
 326           else if (irun < nRuns-1) cEtaPhi->SaveAs("Plots/DATA/EtaPhi.pdf","pdf");
 327           else if (irun == nRuns-1) cEtaPhi->SaveAs("Plots/DATA/EtaPhi.pdf)","pdf");
 328           cEtaPhi->Close();
 329         }
 330
 331       // close the file and clean up before the next iteration
 332       file->Close();
 333       //delete file;
 334       //delete Dir;
 335       //delete hman;
 336
 337       Printf("\n--- Done Processing data from run %i ---\n", runs[irun]);
 338
 339     }
 340
 341   //------------------------------------------------------------
 342   // Plot all the results and save them to the output file     -
 343   //------------------------------------------------------------
 344   PlotAndSave(runs,
 345               nRuns,
 346               useMC,
 347               fout,
 348
 349               TPCnsigMeanTrendPion,
 350               TPCnsigMeanTrendKaon,
 351               TPCnsigMeanTrendProton,
 352               TPCnsigSigmaTrendPion,
 353               TPCnsigSigmaTrendKaon,
 354               TPCnsigSigmaTrendProton,
 355               TOFnsigMeanTrendPion,
 356               TOFnsigMeanTrendKaon,
 357               TOFnsigMeanTrendProton,
 358               TOFnsigSigmaTrendPion,
 359               TOFnsigSigmaTrendKaon,
 360               TOFnsigSigmaTrendProton,
 361
 362               IntegRawYieldAll,
 363               IntegRawYieldPiPlus,
 364               IntegRawYieldKPlus,
 365               IntegRawYieldProton,
 366               IntegRawYieldPiMinus,
 367               IntegRawYieldKMinus,
 368               IntegRawYieldAntiproton,
 369
 370               EfficiencyPiPlus,
 371               EfficiencyKPlus,
 372               EfficiencyProton,
 373               EfficiencyPiMinus,
 374               EfficiencyKMinus,
 375               EfficiencyAntiproton,
 376
 377               MatchEffPos,
 378               MatchEffNeg);
 379
 380   //------------------------------------------------------------------------------
 381   // Find the outliers from the efficiency, raw yield, matching efficiency, and  -
 382   // nsigma fit plots and print the corresponding run numbers to the console.    -
 383   // Also creates a .txt file in the results/ directory with info about outliers -
 384   //------------------------------------------------------------------------------
 385   FindOutliers(runs,
 386                nRuns,
 387                useMC,
 388                icut,
 389                nSigmaCut,
 390                fout,
 391
 392                TPCnsigMeanTrendPion,
 393                TPCnsigMeanTrendKaon,
 394                TPCnsigMeanTrendProton,
 395                TPCnsigSigmaTrendPion,
 396                TPCnsigSigmaTrendKaon,
 397                TPCnsigSigmaTrendProton,
 398                TOFnsigMeanTrendPion,
 399                TOFnsigMeanTrendKaon,
 400                TOFnsigMeanTrendProton,
 401                TOFnsigSigmaTrendPion,
 402                TOFnsigSigmaTrendKaon,
 403                TOFnsigSigmaTrendProton,
 404
 405                IntegRawYieldAll,
 406
 407                EfficiencyPiPlus,
 408                EfficiencyKPlus,
 409                EfficiencyProton,
 410                EfficiencyPiMinus,
 411                EfficiencyKMinus,
 412                EfficiencyAntiproton,
 413
 414                MatchEffPos,
 415                MatchEffNeg);
 416
 417
 418   // close the output file
 419   fout->Close();
 420   // also close the unused efficiency canvas
 421   if (useMC) cEfficienciesAllRuns->Close();
 422
 423   Printf("\n\n--- End of AODQAChecks() ---\n");
 424 }
 425
 426
 427
 428
 429
 430
 431