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