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     fPdg(0),           fFillPileUpHistograms(0),
  48     fhNtracks(0),      fhPt(0),
  49     fhPhiNeg(0),       fhEtaNeg(0),
  50     fhPhiPos(0),       fhEtaPos(0),
  51     fhEtaPhiPos(0),    fhEtaPhiNeg(0),
  52     //MC
  53     fhPtPion(0),       fhPhiPion(0),         fhEtaPion(0),
  54     fhPtProton(0),     fhPhiProton(0),       fhEtaProton(0),
  55     fhPtElectron(0),   fhPhiElectron(0),     fhEtaElectron(0),
  56     fhPtKaon(0),       fhPhiKaon(0),         fhEtaKaon(0),
  57     fhPtUnknown(0),    fhPhiUnknown(0),      fhEtaUnknown(0),
  58     fhTOFSignal(0),    fhTOFSignalPtCut(0),  fhTOFSignalBCOK(0),
  59     fhPtTOFSignal(0),  fhPtTOFStatus0(0),    fhEtaPhiTOFStatus0(0),
  60     fhEtaPhiTOFBC0(0), fhEtaPhiTOFBCPlus(0), fhEtaPhiTOFBCMinus(0),
  61     fhEtaPhiTOFBC0PileUpSPD(0),
  62     fhEtaPhiTOFBCPlusPileUpSPD(0),
  63     fhEtaPhiTOFBCMinusPileUpSPD(0)
  64 {
  65   //Default Ctor
  66
  67   for(Int_t i = 0; i < 7; i++)
  68   {
  69     fhPtPileUp         [i] = 0;
  70     fhPtTOFSignalPileUp[i] = 0;
  71   }
  72
  73   for(Int_t i = 0; i < 3; i++)
  74   {
  75     fhPtDCA               [i] = 0 ;
  76     //fhPtDCAVtxOutBC0      [i] = 0 ;
  77     fhPtDCAPileUp         [i] = 0 ;
  78     //fhPtDCAVtxOutBC0PileUp[i] = 0 ;
  79
  80     fhPtDCATOFBC0         [i] = 0 ;
  81     fhPtDCAPileUpTOFBC0   [i] = 0 ;
  82
  83     fhPtDCANoTOFHit               [i] = 0 ;
  84     //fhPtDCAVtxOutBC0NoTOFHit      [i] = 0 ;
  85     fhPtDCAPileUpNoTOFHit         [i] = 0 ;
  86     //fhPtDCAVtxOutBC0PileUpNoTOFHit[i] = 0 ;
  87   }
  88
  89   //Initialize parameters
  90   InitParameters();
  91
  92 }
  93
  94 //_______________________________________________________
  95 TList *  AliAnaChargedParticles::GetCreateOutputObjects()
  96 {
  97   // Create histograms to be saved in output file and
  98   // store them in fOutputContainer
  99
 100
 101   TList * outputContainer = new TList() ;
 102   outputContainer->SetName("ExampleHistos") ;
 103
 104   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins(); Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
 105   Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();  Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
 106   Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();  Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();  Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
 107
 108   fhNtracks  = new TH1F ("hNtracks","# of tracks", 1000,0,1000);
 109   fhNtracks->SetXTitle("# of tracks");
 110   outputContainer->Add(fhNtracks);
 111
 112   fhPt  = new TH1F ("hPt","p_T distribution", nptbins,ptmin,ptmax);
 113   fhPt->SetXTitle("p_{T} (GeV/c)");
 114   outputContainer->Add(fhPt);
 115
 116   fhPhiNeg  = new TH2F ("hPhiNegative","#phi of negative charges distribution",
 117                         nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 118   fhPhiNeg->SetYTitle("#phi (rad)");
 119   fhPhiNeg->SetXTitle("p_{T} (GeV/c)");
 120   outputContainer->Add(fhPhiNeg);
 121
 122   fhEtaNeg  = new TH2F ("hEtaNegative","#eta of negative charges distribution",
 123                         nptbins,ptmin,ptmax, netabins,etamin,etamax);
 124   fhEtaNeg->SetYTitle("#eta ");
 125   fhEtaNeg->SetXTitle("p_{T} (GeV/c)");
 126   outputContainer->Add(fhEtaNeg);
 127
 128   fhPhiPos  = new TH2F ("hPhiPositive","#phi of positive charges distribution",
 129                         nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 130   fhPhiPos->SetYTitle("#phi (rad)");
 131   fhPhiPos->SetXTitle("p_{T} (GeV/c)");
 132   outputContainer->Add(fhPhiPos);
 133
 134   fhEtaPos  = new TH2F ("hEtaPositive","#eta of positive charges distribution",
 135                         nptbins,ptmin,ptmax, netabins,etamin,etamax);
 136   fhEtaPos->SetYTitle("#eta ");
 137   fhEtaPos->SetXTitle("p_{T} (GeV/c)");
 138   outputContainer->Add(fhEtaPos);
 139
 140   fhEtaPhiPos  = new TH2F ("hEtaPhiPositive","pt/eta/phi of positive charge",netabins,etamin,etamax, nphibins,phimin,phimax);
 141   fhEtaPhiPos->SetXTitle("#eta ");
 142   fhEtaPhiPos->SetYTitle("#phi (rad)");
 143   outputContainer->Add(fhEtaPhiPos);
 144
 145   fhEtaPhiNeg  = new TH2F ("hEtaPhiNegative","eta vs phi of negative charge",netabins,etamin,etamax, nphibins,phimin,phimax);
 146   fhEtaPhiNeg->SetXTitle("#eta ");
 147   fhEtaPhiNeg->SetYTitle("#phi (rad)");
 148   outputContainer->Add(fhEtaPhiNeg);
 149
 150   fhProductionVertexBC      = new TH1F("hProductionVertexBC", "tracks production vertex bunch crossing ", 18 , -9 , 9 ) ;
 151   fhProductionVertexBC->SetYTitle("# tracks");
 152   fhProductionVertexBC->SetXTitle("Bunch crossing");
 153   outputContainer->Add(fhProductionVertexBC);
 154
 155   Int_t ntofbins = 1000;
 156   Int_t mintof = -500;
 157   Int_t maxtof =  500;
 158
 159   fhTOFSignal  = new TH1F ("hTOFSignal","TOF signal", ntofbins,mintof,maxtof);
 160   fhTOFSignal->SetXTitle("TOF signal (ns)");
 161   outputContainer->Add(fhTOFSignal);
 162
 163   fhTOFSignalBCOK  = new TH1F ("hTOFSignalBCOK","TOF signal", ntofbins,mintof,maxtof);
 164   fhTOFSignalBCOK->SetXTitle("TOF signal (ns)");
 165   outputContainer->Add(fhTOFSignalBCOK);
 166
 167   fhTOFSignalPtCut  = new TH1F ("hTOFSignalPtCut","TOF signal", ntofbins,mintof,maxtof);
 168   fhTOFSignalPtCut->SetXTitle("TOF signal (ns)");
 169   outputContainer->Add(fhTOFSignalPtCut);
 170
 171   fhPtTOFSignal  = new TH2F ("hPtTOFSignal","TOF signal", nptbins,ptmin,ptmax,ntofbins,mintof,maxtof);
 172   fhPtTOFSignal->SetYTitle("TOF signal (ns)");
 173   fhPtTOFSignal->SetXTitle("p_{T} (GeV/c)");
 174   outputContainer->Add(fhPtTOFSignal);
 175
 176   if(fFillPileUpHistograms)
 177   {
 178     TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
 179
 180     for(Int_t i = 0 ; i < 7 ; i++)
 181     {
 182       fhPtPileUp[i]  = new TH1F(Form("hPtPileUp%s",pileUpName[i].Data()),
 183                                 Form("Track p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()),
 184                                 nptbins,ptmin,ptmax);
 185       fhPtPileUp[i]->SetXTitle("p_{T} (GeV/c)");
 186       outputContainer->Add(fhPtPileUp[i]);
 187
 188       fhPtTOFSignalPileUp[i]  = new TH2F(Form("hPtTOFSignalPileUp%s",pileUpName[i].Data()),
 189                                          Form("Track TOF vs p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()),
 190                                          nptbins,ptmin,ptmax,ntofbins,mintof,maxtof);
 191       fhPtTOFSignalPileUp[i]->SetXTitle("p_{T} (GeV/c)");
 192       fhPtTOFSignalPileUp[i]->SetXTitle("TOF signal (ns)");
 193       outputContainer->Add(fhPtTOFSignalPileUp[i]);
 194     }
 195
 196     fhEtaPhiTOFBC0  = new TH2F ("hEtaPhiTOFBC0","eta-phi for tracks with hit on TOF, and tof corresponding to BC=0",netabins,etamin,etamax, nphibins,phimin,phimax);
 197     fhEtaPhiTOFBC0->SetXTitle("#eta ");
 198     fhEtaPhiTOFBC0->SetYTitle("#phi (rad)");
 199     outputContainer->Add(fhEtaPhiTOFBC0);
 200
 201     fhEtaPhiTOFBCPlus  = new TH2F ("hEtaPhiTOFBCPlus","eta-phi for tracks with hit on TOF, and tof corresponding to BC>0",netabins,etamin,etamax, nphibins,phimin,phimax);
 202     fhEtaPhiTOFBCPlus->SetXTitle("#eta ");
 203     fhEtaPhiTOFBCPlus->SetYTitle("#phi (rad)");
 204     outputContainer->Add(fhEtaPhiTOFBCPlus);
 205
 206     fhEtaPhiTOFBCMinus  = new TH2F ("hEtaPhiTOFBCMinus","eta-phi for tracks with hit on TOF, and tof corresponding to BC<0",netabins,etamin,etamax, nphibins,phimin,phimax);
 207     fhEtaPhiTOFBCMinus->SetXTitle("#eta ");
 208     fhEtaPhiTOFBCMinus->SetYTitle("#phi (rad)");
 209     outputContainer->Add(fhEtaPhiTOFBCMinus);
 210
 211     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);
 212     fhEtaPhiTOFBC0PileUpSPD->SetXTitle("#eta ");
 213     fhEtaPhiTOFBC0PileUpSPD->SetYTitle("#phi (rad)");
 214     outputContainer->Add(fhEtaPhiTOFBC0PileUpSPD);
 215
 216     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);
 217     fhEtaPhiTOFBCPlusPileUpSPD->SetXTitle("#eta ");
 218     fhEtaPhiTOFBCPlusPileUpSPD->SetYTitle("#phi (rad)");
 219     outputContainer->Add(fhEtaPhiTOFBCPlusPileUpSPD);
 220
 221     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);
 222     fhEtaPhiTOFBCMinusPileUpSPD->SetXTitle("#eta ");
 223     fhEtaPhiTOFBCMinusPileUpSPD->SetYTitle("#phi (rad)");
 224     outputContainer->Add(fhEtaPhiTOFBCMinusPileUpSPD);
 225
 226   }
 227
 228   fhPtTOFStatus0  = new TH1F ("hPtTOFStatus0","p_T distribution of tracks not hitting TOF", nptbins,ptmin,ptmax);
 229   fhPtTOFStatus0->SetXTitle("p_{T} (GeV/c)");
 230   outputContainer->Add(fhPtTOFStatus0);
 231
 232
 233   fhEtaPhiTOFStatus0  = new TH2F ("hEtaPhiTOFStatus0","eta-phi for tracks without hit on TOF",netabins,etamin,etamax, nphibins,phimin,phimax);
 234   fhEtaPhiTOFStatus0->SetXTitle("#eta ");
 235   fhEtaPhiTOFStatus0->SetYTitle("#phi (rad)");
 236   outputContainer->Add(fhEtaPhiTOFStatus0);
 237
 238   TString dcaName[] = {"xy","z","Cons"} ;
 239   Int_t ndcabins = 800;
 240   Int_t mindca = -4;
 241   Int_t maxdca =  4;
 242
 243   for(Int_t i = 0 ; i < 3 ; i++)
 244   {
 245
 246     fhPtDCA[i]  = new TH2F(Form("hPtDCA%s",dcaName[i].Data()),
 247                            Form("Track DCA%s vs p_{T} distribution",dcaName[i].Data()),
 248                            nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 249     fhPtDCA[i]->SetXTitle("p_{T} (GeV/c)");
 250     fhPtDCA[i]->SetXTitle(Form("DCA_{%s}",dcaName[i].Data()));
 251     outputContainer->Add(fhPtDCA[i]);
 252
 253     fhPtDCATOFBC0[i]  = new TH2F(Form("hPtDCA%sTOFBC0",dcaName[i].Data()),
 254                            Form("Track DCA%s vs p_{T} distribution",dcaName[i].Data()),
 255                            nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 256     fhPtDCATOFBC0[i]->SetXTitle("p_{T} (GeV/c)");
 257     fhPtDCATOFBC0[i]->SetXTitle(Form("DCA_{%s}",dcaName[i].Data()));
 258     outputContainer->Add(fhPtDCATOFBC0[i]);
 259
 260     fhPtDCANoTOFHit[i]  = new TH2F(Form("hPtDCA%sNoTOFHit",dcaName[i].Data()),
 261                            Form("Track (no TOF hit) DCA%s vs p_{T} distribution",dcaName[i].Data()),
 262                            nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 263     fhPtDCANoTOFHit[i]->SetXTitle("p_{T} (GeV/c)");
 264     fhPtDCANoTOFHit[i]->SetXTitle(Form("DCA_{%s}",dcaName[i].Data()));
 265     outputContainer->Add(fhPtDCANoTOFHit[i]);
 266
 267 //    fhPtDCAVtxOutBC0[i]  = new TH2F(Form("hPtDCA%sVtxOutBC0",dcaName[i].Data()),
 268 //                           Form("Track DCA%s vs p_{T} distribution, vertex with BC!=0",dcaName[i].Data()),
 269 //                           nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 270 //    fhPtDCAVtxOutBC0[i]->SetXTitle("p_{T} (GeV/c)");
 271 //    fhPtDCAVtxOutBC0[i]->SetXTitle(Form("DCA_{%s}",dcaName[i].Data()));
 272 //    outputContainer->Add(fhPtDCAVtxOutBC0[i]);
 273 //
 274 //    fhPtDCAVtxOutBC0NoTOFHit[i]  = new TH2F(Form("hPtDCA%sVtxOutBC0NoTOFHit",dcaName[i].Data()),
 275 //                                   Form("Track (no TOF hit) DCA%s vs p_{T} distribution, vertex with BC!=0",dcaName[i].Data()),
 276 //                                   nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 277 //    fhPtDCAVtxOutBC0NoTOFHit[i]->SetXTitle("p_{T} (GeV/c)");
 278 //    fhPtDCAVtxOutBC0NoTOFHit[i]->SetXTitle(Form("DCA_{%s}",dcaName[i].Data()));
 279 //    outputContainer->Add(fhPtDCAVtxOutBC0NoTOFHit[i]);
 280
 281     if(fFillPileUpHistograms)
 282     {
 283       fhPtDCAPileUp[i]  = new TH2F(Form("hPtDCA%sPileUp",dcaName[i].Data()),
 284                              Form("Track DCA%s vs p_{T} distribution, SPD Pile-Up",dcaName[i].Data()),
 285                              nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 286       fhPtDCAPileUp[i]->SetXTitle("p_{T} (GeV/c)");
 287       fhPtDCAPileUp[i]->SetXTitle(Form("DCA_{%s}",dcaName[i].Data()));
 288       outputContainer->Add(fhPtDCAPileUp[i]);
 289
 290       fhPtDCAPileUpTOFBC0[i]  = new TH2F(Form("hPtDCA%sPileUpTOFBC0",dcaName[i].Data()),
 291                                    Form("Track DCA%s vs p_{T} distribution",dcaName[i].Data()),
 292                                    nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 293       fhPtDCAPileUpTOFBC0[i]->SetXTitle("p_{T} (GeV/c)");
 294       fhPtDCAPileUpTOFBC0[i]->SetXTitle(Form("DCA_{%s}",dcaName[i].Data()));
 295       outputContainer->Add(fhPtDCAPileUpTOFBC0[i]);
 296
 297       fhPtDCAPileUpNoTOFHit[i]  = new TH2F(Form("hPtDCA%sPileUpNoTOFHit",dcaName[i].Data()),
 298                                      Form("Track (no TOF hit) DCA%s vs p_{T} distribution, SPD Pile-Up, vertex with BC!=0",dcaName[i].Data()),
 299                                      nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 300       fhPtDCAPileUpNoTOFHit[i]->SetXTitle("p_{T} (GeV/c)");
 301       fhPtDCAPileUpNoTOFHit[i]->SetXTitle(Form("DCA_{%s}",dcaName[i].Data()));
 302       outputContainer->Add(fhPtDCAPileUpNoTOFHit[i]);
 303
 304 //      fhPtDCAVtxOutBC0PileUp[i]  = new TH2F(Form("hPtDCA%sPileUpVtxOutBC0",dcaName[i].Data()),
 305 //                                   Form("Track DCA%s vs p_{T} distribution, SPD Pile-Up",dcaName[i].Data()),
 306 //                                   nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 307 //      fhPtDCAVtxOutBC0PileUp[i]->SetXTitle("p_{T} (GeV/c)");
 308 //      fhPtDCAVtxOutBC0PileUp[i]->SetXTitle(Form("DCA_{%s}",dcaName[i].Data()));
 309 //      outputContainer->Add(fhPtDCAVtxOutBC0PileUp[i]);
 310 //
 311 //      fhPtDCAVtxOutBC0PileUpNoTOFHit[i]  = new TH2F(Form("hPtDCA%sVtxOutBC0PileUpNoTOFHit",dcaName[i].Data()),
 312 //                                           Form("Track (no TOF hit) DCA%s vs p_{T} distribution, SPD Pile-Up, vertex with BC!=0",dcaName[i].Data()),
 313 //                                           nptbins,ptmin,ptmax,ndcabins,mindca,maxdca);
 314 //      fhPtDCAVtxOutBC0PileUpNoTOFHit[i]->SetXTitle("p_{T} (GeV/c)");
 315 //      fhPtDCAVtxOutBC0PileUpNoTOFHit[i]->SetXTitle(Form("DCA_{%s}",dcaName[i].Data()));
 316 //      outputContainer->Add(fhPtDCAVtxOutBC0PileUpNoTOFHit[i]);
 317     }
 318   }
 319
 320
 321
 322   if(IsDataMC()){
 323
 324     fhPtPion  = new TH1F ("hPtMCPion","p_T distribution from #pi", nptbins,ptmin,ptmax);
 325     fhPtPion->SetXTitle("p_{T} (GeV/c)");
 326     outputContainer->Add(fhPtPion);
 327
 328     fhPhiPion  = new TH2F ("hPhiMCPion","#phi distribution from #pi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 329     fhPhiPion->SetXTitle("#phi (rad)");
 330     outputContainer->Add(fhPhiPion);
 331
 332     fhEtaPion  = new TH2F ("hEtaMCPion","#eta distribution from #pi",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 333     fhEtaPion->SetXTitle("#eta ");
 334     outputContainer->Add(fhEtaPion);
 335
 336     fhPtProton  = new TH1F ("hPtMCProton","p_T distribution from proton", nptbins,ptmin,ptmax);
 337     fhPtProton->SetXTitle("p_{T} (GeV/c)");
 338     outputContainer->Add(fhPtProton);
 339
 340     fhPhiProton  = new TH2F ("hPhiMCProton","#phi distribution from proton",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 341     fhPhiProton->SetXTitle("#phi (rad)");
 342     outputContainer->Add(fhPhiProton);
 343
 344     fhEtaProton  = new TH2F ("hEtaMCProton","#eta distribution from proton",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 345     fhEtaProton->SetXTitle("#eta ");
 346     outputContainer->Add(fhEtaProton);
 347
 348     fhPtKaon  = new TH1F ("hPtMCKaon","p_T distribution from kaon", nptbins,ptmin,ptmax);
 349     fhPtKaon->SetXTitle("p_{T} (GeV/c)");
 350     outputContainer->Add(fhPtKaon);
 351
 352     fhPhiKaon  = new TH2F ("hPhiMCKaon","#phi distribution from kaon",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 353     fhPhiKaon->SetXTitle("#phi (rad)");
 354     outputContainer->Add(fhPhiKaon);
 355
 356     fhEtaKaon  = new TH2F ("hEtaMCKaon","#eta distribution from kaon",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 357     fhEtaKaon->SetXTitle("#eta ");
 358     outputContainer->Add(fhEtaKaon);
 359
 360     fhPtElectron  = new TH1F ("hPtMCElectron","p_T distribution from electron", nptbins,ptmin,ptmax);
 361     fhPtElectron->SetXTitle("p_{T} (GeV/c)");
 362     outputContainer->Add(fhPtElectron);
 363
 364     fhPhiElectron  = new TH2F ("hPhiMCElectron","#phi distribution from electron",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 365     fhPhiElectron->SetXTitle("#phi (rad)");
 366     outputContainer->Add(fhPhiElectron);
 367
 368     fhEtaElectron  = new TH2F ("hEtaMCElectron","#eta distribution from electron",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 369     fhEtaElectron->SetXTitle("#eta ");
 370     outputContainer->Add(fhEtaElectron);
 371
 372     fhPtUnknown  = new TH1F ("hPtMCUnknown","p_T distribution from unknown", nptbins,ptmin,ptmax);
 373     fhPtUnknown->SetXTitle("p_{T} (GeV/c)");
 374     outputContainer->Add(fhPtUnknown);
 375
 376     fhPhiUnknown  = new TH2F ("hPhiMCUnknown","#phi distribution from unknown",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 377     fhPhiUnknown->SetXTitle("#phi (rad)");
 378     outputContainer->Add(fhPhiUnknown);
 379
 380     fhEtaUnknown  = new TH2F ("hEtaMCUnknown","#eta distribution from unknown",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 381     fhEtaUnknown->SetXTitle("#eta ");
 382     outputContainer->Add(fhEtaUnknown);
 383
 384   }
 385
 386   return outputContainer;
 387
 388 }
 389
 390 //___________________________________________
 391 void AliAnaChargedParticles::InitParameters()
 392 {
 393   //Initialize the parameters of the analysis.
 394   SetOutputAODClassName("AliAODPWG4Particle");
 395   SetOutputAODName("PWG4Particle");
 396
 397   AddToHistogramsName("AnaCharged_");
 398
 399   fPdg = -1; //Select all tracks
 400
 401 }
 402
 403 //____________________________________________________________
 404 void AliAnaChargedParticles::Print(const Option_t * opt) const
 405 {
 406   //Print some relevant parameters set for the analysis
 407   if(! opt)
 408     return;
 409
 410   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 411   AliAnaCaloTrackCorrBaseClass::Print(" ");
 412
 413   printf("Min Pt = %3.2f\n", GetMinPt());
 414   printf("Max Pt = %3.2f\n", GetMaxPt());
 415   printf("Select clusters with pdg %d \n",fPdg);
 416
 417 }
 418
 419 //_________________________________
 420 void AliAnaChargedParticles::Init()
 421 {
 422   //Init
 423   //Do some checks
 424   if(!GetReader()->IsCTSSwitchedOn()){
 425     printf("AliAnaChargedParticles::Init() - STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!\n");
 426     abort();
 427   }
 428
 429 }
 430
 431 //_________________________________________________
 432 void  AliAnaChargedParticles::MakeAnalysisFillAOD()
 433 {
 434   //Do analysis and fill aods
 435   if(!GetCTSTracks() || GetCTSTracks()->GetEntriesFast() == 0) return ;
 436
 437   Int_t ntracks = GetCTSTracks()->GetEntriesFast();
 438   Double_t vert[3] = {0,0,0}; //vertex ;
 439
 440   //Some prints
 441   if(GetDebug() > 0)
 442     printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);
 443
 444   AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (GetReader()->GetInputEvent());
 445   AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (GetReader()->GetInputEvent());
 446
 447   Double_t bz = GetReader()->GetInputEvent()->GetMagneticField();
 448
 449   //Fill AODParticle with CTS aods
 450   TVector3 p3;
 451   Int_t evtIndex = 0;
 452   for(Int_t i = 0; i < ntracks; i++){
 453
 454     AliVTrack * track =  (AliVTrack*) (GetCTSTracks()->At(i));
 455
 456      //Fill AODParticle after some selection
 457     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 458     p3.SetXYZ(mom[0],mom[1],mom[2]);
 459
 460     //TOF
 461     ULong_t status = track->GetStatus();
 462     Bool_t okTOF = (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ;
 463     Double32_t tof = track->GetTOFsignal()*1e-3;
 464     //if( tof < 0) printf("TOF Signal %e, status %d, pt %f\n", tof,status,status2,p3.Pt());
 465
 466     Double_t dcaCons = -999;
 467     Int_t vtxBC = -999;
 468     if     (esdevent) vtxBC = esdevent->GetPrimaryVertex()->GetBC();
 469     else if(aodevent) vtxBC = aodevent->GetPrimaryVertex()->GetBC();
 470
 471     if(vtxBC!=AliVTrack::kTOFBCNA)printf("BC primary %d",vtxBC);
 472
 473     AliAODTrack * aodTrack = dynamic_cast<AliAODTrack*>(track);
 474     if(aodTrack)
 475     {
 476       dcaCons = aodTrack->DCA();
 477       vtxBC   = aodTrack->GetProdVertex()->GetBC();
 478     }
 479
 480     if(vtxBC!=AliVTrack::kTOFBCNA) printf(" - production %d\n",vtxBC);
 481
 482     fhProductionVertexBC->Fill(vtxBC);
 483
 484     Double_t dca[2]   = {1e6,1e6};
 485     Double_t covar[3] = {1e6,1e6,1e6};
 486     track->PropagateToDCA(GetReader()->GetInputEvent()->GetPrimaryVertex(),bz,100.,dca,covar);
 487
 488     if(dcaCons == -999)
 489     {
 490       fhPtDCA[0]->Fill(p3.Pt(), dca[0]);
 491       fhPtDCA[1]->Fill(p3.Pt(), dca[1]);
 492     }
 493     else fhPtDCA[2]->Fill(p3.Pt(), dcaCons);
 494
 495 //    if(TMath::Abs(vtxBC) > 0 && vtxBC!=AliVTrack::kTOFBCNA)
 496 //    {
 497 //      if(dcaCons == -999)
 498 //      {
 499 //        fhPtDCAVtxOutBC0[0]->Fill(p3.Pt(), dca[0]);
 500 //        fhPtDCAVtxOutBC0[1]->Fill(p3.Pt(), dca[1]);
 501 //      }
 502 //      else
 503 //        fhPtDCAVtxOutBC0[2]->Fill(p3.Pt(), dcaCons);
 504 //    }
 505
 506     if(fFillPileUpHistograms && GetReader()->IsPileUpFromSPD())
 507     {
 508       if(dcaCons == -999)
 509       {
 510         fhPtDCAPileUp[0]->Fill(p3.Pt(), dca[0]);
 511         fhPtDCAPileUp[1]->Fill(p3.Pt(), dca[1]);
 512       }
 513       else
 514         fhPtDCAPileUp[2]->Fill(p3.Pt(), dcaCons);
 515
 516 //      if(TMath::Abs(vtxBC) > 0 && vtxBC!=AliVTrack::kTOFBCNA)
 517 //      {
 518 //        if(dcaCons == -999)
 519 //        {
 520 //          fhPtDCAVtxOutBC0PileUp[0]->Fill(p3.Pt(), dca[0]);
 521 //          fhPtDCAVtxOutBC0PileUp[1]->Fill(p3.Pt(), dca[1]);
 522 //        }
 523 //        else fhPtDCAVtxOutBC0PileUp[2]->Fill(p3.Pt(), dcaCons);
 524 //      }
 525     }
 526
 527     if(!okTOF)
 528     {
 529       if(dcaCons == -999)
 530       {
 531         fhPtDCANoTOFHit[0]->Fill(p3.Pt(), dca[0]);
 532         fhPtDCANoTOFHit[1]->Fill(p3.Pt(), dca[1]);
 533       }
 534       else
 535         fhPtDCANoTOFHit[2]->Fill(p3.Pt(), dcaCons);
 536
 537 //      if(TMath::Abs(vtxBC) > 0 && vtxBC!=AliVTrack::kTOFBCNA)
 538 //      {
 539 //        if(dcaCons == -999)
 540 //        {
 541 //          fhPtDCAVtxOutBC0NoTOFHit[0]->Fill(p3.Pt(), dca[0]);
 542 //          fhPtDCAVtxOutBC0NoTOFHit[1]->Fill(p3.Pt(), dca[1]);
 543 //        }
 544 //        else
 545 //          fhPtDCAVtxOutBC0NoTOFHit[2]->Fill(p3.Pt(), dcaCons);
 546 //      }
 547
 548       if(fFillPileUpHistograms && GetReader()->IsPileUpFromSPD())
 549       {
 550         if(dcaCons == -999)
 551         {
 552           fhPtDCAPileUpNoTOFHit[0]->Fill(p3.Pt(), dca[0]);
 553           fhPtDCAPileUpNoTOFHit[1]->Fill(p3.Pt(), dca[1]);
 554         }
 555         else
 556           fhPtDCAPileUpNoTOFHit[2]->Fill(p3.Pt(), dcaCons);
 557
 558 //        if(TMath::Abs(vtxBC) > 0 && vtxBC!=AliVTrack::kTOFBCNA)
 559 //        {
 560 //          if(dcaCons == -999)
 561 //          {
 562 //            fhPtDCAVtxOutBC0PileUpNoTOFHit[0]->Fill(p3.Pt(), dca[0]);
 563 //            fhPtDCAVtxOutBC0PileUpNoTOFHit[1]->Fill(p3.Pt(), dca[1]);
 564 //          }
 565 //          else
 566 //            fhPtDCAVtxOutBC0PileUpNoTOFHit[2]->Fill(p3.Pt(), dcaCons);
 567 //        }
 568       }
 569     }
 570
 571     //printf("track pT %2.2f, DCA Cons %f, DCA1 %f, DCA2 %f, TOFBC %d, oktof %d, tof %f\n",
 572     //      p3.Pt(),dcaCons,dca[0],dca[1],track->GetTOFBunchCrossing(bz),okTOF, tof);
 573
 574     Int_t trackBC = track->GetTOFBunchCrossing(bz);
 575
 576     if(okTOF && trackBC==0)
 577     {
 578       fhTOFSignalBCOK->Fill(tof);
 579
 580       if(dcaCons == -999)
 581       {
 582         fhPtDCATOFBC0[0]->Fill(p3.Pt(), dca[0]);
 583         fhPtDCATOFBC0[1]->Fill(p3.Pt(), dca[1]);
 584       }
 585       else
 586         fhPtDCATOFBC0[2]->Fill(p3.Pt(), dcaCons);
 587
 588       if(fFillPileUpHistograms && GetReader()->IsPileUpFromSPD())
 589       {
 590         if(dcaCons == -999)
 591         {
 592           fhPtDCAPileUpTOFBC0[0]->Fill(p3.Pt(), dca[0]);
 593           fhPtDCAPileUpTOFBC0[1]->Fill(p3.Pt(), dca[1]);
 594         }
 595         else
 596           fhPtDCAPileUpTOFBC0[2]->Fill(p3.Pt(), dcaCons);
 597       }
 598     }
 599
 600     if(okTOF && fFillPileUpHistograms)
 601     {
 602       fhTOFSignal  ->Fill(tof);
 603       fhPtTOFSignal->Fill(p3.Pt(), tof);
 604
 605       if(GetReader()->IsPileUpFromSPD())               fhPtTOFSignalPileUp[0]->Fill(p3.Pt(), tof);
 606       if(GetReader()->IsPileUpFromEMCal())             fhPtTOFSignalPileUp[1]->Fill(p3.Pt(), tof);
 607       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtTOFSignalPileUp[2]->Fill(p3.Pt(), tof);
 608       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtTOFSignalPileUp[3]->Fill(p3.Pt(), tof);
 609       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtTOFSignalPileUp[4]->Fill(p3.Pt(), tof);
 610       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtTOFSignalPileUp[5]->Fill(p3.Pt(), tof);
 611       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtTOFSignalPileUp[6]->Fill(p3.Pt(), tof);
 612
 613       if      (trackBC ==0)  { fhEtaPhiTOFBC0    ->Fill(track->Eta(),track->Phi()); if(GetReader()->IsPileUpFromSPD()) fhEtaPhiTOFBC0PileUpSPD    ->Fill(track->Eta(),track->Phi()); }
 614       else if (trackBC < 0)  { fhEtaPhiTOFBCPlus ->Fill(track->Eta(),track->Phi()); if(GetReader()->IsPileUpFromSPD()) fhEtaPhiTOFBCPlusPileUpSPD ->Fill(track->Eta(),track->Phi()); }
 615       else if (trackBC > 0)  { fhEtaPhiTOFBCMinus->Fill(track->Eta(),track->Phi()); if(GetReader()->IsPileUpFromSPD()) fhEtaPhiTOFBCMinusPileUpSPD->Fill(track->Eta(),track->Phi()); }
 616
 617     }
 618
 619     Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
 620
 621     if(GetDebug() > 1)
 622       printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, eta %2.2f in fiducial cut %d\n",p3.Pt(), p3.Phi(), p3.Eta(),in);
 623
 624     //Acceptance selection
 625     if(IsFiducialCutOn() && ! in ) continue ;
 626
 627     // Momentum selection
 628     if(track->Pt() < GetMinPt() || track->Pt() > GetMaxPt()) continue;
 629
 630     if(okTOF) fhTOFSignalPtCut->Fill(tof);
 631     else
 632     {
 633       fhPtTOFStatus0    ->Fill(track->Pt());
 634       fhEtaPhiTOFStatus0->Fill(track->Eta(),track->Phi());
 635     }
 636
 637     //Keep only particles identified with fPdg
 638     //Selection not done for the moment
 639     //Should be done here.
 640
 641     // Mixed event
 642     if (GetMixedEvent())
 643     {
 644       evtIndex = GetMixedEvent()->EventIndex(track->GetID()) ;
 645     }
 646
 647     GetVertex(vert,evtIndex);
 648     if(TMath::Abs(vert[2])> GetZvertexCut()) return;
 649
 650     AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
 651     tr.SetDetector("CTS");
 652     tr.SetLabel(track->GetLabel());
 653     tr.SetTrackLabel(track->GetID(),-1);
 654     tr.SetChargedBit(track->Charge()>0);
 655
 656     AddAODParticle(tr);
 657
 658   }//loop
 659
 660   if(GetDebug() > 0)
 661     printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
 662 }
 663
 664 //__________________________________________________________________
 665 void  AliAnaChargedParticles::MakeAnalysisFillHistograms()
 666 {
 667   //Do analysis and fill histograms
 668
 669   //Loop on stored AODParticles
 670   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 671
 672   fhNtracks->Fill(GetReader()->GetTrackMultiplicity()) ;
 673
 674   if(GetDebug() > 0)
 675     printf("AliAnaChargedParticles::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 676
 677
 678   for(Int_t iaod = 0; iaod < naod ; iaod++){
 679     AliAODPWG4Particle* tr =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 680
 681     fhPt->Fill(tr->Pt());
 682
 683     if(tr->GetChargedBit()){
 684       fhPhiPos   ->Fill(tr->Pt(), tr->Phi());
 685       fhEtaPos   ->Fill(tr->Pt(), tr->Eta());
 686       fhEtaPhiPos->Fill(tr->Eta(),tr->Phi());
 687     }
 688     else{
 689       fhPhiNeg   ->Fill(tr->Pt(), tr->Phi());
 690       fhEtaNeg   ->Fill(tr->Pt(), tr->Eta());
 691       fhEtaPhiNeg->Fill(tr->Eta(),tr->Phi());
 692     }
 693
 694     if(fFillPileUpHistograms)
 695     {
 696       if(GetReader()->IsPileUpFromSPD())               {fhPtPileUp[0]->Fill(tr->Pt());}
 697       if(GetReader()->IsPileUpFromEMCal())             {fhPtPileUp[1]->Fill(tr->Pt());}
 698       if(GetReader()->IsPileUpFromSPDOrEMCal())        {fhPtPileUp[2]->Fill(tr->Pt());}
 699       if(GetReader()->IsPileUpFromSPDAndEMCal())       {fhPtPileUp[3]->Fill(tr->Pt());}
 700       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    {fhPtPileUp[4]->Fill(tr->Pt());}
 701       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    {fhPtPileUp[5]->Fill(tr->Pt());}
 702       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) {fhPtPileUp[6]->Fill(tr->Pt());}
 703     }
 704
 705
 706     if(IsDataMC())
 707     {
 708       //Play with the MC stack if available
 709       Int_t mompdg = -1;
 710       Int_t label  = tr->GetLabel();
 711       if(label >= 0)
 712       {
 713         if( GetReader()->ReadStack() && label < GetMCStack()->GetNtrack())
 714         {
 715           TParticle * mom = GetMCStack()->Particle(label);
 716           mompdg =TMath::Abs(mom->GetPdgCode());
 717         }
 718         else if(GetReader()->ReadAODMCParticles())
 719         {
 720           AliAODMCParticle * aodmom = 0;
 721           //Get the list of MC particles
 722           aodmom = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(tr->GetInputFileIndex()))->At(label);
 723           mompdg =TMath::Abs(aodmom->GetPdgCode());
 724         }
 725       }
 726
 727       if(mompdg==211)
 728       {
 729         fhPtPion ->Fill(tr->Pt());
 730         fhPhiPion->Fill(tr->Pt(), tr->Phi());
 731         fhEtaPion->Fill(tr->Pt(), tr->Eta());
 732       }
 733       else if(mompdg==2212)
 734       {
 735         fhPtProton ->Fill(tr->Pt());
 736         fhPhiProton->Fill(tr->Pt(), tr->Phi());
 737         fhEtaProton->Fill(tr->Pt(), tr->Eta());
 738       }
 739       else if(mompdg==321)
 740       {
 741         fhPtKaon ->Fill(tr->Pt());
 742         fhPhiKaon->Fill(tr->Pt(), tr->Phi());
 743         fhEtaKaon->Fill(tr->Pt(), tr->Eta());
 744       }
 745       else if(mompdg==11)
 746       {
 747         fhPtElectron ->Fill(tr->Pt());
 748         fhPhiElectron->Fill(tr->Pt(), tr->Phi());
 749         fhEtaElectron->Fill(tr->Pt(), tr->Eta());
 750       }
 751       else {
 752         //printf("unknown pdg %d\n",mompdg);
 753         fhPtUnknown ->Fill(tr->Pt());
 754         fhPhiUnknown->Fill(tr->Pt(), tr->Phi());
 755         fhEtaUnknown->Fill(tr->Pt(), tr->Eta());
 756       }
 757
 758     }//Work with stack also
 759
 760   }// aod branch loop
 761
 762 }