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