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