1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 ///////////////////////////////////////////////////////////////////////////
19 // Analysis for identified charged hadron spectra. //
22 ///////////////////////////////////////////////////////////////////////////
24 #include "Riostream.h"
33 #include "TObjArray.h"
37 #include "AliAnalysisTaskSE.h"
38 #include "AliAnalysisManager.h"
40 #include "AliHeader.h"
41 #include "AliGenPythiaEventHeader.h"
42 #include "AliGenCocktailEventHeader.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"
55 #include "AliMCEventHandler.h"
56 #include "AliMCEvent.h"
61 #include "AliAnalysisDeuteronpA.h"
63 #include "AliAnalysisUtils.h"
66 ClassImp(AliAnalysisDeuteronpA)
68 //________________________________________________________________________
69 AliAnalysisDeuteronpA::AliAnalysisDeuteronpA()
70 : AliAnalysisTaskSE("TaskChargedHadron"), fESD(0), fListHist(0), fESDtrackCuts(0),fESDTrackCutsMult(0),fESDpid(0),
85 fHistVertexResTracks(0)
87 // default Constructor
88 /* fast compilation test
89 gSystem->Load("libANALYSIS");
90 gSystem->Load("libANALYSISalice");
91 .L AliAnalysisDeuteronpA.cxx++
96 //________________________________________________________________________
97 AliAnalysisDeuteronpA::AliAnalysisDeuteronpA(const char *name)
98 : AliAnalysisTaskSE(name), fESD(0), fListHist(0), fESDtrackCuts(0),fESDTrackCutsMult(0),fESDpid(0),
113 fHistVertexResTracks(0)
116 // standard constructur which should be used
118 Printf("*** CONSTRUCTOR CALLED ****");
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;
129 // initialize PID object
131 //fESDpid = new AliESDpid();
135 fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
138 // Output slot #0 writes into a TList container
139 DefineOutput(1, TList::Class());
144 //________________________________________________________________________
145 void AliAnalysisDeuteronpA::Initialize()
148 // updating parameters in case of changes
152 fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE,kTRUE);
153 fESDtrackCuts->SetMaxDCAToVertexXY(3);
154 fESDtrackCuts->SetMaxDCAToVertexZ(2);
155 fESDtrackCuts->SetEtaRange(-0.9,0.9);
162 //________________________________________________________________________
163 void AliAnalysisDeuteronpA::UserCreateOutputObjects()
166 //Create analysis utils for event selection and pileup rejection
167 fUtils = new AliAnalysisUtils();
168 fUtils->SetCutOnZVertexSPD(kFALSE);
173 fListHist = new TList();
174 fListHist->SetOwner(kTRUE);
176 const Int_t kPtBins = 28;
177 const Int_t kMultBins = 11;
178 const Int_t kDcaBins = 38;
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;
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};
197 // provide a reasonable dca/binning
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;
206 // create the histograms with all necessary information --> it is filled 4x for each particle assumption
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
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
217 // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident prim, 3-second weak, 4-second material, 5-misident sec
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);
225 fHistRealTracks->GetAxis(2)->Set(kPtBins, binsPt);
226 fHistRealTracks->GetAxis(8)->Set(kDcaBins, binsDca);
227 fListHist->Add(fHistRealTracks);
230 fHistPidQA = new TH3F("fHistPidQA","PID QA",500,0.1,10,1000,0,1000,2,-2,2);
231 BinLogAxis(fHistPidQA);
232 fListHist->Add(fHistPidQA);
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};
241 // different binning for CODE axis, if we want to save motherPDG
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);
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);
253 fHistMomCorr = new TH2F("fHistMomCorr","momentum correction; pT_{MCtrue}; rel. difference",200,0,3,200,-0.2,0.2);
254 fListHist->Add(fHistMomCorr);
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);
259 fHistVertex = new TH2F("fHistVertex","vertex position; V_{XY} (cm); V_{Z}",2000,-1.,1.,400,-20.,20.);
260 fListHist->Add(fHistVertex);
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);
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);
268 PostData(1, fListHist);
272 //________________________________________________________________________
273 void AliAnalysisDeuteronpA::UserExec(Option_t *)
278 if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
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);
283 //AliLog::SetGlobalLogLevel(AliLog::kError);
285 // Check Monte Carlo information and other access first:
287 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
289 //Printf("ERROR: Could not retrieve MC event handler");
293 AliMCEvent* mcEvent = 0x0;
294 AliStack* stack = 0x0;
295 if (eventHandler) mcEvent = eventHandler->MCEvent();
297 //Printf("ERROR: Could not retrieve MC event");
301 stack = mcEvent->Stack();
305 fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
307 //Printf("ERROR: fESD not available");
311 //check with analysis utils
312 //first event in chunck --> continue
313 if(fUtils->IsFirstEventInChunk(fESD)) {
314 PostData(1, fListHist);
318 if (fUtils->IsPileUpEvent(fESD)){
319 PostData(1, fListHist);
323 Bool_t isVtxOk = kTRUE;
324 if(!fUtils->IsVertexSelected2013pA(fESD)) isVtxOk = kFALSE;
325 //NEED TO PUT THIS IN ACTION SOMEWHERE
328 if (!fESDtrackCuts) {
329 Printf("ERROR: fESDtrackCuts not available");
333 // check if event is selected by physics selection class
336 // monitor vertex position
338 const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
339 if(vertex->GetNContributors()<1) {
341 vertex = fESD->GetPrimaryVertexSPD();
342 if(vertex->GetNContributors()<1) vertex = 0x0;
345 // fill histogram to monitor vertex position
346 if (vertex && isVtxOk) fHistVertex->Fill(TMath::Sqrt(vertex->GetX()*vertex->GetX() + vertex->GetY()*vertex->GetY()),vertex->GetZ());
348 Float_t vertexRes = -100;
349 Int_t nPrimaries = 0;
352 if (fMCtrue && vertex && isVtxOk) {
355 for(Int_t i = 0; i < stack->GetNtrack(); i++) {
356 TParticle * trackMC = stack->Particle(i);
358 if (!trackMC->IsPrimary()) continue;
360 //Double_t xv = trackMC->Vx();
361 //Double_t yv = trackMC->Vy();
362 Double_t zv = trackMC->Vz();
364 vertexRes = vertex->GetZ() - zv;
365 fHistVertexRes->Fill(vertex->GetZ() - zv);
370 for(Int_t i = 0; i < stack->GetNtrack(); i++) {
371 TParticle * trackMC = stack->Particle(i);
373 // if (trackMC->IsPrimary()) {
374 if (stack->IsPhysicalPrimary(i)) {
375 if (trackMC->Eta() > -0.9 && trackMC->Eta() < 0.9 && trackMC->GetPDG()->Charge() != 0){
377 //cout << nPrimaries << endl;
383 fHistVertexResTracks->Fill(nPrimaries, vertexRes);
389 Float_t centrality = -1;
392 Int_t rootS = fESD->GetBeamEnergy() < 1000 ? 0 : 1;
393 if (fESD->GetEventSpecie() == 4) { // PbPb
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");
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");
410 Int_t processCode = 0;
412 // important change: fill generated only after vertex cut in case of heavy-ions
414 if ( fESD->GetEventSpecie() == 4) {
415 if (!vertex || !isVtxOk) {
416 fHistMult->Fill(-1, processCode);
417 PostData(1, fListHist);
420 if (TMath::Abs(vertex->GetZ()) > 10) {
421 fHistMult->Fill(-1, processCode);
422 PostData(1, fListHist);
431 for(Int_t i = 0; i < stack->GetNtrack(); i++) {
432 TParticle * trackMC = stack->Particle(i);
433 Int_t pdg = trackMC->GetPdgCode();
435 Double_t xv = trackMC->Vx();
436 Double_t yv = trackMC->Vy();
437 //Double_t zv = trackMC->Vz();
439 dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
441 //dz = TMath::Abs(zv); // so stupid to avoid warnings
443 // vertex cut - selection of primaries
445 //if (dxy > 3 || dz > 10) continue; // fixed cut at 3cm in r
447 if (!stack->IsPhysicalPrimary(i)) continue;
449 // fill MC histograms here...
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();
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;
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);
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.);
475 fHistMult->Fill(-1, processCode);
476 PostData(1, fListHist);
479 if (TMath::Abs(vertex->GetZ()) > 10) {
480 fHistMult->Fill(-1, processCode);
481 PostData(1, fListHist);
486 // count events after physics selection and after vertex selection
488 fHistMult->Fill(fESDtrackCuts->CountAcceptedTracks(fESD), 0);
489 fHistMult->Fill(fESD->GetNumberOfTracks(), 1);
490 fHistCentrality->Fill(centrality);
492 //***************************************************
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*/;
498 Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
500 for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
502 AliESDtrack *track =fESD->GetTrack(i);
503 if (!track->GetInnerParam()) continue;
505 Double_t ptot = track->GetInnerParam()->GetP(); // momentum for dEdx determination
507 // momentum correction for different mass assumption in tracking
509 Double_t pT = track->Pt()/(1 - 0.333303/TMath::Power(track->Pt() + 0.651111, 5.27268));
510 // -[0]/TMath::Power(x - [1],[2])
512 track->GetImpactParameters(dca, cov);
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;
517 // 2.a) apply some standard track cuts according to general recommendations
519 if (!fESDtrackCuts->AcceptTrack(track)) continue;
521 UInt_t status = track->GetStatus();
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;
530 // calculate rapidities and kinematics
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));
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));
547 rapDeuteron = rapDeuteron + 0.465;
548 rapTriton = rapTriton + 0.465;
549 rapHe3 = rapHe3 + 0.465;
555 Double_t sign = track->GetSign();
556 Double_t tpcSignal = track->GetTPCsignal();
559 // 3.a. calculate expected signals in nsigma
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
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
571 // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident, 3-second weak, 4-second material
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]);
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),
593 // Process TOF information
595 //Float_t time0 = fESDpid->GetTOFResponse().GetTimeZero();
596 Float_t time0 = fESDpid->GetTOFResponse().GetStartTime(track->P());//fESDpid->GetTOFResponse().GetTimeZero();
602 time = track->GetTOFsignal() - time0;
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.
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));
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;
623 if (hasTOFout && hasTOFtime && TMath::Abs(pullsTPC[1]) < 3) fHistTofQA->Fill(ptot*sign, beta);
626 for(Int_t iPart = 0; iPart < 3; iPart++) { // loop over assumed particle type
628 // temporary measure to reduce memory: fill only deuterons
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.
636 // using MC truth for precise efficiencies...
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);
642 if (iPart == 0) assumedPdg = 1000010020;
643 if (iPart == 1) assumedPdg = 1000010030;
644 if (iPart == 2) assumedPdg = 1000020030;
647 TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
648 Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
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()))) {
656 if (pdg == assumedPdg && stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()))) code = 4;
658 // fill momentum correction histogram
660 if (code == 1) fHistMomCorr->Fill(trackMC->Pt(), (pT - trackMC->Pt())/trackMC->Pt());
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());
667 // check TOF mismatch on MC basis with TOF label
669 Int_t tofLabel[3] = {0,0,0};
670 track->GetTOFLabel(tofLabel);
673 if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) {
675 if (TMath::Abs(trackMC->GetFirstMother()) == TMath::Abs(tofLabel[0])) hasTOF = kTRUE;
678 // IMPORTANT BIG PROBLEM HERE THE PROBABLILITY TO HAVE A PID SIGNAL MUST BE IN !!!!!!!!!!!!
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);
687 } // end loop over assumed particle type
690 } // end of track loop
693 PostData(1, fListHist);
698 //________________________________________________________________________
699 void AliAnalysisDeuteronpA::Terminate(Option_t *)
701 // Draw result to the screen
702 // Called once at the end of the query
703 Printf("*** CONSTRUCTOR CALLED ****");
708 //________________________________________________________________________
709 void AliAnalysisDeuteronpA::BinLogAxis(const TH1 *h) {
711 // Method for the correct logarithmic binning of histograms
713 TAxis *axis = const_cast<TAxis*>(h->GetXaxis());
714 int bins = axis->GetNbins();
716 Double_t from = axis->GetXmin();
717 Double_t to = axis->GetXmax();
718 Double_t *newBins = new Double_t[bins + 1];
721 Double_t factor = pow(to/from, 1./bins);
723 for (int i = 1; i <= bins; i++) {
724 newBins[i] = factor * newBins[i-1];
726 axis->Set(bins, newBins);