PWGLF/SPECTRA/PiKaPr/TPCTOFpA/AliAnalysisTPCTOFpA.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 in pPb.                //
  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 "AliAnalysisTPCTOFpA.h"
  62
  63
  64 ClassImp(AliAnalysisTPCTOFpA)
  65
  66 //________________________________________________________________________
  67 AliAnalysisTPCTOFpA::AliAnalysisTPCTOFpA()
  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   fSaveMotherPDG(0),
  74   fSmallTHnSparse(0),
  75   fIspA(0),
  76   fRapCMS(0),
  77   fTPCnSigmaCutLow(0),
  78   fTPCnSigmaCutHigh(0),
  79   fRapidityCutLow(0),
  80   fRapidityCutHigh(0),
  81   fEvenDCAbinning(0),
  82   fAlephParameters(),
  83   fHistRealTracks(0),
  84   fHistMCparticles(0),
  85   fHistPidQA(0),
  86   fHistMult(0),
  87   fHistCentrality(0)
  88 {
  89   // default Constructor
  90   /* fast compilation test
  91      gSystem->Load("libANALYSIS");
  92      gSystem->Load("libANALYSISalice");
  93      .L /d/alice09/akalweit/train/trunk/akalweit_hadronSpectra/AliAnalysisTPCTOFpA.cxx++
  94    */
  95 }
  96
  97
  98 //________________________________________________________________________
  99 AliAnalysisTPCTOFpA::AliAnalysisTPCTOFpA(const char *name)
 100   : AliAnalysisTaskSE(name), fESD(0), fListHist(0), fESDtrackCuts(0),fESDTrackCutsMult(0),fESDpid(0),
 101     fMCtrue(0),
 102     fOnlyQA(0),
 103     fUseHBTmultiplicity(0),
 104     fUseTPConlyTracks(0),
 105     fSaveMotherPDG(0),
 106     fSmallTHnSparse(0),
 107     fIspA(0),
 108     fRapCMS(0),
 109     fTPCnSigmaCutLow(0),
 110     fTPCnSigmaCutHigh(0),
 111     fRapidityCutLow(0),
 112     fRapidityCutHigh(0),
 113     fEvenDCAbinning(0),
 114     fAlephParameters(),
 115     fHistRealTracks(0),
 116     fHistMCparticles(0),
 117     fHistPidQA(0),
 118     fHistMult(0),
 119     fHistCentrality(0)
 120     {
 121   //
 122   // standard constructur which should be used
 123   //
 124   Printf("*** CONSTRUCTOR CALLED ****");
 125   //
 126   fMCtrue = kTRUE;
 127   fOnlyQA = kFALSE;
 128   fUseHBTmultiplicity = kTRUE;
 129   fUseTPConlyTracks = kFALSE;
 130
 131   fSmallTHnSparse = kFALSE;
 132   fTPCnSigmaCutLow = -3.;
 133   fTPCnSigmaCutHigh = 3.;
 134   fRapidityCutLow = -0.2;
 135   fRapidityCutHigh = 0.2;
 136   fEvenDCAbinning = kFALSE;
 137
 138   /* real */
 139   fAlephParameters[0] = 0.0283086;
 140   fAlephParameters[1] = 2.63394e+01;
 141   fAlephParameters[2] = 5.04114e-11;
 142   fAlephParameters[3] = 2.12543e+00;
 143   fAlephParameters[4] = 4.88663e+00;
 144   //
 145   // initialize PID object
 146   //
 147   //fESDpid = new AliESDpid();
 148   //
 149   // create track cuts
 150   //
 151   fESDtrackCuts = new AliESDtrackCuts("AliESDtrackCuts","AliESDtrackCuts");
 152   //
 153   //Initialize();
 154   // Output slot #0 writes into a TList container
 155   DefineOutput(1, TList::Class());
 156
 157 }
 158
 159
 160 //________________________________________________________________________
 161 void AliAnalysisTPCTOFpA::Initialize()
 162 {
 163   //
 164   // updating parameters in case of changes
 165   //
 166   // 1. pid parameters
 167   //
 168   //fESDpid->GetTPCResponse().SetBetheBlochParameters(fAlephParameters[0],fAlephParameters[1],fAlephParameters[2],fAlephParameters[3],fAlephParameters[4]);
 169   //
 170   // 2. track cuts
 171   //
 172   /*
 173   fESDtrackCuts->SetMaxCovDiagonalElements(2, 2, 0.5, 0.5, 2);  // BEWARE STANDARD VALUES ARE: 2, 2, 0.5, 0.5, 2
 174   fESDtrackCuts->SetMaxNsigmaToVertex(3);
 175   fESDtrackCuts->SetRequireSigmaToVertex(kTRUE);
 176   fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
 177   fESDtrackCuts->SetMinNClustersTPC(70);
 178   fESDtrackCuts->SetMaxChi2PerClusterTPC(4);
 179   fESDtrackCuts->SetMaxDCAToVertexXY(3);
 180   fESDtrackCuts->SetMaxDCAToVertexZ(3);
 181   fESDtrackCuts->SetRequireTPCRefit(kTRUE);
 182   fESDtrackCuts->SetRequireITSRefit(kTRUE);
 183   fESDtrackCuts->SetClusterRequirementITS(AliESDtrackCuts::kSPD, AliESDtrackCuts::kAny); //TEMPORARY <-> REMOVE
 184   fESDtrackCuts->SetMinNClustersITS(3);
 185   */
 186   //fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE); // kTRUE = sel. primaries --> patch for the moment, do TFractionFitter later
 187
 188
 189   if (!fUseTPConlyTracks) {
 190           fESDtrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 191           fESDtrackCuts->SetMaxDCAToVertexXY(3);
 192           fESDtrackCuts->SetMaxDCAToVertexZ(3);
 193           fESDtrackCuts->SetEtaRange(-0.9,0.9);
 194   }
 195   else {
 196           //fESDtrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 197           fESDtrackCuts->SetMinNClustersTPC(70);
 198           fESDtrackCuts->SetMaxChi2PerClusterTPC(4);
 199           fESDtrackCuts->SetAcceptKinkDaughters(kFALSE);
 200           fESDtrackCuts->SetRequireTPCRefit(kFALSE);
 201
 202           fESDtrackCuts->SetMaxDCAToVertexXY(15);
 203           fESDtrackCuts->SetMaxDCAToVertexZ(6);
 204           fESDtrackCuts->SetDCAToVertex2D(kFALSE);
 205           fESDtrackCuts->SetRequireSigmaToVertex(kFALSE);
 206
 207           fESDtrackCuts->SetEtaRange(-0.9,0.9);
 208   }
 209   //
 210   //
 211   //
 212   //
 213   fESDTrackCutsMult = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 214   fESDTrackCutsMult->SetEtaRange(-1.2,+1.2);
 215   fESDTrackCutsMult->SetPtRange(0.15,1e10);
 216
 217 }
 218
 219
 220 //________________________________________________________________________
 221 void AliAnalysisTPCTOFpA::UserCreateOutputObjects()
 222 {
 223   // Create histograms
 224   // Called once
 225   fListHist = new TList();
 226   fListHist->SetOwner(kTRUE);
 227   //
 228   const Int_t kPtBins = 35;
 229   const Int_t kMultBins = 11;
 230
 231   Int_t kDcaBinsTemp = 76;
 232   if (fEvenDCAbinning) kDcaBinsTemp = 150;
 233   const Int_t kDcaBins = (const Int_t) kDcaBinsTemp;
 234
 235   const Float_t kDcaBinsTPConlyFactor = 5; //need to change binning of DCA plot for tpconly
 236   // sort pT-bins ..
 237   Double_t binsPt[77] = {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};
 238
 239   Double_t binsDca[77] =  {-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};
 240
 241   // DCA bins borders get multiplied by constant factor for TPConlyTracks
 242   Double_t binsDcaTPConly[77];
 243   for (Int_t i = 0; i< 77; i++) {
 244         binsDcaTPConly[i] = kDcaBinsTPConlyFactor * binsDca[i];
 245   }
 246
 247   //
 248   // create the histograms with all necessary information --> it is filled 4x for each particle assumption
 249   //
 250   // (0.) assumed particle: 0. pion, 1. kaon, 2. proton, 3. deuteron
 251   // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
 252   // (2.) pT
 253   // (3.) sign
 254   // (4.) rapidity --> filled 4x
 255   // (5.)  pull TPC dEx --> filled 4x
 256   // (6.) has valid TOF pid signal
 257   // (7.) nsigma TOF --> filled 4x
 258   // (8..) dca_xy
 259   // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident prim, 3-second weak, 4-second material, 5-misident sec, 6-sec. K0, 7-sec. lambda, 8-sec sigma+
 260   //
 261
 262   //dimensions of standard THnSparse
 263   //                              0,           1,           2,  3,       4,    5,     6,   7,     8
 264   Int_t    binsHistReal[9] = {   3,   kMultBins,     kPtBins,  2,       10,    50,    2,  80,    kDcaBins};
 265   Double_t xminHistReal[9] = {-0.5,        -0.5,           0, -2,      -0.5,   -5, -0.5,  -8,       -3};
 266   Double_t xmaxHistReal[9] = { 2.5,        10.5,           3,  2,       0.5,    5,  1.5,   8,        3};
 267
 268   //dimensions of small THnSparse
 269   //                              0,           1,           2,  3,                        4,   5,       6
 270   Int_t    binsHistRealSm[7] = {   3,   kMultBins,     kPtBins,  2,  /*    10,   50,*/    2,  80, kDcaBins};
 271   Double_t xminHistRealSm[7] = {-0.5,        -0.5,           0, -2,  /*  -0.5,   -5,*/ -0.5,  -8,       -3};
 272   Double_t xmaxHistRealSm[7] = { 2.5,        10.5,           3,  2,  /*   0.5,    5,*/  1.5,   8,        3};
 273
 274   if (!fSmallTHnSparse) fHistRealTracks = new THnSparseF("fHistRealTracks","real tracks",9,binsHistReal,xminHistReal,xmaxHistReal);
 275   else                  fHistRealTracks = new THnSparseF("fHistRealTracks","real tracks",7,binsHistRealSm,xminHistRealSm,xmaxHistRealSm);
 276   //
 277   fHistRealTracks->GetAxis(2)->Set(kPtBins, binsPt);
 278
 279   //different DCAxy binning for TPConlyTracks
 280
 281   Int_t dcaAxisNumber = 8;
 282   if (fSmallTHnSparse) dcaAxisNumber = 6;
 283
 284   if (!fUseTPConlyTracks){
 285     if (fEvenDCAbinning) fHistRealTracks->GetAxis(dcaAxisNumber)->Set(kDcaBins,-3.,3.);
 286     else fHistRealTracks->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDca);
 287   }
 288   else {
 289     if (fEvenDCAbinning) fHistRealTracks->GetAxis(dcaAxisNumber)->Set(kDcaBins,-15.,15.);
 290     else fHistRealTracks->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDcaTPConly);
 291    }
 292   fListHist->Add(fHistRealTracks);
 293   //
 294   //                      0.ptot,1.tpcSig,2.hasTOF, 3. assumed part., 4. nclDedx, 5. nSigmaTPC (4x), 6. nSigmaTOF (4x), 7. centrality
 295   fHistPidQA = new TH3D("fHistPidQA","PID QA",500,0.1,10,1000,0,1000,2,-2,2);
 296   BinLogAxis(fHistPidQA);
 297   fListHist->Add(fHistPidQA);
 298
 299   // dimensions of standard THnSparse
 300   //                            0,            1,           2,  3,      4,   5,    6,   7,       8.,    9.
 301   Int_t    binsHistMC[10] = {   3,    kMultBins,     kPtBins,  2,     10,  50,    2,  80, kDcaBins,    6};
 302   Double_t xminHistMC[10] = {-0.5,         -0.5,           0, -2,   -0.5,  -5, -0.5,  -8,       -3, -0.5};
 303   Double_t xmaxHistMC[10] = { 2.5,         10.5,           3,  2,    0.5,   5,  1.5,   8,        3,  5.5};
 304
 305   //dimensions of small THnSparse
 306   //                              0,           1,           2,  3,                     4,   5,       6.,    7.
 307   Int_t    binsHistMCSm[8] = {   3,    kMultBins,     kPtBins,  2,  /*   10,  50,*/    2,  80, kDcaBins,    6};
 308   Double_t xminHistMCSm[8] = {-0.5,         -0.5,           0, -2,  /* -0.5,  -5,*/ -0.5,  -8,       -3, -0.5};
 309   Double_t xmaxHistMCSm[8] = { 2.5,         10.5,           3,  2,  /*  0.5,   5,*/  1.5,   8,        3,  5.5};
 310
 311   //different binning for CODE axis, if we want to save motherPDG
 312   if (fSaveMotherPDG) {
 313     binsHistMC[9] = 9;
 314     xmaxHistMC[9] = 8.5;
 315     binsHistMCSm[7] = 9;
 316     xmaxHistMCSm[7] = 8.5;
 317   }
 318
 319
 320   if (!fSmallTHnSparse) fHistMCparticles = new THnSparseF("fHistMCparticles","MC histogram",10,binsHistMC,xminHistMC,xmaxHistMC);
 321   else                  fHistMCparticles = new THnSparseF("fHistMCparticles","MC histogram",8,binsHistMCSm,xminHistMCSm,xmaxHistMCSm);
 322
 323
 324   fHistMCparticles->GetAxis(2)->Set(kPtBins, binsPt);
 325
 326   /*
 327   //different DCAxy binning for TPConlyTracks
 328   if (!fEvenDCAbinning){
 329     if (!fUseTPConlyTracks) fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDca);
 330     else fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDcaTPConly);
 331   }
 332   */
 333
 334   if (!fUseTPConlyTracks){
 335     if (fEvenDCAbinning) fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins,-3.,3.);
 336     else fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDca);
 337   }
 338   else {
 339     if (fEvenDCAbinning) fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins,-15.,15.);
 340     else fHistMCparticles->GetAxis(dcaAxisNumber)->Set(kDcaBins, binsDcaTPConly);
 341    }
 342
 343
 344
 345   fListHist->Add(fHistMCparticles);
 346   //
 347   fHistMult = new TH2D("fHistMult", "control histogram to count number of events", 502, -2.5, 499.5,4,-0.5,3.5);
 348   fHistCentrality = new TH1D("fHistCentrality", "control histogram to count number of events", 22, -1.5, 20.5);
 349   fListHist->Add(fHistMult);
 350   fListHist->Add(fHistCentrality);
 351
 352 }
 353
 354 //________________________________________________________________________
 355 void AliAnalysisTPCTOFpA::UserExec(Option_t *)
 356 {
 357   //
 358   // main event loop
 359   //
 360   if (!fESDpid) fESDpid = ((AliESDInputHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->GetESDpid();
 361   if (!fESDpid) {
 362     fESDpid = new AliESDpid(); // HACK FOR MC PBPB --> PLEASE REMOVE AS SOON AS POSSIBLE
 363     fESDpid->GetTPCResponse().SetBetheBlochParameters(1.28778e+00/50., 3.13539e+01, TMath::Exp(-3.16327e+01), 1.87901e+00, 6.41583e+00);
 364   }
 365   //AliLog::SetGlobalLogLevel(AliLog::kError);
 366   //
 367   // Check Monte Carlo information and other access first:
 368   //
 369   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 370   if (!eventHandler) {
 371     //Printf("ERROR: Could not retrieve MC event handler");
 372     fMCtrue = kFALSE;
 373   }
 374   //
 375   AliMCEvent* mcEvent = 0x0;
 376   AliStack* stack = 0x0;
 377   if (eventHandler) mcEvent = eventHandler->MCEvent();
 378   if (!mcEvent) {
 379     //Printf("ERROR: Could not retrieve MC event");
 380     if (fMCtrue) return;
 381   }
 382   if (fMCtrue) {
 383     stack = mcEvent->Stack();
 384     if (!stack) return;
 385   }
 386   //
 387   fESD = dynamic_cast<AliESDEvent*>( InputEvent() );
 388   if (!fESD) {
 389     //Printf("ERROR: fESD not available");
 390     return;
 391   }
 392
 393   if (!fESDtrackCuts) {
 394     Printf("ERROR: fESDtrackCuts not available");
 395     return;
 396   }
 397   //
 398   // check if event is selected by physics selection class
 399   //
 400   Bool_t isSelected = kTRUE; // for reasons of backward compatibility --> check is now in AddTask macro
 401   /*
 402   Bool_t isSelected = kFALSE;
 403   isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB) == AliVEvent::kMB;
 404   */
 405   //
 406   // monitor vertex position
 407   //
 408   const AliESDVertex *vertex = fESD->GetPrimaryVertexTracks();
 409   if(vertex->GetNContributors()<1) {
 410     // SPD vertex
 411     vertex = fESD->GetPrimaryVertexSPD();
 412     if(vertex->GetNContributors()<1) vertex = 0x0;
 413   }
 414   //
 415   // small track loop to determine trigger Pt, multiplicity or centrality
 416   //
 417   Double_t triggerPt = 0;
 418   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
 419     AliESDtrack *track =fESD->GetTrack(i);
 420     if (!fESDTrackCutsMult->AcceptTrack(track)) continue;
 421     if (track->Pt() > triggerPt) triggerPt = track->Pt();
 422   }
 423   Int_t trackCounter = fESDTrackCutsMult->CountAcceptedTracks(fESD);
 424   //
 425   // 2nd multiplicity estimator SPD hits; replaces multiplicity from HBT
 426   //
 427   Float_t spdCorr = -1;
 428   const AliMultiplicity *mult = fESD->GetMultiplicity();
 429   Float_t nClusters[6]={0.0,0.0,0.0,0.0,0.0,0.0};
 430   for(Int_t ilay=0; ilay<6; ilay++) nClusters[ilay] = (Float_t)mult->GetNumberOfITSClusters(ilay);
 431   if (vertex) spdCorr = AliESDUtils::GetCorrSPD2(nClusters[1],vertex->GetZ());
 432   //
 433   Float_t centrality = -1;
 434   //
 435   // IMPORTANT CENTRALITY DEFINITION FOR pp
 436   //
 437   if (!fUseHBTmultiplicity) {
 438     // test binning to check if multiplicity-dependence is due to diffractive events being s
 439     if (trackCounter >= 0  && trackCounter <= 0)  centrality = 0;
 440     if (trackCounter >= 1  && trackCounter <= 1)  centrality = 1;
 441     if (trackCounter >= 2  && trackCounter <= 2)  centrality = 2;
 442     if (trackCounter >= 3  && trackCounter <= 3)  centrality = 3;
 443     if (trackCounter >= 4  && trackCounter <= 4)  centrality = 4;
 444     if (trackCounter >= 5  && trackCounter <= 5)  centrality = 5;
 445     // this was all the original bin 1 being [0..5]
 446     if (trackCounter >= 6  && trackCounter <= 9)  centrality = 6;
 447     if (trackCounter >= 10 && trackCounter <= 14) centrality = 7;
 448     if (trackCounter >= 15 && trackCounter <= 22) centrality = 8;
 449     if (trackCounter >= 23 && trackCounter <= 32) centrality = 9;
 450     if (trackCounter >= 33 && trackCounter <= 42) centrality = 10;
 451     /*
 452     if (trackCounter >= 43 && trackCounter <= 52) centrality = 7
 453     if (trackCounter >= 53 && trackCounter <= 62) centrality = 8;
 454     if (trackCounter >= 63 && trackCounter <= 72) centrality = 9;
 455     if (trackCounter >= 73 && trackCounter <= 82) centrality = 10;
 456     */
 457   } else {
 458     if (spdCorr >= 0  && spdCorr <=  2)  centrality  = 0;
 459     if (spdCorr >= 3  && spdCorr <=  5)  centrality  = 1;
 460     if (spdCorr >= 6  && spdCorr <=  8)  centrality  = 2;
 461     if (spdCorr >= 9  && spdCorr <= 11)  centrality  = 3;
 462     if (spdCorr >= 12 && spdCorr <= 14)  centrality  = 4;
 463     if (spdCorr >= 15 && spdCorr <= 16)  centrality  = 5;
 464     // this was all the original bin 1 being [0..16]
 465     if (spdCorr >= 17 && spdCorr <= 30)  centrality =  6;
 466     if (spdCorr >= 31 && spdCorr <= 45)  centrality =  7;
 467     if (spdCorr >= 46 && spdCorr <= 68)  centrality =  8;
 468     if (spdCorr >= 69 && spdCorr <= 97)  centrality =  9;
 469     if (spdCorr >= 98)                   centrality = 10;
 470     /*
 471     if (spdCorr >= 17 && spdCorr <= 30)  centrality = 2;
 472     if (spdCorr >= 31 && spdCorr <= 45)  centrality = 3;
 473     if (spdCorr >= 46 && spdCorr <= 68)  centrality = 4;
 474     if (spdCorr >= 69 && spdCorr <= 97)  centrality = 5;
 475     if (spdCorr >= 98)                   centrality = 6;
 476     */
 477   }
 478   //
 479   Int_t rootS = fESD->GetBeamEnergy() < 1000 ? 0 : 1;
 480   if (fESD->GetEventSpecie() == 4) { // PbPb
 481     rootS = 2;
 482     AliCentrality *esdCentrality = fESD->GetCentrality();
 483     centrality = esdCentrality->GetCentralityClass10("V0M") + 1; // centrality percentile determined with V0
 484     if (TMath::Abs(centrality - 1) < 1e-5) {
 485       centrality = esdCentrality->GetCentralityClass5("V0M");
 486     }
 487   }
 488
 489   if (fIspA) {
 490     AliCentrality *esdCentrality = fESD->GetCentrality();
 491     Float_t pApercentile = esdCentrality->GetCentralityPercentile("V0M"); // centrality percentile determined with V0A
 492     if (pApercentile >=  0. && pApercentile < 20.) centrality = 0;
 493     if (pApercentile >= 20. && pApercentile < 40.) centrality = 1;
 494     if (pApercentile >= 40. && pApercentile < 60.) centrality = 2;
 495     if (pApercentile >= 60. && pApercentile <= 100.) centrality = 3;
 496     //if (pApercentile >= 80. && pApercentile <= 100.) centrality = 4;
 497     //cout << "*****************ispA switch works***************************" << endl;
 498     //cout << "centrality percentile is:  " << pApercentile << endl;
 499     //cout << "*************************************************************" << endl;
 500   }
 501
 502
 503
 504
 505   Int_t nContributors = 0;
 506   if (fESD->GetPrimaryVertexTPC()) nContributors = fESD->GetPrimaryVertexTPC()->GetNContributors();
 507   //
 508
 509   Int_t processtype = 0;
 510   Int_t processCode = 0;
 511   //
 512   // important change: fill generated only after vertex cut in case of heavy-ions
 513   //
 514
 515
 516   if (!vertex) {
 517     fHistMult->Fill(-1, processCode);
 518     PostData(1, fListHist);
 519     return;
 520   } else {
 521     if (TMath::Abs(vertex->GetZv()) > 10) {
 522       fHistMult->Fill(-1, processCode);
 523       PostData(1, fListHist);
 524       return;
 525     }
 526   }
 527   //
 528
 529
 530
 531   if (fMCtrue) {
 532     //
 533     //
 534     //
 535     /*
 536     AliHeader * header = mcEvent->Header();
 537     processtype = GetPythiaEventProcessType(header);
 538     // non diffractive
 539     if (processtype !=92 && processtype !=93 && processtype != 94) processCode = 1;
 540     // single diffractive
 541     if ((processtype == 92 || processtype == 93)) processCode = 2;
 542     // double diffractive
 543     if (processtype == 94) processCode = 3;
 544     //
 545     */
 546
 547     for(Int_t i = 0; i < stack->GetNtrack(); i++) {
 548       TParticle * trackMC = stack->Particle(i);
 549       Int_t pdg = trackMC->GetPdgCode();
 550       //
 551       Double_t xv = trackMC->Vx();
 552       Double_t yv = trackMC->Vy();
 553       Double_t zv = trackMC->Vz();
 554       Double_t dxy = 0;
 555       dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
 556       Double_t dz = 0;
 557       dz = TMath::Abs(zv); // so stupid to avoid warnings
 558       //
 559       // vertex cut - selection of primaries
 560       //
 561       //if (dxy > 3 || dz > 10) continue; // fixed cut at 3cm in r
 562       //
 563       if (!stack->IsPhysicalPrimary(i)) continue;
 564       //
 565       // fill MC histograms here...
 566       //
 567       Double_t rap = trackMC->Y();
 568       if (fRapCMS) rap = rap + 0.465;
 569       Double_t pT  = trackMC->Pt();
 570       Int_t sign = pdg < 0 ? -1 : 1; // only works for charged pi,K,p !!
 571 //      Double_t transMass = TMath::Sqrt(trackMC->Pt()*trackMC->Pt() + trackMC->GetMass()*trackMC->GetMass()) - trackMC->GetMass();
 572       //
 573       Int_t iPart = -1;
 574       if (TMath::Abs(pdg) == 211)  iPart = 0; // select Pi+/Pi- only
 575       if (TMath::Abs(pdg) == 321)  iPart = 1; // select K+/K- only
 576       if (TMath::Abs(pdg) == 2212) iPart = 2; // select p+/p- only
 577       if (iPart == -1) continue;
 578       //
 579
 580       if (!fSmallTHnSparse){
 581         Double_t vecHistMC[10] = {iPart, centrality,  pT, sign, rap, 0, 1, 0, dxy, 0};
 582         if (!fOnlyQA) fHistMCparticles->Fill(vecHistMC);
 583       }
 584       else{
 585         if (rap>fRapidityCutLow && rap<fRapidityCutHigh){
 586           Double_t vecHistMC[8] = {iPart, centrality,  pT, sign, 1, 0, dxy, 0};
 587           if (!fOnlyQA) fHistMCparticles->Fill(vecHistMC);
 588         }
 589       }
 590     }
 591   }
 592   //
 593   if (!isSelected && !fOnlyQA) {
 594     PostData(1, fListHist);
 595     return;
 596   }
 597   //
 598   if (!vertex) {
 599     fHistMult->Fill(-1, processCode);
 600     PostData(1, fListHist);
 601     return;
 602   } else {
 603     if (TMath::Abs(vertex->GetZv()) > 10) {
 604       fHistMult->Fill(-1, processCode);
 605       PostData(1, fListHist);
 606       return;
 607     }
 608   }
 609   //
 610   // count events after physics selection and after vertex selection
 611   //
 612   //cout << "MULTIPLICITY " << trackCounter << " " << fESD->GetEventNumberInFile() <<endl;
 613   fHistMult->Fill(trackCounter, processCode);
 614   fHistCentrality->Fill(centrality);
 615
 616
 617   //***************************************************
 618   // track loop
 619   //***************************************************
 620   //const Float_t kNsigmaCut = 3;
 621   //const Float_t k2sigmaCorr = 1/(0.5*(TMath::Erf(kNsigmaCut/sqrt(2))-TMath::Erf(-kNsigmaCut/sqrt(2))))/*1/0.9545*/;
 622   //
 623   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
 624   //
 625   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
 626
 627         AliESDtrack *track = 0;
 628         AliESDtrack *trackForTOF = 0; //normal track for all TOF information needed when using tpconly-tracks
 629
 630         //normal tracks, if tpconly flag is set, use tpconlytracks
 631         if (!fUseTPConlyTracks){
 632         track =fESD->GetTrack(i);
 633         }
 634         else {
 635                 track = fESDtrackCuts->GetTPCOnlyTrack(fESD,i);
 636                 if (!track) continue;
 637                 trackForTOF = fESD->GetTrack(i);
 638         }
 639     //
 640     if (!track->GetInnerParam()) {
 641                 if (fUseTPConlyTracks) {delete track; track = 0;} //need to delete tpconlytrack
 642                 continue;
 643         }
 644     Double_t ptot = track->GetInnerParam()->GetP(); // momentum for dEdx determination
 645     Double_t pT = track->Pt();
 646     track->GetImpactParameters(dca, cov);
 647     //
 648     //
 649     // cut for dead regions in the detector
 650     // if (track->Eta() > 0.1 && (track->Eta() < 0.2 && track->Phi() > 0.1 && track->Phi() < 0.1) continue;
 651     //
 652     // 2.a) apply some standard track cuts according to general recommendations
 653     //
 654     if (!fESDtrackCuts->AcceptTrack(track)) {
 655                 if (fUseTPConlyTracks) {delete track; track = 0;} //need to delete tpconlytrack
 656                 continue;
 657         }
 658
 659         UInt_t status = 0;
 660     if (!fUseTPConlyTracks) status = track->GetStatus();
 661         else status = trackForTOF->GetStatus();
 662     Bool_t hasTOFout  = status&AliESDtrack::kTOFout;
 663     Bool_t hasTOFtime = status&AliESDtrack::kTIME;
 664     Bool_t hasTOFpid  = status&AliESDtrack::kTOFpid;
 665     Bool_t hasTOF     = kFALSE;
 666     if (hasTOFout && hasTOFtime && hasTOFpid) hasTOF = kTRUE;
 667         Float_t length = 0.;
 668     if (!fUseTPConlyTracks) length = track->GetIntegratedLength();
 669     else length = trackForTOF->GetIntegratedLength();
 670
 671     if (length < 350.) hasTOF = kFALSE;
 672     //
 673     // calculate rapidities and kinematics
 674     //
 675     //
 676     Double_t pvec[3];
 677     track->GetPxPyPz(pvec);
 678     Double_t energyPion = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kPion)*AliPID::ParticleMass(AliPID::kPion));
 679     Double_t energyKaon = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kKaon)*AliPID::ParticleMass(AliPID::kKaon));
 680     Double_t energyProton = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton));
 681     Double_t energyDeuteron = TMath::Sqrt(track->GetP()*track->GetP() + 4*AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton));
 682     //
 683     Double_t rapPion = 0.5*TMath::Log((energyPion + pvec[2])/(energyPion - pvec[2]));
 684     Double_t rapKaon = 0.5*TMath::Log((energyKaon + pvec[2])/(energyKaon - pvec[2]));
 685     Double_t rapProton = 0.5*TMath::Log((energyProton + pvec[2])/(energyProton - pvec[2]));
 686     Double_t rapDeuteron = 0.5*TMath::Log((energyDeuteron + pvec[2])/(energyDeuteron - pvec[2]));
 687
 688     if (fRapCMS) {
 689       rapPion = rapPion + 0.465;
 690       rapKaon = rapKaon + 0.465;
 691       rapProton = rapProton + 0.465;
 692       rapDeuteron = rapDeuteron + 0.465;
 693     }
 694
 695
 696     //
 697 //    Double_t transMassPion = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kPion)*AliPID::ParticleMass(AliPID::kPion))      -  AliPID::ParticleMass(AliPID::kPion);
 698 //    Double_t transMassKaon = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kKaon)*AliPID::ParticleMass(AliPID::kKaon))     -  AliPID::ParticleMass(AliPID::kKaon);
 699  //   Double_t transMassProton = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))    -  AliPID::ParticleMass(AliPID::kProton);
 700 //    Double_t transMassDeuteron = TMath::Sqrt(track->Pt()*track->Pt() + 4*AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))    -  2*AliPID::ParticleMass(AliPID::kProton);
 701     //
 702     // 3. make the PID
 703     //
 704     Double_t sign = track->GetSign();
 705     Double_t tpcSignal = track->GetTPCsignal();
 706     //
 707     //
 708     // 3.a. calculate expected signals in nsigma
 709     //
 710     //
 711     // (0.) assumed particle: 0. pion, 1. kaon, 2. proton, 3. deuteron
 712     // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
 713     // (2.) pT
 714     // (3.) sign
 715     // (4.) rapidity --> filled 4x
 716     // (5.)  pull TPC dEx --> filled 4x
 717     // (6.) has valid TOF pid signal
 718     // (7.) nsigma TOF --> filled 4x
 719     // (8..) dca_xy
 720     // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident, 3-second weak, 4-second material
 721     //
 722 //    Double_t transMass[4] = {transMassPion,transMassKaon,transMassProton,transMassDeuteron};
 723     Double_t rap[4] = {rapPion,rapKaon,rapProton,rapDeuteron};
 724     Double_t pullsTPC[4] = {fESDpid->NumberOfSigmasTPC(track,AliPID::kPion),
 725                             fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon),
 726                             fESDpid->NumberOfSigmasTPC(track,AliPID::kProton),
 727                             0}; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!
 728     Float_t time0 = fESDpid->GetTOFResponse().GetTimeZero();
 729     //fESDpid->GetTOFResponse().SetTimeResolution(130.);
 730     Double_t pullsTOF[4] ={0.,0.,0.,0.};
 731     if (!fUseTPConlyTracks) {
 732       pullsTOF[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kPion, time0);
 733       pullsTOF[1] = fESDpid->NumberOfSigmasTOF(track,AliPID::kKaon, time0);
 734       pullsTOF[2] = fESDpid->NumberOfSigmasTOF(track,AliPID::kProton, time0);
 735       pullsTOF[3] = 0; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!;
 736     }
 737     else {
 738       pullsTOF[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kPion, time0);
 739       pullsTOF[1] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kKaon, time0);
 740       pullsTOF[2] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kProton, time0);
 741       pullsTOF[3] = 0; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!;
 742     }
 743
 744     //
 745 //    Double_t tpcQA[4] = {fESDpid->NumberOfSigmasTPC(track,AliPID::kElectron),
 746 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kPion),
 747 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon),
 748 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kProton)};
 749
 750     Double_t tofQA[4] = {0.,0.,0.,0.};
 751     if (!fUseTPConlyTracks) {
 752       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kElectron, time0);
 753       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kPion, time0);
 754       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kKaon, time0);
 755       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kProton, time0);
 756     }
 757     else{
 758       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kElectron, time0);
 759       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kPion, time0);
 760       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kKaon, time0);
 761       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kProton, time0);
 762     }
 763
 764
 765     //save information for every particle type  // loop over assumed particle type
 766     for(Int_t iPart = 0; iPart < 3; iPart++) {
 767
 768       if (!fSmallTHnSparse) {
 769         Double_t vecHistReal[9]  = {iPart,  centrality,   pT, sign,  rap[iPart], pullsTPC[iPart], hasTOF, pullsTOF[iPart], dca[0]};
 770         if (!fOnlyQA) fHistRealTracks->Fill(vecHistReal);
 771       }
 772       else {
 773         if (pullsTPC[iPart]>fTPCnSigmaCutLow && pullsTPC[iPart]<fTPCnSigmaCutHigh && rap[iPart]>fRapidityCutLow && rap[iPart]<fRapidityCutHigh) {
 774           Double_t vecHistReal[7]  = {iPart,  centrality,   pT, sign, hasTOF, pullsTOF[iPart], dca[0]};
 775           if (!fOnlyQA) fHistRealTracks->Fill(vecHistReal);
 776         }
 777       }
 778
 779
 780       // using MC truth for precise efficiencies...
 781       //
 782       if (fMCtrue && !fOnlyQA) {
 783         Int_t code = 9; // code: 0-generated, 1-true rec. primaries, 2-misident, 3-second weak, 4-second material
 784         Int_t assumedPdg = 0;//2212(proton); 321(Kaon); 211(pion);
 785         Int_t motherCode = -1;
 786         if (iPart == 0) assumedPdg = 211;
 787         if (iPart == 1) assumedPdg = 321;
 788         if (iPart == 2) assumedPdg = 2212;
 789         //
 790         //
 791         TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 792         Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 793         //
 794         if (pdg != assumedPdg && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 2;
 795         if (pdg != assumedPdg && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) code = 5;
 796         if (pdg == assumedPdg && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 1;
 797         if (pdg == assumedPdg && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) {
 798           code = 3;
 799           if (fSaveMotherPDG){
 800             TParticle *trackMother =  stack->Particle(TMath::Abs(trackMC->GetFirstMother()));
 801             if (trackMother->GetPdgCode() == 310) motherCode = 6; //K0
 802             if (trackMother->GetPdgCode() == 3122) motherCode = 7; //Lambda
 803             if (trackMother->GetPdgCode() == 3222) motherCode = 8; //Sigma+
 804           }
 805         }
 806         if (pdg == assumedPdg && stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()))) code = 4;
 807
 808         //
 809         // muons need special treatment, because they are indistinguishable from pions
 810         //
 811         if (iPart == 0 && pdg == 13  && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 1;
 812         if (iPart == 0 && pdg == 13  && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) code = 3;
 813         //
 814         // check TOF mismatch on MC basis with TOF label
 815         //
 816         Int_t tofLabel[3];
 817         if (!fUseTPConlyTracks) track->GetTOFLabel(tofLabel);
 818         else trackForTOF->GetTOFLabel(tofLabel);
 819         if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) hasTOF = kFALSE;
 820         //
 821         // IMPORTANT BIG PROBLEM HERE THE PROBABLILITY TO HAVE A PID SIGNAL MUST BE IN !!!!!!!!!!!!
 822         //
 823         if (!fSmallTHnSparse){
 824           Double_t vectorHistMC[10] = {iPart,  centrality,  pT, sign,  rap[iPart], pullsTPC[iPart], hasTOF, pullsTOF[iPart], dca[0], code};
 825           if (!fOnlyQA) {
 826             fHistMCparticles->Fill(vectorHistMC);
 827             if (motherCode != -1 && fSaveMotherPDG) { //if mother of weak decay is K0, lambda or sigma+ add track again with this information
 828               Double_t vectorHistMCmother[10] = {iPart,  centrality,  pT, sign,  rap[iPart], pullsTPC[iPart], hasTOF, pullsTOF[iPart], dca[0], motherCode};
 829               fHistMCparticles->Fill(vectorHistMCmother);
 830             }
 831           }
 832         }
 833         else{
 834           if (pullsTPC[iPart]>fTPCnSigmaCutLow && pullsTPC[iPart]<fTPCnSigmaCutHigh && rap[iPart]>fRapidityCutLow && rap[iPart]<fRapidityCutHigh) {
 835             //                              0,           1,   2,    3,           4,               5,      6,               7,      8,   9
 836             Double_t vectorHistMC[8] = {iPart,  centrality,  pT, sign, hasTOF, pullsTOF[iPart], dca[0], code};
 837             if (!fOnlyQA) {
 838               fHistMCparticles->Fill(vectorHistMC);
 839               if (motherCode != -1 && fSaveMotherPDG) { //if mother of weak decay is K0, lambda or sigma+ add track again with this information
 840                 Double_t vectorHistMCmother[8] = {iPart,  centrality,  pT, sign, hasTOF, pullsTOF[iPart], dca[0], motherCode};
 841                 fHistMCparticles->Fill(vectorHistMCmother);
 842               }
 843             }
 844           }
 845         }
 846       }
 847       //
 848       //
 849       Int_t tpcShared = track->GetTPCnclsS();
 850       if (TMath::Abs(track->Eta()) < 0.8 && iPart == 0 && tpcShared < 4) fHistPidQA->Fill(ptot,tpcSignal,sign);
 851     } // end loop over assumed particle type
 852
 853           //need to delete tpconlytrack
 854           if (fUseTPConlyTracks){
 855                 delete track;
 856                 track = 0;
 857           }
 858
 859   } // end of track loop
 860
 861   // Post output data
 862   PostData(1, fListHist);
 863
 864 }
 865
 866
 867 //________________________________________________________________________
 868 void AliAnalysisTPCTOFpA::Terminate(Option_t *)
 869 {
 870   // Draw result to the screen
 871   // Called once at the end of the query
 872   Printf("*** CONSTRUCTOR CALLED ****");
 873
 874 }
 875
 876
 877 //________________________________________________________________________
 878 Bool_t AliAnalysisTPCTOFpA::SelectOnImpPar(AliESDtrack* t) {
 879   //
 880   // cut on transverse impact parameter // DEPRECATED
 881   //
 882   Float_t d0z0[2],covd0z0[3];
 883   t->GetImpactParameters(d0z0,covd0z0);
 884   Float_t sigma= 0.0050+0.0060/TMath::Power(t->Pt(),0.9);
 885   Float_t d0max = 7.*sigma;
 886   //
 887   Float_t sigmaZ = 0.0146+0.0070/TMath::Power(t->Pt(),1.114758);
 888   if (t->Pt() > 1) sigmaZ = 0.0216;
 889   Float_t d0maxZ = 5.*sigmaZ;
 890   //
 891   if(TMath::Abs(d0z0[0]) < d0max && TMath::Abs(d0z0[1]) < d0maxZ) return kTRUE;
 892   return kFALSE;
 893 }
 894
 895
 896 //________________________________________________________________________
 897 void AliAnalysisTPCTOFpA::BinLogAxis(const TH1 *h) {
 898   //
 899   // Method for the correct logarithmic binning of histograms
 900   //
 901   TAxis *axis = h->GetXaxis();
 902   int bins = axis->GetNbins();
 903
 904   Double_t from = axis->GetXmin();
 905   Double_t to = axis->GetXmax();
 906   Double_t *newBins = new Double_t[bins + 1];
 907
 908   newBins[0] = from;
 909   Double_t factor = pow(to/from, 1./bins);
 910
 911   for (int i = 1; i <= bins; i++) {
 912    newBins[i] = factor * newBins[i-1];
 913   }
 914   axis->Set(bins, newBins);
 915   delete [] newBins;
 916
 917 }
 918
 919
 920 //________________________________________________________________________
 921 Int_t AliAnalysisTPCTOFpA::GetPythiaEventProcessType(const AliHeader* aHeader, const Bool_t adebug) const {
 922   //
 923   // get the process type of the event.
 924   //
 925
 926   // can only read pythia headers, either directly or from cocktalil header
 927   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
 928
 929   if (!pythiaGenHeader) {
 930
 931     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
 932     if (!genCocktailHeader) {
 933       //printf("AliAnalysisTPCTOFpA::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
 934       return -1;
 935     }
 936
 937     TList* headerList = genCocktailHeader->GetHeaders();
 938     if (!headerList) {
 939       return -1;
 940     }
 941
 942     for (Int_t i=0; i<headerList->GetEntries(); i++) {
 943       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
 944       if (pythiaGenHeader)
 945         break;
 946     }
 947
 948     if (!pythiaGenHeader) {
 949       //printf("AliAnalysisTPCTOFpA::GetProcessType : Could not find Pythia header. \n");
 950       return -1;
 951     }
 952   }
 953
 954   if (adebug) {
 955     //printf("AliAnalysisTPCTOFpA::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType());
 956   }
 957
 958   return pythiaGenHeader->ProcessType();
 959 }