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