PWGDQ/dielectron/macrosJPSI/RaaPbPb2010/RunDstPbPbJpsiAnalysis.C

   1 // J/psi analysis in PbPb
   2 // author: Ionut-Cristian Arsene, i.c.arsene@gsi.de
   3 //         2012/Apr/11
   4
   5 #include <iostream>
   6 using namespace std;
   7
   8 #include <TTimeStamp.h>
   9
  10 #ifndef ALICORRELATIONREDUCEDEVENT_H
  11 #include "AliCorrelationReducedEvent.h"
  12 #endif
  13
  14 #include "DstCommonMacros.C"
  15 using namespace DstCommonMacros;
  16
  17 // function prototypes
  18 void DefineHistograms(const Char_t* histClasses, const Char_t* output);
  19 Bool_t IsEventSelected(AliCorrelationReducedEvent* event);
  20 Bool_t IsLegSelected(AliCorrelationReducedTrack* track, Bool_t* legQualityMap);
  21 Bool_t IsLegSelected(AliCorrelationReducedTrack* track);
  22 Bool_t IsLegSelectedTRD(AliCorrelationReducedTrack* track, Int_t trdCut=-1, Int_t trdMinNtr=4, Float_t eleCut=0.7);
  23 Bool_t IsLegSelectedTOF(AliCorrelationReducedTrack* track, Int_t tofCut=-1, Float_t lowExcl=0.3, Float_t highExcl=0.93);
  24 Bool_t IsTrackUsedForMixing(UShort_t* idxArr, UShort_t size, UShort_t idx);
  25 Bool_t IsPairSelected(AliCorrelationReducedPair* pair);
  26 void RunDstPbPbJpsiAnalysis(const Char_t* output, const Char_t* inputfilename,
  27                             Int_t howMany/* = 5000*/, Int_t offset/* = 0*/);
  28
  29 TFile* gSaveFile=0x0;
  30 // event plane friend file
  31 const Char_t* gkEPfile="dstTree_VZERO_TPC_recentered.root";
  32
  33 // mixing pool depth
  34 const Int_t gkEventMixingPoolDepth = 100;
  35
  36 // centrality binning for the event mixing
  37 const Int_t gkNCentRanges = 3;
  38 /*Double_t gkEMCentLims[gkNCentRanges+1] = { 0.0,  2.5,  5.0,  7.5, 10.0,
  39                                           12.5, 15.0, 17.5, 20.0, 22.5,
  40                                           25.0, 27.5, 30.0, 32.5, 35.0,
  41                                           37.5, 40.0, 50.0, 60.0, 70.0, 80.0};*/
  42 Double_t gkEMCentLims[gkNCentRanges+1] = { 0.0,  10.0, 40.0, 80.0};
  43 // vertex binning for the event mixing
  44 const Int_t gkNVtxRanges = 1;
  45 //Double_t gkEMVtxLims[gkNVtxRanges+1] = {-10.0, -6.0, -2.0, 2.0, 6.0, 10.0};
  46 Double_t gkEMVtxLims[gkNVtxRanges+1] = {-10.0, 10.0};
  47
  48 // event plane angle binning for the event mixing
  49 const Int_t gkNPhiEPRanges = 1;
  50 Double_t gkEMPhiEPLims[gkNPhiEPRanges+1] = {-1.5708, 1.5708};
  51 //Double_t gkEMPhiEPLims[gkNPhiEPRanges+1] = {-1.5708, -1.2566, -0.94248, -0.62832, -0.31415926, 0.0, 0.31415926, 0.62832, 0.94248, 1.2566, 1.5708};
  52 /*Double_t gkEMPhiEPLims[gkNPhiEPRanges+1] = {-1.5708, -1.4137, -1.2566, -1.0996, -0.94248, -0.78540, -0.62832, -0.47124, -0.31416, -0.15708, 0.0,
  53                                              0.15708, 0.31416, 0.47124, 0.62832, 0.7854,   0.94248,  1.09956,  1.2566,   1.4137,   1.5708};*/
  54 // which event plane to be used
  55 const Int_t gkEPused = kVZERORP + 6*1 + 1;    // VZEROC 2nd harmonic event plane
  56 //const Int_t gkEPused = kVZERORP + 6*0 + 1;    // VZEROA 2nd harmonic event plane
  57 //const Int_t gkEPused = kTPCRP + 1;    // TPC harmonic event plane
  58
  59 // define detector cuts
  60 const Int_t gkNDetCuts = 6;
  61 const Char_t* gkLegDetCutNames[gkNDetCuts] = {
  62   "TPC",
  63   "TPC_TRD4_1_0.7", "TPC_TRD4_1_0.8", "TPC_TRD4_1_0.9",
  64   "TPC_TOF", "TPC_TRD4_1_0.9_TOF"
  65 };
  66
  67 //_________________________________________________________________________________________________
  68 void RunDstPbPbJpsiAnalysis(const Char_t* output, const Char_t* inputfilename,
  69                             Int_t howMany/* = 5000*/, Int_t offset/* = 0*/) {
  70   //
  71   // J/psi analysis in Pb-Pb
  72   //
  73   cout << "start ..." << endl;
  74   TTimeStamp start;
  75   cout << "creating chain ..." << endl;
  76   // create the input chains -----------------------------------------
  77   Long64_t entries=0;
  78   TChain* friendChain=0x0;
  79   if(gkEPfile[0]!='\0')
  80     friendChain = new TChain("DstFriendTree");
  81   TChain* chain = GetChain(inputfilename, howMany, offset, entries, friendChain, gkEPfile);
  82   if(!chain) return;
  83   if(gkEPfile[0]!='\0' && !friendChain) return;
  84   AliCorrelationReducedEvent* event = new AliCorrelationReducedEvent();
  85   chain->SetBranchAddress("Event",&event);
  86   AliCorrelationReducedEventFriend* eventF = 0x0;
  87   if(gkEPfile[0]!='\0') {
  88     eventF = new AliCorrelationReducedEventFriend();
  89     friendChain->SetBranchAddress("Event",&eventF);
  90   }
  91
  92   cout << "initialize event mixing lists ..." << endl;
  93
  94   // define histograms -----------------------------------------------
  95   TString histClasses = "";
  96   histClasses += "Event_BeforeCuts;Event_AfterCuts;VZERORP;TPCRP;";
  97   histClasses += "OfflineTriggers;";
  98   histClasses += "TrackingFlags_Before;TrackQA_JpsiLegs_BeforeCuts_ITS_TPC_TRD_TOF;";
  99   for(Int_t iLegCut=0; iLegCut<gkNDetCuts; ++iLegCut) {
 100     histClasses += Form("TrackQA_JpsiLegs_%s_AfterCuts_ITS_TPC_TRD_TOF;TrackingFlags_%s_AfterCuts;",
 101                         gkLegDetCutNames[iLegCut], gkLegDetCutNames[iLegCut]);
 102     for(Int_t iVtx=0; iVtx<gkNVtxRanges; ++iVtx)
 103       histClasses += Form("PairQA_SE_%s_vtx%.1f_%.1f;", gkLegDetCutNames[iLegCut], gkEMVtxLims[iVtx], gkEMVtxLims[iVtx+1]);
 104   }
 105
 106   // initialize event lists for mixing
 107   //const Int_t kNVarMixingRanges = (gkEventMixingType==0 ? 1 : (gkEventMixingType==1 ? gkNVtxRanges : gkNPhiEPRanges));
 108   gEMCategories = gkNDetCuts*gkNCentRanges*gkNPhiEPRanges*gkNVtxRanges;
 109   TList* list1[gkNDetCuts][gkNCentRanges][gkNVtxRanges][gkNPhiEPRanges];  // event master lists
 110   TList* list2[gkNDetCuts][gkNCentRanges][gkNVtxRanges][gkNPhiEPRanges];  // event master lists
 111   TList* selectedPosLegs[gkNDetCuts][gkNCentRanges][gkNVtxRanges][gkNPhiEPRanges][gkEventMixingPoolDepth];
 112   TList* selectedNegLegs[gkNDetCuts][gkNCentRanges][gkNVtxRanges][gkNPhiEPRanges][gkEventMixingPoolDepth];
 113   Int_t mixingPoolSize[gkNDetCuts][gkNCentRanges][gkNVtxRanges][gkNPhiEPRanges] = {{{{0}}}};
 114   UShort_t tracksUsedForMixing[gkNDetCuts][gkNCentRanges][gkNVtxRanges][gkNPhiEPRanges][100] = {{{{{0}}}}};
 115   UShort_t nTracksUsedForMixing[gkNDetCuts][gkNCentRanges][gkNVtxRanges][gkNPhiEPRanges] = {{{{0}}}};
 116   Bool_t eventHasCandidates[gkNDetCuts][gkNCentRanges][gkNVtxRanges][gkNPhiEPRanges] = {{{{kFALSE}}}};
 117   Bool_t fullListFound = kFALSE;
 118   for(Int_t iDetCut=0; iDetCut<gkNDetCuts; ++iDetCut) {
 119     for(Int_t iVtx=0; iVtx<gkNVtxRanges; ++iVtx)
 120       histClasses += Form("PairQA_ME_%s_vtx%.1f_%.1f;", gkLegDetCutNames[iDetCut], gkEMVtxLims[iVtx], gkEMVtxLims[iVtx+1]);
 121   }
 122   for(Int_t iDetCut=0; iDetCut<gkNDetCuts; ++iDetCut) {
 123     for(Int_t iCent=0; iCent<gkNCentRanges; ++iCent) {
 124       for(Int_t iVtx=0; iVtx<gkNVtxRanges; ++iVtx) {
 125         for(Int_t iPhi=0; iPhi<gkNPhiEPRanges; ++iPhi) {
 126           list1[iDetCut][iCent][iVtx][iPhi] = new TList();
 127           list2[iDetCut][iCent][iVtx][iPhi] = new TList();
 128           for(Int_t i=0; i<gkEventMixingPoolDepth; ++i) {
 129             selectedPosLegs[iDetCut][iCent][iVtx][iPhi][i] = new TList();
 130             selectedPosLegs[iDetCut][iCent][iVtx][iPhi][i]->SetOwner();
 131             selectedNegLegs[iDetCut][iCent][iVtx][iPhi][i] = new TList();
 132             selectedNegLegs[iDetCut][iCent][iVtx][iPhi][i]->SetOwner();
 133           }
 134         }  // end loop over phi bins
 135       }  // end loop over vtx bins
 136     }  // end loop over centrality bins
 137   }  // end loop over leg cuts
 138   DefineHistograms(histClasses.Data(), output);
 139
 140   Float_t values[kNVars];
 141   Int_t trackIdMap[20000] = {-1};
 142   TClonesArray* trackList=0x0;
 143   TClonesArray* pairList=0x0;
 144   Bool_t legsQuality[2][gkNDetCuts] = {{kFALSE}};
 145   Bool_t allLegCutsOr[2] = {kFALSE};
 146   Bool_t pairQuality[gkNDetCuts] = {kFALSE};
 147
 148   Float_t oldVertex = -999.; Float_t oldBC=-999.; Float_t oldCentVZERO=-999.; Float_t oldCentSPD=-999.; Float_t oldCentTPC=-999.;
 149   Double_t mixingTime = 0.0;
 150   TTimeStamp startEventLoop;
 151
 152   for(Int_t ie=0; ie<entries; ++ie) {
 153     chain->GetEntry(ie);
 154     friendChain->GetEntry(ie);
 155
 156     gCurrentEvent = event;
 157     gCurrentEventFriend = eventF;
 158     if(ie%100==0) cout << "event " << ie << endl;
 159
 160     FillEventInfo(event, values, eventF);
 161     FillHistClass("Event_BeforeCuts", values);
 162     // event cuts
 163     if(!IsEventSelected(event)) continue;
 164     // check wheter this event is the same as the previous
 165     if(TMath::Abs(event->Vertex(2)-oldVertex)<0.0001 && TMath::Abs(event->CentralityVZERO()-oldCentVZERO)<0.001) {
 166       cout << "This event is the a copy of the previous event" << endl;
 167       cout << "OLD/NEW:  vtxZ " << oldVertex << "/" << event->Vertex(2)
 168            << "  BC " << oldBC << "/" << event->BC()
 169            << "  CentVZERO " << oldCentVZERO << "/" << event->CentralityVZERO()
 170            << "  CentSPD " << oldCentSPD << "/" << event->CentralitySPD()
 171            << "  CentTPC " << oldCentTPC << "/" << event->CentralityTPC() << endl;
 172       ++ie;
 173       continue;
 174     }
 175     else {
 176       oldVertex = event->Vertex(2); oldBC = event->BC();
 177       oldCentVZERO = event->CentralityVZERO(); oldCentSPD = event->CentralitySPD();
 178       oldCentTPC = event->CentralityTPC();
 179     }
 180
 181     for(UShort_t ibit=0; ibit<64; ++ibit) {
 182       FillEventOfflineTriggers(ibit, values);
 183       FillHistClass("OfflineTriggers", values);
 184     }
 185
 186     // get the event vtx and centrality bins
 187     Float_t centVZERO = values[kCentVZERO];
 188     Float_t vtxZ = values[kVtxZ];
 189     Int_t binCent = -1;
 190     for(Int_t i=0; i<gkNCentRanges; ++i) {
 191       if(centVZERO>=gkEMCentLims[i] && centVZERO<gkEMCentLims[i+1]) {
 192         binCent = i; break;
 193       }
 194     }
 195     Int_t binVtxZ = -1;
 196     for(Int_t i=0; i<gkNVtxRanges; ++i) {
 197       if(vtxZ>=gkEMVtxLims[i] && vtxZ<gkEMVtxLims[i+1]) {
 198         binVtxZ = i; break;
 199       }
 200     }
 201     if(binCent==-1) {
 202       cout << "Warning: Centrality bin for this event is -1!! Something went wrong, check it out!" << endl;
 203       cout << "centVZERO = " << centVZERO << endl;
 204       continue;
 205     }
 206     if(binVtxZ==-1) {
 207       cout << "Warning: Vertex bin for this event is -1!! Something went wrong, check it out!" << endl;
 208       cout << "vtxZ = " << vtxZ << endl;
 209       continue;
 210     }
 211
 212     Float_t phiEP = values[gkEPused];
 213     Int_t binPhiEP = -1;
 214     for(Int_t i=0; i<gkNPhiEPRanges; ++i) {
 215       if(values[gkEPused]>=gkEMPhiEPLims[i] &&
 216         values[gkEPused]<gkEMPhiEPLims[i+1]) {
 217           binPhiEP = i; break;
 218       }
 219     }
 220     if(binPhiEP==-1) {
 221       cout << "Warning: EP Phi bin for this event is -1!! Something went wrong, check it out!" << endl;
 222       cout << "phi = " << values[gkEPused] << endl;
 223       continue;
 224     }
 225
 226     // Do the mixing if the events in the pool reached the maximum number of events
 227     if(fullListFound) {
 228       TTimeStamp startMix;
 229       for(Int_t iCut=0; iCut<gkNDetCuts; ++iCut) {
 230         for(Int_t iCent=0; iCent<gkNCentRanges; ++iCent) {
 231           for(Int_t iVtx=0; iVtx<gkNVtxRanges; ++iVtx) {
 232             for(Int_t iPhi=0; iPhi<gkNPhiEPRanges; ++iPhi) {
 233               if(mixingPoolSize[iCut][iCent][iVtx][iPhi]>=gkEventMixingPoolDepth) {
 234                 values[kCentVZERO] = 0.5*(gkEMCentLims[iCent]+gkEMCentLims[iCent+1]);
 235                 values[kVtxZ] = 0.5*(gkEMVtxLims[iVtx]+gkEMVtxLims[iVtx+1]);
 236                 values[gkEPused] = 0.5*(gkEMPhiEPLims[iPhi]+gkEMPhiEPLims[iPhi+1]);
 237                 DoEventMixing(list1[iCut][iCent][iVtx][iPhi], list2[iCut][iCent][iVtx][iPhi], values,
 238                               AliCorrelationReducedPair::kJpsiToEE,
 239                               Form("PairQA_ME_%s_vtx%.1f_%.1f", gkLegDetCutNames[iCut], gkEMVtxLims[iVtx], gkEMVtxLims[iVtx+1]));
 240                 mixingPoolSize[iCut][iCent][iVtx][iPhi] = 0;   // reset the mixing events number
 241               }
 242             }  // end loop over phi ranges
 243           }   // end loop over vertex ranges
 244         }   // end loop over centralities ranges
 245         TTimeStamp stopMix;
 246         mixingTime += Double_t(stopMix.GetSec())+Double_t(stopMix.GetNanoSec())/1.0e+9 -
 247                      (Double_t(startMix.GetSec())+Double_t(startMix.GetNanoSec())/1.0e+9);
 248       }   // end loop over detector cuts
 249       fullListFound = kFALSE;
 250     }  // end if (fullListFound)
 251     // reset the correct centrality, vtx and phi for this event
 252     values[kCentVZERO] = centVZERO;
 253     values[kVtxZ] = vtxZ;
 254     values[gkEPused] = phiEP;
 255
 256     // make an index map for track access optimization
 257     AliCorrelationReducedTrack* track = 0x0;
 258     trackList = event->GetTracks();
 259     TIter nextTrack(trackList);
 260     for(Int_t i=0; i<20000; ++i) trackIdMap[i] = -1;
 261     for(Int_t it=0; it<event->NTracks(); ++it) {
 262       track = (AliCorrelationReducedTrack*)nextTrack();
 263       if(track)
 264         trackIdMap[track->TrackId()] = it;
 265     }
 266
 267     pairList = event->GetPairs();
 268     TIter nextPair(pairList);
 269     AliCorrelationReducedPair* pair = 0x0;
 270     AliCorrelationReducedTrack* leg1 = 0x0;
 271     AliCorrelationReducedTrack* leg2 = 0x0;
 272     values[kNpairsSelected] = 0.0;
 273
 274     // reset arrays
 275     for(Int_t iCut=0; iCut<gkNDetCuts; ++iCut) {
 276       for(Int_t iCent=0; iCent<gkNCentRanges; ++iCent) {
 277         for(Int_t iVtx=0; iVtx<gkNVtxRanges; ++iVtx) {
 278           for(Int_t iPhi=0; iPhi<gkNPhiEPRanges; ++iPhi) {
 279             nTracksUsedForMixing[iCut][iCent][iVtx][iPhi] = 0;
 280             eventHasCandidates[iCut][iCent][iVtx][iPhi] = kFALSE;
 281           }  // end loop over phi bins
 282         }  // end loop over vtx bins
 283       }  // end loop over centrality bins
 284     }  // end loop over detector cuts
 285
 286     // loop over pairs
 287     for(Int_t ip=0; ip<event->NDielectrons(); ++ip) {
 288       pair = event->GetDielectronPair(ip);
 289       // pair cuts
 290       if(!IsPairSelected(pair)) continue;
 291
 292       // reset flag arrays
 293       for(Int_t iCut=0; iCut<gkNDetCuts; ++iCut) {
 294         legsQuality[0][iCut] = kFALSE; legsQuality[1][iCut] = kFALSE;
 295         pairQuality[iCut] = kFALSE;
 296       }
 297       allLegCutsOr[0] = kFALSE; allLegCutsOr[1] = kFALSE;
 298       // apply cuts on the first leg
 299       if(trackIdMap[pair->LegId(0)]!=-1) {
 300         leg1 = event->GetTrack(trackIdMap[pair->LegId(0)]);
 301         FillTrackInfo(leg1, values);
 302         for(UShort_t iflag=0; iflag<kNTrackingFlags; ++iflag) {
 303           FillTrackingFlag(leg1, iflag, values);
 304           FillHistClass("TrackingFlags_Before", values);
 305         }
 306         FillHistClass("TrackQA_JpsiLegs_BeforeCuts_ITS_TPC_TRD_TOF", values);
 307         allLegCutsOr[0] = IsLegSelected(leg1, legsQuality[0]);
 308       }   // end if
 309
 310       // apply cuts on the second leg
 311       if(trackIdMap[pair->LegId(1)]!=-1) {
 312         leg2 = event->GetTrack(trackIdMap[pair->LegId(1)]);
 313         FillTrackInfo(leg2, values);
 314         for(UShort_t iflag=0; iflag<kNTrackingFlags; ++iflag) {
 315           FillTrackingFlag(leg2, iflag, values);
 316           FillHistClass("TrackingFlags_Before", values);
 317         }
 318         FillHistClass("TrackQA_JpsiLegs_BeforeCuts_ITS_TPC_TRD_TOF", values);
 319         allLegCutsOr[1] = IsLegSelected(leg2, legsQuality[1]);
 320       }   // end if
 321       if(!allLegCutsOr[0] || !allLegCutsOr[1]) continue;
 322       Int_t nCutsPassed = 0;
 323       for(Int_t iCut=0; iCut<gkNDetCuts; ++iCut) {
 324         pairQuality[iCut] = legsQuality[0][iCut] && legsQuality[1][iCut];
 325         if(pairQuality[iCut]) ++nCutsPassed;
 326       }
 327       if(nCutsPassed==0) continue;
 328
 329       FillPairInfo(pair, values);
 330
 331       // fill track leg histograms
 332       FillTrackInfo(leg1, values);
 333       for(Int_t iCut=0; iCut<gkNDetCuts; ++iCut) {
 334         if(pairQuality[iCut]) {
 335           if(!IsTrackUsedForMixing(tracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP],
 336                                    nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP],
 337                                    leg1->TrackId())) {
 338             FillHistClass(Form("TrackQA_JpsiLegs_%s_AfterCuts_ITS_TPC_TRD_TOF", gkLegDetCutNames[iCut]), values);
 339             for(UShort_t iflag=0; iflag<kNTrackingFlags; ++iflag) {
 340               FillTrackingFlag(leg1, iflag, values);
 341               FillHistClass(Form("TrackingFlags_%s_AfterCuts", gkLegDetCutNames[iCut]), values);
 342             }
 343           }
 344         }
 345       }
 346       FillTrackInfo(leg2, values);
 347       for(Int_t iCut=0; iCut<gkNDetCuts; ++iCut) {
 348         if(pairQuality[iCut]) {
 349           if(!IsTrackUsedForMixing(tracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP],
 350                                    nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP],
 351                                    leg2->TrackId())) {
 352             FillHistClass(Form("TrackQA_JpsiLegs_%s_AfterCuts_ITS_TPC_TRD_TOF", gkLegDetCutNames[iCut]), values);
 353             for(UShort_t iflag=0; iflag<kNTrackingFlags; ++iflag) {
 354               FillTrackingFlag(leg2, iflag, values);
 355               FillHistClass(Form("TrackingFlags_%s_AfterCuts", gkLegDetCutNames[iCut]), values);
 356             }
 357           }
 358         }
 359       }
 360
 361       // fill single-event pair histograms
 362       for(Int_t iCut=0; iCut<gkNDetCuts; ++iCut) {
 363         if(pairQuality[iCut]) {
 364           if(IsPairSelectedEM(values))
 365             FillHistClass(Form("PairQA_SE_%s_vtx%.1f_%.1f", gkLegDetCutNames[iCut], gkEMVtxLims[binVtxZ], gkEMVtxLims[binVtxZ+1]), values);
 366         }
 367       }
 368
 369       values[kNpairsSelected] += 1.0;
 370
 371       // add legs to the mixing lists
 372       for(Int_t iCut=0; iCut<gkNDetCuts; ++iCut) {
 373         if(pairQuality[iCut]) {
 374           if(!IsTrackUsedForMixing(tracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP],
 375                                    nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP],
 376                                    leg1->TrackId())) {
 377             if(leg1->Charge()>0) selectedPosLegs[iCut][binCent][binVtxZ][binPhiEP][mixingPoolSize[iCut][binCent][binVtxZ][binPhiEP]]->Add(leg1->Clone());
 378             else selectedNegLegs[iCut][binCent][binVtxZ][binPhiEP][mixingPoolSize[iCut][binCent][binVtxZ][binPhiEP]]->Add(leg1->Clone());
 379             tracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP][nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP]] = leg1->TrackId();
 380             nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP] += 1;
 381           }
 382           if(!IsTrackUsedForMixing(tracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP],
 383                                    nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP],
 384                                    leg2->TrackId())) {
 385             if(leg2->Charge()>0) selectedPosLegs[iCut][binCent][binVtxZ][binPhiEP][mixingPoolSize[iCut][binCent][binVtxZ][binPhiEP]]->Add(leg2->Clone());
 386             else selectedNegLegs[iCut][binCent][binVtxZ][binPhiEP][mixingPoolSize[iCut][binCent][binVtxZ][binPhiEP]]->Add(leg2->Clone());
 387             tracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP][nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP]] = leg2->TrackId();
 388             nTracksUsedForMixing[iCut][binCent][binVtxZ][binPhiEP] += 1;
 389           }
 390           eventHasCandidates[iCut][binCent][binVtxZ][binPhiEP] = kTRUE;
 391         }  // end if(pairQuality)
 392       }
 393     }  // end loop over pairs
 394
 395     // If this event has candidates, then store the event in the list for event mixing
 396     for(Int_t iCut=0; iCut<gkNDetCuts; ++iCut) {
 397       for(Int_t iCent=0; iCent<gkNCentRanges; ++iCent) {
 398         for(Int_t iVtx=0; iVtx<gkNVtxRanges; ++iVtx) {
 399           for(Int_t iPhi=0; iPhi<gkNPhiEPRanges; ++iPhi) {
 400             if(!eventHasCandidates[iCut][iCent][iVtx][iPhi]) continue;
 401             list1[iCut][iCent][iVtx][iPhi]->Add(selectedPosLegs[iCut][iCent][iVtx][iPhi][mixingPoolSize[iCut][iCent][iVtx][iPhi]]);
 402             list2[iCut][iCent][iVtx][iPhi]->Add(selectedNegLegs[iCut][iCent][iVtx][iPhi][mixingPoolSize[iCut][iCent][iVtx][iPhi]]);
 403             /*
 404             cout << "event mixing category (cut/cent/vtx) " << iCut << "/" << iCent << "/" << iVtx << ", pool size (pos/neg): "
 405                  << list1[iCut][iCent][iVtx]->GetEntries()
 406                  << "/" << list2[iCut][iCent][iVtx]->GetEntries() << "; ntracks this entry (pos/neg): "
 407                  << selectedPosLegs[iCut][iCent][iVtx][mixingPoolSize[iCut][iCent][iVtx]]->GetEntries() << "/"
 408                  << selectedNegLegs[iCut][iCent][iVtx][mixingPoolSize[iCut][iCent][iVtx]]->GetEntries() << endl;
 409             */
 410             mixingPoolSize[iCut][iCent][iVtx][iPhi] += 1;
 411             if(mixingPoolSize[iCut][iCent][iVtx][iPhi]==gkEventMixingPoolDepth)
 412               fullListFound = kTRUE;
 413           }  // end loop over phi ranges
 414         }  // end loop over vertex ranges
 415       }  // end loop over centrality ranges
 416     }  // end loop over detector cuts
 417
 418     FillHistClass("Event_AfterCuts", values);
 419     FillHistClass("TPCRP", values);
 420     FillHistClass("VZERORP", values);
 421   }  // end loop over events
 422
 423   // Mixing the leftover events
 424   cout << "Leftover mixing ..." << endl;
 425   TTimeStamp startMix;
 426   for(Int_t iCut=0; iCut<gkNDetCuts; ++iCut) {
 427     for(Int_t iCent=0; iCent<gkNCentRanges; ++iCent) {
 428       for(Int_t iVtx=0; iVtx<gkNVtxRanges; ++iVtx) {
 429         for(Int_t iPhi=0; iPhi<gkNPhiEPRanges; ++iPhi) {
 430           values[kCentVZERO] = 0.5*(gkEMCentLims[iCent]+gkEMCentLims[iCent+1]);
 431           values[kVtxZ] = 0.5*(gkEMVtxLims[iVtx]+gkEMVtxLims[iVtx+1]);
 432           values[gkEPused] = 0.5*(gkEMPhiEPLims[iPhi]+gkEMPhiEPLims[iPhi+1]);
 433           DoEventMixing(list1[iCut][iCent][iVtx][iPhi], list2[iCut][iCent][iVtx][iPhi], values,
 434                         AliCorrelationReducedPair::kJpsiToEE,
 435                         Form("PairQA_ME_%s_vtx%.1f_%.1f", gkLegDetCutNames[iCut], gkEMVtxLims[iVtx], gkEMVtxLims[iVtx+1])); // after the event mixing the lists will be cleared
 436           mixingPoolSize[iCut][iCent][iVtx][iPhi] = 0;   // reset the mixing events number
 437         }  // end loop over phi ranges
 438       }   // end loop over vertex ranges
 439     }   // end loop over centralities ranges
 440   }   // end loop over detector cuts
 441   TTimeStamp stopMix;
 442   mixingTime += Double_t(stopMix.GetSec())+Double_t(stopMix.GetNanoSec())/1.0e+9 -
 443                (Double_t(startMix.GetSec())+Double_t(startMix.GetNanoSec())/1.0e+9);
 444
 445   TTimeStamp stopEventLoop;
 446
 447   WriteOutput(gSaveFile);
 448
 449   cout << "Initialization time: " << startEventLoop.GetSec() - start.GetSec() << " seconds" << endl;
 450   cout << "Total looping time: " << stopEventLoop.GetSec() - startEventLoop.GetSec() << " seconds" << endl;
 451   cout << "Time spent in mixing: " << mixingTime << " seconds" << endl;
 452   cout << "Overall speed       : " << Double_t(stopEventLoop.GetSec() - startEventLoop.GetSec())/Double_t(entries) << " sec./event" << endl;
 453 }
 454
 455
 456 //__________________________________________________________________
 457 Bool_t IsTrackUsedForMixing(UShort_t* idxArray, UShort_t size, UShort_t idx) {
 458   //
 459   // check whether track with idx was already added to the mixing list
 460   //
 461   for(Int_t i=0; i<size; ++i) {
 462     if(idxArray[i]==idx) return kTRUE;
 463   }
 464   return kFALSE;
 465 }
 466
 467
 468 //__________________________________________________________________
 469 Bool_t IsEventSelected(AliCorrelationReducedEvent* event) {
 470   //
 471   // event selection
 472   //
 473   if(!event->IsPhysicsSelection()) return kFALSE;
 474   if(event->VertexNContributors()<1) return kFALSE;
 475   if(event->CentralityQuality()!=0) return kFALSE;
 476   if(event->CentralityVZERO()<-0.0001) return kFALSE;
 477   if(event->CentralityVZERO()>80.0) return kFALSE;
 478   if(TMath::Abs(event->Vertex(2))>10.0) return kFALSE;
 479   return kTRUE;
 480 }
 481
 482
 483 //__________________________________________________________________
 484 Bool_t IsLegSelected(AliCorrelationReducedTrack* track, Bool_t* legQuality) {
 485   //
 486   //  track leg selection
 487   //
 488   legQuality[0] = IsLegSelected(track);
 489   legQuality[1] = legQuality[0] && IsLegSelectedTRD(track, 1, 4, 0.7);
 490   legQuality[2] = legQuality[0] && IsLegSelectedTRD(track, 1, 4, 0.80);
 491   legQuality[3] = legQuality[0] && IsLegSelectedTRD(track, 1, 4, 0.90);
 492   legQuality[4] = legQuality[0] && IsLegSelectedTOF(track, 2);
 493   legQuality[5] = legQuality[3] && legQuality[4];
 494   Bool_t globalOr = kFALSE;
 495   for(Int_t i=0; i<gkNDetCuts; ++i) globalOr = globalOr || legQuality[i];
 496   return globalOr;
 497 }
 498
 499
 500 //__________________________________________________________________
 501 Bool_t IsLegSelected(AliCorrelationReducedTrack* track) {
 502   //
 503   // track leg selection
 504   //
 505   if(!track->CheckTrackStatus(kTPCrefit)) return kFALSE;
 506   if(!track->CheckTrackStatus(kITSrefit)) return kFALSE;
 507   if(!(track->ITSLayerHit(0) || track->ITSLayerHit(1))) return kFALSE;
 508
 509   if(track->Pt()<0.7) return kFALSE;
 510   if(TMath::Abs(track->Eta())>0.9) return kFALSE;
 511
 512   if(track->TPCncls()<70) return kFALSE;
 513
 514   if(track->TPCnSig(kElectron)>3.0) return kFALSE;
 515   if(track->TPCnSig(kElectron)<-2.0) return kFALSE;
 516
 517   if(TMath::Abs(track->DCAz())>3.0) return kFALSE;
 518   if(TMath::Abs(track->DCAxy())>1.0) return kFALSE;
 519
 520   //if(track->Pin()<2.5) {
 521   if(TMath::Abs(track->TPCnSig(kProton))<3.5) return kFALSE;
 522   if(TMath::Abs(track->TPCnSig(kPion))<3.5) return kFALSE;
 523   if(track->TPCnSig(kProton)<-3.5) return kFALSE;
 524   //if(track->TPCnSig(kProton)<-3.5) return kFALSE;
 525   //}
 526
 527   return kTRUE;
 528 }
 529
 530
 531 //__________________________________________________________________
 532 Bool_t IsLegSelectedTRD(AliCorrelationReducedTrack* track,
 533                         Int_t trdCut/*=-1*/, Int_t trdMinNtr/*=4*/, Float_t eleCut/*=0.7*/) {
 534   //
 535   // TRD leg selection
 536   //
 537   if(trdCut==1 && track->TRDntracklets(0)>=trdMinNtr) {
 538     if(track->TRDpid(0)<eleCut) return kFALSE;
 539   }
 540   if(trdCut==2 && track->TRDntracklets(0)>=trdMinNtr && track->Pin()>1.0) {
 541     if(track->TRDpid(0)<eleCut) return kFALSE;
 542   }
 543   if(trdCut==3 && track->CheckTrackStatus(kTRDpid) && track->TRDntracklets(0)>=trdMinNtr) {
 544     if(track->TRDpid(0)<eleCut) return kFALSE;
 545   }
 546   if(trdCut==4 && track->CheckTrackStatus(kTRDpid) && track->TRDntracklets(1)>=trdMinNtr && track->Pin()>1.0) {
 547     if(track->TRDpid(0)<eleCut) return kFALSE;
 548   }
 549   return kTRUE;
 550 }
 551
 552
 553 //__________________________________________________________________
 554 Bool_t IsLegSelectedTOF(AliCorrelationReducedTrack* track,
 555                         Int_t tofCut/*=-1*/, Float_t lowExcl/*=0.3*/, Float_t highExcl/*=0.93*/) {
 556   //
 557   // TOF leg selection
 558   //
 559   if(tofCut==1 && track->TOFbeta()>lowExcl && track->TOFbeta()<highExcl) return kFALSE;
 560   if(tofCut==2 && track->CheckTrackStatus(kTOFpid) && TMath::Abs(track->TOFnSig(kElectron))>3.5) return kFALSE;
 561   return kTRUE;
 562 }
 563
 564
 565 //__________________________________________________________________
 566 Bool_t IsPairSelected(AliCorrelationReducedPair* pair) {
 567   //
 568   // pair selection
 569   //
 570   if(!pair) return kFALSE;
 571   if(pair->CandidateId()!=AliCorrelationReducedPair::kJpsiToEE) return kFALSE;
 572   if(pair->PairType()!=1) return kFALSE;
 573   if(TMath::Abs(pair->Rapidity())>0.9) return kFALSE;
 574   return kTRUE;
 575 }
 576
 577
 578 //__________________________________________________________________
 579 void DefineHistograms(const Char_t* histClasses, const Char_t* output) {
 580   //
 581   // define the histograms
 582   //
 583   cout << "Defining histograms ..." << flush;
 584   cout << "histogram classes: " << histClasses << endl;
 585
 586   gSaveFile=new TFile(output,"RECREATE");
 587
 588   TString classesStr(histClasses);
 589   TObjArray* arr=classesStr.Tokenize(";");
 590
 591   const Int_t kNRunBins = 3000;
 592   Double_t runHistRange[2] = {137000.,140000.};
 593
 594   for(Int_t iclass=0; iclass<arr->GetEntries(); ++iclass) {
 595     TString classStr = arr->At(iclass)->GetName();
 596     cout << "hist class: " << classStr.Data() << endl;
 597
 598     // Event wise histograms
 599     if(classStr.Contains("Event")) {
 600       cout << "Event" << endl;
 601       AddHistClass(classStr.Data());
 602       AddHistogram(classStr.Data(),"RunNo","Run numbers;Run",kFALSE, kNRunBins, runHistRange[0], runHistRange[1], kRunNo);
 603       AddHistogram(classStr.Data(),"BC","Bunch crossing;BC",kFALSE,3000,0.,3000.,kBC);
 604       AddHistogram(classStr.Data(),"IsPhysicsSelection","Physics selection flag;;",kFALSE,
 605                    2,-0.5,1.5,kIsPhysicsSelection, 0,0.0,0.0,kNothing, 0,0.0,0.0,kNothing, "off;on");
 606
 607       AddHistogram(classStr.Data(),"VtxZ","Vtx Z;vtx Z (cm)",kFALSE,300,-15.,15.,kVtxZ);
 608       AddHistogram(classStr.Data(),"VtxX","Vtx X;vtx X (cm)",kFALSE,300,-0.4,0.4,kVtxX);
 609       AddHistogram(classStr.Data(),"VtxY","Vtx Y;vtx Y (cm)",kFALSE,300,-0.4,0.4,kVtxY);
 610
 611       AddHistogram(classStr.Data(),"VtxZ_Run_prof", "<VtxZ> vs run; Run; vtxZ (cm)", kTRUE,
 612                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, -20, 20., 1000., kVtxZ);
 613       AddHistogram(classStr.Data(),"VtxX_Run_prof", "<VtxX> vs run; Run; vtxX (cm)", kTRUE,
 614                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, -20, 20., 1000., kVtxX);
 615       AddHistogram(classStr.Data(),"VtxY_Run_prof", "<VtxY> vs run; Run; vtxY (cm)", kTRUE,
 616                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, -20, 20., 1000., kVtxY);
 617
 618       AddHistogram(classStr.Data(),"CentVZERO","Centrality(VZERO);centrality VZERO (percents)",kFALSE,
 619                    100, 0.0, 100.0, kCentVZERO);
 620       AddHistogram(classStr.Data(),"CentQuality","Centrality quality;centrality quality",kFALSE,
 621                    100, -50.5, 49.5, kCentQuality);
 622       AddHistogram(classStr.Data(),"CentVZERO_Run_prof","<Centrality(VZERO)> vs run;Run; centrality VZERO (%)",kTRUE,
 623                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0.0, 100.0, kCentVZERO);
 624
 625       AddHistogram(classStr.Data(),"NPairs","Number of candidates per event;# pairs",kFALSE,
 626                    5000,0.,5000.,kNdielectrons);
 627       AddHistogram(classStr.Data(),"NPairsSelected", "Number of selected pairs per event; #pairs", kFALSE,
 628                    5000,0.,5000.,kNpairsSelected);
 629       AddHistogram(classStr.Data(),"NTracksTotal","Number of total tracks per event;# tracks",kFALSE,
 630                    1000,0.,30000.,kNtracksTotal);
 631       AddHistogram(classStr.Data(),"NTracksSelected","Number of selected tracks per event;# tracks",kFALSE,
 632                    1000,0.,30000.,kNtracksSelected);
 633       AddHistogram(classStr.Data(),"SPDntracklets", "SPD #tracklets in |#eta|<1.0; tracklets", kFALSE,
 634                    3000, -0.5, 2999.5, kSPDntracklets);
 635
 636       AddHistogram(classStr.Data(),"Ndielectrons_Run_prof", "<Number of dielectrons> per run; Run; #tracks", kTRUE,
 637                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0., 10000., kNdielectrons);
 638       AddHistogram(classStr.Data(),"NpairsSelected_Run_prof", "<Number of selected pairs> per run; Run; #tracks", kTRUE,
 639                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0., 10000., kNpairsSelected);
 640       AddHistogram(classStr.Data(),"NTracksTotal_Run_prof", "<Number of tracks> per run; Run; #tracks", kTRUE,
 641                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0., 10000., kNtracksTotal);
 642       AddHistogram(classStr.Data(),"NTracksSelected_Run_prof", "<Number of selected tracks> per run; Run; #tracks", kTRUE,
 643                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0., 10000., kNtracksSelected);
 644       AddHistogram(classStr.Data(),"SPDntracklets_Run_prof", "<SPD ntracklets> per run; Run; #tracks", kTRUE,
 645                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, 0., 10000., kSPDntracklets);
 646
 647       AddHistogram(classStr.Data(),"VtxZ_CentVZERO","Centrality(VZERO) vs vtx. Z;vtx Z (cm); centrality VZERO (%)",kFALSE,
 648                    300,-15.,15.,kVtxZ, 100, 0.0, 100.0, kCentVZERO);
 649       AddHistogram(classStr.Data(),"VtxZ_CentSPD","Centrality(SPD) vs vtx. Z;vtx Z (cm); centrality SPD (%)",kFALSE,
 650                    300,-15.,15.,kVtxZ, 100, 0.0, 100.0, kCentSPD);
 651       AddHistogram(classStr.Data(),"VtxZ_CentTPC","Centrality(TPC) vs vtx. Z;vtx Z (cm); centrality TPC (%)",kFALSE,
 652                    300,-15.,15.,kVtxZ, 100, 0.0, 100.0, kCentTPC);
 653       continue;
 654     }  // end if className contains "Event"
 655
 656     // Offline trigger histograms
 657     if(classStr.Contains("OfflineTriggers")) {
 658       cout << "OfflineTriggers" << endl;
 659       AddHistClass(classStr.Data());
 660
 661       TString triggerNames = "";
 662       for(Int_t i=0; i<64; ++i) {triggerNames += gkOfflineTriggerNames[i]; triggerNames+=";";}
 663
 664       AddHistogram(classStr.Data(), "Triggers", "Offline triggers fired; ; ;", kFALSE,
 665                    64, -0.5, 63.5, kOfflineTrigger, 2, -0.5, 1.5, kOfflineTriggerFired, 0, 0.0, 0.0, kNothing, triggerNames.Data(), "off;on");
 666       AddHistogram(classStr.Data(), "Triggers2", "Offline triggers fired; ; ;", kFALSE,
 667                    64, -0.5, 63.5, kOfflineTriggerFired2, 0, 0.0, 0.0, kNothing, 0, 0.0, 0.0, kNothing, triggerNames.Data());
 668       AddHistogram(classStr.Data(), "CentVZERO_Triggers2", "Offline triggers fired vs centrality VZERO; ; centrality VZERO;", kFALSE,
 669                    64, -0.5, 63.5, kOfflineTriggerFired2, 20, 0.0, 100.0, kCentVZERO, 0, 0.0, 0.0, kNothing, triggerNames.Data());
 670       AddHistogram(classStr.Data(), "VtxZ_Triggers2", "Offline triggers fired vs vtxZ; ; vtx Z (cm.);", kFALSE,
 671                    64, -0.5, 63.5, kOfflineTriggerFired2, 200, -20.0, 20.0, kVtxZ, 0, 0.0, 0.0, kNothing, triggerNames.Data());
 672       continue;
 673     }
 674
 675     if(classStr.Contains("VZERORP")) {
 676       cout << "VZERORP" << endl;
 677       AddHistClass(classStr.Data());
 678       const Char_t* sname[3] = {"A","C","A&C"};
 679       for(Int_t ih=1; ih<2; ++ih) {
 680         for(Int_t iS=0; iS<3; ++iS) {
 681           AddHistogram(classStr.Data(), Form("QvecX_side%s_h%d_CentSPDVtxZ_prof",sname[iS],ih+1),
 682                        Form("Q_{x}, side %s, harmonic %d, vs centSPD and vtxZ; Centrality (percents); VtxZ (cm); <Q_{x}>",sname[iS],ih+1), kTRUE,
 683                        20, 0.0, 100.0, kCentSPD, 24, -12.0, 12.0, kVtxZ, 500, -10000.0, 10000.0, kVZEROQvecX+iS*6+ih);
 684           AddHistogram(classStr.Data(), Form("QvecY_side%s_h%d_CentSPDVtxZ_prof",sname[iS],ih+1),
 685                        Form("Q_{y}, side %s, harmonic %d, vs centSPD and vtxZ; Centrality (percents); VtxZ (cm); <Q_{y}>",sname[iS],ih+1), kTRUE,
 686                        20, 0.0, 100.0, kCentSPD, 24, -12.0, 12.0, kVtxZ, 500, -10000.0, 10000.0, kVZEROQvecY+iS*6+ih);
 687           AddHistogram(classStr.Data(), Form("QvecX_side%s_h%d_Run_prof", sname[iS], ih+1),
 688                        Form("<Q_{x}>, VZERO side %s, harmonic %d, vs run; Run; <Q_{x}>", sname[iS], ih+1), kTRUE,
 689                        kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, -100., 100., kVZEROQvecX+iS*6+ih);
 690           AddHistogram(classStr.Data(), Form("QvecY_side%s_h%d_Run_prof", sname[iS], ih+1),
 691                        Form("<Q_{y}>, VZERO side %s, harmonic %d, vs run; Run; <Q_{y}>", sname[iS], ih+1), kTRUE,
 692                        kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, -100., 100., kVZEROQvecY+iS*6+ih);
 693           AddHistogram(classStr.Data(), Form("RP_side%s_h%d_CentSPD",sname[iS],ih+1),
 694                        Form("VZERO reaction plane, side %s, harmonic %d, vs centrality SPD; #Psi (rad.); centSPD (percents)",sname[iS],ih+1), kFALSE,
 695                        400, -4.0/Double_t(ih+1), 4.0/Double_t(ih+1), kVZERORP+iS*6+ih, 20, 0.0, 100.0, kCentSPD);
 696           AddHistogram(classStr.Data(), Form("RP_side%s_h%d_VtxZ",sname[iS],ih+1),
 697                        Form("VZERO reaction plane, side %s, harmonic %d, vs vtxZ; #Psi (rad.); vtxZ (cm)",sname[iS],ih+1), kFALSE,
 698                        400, -4.0/Double_t(ih+1), 4.0/Double_t(ih+1), kVZERORP+iS*6+ih, 24, -12.0, +12.0, kVtxZ);
 699         }   // end loop over VZERO sides
 700       }   // end loop over harmonics
 701       continue;
 702     }   // end if for the VZERO reaction plane histograms
 703
 704     if(classStr.Contains("TPCRP")) {
 705       cout << "TPCRP" << endl;
 706       AddHistClass(classStr.Data());
 707       for(Int_t ih=1; ih<2; ++ih) {
 708         AddHistogram(classStr.Data(), Form("QvecX_TPC_h%d_Run_prof", ih+1),
 709                      Form("<Q_{x}>, TPC, harmonic %d, vs run; Run; <Q_{x}>", ih+1), kTRUE,
 710                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, -100., 100., kTPCQvecX+ih);
 711         AddHistogram(classStr.Data(), Form("QvecY_TPC_h%d_Run_prof", ih+1),
 712                      Form("<Q_{y}>, TPC, harmonic %d, vs run; Run; <Q_{y}>", ih+1), kTRUE,
 713                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 100, -100., 100., kTPCQvecY+ih);
 714         AddHistogram(classStr.Data(), Form("TPCRP_h%d", ih+1),
 715                      Form("TPC event plane, harmonic %d; #Psi (rad.)", ih+1), kFALSE,
 716                      400, -4.0/Double_t(ih+1), 4.0/Double_t(ih+1), kTPCRP+ih);
 717         AddHistogram(classStr.Data(), Form("TPCQvecX_h%d_CentSPDVtxZ_prof", ih+1),
 718                      Form("TPC Q_{x}, harmonic %d, vs centSPD and vtxZ; Centrality (percents); VtxZ (cm); <Q_{x}>", ih+1), kTRUE,
 719                      20, 0.0, 100.0, kCentSPD, 24, -12.0, 12.0, kVtxZ, 500, -600.0, 600.0, kTPCQvecX+ih);
 720         AddHistogram(classStr.Data(), Form("TPCQvecY_h%d_CentSPDVtxZ_prof", ih+1),
 721                      Form("TPC Q_{y}, harmonic %d, vs centSPD and vtxZ; Centrality (percents); VtxZ (cm); <Q_{y}>", ih+1), kTRUE,
 722                      20, 0.0, 100.0, kCentSPD, 24, -12.0, 12.0, kVtxZ, 500, -600.0, 600.0, kTPCQvecY+ih);
 723       }    // end loop over harmonics
 724       continue;
 725     }    // end if for the TPC reaction plane histograms
 726
 727     TString trkFlagNames = "";
 728     for(Int_t iflag=0; iflag<kNTrackingFlags; ++iflag) {
 729       trkFlagNames += gkTrackingFlagNames[iflag];
 730       trkFlagNames += ";";
 731     }
 732
 733     // Track histograms
 734     if(classStr.Contains("TrackingFlags")) {
 735       AddHistClass(classStr.Data());
 736
 737       AddHistogram(classStr.Data(), "TrackingFlags", "Tracking flags;;", kFALSE,
 738                    kNTrackingFlags, -0.5, kNTrackingFlags-0.5, kTrackingFlag, 0, 0.0, 0.0, kNothing, 0, 0.0, 0.0, kNothing, trkFlagNames.Data());
 739       AddHistogram(classStr.Data(), "TrackingFlags_Pt", "Tracking flags vs p_{T};p_{T} (GeV/c);", kFALSE,
 740                    100, 0.0, 20.0, kPt, kNTrackingFlags, -0.5, kNTrackingFlags-0.5, kTrackingFlag, 0, 0.0, 0.0, kNothing, "", trkFlagNames.Data());
 741       AddHistogram(classStr.Data(), "TrackingFlags_Eta", "Tracking flags vs #eta;#eta;", kFALSE,
 742                    30, -1.5, 1.5, kEta, kNTrackingFlags, -0.5, kNTrackingFlags-0.5, kTrackingFlag, 0, 0.0, 0.0, kNothing, "", trkFlagNames.Data());
 743       AddHistogram(classStr.Data(), "TrackingFlags_Phi", "Tracking flags vs #varphi;#varphi (rad.);", kFALSE,
 744                    60, 0.0, 6.3, kPhi, kNTrackingFlags, -0.5, kNTrackingFlags-0.5, kTrackingFlag, 0, 0.0, 0.0, kNothing, "", trkFlagNames.Data());
 745       AddHistogram(classStr.Data(), "TrackingFlags_CentVZERO", "Tracking flags vs centrality VZERO; centrality (%);", kFALSE,
 746                    20, 0.0, 100.0, kCentVZERO, kNTrackingFlags, -0.5, kNTrackingFlags-0.5, kTrackingFlag, 0, 0.0, 0.0, kNothing, "", trkFlagNames.Data());
 747
 748       AddHistogram(classStr.Data(), "TrackingFlags_TRDntracklets", "Tracking flags vs TRD #tracklets; TRD #tracklets;", kFALSE,
 749                    7, -0.5, 6.5, kTRDntracklets, kNTrackingFlags, -0.5, kNTrackingFlags-0.5, kTrackingFlag, 0, 0.0, 0.0, kNothing, "", trkFlagNames.Data());
 750       AddHistogram(classStr.Data(), "TrackingFlags_TRDntrackletsPID", "Tracking flags vs TRD # pid tracklets; TRD #tracklets;", kFALSE,
 751                    7, -0.5, 6.5, kTRDntrackletsPID, kNTrackingFlags, -0.5, kNTrackingFlags-0.5, kTrackingFlag, 0, 0.0, 0.0, kNothing, "", trkFlagNames.Data());
 752       AddHistogram(classStr.Data(), "TrackingFlags_TRDeleProbab", "Tracking flags vs TRD electron probability; TRD probab.;", kFALSE,
 753                    50, 0.0, 1.0, kTRDpidProbabilities, kNTrackingFlags, -0.5, kNTrackingFlags-0.5, kTrackingFlag, 0, 0.0, 0.0, kNothing, "", trkFlagNames.Data());
 754       AddHistogram(classStr.Data(), "TrackingFlags_TOFbeta", "Tracking flags vs TOF #beta; TOF #beta;", kFALSE,
 755                    50, 0.0, 1.0, kTOFbeta, kNTrackingFlags, -0.5, kNTrackingFlags-0.5, kTrackingFlag, 0, 0.0, 0.0, kNothing, "", trkFlagNames.Data());
 756     }
 757
 758     if(classStr.Contains("TrackQA")) {
 759       AddHistClass(classStr.Data());
 760
 761       AddHistogram(classStr.Data(), "Pt", "p_{T} distribution; p_{T} (GeV/c^{2});", kFALSE,
 762                    1000, 0.0, 50.0, kPt);
 763       AddHistogram(classStr.Data(), "Eta", "#eta illumination; #eta;", kFALSE,
 764                    1000, -1.5, 1.5, kEta);
 765       AddHistogram(classStr.Data(), "Phi", "#varphi illumination; #varphi;", kFALSE,
 766                    1000, 0.0, 6.3, kPhi);
 767       AddHistogram(classStr.Data(), "DCAxy", "DCAxy; DCAxy (cm.)", kFALSE,
 768                    1000, -10.0, 10.0, kDcaXY);
 769       AddHistogram(classStr.Data(), "DCAz", "DCAz; DCAz (cm.)", kFALSE,
 770                    1000, -10.0, 10.0, kDcaZ);
 771
 772       // run dependence
 773       AddHistogram(classStr.Data(), "Pt_Run", "<p_{T}> vs run; run;", kTRUE,
 774                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, 0.0, 50.0, kPt);
 775       AddHistogram(classStr.Data(), "Eta_Run", "<#eta> vs run; run;", kTRUE,
 776                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, -1.5, 1.5, kEta);
 777       AddHistogram(classStr.Data(), "Phi_Run", "<#varphi> vs run; run;", kTRUE,
 778                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, 0.0, 6.3, kPhi);
 779       AddHistogram(classStr.Data(), "DCAxy_Run", "<DCAxy> vs run; run;", kTRUE,
 780                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, -10.0, 10.0, kDcaXY);
 781       AddHistogram(classStr.Data(), "DCAz_Run", "<DCAz> vs run; run;", kTRUE,
 782                    kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, -10.0, 10.0, kDcaZ);
 783
 784       // correlations between parameters
 785       AddHistogram(classStr.Data(), "Eta_Pt_prof", "<p_{T}> vs #eta; #eta; p_{T} (GeV/c);", kTRUE,
 786                    300, -1.5, +1.5, kEta, 100, 0.0, 10.0, kPt);
 787       AddHistogram(classStr.Data(), "Phi_Pt_prof", "<p_{T}> vs #varphi; #varphi (rad.); p_{T} (GeV/c)", kTRUE,
 788                    300, 0.0, 6.3, kPhi, 100, 0.0, 10.0, kPt);
 789       AddHistogram(classStr.Data(), "Eta_Phi", "#varphi vs #eta; #eta; #varphi (rad.);", kFALSE,
 790                    200, -1.0, +1.0, kEta, 100, 0.0, 6.3, kPhi);
 791       AddHistogram(classStr.Data(), "Pt_DCAxy", "DCAxy vs p_{T}; p_{T} (GeV/c); DCA_{xy} (cm)", kFALSE,
 792                    100, 0.0, 10.0, kPt, 500, -2.0, 2.0, kDcaXY);
 793       AddHistogram(classStr.Data(), "Pt_DCAz", "DCAz vs p_{T}; p_{T} (GeV/c); DCA_{z} (cm)", kFALSE,
 794                    100, 0.0, 10.0, kPt, 500, -2.0, 2.0, kDcaZ);
 795       AddHistogram(classStr.Data(), "Eta_DCAxy", "DCAxy vs #eta; #eta; DCA_{xy} (cm)", kFALSE,
 796                    100, -1.0, 1.0, kEta, 500, -2.0, 2.0, kDcaXY);
 797       AddHistogram(classStr.Data(), "Eta_DCAz", "DCAz vs #eta; #eta; DCA_{z} (cm)", kFALSE,
 798                    100, -1.0, 1.0, kEta, 500, -2.0, 2.0, kDcaZ);
 799
 800       if(classStr.Contains("ITS")) {
 801         AddHistogram(classStr.Data(),"ITSncls", "ITS nclusters;# clusters ITS", kFALSE,
 802                      7,-0.5,6.5,kITSncls);
 803         AddHistogram(classStr.Data(),"ITSncls_Run", "ITS <nclusters> vs run;run;# clusters ITS", kTRUE,
 804                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 7,-0.5,6.5, kITSncls);
 805         AddHistogram(classStr.Data(),"ITSncls_Cent_prof", "ITS <nclusters> vs centrality; centrality;# clusters ITS", kTRUE,
 806                      20, 0.0, 100.0, kCentVZERO, 7,-0.5,6.5, kITSncls);
 807         AddHistogram(classStr.Data(),"Pt_ITSncls","ITS nclusters vs p_{T};p_{T} (GeV/c);# clusters ITS",kFALSE,
 808                      100, 0.0, 10.0, kPt, 7,-0.5,6.5, kITSncls);
 809         AddHistogram(classStr.Data(),"Eta_Phi_ITSncls_prof","ITS <nclusters> vs (#eta,#phi);#eta;#phi (rad.);# clusters ITS",kTRUE,
 810                      192, -1.2, 1.2, kEta, 126, 0.0, 6.3, kPhi, 7, -0.5, 6.5, kITSncls);
 811       }  // end if ITS histograms
 812
 813       if(classStr.Contains("TPC")) {
 814         AddHistogram(classStr.Data(),"TPCncls","TPC nclusters;# clusters TPC",kFALSE,
 815                      160,-0.5,159.5,kTPCncls);
 816         AddHistogram(classStr.Data(),"TPCncls_Run","TPC nclusters vs run;run;# clusters TPC",kTRUE,
 817                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 160,-0.5,159.5, kTPCncls);
 818         AddHistogram(classStr.Data(),"TPCncls_CentVZERO","TPC nclusters vs centrality;centrality;# clusters TPC",kFALSE,
 819                      20, 0.0, 100.0, kCentVZERO, 160,-0.5,159.5, kTPCncls);
 820         AddHistogram(classStr.Data(),"TPCncls_Pin","TPC nclusters vs inner param P; P (GeV/c); # clusters TPC",kFALSE,
 821                      200,0.0,20,kPin,160,-0.5,159.5,kTPCncls);
 822         AddHistogram(classStr.Data(),"Eta_Phi_TPCncls_prof","TPC <nclusters> vs (#eta,#phi);#eta;#phi (rad.);# clusters TPC",kTRUE,
 823                      192, -1.2, 1.2, kEta, 126, 0.0, 6.3, kPhi, 160, -0.5, 159.5, kTPCncls);
 824         AddHistogram(classStr.Data(),"TPCsignal_Pin","TPC dE/dx vs. inner param P;P (GeV/c);TPC dE/dx",kFALSE,
 825                      1000,0.0,20.0,kPin,151,-0.5,150.5,kTPCsignal);
 826         AddHistogram(classStr.Data(),"TPCsignal_TPCclusters","TPC dE/dx vs. TPC nclusters;TPC #clusters;TPC dE/dx",kFALSE,
 827                      160,-0.5,159.5,kTPCncls,151,-0.5,150.5,kTPCsignal);
 828         AddHistogram(classStr.Data(),"TPCnsigElectron_Pin","TPC N_{#sigma} electron vs. inner param P;P (GeV/c);N_{#sigma}",kFALSE,
 829                      1000,0.0,20.0,kPin,100,-5.0,5.0,kTPCnSig+kElectron);
 830         AddHistogram(classStr.Data(),"TPCnsigElectron_CentEta_prof","TPC N_{#sigma} electron vs. (#eta,centrality); #eta; centrality VZERO;N_{#sigma}",kTRUE,
 831                      100,-1.0,1.0,kEta, 20, 0.0, 100.0, kCentVZERO, 100,-5.0,5.0,kTPCnSig+kElectron);
 832         AddHistogram(classStr.Data(),"TPCnsigElectron_Run","TPC N_{#sigma} electron vs. run; run; N_{#sigma}",kTRUE,
 833                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo,100,-5.0,5.0,kTPCnSig+kElectron);
 834       }      // end if TPC histograms
 835
 836       if(classStr.Contains("TRD")) {
 837         AddHistogram(classStr.Data(),"TRDntracklets","TRD ntracklets; #tracklets TRD",kFALSE,
 838                      7,-0.5,6.5,kTRDntracklets);
 839         AddHistogram(classStr.Data(),"TRDntrackletsPID","TRD ntracklets PID; #tracklets TRD",kFALSE,
 840                      7,-0.5,6.5,kTRDntrackletsPID);
 841         AddHistogram(classStr.Data(),"TRDprobabElectron","TRD electron probability; probability",kFALSE,
 842                      500,0.0,1.0,kTRDpidProbabilities);
 843         AddHistogram(classStr.Data(),"TRDprobabPion","TRD pion probability; probability",kFALSE,
 844                      500,0.0,1.0,kTRDpidProbabilities+1);
 845         AddHistogram(classStr.Data(),"TRDntracklets_Run","TRD <ntracklets> vs run; run; #tracklets TRD",kTRUE,
 846                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 7,-0.5,6.5,kTRDntracklets);
 847         AddHistogram(classStr.Data(),"TRDntrackletsPID_Run","TRD <ntracklets PID> vs run; run; #tracklets TRD",kTRUE,
 848                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 7,-0.5,6.5,kTRDntrackletsPID);
 849         AddHistogram(classStr.Data(),"TRDprobabEle_Run","TRD <electron probability> vs run; run; probability",kTRUE,
 850                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 10,0.0,1.0,kTRDpidProbabilities);
 851         AddHistogram(classStr.Data(),"TRDprobabPion_Run","TRD <pion probability> vs run; run; probability",kTRUE,
 852                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 10,0.0,1.0,kTRDpidProbabilities+1);
 853         AddHistogram(classStr.Data(),"Eta_Phi_TRDntracklets_prof","TRD <ntracklets> vs (#eta,#phi);#eta;#phi (rad.);#tracklets TRD",kTRUE,
 854                      192, -1.2, 1.2, kEta, 126, 0.0, 6.3, kPhi, 7, -0.5, 6.5, kTRDntracklets);
 855         AddHistogram(classStr.Data(),"TRDntracklets_cent","TRD ntracklets vs centrality; #tracklets TRD; centrality (%)",kFALSE,
 856                      7,-0.5,6.5,kTRDntracklets,20,0.0, 100.0, kCentSPD);
 857         AddHistogram(classStr.Data(),"Eta_Phi_TRDntrackletsPID_prof","TRD <ntracklets PID> vs (#eta,#phi);#eta;#phi (rad.);#tracklets TRD",kTRUE,
 858                      192, -1.2, 1.2, kEta, 126, 0.0, 6.3, kPhi, 7, -0.5, 6.5, kTRDntrackletsPID);
 859         AddHistogram(classStr.Data(),"TRDntrackletsPID_cent","TRD ntracklets PID vs centrality; #tracklets TRD; centrality (%)",kFALSE,
 860                      7,-0.5,6.5,kTRDntrackletsPID,20,0.0, 100.0, kCentSPD);
 861         AddHistogram(classStr.Data(),"TRDntracklets_TRDntrackletsPID","TRD ntracklets vs TRD ntracklets PID; #tracklets TRD; #tracklets TRD pid",kFALSE,
 862                      7,-0.5,6.5,kTRDntracklets,7,-0.5, 6.5, kTRDntrackletsPID);
 863       }   // end if TRD histograms
 864
 865       if(classStr.Contains("TOF")) {
 866         AddHistogram(classStr.Data(),"TOFbeta_P","TOF #beta vs P;P (GeV/c);#beta",kFALSE,
 867                      200,0.0,20.0,kP, 220,0.0,1.1,kTOFbeta);
 868         AddHistogram(classStr.Data(),"Eta_Phi_TOFbeta_prof","TOF <#beta> vs (#eta,#phi);#eta;#phi (rad.);TOF #beta",kTRUE,
 869                      192, -1.2, 1.2, kEta, 126, 0.0, 6.3, kPhi, 160, -0.5, 159.5, kTOFbeta);
 870         AddHistogram(classStr.Data(),"TOFnsigElectron_P","TOF N_{#sigma} electron vs. P;P (GeV/c);N_{#sigma}",kFALSE,
 871                      200,0.0,20.0,kP, 100,-5.0,5.0,kTOFnSig+kElectron);
 872         AddHistogram(classStr.Data(),"TOFnSigElectron_Run","TOF <n-#sigma_{e}> vs run number;run;TOF n-#sigma_{e}",kTRUE,
 873                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 160, -0.5, 159.5, kTOFnSig+kElectron);
 874       }    // end if TOF histograms
 875       continue;
 876     }  // end if "TrackQA"
 877
 878     Double_t massBinWidth = 0.04;     // *GeV/c^2
 879     Double_t massRange[2] = {0.0,6.0};
 880     const Int_t nMassBins = TMath::Nint((massRange[1]-massRange[0])/massBinWidth);
 881     Double_t massBinLims[nMassBins+1]; for(Int_t im=0; im<=nMassBins; ++im) massBinLims[im] = massBinWidth*im;
 882
 883     // Histograms for pairs
 884     if(classStr.Contains("PairQA")) {
 885       AddHistClass(classStr.Data());
 886
 887       AddHistogram(classStr.Data(), "Mass_CentVZEROCep",
 888                    "Inv. mass vs (centVZERO, #Psi^{2}); centrality VZERO; #Psi^{2}; m (GeV/c^{2})",
 889                    kFALSE, gkNCentRanges, gkEMCentLims, kCentVZERO, gkNPhiEPRanges, gkEMPhiEPLims, gkEPused, nMassBins, massBinLims, kMass);
 890       AddHistogram(classStr.Data(), "Pt", "Pt; p_{T} (GeV/c)", kFALSE,
 891                      1000, 0.0, 10.0, kPairPt);
 892
 893       if(classStr.Contains("SE")) {
 894         AddHistogram(classStr.Data(), "CandidateId", "Candidate id; p_{T} (GeV/c)", kFALSE,
 895                      AliCorrelationReducedPair::kNMaxCandidateTypes+1, -0.5, Double_t(AliCorrelationReducedPair::kNMaxCandidateTypes)+0.5, kCandidateId);
 896         AddHistogram(classStr.Data(), "PairType", "Pair type; pair type", kFALSE,
 897                      4, -0.5, 3.5, kPairType);
 898         AddHistogram(classStr.Data(), "Mass", "Invariant mass; m_{inv} (GeV/c^{2})", kFALSE,
 899                      nMassBins, massRange[0], massRange[1], kMass);
 900         AddHistogram(classStr.Data(), "Rapidity", "Rapidity; y", kFALSE,
 901                      240, -1.2, 1.2, kPairRap);
 902         AddHistogram(classStr.Data(), "Eta", "Pseudo-rapidity #eta; #eta", kFALSE,
 903                      240, -2.0, 2.0, kPairEta);
 904         AddHistogram(classStr.Data(), "Phi", "Azimuthal distribution; #phi (rad.)", kFALSE,
 905                      315, 0.0, 6.3, kPairPhi);
 906         AddHistogram(classStr.Data(), "OpeningAngle", "Opening angle; op. angle (rad)", kFALSE,
 907                      1000, 0.0, 3.2, kPairOpeningAngle);
 908         AddHistogram(classStr.Data(), "Mass_Pt", "Pt vs invariant mass; m_{inv} (GeV/c^{2}); p_{T} (GeV/c)", kFALSE,
 909                      nMassBins, massRange[0], massRange[1], kMass, 100, 0.0, 10.0, kPairPt);
 910         AddHistogram(classStr.Data(), "Mass_Run_prof", "<Invariant mass> vs run; run;m_{inv} (GeV/c^{2})", kTRUE,
 911                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo, nMassBins, massRange[0], massRange[1], kMass);
 912         AddHistogram(classStr.Data(), "Pt_Run_prof", "<p_{T}> vs run; run;p_{T} (GeV/c)", kTRUE,
 913                      kNRunBins, runHistRange[0], runHistRange[1], kRunNo, 1000, 0.0, 10.0, kPairPt);
 914         AddHistogram(classStr.Data(), "VZEROflow_sideA_v2_Mass",
 915                      "Pair v2 coefficient using VZERO-A reaction plane; inv.mass (GeV/c^{2});v_{2}(#Psi_{2}^{VZERO-A})", kTRUE,
 916                      nMassBins, massRange[0], massRange[1], kMass, 1000, -1., 1., kPairVZEROFlowVn+0*6+1);
 917         AddHistogram(classStr.Data(), "VZEROflow_sideC_v2_Mass",
 918                      "Pair v2 coefficient using VZERO-C reaction plane; inv.mass (GeV/c^{2});v_{2}(#Psi_{2}^{VZERO-C})", kTRUE,
 919                      nMassBins, massRange[0], massRange[1], kMass, 1000, -1., 1., kPairVZEROFlowVn+1*6+1);
 920         AddHistogram(classStr.Data(), "TPCflow_v2_Mass",
 921                      "Pair v2 coefficient using TPC reaction plane; inv.mass (GeV/c^{2});v_{2}(#Psi_{2}^{TPC})", kTRUE,
 922                      nMassBins, massRange[0], massRange[1], kMass, 1000, -1., 1., kPairTPCFlowVn+1);
 923       }
 924       continue;
 925     }   // end if for Pair classes of histograms
 926
 927   }   // end loop over histogram classes
 928 }