PWG4/PartCorrDep/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 /* $Id: $ */
  16
  17 //_________________________________________________________________________
  18 //
  19 // Class for track selection and identification (not done now)
  20 // Tracks from the CTS are kept in the AOD.
  21 // Few histograms produced.
  22 //
  23 //-- Author: Gustavo Conesa (INFN-LNF)
  24 //_________________________________________________________________________
  25
  26
  27 // --- ROOT system ---
  28 #include "TParticle.h"
  29 #include "TH2F.h"
  30 #include "TH3D.h"
  31 //---- AliRoot system ----
  32 #include "AliAnaChargedParticles.h"
  33 #include "AliCaloTrackReader.h"
  34 #include "AliAODPWG4Particle.h"
  35 #include "AliStack.h"
  36 #include "AliFiducialCut.h"
  37 #include "AliAODTrack.h"
  38 #include "AliAODMCParticle.h"
  39
  40 ClassImp(AliAnaChargedParticles)
  41
  42 //____________________________________________________________________________
  43   AliAnaChargedParticles::AliAnaChargedParticles() :
  44     AliAnaPartCorrBaseClass(),fPdg(0), fhNtracks(0),fhVertex(0), fhPt(0),fhPhi(0),fhEta(0),
  45     fhPtEtaPhiPos(0), fhPtEtaPhiNeg(0),
  46     fhPtPion(0),fhPhiPion(0),fhEtaPion(0),
  47     fhPtProton(0),fhPhiProton(0),fhEtaProton(0),
  48     fhPtElectron(0),fhPhiElectron(0),fhEtaElectron(0),
  49     fhPtKaon(0),fhPhiKaon(0),fhEtaKaon(0),
  50     fhPtUnknown(0),fhPhiUnknown(0),fhEtaUnknown(0)
  51
  52 {
  53   //Default Ctor
  54
  55   //Initialize parameters
  56   InitParameters();
  57 }
  58
  59 /*
  60 //____________________________________________________________________________
  61 AliAnaChargedParticles::AliAnaChargedParticles(const AliAnaChargedParticles & ch) :
  62   AliAnaPartCorrBaseClass(ch), fPdg(ch.fPdg),  fhPt(ch.fhPt),  fhPhi(ch.fhPhi),fhEta(ch.fhEta),
  63   fhPtPion(ch. fhPtPion),fhPhiPion(ch.fhPhiPion),fhEtaPion(ch.fhEtaPion),
  64   fhPtProton(ch.fhPtProton),fhPhiProton(ch.fhPhiProton),fhEtaProton(ch.fhEtaProton),
  65   fhPtElectron(ch. fhPtElectron),fhPhiElectron(ch.fhPhiElectron),fhEtaElectron(ch.fhEtaElectron),
  66   fhPtKaon(ch. fhPtKaon),fhPhiKaon(ch.fhPhiKaon),fhEtaKaon(ch.fhEtaKaon),
  67   fhPtUnknown(ch.fhPtUnknown),fhPhiUnknown(ch.fhPhiUnknown),fhEtaUnknown(ch.fhEtaUnknown)
  68
  69 {
  70   // cpy ctor
  71
  72 }
  73
  74 //_________________________________________________________________________
  75 AliAnaChargedParticles & AliAnaChargedParticles::operator = (const AliAnaChargedParticles & ch)
  76 {
  77   // assignment operator
  78
  79   if(this == &ch)return *this;
  80   ((AliAnaPartCorrBaseClass *)this)->operator=(ch);
  81
  82   fPdg = ch.fPdg;
  83   fhPt = ch.fhPt;
  84   fhPhi = ch.fhPhi;
  85   fhEta = ch.fhEta;
  86
  87   fhPtPion = ch. fhPtPion; fhPhiPion = ch.fhPhiPion; fhEtaPion = ch.fhEtaPion;
  88   fhPtProton = ch.fhPtProton; fhPhiProton = ch.fhPhiProton; fhEtaProton = ch.fhEtaProton;
  89   fhPtElectron = ch. fhPtElectron; fhPhiElectron = ch.fhPhiElectron; fhEtaElectron = ch.fhEtaElectron;
  90   fhPtKaon = ch. fhPtKaon; fhPhiKaon = ch.fhPhiKaon; fhEtaKaon = ch.fhEtaKaon;
  91   fhPtUnknown = ch.fhPtUnknown; fhPhiUnknown = ch.fhPhiUnknown; fhEtaUnknown = ch.fhEtaUnknown;
  92
  93   return *this;
  94
  95 }
  96 */
  97
  98 //________________________________________________________________________
  99 TList *  AliAnaChargedParticles::GetCreateOutputObjects()
 100 {
 101   // Create histograms to be saved in output file and
 102   // store them in fOutputContainer
 103
 104
 105   TList * outputContainer = new TList() ;
 106   outputContainer->SetName("ExampleHistos") ;
 107
 108   Int_t nptbins  = GetHistoPtBins();
 109   Int_t nphibins = GetHistoPhiBins();
 110   Int_t netabins = GetHistoEtaBins();
 111   Float_t ptmax  = GetHistoPtMax();
 112   Float_t phimax = GetHistoPhiMax();
 113   Float_t etamax = GetHistoEtaMax();
 114   Float_t ptmin  = GetHistoPtMin();
 115   Float_t phimin = GetHistoPhiMin();
 116   Float_t etamin = GetHistoEtaMin();
 117
 118   fhNtracks  = new TH1F ("hNtracks","# of tracks", 1000,0,1000);
 119   fhNtracks->SetXTitle("# of tracks");
 120   outputContainer->Add(fhNtracks);
 121
 122   fhVertex  = new TH3D ("Vertex","vertex position", 100,-50.,50., 100,-50.,50., 100,-50.,50.);
 123   fhVertex->SetXTitle("X");
 124   fhVertex->SetYTitle("Y");
 125   fhVertex->SetZTitle("Z");
 126   outputContainer->Add(fhVertex);
 127
 128   fhPt  = new TH1F ("hPtCharged","p_T distribution", nptbins,ptmin,ptmax);
 129   fhPt->SetXTitle("p_{T} (GeV/c)");
 130   outputContainer->Add(fhPt);
 131
 132   fhPhi  = new TH2F ("hPhiCharged","#phi distribution",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 133   fhPhi->SetXTitle("#phi (rad)");
 134   outputContainer->Add(fhPhi);
 135
 136   fhEta  = new TH2F ("hEtaCharged","#eta distribution",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 137   fhEta->SetXTitle("#eta ");
 138   outputContainer->Add(fhEta);
 139
 140   fhPtEtaPhiPos  = new TH3D ("hPtEtaPhiPos","pt/eta/phi of positive charge",nptbins,ptmin,ptmax, netabins,etamin,etamax, nphibins,phimin,phimax);
 141   fhPtEtaPhiPos->SetXTitle("p_{T}^{h^{+}} (GeV/c)");
 142   fhPtEtaPhiPos->SetYTitle("#eta ");
 143   fhPtEtaPhiPos->SetZTitle("#phi (rad)");
 144   outputContainer->Add(fhPtEtaPhiPos);
 145
 146   fhPtEtaPhiNeg  = new TH3D ("hPtEtaPhiNeg","pt/eta/phi of negative charge",nptbins,ptmin,ptmax, netabins,etamin,etamax, nphibins,phimin,phimax);
 147   fhPtEtaPhiNeg->SetXTitle("p_{T}^{h^{-}} (GeV/c)");
 148   fhPtEtaPhiNeg->SetYTitle("#eta ");
 149   fhPtEtaPhiNeg->SetZTitle("#phi (rad)");
 150   outputContainer->Add(fhPtEtaPhiNeg);
 151
 152   if(IsDataMC()){
 153
 154     fhPtPion  = new TH1F ("hPtMCPion","p_T distribution from #pi", nptbins,ptmin,ptmax);
 155     fhPtPion->SetXTitle("p_{T} (GeV/c)");
 156     outputContainer->Add(fhPtPion);
 157
 158     fhPhiPion  = new TH2F ("hPhiMCPion","#phi distribution from #pi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 159     fhPhiPion->SetXTitle("#phi (rad)");
 160     outputContainer->Add(fhPhiPion);
 161
 162     fhEtaPion  = new TH2F ("hEtaMCPion","#eta distribution from #pi",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 163     fhEtaPion->SetXTitle("#eta ");
 164     outputContainer->Add(fhEtaPion);
 165
 166     fhPtProton  = new TH1F ("hPtMCProton","p_T distribution from proton", nptbins,ptmin,ptmax);
 167     fhPtProton->SetXTitle("p_{T} (GeV/c)");
 168     outputContainer->Add(fhPtProton);
 169
 170     fhPhiProton  = new TH2F ("hPhiMCProton","#phi distribution from proton",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 171     fhPhiProton->SetXTitle("#phi (rad)");
 172     outputContainer->Add(fhPhiProton);
 173
 174     fhEtaProton  = new TH2F ("hEtaMCProton","#eta distribution from proton",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 175     fhEtaProton->SetXTitle("#eta ");
 176     outputContainer->Add(fhEtaProton);
 177
 178     fhPtKaon  = new TH1F ("hPtMCKaon","p_T distribution from kaon", nptbins,ptmin,ptmax);
 179     fhPtKaon->SetXTitle("p_{T} (GeV/c)");
 180     outputContainer->Add(fhPtKaon);
 181
 182     fhPhiKaon  = new TH2F ("hPhiMCKaon","#phi distribution from kaon",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 183     fhPhiKaon->SetXTitle("#phi (rad)");
 184     outputContainer->Add(fhPhiKaon);
 185
 186     fhEtaKaon  = new TH2F ("hEtaMCKaon","#eta distribution from kaon",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 187     fhEtaKaon->SetXTitle("#eta ");
 188     outputContainer->Add(fhEtaKaon);
 189
 190     fhPtElectron  = new TH1F ("hPtMCElectron","p_T distribution from electron", nptbins,ptmin,ptmax);
 191     fhPtElectron->SetXTitle("p_{T} (GeV/c)");
 192     outputContainer->Add(fhPtElectron);
 193
 194     fhPhiElectron  = new TH2F ("hPhiMCElectron","#phi distribution from electron",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 195     fhPhiElectron->SetXTitle("#phi (rad)");
 196     outputContainer->Add(fhPhiElectron);
 197
 198     fhEtaElectron  = new TH2F ("hEtaMCElectron","#eta distribution from electron",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 199     fhEtaElectron->SetXTitle("#eta ");
 200     outputContainer->Add(fhEtaElectron);
 201
 202     fhPtUnknown  = new TH1F ("hPtMCUnknown","p_T distribution from unknown", nptbins,ptmin,ptmax);
 203     fhPtUnknown->SetXTitle("p_{T} (GeV/c)");
 204     outputContainer->Add(fhPtUnknown);
 205
 206     fhPhiUnknown  = new TH2F ("hPhiMCUnknown","#phi distribution from unknown",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 207     fhPhiUnknown->SetXTitle("#phi (rad)");
 208     outputContainer->Add(fhPhiUnknown);
 209
 210     fhEtaUnknown  = new TH2F ("hEtaMCUnknown","#eta distribution from unknown",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 211     fhEtaUnknown->SetXTitle("#eta ");
 212     outputContainer->Add(fhEtaUnknown);
 213
 214   }
 215
 216   return outputContainer;
 217
 218 }
 219
 220 //__________________________________________________
 221 void AliAnaChargedParticles::InitParameters()
 222 {
 223   //Initialize the parameters of the analysis.
 224   SetOutputAODClassName("AliAODPWG4Particle");
 225   SetOutputAODName("PWG4Particle");
 226
 227   AddToHistogramsName("AnaCharged_");
 228
 229   fPdg = -1; //Select all tracks
 230
 231 }
 232
 233 //__________________________________________________________________
 234 void AliAnaChargedParticles::Print(const Option_t * opt) const
 235 {
 236   //Print some relevant parameters set for the analysis
 237   if(! opt)
 238     return;
 239
 240   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 241   AliAnaPartCorrBaseClass::Print(" ");
 242
 243   printf("Min Pt = %3.2f\n", GetMinPt());
 244   printf("Max Pt = %3.2f\n", GetMaxPt());
 245   printf("Select clusters with pdg %d \n",fPdg);
 246
 247 }
 248
 249 //____________________________________________________________________________
 250 void AliAnaChargedParticles::Init()
 251 {
 252   //Init
 253   //Do some checks
 254   if(!GetReader()->IsCTSSwitchedOn()){
 255     printf("AliAnaChargedParticles::Init() - STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!\n");
 256     abort();
 257   }
 258
 259 }
 260
 261 //__________________________________________________________________
 262 void  AliAnaChargedParticles::MakeAnalysisFillAOD()
 263 {
 264   //Do analysis and fill aods
 265   if(!GetAODCTS() || GetAODCTS()->GetEntriesFast() == 0) return ;
 266   Int_t ntracks = GetAODCTS()->GetEntriesFast();
 267   Double_t vert[3] = {0,0,0}; //vertex ;
 268   //Some prints
 269   if(GetDebug() > 0)
 270     printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);
 271
 272   //Fill AODParticle with CTS aods
 273   TVector3 p3;
 274   Int_t evtIndex = 0;
 275   for(Int_t i = 0; i < ntracks; i++){
 276
 277     AliAODTrack * track =  (AliAODTrack*) (GetAODCTS()->At(i));
 278
 279     //Fill AODParticle after some selection
 280     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 281     p3.SetXYZ(mom[0],mom[1],mom[2]);
 282
 283     //Acceptance selection
 284     Bool_t in =  GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
 285     if(GetDebug() > 1) printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, in fiducial cut %d\n",p3.Pt(), p3.Phi(), in);
 286     if(p3.Pt() > GetMinPt() && in) {
 287       //Keep only particles identified with fPdg
 288       //Selection not done for the moment
 289       //Should be done here.
 290       if (GetMixedEvent()){
 291         evtIndex = GetMixedEvent()->EventIndex(track->GetID()) ;
 292       }
 293       GetVertex(vert,evtIndex);
 294       if(TMath::Abs(vert[2])> GetZvertexCut()) return;
 295       AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
 296       tr.SetDetector("CTS");
 297       tr.SetLabel(track->GetLabel());
 298       tr.SetTrackLabel(track->GetID(),-1);
 299     //  tr.SetChargedBit(track->Charge());
 300           //Input from second AOD?
 301           //if(GetReader()->GetAODCTSNormalInputEntries() <= i) tr.SetInputFileIndex(1);
 302
 303       AddAODParticle(tr);
 304     }//selection
 305   }//loop
 306
 307   if(GetDebug() > 0)
 308     printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
 309 }
 310
 311 //__________________________________________________________________
 312 void  AliAnaChargedParticles::MakeAnalysisFillHistograms()
 313 {
 314   //Do analysis and fill histograms
 315
 316   //Loop on stored AODParticles
 317   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 318   Double_t v[3] = {0,0,0}; //vertex ;
 319   GetReader()->GetVertex(v);
 320   fhVertex->Fill(v[0],v[1],v[2]);
 321   if(TMath::Abs(v[2]) >GetZvertexCut()) return ;
 322   fhNtracks->Fill(GetReader()->GetTrackMultiplicity()) ;
 323   if(GetDebug() > 0) printf("AliAnaChargedParticles::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 324   for(Int_t iaod = 0; iaod < naod ; iaod++){
 325     AliAODPWG4Particle* tr =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 326
 327     fhPt->Fill(tr->Pt());
 328     fhPhi->Fill(tr->Pt(), tr->Phi());
 329     fhEta->Fill(tr->Pt(), tr->Eta());
 330     //for charge information
 331     AliAODTrack * track =  (AliAODTrack*) (GetAODCTS()->At(iaod));
 332     if(track->Charge()>0)fhPtEtaPhiPos->Fill(tr->Pt(), tr->Eta(),tr->Phi());
 333     if(track->Charge()<0)fhPtEtaPhiNeg->Fill(tr->Pt(), tr->Eta(),tr->Phi());
 334
 335     if(IsDataMC()){
 336       //Play with the MC stack if available
 337       Int_t mompdg = -1;
 338       Int_t label  = TMath::Abs(tr->GetLabel());
 339       if(GetReader()->ReadStack()){
 340         TParticle * mom = GetMCStack()->Particle(label);
 341         mompdg =TMath::Abs(mom->GetPdgCode());
 342       }
 343       else if(GetReader()->ReadAODMCParticles()){
 344         AliAODMCParticle * aodmom = 0;
 345         //Get the list of MC particles
 346         aodmom = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(tr->GetInputFileIndex()))->At(label);
 347         mompdg =TMath::Abs(aodmom->GetPdgCode());
 348       }
 349
 350
 351       if(mompdg==211){
 352         fhPtPion->Fill(tr->Pt());
 353         fhPhiPion->Fill(tr->Pt(), tr->Phi());
 354         fhEtaPion->Fill(tr->Pt(), tr->Eta());
 355       }
 356       else if(mompdg==2212){
 357         fhPtProton->Fill(tr->Pt());
 358         fhPhiProton->Fill(tr->Pt(), tr->Phi());
 359         fhEtaProton->Fill(tr->Pt(), tr->Eta());
 360       }
 361       else if(mompdg==321){
 362         fhPtKaon->Fill(tr->Pt());
 363         fhPhiKaon->Fill(tr->Pt(), tr->Phi());
 364         fhEtaKaon->Fill(tr->Pt(), tr->Eta());
 365       }
 366       else if(mompdg==11){
 367         fhPtElectron->Fill(tr->Pt());
 368         fhPhiElectron->Fill(tr->Pt(), tr->Phi());
 369         fhEtaElectron->Fill(tr->Pt(), tr->Eta());
 370       }
 371       else {
 372         //printf("unknown pdg %d\n",mompdg);
 373         fhPtUnknown->Fill(tr->Pt());
 374         fhPhiUnknown->Fill(tr->Pt(), tr->Phi());
 375         fhEtaUnknown->Fill(tr->Pt(), tr->Eta());
 376       }
 377     }//Work with stack also
 378
 379   }// aod branch loop
 380
 381 }