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