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
  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     fhTOFSignal(0),   fhTOFSignalPtCut(0),
  56     fhPtTOFSignal(0), fhPtTOFStatus0(0), fhEtaPhiTOFStatus0(0)
  57 {
  58   //Default Ctor
  59
  60   for(Int_t i = 0; i < 7; i++) fhPtTOFSignalPileUp[i] = 0;
  61
  62   //Initialize parameters
  63   InitParameters();
  64
  65 }
  66
  67 //_______________________________________________________
  68 TList *  AliAnaChargedParticles::GetCreateOutputObjects()
  69 {
  70   // Create histograms to be saved in output file and
  71   // store them in fOutputContainer
  72
  73
  74   TList * outputContainer = new TList() ;
  75   outputContainer->SetName("ExampleHistos") ;
  76
  77   Int_t nptbins  = GetHistogramRanges()->GetHistoPtBins(); Int_t nphibins = GetHistogramRanges()->GetHistoPhiBins(); Int_t netabins = GetHistogramRanges()->GetHistoEtaBins();
  78   Float_t ptmax  = GetHistogramRanges()->GetHistoPtMax();  Float_t phimax = GetHistogramRanges()->GetHistoPhiMax();  Float_t etamax = GetHistogramRanges()->GetHistoEtaMax();
  79   Float_t ptmin  = GetHistogramRanges()->GetHistoPtMin();  Float_t phimin = GetHistogramRanges()->GetHistoPhiMin();  Float_t etamin = GetHistogramRanges()->GetHistoEtaMin();
  80
  81   fhNtracks  = new TH1F ("hNtracks","# of tracks", 1000,0,1000);
  82   fhNtracks->SetXTitle("# of tracks");
  83   outputContainer->Add(fhNtracks);
  84
  85   fhPt  = new TH1F ("hPt","p_T distribution", nptbins,ptmin,ptmax);
  86   fhPt->SetXTitle("p_{T} (GeV/c)");
  87   outputContainer->Add(fhPt);
  88
  89   fhPhiNeg  = new TH2F ("hPhiNegative","#phi of negative charges distribution",
  90                         nptbins,ptmin,ptmax, nphibins,phimin,phimax);
  91   fhPhiNeg->SetYTitle("#phi (rad)");
  92   fhPhiNeg->SetXTitle("p_{T} (GeV/c)");
  93   outputContainer->Add(fhPhiNeg);
  94
  95   fhEtaNeg  = new TH2F ("hEtaNegative","#eta of negative charges distribution",
  96                         nptbins,ptmin,ptmax, netabins,etamin,etamax);
  97   fhEtaNeg->SetYTitle("#eta ");
  98   fhEtaNeg->SetXTitle("p_{T} (GeV/c)");
  99   outputContainer->Add(fhEtaNeg);
 100
 101   fhPhiPos  = new TH2F ("hPhiPositive","#phi of positive charges distribution",
 102                         nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 103   fhPhiPos->SetYTitle("#phi (rad)");
 104   fhPhiPos->SetXTitle("p_{T} (GeV/c)");
 105   outputContainer->Add(fhPhiPos);
 106
 107   fhEtaPos  = new TH2F ("hEtaPositive","#eta of positive charges distribution",
 108                         nptbins,ptmin,ptmax, netabins,etamin,etamax);
 109   fhEtaPos->SetYTitle("#eta ");
 110   fhEtaPos->SetXTitle("p_{T} (GeV/c)");
 111   outputContainer->Add(fhEtaPos);
 112
 113   fhEtaPhiPos  = new TH2F ("hEtaPhiPositive","pt/eta/phi of positive charge",netabins,etamin,etamax, nphibins,phimin,phimax);
 114   fhEtaPhiPos->SetXTitle("#eta ");
 115   fhEtaPhiPos->SetYTitle("#phi (rad)");
 116   outputContainer->Add(fhEtaPhiPos);
 117
 118   fhEtaPhiNeg  = new TH2F ("hEtaPhiNegative","eta vs phi of negative charge",netabins,etamin,etamax, nphibins,phimin,phimax);
 119   fhEtaPhiNeg->SetXTitle("#eta ");
 120   fhEtaPhiNeg->SetYTitle("#phi (rad)");
 121   outputContainer->Add(fhEtaPhiNeg);
 122
 123   Int_t ntofbins = 1000;
 124   Int_t mintof = -500;
 125   Int_t maxtof =  500;
 126
 127   fhTOFSignal  = new TH1F ("hTOFSignal","TOF signal", ntofbins,mintof,maxtof);
 128   fhTOFSignal->SetXTitle("TOF signal #times (ns?)");
 129   outputContainer->Add(fhTOFSignal);
 130
 131   fhTOFSignalPtCut  = new TH1F ("hTOFSignalPtCut","TOF signal", ntofbins,mintof,maxtof);
 132   fhTOFSignalPtCut->SetXTitle("TOF signal (ns?)");
 133   outputContainer->Add(fhTOFSignalPtCut);
 134
 135   fhPtTOFSignal  = new TH2F ("hPtTOFSignal","TOF signal", nptbins,ptmin,ptmax,ntofbins,mintof,maxtof);
 136   fhPtTOFSignal->SetYTitle("TOF signal (ns?)");
 137   fhPtTOFSignal->SetXTitle("p_{T} (GeV/c)");
 138   outputContainer->Add(fhPtTOFSignal);
 139
 140   TString pileUpName[] = {"SPD","EMCAL","SPDOrEMCAL","SPDAndEMCAL","SPDAndNotEMCAL","EMCALAndNotSPD","NotSPDAndNotEMCAL"} ;
 141
 142   for(Int_t i = 0 ; i < 7 ; i++)
 143   {
 144     fhPtTOFSignalPileUp[i]  = new TH2F(Form("hPtTOFSignalPileUp%s",pileUpName[i].Data()),
 145                                        Form("Track TOF vs p_{T} distribution, %s Pile-Up event",pileUpName[i].Data()),
 146                                        nptbins,ptmin,ptmax,ntofbins,mintof,maxtof);
 147     fhPtTOFSignalPileUp[i]->SetXTitle("p_{T} (GeV/c)");
 148     fhPtTOFSignalPileUp[i]->SetXTitle("TOF signal (ns?)");
 149     outputContainer->Add(fhPtTOFSignalPileUp[i]);
 150   }
 151
 152   fhPtTOFStatus0  = new TH1F ("hPtTOFStatus0","p_T distribution of tracks not hitting TOF", nptbins,ptmin,ptmax);
 153   fhPtTOFStatus0->SetXTitle("p_{T} (GeV/c)");
 154   outputContainer->Add(fhPtTOFStatus0);
 155
 156
 157   fhEtaPhiTOFStatus0  = new TH2F ("hEtaPhiTOFStatus0","pt/eta/phi for tracks without hit on TOF",netabins,etamin,etamax, nphibins,phimin,phimax);
 158   fhEtaPhiTOFStatus0->SetXTitle("#eta ");
 159   fhEtaPhiTOFStatus0->SetYTitle("#phi (rad)");
 160   outputContainer->Add(fhEtaPhiTOFStatus0);
 161
 162
 163   if(IsDataMC()){
 164
 165     fhPtPion  = new TH1F ("hPtMCPion","p_T distribution from #pi", nptbins,ptmin,ptmax);
 166     fhPtPion->SetXTitle("p_{T} (GeV/c)");
 167     outputContainer->Add(fhPtPion);
 168
 169     fhPhiPion  = new TH2F ("hPhiMCPion","#phi distribution from #pi",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 170     fhPhiPion->SetXTitle("#phi (rad)");
 171     outputContainer->Add(fhPhiPion);
 172
 173     fhEtaPion  = new TH2F ("hEtaMCPion","#eta distribution from #pi",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 174     fhEtaPion->SetXTitle("#eta ");
 175     outputContainer->Add(fhEtaPion);
 176
 177     fhPtProton  = new TH1F ("hPtMCProton","p_T distribution from proton", nptbins,ptmin,ptmax);
 178     fhPtProton->SetXTitle("p_{T} (GeV/c)");
 179     outputContainer->Add(fhPtProton);
 180
 181     fhPhiProton  = new TH2F ("hPhiMCProton","#phi distribution from proton",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 182     fhPhiProton->SetXTitle("#phi (rad)");
 183     outputContainer->Add(fhPhiProton);
 184
 185     fhEtaProton  = new TH2F ("hEtaMCProton","#eta distribution from proton",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 186     fhEtaProton->SetXTitle("#eta ");
 187     outputContainer->Add(fhEtaProton);
 188
 189     fhPtKaon  = new TH1F ("hPtMCKaon","p_T distribution from kaon", nptbins,ptmin,ptmax);
 190     fhPtKaon->SetXTitle("p_{T} (GeV/c)");
 191     outputContainer->Add(fhPtKaon);
 192
 193     fhPhiKaon  = new TH2F ("hPhiMCKaon","#phi distribution from kaon",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 194     fhPhiKaon->SetXTitle("#phi (rad)");
 195     outputContainer->Add(fhPhiKaon);
 196
 197     fhEtaKaon  = new TH2F ("hEtaMCKaon","#eta distribution from kaon",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 198     fhEtaKaon->SetXTitle("#eta ");
 199     outputContainer->Add(fhEtaKaon);
 200
 201     fhPtElectron  = new TH1F ("hPtMCElectron","p_T distribution from electron", nptbins,ptmin,ptmax);
 202     fhPtElectron->SetXTitle("p_{T} (GeV/c)");
 203     outputContainer->Add(fhPtElectron);
 204
 205     fhPhiElectron  = new TH2F ("hPhiMCElectron","#phi distribution from electron",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 206     fhPhiElectron->SetXTitle("#phi (rad)");
 207     outputContainer->Add(fhPhiElectron);
 208
 209     fhEtaElectron  = new TH2F ("hEtaMCElectron","#eta distribution from electron",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 210     fhEtaElectron->SetXTitle("#eta ");
 211     outputContainer->Add(fhEtaElectron);
 212
 213     fhPtUnknown  = new TH1F ("hPtMCUnknown","p_T distribution from unknown", nptbins,ptmin,ptmax);
 214     fhPtUnknown->SetXTitle("p_{T} (GeV/c)");
 215     outputContainer->Add(fhPtUnknown);
 216
 217     fhPhiUnknown  = new TH2F ("hPhiMCUnknown","#phi distribution from unknown",nptbins,ptmin,ptmax, nphibins,phimin,phimax);
 218     fhPhiUnknown->SetXTitle("#phi (rad)");
 219     outputContainer->Add(fhPhiUnknown);
 220
 221     fhEtaUnknown  = new TH2F ("hEtaMCUnknown","#eta distribution from unknown",nptbins,ptmin,ptmax, netabins,etamin,etamax);
 222     fhEtaUnknown->SetXTitle("#eta ");
 223     outputContainer->Add(fhEtaUnknown);
 224
 225   }
 226
 227   return outputContainer;
 228
 229 }
 230
 231 //___________________________________________
 232 void AliAnaChargedParticles::InitParameters()
 233 {
 234   //Initialize the parameters of the analysis.
 235   SetOutputAODClassName("AliAODPWG4Particle");
 236   SetOutputAODName("PWG4Particle");
 237
 238   AddToHistogramsName("AnaCharged_");
 239
 240   fPdg = -1; //Select all tracks
 241
 242 }
 243
 244 //____________________________________________________________
 245 void AliAnaChargedParticles::Print(const Option_t * opt) const
 246 {
 247   //Print some relevant parameters set for the analysis
 248   if(! opt)
 249     return;
 250
 251   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 252   AliAnaCaloTrackCorrBaseClass::Print(" ");
 253
 254   printf("Min Pt = %3.2f\n", GetMinPt());
 255   printf("Max Pt = %3.2f\n", GetMaxPt());
 256   printf("Select clusters with pdg %d \n",fPdg);
 257
 258 }
 259
 260 //_________________________________
 261 void AliAnaChargedParticles::Init()
 262 {
 263   //Init
 264   //Do some checks
 265   if(!GetReader()->IsCTSSwitchedOn()){
 266     printf("AliAnaChargedParticles::Init() - STOP!: You want to use CTS tracks in analysis but not read!! \n!!Check the configuration file!!\n");
 267     abort();
 268   }
 269
 270 }
 271
 272 //_________________________________________________
 273 void  AliAnaChargedParticles::MakeAnalysisFillAOD()
 274 {
 275   //Do analysis and fill aods
 276   if(!GetCTSTracks() || GetCTSTracks()->GetEntriesFast() == 0) return ;
 277
 278   Int_t ntracks = GetCTSTracks()->GetEntriesFast();
 279   Double_t vert[3] = {0,0,0}; //vertex ;
 280
 281   //Some prints
 282   if(GetDebug() > 0)
 283     printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - In CTS aod entries %d\n", ntracks);
 284
 285   //Fill AODParticle with CTS aods
 286   TVector3 p3;
 287   Int_t evtIndex = 0;
 288   for(Int_t i = 0; i < ntracks; i++){
 289
 290     AliVTrack * track =  (AliVTrack*) (GetCTSTracks()->At(i));
 291
 292      //Fill AODParticle after some selection
 293     Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 294     p3.SetXYZ(mom[0],mom[1],mom[2]);
 295
 296     //TOF
 297     ULong_t status = track->GetStatus();
 298     Bool_t okTOF = ( (status & AliVTrack::kTOFout) == AliVTrack::kTOFout ) && ( (status & AliVTrack::kTIME) == AliVTrack::kTIME );
 299     Double32_t tof = track->GetTOFsignal()*1e-3;
 300     //if( tof < 0) printf("TOF Signal %e, status %d, pt %f\n", tof,status,status2,p3.Pt());
 301
 302     if(okTOF)
 303     {
 304       fhTOFSignal  ->Fill(tof);
 305       fhPtTOFSignal->Fill(p3.Pt(), tof);
 306
 307       if(GetReader()->IsPileUpFromSPD())               fhPtTOFSignalPileUp[0]->Fill(p3.Pt(), tof);
 308       if(GetReader()->IsPileUpFromEMCal())             fhPtTOFSignalPileUp[1]->Fill(p3.Pt(), tof);
 309       if(GetReader()->IsPileUpFromSPDOrEMCal())        fhPtTOFSignalPileUp[2]->Fill(p3.Pt(), tof);
 310       if(GetReader()->IsPileUpFromSPDAndEMCal())       fhPtTOFSignalPileUp[3]->Fill(p3.Pt(), tof);
 311       if(GetReader()->IsPileUpFromSPDAndNotEMCal())    fhPtTOFSignalPileUp[4]->Fill(p3.Pt(), tof);
 312       if(GetReader()->IsPileUpFromEMCalAndNotSPD())    fhPtTOFSignalPileUp[5]->Fill(p3.Pt(), tof);
 313       if(GetReader()->IsPileUpFromNotSPDAndNotEMCal()) fhPtTOFSignalPileUp[6]->Fill(p3.Pt(), tof);
 314     }
 315
 316     Bool_t in = GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
 317
 318     if(GetDebug() > 1)
 319       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);
 320
 321     //Acceptance selection
 322     if(IsFiducialCutOn() && ! in ) continue ;
 323
 324     // Momentum selection
 325     if(track->Pt() < GetMinPt() || track->Pt() > GetMaxPt()) continue;
 326
 327     if(okTOF) fhTOFSignalPtCut->Fill(tof);
 328     else
 329     {
 330       fhPtTOFStatus0    ->Fill(track->Pt());
 331       fhEtaPhiTOFStatus0->Fill(track->Eta(),track->Phi());
 332     }
 333
 334     //Keep only particles identified with fPdg
 335     //Selection not done for the moment
 336     //Should be done here.
 337
 338     // Mixed event
 339     if (GetMixedEvent())
 340     {
 341       evtIndex = GetMixedEvent()->EventIndex(track->GetID()) ;
 342     }
 343
 344     GetVertex(vert,evtIndex);
 345     if(TMath::Abs(vert[2])> GetZvertexCut()) return;
 346
 347     AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
 348     tr.SetDetector("CTS");
 349     tr.SetLabel(track->GetLabel());
 350     tr.SetTrackLabel(track->GetID(),-1);
 351     tr.SetChargedBit(track->Charge()>0);
 352
 353     AddAODParticle(tr);
 354
 355   }//loop
 356
 357   if(GetDebug() > 0)
 358     printf("AliAnaChargedParticles::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
 359 }
 360
 361 //__________________________________________________________________
 362 void  AliAnaChargedParticles::MakeAnalysisFillHistograms()
 363 {
 364   //Do analysis and fill histograms
 365
 366   //Loop on stored AODParticles
 367   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 368
 369   fhNtracks->Fill(GetReader()->GetTrackMultiplicity()) ;
 370
 371   if(GetDebug() > 0)
 372     printf("AliAnaChargedParticles::MakeAnalysisFillHistograms() - aod branch entries %d\n", naod);
 373
 374   for(Int_t iaod = 0; iaod < naod ; iaod++){
 375     AliAODPWG4Particle* tr =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 376
 377     fhPt->Fill(tr->Pt());
 378
 379     if(tr->GetChargedBit()){
 380       fhPhiPos   ->Fill(tr->Pt(), tr->Phi());
 381       fhEtaPos   ->Fill(tr->Pt(), tr->Eta());
 382       fhEtaPhiPos->Fill(tr->Eta(),tr->Phi());
 383     }
 384     else{
 385       fhPhiNeg   ->Fill(tr->Pt(), tr->Phi());
 386       fhEtaNeg   ->Fill(tr->Pt(), tr->Eta());
 387       fhEtaPhiNeg->Fill(tr->Eta(),tr->Phi());
 388     }
 389
 390     if(IsDataMC())
 391     {
 392       //Play with the MC stack if available
 393       Int_t mompdg = -1;
 394       Int_t label  = tr->GetLabel();
 395       if(label >= 0)
 396       {
 397         if( GetReader()->ReadStack() && label < GetMCStack()->GetNtrack())
 398         {
 399           TParticle * mom = GetMCStack()->Particle(label);
 400           mompdg =TMath::Abs(mom->GetPdgCode());
 401         }
 402         else if(GetReader()->ReadAODMCParticles())
 403         {
 404           AliAODMCParticle * aodmom = 0;
 405           //Get the list of MC particles
 406           aodmom = (AliAODMCParticle*) (GetReader()->GetAODMCParticles(tr->GetInputFileIndex()))->At(label);
 407           mompdg =TMath::Abs(aodmom->GetPdgCode());
 408         }
 409       }
 410
 411       if(mompdg==211)
 412       {
 413         fhPtPion ->Fill(tr->Pt());
 414         fhPhiPion->Fill(tr->Pt(), tr->Phi());
 415         fhEtaPion->Fill(tr->Pt(), tr->Eta());
 416       }
 417       else if(mompdg==2212)
 418       {
 419         fhPtProton ->Fill(tr->Pt());
 420         fhPhiProton->Fill(tr->Pt(), tr->Phi());
 421         fhEtaProton->Fill(tr->Pt(), tr->Eta());
 422       }
 423       else if(mompdg==321)
 424       {
 425         fhPtKaon ->Fill(tr->Pt());
 426         fhPhiKaon->Fill(tr->Pt(), tr->Phi());
 427         fhEtaKaon->Fill(tr->Pt(), tr->Eta());
 428       }
 429       else if(mompdg==11)
 430       {
 431         fhPtElectron ->Fill(tr->Pt());
 432         fhPhiElectron->Fill(tr->Pt(), tr->Phi());
 433         fhEtaElectron->Fill(tr->Pt(), tr->Eta());
 434       }
 435       else {
 436         //printf("unknown pdg %d\n",mompdg);
 437         fhPtUnknown ->Fill(tr->Pt());
 438         fhPhiUnknown->Fill(tr->Pt(), tr->Phi());
 439         fhEtaUnknown->Fill(tr->Pt(), tr->Eta());
 440       }
 441
 442     }//Work with stack also
 443
 444   }// aod branch loop
 445
 446 }