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