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