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