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 = -99;
 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     else if (pApercentile >=  5. && pApercentile < 10.) centrality = 0;
 542     else if (pApercentile >= 10. && pApercentile < 20.) centrality = 1;
 543     else if (pApercentile >= 20. && pApercentile < 30.) centrality = 2;
 544     else if (pApercentile >= 30. && pApercentile < 40.) centrality = 3;
 545     else if (pApercentile >= 40. && pApercentile < 50.) centrality = 4;
 546     else if (pApercentile >= 50. && pApercentile < 60.) centrality = 5;
 547     else if (pApercentile >= 60. && pApercentile < 70.) centrality = 6;
 548     else if (pApercentile >= 70. && pApercentile < 80.) centrality = 7;
 549     else if (pApercentile >= 80. && pApercentile < 90.) centrality = 8;
 550     else if (pApercentile >= 90. && pApercentile <= 100.) centrality = 9;
 551     else centrality = -99;
 552
 553     /*
 554     cout << "*****************ispA switch works***************************" << endl;
 555     cout << "centrality estimator  is:  " << fCentEst.Data() << endl;
 556     cout << "centrality percentile is:  " << pApercentile << endl;
 557     cout << "*************************************************************" << endl;
 558     */
 559   }
 560
 561
 562
 563
 564   Int_t nContributors = 0;
 565   if (fESD->GetPrimaryVertexTPC()) nContributors = fESD->GetPrimaryVertexTPC()->GetNContributors();
 566   //
 567
 568   //  Int_t processtype = 0;
 569   Int_t processCode = 0;
 570   //
 571   // important change: fill generated only after vertex cut in case of heavy-ions
 572   //
 573
 574
 575   if (!vertex || !isVertexOk) {
 576     fHistMult->Fill(-1, processCode);
 577     PostData(1, fListHist);
 578     return;
 579   } else {
 580     if (TMath::Abs(vertex->GetZv()) > 10) {
 581       fHistMult->Fill(-1, processCode);
 582       PostData(1, fListHist);
 583       return;
 584     }
 585   }
 586   //
 587
 588
 589
 590   if (fMCtrue) {
 591     //
 592     //
 593     //
 594     /*
 595     AliHeader * header = mcEvent->Header();
 596     processtype = GetPythiaEventProcessType(header);
 597     // non diffractive
 598     if (processtype !=92 && processtype !=93 && processtype != 94) processCode = 1;
 599     // single diffractive
 600     if ((processtype == 92 || processtype == 93)) processCode = 2;
 601     // double diffractive
 602     if (processtype == 94) processCode = 3;
 603     //
 604     */
 605
 606     for(Int_t i = 0; i < stack->GetNtrack(); i++) {
 607       TParticle * trackMC = stack->Particle(i);
 608       Int_t pdg = trackMC->GetPdgCode();
 609       //
 610       Double_t xv = trackMC->Vx();
 611       Double_t yv = trackMC->Vy();
 612       Double_t zv = trackMC->Vz();
 613       Double_t dxy = 0;
 614       dxy = TMath::Sqrt(xv*xv + yv*yv); // so stupid to avoid warnings
 615       Double_t dz = 0;
 616       dz = TMath::Abs(zv); // so stupid to avoid warnings
 617       //
 618       // vertex cut - selection of primaries
 619       //
 620       //if (dxy > 3 || dz > 10) continue; // fixed cut at 3cm in r
 621       //
 622       if (!stack->IsPhysicalPrimary(i)) continue;
 623       //
 624       // fill MC histograms here...
 625       //
 626       Double_t rap = trackMC->Y();
 627       if (fRapCMS) rap = rap + 0.465;
 628       Double_t pT  = trackMC->Pt();
 629       Int_t sign = pdg < 0 ? -1 : 1; // only works for charged pi,K,p !!
 630 //      Double_t transMass = TMath::Sqrt(trackMC->Pt()*trackMC->Pt() + trackMC->GetMass()*trackMC->GetMass()) - trackMC->GetMass();
 631       //
 632       Int_t iPart = -1;
 633       if (TMath::Abs(pdg) == 211)  iPart = 0; // select Pi+/Pi- only
 634       if (TMath::Abs(pdg) == 321)  iPart = 1; // select K+/K- only
 635       if (TMath::Abs(pdg) == 2212) iPart = 2; // select p+/p- only
 636       if (iPart == -1) continue;
 637       //
 638
 639       if (!fSmallTHnSparse){
 640         Double_t vecHistMC[10] = {iPart, centrality,  pT, sign, rap, 0, 1, 0, dxy, 0};
 641         if (!fOnlyQA) fHistMCparticles->Fill(vecHistMC);
 642       }
 643       else{
 644         if (rap>fRapidityCutLow && rap<fRapidityCutHigh){
 645           Double_t vecHistMC[8] = {iPart, centrality,  pT, sign, 1, 0, dxy, 0};
 646           if (!fOnlyQA) fHistMCparticles->Fill(vecHistMC);
 647         }
 648       }
 649     }
 650   }
 651   //
 652   if (!isSelected && !fOnlyQA) {
 653     PostData(1, fListHist);
 654     return;
 655   }
 656   //
 657   if (!vertex || !isVertexOk) {
 658     fHistMult->Fill(-1, processCode);
 659     PostData(1, fListHist);
 660     return;
 661   } else {
 662     if (TMath::Abs(vertex->GetZv()) > 10) {
 663       fHistMult->Fill(-1, processCode);
 664       PostData(1, fListHist);
 665       return;
 666     }
 667   }
 668   //
 669   // count events after physics selection and after vertex selection
 670   //
 671   //cout << "MULTIPLICITY " << trackCounter << " " << fESD->GetEventNumberInFile() <<endl;
 672   fHistMult->Fill(trackCounter, processCode);
 673   fHistCentrality->Fill(centrality);
 674
 675
 676   //***************************************************
 677   // track loop
 678   //***************************************************
 679   //const Float_t kNsigmaCut = 3;
 680   //const Float_t k2sigmaCorr = 1/(0.5*(TMath::Erf(kNsigmaCut/sqrt(2))-TMath::Erf(-kNsigmaCut/sqrt(2))))/*1/0.9545*/;
 681   //
 682   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z for the vertex cut
 683   //
 684   for (Int_t i=0;i<fESD->GetNumberOfTracks();++i) {
 685
 686         AliESDtrack *track = 0;
 687         AliESDtrack *trackForTOF = 0; //normal track for all TOF information needed when using tpconly-tracks
 688
 689         //normal tracks, if tpconly flag is set, use tpconlytracks
 690         if (!fUseTPConlyTracks){
 691         track =fESD->GetTrack(i);
 692         }
 693         else {
 694                 track = fESDtrackCuts->GetTPCOnlyTrack(fESD,i);
 695                 if (!track) continue;
 696                 trackForTOF = fESD->GetTrack(i);
 697         }
 698     //
 699     if (!track->GetInnerParam()) {
 700                 if (fUseTPConlyTracks) {delete track; track = 0;} //need to delete tpconlytrack
 701                 continue;
 702         }
 703     Double_t ptot = track->GetInnerParam()->GetP(); // momentum for dEdx determination
 704     Double_t pT = track->Pt();
 705     track->GetImpactParameters(dca, cov);
 706     //
 707     //
 708     // cut for dead regions in the detector
 709     // if (track->Eta() > 0.1 && (track->Eta() < 0.2 && track->Phi() > 0.1 && track->Phi() < 0.1) continue;
 710     //
 711     // 2.a) apply some standard track cuts according to general recommendations
 712     //
 713     if (!fESDtrackCuts->AcceptTrack(track)) {
 714                 if (fUseTPConlyTracks) {delete track; track = 0;} //need to delete tpconlytrack
 715                 continue;
 716         }
 717
 718         UInt_t status = 0;
 719     if (!fUseTPConlyTracks) status = track->GetStatus();
 720         else status = trackForTOF->GetStatus();
 721     Bool_t hasTOFout  = status&AliESDtrack::kTOFout;
 722     Bool_t hasTOFtime = status&AliESDtrack::kTIME;
 723     Bool_t hasTOFpid  = status&AliESDtrack::kTOFpid;
 724     //Bool_t hasTOFmismatch  = status&AliESDtrack::kTOFmismatch;
 725     Bool_t hasTOF     = kFALSE;
 726     if (hasTOFout && hasTOFtime && hasTOFpid) hasTOF = kTRUE;
 727
 728
 729     //check if TRD contributed to tracking and throw track away if  fTRDinReject flag is set
 730     Bool_t hasTRDin = status&AliESDtrack::kTRDin;
 731     if (hasTRDin && fTRDinReject) {
 732       //hasTOF = kFALSE;
 733       if (fUseTPConlyTracks) {delete track; track = 0;} //need to delete tpconlytrack
 734       continue;
 735     }
 736
 737
 738     //check TOF window
 739     Float_t dxTOF = track->GetTOFsignalDx();
 740     Float_t dzTOF = track->GetTOFsignalDz();
 741
 742     if (hasTOF) fHistTOFwindow->Fill(dxTOF,dzTOF);
 743
 744     //******************************************
 745     //*******NEEDS PROPER CUT SETTER************
 746     //******************************************
 747     //cut on TOF window here
 748     if (TMath::Abs(dxTOF) > fTOFwindow || TMath::Abs(dzTOF) > fTOFwindow) hasTOF = kFALSE;
 749
 750
 751
 752     Float_t length = 0.;
 753     if (!fUseTPConlyTracks) length = track->GetIntegratedLength();
 754     else length = trackForTOF->GetIntegratedLength();
 755
 756     if (length < 350.) hasTOF = kFALSE;
 757
 758     //
 759     // calculate rapidities and kinematics
 760     //
 761     //
 762     Double_t pvec[3];
 763     track->GetPxPyPz(pvec);
 764     Double_t energyPion = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kPion)*AliPID::ParticleMass(AliPID::kPion));
 765     Double_t energyKaon = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kKaon)*AliPID::ParticleMass(AliPID::kKaon));
 766     Double_t energyProton = TMath::Sqrt(track->GetP()*track->GetP() + AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton));
 767     Double_t energyDeuteron = TMath::Sqrt(track->GetP()*track->GetP() + 4*AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton));
 768     //
 769     Double_t rapPion = 0.5*TMath::Log((energyPion + pvec[2])/(energyPion - pvec[2]));
 770     Double_t rapKaon = 0.5*TMath::Log((energyKaon + pvec[2])/(energyKaon - pvec[2]));
 771     Double_t rapProton = 0.5*TMath::Log((energyProton + pvec[2])/(energyProton - pvec[2]));
 772     Double_t rapDeuteron = 0.5*TMath::Log((energyDeuteron + pvec[2])/(energyDeuteron - pvec[2]));
 773
 774     if (fRapCMS) {
 775       rapPion = rapPion + 0.465;
 776       rapKaon = rapKaon + 0.465;
 777       rapProton = rapProton + 0.465;
 778       rapDeuteron = rapDeuteron + 0.465;
 779     }
 780
 781
 782     //
 783 //    Double_t transMassPion = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kPion)*AliPID::ParticleMass(AliPID::kPion))      -  AliPID::ParticleMass(AliPID::kPion);
 784 //    Double_t transMassKaon = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kKaon)*AliPID::ParticleMass(AliPID::kKaon))     -  AliPID::ParticleMass(AliPID::kKaon);
 785  //   Double_t transMassProton = TMath::Sqrt(track->Pt()*track->Pt() + AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))    -  AliPID::ParticleMass(AliPID::kProton);
 786 //    Double_t transMassDeuteron = TMath::Sqrt(track->Pt()*track->Pt() + 4*AliPID::ParticleMass(AliPID::kProton)*AliPID::ParticleMass(AliPID::kProton))    -  2*AliPID::ParticleMass(AliPID::kProton);
 787     //
 788     // 3. make the PID
 789     //
 790     Double_t sign = track->GetSign();
 791     Double_t tpcSignal = track->GetTPCsignal();
 792     //
 793     //
 794     // 3.a. calculate expected signals in nsigma
 795     //
 796     //
 797     // (0.) assumed particle: 0. pion, 1. kaon, 2. proton, 3. deuteron
 798     // (1.) multiplicity or centrality -- number of accepted ESD tracks per events (deprecated), but now classes from 1 to 10, 0: Min. Bias
 799     // (2.) pT
 800     // (3.) sign
 801     // (4.) rapidity --> filled 4x
 802     // (5.)  pull TPC dEx --> filled 4x
 803     // (6.) has valid TOF pid signal
 804     // (7.) nsigma TOF --> filled 4x
 805     // (8..) dca_xy
 806     // (9.) CODE -- only MC 0-generated, 1-true rec. primaries, 2-misident, 3-second weak, 4-second material
 807     //
 808 //    Double_t transMass[4] = {transMassPion,transMassKaon,transMassProton,transMassDeuteron};
 809     Double_t rap[4] = {rapPion,rapKaon,rapProton,rapDeuteron};
 810     Double_t pullsTPC[4] = {fESDpid->NumberOfSigmasTPC(track,AliPID::kPion),
 811                             fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon),
 812                             fESDpid->NumberOfSigmasTPC(track,AliPID::kProton),
 813                             0}; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!
 814
 815
 816     Float_t time0 = fESDpid->GetTOFResponse().GetStartTime(track->P());
 817     //Float_t time0 = fESDpid->GetTOFResponse().GetTimeZero(); //old way of getting time0
 818     //fESDpid->GetTOFResponse().SetTimeResolution(130.);
 819     Double_t pullsTOF[4] ={0.,0.,0.,0.};
 820     if (!fUseTPConlyTracks) {
 821       pullsTOF[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kPion, time0);
 822       pullsTOF[1] = fESDpid->NumberOfSigmasTOF(track,AliPID::kKaon, time0);
 823       pullsTOF[2] = fESDpid->NumberOfSigmasTOF(track,AliPID::kProton, time0);
 824       pullsTOF[3] = 0; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!;
 825     }
 826     else {
 827       pullsTOF[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kPion, time0);
 828       pullsTOF[1] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kKaon, time0);
 829       pullsTOF[2] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kProton, time0);
 830       pullsTOF[3] = 0; // ASK FOR PUTTING THE DEUTERON TO AliPID !!!!!!!!!!!!!!;
 831     }
 832
 833     //
 834 //    Double_t tpcQA[4] = {fESDpid->NumberOfSigmasTPC(track,AliPID::kElectron),
 835 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kPion),
 836 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kKaon),
 837 //                       fESDpid->NumberOfSigmasTPC(track,AliPID::kProton)};
 838
 839     Double_t tofQA[4] = {0.,0.,0.,0.};
 840     if (!fUseTPConlyTracks) {
 841       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kElectron, time0);
 842       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kPion, time0);
 843       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kKaon, time0);
 844       tofQA[0] = fESDpid->NumberOfSigmasTOF(track,AliPID::kProton, time0);
 845     }
 846     else{
 847       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kElectron, time0);
 848       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kPion, time0);
 849       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kKaon, time0);
 850       tofQA[0] = fESDpid->NumberOfSigmasTOF(trackForTOF,AliPID::kProton, time0);
 851     }
 852
 853
 854     //save information for every particle type  // loop over assumed particle type
 855     for(Int_t iPart = 0; iPart < 3; iPart++) {
 856
 857       if (!fSmallTHnSparse) {
 858         Double_t vecHistReal[9]  = {iPart,  centrality,   pT, sign,  rap[iPart], pullsTPC[iPart], hasTOF, pullsTOF[iPart], dca[0]};
 859         if (!fOnlyQA) fHistRealTracks->Fill(vecHistReal);
 860       }
 861       else {
 862         if (pullsTPC[iPart]>fTPCnSigmaCutLow && pullsTPC[iPart]<fTPCnSigmaCutHigh && rap[iPart]>fRapidityCutLow && rap[iPart]<fRapidityCutHigh) {
 863           Double_t vecHistReal[7]  = {iPart,  centrality,   pT, sign, hasTOF, pullsTOF[iPart], dca[0]};
 864           if (!fOnlyQA) fHistRealTracks->Fill(vecHistReal);
 865         }
 866       }
 867
 868
 869       // using MC truth for precise efficiencies...
 870       //
 871       if (fMCtrue && !fOnlyQA) {
 872         Int_t code = 9; // code: 0-generated, 1-true rec. primaries, 2-misident, 3-second weak, 4-second material
 873         Int_t assumedPdg = 0;//2212(proton); 321(Kaon); 211(pion);
 874         Int_t motherCode = -1;
 875         if (iPart == 0) assumedPdg = 211;
 876         if (iPart == 1) assumedPdg = 321;
 877         if (iPart == 2) assumedPdg = 2212;
 878         //
 879         //
 880         TParticle *trackMC = stack->Particle(TMath::Abs(track->GetLabel()));
 881         Int_t pdg = TMath::Abs(trackMC->GetPdgCode());
 882         //
 883         if (pdg != assumedPdg) code = 2;
 884         //if (pdg != assumedPdg && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) code = 5;
 885         if (pdg == assumedPdg && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 1;
 886         if (pdg == assumedPdg && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) {
 887           code = 3;
 888           if (fSaveMotherPDG){
 889             TParticle *trackMother =  stack->Particle(TMath::Abs(trackMC->GetFirstMother()));
 890             if (trackMother->GetPdgCode() == 310) motherCode = 6; //K0
 891             if (trackMother->GetPdgCode() == 3122) motherCode = 7; //Lambda
 892             if (trackMother->GetPdgCode() == 3222) motherCode = 8; //Sigma+
 893           }
 894         }
 895
 896
 897         //FILL MATERIAL TEMPLATE FOR KAONS WITH ELECTRONS
 898         if (iPart != 1){
 899           if (pdg == assumedPdg && stack->IsSecondaryFromMaterial(TMath::Abs(track->GetLabel()))) code = 4;
 900         }
 901         else {
 902           if (pdg == 11) code = 4;
 903           //cout << "got an electron for  kaons!" << endl;
 904         }
 905
 906         //
 907         // muons need special treatment, because they are indistinguishable from pions
 908         //
 909         if (iPart == 0 && pdg == 13  && stack->IsPhysicalPrimary(TMath::Abs(track->GetLabel()))) code = 1;
 910         if (iPart == 0 && pdg == 13  && stack->IsSecondaryFromWeakDecay(TMath::Abs(track->GetLabel()))) code = 3;
 911         //
 912         // check TOF mismatch on MC basis with TOF label
 913         //
 914         Int_t tofLabel[3];
 915         if (!fUseTPConlyTracks) track->GetTOFLabel(tofLabel);
 916         else trackForTOF->GetTOFLabel(tofLabel);
 917
 918
 919         //three options:
 920         //0: do NOT check at all
 921         //1: do check
 922         //2: in case of decays, check if mother label matches --> if yes, assign hasTOF = kTRUE
 923
 924         if (fTOFmisMatch == 0) {
 925           //if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) hasTOF = kFALSE;
 926         }
 927         if (fTOFmisMatch == 1) {
 928           if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) hasTOF = kFALSE;
 929         }
 930         if (fTOFmisMatch == 2) {
 931           if (TMath::Abs(track->GetLabel()) != TMath::Abs(tofLabel[0]) || tofLabel[1] > 0) hasTOF = kFALSE;
 932           TParticle *matchedTrack = stack->Particle(TMath::Abs(tofLabel[0]));
 933           if (TMath::Abs(matchedTrack->GetFirstMother()) == TMath::Abs(track->GetLabel())) hasTOF = kTRUE;
 934         }
 935
 936
 937
 938           //
 939         // IMPORTANT BIG PROBLEM HERE THE PROBABLILITY TO HAVE A PID SIGNAL MUST BE IN !!!!!!!!!!!!
 940         //
 941         if (!fSmallTHnSparse){
 942           Double_t vectorHistMC[10] = {iPart,  centrality,  pT, sign,  rap[iPart], pullsTPC[iPart], hasTOF, pullsTOF[iPart], dca[0], code};
 943           if (!fOnlyQA) {
 944             fHistMCparticles->Fill(vectorHistMC);
 945             if (motherCode != -1 && fSaveMotherPDG) { //if mother of weak decay is K0, lambda or sigma+ add track again with this information
 946               Double_t vectorHistMCmother[10] = {iPart,  centrality,  pT, sign,  rap[iPart], pullsTPC[iPart], hasTOF, pullsTOF[iPart], dca[0], motherCode};
 947               fHistMCparticles->Fill(vectorHistMCmother);
 948             }
 949           }
 950         }
 951         else{
 952           if (pullsTPC[iPart]>fTPCnSigmaCutLow && pullsTPC[iPart]<fTPCnSigmaCutHigh && rap[iPart]>fRapidityCutLow && rap[iPart]<fRapidityCutHigh) {
 953             //                              0,           1,   2,    3,           4,               5,      6,               7,      8,   9
 954             Double_t vectorHistMC[8] = {iPart,  centrality,  pT, sign, hasTOF, pullsTOF[iPart], dca[0], code};
 955             if (!fOnlyQA) {
 956               fHistMCparticles->Fill(vectorHistMC);
 957               if (motherCode != -1 && fSaveMotherPDG) { //if mother of weak decay is K0, lambda or sigma+ add track again with this information
 958                 Double_t vectorHistMCmother[8] = {iPart,  centrality,  pT, sign, hasTOF, pullsTOF[iPart], dca[0], motherCode};
 959                 fHistMCparticles->Fill(vectorHistMCmother);
 960               }
 961             }
 962           }
 963         }
 964       }
 965       //
 966       //
 967       Int_t tpcShared = track->GetTPCnclsS();
 968       if (TMath::Abs(track->Eta()) < 0.8 && iPart == 0 && tpcShared < 4) fHistPidQA->Fill(ptot,tpcSignal,sign);
 969     } // end loop over assumed particle type
 970
 971           //need to delete tpconlytrack
 972           if (fUseTPConlyTracks){
 973                 delete track;
 974                 track = 0;
 975           }
 976
 977   } // end of track loop
 978
 979   // Post output data
 980   PostData(1, fListHist);
 981
 982 }
 983
 984
 985 //________________________________________________________________________
 986 void AliAnalysisTPCTOFpA::Terminate(Option_t *)
 987 {
 988   // Draw result to the screen
 989   // Called once at the end of the query
 990   Printf("*** CONSTRUCTOR CALLED ****");
 991
 992 }
 993
 994
 995 //________________________________________________________________________
 996 Bool_t AliAnalysisTPCTOFpA::SelectOnImpPar(AliESDtrack* t) {
 997   //
 998   // cut on transverse impact parameter // DEPRECATED
 999   //
1000   Float_t d0z0[2],covd0z0[3];
1001   t->GetImpactParameters(d0z0,covd0z0);
1002   Float_t sigma= 0.0050+0.0060/TMath::Power(t->Pt(),0.9);
1003   Float_t d0max = 7.*sigma;
1004   //
1005   Float_t sigmaZ = 0.0146+0.0070/TMath::Power(t->Pt(),1.114758);
1006   if (t->Pt() > 1) sigmaZ = 0.0216;
1007   Float_t d0maxZ = 5.*sigmaZ;
1008   //
1009   if(TMath::Abs(d0z0[0]) < d0max && TMath::Abs(d0z0[1]) < d0maxZ) return kTRUE;
1010   return kFALSE;
1011 }
1012
1013
1014 //________________________________________________________________________
1015 void AliAnalysisTPCTOFpA::BinLogAxis(const TH1 *h) {
1016   //
1017   // Method for the correct logarithmic binning of histograms
1018   //
1019   TAxis *axis = h->GetXaxis();
1020   int bins = axis->GetNbins();
1021
1022   Double_t from = axis->GetXmin();
1023   Double_t to = axis->GetXmax();
1024   Double_t *newBins = new Double_t[bins + 1];
1025
1026   newBins[0] = from;
1027   Double_t factor = pow(to/from, 1./bins);
1028
1029   for (int i = 1; i <= bins; i++) {
1030    newBins[i] = factor * newBins[i-1];
1031   }
1032   axis->Set(bins, newBins);
1033   delete [] newBins;
1034
1035 }
1036
1037
1038 //________________________________________________________________________
1039 Int_t AliAnalysisTPCTOFpA::GetPythiaEventProcessType(const AliHeader* aHeader, const Bool_t adebug) const {
1040   //
1041   // get the process type of the event.
1042   //
1043
1044   // can only read pythia headers, either directly or from cocktalil header
1045   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader());
1046
1047   if (!pythiaGenHeader) {
1048
1049     AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader());
1050     if (!genCocktailHeader) {
1051       //printf("AliAnalysisTPCTOFpA::GetProcessType : Unknown header type (not Pythia or Cocktail). \n");
1052       return -1;
1053     }
1054
1055     TList* headerList = genCocktailHeader->GetHeaders();
1056     if (!headerList) {
1057       return -1;
1058     }
1059
1060     for (Int_t i=0; i<headerList->GetEntries(); i++) {
1061       pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i));
1062       if (pythiaGenHeader)
1063         break;
1064     }
1065
1066     if (!pythiaGenHeader) {
1067       //printf("AliAnalysisTPCTOFpA::GetProcessType : Could not find Pythia header. \n");
1068       return -1;
1069     }
1070   }
1071
1072   if (adebug) {
1073     //printf("AliAnalysisTPCTOFpA::GetProcessType : Pythia process type found: %d \n",pythiaGenHeader->ProcessType());
1074   }
1075
1076   return pythiaGenHeader->ProcessType();
1077 }