Added possibility to do cut on TPCnSigma and rapidity in the task.
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / PiKaPr / TPCTOF / AliAnalysisCombinedHadronSpectra.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, proviyaded 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 purapose. It is         *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 ///////////////////////////////////////////////////////////////////////////
17 //                                                                       //
18 //                                                                       //
19 // Analysis for identified charged hadron spectra.                       //
20 //                                                                       //
21 //                                                                       //
22 ///////////////////////////////////////////////////////////////////////////
23
24 #include "Riostream.h"
25 #include "TChain.h"
26 #include "TTree.h"
27 #include "TH1F.h"
28 #include "TH2D.h"
29 #include "TH3F.h"
30 #include "TList.h"
31 #include "TMath.h"
32 #include "TCanvas.h"
33 #include "TObjArray.h"
34 #include "TF1.h"
35 #include "TFile.h"
36
37 #include "AliAnalysisTaskSE.h"
38 #include "AliAnalysisManager.h"
39
40 #include "AliHeader.h"
41 #include "AliGenPythiaEventHeader.h"
42 #include "AliGenCocktailEventHeader.h"
43
44 #include "AliPID.h"
45 #include "AliESDtrackCuts.h"
46 #include "AliESDVertex.h"
47 #include "AliESDEvent.h"
48 #include "AliESDInputHandler.h"
49 #include "AliESDtrack.h"
50 #include "AliESDpid.h"
51 #include "AliCentrality.h"
52 #include "AliESDUtils.h"
53 #include "AliMultiplicity.h"
54
55 #include "AliMCEventHandler.h"
56 #include "AliMCEvent.h"
57 #include "AliStack.h"
58
59 #include "AliLog.h"
60
61 #include "AliAnalysisCombinedHadronSpectra.h"
62
63
64 ClassImp(AliAnalysisCombinedHadronSpectra)
65
66 //________________________________________________________________________
67 AliAnalysisCombinedHadronSpectra::AliAnalysisCombinedHadronSpectra() 
68   : AliAnalysisTaskSE("TaskChargedHadron"), fESD(0), fListHist(0), fESDtrackCuts(0),fESDTrackCutsMult(0),fESDpid(0),
69     fMCtrue(0),
70     fOnlyQA(0),
71     fUseHBTmultiplicity(0),
72     fUseTPConlyTracks(0),
73     fSaveMotherPDG(0),
74   fSmallTHnSparse(0),
75   fTPCnSigmaCutLow(0),
76   fTPCnSigmaCutHigh(0),
77   fRapidityCutLow(0),
78   fRapidityCutHigh(0),
79     fAlephParameters(),
80     fHistRealTracks(0),
81     fHistMCparticles(0),
82     fHistPidQA(0),
83     fHistMult(0),
84     fHistCentrality(0)
85 {
86   // default Constructor
87   /* fast compilation test
88      gSystem->Load("libANALYSIS");
89      gSystem->Load("libANALYSISalice");
90      .L /d/alice09/akalweit/train/trunk/akalweit_hadronSpectra/AliAnalysisCombinedHadronSpectra.cxx++
91    */
92 }
93
94
95 //________________________________________________________________________
96 AliAnalysisCombinedHadronSpectra::AliAnalysisCombinedHadronSpectra(const char *name) 
97   : AliAnalysisTaskSE(name), fESD(0), fListHist(0), fESDtrackCuts(0),fESDTrackCutsMult(0),fESDpid(0),
98     fMCtrue(0),
99     fOnlyQA(0),
100     fUseHBTmultiplicity(0),
101     fUseTPConlyTracks(0),
102     fSaveMotherPDG(0),
103   fSmallTHnSparse(0),
104   fTPCnSigmaCutLow(0),
105   fTPCnSigmaCutHigh(0),
106   fRapidityCutLow(0),
107   fRapidityCutHigh(0),
108     fAlephParameters(),
109     fHistRealTracks(0),
110     fHistMCparticles(0),
111     fHistPidQA(0),
112     fHistMult(0),
113     fHistCentrality(0)
114     {
115   //
116   // standard constructur which should be used
117   //
118   Printf("*** CONSTRUCTOR CALLED ****");
119   //
120   fMCtrue = kTRUE; 
121   fOnlyQA = kFALSE;
122   fUseHBTmultiplicity = kTRUE;
123   fUseTPConlyTracks = kFALSE;
124
125   fSmallTHnSparse = kFALSE;
126   fTPCnSigmaCutLow = -3.;
127   fTPCnSigmaCutHigh = 3.;
128   fRapidityCutLow = -0.2;
129   fRapidityCutHigh = 0.2;
130   /* real */
131   fAlephParameters[0] = 0.0283086;
132   fAlephParameters[1] = 2.63394e+01;
133   fAlephParameters[2] = 5.04114e-11;
134   fAlephParameters[3] = 2.12543e+00;
135   fAlephParameters[4] = 4.88663e+00;
136   //
137   // initialize PID object
138   //
139   //fESDpid = new AliESDpid();
140   //
141   // create track cuts
142   //
143   fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
144   //
145   //Initialize();
146   // Output slot #0 writes into a TList container
147   DefineOutput(1, TList::Class());
148
149 }
150
151
152 //________________________________________________________________________
153 void AliAnalysisCombinedHadronSpectra::Initialize()
154 {
155   //
156   // updating parameters in case of changes
157   //
158   // 1. pid parameters
159   //
160   //fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
161   //
162   // 2. track cuts
163   //
164   /*
165   fESDtrackCuts->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2);  // BEWARE STANDARD VALUES ARE: 2, 2, 0.5, 0.5, 2
166   fESDtrackCuts->SetMaxNsigmaToVertex(3);
167   fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
168   fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
169   fESDtrackCuts->SetMinNClustersTPC(70);
170   fESDtrackCuts->SetMaxChi2PerClusterTPC(4);
171   fESDtrackCuts->SetMaxDCAToVertexXY(3);
172   fESDtrackCuts->SetMaxDCAToVertexZ(3);
173   fESDtrackCuts->SetRequireTPCRefit(kTRUE);
174   fESDtrackCuts->SetRequireITSRefit(kTRUE);
175   fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny); //TEMPORARY <-> REMOVE
176   fESDtrackCuts->SetMinNClustersITS(3);
177   */
178   //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); // kTRUE = sel. primaries --> patch for the moment, do TFractionFitter later
179
180
181   if (!fUseTPConlyTracks) {
182           fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kFALSE);
183           fESDtrackCuts->SetMaxDCAToVertexXY(3);
184           fESDtrackCuts->SetMaxDCAToVertexZ(2);
185           fESDtrackCuts->SetEtaRange(-0.9,0.9);
186   }
187   else {
188           //fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
189           fESDtrackCuts->SetMinNClustersTPC(70);
190           fESDtrackCuts->SetMaxChi2PerClusterTPC(4);
191           fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
192           fESDtrackCuts->SetRequireTPCRefit(kFALSE);
193
194           fESDtrackCuts->SetMaxDCAToVertexXY(15);
195           fESDtrackCuts->SetMaxDCAToVertexZ(6);
196           fESDtrackCuts->SetDCAToVertex2D(kFALSE);
197           fESDtrackCuts->SetRequireSigmaToVertex(kFALSE);
198
199           fESDtrackCuts->SetEtaRange(-0.9,0.9);
200   }
201   //
202   //
203   //
204   //
205   fESDTrackCutsMult = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
206   fESDTrackCutsMult->SetEtaRange(-1.2,+1.2);
207   fESDTrackCutsMult->SetPtRange(0.15,1e10);
208
209 }
210
211
212 //________________________________________________________________________
213 void AliAnalysisCombinedHadronSpectra::UserCreateOutputObjects() 
214 {
215   // Create histograms
216   // Called once
217   fListHist = new TList();
218   fListHist->SetOwner(kTRUE);
219   //
220   const Int_t kPtBins = 35;
221   const Int_t kMultBins = 11;
222   const Int_t kDcaBins = 76;
223   const Float_t kDcaBinsTPConlyFactor = 5; //need to change binning of DCA plot for tpconly
224   // sort pT-bins ..
225   Double_t binsPt[kPtBins+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};
226
227   Double_t binsDca[kDcaBins+1] =  {-3,-2.85,-2.7,-2.55,-2.4,-2.25,-2.1,-1.95,-1.8,-1.65,-1.5,-1.35,-1.2,-1.05,-0.9,-0.75,-0.6,-0.45,-0.3,-0.285,-0.27,-0.255,-0.24,-0.225,-0.21,-0.195,-0.18,-0.165,-0.15,-0.135,-0.12,-0.105,-0.09,-0.075,-0.06,-0.045,-0.03,-0.015,0,0.015,0.03,0.045,0.06,0.075,0.09,0.105,0.12,0.135,0.15,0.165,0.18,0.195,0.21,0.225,0.24,0.255,0.27,0.285,0.3,0.45,0.6,0.75,0.9,1.05,1.2,1.35,1.5,1.65,1.8,1.95,2.1,2.25,2.4,2.55,2.7,2.85,3};
228
229   // DCA bins borders get multiplied by constant factor for TPConlyTracks
230   Double_t binsDcaTPConly[kDcaBins+1];
231   for (Int_t i = 0; i< kDcaBins+1; i++) {
232         binsDcaTPConly[i] = kDcaBinsTPConlyFactor * binsDca[i];
233   }
234
235   //
236   // create the histograms with all necessary information --> it is filled 4x for each particle assumption
237   //
238   // (0.) assumed particle: 0. pion, 1. kaon, 2. proton, 3. deuteron
239   // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
240   // (2.) pT
241   // (3.) sign
242   // (4.) rapidity --> filled 4x
243   // (5.)  pull TPC dEx --> filled 4x
244   // (6.) has valid TOF pid signal
245   // (7.) nsigma TOF --> filled 4x
246   // (8..) dca_xy
247   // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident prim, 3-second weak, 4-second material, 5-misident sec, 6-sec. K0, 7-sec. lambda, 8-sec sigma+
248   //
249   
250   //dimensions of standard THnSparse
251   //                              0,           1,           2,  3,       4,    5,     6,   7,     8
252   Int_t    binsHistReal[9] = {   3,   kMultBins,     kPtBins,  2,       10,    50,    2,  80,    kDcaBins};
253   Double_t xminHistReal[9] = {-0.5,        -0.5,           0, -2,      -0.5,   -5, -0.5,  -8,       -3};
254   Double_t xmaxHistReal[9] = { 2.5,        10.5,           3,  2,       0.5,    5,  1.5,   8,        3};
255
256   //dimensions of small THnSparse
257   //                              0,           1,           2,  3,                        4,   5,       6
258   Int_t    binsHistRealSm[7] = {   3,   kMultBins,     kPtBins,  2,  /*    10,   50,*/    2,  80, kDcaBins};
259   Double_t xminHistRealSm[7] = {-0.5,        -0.5,           0, -2,  /*  -0.5,   -5,*/ -0.5,  -8,       -3};
260   Double_t xmaxHistRealSm[7] = { 2.5,        10.5,           3,  2,  /*   0.5,    5,*/  1.5,   8,        3};
261
262   if (!fSmallTHnSparse) fHistRealTracks = new THnSparseF("fHistRealTracks","real tracks",9,binsHistReal,xminHistReal,xmaxHistReal);
263   else                  fHistRealTracks = new THnSparseF("fHistRealTracks","real tracks",7,binsHistRealSm,xminHistRealSm,xmaxHistRealSm);
264   //
265   fHistRealTracks->GetAxis(2)->Set(kPtBins, binsPt);
266
267   //different DCAxy binning for TPConlyTracks
268
269   if (!fUseTPConlyTracks){
270     if (!fSmallTHnSparse) fHistRealTracks->GetAxis(8)->Set(kDcaBins, binsDca);
271     else                  fHistRealTracks->GetAxis(6)->Set(kDcaBins, binsDca);
272   }
273   else {
274     if (!fSmallTHnSparse) fHistRealTracks->GetAxis(8)->Set(kDcaBins, binsDcaTPConly);
275     else                  fHistRealTracks->GetAxis(6)->Set(kDcaBins, binsDcaTPConly);
276   }
277   fListHist->Add(fHistRealTracks);
278   //
279   //                      0.ptot,1.tpcSig,2.hasTOF, 3. assumed part., 4. nclDedx, 5. nSigmaTPC (4x), 6. nSigmaTOF (4x), 7. centrality
280   fHistPidQA = new TH3D("fHistPidQA","PID QA",500,0.1,10,1000,0,1000,2,-2,2);
281   BinLogAxis(fHistPidQA);
282   fListHist->Add(fHistPidQA);
283   
284   // dimensions of standard THnSparse
285   //                            0,            1,           2,  3,      4,   5,    6,   7,       8.,    9.
286   Int_t    binsHistMC[10] = {   3,    kMultBins,     kPtBins,  2,     10,  50,    2,  80, kDcaBins,    6};
287   Double_t xminHistMC[10] = {-0.5,         -0.5,           0, -2,   -0.5,  -5, -0.5,  -8,       -3, -0.5};
288   Double_t xmaxHistMC[10] = { 2.5,         10.5,           3,  2,    0.5,   5,  1.5,   8,        3,  5.5};
289
290   //dimensions of small THnSparse
291   //                              0,           1,           2,  3,                     4,   5,       6.,    7.
292   Int_t    binsHistMCSm[8] = {   3,    kMultBins,     kPtBins,  2,  /*   10,  50,*/    2,  80, kDcaBins,    6};
293   Double_t xminHistMCSm[8] = {-0.5,         -0.5,           0, -2,  /* -0.5,  -5,*/ -0.5,  -8,       -3, -0.5};
294   Double_t xmaxHistMCSm[8] = { 2.5,         10.5,           3,  2,  /*  0.5,   5,*/  1.5,   8,        3,  5.5};
295
296   //different binning for CODE axis, if we want to save motherPDG
297   if (fSaveMotherPDG) {
298     binsHistMC[9] = 9;
299     xmaxHistMC[9] = 8.5;
300     binsHistMCSm[7] = 9;
301     xmaxHistMCSm[7] = 8.5;
302
303   }
304
305   if (!fSmallTHnSparse) fHistMCparticles = new THnSparseF("fHistMCparticles","MC histogram",10,binsHistMC,xminHistMC,xmaxHistMC);
306   else                  fHistMCparticles = new THnSparseF("fHistMCparticles","MC histogram",8,binsHistMCSm,xminHistMCSm,xmaxHistMCSm);
307
308
309   fHistMCparticles->GetAxis(2)->Set(kPtBins, binsPt);
310
311   //different DCAxy binning for TPConlyTracks
312   if (!fUseTPConlyTracks) {
313     if (!fSmallTHnSparse) fHistMCparticles->GetAxis(8)->Set(kDcaBins, binsDca);
314     else                  fHistMCparticles->GetAxis(6)->Set(kDcaBins, binsDca);
315   }
316   else {
317     if (!fSmallTHnSparse) fHistMCparticles->GetAxis(8)->Set(kDcaBins, binsDcaTPConly);
318     else                  fHistMCparticles->GetAxis(6)->Set(kDcaBins, binsDcaTPConly);
319   }
320
321
322   fListHist->Add(fHistMCparticles);
323   //
324   fHistMult = new TH2D("fHistMult", "control histogram to count number of events", 502, -2.5, 499.5,4,-0.5,3.5);
325   fHistCentrality = new TH1D("fHistCentrality", "control histogram to count number of events", 22, -1.5, 20.5);
326   fListHist->Add(fHistMult);
327   fListHist->Add(fHistCentrality);
328   
329 }
330
331 //________________________________________________________________________
332 void AliAnalysisCombinedHadronSpectra::UserExec(Option_t *) 
333 {
334   //
335   // main event loop
336   //
337   if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
338   if (!fESDpid) {
339     fESDpid = new AliESDpid(); // HACK FOR MC PBPB --> PLEASE REMOVE AS SOON AS POSSIBLE
340     fESDpid->GetTPCResponse().SetBetheBlochParameters(1.28778e+00/50., 3.13539e+01, TMath::Exp(-3.16327e+01), 1.87901e+00, 6.41583e+00);
341   }
342   //AliLog::SetGlobalLogLevel(AliLog::kError);
343   //
344   // Check Monte Carlo information and other access first:
345   //
346   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
347   if (!eventHandler) {
348     //Printf("ERROR: Could not retrieve MC event handler");
349     fMCtrue = kFALSE;
350   }
351   //
352   AliMCEvent* mcEvent = 0x0;
353   AliStack* stack = 0x0;
354   if (eventHandler) mcEvent = eventHandler->MCEvent();
355   if (!mcEvent) {
356     //Printf("ERROR: Could not retrieve MC event");
357     if (fMCtrue) return;
358   }
359   if (fMCtrue) {
360     stack = mcEvent->Stack();
361     if (!stack) return;
362   }
363   //
364   fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
365   if (!fESD) {
366     //Printf("ERROR: fESD not available");
367     return;
368   }
369   
370   if (!fESDtrackCuts) {
371     Printf("ERROR: fESDtrackCuts not available");
372     return;
373   }
374   //
375   // check if event is selected by physics selection class
376   //
377   Bool_t isSelected = kTRUE; // for reasons of backward compatibility --> check is now in AddTask macro
378   /*
379   Bool_t isSelected = kFALSE;
380   isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) == AliVEvent::kMB;
381   */
382   //
383   // monitor vertex position
384   //
385   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
386   if(vertex->GetNContributors()<1) {
387     // SPD vertex
388     vertex = fESD->GetPrimaryVertexSPD();
389     if(vertex->GetNContributors()<1) vertex = 0x0;
390   }  
391   //
392   // small track loop to determine trigger Pt, multiplicity or centrality
393   //
394   Double_t triggerPt = 0;
395   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
396     AliESDtrack *track =fESD->GetTrack(i);
397     if (!fESDTrackCutsMult->AcceptTrack(track)) continue;
398     if (track->Pt() > triggerPt) triggerPt = track->Pt();
399   }
400   Int_t trackCounter = fESDTrackCutsMult->CountAcceptedTracks(fESD);
401   //
402   // 2nd multiplicity estimator SPD hits; replaces multiplicity from HBT
403   //
404   Float_t spdCorr = -1;
405   const AliMultiplicity *mult = fESD->GetMultiplicity();
406   Float_t nClusters[6]={0.0,0.0,0.0,0.0,0.0,0.0};
407   for(Int_t ilay=0; ilay<6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
408   if (vertex) spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],vertex->GetZ());
409   //
410   Float_t centrality = -1;
411   //
412   // IMPORTANT CENTRALITY DEFINITION FOR pp
413   //
414   if (!fUseHBTmultiplicity) {
415     // test binning to check if multiplicity-dependence is due to diffractive events being s
416     if (trackCounter >= 0  && trackCounter <= 0)  centrality = 0;
417     if (trackCounter >= 1  && trackCounter <= 1)  centrality = 1;
418     if (trackCounter >= 2  && trackCounter <= 2)  centrality = 2;
419     if (trackCounter >= 3  && trackCounter <= 3)  centrality = 3;
420     if (trackCounter >= 4  && trackCounter <= 4)  centrality = 4;
421     if (trackCounter >= 5  && trackCounter <= 5)  centrality = 5;
422     // this was all the original bin 1 being [0..5]
423     if (trackCounter >= 6  && trackCounter <= 9)  centrality = 6;
424     if (trackCounter >= 10 && trackCounter <= 14) centrality = 7;
425     if (trackCounter >= 15 && trackCounter <= 22) centrality = 8;
426     if (trackCounter >= 23 && trackCounter <= 32) centrality = 9;
427     if (trackCounter >= 33 && trackCounter <= 42) centrality = 10;
428     /*
429     if (trackCounter >= 43 && trackCounter <= 52) centrality = 7
430     if (trackCounter >= 53 && trackCounter <= 62) centrality = 8;
431     if (trackCounter >= 63 && trackCounter <= 72) centrality = 9;
432     if (trackCounter >= 73 && trackCounter <= 82) centrality = 10;
433     */
434   } else {
435     if (spdCorr >= 0  && spdCorr <=  2)  centrality  = 0;
436     if (spdCorr >= 3  && spdCorr <=  5)  centrality  = 1;
437     if (spdCorr >= 6  && spdCorr <=  8)  centrality  = 2;
438     if (spdCorr >= 9  && spdCorr <= 11)  centrality  = 3;
439     if (spdCorr >= 12 && spdCorr <= 14)  centrality  = 4;
440     if (spdCorr >= 15 && spdCorr <= 16)  centrality  = 5;
441     // this was all the original bin 1 being [0..16]
442     if (spdCorr >= 17 && spdCorr <= 30)  centrality =  6;
443     if (spdCorr >= 31 && spdCorr <= 45)  centrality =  7;
444     if (spdCorr >= 46 && spdCorr <= 68)  centrality =  8;
445     if (spdCorr >= 69 && spdCorr <= 97)  centrality =  9;
446     if (spdCorr >= 98)                   centrality = 10;
447     /*
448     if (spdCorr >= 17 && spdCorr <= 30)  centrality = 2;
449     if (spdCorr >= 31 && spdCorr <= 45)  centrality = 3;
450     if (spdCorr >= 46 && spdCorr <= 68)  centrality = 4;
451     if (spdCorr >= 69 && spdCorr <= 97)  centrality = 5;
452     if (spdCorr >= 98)                   centrality = 6;
453     */
454   }
455   //
456   Int_t rootS = fESD->GetBeamEnergy() < 1000 ? 0 : 1;
457   if (fESD->GetEventSpecie() == 4) { // PbPb
458     rootS = 2;
459     AliCentrality *esdCentrality = fESD->GetCentrality();
460     centrality = esdCentrality->GetCentralityClass10("V0M") + 1; // centrality percentile determined with V0
461     if (TMath::Abs(centrality - 1) < 1e-5) {
462       centrality = esdCentrality->GetCentralityClass5("V0M");
463     }
464   }
465   Int_t nContributors = 0;
466   if (fESD->GetPrimaryVertexTPC()) nContributors = fESD->GetPrimaryVertexTPC()->GetNContributors();
467   //
468   Int_t processtype = 0;
469   Int_t processCode = 0;
470   //
471   // important change: fill generated only after vertex cut in case of heavy-ions
472   //
473   if (!vertex && fESD->GetEventSpecie() == 4) {
474     fHistMult->Fill(-1, processCode);
475     PostData(1, fListHist);
476     return;
477   } else {
478     if (TMath::Abs(vertex->GetZv()) > 10 && fESD->GetEventSpecie() == 4) {
479       fHistMult->Fill(-1, processCode);
480       PostData(1, fListHist);
481       return;
482     }
483   }
484   //
485   if (fMCtrue) {
486     //
487     //
488     //
489     AliHeader * header = mcEvent->Header();
490     processtype = GetPythiaEventProcessType(header);
491     // non diffractive
492     if (processtype !=92 && processtype !=93 && processtype != 94) processCode = 1;
493     // single diffractive
494     if ((processtype == 92 || processtype == 93)) processCode = 2;
495     // double diffractive
496     if (processtype == 94) processCode = 3;
497     //
498     for(Int_t i = 0; i < stack->GetNtrack(); i++) {
499       TParticle * trackMC = stack->Particle(i);
500       Int_t pdg = trackMC->GetPdgCode();
501       //
502       Double_t xv = trackMC->Vx();
503       Double_t yv = trackMC->Vy();
504       Double_t zv = trackMC->Vz();
505       Double_t dxy = 0;
506       dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
507       Double_t dz = 0;
508       dz = TMath::Abs(zv); // so stupid to avoid warnings
509       //
510       // vertex cut - selection of primaries
511       //
512       //if (dxy > 3 || dz > 10) continue; // fixed cut at 3cm in r
513       //
514       if (!stack->IsPhysicalPrimary(i)) continue;
515       //
516       // fill MC histograms here...
517       // 
518       Double_t rap = trackMC->Y();
519       Double_t pT  = trackMC->Pt();
520       Int_t sign = pdg < 0 ? -1 : 1; // only works for charged pi,K,p !!
521 //      Double_t transMass = TMath::Sqrt(trackMC->Pt()*trackMC->Pt() + trackMC->GetMass()*trackMC->GetMass()) - trackMC->GetMass();
522       //
523       Int_t iPart = -1;
524       if (TMath::Abs(pdg) == 211)  iPart = 0; // select Pi+/Pi- only
525       if (TMath::Abs(pdg) == 321)  iPart = 1; // select K+/K- only
526       if (TMath::Abs(pdg) == 2212) iPart = 2; // select p+/p- only
527       if (iPart == -1) continue;
528       //
529
530       if (!fSmallTHnSparse){
531         Double_t vecHistMC[10] = {iPart, centrality,  pT, sign, rap, 0, 1, 0, dxy, 0};
532         if (!fOnlyQA) fHistMCparticles->Fill(vecHistMC);
533       }
534       else{
535         if (rap>fRapidityCutLow && rap<fRapidityCutHigh){
536           Double_t vecHistMC[8] = {iPart, centrality,  pT, sign, 1, 0, dxy, 0};
537           if (!fOnlyQA) fHistMCparticles->Fill(vecHistMC);
538         }
539       }
540     }
541   }
542   //
543   if (!isSelected && !fOnlyQA) {
544     PostData(1, fListHist);
545     return;
546   }
547   //
548   if (!vertex) {
549     fHistMult->Fill(-1, processCode);
550     PostData(1, fListHist);
551     return;
552   } else {
553     if (TMath::Abs(vertex->GetZv()) > 10) {
554       fHistMult->Fill(-1, processCode);
555       PostData(1, fListHist);
556       return;
557     }
558   }
559   //
560   // count events after physics selection and after vertex selection
561   //
562   //cout << "MULTIPLICITY " << trackCounter << " " << fESD->GetEventNumberInFile() <<endl;
563   fHistMult->Fill(trackCounter, processCode);
564   fHistCentrality->Fill(centrality);
565
566
567   //***************************************************
568   // track loop
569   //***************************************************
570   //const Float_t kNsigmaCut = 3;
571   //const Float_t k2sigmaCorr = 1/(0.5*(TMath::Erf(kNsigmaCut/sqrt(2))-TMath::Erf(-kNsigmaCut/sqrt(2))))/*1/0.9545*/;
572   //
573   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
574   //
575   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
576         
577         AliESDtrack *track = 0;
578         AliESDtrack *trackForTOF = 0; //normal track for all TOF information needed when using tpconly-tracks
579         
580         //normal tracks, if tpconly flag is set, use tpconlytracks
581         if (!fUseTPConlyTracks){
582         track =fESD->GetTrack(i); 
583         }
584         else {
585                 track = fESDtrackCuts->GetTPCOnlyTrack(fESD,i);
586                 if (!track) continue;
587                 trackForTOF = fESD->GetTrack(i);
588         }
589     //
590     if (!track->GetInnerParam()) {
591                 if (fUseTPConlyTracks) {delete track; track = 0;} //need to delete tpconlytrack
592                 continue;
593         }
594     Double_t ptot = track->GetInnerParam()->GetP(); // momentum for dEdx determination
595     Double_t pT = track->Pt();
596     track->GetImpactParameters(dca, cov);
597     //
598     //
599     // cut for dead regions in the detector
600     // if (track->Eta() > 0.1 && (track->Eta() < 0.2 && track->Phi() > 0.1 && track->Phi() < 0.1) continue;
601     //
602     // 2.a) apply some standard track cuts according to general recommendations
603     //
604     if (!fESDtrackCuts->AcceptTrack(track)) {
605                 if (fUseTPConlyTracks) {delete track; track = 0;} //need to delete tpconlytrack
606                 continue;
607         }
608
609         UInt_t status = 0;
610     if (!fUseTPConlyTracks) status = track->GetStatus();
611         else status = trackForTOF->GetStatus();
612     Bool_t hasTOFout  = status&AliESDtrack::kTOFout; 
613     Bool_t hasTOFtime = status&AliESDtrack::kTIME;
614     Bool_t hasTOFpid  = status&AliESDtrack::kTOFpid;
615     Bool_t hasTOF     = kFALSE;
616     if (hasTOFout && hasTOFtime && hasTOFpid) hasTOF = kTRUE;
617         Float_t length = 0.;
618     if (!fUseTPConlyTracks) length = track->GetIntegratedLength(); 
619     else length = trackForTOF->GetIntegratedLength();
620
621     if (length < 350.) hasTOF = kFALSE;
622     //
623     // calculate rapidities and kinematics
624     // 
625     //
626     Double_t pvec[3];
627     track->GetPxPyPz(pvec);
628     Double_t energyPion = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kPion)*AliPID::ParticleMass(AliPID::kPion));
629     Double_t energyKaon = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kKaon)*AliPID::ParticleMass(AliPID::kKaon));
630     Double_t energyProton = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton));
631     Double_t energyDeuteron = TMath::Sqrt(track->GetP()*track->GetP() + 4*AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton));
632     //
633     Double_t rapPion = 0.5*TMath::Log((energyPion + pvec[2])/(energyPion - pvec[2]));
634     Double_t rapKaon = 0.5*TMath::Log((energyKaon + pvec[2])/(energyKaon - pvec[2]));
635     Double_t rapProton = 0.5*TMath::Log((energyProton + pvec[2])/(energyProton - pvec[2]));
636     Double_t rapDeuteron = 0.5*TMath::Log((energyDeuteron + pvec[2])/(energyDeuteron - pvec[2]));
637     //
638 //    Double_t transMassPion = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kPion)*AliPID::ParticleMass(AliPID::kPion))      -  AliPID::ParticleMass(AliPID::kPion);
639 //    Double_t transMassKaon = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kKaon)*AliPID::ParticleMass(AliPID::kKaon))     -  AliPID::ParticleMass(AliPID::kKaon);
640  //   Double_t transMassProton = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))    -  AliPID::ParticleMass(AliPID::kProton);
641 //    Double_t transMassDeuteron = TMath::Sqrt(track->Pt()*track->Pt() + 4*AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))    -  2*AliPID::ParticleMass(AliPID::kProton);
642     //
643     // 3. make the PID
644     //
645     Double_t sign = track->GetSign();   
646     Double_t tpcSignal = track->GetTPCsignal();
647     //
648     //
649     // 3.a. calculate expected signals in nsigma
650     //
651     //  
652     // (0.) assumed particle: 0. pion, 1. kaon, 2. proton, 3. deuteron
653     // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
654     // (2.) pT
655     // (3.) sign
656     // (4.) rapidity --> filled 4x
657     // (5.)  pull TPC dEx --> filled 4x
658     // (6.) has valid TOF pid signal
659     // (7.) nsigma TOF --> filled 4x
660     // (8..) dca_xy
661     // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident, 3-second weak, 4-second material
662     //
663 //    Double_t transMass[4] = {transMassPion,transMassKaon,transMassProton,transMassDeuteron};
664     Double_t rap[4] = {rapPion,rapKaon,rapProton,rapDeuteron};
665     Double_t pullsTPC[4] = {fESDpid->NumberOfSigmasTPC(track,AliPID::kPion),
666                             fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon),
667                             fESDpid->NumberOfSigmasTPC(track,AliPID::kProton),
668                             0}; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!
669     Float_t time0 = fESDpid->GetTOFResponse().GetTimeZero();
670     //fESDpid->GetTOFResponse().SetTimeResolution(130.);
671     Double_t pullsTOF[4] ={0.,0.,0.,0.};
672     if (!fUseTPConlyTracks) {
673       pullsTOF[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kPion, time0);
674       pullsTOF[1] = fESDpid->NumberOfSigmasTOF(track,AliPID::kKaon, time0);
675       pullsTOF[2] = fESDpid->NumberOfSigmasTOF(track,AliPID::kProton, time0);
676       pullsTOF[3] = 0; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!;
677     }
678     else {
679       pullsTOF[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kPion, time0);
680       pullsTOF[1] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kKaon, time0);
681       pullsTOF[2] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kProton, time0);
682       pullsTOF[3] = 0; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!;
683     }
684
685     //
686 //    Double_t tpcQA[4] = {fESDpid->NumberOfSigmasTPC(track,AliPID::kElectron),
687 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kPion),
688 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon),
689 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kProton)};
690
691     Double_t tofQA[4] = {0.,0.,0.,0.}; 
692     if (!fUseTPConlyTracks) {
693       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kElectron, time0);
694       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kPion, time0);
695       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kKaon, time0);
696       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kProton, time0);
697     }
698     else{
699       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kElectron, time0);
700       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kPion, time0);
701       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kKaon, time0);
702       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kProton, time0);
703     }
704
705
706     //save information for every particle type  // loop over assumed particle type
707     for(Int_t iPart = 0; iPart < 3; iPart++) {
708
709       if (!fSmallTHnSparse) {
710         Double_t vecHistReal[9]  = {iPart,  centrality,   pT, sign,  rap[iPart], pullsTPC[iPart], hasTOF, pullsTOF[iPart], dca[0]};
711         if (!fOnlyQA) fHistRealTracks->Fill(vecHistReal);
712       }
713       else {
714         if (pullsTPC[iPart]>fTPCnSigmaCutLow && pullsTPC[iPart]<fTPCnSigmaCutHigh && rap[iPart]>fRapidityCutLow && rap[iPart]<fRapidityCutHigh) {
715           Double_t vecHistReal[7]  = {iPart,  centrality,   pT, sign, hasTOF, pullsTOF[iPart], dca[0]};
716           if (!fOnlyQA) fHistRealTracks->Fill(vecHistReal);
717         }
718       }
719
720
721       // using MC truth for precise efficiencies...
722       //
723       if (fMCtrue && !fOnlyQA) {
724         Int_t code = 9; // code: 0-generated, 1-true rec. primaries, 2-misident, 3-second weak, 4-second material
725         Int_t assumedPdg = 0;//2212(proton); 321(Kaon); 211(pion);
726         Int_t motherCode = -1;
727         if (iPart == 0) assumedPdg = 211;
728         if (iPart == 1) assumedPdg = 321;
729         if (iPart == 2) assumedPdg = 2212;
730         //
731         //
732         TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
733         Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
734         //
735         if (pdg != assumedPdg && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 2;
736         if (pdg != assumedPdg && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) code = 5;
737         if (pdg == assumedPdg && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 1;
738         if (pdg == assumedPdg && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) {
739           code = 3;
740           if (fSaveMotherPDG){
741             TParticle *trackMother =  stack->Particle(TMath::Abs(trackMC->GetFirstMother()));
742             if (trackMother->GetPdgCode() == 310) motherCode = 6; //K0
743             if (trackMother->GetPdgCode() == 3122) motherCode = 7; //Lambda
744             if (trackMother->GetPdgCode() == 3222) motherCode = 8; //Sigma+
745           }
746         }
747         if (pdg == assumedPdg && stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()))) code = 4;
748         
749         //
750         // muons need special treatment, because they are indistinguishable from pions
751         //
752         if (iPart == 0 && pdg == 13  && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 1;
753         if (iPart == 0 && pdg == 13  && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) code = 3;
754         //
755         // check TOF mismatch on MC basis with TOF label
756         //
757         Int_t tofLabel[3];
758         if (!fUseTPConlyTracks) track->GetTOFLabel(tofLabel);
759         else trackForTOF->GetTOFLabel(tofLabel);
760         if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0])) hasTOF = kFALSE;
761         //
762         // IMPORTANT BIG PROBLEM HERE THE PROBABLILITY TO HAVE A PID SIGNAL MUST BE IN !!!!!!!!!!!!
763         //
764         if (!fSmallTHnSparse){
765           Double_t vectorHistMC[10] = {iPart,  centrality,  pT, sign,  rap[iPart], pullsTPC[iPart], hasTOF, pullsTOF[iPart], dca[0], code};
766           if (!fOnlyQA) { 
767             fHistMCparticles->Fill(vectorHistMC);
768             if (motherCode != -1 && fSaveMotherPDG) { //if mother of weak decay is K0, lambda or sigma+ add track again with this information
769               Double_t vectorHistMCmother[10] = {iPart,  centrality,  pT, sign,  rap[iPart], pullsTPC[iPart], hasTOF, pullsTOF[iPart], dca[0], motherCode};
770               fHistMCparticles->Fill(vectorHistMCmother);
771             }
772           }
773         }
774         else{
775           if (pullsTPC[iPart]>fTPCnSigmaCutLow && pullsTPC[iPart]<fTPCnSigmaCutHigh && rap[iPart]>fRapidityCutLow && rap[iPart]<fRapidityCutHigh) {
776             //                              0,           1,   2,    3,           4,               5,      6,               7,      8,   9
777             Double_t vectorHistMC[8] = {iPart,  centrality,  pT, sign, hasTOF, pullsTOF[iPart], dca[0], code};
778             if (!fOnlyQA) { 
779               fHistMCparticles->Fill(vectorHistMC);
780               if (motherCode != -1 && fSaveMotherPDG) { //if mother of weak decay is K0, lambda or sigma+ add track again with this information
781                 Double_t vectorHistMCmother[8] = {iPart,  centrality,  pT, sign, hasTOF, pullsTOF[iPart], dca[0], motherCode};
782                 fHistMCparticles->Fill(vectorHistMCmother);
783               }
784             }
785           }
786         }
787       }
788       //
789       //
790       Int_t tpcShared = track->GetTPCnclsS();
791       if (TMath::Abs(track->Eta()) < 0.8 && iPart == 0 && tpcShared < 4) fHistPidQA->Fill(ptot,tpcSignal,sign);
792     } // end loop over assumed particle type
793
794           //need to delete tpconlytrack
795           if (fUseTPConlyTracks){
796                 delete track; 
797                 track = 0;     
798           }
799
800   } // end of track loop
801   
802   // Post output data  
803   PostData(1, fListHist);
804   
805 }      
806
807
808 //________________________________________________________________________
809 void AliAnalysisCombinedHadronSpectra::Terminate(Option_t *) 
810 {
811   // Draw result to the screen
812   // Called once at the end of the query
813   Printf("*** CONSTRUCTOR CALLED ****");
814
815 }
816
817
818 //________________________________________________________________________
819 Bool_t AliAnalysisCombinedHadronSpectra::SelectOnImpPar(AliESDtrack* t) {
820   //
821   // cut on transverse impact parameter // DEPRECATED
822   //
823   Float_t d0z0[2],covd0z0[3];
824   t->GetImpactParameters(d0z0,covd0z0);
825   Float_t sigma= 0.0050+0.0060/TMath::Power(t->Pt(),0.9);
826   Float_t d0max = 7.*sigma;
827   //
828   Float_t sigmaZ = 0.0146+0.0070/TMath::Power(t->Pt(),1.114758);
829   if (t->Pt() > 1) sigmaZ = 0.0216;
830   Float_t d0maxZ = 5.*sigmaZ;
831   //
832   if(TMath::Abs(d0z0[0]) < d0max && TMath::Abs(d0z0[1]) < d0maxZ) return kTRUE;
833   return kFALSE;
834 }
835
836
837 //________________________________________________________________________
838 void AliAnalysisCombinedHadronSpectra::BinLogAxis(const TH1 *h) {
839   //
840   // Method for the correct logarithmic binning of histograms
841   //
842   TAxis *axis = h->GetXaxis();
843   int bins = axis->GetNbins();
844
845   Double_t from = axis->GetXmin();
846   Double_t to = axis->GetXmax();
847   Double_t *newBins = new Double_t[bins + 1];
848    
849   newBins[0] = from;
850   Double_t factor = pow(to/from, 1./bins);
851   
852   for (int i = 1; i <= bins; i++) {
853    newBins[i] = factor * newBins[i-1];
854   }
855   axis->Set(bins, newBins);
856   delete [] newBins;
857   
858 }
859
860
861 //________________________________________________________________________
862 Int_t AliAnalysisCombinedHadronSpectra::GetPythiaEventProcessType(const AliHeader* aHeader, const Bool_t adebug) const {
863   //
864   // get the process type of the event.
865   //
866
867   // can only read pythia headers, either directly or from cocktalil header
868   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
869
870   if (!pythiaGenHeader) {
871
872     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
873     if (!genCocktailHeader) {
874       //printf("AliAnalysisCombinedHadronSpectra::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
875       return -1;
876     }
877
878     TList* headerList = genCocktailHeader->GetHeaders();
879     if (!headerList) {
880       return -1;
881     }
882
883     for (Int_t i=0; i<headerList->GetEntries(); i++) {
884       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
885       if (pythiaGenHeader)
886         break;
887     }
888
889     if (!pythiaGenHeader) {
890       //printf("AliAnalysisCombinedHadronSpectra::GetProcessType : Could not find Pythia header. \n");
891       return -1;
892     }
893   }
894
895   if (adebug) {
896     //printf("AliAnalysisCombinedHadronSpectra::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType());
897   }
898
899   return pythiaGenHeader->ProcessType();
900 }