ce877a946f90060186627043ad00d0ecaa21a289
[u/mrichter/AliRoot.git] / PWGLF / SPECTRA / Nuclei / deuteronpA / AliAnalysisDeuteronpA.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 "TH2F.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 "AliAnalysisDeuteronpA.h"
62
63 #include "AliAnalysisUtils.h"
64
65
66 ClassImp(AliAnalysisDeuteronpA)
67
68 //________________________________________________________________________
69 AliAnalysisDeuteronpA::AliAnalysisDeuteronpA() 
70 : AliAnalysisTaskSE("TaskChargedHadron"), fESD(0), fListHist(0), fESDtrackCuts(0),fESDTrackCutsMult(0),fESDpid(0),
71   fMCtrue(0),
72   fRapCMSpA(0),
73   fAlephParameters(),
74   fUtils(0),
75   fHistRealTracks(0),
76   fHistMCparticles(0),
77   fHistPidQA(0),
78   fHistTofQA(0),
79   fHistMult(0),
80   fHistCentrality(0),
81   fHistMomCorr(0),
82   fHistEtaPtGen(0),
83   fHistVertex(0),
84   fHistVertexRes(0),
85   fHistVertexResTracks(0)
86 {
87   // default Constructor
88   /* fast compilation test
89      gSystem->Load("libANALYSIS");
90      gSystem->Load("libANALYSISalice");
91      .L AliAnalysisDeuteronpA.cxx++
92    */
93 }
94
95
96 //________________________________________________________________________
97 AliAnalysisDeuteronpA::AliAnalysisDeuteronpA(const char *name) 
98   : AliAnalysisTaskSE(name), fESD(0), fListHist(0), fESDtrackCuts(0),fESDTrackCutsMult(0),fESDpid(0),
99     fMCtrue(0),
100     fRapCMSpA(0),
101     fAlephParameters(),
102     fUtils(0),
103     fHistRealTracks(0),
104     fHistMCparticles(0),
105     fHistPidQA(0),
106     fHistTofQA(0),
107     fHistMult(0),
108     fHistCentrality(0),
109     fHistMomCorr(0),
110     fHistEtaPtGen(0),
111     fHistVertex(0),
112     fHistVertexRes(0),
113     fHistVertexResTracks(0)
114 {
115   //
116   // standard constructur which should be used
117   //
118   Printf("*** CONSTRUCTOR CALLED ****");
119   //
120   fMCtrue = kTRUE; 
121   fRapCMSpA = kTRUE;
122   /* real */
123   fAlephParameters[0] = 0.0283086;
124   fAlephParameters[1] = 2.63394e+01;
125   fAlephParameters[2] = 5.04114e-11;
126   fAlephParameters[3] = 2.12543e+00;
127   fAlephParameters[4] = 4.88663e+00;
128   //
129   // initialize PID object
130   //
131   //fESDpid = new AliESDpid();
132   //
133   // create track cuts
134   //
135   fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
136   //
137   Initialize();
138   // Output slot #0 writes into a TList container
139   DefineOutput(1, TList::Class());
140
141 }
142
143
144 //________________________________________________________________________
145 void AliAnalysisDeuteronpA::Initialize()
146 {
147   //
148   // updating parameters in case of changes
149   //
150
151
152   fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE,kTRUE);
153   fESDtrackCuts->SetMaxDCAToVertexXY(3);
154   fESDtrackCuts->SetMaxDCAToVertexZ(2);
155   fESDtrackCuts->SetEtaRange(-0.9,0.9);
156   //
157   //
158
159 }
160
161
162 //________________________________________________________________________
163 void AliAnalysisDeuteronpA::UserCreateOutputObjects() 
164 {
165
166   //Create analysis utils for event selection and pileup rejection
167   fUtils = new AliAnalysisUtils();
168   fUtils->SetCutOnZVertexSPD(kFALSE);
169   //
170   //
171   // Create histograms
172   // Called once
173   fListHist = new TList();
174   fListHist->SetOwner(kTRUE);
175   //
176   const Int_t kPtBins = 28;
177   const Int_t kMultBins = 11;
178   const Int_t kDcaBins = 38;
179
180
181   //different binning for YCMS in pA to cover detector acceptance
182   Double_t kYlowBorder = -0.6;
183   Double_t kYhighBorder = 0.6;
184   Double_t kYBins = 12;
185
186   if (fRapCMSpA){
187     kYlowBorder = -0.1;
188     kYhighBorder = 1.1;
189     kYBins = 12;
190   }
191
192   //
193   // sort pT-bins ..
194   //
195   Double_t binsPt[kPtBins+1] = {0., 0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0, 5.2, 5.4, 5.6};
196   //
197   // provide a reasonable dca/binning
198   //
199   Double_t binsDca[kDcaBins+1];
200   for(Int_t i = 0; i< kDcaBins+1; i++) {
201     binsDca[i] = 1.5*(1./TMath::Exp(-TMath::Abs((i - kDcaBins/2.))/(kDcaBins/2)) -1);
202     if ((i - kDcaBins/2.) < 0) binsDca[i] *= -1;
203     //cout <<  binsDca[i] << endl;
204   }
205   //
206   // create the histograms with all necessary information --> it is filled 4x for each particle assumption
207   //
208   // (0.) assumed particle: 0. deuteron, 1. triton, 2. He-3
209   // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
210   // (2.) pT
211   // (3.) sign
212   // (4.) rapidity --> filled 4x
213   // (5.)  pull TPC dEx --> filled 4x
214   // (6.) has valid TOF pid signal
215   // (7.) nsigma TOF --> filled 4x XXXXXXXXX no mass*mass
216   // (8..) dca_xy
217   // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident prim, 3-second weak, 4-second material, 5-misident sec
218   //
219   //                              0,           1,           2,  3,      4,              5,    6,    7,       8
220   Int_t    binsHistReal[9] = {   3,   kMultBins,     kPtBins,  2,      static_cast<Int_t>(kYBins)      ,   50,    2,  100, kDcaBins};
221   Double_t xminHistReal[9] = {-0.5,        -0.5,           0, -2,      kYlowBorder ,   -5,- 0.5, -2.5,       -3};
222   Double_t xmaxHistReal[9] = { 2.5,        10.5,           3,  2,      kYhighBorder,    5,  1.5,  2.5,        3};
223   fHistRealTracks = new THnSparseF("fHistRealTracks","real tracks",9,binsHistReal,xminHistReal,xmaxHistReal);
224   //
225   fHistRealTracks->GetAxis(2)->Set(kPtBins, binsPt);
226   fHistRealTracks->GetAxis(8)->Set(kDcaBins, binsDca);
227   fListHist->Add(fHistRealTracks);
228   //
229   //
230   fHistPidQA = new TH3F("fHistPidQA","PID QA",500,0.1,10,1000,0,1000,2,-2,2);
231   BinLogAxis(fHistPidQA);
232   fListHist->Add(fHistPidQA);
233   //
234   fHistTofQA = new TH2F("fHistTofQA","TOF-QA",200,-4.,4.,400,0,1.1);
235   fListHist->Add(fHistTofQA);
236   //                            0,            1,           2,  3,      4,            5,    6,    7,        8,    9
237   Int_t    binsHistMC[10] = {   3,    kMultBins,     kPtBins,  2,    static_cast<Int_t>(kYBins)      ,  50,    2,  100, kDcaBins,    6};
238   Double_t xminHistMC[10] = {-0.5,         -0.5,           0, -2,    kYlowBorder ,  -5,- 0.5, -2.5,       -3, -0.5};
239   Double_t xmaxHistMC[10] = { 2.5,         10.5,           3,  2,    kYhighBorder,   5,  1.5,  2.5,        3,  5.5};
240   //
241   // different binning for CODE axis, if we want to save motherPDG
242   //
243   fHistMCparticles = new THnSparseF("fHistMCparticles","MC histogram",10,binsHistMC,xminHistMC,xmaxHistMC);
244   fHistMCparticles->GetAxis(2)->Set(kPtBins, binsPt);
245   fHistMCparticles->GetAxis(8)->Set(kDcaBins, binsDca);
246   fListHist->Add(fHistMCparticles);
247   //
248   fHistMult = new TH2F("fHistMult", "control histogram to count number of events", 502, -2.5, 499.5,4,-0.5,3.5);
249   fHistCentrality = new TH1F("fHistCentrality", "control histogram to count number of events", 22, -1.5, 20.5);
250   fListHist->Add(fHistMult);
251   fListHist->Add(fHistCentrality);
252   //
253   fHistMomCorr = new TH2F("fHistMomCorr","momentum correction; pT_{MCtrue}; rel. difference",200,0,3,200,-0.2,0.2);
254   fListHist->Add(fHistMomCorr);
255   //
256   fHistEtaPtGen = new TH3F("fHistPtEtaGen","rapidity correction; p_{T} (GeV/c); rapidity y; pid",200,0,12, 200,-1.2,1.2, 3,-0.5,2.5);
257   fListHist->Add(fHistEtaPtGen);
258   //
259   fHistVertex = new TH2F("fHistVertex","vertex position; V_{XY} (cm); V_{Z}",2000,-1.,1.,400,-20.,20.);
260   fListHist->Add(fHistVertex);
261   //
262   fHistVertexRes = new TH1F("fHistVertexRes","vertex position difference MC truth and rec; V_{Z,rec} - V_{Z,truth} (cm)",2000,-0.5,0.5);
263   fListHist->Add(fHistVertexRes);
264   //
265   fHistVertexResTracks = new TH2F("fHistVertexResTracks","vertex res vs no of tracks; no. of ch. primaries |#eta| < 0.9; V_{Z,rec} - V_{Z,truth} (cm)",50,-0.5,499.5,200,-0.05,0.05);
266   fListHist->Add(fHistVertexResTracks);
267
268   PostData(1, fListHist);
269
270 }
271
272 //________________________________________________________________________
273 void AliAnalysisDeuteronpA::UserExec(Option_t *) 
274 {
275   //
276   // main event loop
277   //
278   if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
279   if (!fESDpid) {
280     fESDpid = new AliESDpid(); // HACK FOR MC PBPB --> PLEASE REMOVE AS SOON AS POSSIBLE
281     fESDpid->GetTPCResponse().SetBetheBlochParameters(1.28778e+00/50., 3.13539e+01, TMath::Exp(-3.16327e+01), 1.87901e+00, 6.41583e+00);
282   }
283   //AliLog::SetGlobalLogLevel(AliLog::kError);
284   //
285   // Check Monte Carlo information and other access first:
286   //
287   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
288   if (!eventHandler) {
289     //Printf("ERROR: Could not retrieve MC event handler");
290     fMCtrue = kFALSE;
291   }
292   //
293   AliMCEvent* mcEvent = 0x0;
294   AliStack* stack = 0x0;
295   if (eventHandler) mcEvent = eventHandler->MCEvent();
296   if (!mcEvent) {
297     //Printf("ERROR: Could not retrieve MC event");
298     if (fMCtrue) return;
299   }
300   if (fMCtrue) {
301     stack = mcEvent->Stack();
302     if (!stack) return;
303   }
304   //
305   fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
306   if (!fESD) {
307     //Printf("ERROR: fESD not available");
308     return;
309   }
310   
311   //check with analysis utils
312   //first event in chunck --> continue
313   if(fUtils->IsFirstEventInChunk(fESD)) {
314     PostData(1, fListHist);
315     return;
316   }
317   //check for pileup
318   if (fUtils->IsPileUpEvent(fESD)){
319     PostData(1, fListHist);
320     return;
321   }
322   //vetex cuts
323   Bool_t isVtxOk = kTRUE;
324   if(!fUtils->IsVertexSelected2013pA(fESD)) isVtxOk = kFALSE;
325   //NEED TO PUT THIS IN ACTION SOMEWHERE
326
327
328   if (!fESDtrackCuts) {
329     Printf("ERROR: fESDtrackCuts not available");
330     return;
331   }
332   //
333   // check if event is selected by physics selection class
334   //
335   //
336   // monitor vertex position
337   //
338   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
339   if(vertex->GetNContributors()<1) {
340     // SPD vertex
341     vertex = fESD->GetPrimaryVertexSPD();
342     if(vertex->GetNContributors()<1) vertex = 0x0;
343   }  
344   //
345   // fill histogram to monitor vertex position
346   if (vertex && isVtxOk) fHistVertex->Fill(TMath::Sqrt(vertex->GetXv()*vertex->GetXv() + vertex->GetYv()*vertex->GetYv()),vertex->GetZv());
347
348   Float_t vertexRes = -100;
349   Int_t nPrimaries = 0;
350
351
352   if (fMCtrue && vertex && isVtxOk) {
353     //
354     //
355     for(Int_t i = 0; i < stack->GetNtrack(); i++) {
356       TParticle * trackMC = stack->Particle(i);
357       
358       if (!trackMC->IsPrimary()) continue;
359       //
360       //Double_t xv = trackMC->Vx();
361       //Double_t yv = trackMC->Vy();
362       Double_t zv = trackMC->Vz();
363
364       vertexRes = vertex->GetZv() - zv;
365       fHistVertexRes->Fill(vertex->GetZv() - zv);
366
367       break;
368     }
369
370     for(Int_t i = 0; i < stack->GetNtrack(); i++) {
371       TParticle * trackMC = stack->Particle(i);
372       
373       //      if (trackMC->IsPrimary()) {
374       if (stack->IsPhysicalPrimary(i)) {
375         if (trackMC->Eta() > -0.9 && trackMC->Eta() < 0.9 && trackMC->GetPDG()->Charge() != 0){
376           nPrimaries++;
377           //cout << nPrimaries << endl;
378         }
379       }
380     }
381   }
382
383   fHistVertexResTracks->Fill(nPrimaries, vertexRes);
384
385
386
387
388   //
389   Float_t centrality = -1;
390   //
391   /*
392   Int_t rootS = fESD->GetBeamEnergy() < 1000 ? 0 : 1;
393   if (fESD->GetEventSpecie() == 4) { // PbPb
394     rootS = 2;
395     AliCentrality *esdCentrality = fESD->GetCentrality();
396     centrality = esdCentrality->GetCentralityClass10("V0M") + 1; // centrality percentile determined with V0
397     if (TMath::Abs(centrality - 1) < 1e-5) {
398       centrality = esdCentrality->GetCentralityClass5("V0M");
399     }
400   }
401   */
402   AliCentrality *esdCentrality = fESD->GetCentrality();
403   centrality = esdCentrality->GetCentralityClass10("V0A") + 1; // centrality percentile determined with V0A
404   if (TMath::Abs(centrality - 1) < 1e-5) {
405     centrality = esdCentrality->GetCentralityClass5("V0A");
406   }
407
408
409   //
410   Int_t processCode = 0;
411   //
412   // important change: fill generated only after vertex cut in case of heavy-ions
413   //
414   if ( fESD->GetEventSpecie() == 4) {
415     if (!vertex || !isVtxOk) {
416       fHistMult->Fill(-1, processCode);
417       PostData(1, fListHist);
418       return;
419     } else {
420       if (TMath::Abs(vertex->GetZv()) > 10) {
421         fHistMult->Fill(-1, processCode);
422         PostData(1, fListHist);
423         return;
424       }
425     }
426   }
427   //
428   if (fMCtrue) {
429     //
430     //
431     for(Int_t i = 0; i < stack->GetNtrack(); i++) {
432       TParticle * trackMC = stack->Particle(i);
433       Int_t pdg = trackMC->GetPdgCode();
434       //
435       Double_t xv = trackMC->Vx();
436       Double_t yv = trackMC->Vy();
437       //Double_t zv = trackMC->Vz();
438       Double_t dxy = 0;
439       dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
440       //Double_t dz = 0;
441       //dz = TMath::Abs(zv); // so stupid to avoid warnings
442       //
443       // vertex cut - selection of primaries
444       //
445       //if (dxy > 3 || dz > 10) continue; // fixed cut at 3cm in r
446       //
447       if (!stack->IsPhysicalPrimary(i)) continue;
448       //
449       // fill MC histograms here...
450       // 
451       Double_t rap = trackMC->Y();
452       if (fRapCMSpA) rap = rap + 0.465;
453       Double_t pT  = trackMC->Pt();
454       Int_t sign = pdg < 0 ? -1 : 1; // only works for charged pi,K,p !!
455       //      Double_t transMass = TMath::Sqrt(trackMC->Pt()*trackMC->Pt() + trackMC->GetMass()*trackMC->GetMass()) - trackMC->GetMass();
456       //
457       Int_t iPart = -1;
458       if (TMath::Abs(pdg) == 1000010020)  iPart = 0; // select d+/d- only
459       if (TMath::Abs(pdg) == 1000010030)  iPart = 1; // select t+/t- only
460       if (TMath::Abs(pdg) == 1000020030)  iPart = 2; // select He+/He- only
461       if (iPart == -1) continue;
462       //
463       Double_t vecHistMC[10] = {static_cast<Double_t>(iPart), centrality,  pT, static_cast<Double_t>(sign), rap, 0, 1, 0, dxy, 0};
464       fHistMCparticles->Fill(vecHistMC);
465       //
466       if (iPart == 0) fHistEtaPtGen->Fill(pT, trackMC->Y(), 0.);
467       if (iPart == 1) fHistEtaPtGen->Fill(pT, trackMC->Y(), 1.);
468       if (iPart == 2) fHistEtaPtGen->Fill(pT, trackMC->Y(), 2.);
469     }
470   }
471   //
472   // check vertex
473   //
474   if (!vertex) {
475     fHistMult->Fill(-1, processCode);
476     PostData(1, fListHist);
477     return;
478   } else {
479     if (TMath::Abs(vertex->GetZv()) > 10) {
480       fHistMult->Fill(-1, processCode);
481       PostData(1, fListHist);
482       return;
483     }
484   }
485   //
486   // count events after physics selection and after vertex selection
487   //
488   fHistMult->Fill(fESDtrackCuts->CountAcceptedTracks(fESD), 0);
489   fHistMult->Fill(fESD->GetNumberOfTracks(), 1);
490   fHistCentrality->Fill(centrality);
491   //
492   //***************************************************
493   // track loop
494   //***************************************************
495   //const Float_t kNsigmaCut = 3;
496   //const Float_t k2sigmaCorr = 1/(0.5*(TMath::Erf(kNsigmaCut/sqrt(2))-TMath::Erf(-kNsigmaCut/sqrt(2))))/*1/0.9545*/;
497   //
498   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
499   //
500   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
501     //
502     AliESDtrack *track  =fESD->GetTrack(i); 
503     if (!track->GetInnerParam()) continue;
504     //
505     Double_t ptot = track->GetInnerParam()->GetP(); // momentum for dEdx determination
506     //
507     // momentum correction for different mass assumption in tracking
508     //
509     Double_t pT = track->Pt()/(1 - 0.333303/TMath::Power(track->Pt() + 0.651111, 5.27268));
510     // -[0]/TMath::Power(x - [1],[2])
511     //
512     track->GetImpactParameters(dca, cov);
513     //
514     // cut for dead regions in the detector
515     // if (track->Eta() > 0.1 && (track->Eta() < 0.2 && track->Phi() > 0.1 && track->Phi() < 0.1) continue;
516     //
517     // 2.a) apply some standard track cuts according to general recommendations
518     //
519     if (!fESDtrackCuts->AcceptTrack(track)) continue;
520     //
521     UInt_t status = track->GetStatus();
522     //
523     Bool_t hasTOFout  = status&AliESDtrack::kTOFout; 
524     Bool_t hasTOFtime = status&AliESDtrack::kTIME;
525     Bool_t hasTOF     = kFALSE;
526     if (hasTOFout) hasTOF = kTRUE;
527     Float_t length = track->GetIntegratedLength(); 
528     if (length < 350.) hasTOF = kFALSE;
529     //
530     // calculate rapidities and kinematics
531     // 
532     //
533     Double_t pvec[3];
534     track->GetPxPyPz(pvec);
535     Double_t energyDeuteron = TMath::Sqrt(pT*pT + pvec[2]*pvec[2] +
536                                           AliPID::ParticleMass(AliPID::kDeuteron)*AliPID::ParticleMass(AliPID::kDeuteron));
537     Double_t energyTriton = TMath::Sqrt(track->GetP()*track->GetP() + 
538                                         AliPID::ParticleMass(AliPID::kTriton)*AliPID::ParticleMass(AliPID::kTriton));
539     Double_t energyHe3 = TMath::Sqrt(track->GetP()*2*track->GetP()*2 + 
540                                           AliPID::ParticleMass(AliPID::kHe3)*AliPID::ParticleMass(AliPID::kHe3));
541     //
542     Double_t rapDeuteron = 0.5*TMath::Log((energyDeuteron + pvec[2])/(energyDeuteron - pvec[2]));
543     Double_t rapTriton = 0.5*TMath::Log((energyTriton + pvec[2])/(energyTriton - pvec[2]));
544     Double_t rapHe3 = 0.5*TMath::Log((energyHe3 + pvec[2]*2)/(energyHe3 - pvec[2]*2));
545     //
546     if (fRapCMSpA) {
547       rapDeuteron = rapDeuteron + 0.465;
548       rapTriton = rapTriton + 0.465;
549       rapHe3 = rapHe3 + 0.465;
550     }
551     //
552     //
553     // 3. make the PID
554     //
555     Double_t sign = track->GetSign();   
556     Double_t tpcSignal = track->GetTPCsignal();
557     //
558     //
559     // 3.a. calculate expected signals in nsigma
560     //
561     //  
562     // (0.) assumed particle: 0. pion, 1. kaon, 2. proton, 3. deuteron
563     // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
564     // (2.) pT
565     // (3.) sign
566     // (4.) rapidity --> filled 4x
567     // (5.)  pull TPC dEx --> filled 4x
568     // (6.) has valid TOF pid signal
569     // (7.) nsigma TOF --> filled 4x
570     // (8..) dca_xy
571     // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident, 3-second weak, 4-second material
572     //
573     //
574     Float_t deutExp = AliExternalTrackParam::BetheBlochAleph(ptot/(0.938*2),1.45802,27.4992,4.00313e-15,2.48485,8.31768);
575     Float_t tritExp = AliExternalTrackParam::BetheBlochAleph(ptot/(0.938*3),1.45802,27.4992,4.00313e-15,2.48485,8.31768);
576     Float_t hel3Exp = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/(0.938*3),1.74962,27.4992,4.00313e-15,2.42485,8.31768);
577     if (fMCtrue && ptot > 0.3) {
578       //Double_t parMC[5] = {1.17329, 27.4992, 4.00313e-15, 2.35563, 9.47569}; // OLD FOR LHC11b9_1 !!
579       //Double_t parMC[5] = {1.17329, 27.4992, 4.00313e-15, 2.1204316, 4.1373729};  // NEW!!! PbPb
580       Double_t parMC[5] = {20.1533, 2.58127, 0.00114169, 2.0373, 0.502123}; //for LHC13b2efix enhanced
581       deutExp = AliExternalTrackParam::BetheBlochAleph(ptot/(0.938*2),parMC[0],parMC[1],parMC[2],parMC[3],parMC[4]);
582       tritExp = AliExternalTrackParam::BetheBlochAleph(ptot/(0.938*3),parMC[0],parMC[1],parMC[2],parMC[3],parMC[4]);
583       hel3Exp = 4*AliExternalTrackParam::BetheBlochAleph(2*ptot/(0.938*3),parMC[0],parMC[1],parMC[2],parMC[3],parMC[4]);
584     }
585     //
586     //
587     Double_t rap[4] = {rapDeuteron, rapTriton, rapHe3,0};
588     Double_t pullsTPC[4] = {(tpcSignal - deutExp)/(0.07*deutExp),
589                             (tpcSignal - tritExp)/(0.07*tritExp),
590                             (tpcSignal - hel3Exp)/(0.07*hel3Exp),
591                             0};
592     //
593     // Process TOF information
594     //
595     //Float_t time0 = fESDpid->GetTOFResponse().GetTimeZero();
596     Float_t time0 = fESDpid->GetTOFResponse().GetStartTime(track->P());//fESDpid->GetTOFResponse().GetTimeZero();
597     Float_t mass = 0;
598     Float_t time = -1; 
599     Float_t beta = 0;
600     //
601     if (length > 350.) {
602       time = track->GetTOFsignal() - time0;
603       if (time > 0) {
604         beta = length / (2.99792457999999984e-02 * time);
605         Float_t gamma = 1/TMath::Sqrt(1 - beta*beta);
606         mass = ptot/TMath::Sqrt(gamma*gamma - 1); // using inner TPC mom. as approx.
607       }
608     }
609     //
610     /*
611     Double_t timeDeuteron = (length/2.99792457999999984e-02) * TMath::Sqrt(1 + 1.876*1.876/(ptot*ptot));
612     Double_t timeTriton   = (length/2.99792457999999984e-02) * TMath::Sqrt(1 + 2.809*2.809/(ptot*ptot));
613     Double_t timeHe3      = (length/2.99792457999999984e-02) * TMath::Sqrt(1 + 2.809*2.809/(2*ptot*2*ptot));
614     */
615     //
616     Double_t pullsTOF[4] ={0.,0.,0.,0.};
617     pullsTOF[0] = mass*mass - 1.876*1.876; // assuming 130.ps time resolution
618     pullsTOF[1] = mass*mass - 2.809*2.809;
619     pullsTOF[2] = mass*mass - 2.809*2.809;
620     pullsTOF[3] = 0;
621     //
622     //
623     if (hasTOFout && hasTOFtime && TMath::Abs(pullsTPC[1]) < 3) fHistTofQA->Fill(ptot*sign, beta);
624     //
625     //
626     for(Int_t iPart = 0; iPart < 3; iPart++) { // loop over assumed particle type
627       //
628       // temporary measure to reduce memory: fill only deuterons
629       //
630       if (iPart != 0) continue;
631       //                              0,           1,    2,    3,           4,               5,      6,              7,     8
632       Double_t vecHistReal[9]  = {static_cast<Double_t>(iPart),  centrality,   pT, static_cast<Double_t>(sign),  rap[iPart], pullsTPC[iPart], static_cast<Double_t>(hasTOF), pullsTOF[iPart], dca[0]};
633       fHistRealTracks->Fill(vecHistReal);
634       if (TMath::Abs(rap[1]) < 0.5) fHistPidQA->Fill(ptot,tpcSignal,sign); // has to be used for the very difficult deuterons.
635       //
636       // using MC truth for precise efficiencies...
637       //
638       if (fMCtrue) {
639         Int_t code = 9; // code: 0-generated, 1-true rec. primaries, 2-misident, 3-second weak, 4-second material
640         Int_t assumedPdg = 0;//2212(proton); 321(Kaon); 211(pion);
641         //
642         if (iPart == 0) assumedPdg = 1000010020;
643         if (iPart == 1) assumedPdg = 1000010030;
644         if (iPart == 2) assumedPdg = 1000020030;
645         //
646         //
647         TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
648         Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
649         //
650         if (pdg != assumedPdg && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 2;
651         if (pdg != assumedPdg && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) code = 5;
652         if (pdg == assumedPdg && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 1;
653         if (pdg == assumedPdg && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) {
654           code = 3;
655         }
656         if (pdg == assumedPdg && stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()))) code = 4;
657         //
658         // fill momentum correction histogram
659         //
660         if (code == 1) fHistMomCorr->Fill(trackMC->Pt(), (pT - trackMC->Pt())/trackMC->Pt());
661         //
662         //here cout rap[ipart] and trackMC->Y() and pT
663         //if (trackMC->Pt()>1.5 && TMath::Abs(rap[iPart] - trackMC->Y()) > 0.1) {
664         //    printf("Y_rec: %4.3f    Y_MC: %4.3f    p_T: %5.3f \n", rap[iPart], trackMC->Y(), trackMC->Pt());
665         //}
666          //
667         // check TOF mismatch on MC basis with TOF label
668         //
669         Int_t tofLabel[3] = {0,0,0};
670         track->GetTOFLabel(tofLabel);
671         //      
672         //
673         if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) {
674           hasTOF = kFALSE;
675           if (TMath::Abs(trackMC->GetFirstMother()) == TMath::Abs(tofLabel[0])) hasTOF = kTRUE;
676         }
677         //
678         // IMPORTANT BIG PROBLEM HERE THE PROBABLILITY TO HAVE A PID SIGNAL MUST BE IN !!!!!!!!!!!!
679         //
680         //                              0,           1,   2,    3,           4,               5,      6,               7,      8,   9
681         Double_t vectorHistMC[10] = {static_cast<Double_t>(iPart),  centrality,  pT, static_cast<Double_t>(sign),  rap[iPart], pullsTPC[iPart], static_cast<Double_t>(hasTOF), pullsTOF[iPart], dca[0], static_cast<Double_t>(code)};
682         fHistMCparticles->Fill(vectorHistMC);
683
684       }
685       //
686       //
687     } // end loop over assumed particle type
688
689
690   } // end of track loop
691   
692   // Post output data  
693   PostData(1, fListHist);
694   
695 }      
696
697
698 //________________________________________________________________________
699 void AliAnalysisDeuteronpA::Terminate(Option_t *) 
700 {
701   // Draw result to the screen
702   // Called once at the end of the query
703   Printf("*** CONSTRUCTOR CALLED ****");
704
705 }
706
707
708 //________________________________________________________________________
709 void AliAnalysisDeuteronpA::BinLogAxis(const TH1 *h) {
710   //
711   // Method for the correct logarithmic binning of histograms
712   //
713   TAxis *axis = h->GetXaxis();
714   int bins = axis->GetNbins();
715
716   Double_t from = axis->GetXmin();
717   Double_t to = axis->GetXmax();
718   Double_t *newBins = new Double_t[bins + 1];
719    
720   newBins[0] = from;
721   Double_t factor = pow(to/from, 1./bins);
722   
723   for (int i = 1; i <= bins; i++) {
724    newBins[i] = factor * newBins[i-1];
725   }
726   axis->Set(bins, newBins);
727   delete [] newBins;
728   
729 }
730