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