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