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