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