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