PWG0/dNdPt/AlidNdPtCutAnalysisPbPb.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 //------------------------------------------------------------------------------\r
  16 // AlidNdPtCutAnalysisPbPb class. \r
  17 //\r
  18 // a. functionality:\r
  19 // - fills generic cut histograms\r
  20 // - generates cuts (selection criteria)\r
  21 //\r
  22 // b. data members:\r
  23 // - generic cut histograms\r
  24 // - control histograms\r
  25 //\r
  26 // Author: J.Otwinowski 04/11/2008 \r
  27 //------------------------------------------------------------------------------\r
  28 #include "TH1.h"\r
  29 #include "TH2.h"\r
  30 \r
  31 #include "AliHeader.h"  \r
  32 #include "AliGenEventHeader.h"  \r
  33 #include "AliStack.h"  \r
  34 #include "AliESDEvent.h"  \r
  35 #include "AliMCEvent.h"  \r
  36 #include "AliESDtrackCuts.h"  \r
  37 #include "AliLog.h" \r
  38 #include "AliTracker.h" \r
  39 \r
  40 #include "AlidNdPtEventCuts.h"\r
  41 #include "AlidNdPtAcceptanceCuts.h"\r
  42 #include "AlidNdPtBackgroundCuts.h"\r
  43 #include "AlidNdPtAnalysis.h"\r
  44 #include "AliPhysicsSelection.h"\r
  45 \r
  46 #include "AliPWG0Helper.h"\r
  47 #include "AlidNdPtHelper.h"\r
  48 #include "AlidNdPtCutAnalysisPbPb.h"\r
  49 \r
  50 using namespace std;\r
  51 \r
  52 ClassImp(AlidNdPtCutAnalysisPbPb)\r
  53 \r
  54 //_____________________________________________________________________________\r
  55   AlidNdPtCutAnalysisPbPb::AlidNdPtCutAnalysisPbPb(): AlidNdPt(),\r
  56   fAnalysisFolder(0),\r
  57   fEventCount(0),\r
  58   fRecEventHist(0),\r
  59   fMCEventHist(0),\r
  60   fRecMCEventHist(0),\r
  61   fRecMCTrackHist(0)\r
  62 {\r
  63   // default constructor\r
  64   Init();\r
  65 }\r
  66 \r
  67 //_____________________________________________________________________________\r
  68 AlidNdPtCutAnalysisPbPb::AlidNdPtCutAnalysisPbPb(Char_t* name, Char_t* title): AlidNdPt(name,title),\r
  69   fAnalysisFolder(0),\r
  70   fEventCount(0),\r
  71   fRecEventHist(0),\r
  72   fMCEventHist(0),\r
  73   fRecMCEventHist(0),\r
  74   fRecMCTrackHist(0)\r
  75 {\r
  76   // constructor\r
  77   Init();\r
  78 }\r
  79 \r
  80 //_____________________________________________________________________________\r
  81 AlidNdPtCutAnalysisPbPb::~AlidNdPtCutAnalysisPbPb() {\r
  82   // \r
  83   if(fEventCount) delete fEventCount; fEventCount=0;\r
  84   if(fRecEventHist) delete fRecEventHist; fRecEventHist=0;\r
  85   if(fMCEventHist) delete fMCEventHist; fMCEventHist=0;\r
  86   if(fRecMCEventHist) delete fRecMCEventHist; fRecMCEventHist=0;\r
  87   if(fRecMCTrackHist) delete fRecMCTrackHist; fRecMCTrackHist=0;\r
  88 \r
  89   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;\r
  90 }\r
  91 \r
  92 //_____________________________________________________________________________\r
  93 void AlidNdPtCutAnalysisPbPb::Init(){\r
  94   //\r
  95   // Init histograms\r
  96   //\r
  97   const Int_t ptNbins = 58; \r
  98   const Double_t ptMin = 0.; \r
  99   const Double_t ptMax = 20.; \r
 100 \r
 101   Double_t binsPt[ptNbins+1] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8,0.85,0.9,0.95,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.2,2.4,2.6,2.8,3.0,3.2,3.4,3.6,3.8,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,9.0,10.0,11.0,12.0,13.0,14.0,15.0,16.0,18.0, 20.};\r
 102 \r
 103   // \r
 104   Int_t binsEventCount[2]={2,2};\r
 105   Double_t minEventCount[2]={0,0}; \r
 106   Double_t maxEventCount[2]={2,2}; \r
 107   fEventCount = new THnSparseF("fEventCount","trig vs trig+vertex",2,binsEventCount,minEventCount,maxEventCount);\r
 108   fEventCount->GetAxis(0)->SetTitle("trig");\r
 109   fEventCount->GetAxis(1)->SetTitle("trig+vert");\r
 110   fEventCount->Sumw2();\r
 111 \r
 112   //Xv:Yv:Zv:ResZv:Mult\r
 113   Double_t kFact = 1.0;\r
 114 \r
 115   Int_t binsRecEventHist[5]={80,80,100,80,150};\r
 116   Double_t minRecEventHist[5]={-3.*kFact,-3.*kFact,-35.,0.,0.}; \r
 117   Double_t maxRecEventHist[5]={3.*kFact,3.*kFact,35.,10.,150.}; \r
 118   fRecEventHist = new THnSparseF("fRecEventHist","Xv:Yv:Zv:ResZv:Mult",5,binsRecEventHist,minRecEventHist,maxRecEventHist);\r
 119   fRecEventHist->GetAxis(0)->SetTitle("Xv (cm)");\r
 120   fRecEventHist->GetAxis(1)->SetTitle("Yv (cm)");\r
 121   fRecEventHist->GetAxis(2)->SetTitle("Zv (cm)");\r
 122   fRecEventHist->GetAxis(3)->SetTitle("ResZv (cm)");\r
 123   fRecEventHist->GetAxis(4)->SetTitle("Mult");\r
 124   fRecEventHist->Sumw2();\r
 125 \r
 126   //Xv:Yv:Zv\r
 127   Int_t binsMCEventHist[3]={80,80,100};\r
 128   Double_t minMCEventHist[3]={-0.1,-0.1,-35.}; \r
 129   Double_t maxMCEventHist[3]={0.1,0.1,35.}; \r
 130   fMCEventHist = new THnSparseF("fMCEventHist","mcXv:mcYv:mcZv",3,binsMCEventHist,minMCEventHist,maxMCEventHist);\r
 131   fMCEventHist->GetAxis(0)->SetTitle("mcXv (cm)");\r
 132   fMCEventHist->GetAxis(1)->SetTitle("mcYv (cm)");\r
 133   fMCEventHist->GetAxis(2)->SetTitle("mcZv (cm)");\r
 134   fMCEventHist->Sumw2();\r
 135 \r
 136   //Xv-mcXv:Yv-mcYv:Zv-mcZv:Mult\r
 137   Int_t binsRecMCEventHist[4]={100,100,100,150};\r
 138   Double_t minRecMCEventHist[4]={-1.0*kFact,-1.0*kFact,-1.0*kFact,0.}; \r
 139   Double_t maxRecMCEventHist[4]={1.0*kFact,1.0*kFact,1.0*kFact,150.}; \r
 140   fRecMCEventHist = new THnSparseF("fRecMCEventHist","Xv-mcXv:Yv-mcYv:Zv-mcZv:Mult",4,binsRecMCEventHist,minRecMCEventHist,maxRecMCEventHist);\r
 141   fRecMCEventHist->GetAxis(0)->SetTitle("Xv-mcXv (cm)");\r
 142   fRecMCEventHist->GetAxis(1)->SetTitle("Yv-mcYv (cm)");\r
 143   fRecMCEventHist->GetAxis(2)->SetTitle("Zv-mcZv (cm)");\r
 144   fRecMCEventHist->GetAxis(3)->SetTitle("Mult");\r
 145   fRecMCEventHist->Sumw2();\r
 146 \r
 147   //\r
 148   // THnSparse track histograms\r
 149   //\r
 150   //nClust:chi2PerClust:nClust/nFindableClust:DCAy:DCAz:eta:phi:pt:hasStrangeMother:isFromMaterial:isPrim:charge\r
 151   Int_t binsRecMCTrackHist[12]={160,20,40,100,100,30,90,ptNbins, 2, 2, 2, 3};\r
 152   //Double_t minRecMCTrackHist[12]={0., 0., 0., -10.,-10.,-1.5, 0., ptMin, 0., 0., 0., -1.};\r
 153   //Double_t maxRecMCTrackHist[12]={160.,10.,1.2, 10.,10.,1.5, 2.*TMath::Pi(), ptMax, 2., 2., 2., 2.};\r
 154   Double_t minRecMCTrackHist[12]={0., 0., 0., -1.,-1.,-1.5, 0., ptMin, 0., 0., 0., -1.};\r
 155   Double_t maxRecMCTrackHist[12]={160.,10.,1.2, 1.,1.,1.5, 2.*TMath::Pi(), ptMax, 2., 2., 2., 2.};\r
 156 \r
 157   fRecMCTrackHist = new THnSparseF("fRecMCTrackHist","nClust:chi2PerClust:nClust/nFindableClust:DCAy:DCAz:eta:phi:pt:hasStrangeMother:isFromMaterial:isPrim:charge",12,binsRecMCTrackHist,minRecMCTrackHist,maxRecMCTrackHist);\r
 158   fRecMCTrackHist->SetBinEdges(7,binsPt);\r
 159 \r
 160   fRecMCTrackHist->GetAxis(0)->SetTitle("nClust");\r
 161   fRecMCTrackHist->GetAxis(1)->SetTitle("chi2PerClust");\r
 162   fRecMCTrackHist->GetAxis(2)->SetTitle("nClust/nFindableClust");\r
 163   fRecMCTrackHist->GetAxis(3)->SetTitle("DCAy (cm)");\r
 164   fRecMCTrackHist->GetAxis(4)->SetTitle("DCAz (cm)");\r
 165   fRecMCTrackHist->GetAxis(5)->SetTitle("#eta");\r
 166   fRecMCTrackHist->GetAxis(6)->SetTitle("#phi (rad)");\r
 167   fRecMCTrackHist->GetAxis(7)->SetTitle("p_{T} (GeV/c)");\r
 168   fRecMCTrackHist->GetAxis(8)->SetTitle("hasStrangeMother");\r
 169   fRecMCTrackHist->GetAxis(9)->SetTitle("isFromMaterial");\r
 170   fRecMCTrackHist->GetAxis(10)->SetTitle("isPrim");\r
 171   fRecMCTrackHist->GetAxis(11)->SetTitle("charge");\r
 172   fRecMCTrackHist->Sumw2();\r
 173 \r
 174   // init output folder\r
 175   fAnalysisFolder = CreateFolder("folderdNdPt","Analysis dNdPt Folder");\r
 176 \r
 177 }\r
 178 \r
 179 //_____________________________________________________________________________\r
 180 void AlidNdPtCutAnalysisPbPb::Process(AliESDEvent *const esdEvent, AliMCEvent * const mcEvent)\r
 181 {\r
 182   //\r
 183   // Process real and/or simulated events\r
 184   //\r
 185   if(!esdEvent) {\r
 186     AliDebug(AliLog::kError, "esdEvent not available");\r
 187     return;\r
 188   }\r
 189 \r
 190   // get selection cuts\r
 191   AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
 192   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
 193   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
 194 \r
 195   if(!evtCuts || !accCuts  || !esdTrackCuts) {\r
 196     AliDebug(AliLog::kError, "cuts not available");\r
 197     return;\r
 198   }\r
 199 \r
 200   // trigger selection\r
 201   Bool_t isEventTriggered = kTRUE;\r
 202   AliPhysicsSelection *trigSel = NULL;\r
 203   AliTriggerAnalysis *trigAna = NULL; // needed for andV0\r
 204 \r
 205   if(evtCuts->IsTriggerRequired())  \r
 206   {\r
 207     //\r
 208     trigSel = GetPhysicsTriggerSelection();\r
 209     if(!trigSel) {\r
 210       printf("cannot get trigSel \n");\r
 211       return;\r
 212     }\r
 213 \r
 214     //\r
 215     if(IsUseMCInfo()) \r
 216     { \r
 217       trigSel->SetAnalyzeMC();\r
 218 \r
 219       \r
 220         isEventTriggered = trigSel->IsCollisionCandidate(esdEvent);\r
 221         \r
 222         if(GetTrigger() == AliTriggerAnalysis::kV0AND) \r
 223         {\r
 224           trigAna = trigSel->GetTriggerAnalysis();\r
 225           if(!trigAna) \r
 226             return;\r
 227 \r
 228           isEventTriggered = trigAna->IsOfflineTriggerFired(esdEvent, GetTrigger());\r
 229         }//if(GetTrigger() == AliTriggerAnalysis::kV0AND)\r
 230      }//if(IsUseMCInfo())\r
 231   }//if(evtCuts->IsTriggerRequired())\r
 232 \r
 233   // use MC information\r
 234   AliHeader* header = 0;\r
 235   AliGenEventHeader* genHeader = 0;\r
 236   AliStack* stack = 0;\r
 237   TArrayF vtxMC(3);\r
 238 \r
 239   Int_t multMCTrueTracks = 0;\r
 240   if(IsUseMCInfo())\r
 241   {\r
 242     //\r
 243     if(!mcEvent) {\r
 244       AliDebug(AliLog::kError, "mcEvent not available");\r
 245       return;\r
 246     }\r
 247     // get MC event header\r
 248     header = mcEvent->Header();\r
 249     if (!header) {\r
 250       AliDebug(AliLog::kError, "Header not available");\r
 251       return;\r
 252     }\r
 253     // MC particle stack\r
 254     stack = mcEvent->Stack();\r
 255     if (!stack) {\r
 256       AliDebug(AliLog::kError, "Stack not available");\r
 257       return;\r
 258     }\r
 259 \r
 260     // get MC vertex\r
 261     genHeader = header->GenEventHeader();\r
 262     if (!genHeader) {\r
 263       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");\r
 264       return;\r
 265     }\r
 266     genHeader->PrimaryVertex(vtxMC);\r
 267 \r
 268     Double_t vMCEventHist[3]={vtxMC[0],vtxMC[1],vtxMC[2]};\r
 269     fMCEventHist->Fill(vMCEventHist);\r
 270 \r
 271     // multipliticy of all MC primary tracks\r
 272     // in Zv, pt and eta ranges)\r
 273     multMCTrueTracks = AlidNdPtHelper::GetMCTrueTrackMult(mcEvent,evtCuts,accCuts);\r
 274 \r
 275   } // end bUseMC\r
 276 \r
 277   // get reconstructed vertex  \r
 278   const AliESDVertex* vtxESD = 0; \r
 279   if(evtCuts->IsRecVertexRequired()) \r
 280   {\r
 281      if(GetAnalysisMode() == AlidNdPtHelper::kTPC) {\r
 282         vtxESD = esdEvent->GetPrimaryVertexTPC();\r
 283     }\r
 284     else if(GetAnalysisMode() == AlidNdPtHelper::kTPCITS) {\r
 285       vtxESD = esdEvent->GetPrimaryVertexTracks();\r
 286     }\r
 287     else {\r
 288         return;\r
 289     }\r
 290   }\r
 291   Bool_t isEventOK = evtCuts->AcceptEvent(esdEvent,mcEvent,vtxESD); \r
 292   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK, isEventTriggered);\r
 293   //printf("GetAnalysisMode() %d \n",GetAnalysisMode());\r
 294 \r
 295   TObjArray *allChargedTracks=0;\r
 296   Int_t multAll=0;\r
 297   \r
 298   //\r
 299   // event counter\r
 300   // \r
 301   //printf("isEventOK %d, isEventTriggered %d \n",isEventOK,isEventTriggered);\r
 302 \r
 303   Bool_t isTrigAndVertex = isEventTriggered && isEventOK;\r
 304   Double_t vEventCount[2] = { isEventTriggered, isTrigAndVertex};\r
 305   fEventCount->Fill(vEventCount);\r
 306 \r
 307 \r
 308   // check event cuts\r
 309   if(isEventOK && isEventTriggered)\r
 310   {\r
 311     // get all charged tracks\r
 312     allChargedTracks = AlidNdPtHelper::GetAllChargedTracks(esdEvent,GetAnalysisMode());\r
 313     if(!allChargedTracks) return;\r
 314 \r
 315     Int_t entries = allChargedTracks->GetEntries();\r
 316     for(Int_t i=0; i<entries;++i) \r
 317     {\r
 318       AliESDtrack *track = (AliESDtrack*)allChargedTracks->At(i);\r
 319       if(!track) continue;\r
 320 \r
 321       if(!esdTrackCuts->AcceptTrack(track)) continue;\r
 322       FillHistograms(track, stack);\r
 323       multAll++;\r
 324     }\r
 325 \r
 326     Double_t vRecEventHist[5] = {vtxESD->GetXv(),vtxESD->GetYv(),vtxESD->GetZv(),vtxESD->GetZRes(),multAll};\r
 327     fRecEventHist->Fill(vRecEventHist);\r
 328 \r
 329     if(IsUseMCInfo()) {\r
 330       Double_t vRecMCEventHist[5] = {vtxESD->GetXv()-vtxMC[0],vtxESD->GetYv()-vtxMC[1],vtxESD->GetZv()-vtxMC[2],multAll};\r
 331       fRecMCEventHist->Fill(vRecMCEventHist);\r
 332     }\r
 333   }\r
 334 \r
 335   if(allChargedTracks) delete allChargedTracks; allChargedTracks = 0;\r
 336 \r
 337 }\r
 338 \r
 339 //_____________________________________________________________________________\r
 340 void AlidNdPtCutAnalysisPbPb::FillHistograms(AliESDtrack *const esdTrack, AliStack *const stack) const\r
 341 {\r
 342   //\r
 343   // Fill ESD track and MC histograms \r
 344   //\r
 345   if(!esdTrack) return;\r
 346   if(esdTrack->Charge() == 0.) return;\r
 347 \r
 348   Float_t pt = esdTrack->Pt();\r
 349   Float_t eta = esdTrack->Eta();\r
 350   Float_t phi = esdTrack->Phi();\r
 351 \r
 352   Int_t nClust = 0;\r
 353   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) { \r
 354     nClust = esdTrack->GetTPCNclsIter1();\r
 355   } else {\r
 356     nClust = esdTrack->GetTPCclusters(0);\r
 357   }\r
 358 \r
 359   Float_t chi2PerCluster = 0.;\r
 360   if(GetAnalysisMode() == AlidNdPtHelper::kTPC) { \r
 361     if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2Iter1()/Float_t(nClust);\r
 362   } else {\r
 363     chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);\r
 364   }\r
 365 \r
 366   Float_t clustPerFindClust = 0.;\r
 367   Int_t nFindableClust = esdTrack->GetTPCNclsF();\r
 368   if(nFindableClust>0.) clustPerFindClust = Float_t(nClust)/nFindableClust;\r
 369 \r
 370   Float_t b[2], bCov[3];\r
 371   esdTrack->GetImpactParameters(b,bCov);\r
 372 \r
 373   // kink idx\r
 374   Int_t kinkIdx = 0;\r
 375   //if(esdTrack->GetKinkIndex(0) > 0.)   isKink  = kTRUE;\r
 376   if(esdTrack->GetKinkIndex(0) > 0)      kinkIdx = 1;   // kink daughter\r
 377   else if(esdTrack->GetKinkIndex(0) < 0) kinkIdx = -1;  // kink mother\r
 378   else kinkIdx = 0; // not kink\r
 379 \r
 380   //printf("esdTrack->GetKinkIndex(0) %d \n", esdTrack->GetKinkIndex(0));\r
 381   //printf("esdTrack->GetKinkIndex(1) %d \n", esdTrack->GetKinkIndex(1));\r
 382   //printf("esdTrack->GetKinkIndex(2) %d \n", esdTrack->GetKinkIndex(2));\r
 383   //printf("kinkIdx %d \n", kinkIdx);\r
 384 \r
 385   //\r
 386   // Fill rec vs MC information\r
 387   //\r
 388   Bool_t isPrim = kTRUE;\r
 389   Bool_t hasStrangeMother = kFALSE;\r
 390   Bool_t isFromMaterial = kFALSE;\r
 391 \r
 392   if(IsUseMCInfo()) {\r
 393     if(!stack) return;\r
 394     Int_t label = TMath::Abs(esdTrack->GetLabel()); \r
 395     TParticle* particle = stack->Particle(label);\r
 396     if(!particle) return;\r
 397     if(particle->GetPDG() && particle->GetPDG()->Charge()==0.) return;\r
 398     isPrim = stack->IsPhysicalPrimary(label);\r
 399 \r
 400     // check whether has stange mother\r
 401     //\r
 402     Int_t motherPdg = -1; \r
 403     TParticle* mother = 0; \r
 404        \r
 405     Int_t motherLabel = particle->GetMother(0);  \r
 406     if(motherLabel>0) mother = stack->Particle(motherLabel); \r
 407     if(mother) motherPdg = TMath::Abs(mother->GetPdgCode()); // take abs for visualisation only \r
 408     Int_t mech = particle->GetUniqueID(); // production mechanism \r
 409 \r
 410     if( (motherPdg == 3122) || (motherPdg == -3122) || (motherPdg == 310)) // lambda, antilambda, k0s\r
 411     {\r
 412       if( (mech == 4) || (mech == 5) ) hasStrangeMother = kTRUE;\r
 413     } \r
 414     else {\r
 415       //if(isPrim==0 && mech == 13)   \r
 416       //printf("mech %d \n", mech);\r
 417       if(!isPrim) isFromMaterial = kTRUE; \r
 418     }\r
 419   }\r
 420 \r
 421   \r
 422   // fill histo\r
 423   Int_t charge = esdTrack->Charge();\r
 424   Double_t vRecMCTrackHist[12] = { nClust,chi2PerCluster,clustPerFindClust,b[0],b[1],eta,phi,pt,hasStrangeMother,isFromMaterial,isPrim,charge }; \r
 425   fRecMCTrackHist->Fill(vRecMCTrackHist);\r
 426 }\r
 427 \r
 428 //_____________________________________________________________________________\r
 429 Long64_t AlidNdPtCutAnalysisPbPb::Merge(TCollection* const list) \r
 430 {\r
 431   // Merge list of objects (needed by PROOF)\r
 432 \r
 433   if (!list)\r
 434   return 0;\r
 435 \r
 436   if (list->IsEmpty())\r
 437   return 1;\r
 438 \r
 439   TIterator* iter = list->MakeIterator();\r
 440   TObject* obj = 0;\r
 441 \r
 442   //TList *collPhysSelection = new TList;\r
 443 \r
 444   // collection of generated histograms\r
 445   Int_t count=0;\r
 446   while((obj = iter->Next()) != 0) {\r
 447     AlidNdPtCutAnalysisPbPb* entry = dynamic_cast<AlidNdPtCutAnalysisPbPb*>(obj);\r
 448     if (entry == 0) continue; \r
 449   \r
 450     // event histo\r
 451     fEventCount->Add(entry->fEventCount);\r
 452     fRecEventHist->Add(entry->fRecEventHist);\r
 453     fRecMCEventHist->Add(entry->fRecMCEventHist);\r
 454     fMCEventHist->Add(entry->fMCEventHist);\r
 455 \r
 456     // track histo\r
 457     fRecMCTrackHist->Add(entry->fRecMCTrackHist);\r
 458 \r
 459     // physics selection\r
 460     //collPhysSelection->Add(entry->GetPhysicsTriggerSelection());\r
 461     \r
 462   count++;\r
 463   }\r
 464 \r
 465   //AliPhysicsSelection *trigSelection = GetPhysicsTriggerSelection();\r
 466   //trigSelection->Merge(collPhysSelection);\r
 467 \r
 468   //if(collPhysSelection) delete collPhysSelection;\r
 469 \r
 470 return count;\r
 471 }\r
 472 \r
 473 //_____________________________________________________________________________\r
 474 void AlidNdPtCutAnalysisPbPb::Analyse() \r
 475 {\r
 476   //\r
 477   // Analyse histograms\r
 478   //\r
 479   TH1::AddDirectory(kFALSE);\r
 480   TObjArray *aFolderObj = new TObjArray;\r
 481   TH1D *h1D = 0; \r
 482   TH2D *h2D = 0; \r
 483 \r
 484 \r
 485   //\r
 486   // get cuts\r
 487   //\r
 488   AlidNdPtEventCuts *evtCuts = GetEventCuts(); \r
 489   AlidNdPtAcceptanceCuts *accCuts = GetAcceptanceCuts(); \r
 490   AliESDtrackCuts *esdTrackCuts = GetTrackCuts(); \r
 491 \r
 492   if(!evtCuts || !accCuts || !esdTrackCuts) {\r
 493     Error("AlidNdPtCutAnalysisPbPb::Analyse()", "cuts not available");\r
 494     return;\r
 495   }\r
 496 \r
 497   //\r
 498   // set min and max values\r
 499   //\r
 500   Double_t minPt = accCuts->GetMinPt();\r
 501   Double_t maxPt = accCuts->GetMaxPt();\r
 502   Double_t minEta = accCuts->GetMinEta();\r
 503   Double_t maxEta = accCuts->GetMaxEta()-0.00001;\r
 504 \r
 505   Double_t maxDCAr = accCuts->GetMaxDCAr();\r
 506 \r
 507   //\r
 508   // Event counters\r
 509   //\r
 510   h2D = (TH2D*)fEventCount->Projection(0,1);\r
 511   h2D->SetName("trig_vs_trigANDvertex");\r
 512   aFolderObj->Add(h2D);\r
 513 \r
 514   fEventCount->GetAxis(0)->SetRange(2,2); // triggered\r
 515   h1D = (TH1D*)fEventCount->Projection(1);\r
 516   h1D->SetTitle("rec. vertex for triggered events");\r
 517   h1D->SetName("trigANDvertex");\r
 518   aFolderObj->Add(h1D);\r
 519 \r
 520   //\r
 521   // Create rec. event histograms\r
 522   //\r
 523   h1D = (TH1D *)fRecEventHist->Projection(0);\r
 524   h1D->SetName("rec_xv");\r
 525   aFolderObj->Add(h1D);\r
 526 \r
 527   h1D = (TH1D *)fRecEventHist->Projection(1);\r
 528   h1D->SetName("rec_yv");\r
 529   aFolderObj->Add(h1D);\r
 530 \r
 531   h1D = (TH1D *)fRecEventHist->Projection(2);\r
 532   h1D->SetName("rec_zv");\r
 533   aFolderObj->Add(h1D);\r
 534 \r
 535   h2D = (TH2D *)fRecEventHist->Projection(3,4);\r
 536   h2D->SetName("rec_resZv_vs_Mult");\r
 537   aFolderObj->Add(h2D);\r
 538 \r
 539   h2D = (TH2D *)fRecEventHist->Projection(0,1);\r
 540   h2D->SetName("rec_xv_vs_yv");\r
 541   aFolderObj->Add(h2D);\r
 542 \r
 543   h2D = (TH2D *)fRecEventHist->Projection(0,2);\r
 544   h2D->SetName("rec_xv_vs_zv");\r
 545   aFolderObj->Add(h2D);\r
 546 \r
 547   h2D = (TH2D *)fRecEventHist->Projection(3,4);\r
 548   h2D->SetName("rec_resZv_vs_Mult");\r
 549   aFolderObj->Add(h2D);\r
 550 \r
 551   //\r
 552   // MC available\r
 553   //\r
 554   if(IsUseMCInfo()) {\r
 555 \r
 556   //\r
 557   // Create mc event histograms\r
 558   //\r
 559   h2D = (TH2D *)fMCEventHist->Projection(0,1);\r
 560   h2D->SetName("mc_xv_vs_yv");\r
 561   aFolderObj->Add(h2D);\r
 562 \r
 563   h2D = (TH2D *)fMCEventHist->Projection(0,2);\r
 564   h2D->SetName("mc_xv_vs_zv");\r
 565   aFolderObj->Add(h2D);\r
 566 \r
 567   //\r
 568   // Create rec-mc event histograms\r
 569   //\r
 570   h2D = (TH2D *)fRecMCEventHist->Projection(0,3);\r
 571   h2D->SetName("rec_mc_deltaXv_vs_mult");\r
 572   aFolderObj->Add(h2D);\r
 573 \r
 574   h2D = (TH2D *)fRecMCEventHist->Projection(1,3);\r
 575   h2D->SetName("rec_mc_deltaYv_vs_mult");\r
 576   aFolderObj->Add(h2D);\r
 577 \r
 578   h2D = (TH2D *)fRecMCEventHist->Projection(2,3);\r
 579   h2D->SetName("rec_mc_deltaZv_vs_mult");\r
 580   aFolderObj->Add(h2D);\r
 581 \r
 582   } // end use MC info \r
 583 \r
 584 \r
 585 \r
 586   //\r
 587   // Create rec-mc track track histograms \r
 588   //\r
 589 \r
 590   // DCA cuts\r
 591   fRecMCTrackHist->GetAxis(3)->SetRangeUser(-maxDCAr,maxDCAr);\r
 592   fRecMCTrackHist->GetAxis(4)->SetRangeUser(-maxDCAr,maxDCAr);\r
 593 \r
 594   h2D = (TH2D *)fRecMCTrackHist->Projection(7,5);\r
 595   h2D->SetName("pt_vs_eta");\r
 596   aFolderObj->Add(h2D);\r
 597 \r
 598   fRecMCTrackHist->GetAxis(7)->SetRangeUser(minPt,maxPt);  \r
 599 \r
 600   h2D = (TH2D *)fRecMCTrackHist->Projection(0,5);\r
 601   h2D->SetName("nClust_vs_eta");\r
 602   aFolderObj->Add(h2D);\r
 603 \r
 604   h2D = (TH2D *)fRecMCTrackHist->Projection(1,5);\r
 605   h2D->SetName("chi2PerClust_vs_eta");\r
 606   aFolderObj->Add(h2D);\r
 607 \r
 608   h2D = (TH2D *)fRecMCTrackHist->Projection(2,5);\r
 609   h2D->SetName("ratio_nClust_nFindableClust_vs_eta");\r
 610   aFolderObj->Add(h2D);\r
 611 \r
 612   h2D = (TH2D *)fRecMCTrackHist->Projection(5,6);\r
 613   h2D->SetName("eta_vs_phi");\r
 614   aFolderObj->Add(h2D);\r
 615 \r
 616   //\r
 617   fRecMCTrackHist->GetAxis(5)->SetRangeUser(minEta,maxEta);  \r
 618 \r
 619   h2D = (TH2D *)fRecMCTrackHist->Projection(0,6);\r
 620   h2D->SetName("nClust_vs_phi");\r
 621   aFolderObj->Add(h2D);\r
 622 \r
 623   h2D = (TH2D *)fRecMCTrackHist->Projection(1,6);\r
 624   h2D->SetName("chi2PerClust_vs_phi");\r
 625   aFolderObj->Add(h2D);\r
 626 \r
 627   h2D = (TH2D *)fRecMCTrackHist->Projection(2,6);\r
 628   h2D->SetName("ratio_nClust_nFindableClust_vs_phi");\r
 629   aFolderObj->Add(h2D);\r
 630 \r
 631   //\r
 632   fRecMCTrackHist->GetAxis(7)->SetRangeUser(0.0,maxPt);  \r
 633 \r
 634   h2D = (TH2D *)fRecMCTrackHist->Projection(0,7);\r
 635   h2D->SetName("nClust_vs_pt");\r
 636   aFolderObj->Add(h2D);\r
 637 \r
 638   h2D = (TH2D *)fRecMCTrackHist->Projection(1,7);\r
 639   h2D->SetName("chi2PerClust_vs_pt");\r
 640   aFolderObj->Add(h2D);\r
 641 \r
 642   h2D = (TH2D *)fRecMCTrackHist->Projection(2,7);\r
 643   h2D->SetName("ratio_nClust_nFindableClust_vs_pt");\r
 644   aFolderObj->Add(h2D);\r
 645 \r
 646   h2D = (TH2D *)fRecMCTrackHist->Projection(6,7);\r
 647   h2D->SetName("phi_vs_pt");\r
 648   aFolderObj->Add(h2D);\r
 649 \r
 650 \r
 651   // fiducial volume\r
 652   fRecMCTrackHist->GetAxis(5)->SetRangeUser(minEta,maxEta);  \r
 653   fRecMCTrackHist->GetAxis(7)->SetRangeUser(minPt,maxPt);  \r
 654 \r
 655   // DCA cuts\r
 656   fRecMCTrackHist->GetAxis(3)->SetRangeUser(-maxDCAr,maxDCAr);\r
 657   fRecMCTrackHist->GetAxis(4)->SetRangeUser(-maxDCAr,maxDCAr);\r
 658 \r
 659   h2D = (TH2D *)fRecMCTrackHist->Projection(0,1);\r
 660   h2D->SetName("nClust_vs_chi2PerClust");\r
 661   aFolderObj->Add(h2D);\r
 662 \r
 663   h2D = (TH2D *)fRecMCTrackHist->Projection(0,2);\r
 664   h2D->SetName("nClust_vs_ratio_nClust_nFindableClust");\r
 665   aFolderObj->Add(h2D);\r
 666 \r
 667   //\r
 668   // DCAy cuts\r
 669   //\r
 670   fRecMCTrackHist->GetAxis(0)->SetRange(50,160); // nClust/track > 50\r
 671   fRecMCTrackHist->GetAxis(1)->SetRangeUser(0.,3.9999); // chi2/cluster < 4.0\r
 672   fRecMCTrackHist->GetAxis(3)->SetRange(1,fRecMCTrackHist->GetAxis(3)->GetNbins());\r
 673   //fRecMCTrackHist->GetAxis(4)->SetRangeUser(-1.0,1.0);\r
 674   fRecMCTrackHist->GetAxis(4)->SetRange(1,fRecMCTrackHist->GetAxis(4)->GetNbins());\r
 675 \r
 676   // sec\r
 677   fRecMCTrackHist->GetAxis(9)->SetRange(1,1);\r
 678   h1D = (TH1D *)fRecMCTrackHist->Projection(3);\r
 679   h1D->SetName("dcay_sec");\r
 680   aFolderObj->Add(h1D);\r
 681 \r
 682   // prim\r
 683   fRecMCTrackHist->GetAxis(9)->SetRange(2,2);\r
 684   h1D = (TH1D *)fRecMCTrackHist->Projection(3);\r
 685   h1D->SetName("dcay_prim");\r
 686   aFolderObj->Add(h1D);\r
 687 \r
 688   // DCAz cuts\r
 689   //fRecMCTrackHist->GetAxis(3)->SetRangeUser(-1.0,1.0);\r
 690   fRecMCTrackHist->GetAxis(4)->SetRange(1,fRecMCTrackHist->GetAxis(4)->GetNbins());\r
 691 \r
 692   // sec\r
 693   fRecMCTrackHist->GetAxis(9)->SetRange(1,1);\r
 694   h1D = (TH1D *)fRecMCTrackHist->Projection(4);\r
 695   h1D->SetName("dcaz_sec");\r
 696   aFolderObj->Add(h1D);\r
 697 \r
 698   // prim\r
 699   fRecMCTrackHist->GetAxis(9)->SetRange(2,2);\r
 700   h1D = (TH1D *)fRecMCTrackHist->Projection(4);\r
 701   h1D->SetName("dcaz_prim");\r
 702   aFolderObj->Add(h1D);\r
 703 \r
 704 \r
 705   // export objects to analysis folder\r
 706   fAnalysisFolder = ExportToFolder(aFolderObj);\r
 707 \r
 708   // delete only TObjArray\r
 709   if(aFolderObj) delete aFolderObj;\r
 710 }\r
 711 \r
 712 //_____________________________________________________________________________\r
 713 TFolder* AlidNdPtCutAnalysisPbPb::ExportToFolder(TObjArray * const array) \r
 714 {\r
 715   // recreate folder avery time and export objects to new one\r
 716   //\r
 717   AlidNdPtCutAnalysisPbPb * comp=this;\r
 718   TFolder *folder = comp->GetAnalysisFolder();\r
 719 \r
 720   TString name, title;\r
 721   TFolder *newFolder = 0;\r
 722   Int_t i = 0;\r
 723   Int_t size = array->GetSize();\r
 724 \r
 725   if(folder) { \r
 726      // get name and title from old folder\r
 727      name = folder->GetName();  \r
 728      title = folder->GetTitle();  \r
 729 \r
 730          // delete old one\r
 731      delete folder;\r
 732 \r
 733          // create new one\r
 734      newFolder = CreateFolder(name.Data(),title.Data());\r
 735      newFolder->SetOwner();\r
 736 \r
 737          // add objects to folder\r
 738      while(i < size) {\r
 739            newFolder->Add(array->At(i));\r
 740            i++;\r
 741          }\r
 742   }\r
 743 \r
 744 return newFolder;\r
 745 }\r
 746 \r
 747 //_____________________________________________________________________________\r
 748 TFolder* AlidNdPtCutAnalysisPbPb::CreateFolder(TString name,TString title) { \r
 749 // create folder for analysed histograms\r
 750 //\r
 751 TFolder *folder = 0;\r
 752   folder = new TFolder(name.Data(),title.Data());\r
 753 \r
 754   return folder;\r
 755 }\r