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