PWGGA/CaloTrackCorrelations/AliAnaChargedParticles.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, provided 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 purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //_________________________________________________________________________
  17 //
  18 // Class for track selection and identification (not done now)
  19 // Tracks from the CTS are kept in the AOD.
  20 // Few histograms produced.
  21 //
  22 //-- Author: Gustavo Conesa (INFN-LNF)
  23 //_________________________________________________________________________
  24
  25
  26 // --- ROOT system ---
  27 #include "TParticle.h"
  28 #include "TH2F.h"
  29 #include "TDatabasePDG.h"
  30
  31 //---- AliRoot system ----
  32 #include "AliAnaChargedParticles.h"
  33 #include "AliCaloTrackReader.h"
  34 #include "AliAODPWG4Particle.h"
  35 #include "AliStack.h"
  36 #include "AliFiducialCut.h"
  37 #include "AliVTrack.h"
  38 #include "AliAODMCParticle.h"
  39 #include "AliAODTrack.h"
  40 #include "AliAODEvent.h"
  41 #include "AliESDEvent.h"
  42
  43 ClassImp(AliAnaChargedParticles)
  44
  45 //__________________________________________________
  46   AliAnaChargedParticles::AliAnaChargedParticles() :
  47     AliAnaCaloTrackCorrBaseClass(),
  48     fFillTrackBCHistograms(0), fFillVertexBC0Histograms(0),
  49     fMomentum(),
  50     //Histograms
  51     fhNtracks(0),      fhPt(0),            fhPtNoCut(0),
  52     fhPtCutDCA(0),     fhPtCutDCABCOK(0),
  53     fhPhiNeg(0),       fhEtaNeg(0),
  54     fhPhiPos(0),       fhEtaPos(0),
  55     fhEtaPhiPos(0),    fhEtaPhiNeg(0),
  56     fhPtVtxOutBC0(0),  fhEtaPhiVtxOutBC0(0),
  57     fhPtVtxInBC0(0),   fhEtaPhiVtxInBC0(0),
  58     fhPtSPDRefit(0),         fhPtNoSPDRefit(0),         fhPtNoSPDNoRefit(0),
  59     fhEtaPhiSPDRefitPt02(0), fhEtaPhiNoSPDRefitPt02(0), fhEtaPhiNoSPDNoRefitPt02(0),
  60     fhEtaPhiSPDRefitPt3(0),  fhEtaPhiNoSPDRefitPt3(0),  fhEtaPhiNoSPDNoRefitPt3(0),
  61
  62     //TOF
  63     fhTOFSignal(0),    fhTOFSignalPtCut(0),  fhTOFSignalBCOK(0),
  64     fhPtTOFSignal(0),  fhPtTOFSignalDCACut(0),
  65     fhPtTOFSignalVtxOutBC0(0), fhPtTOFSignalVtxInBC0(0),
  66     fhPtTOFStatus0(0), fhEtaPhiTOFStatus0(0),
  67     fhEtaPhiTOFBC0(0), fhEtaPhiTOFBCPlus(0), fhEtaPhiTOFBCMinus(0),
  68     fhEtaPhiTOFBC0PileUpSPD(0),
  69     fhEtaPhiTOFBCPlusPileUpSPD(0),
  70     fhEtaPhiTOFBCMinusPileUpSPD(0),
  71     fhProductionVertexBC(0),
  72     fhPtNPileUpSPDVtx(0),    fhPtNPileUpTrkVtx(0),
  73     fhPtNPileUpSPDVtxBC0(0), fhPtNPileUpTrkVtxBC0(0)
  74
  75 {
  76   //Default Ctor
  77
  78   for(Int_t i = 0; i < 7; i++)
  79   {
  80     fhPtPileUp         [i] = 0;
  81     fhPtTOFSignalPileUp[i] = 0;
  82     fhPtTOFSignalVtxOutBC0PileUp[i] = 0;
  83     fhPtTOFSignalVtxInBC0PileUp [i] = 0;
  84     fhProductionVertexBCPileUp  [i] = 0;
  85   }
  86
  87   for(Int_t i = 0; i < 3; i++)
  88   {
  89     fhPtDCA               [i] = 0 ;
  90
  91     fhPtDCASPDRefit       [i] = 0 ;
  92     fhPtDCANoSPDRefit     [i] = 0 ;
  93     fhPtDCANoSPDNoRefit   [i] = 0 ;
  94
  95     fhPtDCAPileUp         [i] = 0 ;
  96     fhPtDCATOFBC0         [i] = 0 ;
  97     fhPtDCATOFBCOut       [i] = 0 ;
  98     fhPtDCAPileUpTOFBC0   [i] = 0 ;
  99     fhPtDCANoTOFHit       [i] = 0 ;
 100     fhPtDCAPileUpNoTOFHit [i] = 0 ;
 101
 102     fhPtDCAVtxOutBC0      [i] = 0 ;
 103     fhPtDCAVtxInBC0       [i] = 0 ;
 104     fhPtDCAVtxOutBC0PileUp[i] = 0 ;
 105     fhPtDCAVtxInBC0PileUp [i] = 0 ;
 106     fhPtDCAVtxOutBC0NoTOFHit[i] = 0 ;
 107     fhPtDCAVtxInBC0NoTOFHit [i] = 0 ;
 108     fhPtDCAVtxOutBC0PileUpNoTOFHit[i] = 0 ;
 109     fhPtDCAVtxInBC0PileUpNoTOFHit [i] = 0 ;
 110   }
 111
 112   // MC
 113   for(Int_t imcPart = 0; imcPart < 6; imcPart++)
 114   {
 115     fhPtMCPart     [imcPart] = 0;
 116     fhPtMCPrimPart [imcPart] = 0;
 117     fhPhiMCPart    [imcPart] = 0;
 118     fhPhiMCPrimPart[imcPart] = 0;
 119     fhEtaMCPart    [imcPart] = 0;
 120     fhEtaMCPrimPart[imcPart] = 0;
 121   }
 122
 123   //Initialize parameters
 124   InitParameters();
 125
 126 }
 127
 128 //__________________________________________________
 129 void AliAnaChargedParticles::FillPrimaryHistograms()
 130 {
 131   // Fill primary generated particles histograms if MC data is available
 132
 133   Int_t    pdg       =  0 ;
 134   Int_t    nprim     =  0 ;
 135
 136   TParticle        * primStack = 0;
 137   AliAODMCParticle * primAOD   = 0;
 138
 139   // Get the ESD MC particles container
 140   AliStack * stack = 0;
 141   if( GetReader()->ReadStack() )
 142   {
 143     stack = GetMCStack();
 144     if(!stack ) return;
 145     nprim = stack->GetNtrack();
 146   }
 147
 148   // Get the AOD MC particles container
 149   TClonesArray * mcparticles = 0;
 150   if( GetReader()->ReadAODMCParticles() )
 151   {
 152     mcparticles = GetReader()->GetAODMCParticles();
 153     if( !mcparticles ) return;
 154     nprim = mcparticles->GetEntriesFast();
 155   }
 156
 157   for(Int_t i=0 ; i < nprim; i++)
 158   {
 159     if(GetReader()->AcceptOnlyHIJINGLabels() && !GetReader()->IsHIJINGLabel(i)) continue ;
 160
 161     if(GetReader()->ReadStack())
 162     {
 163       primStack = stack->Particle(i) ;
 164       if(!primStack)
 165       {
 166         printf("AliAnaChargedParticles::FillPrimaryMCHistograms() - ESD primaries pointer not available!!\n");
 167         continue;
 168       }
 169
 170       if( primStack->GetStatusCode() != 1 ) continue;
 171
 172       Int_t charge = (Int_t )TDatabasePDG::Instance()->GetParticle(primStack->GetPdgCode())->Charge();
 173       if( TMath::Abs(charge) == 0 ) continue;
 174
 175       pdg  = TMath::Abs(primStack->GetPdgCode());
 176
 177       if(primStack->Energy() == TMath::Abs(primStack->Pz()))  continue ; //Protection against floating point exception
 178
 179       //printf("i %d, %s %d  %s %d \n",i, stack->Particle(i)->GetName(), stack->Particle(i)->GetPdgCode(),
 180       //       prim->GetName(), prim->GetPdgCode());
 181
 182       //Charged kinematics
 183       primStack->Momentum(fMomentum);
 184     }
 185     else
 186     {
 187       primAOD = (AliAODMCParticle *) mcparticles->At(i);
 188       if(!primAOD)
 189       {
 190         printf("AliAnaChargedParticles::FillPrimaryHistograms() - AOD primaries pointer not available!!\n");
 191         continue;
 192       }
 193
 194       if( primAOD->GetStatus() != 1 ) continue;
 195
 196       if(TMath::Abs(primAOD->Charge()) == 0 ) continue;
 197
 198       pdg = TMath::Abs(primAOD->GetPdgCode());
 199
 200       if(primAOD->E() == TMath::Abs(primAOD->Pz()))  continue ; //Protection against floating point exception
 201
 202       //Charged kinematics
 203       fMomentum.SetPxPyPzE(primAOD->Px(),primAOD->Py(),primAOD->Pz(),primAOD->E());
 204     }
 205
 206     Int_t mcType = kmcUnknown;
 207     if     (pdg==211 ) mcType = kmcPion;
 208     else if(pdg==2212) mcType = kmcProton;
 209     else if(pdg==321 ) mcType = kmcKaon;
 210     else if(pdg==11  ) mcType = kmcElectron;
 211     else if(pdg==13  ) mcType = kmcMuon;
 212
 213     //printf("pdg %d, charge %f, mctype %d\n",pdg, charge, mcType);
 214
 215     Float_t  ptPrim = fMomentum.Pt();
 216     Float_t etaPrim = fMomentum.Eta();
 217     Float_t phiPrim = fMomentum.Phi();
 218     if(phiPrim < 0) phiPrim+=TMath::TwoPi();
 219
 220     Bool_t in = GetFiducialCut()->IsInFiducialCut(etaPrim,phiPrim,kCTS) ;
 221
 222     if(in)
 223       fhPtMCPrimPart[mcType]->Fill(ptPrim);
 224     fhEtaMCPrimPart[mcType]->Fill(ptPrim,etaPrim);
 225     fhPhiMCPrimPart[mcType]->Fill(ptPrim,phiPrim);
 226   }
 227 }
 228
 229 //_______________________________________________________
 230 TList *  AliAnaChargedParticles::GetCreateOutputObjects()
 231 {
 232   // Create histograms to be saved in output file and
 233   // store them in fOutputContainer
 234
 235
 236   TList * outputContainer = new TList() ;
 237   outputContainer->SetName("ChargedParticleHistos") ;
 238
 239   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins(); Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
 240   Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();  Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
 241   Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();  Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();  Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
 242
 243   fhNtracks  = new TH1F ("hNtracks","# of tracks", 1000,0,1000);
 244   fhNtracks->SetXTitle("# of tracks");
 245   outputContainer->Add(fhNtracks);
 246
 247   fhPt  = new TH1F ("hPt","#it{p}_{T} distribution", nptbins,ptmin,ptmax);
 248   fhPt->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 249   outputContainer->Add(fhPt);
 250
 251   fhPtNoCut  = new TH1F ("hPtNoCut","#it{p}_{T} distribution, raw tracks", nptbins,ptmin,ptmax);
 252   fhPtNoCut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 253   outputContainer->Add(fhPtNoCut);
 254
 255   if(!GetReader()->IsDCACutOn())
 256   {
 257     fhPtCutDCA  = new TH1F ("hPtCutDCA","#it{p}_{T} distribution, cut DCA", nptbins,ptmin,ptmax);
 258     fhPtCutDCA->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 259     outputContainer->Add(fhPtCutDCA);
 260   }
 261
 262   if(fFillVertexBC0Histograms && !GetReader()->IsDCACutOn())
 263   {
 264     fhPtCutDCABCOK  = new TH1F ("hPtCutDCABCOK","#it{p}_{T} distribution, DCA cut, track BC=0 or -100", nptbins,ptmin,ptmax);
 265     fhPtCutDCABCOK->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 266     outputContainer->Add(fhPtCutDCABCOK);
 267   }
 268
 269   fhPhiNeg  = new TH2F ("hPhiNegative","#phi of negative charges distribution",
 270                         nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 271   fhPhiNeg->SetYTitle("#phi (rad)");
 272   fhPhiNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 273   outputContainer->Add(fhPhiNeg);
 274
 275   fhEtaNeg  = new TH2F ("hEtaNegative","#eta of negative charges distribution",
 276                         nptbins,ptmin,ptmax, netabins,etamin,etamax);
 277   fhEtaNeg->SetYTitle("#eta ");
 278   fhEtaNeg->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 279   outputContainer->Add(fhEtaNeg);
 280
 281   fhPhiPos  = new TH2F ("hPhiPositive","#phi of positive charges distribution",
 282                         nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 283   fhPhiPos->SetYTitle("#phi (rad)");
 284   fhPhiPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 285   outputContainer->Add(fhPhiPos);
 286
 287   fhEtaPos  = new TH2F ("hEtaPositive","#eta of positive charges distribution",
 288                         nptbins,ptmin,ptmax, netabins,etamin,etamax);
 289   fhEtaPos->SetYTitle("#eta ");
 290   fhEtaPos->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 291   outputContainer->Add(fhEtaPos);
 292
 293   fhEtaPhiPos  = new TH2F ("hEtaPhiPositive","pt/eta/phi of positive charge",netabins,etamin,etamax, nphibins,phimin,phimax);
 294   fhEtaPhiPos->SetXTitle("#eta ");
 295   fhEtaPhiPos->SetYTitle("#phi (rad)");
 296   outputContainer->Add(fhEtaPhiPos);
 297
 298   fhEtaPhiNeg  = new TH2F ("hEtaPhiNegative","#eta vs #phi of negative charge",netabins,etamin,etamax, nphibins,phimin,phimax);
 299   fhEtaPhiNeg->SetXTitle("#eta ");
 300   fhEtaPhiNeg->SetYTitle("#phi (rad)");
 301   outputContainer->Add(fhEtaPhiNeg);
 302
 303   if(fFillVertexBC0Histograms)
 304   {
 305     fhPtVtxOutBC0  = new TH1F ("hPtVtxOutBC0","#it{p}_{T} distribution, vertex in BC=0", nptbins,ptmin,ptmax);
 306     fhPtVtxOutBC0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 307     outputContainer->Add(fhPtVtxOutBC0);
 308
 309     fhEtaPhiVtxOutBC0  = new TH2F ("hEtaPhiVtxOutBC0","#eta vs #phi of all charges with vertex in BC=0",netabins,etamin,etamax, nphibins,phimin,phimax);
 310     fhEtaPhiVtxOutBC0->SetXTitle("#eta ");
 311     fhEtaPhiVtxOutBC0->SetYTitle("#phi (rad)");
 312     outputContainer->Add(fhEtaPhiVtxOutBC0);
 313
 314     fhPtVtxInBC0  = new TH1F ("hPtVtxInBC0","#it{p}_{T} distribution, vertex in BC=0", nptbins,ptmin,ptmax);
 315     fhPtVtxInBC0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 316     outputContainer->Add(fhPtVtxInBC0);
 317
 318     fhEtaPhiVtxInBC0  = new TH2F ("hEtaPhiVtxInBC0","#eta vs #phi of all charges with vertex in BC=0",netabins,etamin,etamax, nphibins,phimin,phimax);
 319     fhEtaPhiVtxInBC0->SetXTitle("#eta ");
 320     fhEtaPhiVtxInBC0->SetYTitle("#phi (rad)");
 321     outputContainer->Add(fhEtaPhiVtxInBC0);
 322   }
 323
 324   fhPtSPDRefit  = new TH1F ("hPtSPDRefit","#it{p}_{T} distribution of tracks with SPD and ITS refit", nptbins,ptmin,ptmax);
 325   fhPtSPDRefit->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 326   outputContainer->Add(fhPtSPDRefit);
 327
 328   fhEtaPhiSPDRefitPt02  = new TH2F ("hEtaPhiSPDRefitPt02","#eta vs #phi of tracks with SPD and ITS refit, #it{p}_{T}< 2 GeV/#it{c}",
 329                                     netabins,etamin,etamax, nphibins,phimin,phimax);
 330   fhEtaPhiSPDRefitPt02->SetXTitle("#eta ");
 331   fhEtaPhiSPDRefitPt02->SetYTitle("#phi (rad)");
 332   outputContainer->Add(fhEtaPhiSPDRefitPt02);
 333
 334   fhEtaPhiSPDRefitPt3  = new TH2F ("hEtaPhiSPDRefitPt3","#eta vs #phi of tracks with SPD and ITS refit, #it{p}_{T}> 3 GeV/#it{c}",
 335                                    netabins,etamin,etamax, nphibins,phimin,phimax);
 336   fhEtaPhiSPDRefitPt3->SetXTitle("#eta ");
 337   fhEtaPhiSPDRefitPt3->SetYTitle("#phi (rad)");
 338   outputContainer->Add(fhEtaPhiSPDRefitPt3);
 339
 340   fhPtNoSPDRefit  = new TH1F ("hPtNoSPDRefit","#it{p}_{T} distribution of constrained tracks no SPD and with ITSRefit",
 341                               nptbins,ptmin,ptmax);
 342   fhPtNoSPDRefit->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 343   outputContainer->Add(fhPtNoSPDRefit);
 344
 345   fhEtaPhiNoSPDRefitPt02  = new TH2F ("hEtaPhiNoSPDRefitPt02","#eta vs #phi of constrained tracks no SPD and with ITSRefit, #it{p}_{T}< 2 GeV/#it{c}",
 346                                       netabins,etamin,etamax, nphibins,phimin,phimax);
 347   fhEtaPhiNoSPDRefitPt02->SetXTitle("#eta ");
 348   fhEtaPhiNoSPDRefitPt02->SetYTitle("#phi (rad)");
 349   outputContainer->Add(fhEtaPhiNoSPDRefitPt02);
 350
 351   fhEtaPhiNoSPDRefitPt3  = new TH2F ("hEtaPhiNoSPDRefitPt3","#eta vs #phi of of constrained tracks no SPD and with ITSRefit, #it{p}_{T}> 3 GeV/#it{c}",
 352                                      netabins,etamin,etamax, nphibins,phimin,phimax);
 353   fhEtaPhiNoSPDRefitPt3->SetXTitle("#eta ");
 354   fhEtaPhiNoSPDRefitPt3->SetYTitle("#phi (rad)");
 355   outputContainer->Add(fhEtaPhiNoSPDRefitPt3);
 356
 357   fhPtNoSPDNoRefit  = new TH1F ("hPtNoSPDNoRefit","#it{p}_{T} distribution of constrained tracks with no SPD requierement and without ITSRefit",
 358                                 nptbins,ptmin,ptmax);
 359   fhPtNoSPDNoRefit->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 360   outputContainer->Add(fhPtNoSPDNoRefit);
 361
 362   fhEtaPhiNoSPDNoRefitPt02  = new TH2F ("hEtaPhiNoSPDNoRefitPt02",
 363                                         "#eta vs #phi of constrained tracks with no SPD requierement and without ITSRefit, #it{p}_{T}< 2 GeV/#it{c}",
 364                                         netabins,etamin,etamax, nphibins,phimin,phimax);
 365   fhEtaPhiNoSPDNoRefitPt02->SetXTitle("#eta ");
 366   fhEtaPhiNoSPDNoRefitPt02->SetYTitle("#phi (rad)");
 367   outputContainer->Add(fhEtaPhiNoSPDNoRefitPt02);
 368
 369   fhEtaPhiNoSPDNoRefitPt3  = new TH2F ("hEtaPhiNoSPDNoRefitPt3",
 370                                        "#eta vs #phi of constrained tracks with no SPD requierement and without ITSRefit, #it{p}_{T}> 3 GeV/#it{c}",
 371                                        netabins,etamin,etamax, nphibins,phimin,phimax);
 372   fhEtaPhiNoSPDNoRefitPt3->SetXTitle("#eta ");
 373   fhEtaPhiNoSPDNoRefitPt3->SetYTitle("#phi (rad)");
 374   outputContainer->Add(fhEtaPhiNoSPDNoRefitPt3);
 375
 376   if(fFillVertexBC0Histograms)
 377   {
 378     fhProductionVertexBC      = new TH1F("hProductionVertexBC", "tracks production vertex bunch crossing ", 41 , -20 , 20 ) ;
 379     fhProductionVertexBC->SetYTitle("# tracks");
 380     fhProductionVertexBC->SetXTitle("Bunch crossing");
 381     outputContainer->Add(fhProductionVertexBC);
 382   }
 383
 384   Int_t ntofbins = 1000;
 385   Int_t mintof = -500;
 386   Int_t maxtof =  500;
 387
 388   fhTOFSignal  = new TH1F ("hTOFSignal","TOF signal", ntofbins,mintof,maxtof);
 389   fhTOFSignal->SetXTitle("TOF signal (ns)");
 390   outputContainer->Add(fhTOFSignal);
 391
 392   if(fFillTrackBCHistograms)
 393   {
 394     fhTOFSignalBCOK  = new TH1F ("hTOFSignalBCOK","TOF signal", ntofbins,mintof,maxtof);
 395     fhTOFSignalBCOK->SetXTitle("TOF signal (ns)");
 396     outputContainer->Add(fhTOFSignalBCOK);
 397   }
 398
 399   fhTOFSignalPtCut  = new TH1F ("hTOFSignalPtCut","TOF signal", ntofbins,mintof,maxtof);
 400   fhTOFSignalPtCut->SetXTitle("TOF signal (ns)");
 401   outputContainer->Add(fhTOFSignalPtCut);
 402
 403   fhPtTOFSignal  = new TH2F ("hPtTOFSignal","TOF signal", nptbins,ptmin,ptmax,ntofbins,mintof,maxtof);
 404   fhPtTOFSignal->SetYTitle("TOF signal (ns)");
 405   fhPtTOFSignal->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 406   outputContainer->Add(fhPtTOFSignal);
 407
 408   if(!GetReader()->IsDCACutOn())
 409   {
 410     fhPtTOFSignalDCACut  = new TH2F ("hPtTOFSignalDCACut","TOF signal after DCA cut", nptbins,ptmin,ptmax,ntofbins,mintof,maxtof);
 411     fhPtTOFSignalDCACut->SetYTitle("TOF signal (ns)");
 412     fhPtTOFSignalDCACut->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 413     outputContainer->Add(fhPtTOFSignalDCACut);
 414   }
 415
 416   if(fFillVertexBC0Histograms)
 417   {
 418     fhPtTOFSignalVtxOutBC0  = new TH2F ("hPtTOFSignalVtxOutBC0","TOF signal, vtx BC!=0", nptbins,ptmin,ptmax,ntofbins,mintof,maxtof);
 419     fhPtTOFSignalVtxOutBC0->SetYTitle("TOF signal (ns)");
 420     fhPtTOFSignalVtxOutBC0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 421     outputContainer->Add(fhPtTOFSignalVtxOutBC0);
 422
 423     fhPtTOFSignalVtxInBC0  = new TH2F ("hPtTOFSignalVtxInBC0","TOF signal, vtx BC=0", nptbins,ptmin,ptmax,ntofbins,mintof,maxtof);
 424     fhPtTOFSignalVtxInBC0->SetYTitle("TOF signal (ns)");
 425     fhPtTOFSignalVtxInBC0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 426     outputContainer->Add(fhPtTOFSignalVtxInBC0);
 427   }
 428
 429   if(IsPileUpAnalysisOn())
 430   {
 431     TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
 432
 433     for(Int_t i = 0 ; i < 7 ; i++)
 434     {
 435       fhPtPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
 436                                 Form("Track #it{p}_{T}distribution, %s Pile-Up event",pileUpName[i].Data()),
 437                                 nptbins,ptmin,ptmax);
 438       fhPtPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 439       outputContainer->Add(fhPtPileUp[i]);
 440
 441       fhPtTOFSignalPileUp[i]  = new TH2F(Form("hPtTOFSignalPileUp%s",pileUpName[i].Data()),
 442                                          Form("Track TOF vs #it{p}_{T}distribution, %s Pile-Up event",pileUpName[i].Data()),
 443                                          nptbins,ptmin,ptmax,ntofbins,mintof,maxtof);
 444       fhPtTOFSignalPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 445       fhPtTOFSignalPileUp[i]->SetXTitle("TOF signal (ns)");
 446       outputContainer->Add(fhPtTOFSignalPileUp[i]);
 447
 448       if(fFillVertexBC0Histograms)
 449       {
 450         fhPtTOFSignalVtxOutBC0PileUp[i]  = new TH2F(Form("hPtTOFSignalVtxOutBC0PileUp%s",pileUpName[i].Data()),
 451                                            Form("Track TOF vs #it{p}_{T}distribution, %s Pile-Up event, vtx BC!=0",pileUpName[i].Data()),
 452                                            nptbins,ptmin,ptmax,ntofbins,mintof,maxtof);
 453         fhPtTOFSignalVtxOutBC0PileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 454         fhPtTOFSignalVtxOutBC0PileUp[i]->SetXTitle("TOF signal (ns)");
 455         outputContainer->Add(fhPtTOFSignalVtxOutBC0PileUp[i]);
 456
 457         fhPtTOFSignalVtxInBC0PileUp[i]  = new TH2F(Form("hPtTOFSignalVtxInBC0PileUp%s",pileUpName[i].Data()),
 458                                                     Form("Track TOF vs #it{p}_{T}distribution, %s Pile-Up event, vtx BC=0",pileUpName[i].Data()),
 459                                                     nptbins,ptmin,ptmax,ntofbins,mintof,maxtof);
 460         fhPtTOFSignalVtxInBC0PileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 461         fhPtTOFSignalVtxInBC0PileUp[i]->SetXTitle("TOF signal (ns)");
 462         outputContainer->Add(fhPtTOFSignalVtxInBC0PileUp[i]);
 463       }
 464
 465       if(fFillVertexBC0Histograms)
 466       {
 467         fhProductionVertexBCPileUp[i]      = new TH1F(Form("hProductionVertexBCPileUp%s",pileUpName[i].Data()),
 468                                                 Form("tracks production vertex bunch crossing, %s Pile-Up event",pileUpName[i].Data()),
 469                                                 41 , -20 , 20 ) ;
 470         fhProductionVertexBCPileUp[i]->SetYTitle("# tracks");
 471         fhProductionVertexBCPileUp[i]->SetXTitle("Bunch crossing");
 472         outputContainer->Add(fhProductionVertexBCPileUp[i]);
 473       }
 474     }
 475
 476
 477     if(fFillTrackBCHistograms)
 478     {
 479       fhEtaPhiTOFBC0  = new TH2F ("hEtaPhiTOFBC0","eta-phi for tracks with hit on TOF, and tof corresponding to BC=0",netabins,etamin,etamax, nphibins,phimin,phimax);
 480       fhEtaPhiTOFBC0->SetXTitle("#eta ");
 481       fhEtaPhiTOFBC0->SetYTitle("#phi (rad)");
 482       outputContainer->Add(fhEtaPhiTOFBC0);
 483
 484       fhEtaPhiTOFBCPlus  = new TH2F ("hEtaPhiTOFBCPlus","eta-phi for tracks with hit on TOF, and tof corresponding to BC>0",netabins,etamin,etamax, nphibins,phimin,phimax);
 485       fhEtaPhiTOFBCPlus->SetXTitle("#eta ");
 486       fhEtaPhiTOFBCPlus->SetYTitle("#phi (rad)");
 487       outputContainer->Add(fhEtaPhiTOFBCPlus);
 488
 489       fhEtaPhiTOFBCMinus  = new TH2F ("hEtaPhiTOFBCMinus","eta-phi for tracks with hit on TOF, and tof corresponding to BC<0",netabins,etamin,etamax, nphibins,phimin,phimax);
 490       fhEtaPhiTOFBCMinus->SetXTitle("#eta ");
 491       fhEtaPhiTOFBCMinus->SetYTitle("#phi (rad)");
 492       outputContainer->Add(fhEtaPhiTOFBCMinus);
 493
 494       if(IsPileUpAnalysisOn())
 495       {
 496         fhEtaPhiTOFBC0PileUpSPD  = new TH2F ("hEtaPhiTOFBC0PileUpSPD","eta-phi for tracks with hit on TOF, and tof corresponding to BC=0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
 497         fhEtaPhiTOFBC0PileUpSPD->SetXTitle("#eta ");
 498         fhEtaPhiTOFBC0PileUpSPD->SetYTitle("#phi (rad)");
 499         outputContainer->Add(fhEtaPhiTOFBC0PileUpSPD);
 500
 501         fhEtaPhiTOFBCPlusPileUpSPD  = new TH2F ("hEtaPhiTOFBCPlusPileUpSPD","eta-phi for tracks with hit on TOF, and tof corresponding to BC>0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
 502         fhEtaPhiTOFBCPlusPileUpSPD->SetXTitle("#eta ");
 503         fhEtaPhiTOFBCPlusPileUpSPD->SetYTitle("#phi (rad)");
 504         outputContainer->Add(fhEtaPhiTOFBCPlusPileUpSPD);
 505
 506         fhEtaPhiTOFBCMinusPileUpSPD  = new TH2F ("hEtaPhiTOFBCMinusPileUpSPD","eta-phi for tracks with hit on TOF, and tof corresponding to BC<0, SPD pile-up",netabins,etamin,etamax, nphibins,phimin,phimax);
 507         fhEtaPhiTOFBCMinusPileUpSPD->SetXTitle("#eta ");
 508         fhEtaPhiTOFBCMinusPileUpSPD->SetYTitle("#phi (rad)");
 509         outputContainer->Add(fhEtaPhiTOFBCMinusPileUpSPD);
 510       }
 511     }
 512
 513   }
 514
 515   fhPtTOFStatus0  = new TH1F ("hPtTOFStatus0","#it{p}_{T} distribution of tracks not hitting TOF", nptbins,ptmin,ptmax);
 516   fhPtTOFStatus0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 517   outputContainer->Add(fhPtTOFStatus0);
 518
 519
 520   fhEtaPhiTOFStatus0  = new TH2F ("hEtaPhiTOFStatus0","eta-phi for tracks without hit on TOF",netabins,etamin,etamax, nphibins,phimin,phimax);
 521   fhEtaPhiTOFStatus0->SetXTitle("#eta ");
 522   fhEtaPhiTOFStatus0->SetYTitle("#phi (rad)");
 523   outputContainer->Add(fhEtaPhiTOFStatus0);
 524
 525   TString dcaName[] = {"xy","z","Cons"} ;
 526   Int_t ndcabins = 800;
 527   Int_t mindca = -4;
 528   Int_t maxdca =  4;
 529
 530   for(Int_t i = 0 ; i < 3 ; i++)
 531   {
 532
 533     fhPtDCA[i]  = new TH2F(Form("hPtDCA%s",dcaName[i].Data()),
 534                            Form("Track DCA%s vs #it{p}_{T}distribution",dcaName[i].Data()),
 535                            nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 536     fhPtDCA[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 537     fhPtDCA[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 538     outputContainer->Add(fhPtDCA[i]);
 539
 540     fhPtDCASPDRefit[i]  = new TH2F(Form("hPtDCA%sSPDRefit",dcaName[i].Data()),
 541                                         Form("Track DCA%s vs #it{p}_{T}distribution of tracks with SPD and ITS refit",dcaName[i].Data()),
 542                                         nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 543     fhPtDCASPDRefit[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 544     fhPtDCASPDRefit[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 545     outputContainer->Add(fhPtDCASPDRefit[i]);
 546
 547     fhPtDCANoSPDRefit[i]  = new TH2F(Form("hPtDCA%sNoSPDRefit",dcaName[i].Data()),
 548                                  Form("Track DCA%s vs #it{p}_{T}distributionof constrained tracks no SPD and with ITSRefit",dcaName[i].Data()),
 549                                  nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 550     fhPtDCANoSPDRefit[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 551     fhPtDCANoSPDRefit[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 552     outputContainer->Add(fhPtDCANoSPDRefit[i]);
 553
 554     fhPtDCANoSPDNoRefit[i]  = new TH2F(Form("hPtDCA%sNoSPDNoRefit",dcaName[i].Data()),
 555                            Form("Track DCA%s vs #it{p}_{T}distribution, constrained tracks with no SPD requierement and without ITSRefit",dcaName[i].Data()),
 556                            nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 557     fhPtDCANoSPDNoRefit[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 558     fhPtDCANoSPDNoRefit[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 559     outputContainer->Add(fhPtDCANoSPDNoRefit[i]);
 560
 561     if(fFillTrackBCHistograms)
 562     {
 563       fhPtDCATOFBC0[i]  = new TH2F(Form("hPtDCA%sTOFBC0",dcaName[i].Data()),
 564                                    Form("Track DCA%s vs #it{p}_{T}distribution, BC=0",dcaName[i].Data()),
 565                                    nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 566       fhPtDCATOFBC0[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 567       fhPtDCATOFBC0[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 568       outputContainer->Add(fhPtDCATOFBC0[i]);
 569
 570       fhPtDCATOFBCOut[i]  = new TH2F(Form("hPtDCA%sTOFBCOut",dcaName[i].Data()),
 571                                      Form("Track DCA%s vs #it{p}_{T}distribution, BC!=0",dcaName[i].Data()),
 572                                      nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 573       fhPtDCATOFBCOut[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 574       fhPtDCATOFBCOut[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 575       outputContainer->Add(fhPtDCATOFBCOut[i]);
 576     }
 577
 578     fhPtDCANoTOFHit[i]  = new TH2F(Form("hPtDCA%sNoTOFHit",dcaName[i].Data()),
 579                            Form("Track (no TOF hit) DCA%s vs #it{p}_{T}distribution",dcaName[i].Data()),
 580                            nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 581     fhPtDCANoTOFHit[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 582     fhPtDCANoTOFHit[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 583     outputContainer->Add(fhPtDCANoTOFHit[i]);
 584
 585     if(fFillVertexBC0Histograms)
 586     {
 587       fhPtDCAVtxOutBC0[i]  = new TH2F(Form("hPtDCA%sVtxOutBC0",dcaName[i].Data()),
 588                                       Form("Track DCA%s vs #it{p}_{T}distribution, vertex with BC!=0",dcaName[i].Data()),
 589                                       nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 590       fhPtDCAVtxOutBC0[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 591       fhPtDCAVtxOutBC0[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 592       outputContainer->Add(fhPtDCAVtxOutBC0[i]);
 593
 594       fhPtDCAVtxOutBC0NoTOFHit[i]  = new TH2F(Form("hPtDCA%sVtxOutBC0NoTOFHit",dcaName[i].Data()),
 595                                               Form("Track (no TOF hit) DCA%s vs #it{p}_{T}distribution, vertex with BC!=0",dcaName[i].Data()),
 596                                               nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 597       fhPtDCAVtxOutBC0NoTOFHit[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 598       fhPtDCAVtxOutBC0NoTOFHit[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 599       outputContainer->Add(fhPtDCAVtxOutBC0NoTOFHit[i]);
 600
 601       fhPtDCAVtxInBC0[i]  = new TH2F(Form("hPtDCA%sVtxInBC0",dcaName[i].Data()),
 602                                       Form("Track DCA%s vs #it{p}_{T}distribution, vertex with BC==0",dcaName[i].Data()),
 603                                       nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 604       fhPtDCAVtxInBC0[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 605       fhPtDCAVtxInBC0[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 606       outputContainer->Add(fhPtDCAVtxInBC0[i]);
 607
 608       fhPtDCAVtxInBC0NoTOFHit[i]  = new TH2F(Form("hPtDCA%sVtxInBC0NoTOFHit",dcaName[i].Data()),
 609                                               Form("Track (no TOF hit) DCA%s vs #it{p}_{T}distribution, vertex with BC==0",dcaName[i].Data()),
 610                                               nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 611       fhPtDCAVtxInBC0NoTOFHit[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 612       fhPtDCAVtxInBC0NoTOFHit[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 613       outputContainer->Add(fhPtDCAVtxInBC0NoTOFHit[i]);
 614     }
 615
 616     if(IsPileUpAnalysisOn())
 617     {
 618       fhPtDCAPileUp[i]  = new TH2F(Form("hPtDCA%sPileUp",dcaName[i].Data()),
 619                              Form("Track DCA%s vs #it{p}_{T}distribution, SPD Pile-Up",dcaName[i].Data()),
 620                              nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 621       fhPtDCAPileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 622       fhPtDCAPileUp[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 623       outputContainer->Add(fhPtDCAPileUp[i]);
 624
 625       if(fFillTrackBCHistograms)
 626       {
 627         fhPtDCAPileUpTOFBC0[i]  = new TH2F(Form("hPtDCA%sPileUpTOFBC0",dcaName[i].Data()),
 628                                            Form("Track DCA%s vs #it{p}_{T}distribution",dcaName[i].Data()),
 629                                            nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 630         fhPtDCAPileUpTOFBC0[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 631         fhPtDCAPileUpTOFBC0[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 632         outputContainer->Add(fhPtDCAPileUpTOFBC0[i]);
 633       }
 634
 635       fhPtDCAPileUpNoTOFHit[i]  = new TH2F(Form("hPtDCA%sPileUpNoTOFHit",dcaName[i].Data()),
 636                                      Form("Track (no TOF hit) DCA%s vs #it{p}_{T}distribution, SPD Pile-Up, vertex with BC!=0",dcaName[i].Data()),
 637                                      nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 638       fhPtDCAPileUpNoTOFHit[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 639       fhPtDCAPileUpNoTOFHit[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 640       outputContainer->Add(fhPtDCAPileUpNoTOFHit[i]);
 641
 642       if(fFillVertexBC0Histograms)
 643       {
 644         fhPtDCAVtxOutBC0PileUp[i]  = new TH2F(Form("hPtDCA%sPileUpVtxOutBC0",dcaName[i].Data()),
 645                                               Form("Track DCA%s vs #it{p}_{T}distribution, SPD Pile-Up, vertex with BC!=0",dcaName[i].Data()),
 646                                               nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 647         fhPtDCAVtxOutBC0PileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 648         fhPtDCAVtxOutBC0PileUp[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 649         outputContainer->Add(fhPtDCAVtxOutBC0PileUp[i]);
 650
 651         fhPtDCAVtxOutBC0PileUpNoTOFHit[i]  = new TH2F(Form("hPtDCA%sVtxOutBC0PileUpNoTOFHit",dcaName[i].Data()),
 652                                                       Form("Track (no TOF hit) DCA%s vs #it{p}_{T}distribution, SPD Pile-Up, vertex with BC!=0",dcaName[i].Data()),
 653                                                       nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 654         fhPtDCAVtxOutBC0PileUpNoTOFHit[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 655         fhPtDCAVtxOutBC0PileUpNoTOFHit[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 656         outputContainer->Add(fhPtDCAVtxOutBC0PileUpNoTOFHit[i]);
 657
 658         fhPtDCAVtxInBC0PileUp[i]  = new TH2F(Form("hPtDCA%sPileUpVtxInBC0",dcaName[i].Data()),
 659                                               Form("Track DCA%s vs #it{p}_{T}distribution, SPD Pile-Up,vertex with BC==0",dcaName[i].Data()),
 660                                               nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 661         fhPtDCAVtxInBC0PileUp[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 662         fhPtDCAVtxInBC0PileUp[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 663         outputContainer->Add(fhPtDCAVtxInBC0PileUp[i]);
 664
 665         fhPtDCAVtxInBC0PileUpNoTOFHit[i]  = new TH2F(Form("hPtDCA%sVtxInBC0PileUpNoTOFHit",dcaName[i].Data()),
 666                                                       Form("Track (no TOF hit) DCA%s vs #it{p}_{T}distribution, SPD Pile-Up, vertex with BC==0",dcaName[i].Data()),
 667                                                       nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 668         fhPtDCAVtxInBC0PileUpNoTOFHit[i]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 669         fhPtDCAVtxInBC0PileUpNoTOFHit[i]->SetYTitle(Form("DCA_{%s}",dcaName[i].Data()));
 670         outputContainer->Add(fhPtDCAVtxInBC0PileUpNoTOFHit[i]);
 671
 672       }
 673     }
 674   }
 675
 676   if(IsDataMC())
 677   {
 678     //enum mvType{kmcPion = 0, kmcProton = 1, kmcKaon = 2, kmcMuon = 3, kmcElectron = 4, kmcUnknown = 4 };
 679
 680     TString histoName[] = {"Pion","Proton","Kaon","Muon","Electron","Unknown"};
 681     TString titleName[] = {"#pi^{#pm}","p^{#pm}","K^{#pm}","#mu^{#pm}","e^{#pm}","x^{#pm}"};
 682
 683     for(Int_t imcPart = 0; imcPart < 6; imcPart++)
 684     {
 685       fhPtMCPart[imcPart]  = new TH1F (Form("hPtMC%s",histoName[imcPart].Data()),
 686                                        Form("reconstructed #it{p}_{T} distribution from %s",titleName[imcPart].Data()),
 687                                        nptbins,ptmin,ptmax);
 688       fhPtMCPart[imcPart]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 689       outputContainer->Add(fhPtMCPart[imcPart]);
 690
 691       fhPhiMCPart[imcPart]  = new TH2F (Form("hPhiMC%s",histoName[imcPart].Data()),
 692                                         Form("reconstructed #phi vs #it{p}_{T} distribution from %s",titleName[imcPart].Data()),
 693                                         nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 694       fhPhiMCPart[imcPart]->SetYTitle("#phi (rad)");
 695       fhPhiMCPart[imcPart]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 696       outputContainer->Add(fhPhiMCPart[imcPart]);
 697
 698       fhEtaMCPart[imcPart]  = new TH2F (Form("hEtaMC%s",histoName[imcPart].Data()),
 699                                         Form("reconstructed #eta vs #it{p}_{T} distribution from %s",titleName[imcPart].Data()),
 700                                         nptbins,ptmin,ptmax, netabins,etamin,etamax);
 701       fhEtaMCPart[imcPart]->SetYTitle("#eta ");
 702       fhEtaMCPart[imcPart]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 703       outputContainer->Add(fhEtaMCPart[imcPart]);
 704
 705       fhPtMCPrimPart[imcPart]  = new TH1F (Form("hPtMCPrimary%s",histoName[imcPart].Data()),
 706                                            Form("generated #it{p}_{T} distribution from %s",titleName[imcPart].Data()),
 707                                            nptbins,ptmin,ptmax);
 708       fhPtMCPrimPart[imcPart]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 709       outputContainer->Add(fhPtMCPrimPart[imcPart]);
 710
 711       fhPhiMCPrimPart[imcPart]  = new TH2F (Form("hPhiMCPrimary%s",histoName[imcPart].Data()),
 712                                             Form("generated #phi vs #it{p}_{T} distribution from %s",titleName[imcPart].Data()),
 713                                             nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 714       fhPhiMCPrimPart[imcPart]->SetYTitle("#phi (rad)");
 715       fhPhiMCPrimPart[imcPart]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 716       outputContainer->Add(fhPhiMCPrimPart[imcPart]);
 717
 718       fhEtaMCPrimPart[imcPart]  = new TH2F (Form("hEtaMCPrimary%s",histoName[imcPart].Data()),
 719                                             Form("generated #eta vs #it{p}_{T} distribution from %s",titleName[imcPart].Data()),
 720                                             nptbins,ptmin,ptmax, netabins,etamin,etamax);
 721       fhEtaMCPrimPart[imcPart]->SetYTitle("#eta ");
 722       fhEtaMCPrimPart[imcPart]->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 723       outputContainer->Add(fhEtaMCPrimPart[imcPart]);
 724     }
 725
 726   }
 727
 728   fhPtNPileUpSPDVtx  = new TH2F ("hPt_NPileUpVertSPD","pT of cluster vs N pile-up SPD vertex",
 729                                  nptbins,ptmin,ptmax,20,0,20);
 730   fhPtNPileUpSPDVtx->SetYTitle("# vertex ");
 731   fhPtNPileUpSPDVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 732   outputContainer->Add(fhPtNPileUpSPDVtx);
 733
 734   fhPtNPileUpTrkVtx  = new TH2F ("hPt_NPileUpVertTracks","pT of cluster vs N pile-up Tracks vertex",
 735                                  nptbins,ptmin,ptmax, 20,0,20 );
 736   fhPtNPileUpTrkVtx->SetYTitle("# vertex ");
 737   fhPtNPileUpTrkVtx->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 738   outputContainer->Add(fhPtNPileUpTrkVtx);
 739
 740   if(fFillVertexBC0Histograms)
 741   {
 742     fhPtNPileUpSPDVtxBC0  = new TH2F ("hPt_NPileUpVertSPD_BC0","pT of cluster vs N pile-up SPD vertex",
 743                                    nptbins,ptmin,ptmax,20,0,20);
 744     fhPtNPileUpSPDVtxBC0->SetYTitle("# vertex ");
 745     fhPtNPileUpSPDVtxBC0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 746     outputContainer->Add(fhPtNPileUpSPDVtxBC0);
 747
 748     fhPtNPileUpTrkVtxBC0  = new TH2F ("hPt_NPileUpVertTracks_BC0","pT of cluster vs N pile-up Tracks vertex",
 749                                    nptbins,ptmin,ptmax, 20,0,20 );
 750     fhPtNPileUpTrkVtxBC0->SetYTitle("# vertex ");
 751     fhPtNPileUpTrkVtxBC0->SetXTitle("#it{p}_{T} (GeV/#it{c})");
 752     outputContainer->Add(fhPtNPileUpTrkVtxBC0);
 753   }
 754
 755   return outputContainer;
 756
 757 }
 758
 759
 760 //___________________________________________
 761 void AliAnaChargedParticles::InitParameters()
 762 {
 763   //Initialize the parameters of the analysis.
 764   SetOutputAODClassName("AliAODPWG4Particle");
 765   SetOutputAODName("PWG4Particle");
 766
 767   AddToHistogramsName("AnaCharged_");
 768
 769 }
 770
 771 //____________________________________________________________
 772 void AliAnaChargedParticles::Print(const Option_t * opt) const
 773 {
 774   //Print some relevant parameters set for the analysis
 775   if(! opt)
 776     return;
 777
 778   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 779   AliAnaCaloTrackCorrBaseClass::Print(" ");
 780
 781   printf("Min Pt = %3.2f\n", GetMinPt());
 782   printf("Max Pt = %3.2f\n", GetMaxPt());
 783
 784 }
 785
 786 //_________________________________
 787 void AliAnaChargedParticles::Init()
 788 {
 789   //Init
 790   //Do some checks
 791
 792   if(!GetReader()->IsCTSSwitchedOn())
 793     AliFatal("STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!\n");
 794
 795
 796 }
 797
 798 //_________________________________________________
 799 void  AliAnaChargedParticles::MakeAnalysisFillAOD()
 800 {
 801   //Do analysis and fill aods
 802   if(!GetCTSTracks() || GetCTSTracks()->GetEntriesFast() == 0) return ;
 803
 804   Int_t ntracks = GetCTSTracks()->GetEntriesFast();
 805   Double_t vert[3] = {0,0,0}; //vertex ;
 806
 807   //Some prints
 808   if(GetDebug() > 0)
 809     printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);
 810
 811   AliVEvent  * event = GetReader()->GetInputEvent();
 812   AliESDEvent* esdEv = dynamic_cast<AliESDEvent*> (event);
 813   AliAODEvent* aodEv = dynamic_cast<AliAODEvent*> (event);
 814
 815   Int_t vtxBC = GetReader()->GetVertexBC();
 816   if(!GetReader()->IsDCACutOn()) vtxBC = GetReader()->GetVertexBC(event->GetPrimaryVertex());
 817
 818   if(fFillVertexBC0Histograms)
 819   {
 820     fhProductionVertexBC->Fill(vtxBC);
 821     if(IsPileUpAnalysisOn())
 822     {
 823       if(GetReader()->IsPileUpFromSPD())               fhProductionVertexBCPileUp[0]->Fill(vtxBC);
 824       if(GetReader()->IsPileUpFromEMCal())             fhProductionVertexBCPileUp[1]->Fill(vtxBC);
 825       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhProductionVertexBCPileUp[2]->Fill(vtxBC);
 826       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhProductionVertexBCPileUp[3]->Fill(vtxBC);
 827       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhProductionVertexBCPileUp[4]->Fill(vtxBC);
 828       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhProductionVertexBCPileUp[5]->Fill(vtxBC);
 829       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhProductionVertexBCPileUp[6]->Fill(vtxBC);
 830     }
 831   }
 832
 833   // N pile up vertices
 834   Int_t nVtxSPD = -1;
 835   Int_t nVtxTrk = -1;
 836
 837   if      (esdEv)
 838   {
 839     nVtxSPD = esdEv->GetNumberOfPileupVerticesSPD();
 840     nVtxTrk = esdEv->GetNumberOfPileupVerticesTracks();
 841
 842   }//ESD
 843   else if (aodEv)
 844   {
 845     nVtxSPD = aodEv->GetNumberOfPileupVerticesSPD();
 846     nVtxTrk = aodEv->GetNumberOfPileupVerticesTracks();
 847   }//AOD
 848
 849
 850   //printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - primary vertex BC %d\n",vtxBC);
 851
 852   Double_t bz = event->GetMagneticField();
 853
 854   //Fill AODParticle with CTS aods
 855   Float_t pt  = 0;
 856   Float_t phi = 0;
 857   Float_t eta = 0;
 858   Int_t evtIndex = 0;
 859   for(Int_t i = 0; i < ntracks; i++)
 860   {
 861     AliVTrack * track =  (AliVTrack*) (GetCTSTracks()->At(i));
 862
 863     pt  = track->Pt();
 864     eta = track->Eta();
 865     phi = track->Phi();
 866
 867     fhPtNoCut->Fill(pt);
 868
 869     fhPtNPileUpSPDVtx->Fill(pt,nVtxSPD);
 870     fhPtNPileUpTrkVtx->Fill(pt,nVtxTrk);
 871
 872     AliAODTrack * aodTrack = dynamic_cast<AliAODTrack*>(track);
 873     AliESDtrack * esdTrack = dynamic_cast<AliESDtrack*>(track);
 874
 875     //TOF
 876     ULong_t status = track->GetStatus();
 877     Bool_t okTOF = (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ;
 878     Double32_t tof = track->GetTOFsignal()*1e-3;
 879
 880     //DCA
 881     Double_t dcaCons = -999;
 882     if(aodTrack)
 883     {
 884       dcaCons = aodTrack->DCA();
 885       //vtxBC   = aodTrack->GetProdVertex()->GetBC();
 886     }
 887
 888     Double_t dca[2]   = {1e6,1e6};
 889     Double_t covar[3] = {1e6,1e6,1e6};
 890     track->PropagateToDCA(GetReader()->GetInputEvent()->GetPrimaryVertex(),bz,100.,dca,covar);
 891
 892     Float_t trackDCA = dca[0];
 893
 894     if(dcaCons == -999)
 895     {
 896       fhPtDCA[0]->Fill(pt, dca[0]);
 897       fhPtDCA[1]->Fill(pt, dca[1]);
 898     }
 899     else
 900     {
 901       trackDCA = dcaCons;
 902       fhPtDCA[2]->Fill(pt, dcaCons);
 903     }
 904
 905     if(GetReader()->AcceptDCA(pt,trackDCA) && !GetReader()->IsDCACutOn() ) fhPtCutDCA->Fill(pt);
 906
 907     if(fFillVertexBC0Histograms)
 908     {
 909       if(TMath::Abs(vtxBC) > 0 && vtxBC!=AliVTrack::kTOFBCNA)
 910       {
 911         fhPtVtxOutBC0->Fill(pt);
 912         fhEtaPhiVtxOutBC0->Fill(eta,phi);
 913
 914         if(dcaCons == -999)
 915         {
 916           fhPtDCAVtxOutBC0[0]->Fill(pt, dca[0]);
 917           fhPtDCAVtxOutBC0[1]->Fill(pt, dca[1]);
 918         }
 919         else
 920           fhPtDCAVtxOutBC0[2]->Fill(pt, dcaCons);
 921       }
 922       else
 923       {
 924         fhPtVtxInBC0->Fill(pt);
 925         fhEtaPhiVtxInBC0->Fill(eta,phi);
 926
 927         fhPtNPileUpSPDVtxBC0->Fill(pt,nVtxSPD);
 928         fhPtNPileUpTrkVtxBC0->Fill(pt,nVtxTrk);
 929
 930         if(GetReader()->AcceptDCA(pt,trackDCA) && !GetReader()->IsDCACutOn() ) fhPtCutDCABCOK->Fill(pt);
 931
 932         if(dcaCons == -999)
 933         {
 934           fhPtDCAVtxInBC0[0]->Fill(pt, dca[0]);
 935           fhPtDCAVtxInBC0[1]->Fill(pt, dca[1]);
 936         }
 937         else
 938           fhPtDCAVtxInBC0[2]->Fill(pt, dcaCons);
 939
 940       }
 941     }
 942
 943     if(IsPileUpAnalysisOn() && GetReader()->IsPileUpFromSPD())
 944     {
 945       if(dcaCons == -999)
 946       {
 947         fhPtDCAPileUp[0]->Fill(pt, dca[0]);
 948         fhPtDCAPileUp[1]->Fill(pt, dca[1]);
 949       }
 950       else
 951         fhPtDCAPileUp[2]->Fill(pt, dcaCons);
 952
 953       if(fFillVertexBC0Histograms)
 954       {
 955         if(TMath::Abs(vtxBC) > 0 && vtxBC!=AliVTrack::kTOFBCNA)
 956         {
 957           if(dcaCons == -999)
 958           {
 959             fhPtDCAVtxOutBC0PileUp[0]->Fill(pt, dca[0]);
 960             fhPtDCAVtxOutBC0PileUp[1]->Fill(pt, dca[1]);
 961           }
 962           else fhPtDCAVtxOutBC0PileUp[2]->Fill(pt, dcaCons);
 963         }
 964         else
 965         {
 966           if(dcaCons == -999)
 967           {
 968             fhPtDCAVtxInBC0PileUp[0]->Fill(pt, dca[0]);
 969             fhPtDCAVtxInBC0PileUp[1]->Fill(pt, dca[1]);
 970           }
 971           else fhPtDCAVtxInBC0PileUp[2]->Fill(pt, dcaCons);
 972         }
 973       }
 974     }
 975
 976     if(!okTOF)
 977     {
 978       if(dcaCons == -999)
 979       {
 980         fhPtDCANoTOFHit[0]->Fill(pt, dca[0]);
 981         fhPtDCANoTOFHit[1]->Fill(pt, dca[1]);
 982       }
 983       else
 984         fhPtDCANoTOFHit[2]->Fill(pt, dcaCons);
 985
 986       if(fFillVertexBC0Histograms)
 987       {
 988         if(TMath::Abs(vtxBC) > 0 && vtxBC!=AliVTrack::kTOFBCNA)
 989         {
 990           if(dcaCons == -999)
 991           {
 992             fhPtDCAVtxOutBC0NoTOFHit[0]->Fill(pt, dca[0]);
 993             fhPtDCAVtxOutBC0NoTOFHit[1]->Fill(pt, dca[1]);
 994           }
 995           else
 996             fhPtDCAVtxOutBC0NoTOFHit[2]->Fill(pt, dcaCons);
 997         }
 998         else
 999         {
1000           if(dcaCons == -999)
1001           {
1002             fhPtDCAVtxInBC0NoTOFHit[0]->Fill(pt, dca[0]);
1003             fhPtDCAVtxInBC0NoTOFHit[1]->Fill(pt, dca[1]);
1004           }
1005           else
1006             fhPtDCAVtxInBC0NoTOFHit[2]->Fill(pt, dcaCons);
1007
1008         }
1009       }
1010
1011       if(IsPileUpAnalysisOn() && GetReader()->IsPileUpFromSPD())
1012       {
1013         if(dcaCons == -999)
1014         {
1015           fhPtDCAPileUpNoTOFHit[0]->Fill(pt, dca[0]);
1016           fhPtDCAPileUpNoTOFHit[1]->Fill(pt, dca[1]);
1017         }
1018         else
1019           fhPtDCAPileUpNoTOFHit[2]->Fill(pt, dcaCons);
1020
1021         if(fFillVertexBC0Histograms)
1022         {
1023           if(TMath::Abs(vtxBC) > 0 && vtxBC!=AliVTrack::kTOFBCNA)
1024           {
1025             if(dcaCons == -999)
1026             {
1027               fhPtDCAVtxOutBC0PileUpNoTOFHit[0]->Fill(pt, dca[0]);
1028               fhPtDCAVtxOutBC0PileUpNoTOFHit[1]->Fill(pt, dca[1]);
1029             }
1030             else
1031               fhPtDCAVtxOutBC0PileUpNoTOFHit[2]->Fill(pt, dcaCons);
1032           }
1033           else
1034           {
1035             if(dcaCons == -999)
1036             {
1037               fhPtDCAVtxInBC0PileUpNoTOFHit[0]->Fill(pt, dca[0]);
1038               fhPtDCAVtxInBC0PileUpNoTOFHit[1]->Fill(pt, dca[1]);
1039             }
1040             else
1041               fhPtDCAVtxInBC0PileUpNoTOFHit[2]->Fill(pt, dcaCons);
1042
1043           }
1044         }
1045       }
1046     }
1047
1048     //printf("track pT %2.2f, DCA Cons %f, DCA1 %f, DCA2 %f, TOFBC %d, oktof %d, tof %f\n",
1049     //      pt,dcaCons,dca[0],dca[1],track->GetTOFBunchCrossing(bz),okTOF, tof);
1050
1051
1052 //    if( vtxBC == 0 && trackBC !=0 && trackBC!=AliVTrack::kTOFBCNA)
1053 //      printf("TOF Signal %e, BC %d, pt %f, dca_xy %f, dca_z %f, dca_tpc %f \n", tof,trackBC, pt,dca[0],dca[1],dcaCons);
1054
1055
1056     if(okTOF)
1057     {
1058       fhTOFSignal  ->Fill(tof);
1059       fhPtTOFSignal->Fill(pt, tof);
1060       if(GetReader()->AcceptDCA(pt,trackDCA) && !GetReader()->IsDCACutOn() ) fhPtTOFSignalDCACut->Fill(pt, tof);
1061
1062       if(fFillVertexBC0Histograms)
1063       {
1064         if(TMath::Abs(vtxBC) > 0 && vtxBC!=AliVTrack::kTOFBCNA)
1065           fhPtTOFSignalVtxOutBC0->Fill(pt, tof);
1066         else
1067           fhPtTOFSignalVtxInBC0->Fill(pt, tof);
1068       }
1069
1070       Int_t trackBC = 1000;
1071
1072       if(fFillTrackBCHistograms)
1073       {
1074         trackBC = track->GetTOFBunchCrossing(bz);
1075
1076         if(trackBC==0)
1077         {
1078           fhTOFSignalBCOK->Fill(tof);
1079
1080           if(dcaCons == -999)
1081           {
1082             fhPtDCATOFBC0[0]->Fill(pt, dca[0]);
1083             fhPtDCATOFBC0[1]->Fill(pt, dca[1]);
1084           }
1085           else
1086             fhPtDCATOFBC0[2]->Fill(pt, dcaCons);
1087
1088           if(IsPileUpAnalysisOn() && GetReader()->IsPileUpFromSPD())
1089           {
1090             if(dcaCons == -999)
1091             {
1092               fhPtDCAPileUpTOFBC0[0]->Fill(pt, dca[0]);
1093               fhPtDCAPileUpTOFBC0[1]->Fill(pt, dca[1]);
1094             }
1095             else
1096               fhPtDCAPileUpTOFBC0[2]->Fill(pt, dcaCons);
1097           }
1098         }
1099         else if(trackBC!=AliVTrack::kTOFBCNA)
1100         {
1101           if(dcaCons == -999)
1102           {
1103             fhPtDCATOFBCOut[0]->Fill(pt, dca[0]);
1104             fhPtDCATOFBCOut[1]->Fill(pt, dca[1]);
1105           }
1106           else
1107             fhPtDCATOFBCOut[2]->Fill(pt, dcaCons);
1108
1109         }
1110       }
1111
1112       if(IsPileUpAnalysisOn())
1113       {
1114         if(GetReader()->IsPileUpFromSPD())               fhPtTOFSignalPileUp[0]->Fill(pt, tof);
1115         if(GetReader()->IsPileUpFromEMCal())             fhPtTOFSignalPileUp[1]->Fill(pt, tof);
1116         if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtTOFSignalPileUp[2]->Fill(pt, tof);
1117         if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtTOFSignalPileUp[3]->Fill(pt, tof);
1118         if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtTOFSignalPileUp[4]->Fill(pt, tof);
1119         if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtTOFSignalPileUp[5]->Fill(pt, tof);
1120         if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtTOFSignalPileUp[6]->Fill(pt, tof);
1121
1122         if(fFillTrackBCHistograms)
1123         {
1124           if      (trackBC ==0)  { fhEtaPhiTOFBC0    ->Fill(eta,phi); if(IsPileUpAnalysisOn() && GetReader()->IsPileUpFromSPD()) fhEtaPhiTOFBC0PileUpSPD    ->Fill(eta,phi); }
1125           else if (trackBC < 0)  { fhEtaPhiTOFBCPlus ->Fill(eta,phi); if(IsPileUpAnalysisOn() && GetReader()->IsPileUpFromSPD()) fhEtaPhiTOFBCPlusPileUpSPD ->Fill(eta,phi); }
1126           else if (trackBC > 0)  { fhEtaPhiTOFBCMinus->Fill(eta,phi); if(IsPileUpAnalysisOn() && GetReader()->IsPileUpFromSPD()) fhEtaPhiTOFBCMinusPileUpSPD->Fill(eta,phi); }
1127         }
1128
1129         if(fFillVertexBC0Histograms)
1130         {
1131           if(TMath::Abs(vtxBC) > 0 && vtxBC!=AliVTrack::kTOFBCNA)
1132           {
1133             if(GetReader()->IsPileUpFromSPD())               fhPtTOFSignalVtxOutBC0PileUp[0]->Fill(pt, tof);
1134             if(GetReader()->IsPileUpFromEMCal())             fhPtTOFSignalVtxOutBC0PileUp[1]->Fill(pt, tof);
1135             if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtTOFSignalVtxOutBC0PileUp[2]->Fill(pt, tof);
1136             if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtTOFSignalVtxOutBC0PileUp[3]->Fill(pt, tof);
1137             if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtTOFSignalVtxOutBC0PileUp[4]->Fill(pt, tof);
1138             if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtTOFSignalVtxOutBC0PileUp[5]->Fill(pt, tof);
1139             if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtTOFSignalVtxOutBC0PileUp[6]->Fill(pt, tof);
1140           }
1141           else
1142           {
1143             if(GetReader()->IsPileUpFromSPD())               fhPtTOFSignalVtxInBC0PileUp[0]->Fill(pt, tof);
1144             if(GetReader()->IsPileUpFromEMCal())             fhPtTOFSignalVtxInBC0PileUp[1]->Fill(pt, tof);
1145             if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtTOFSignalVtxInBC0PileUp[2]->Fill(pt, tof);
1146             if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtTOFSignalVtxInBC0PileUp[3]->Fill(pt, tof);
1147             if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtTOFSignalVtxInBC0PileUp[4]->Fill(pt, tof);
1148             if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtTOFSignalVtxInBC0PileUp[5]->Fill(pt, tof);
1149             if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtTOFSignalVtxInBC0PileUp[6]->Fill(pt, tof);
1150           }
1151         }
1152       }
1153     }
1154
1155     //Fill AODParticle after some selection
1156     fMomentum.SetPxPyPzE(track->Px(),track->Py(),track->Pz(),0);
1157     Bool_t in = GetFiducialCut()->IsInFiducialCut(fMomentum.Eta(),fMomentum.Phi(),kCTS) ;
1158
1159     if(GetDebug() > 1)
1160       printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - Track pt %2.2f, eta %2.2f, phi %2.2f in fiducial cut %d\n",pt,eta,phi,in);
1161
1162     //Acceptance selection
1163     if(IsFiducialCutOn() && ! in ) continue ;
1164
1165     // Momentum selection
1166     if(pt < GetMinPt() || pt > GetMaxPt()) continue;
1167
1168     if(okTOF) fhTOFSignalPtCut->Fill(tof);
1169     else
1170     {
1171       fhPtTOFStatus0    ->Fill(pt);
1172       fhEtaPhiTOFStatus0->Fill(eta,phi);
1173     }
1174
1175     Bool_t bITSRefit    = (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit;
1176     Bool_t bConstrained = kFALSE;
1177     if     (aodTrack) bConstrained = aodTrack->IsGlobalConstrained();
1178     else if(esdTrack) bConstrained = (!esdTrack->HasPointOnITSLayer(0) && !esdTrack->HasPointOnITSLayer(1));
1179     //printf("Track %d, pt %2.2f, eta %2.2f, phi %2.2f, SPDRefit %d, refit %d, dcaCons %2.2f\n",
1180     //       i, pt, eta, phi, bConstrained, bITSRefit, dcaCons);
1181
1182     if(bConstrained)
1183     {
1184       if(bITSRefit)
1185       {
1186         fhPtNoSPDRefit->Fill(pt);
1187         if(pt < 2)fhEtaPhiNoSPDRefitPt02->Fill(eta,phi);
1188         if(pt > 3)fhEtaPhiNoSPDRefitPt3 ->Fill(eta,phi);
1189
1190         if(dcaCons == -999)
1191         {
1192           fhPtDCANoSPDRefit[0]->Fill(pt, dca[0]);
1193           fhPtDCANoSPDRefit[1]->Fill(pt, dca[1]);
1194         }
1195         else
1196           fhPtDCANoSPDRefit[2]->Fill(pt, dcaCons);
1197
1198       }
1199       else
1200       {
1201         fhPtNoSPDNoRefit->Fill(pt);
1202         if(pt < 2)fhEtaPhiNoSPDNoRefitPt02->Fill(eta,phi);
1203         if(pt > 3)fhEtaPhiNoSPDNoRefitPt3 ->Fill(eta,phi);
1204         if(dcaCons == -999)
1205         {
1206           fhPtDCANoSPDNoRefit[0]->Fill(pt, dca[0]);
1207           fhPtDCANoSPDNoRefit[1]->Fill(pt, dca[1]);
1208         }
1209         else
1210           fhPtDCANoSPDNoRefit[2]->Fill(pt, dcaCons);
1211
1212       }
1213     }
1214     else
1215     {
1216       fhPtSPDRefit->Fill(pt);
1217       if(pt < 2)fhEtaPhiSPDRefitPt02->Fill(eta,phi);
1218       if(pt > 3)fhEtaPhiSPDRefitPt3 ->Fill(eta,phi);
1219       if(dcaCons == -999)
1220       {
1221         fhPtDCASPDRefit[0]->Fill(pt, dca[0]);
1222         fhPtDCASPDRefit[1]->Fill(pt, dca[1]);
1223       }
1224       else
1225         fhPtDCASPDRefit[2]->Fill(pt, dcaCons);
1226     }
1227
1228     // Mixed event
1229     if (GetMixedEvent())
1230     {
1231       evtIndex = GetMixedEvent()->EventIndex(track->GetID()) ;
1232     }
1233
1234     GetVertex(vert,evtIndex);
1235     if(TMath::Abs(vert[2])> GetZvertexCut()) return;
1236
1237     AliAODPWG4Particle tr = AliAODPWG4Particle(track->Px(),track->Py(),track->Pz(),0);
1238     tr.SetDetectorTag(kCTS);
1239     tr.SetLabel(track->GetLabel());
1240     tr.SetTrackLabel(track->GetID(),-1);
1241     tr.SetChargedBit(track->Charge()>0);
1242
1243     AddAODParticle(tr);
1244
1245   }//loop
1246
1247   if(GetDebug() > 0)
1248     printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
1249 }
1250
1251 //__________________________________________________________________
1252 void  AliAnaChargedParticles::MakeAnalysisFillHistograms()
1253 {
1254   //Do analysis and fill histograms
1255
1256   if(IsDataMC()) FillPrimaryHistograms();
1257
1258   //Loop on stored AODParticles
1259   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
1260
1261   fhNtracks->Fill(GetReader()->GetTrackMultiplicity()) ;
1262
1263   if(GetDebug() > 0)
1264     printf("AliAnaChargedParticles::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
1265
1266   Float_t pt  = 0;
1267   Float_t phi = 0;
1268   Float_t eta = 0;
1269
1270   for(Int_t iaod = 0; iaod < naod ; iaod++)
1271   {
1272     AliAODPWG4Particle* track =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
1273
1274     pt  = track->Pt();
1275     eta = track->Eta();
1276     phi = track->Phi();
1277
1278     fhPt->Fill(pt);
1279
1280     if(track->GetChargedBit())
1281     {
1282       fhPhiPos   ->Fill(pt, phi);
1283       fhEtaPos   ->Fill(pt, eta);
1284       fhEtaPhiPos->Fill(eta,phi);
1285     }
1286     else
1287     {
1288       fhPhiNeg   ->Fill(pt, phi);
1289       fhEtaNeg   ->Fill(pt, eta);
1290       fhEtaPhiNeg->Fill(eta,phi);
1291     }
1292
1293     if(IsPileUpAnalysisOn())
1294     {
1295       if(GetReader()->IsPileUpFromSPD())               {fhPtPileUp[0]->Fill(pt);}
1296       if(GetReader()->IsPileUpFromEMCal())             {fhPtPileUp[1]->Fill(pt);}
1297       if(GetReader()->IsPileUpFromSPDOrEMCal())        {fhPtPileUp[2]->Fill(pt);}
1298       if(GetReader()->IsPileUpFromSPDAndEMCal())       {fhPtPileUp[3]->Fill(pt);}
1299       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    {fhPtPileUp[4]->Fill(pt);}
1300       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    {fhPtPileUp[5]->Fill(pt);}
1301       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtPileUp[6]->Fill(pt);}
1302     }
1303
1304
1305     if(IsDataMC())
1306     {
1307       //Play with the MC stack if available
1308       Int_t mompdg = -1;
1309       Int_t label  = track->GetLabel();
1310
1311       if(label >= 0)
1312       {
1313         if( GetReader()->ReadStack() && label < GetMCStack()->GetNtrack())
1314         {
1315           TParticle * mom = GetMCStack()->Particle(label);
1316           mompdg =TMath::Abs(mom->GetPdgCode());
1317         }
1318         else if(GetReader()->ReadAODMCParticles())
1319         {
1320           AliAODMCParticle * aodmom = 0;
1321           //Get the list of MC particles
1322           aodmom = (AliAODMCParticle*) (GetReader()->GetAODMCParticles())->At(label);
1323           mompdg =TMath::Abs(aodmom->GetPdgCode());
1324         }
1325       }
1326
1327       Int_t mcType = kmcUnknown;
1328       if     (mompdg==211 ) mcType = kmcPion;
1329       else if(mompdg==2212) mcType = kmcProton;
1330       else if(mompdg==321 ) mcType = kmcKaon;
1331       else if(mompdg==11  ) mcType = kmcElectron;
1332       else if(mompdg==13  ) mcType = kmcMuon;
1333
1334       fhPtMCPart [mcType]->Fill(pt);
1335       fhEtaMCPart[mcType]->Fill(pt,eta);
1336       fhPhiMCPart[mcType]->Fill(pt,phi);
1337
1338     }//Work with stack also
1339
1340   }// aod branch loop
1341
1342 }