PWG4/PartCorrDep/AliAnaExample.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 // Example class on how to read CaloClusters, and AODTracks and how
  19 // fill AODs with PWG4PartCorr analysis frame
  20 // Select the type of detector information that you want to analyze, CTS (tracking), PHOS or EMCAL
  21 // Select the PID custer type of the calorimeters
  22 // Set min momentum of the cluster/tracks
  23 // Fill few histograms
  24 //
  25 //-- Author: Gustavo Conesa (INFN-LNF)
  26 //_________________________________________________________________________
  27
  28
  29 // --- ROOT system ---
  30 //#include "Riostream.h"
  31 #include "TObjArray.h"
  32 #include "TParticle.h"
  33 #include "TCanvas.h"
  34 #include "TPad.h"
  35 #include "TROOT.h"
  36 #include "TH2F.h"
  37
  38 //---- AliRoot system ----
  39 #include "AliAnaExample.h"
  40 #include "AliCaloTrackReader.h"
  41 #include "AliAODPWG4Particle.h"
  42 #include "AliStack.h"
  43 #include "AliCaloPID.h"
  44 #include "AliFiducialCut.h"
  45 #include "AliVCluster.h"
  46 #include "AliVTrack.h"
  47
  48 ClassImp(AliAnaExample)
  49
  50 //____________________________________________________________________________
  51   AliAnaExample::AliAnaExample() :
  52     AliAnaPartCorrBaseClass(),fPdg(0),  fDetector(""), fhPt(0),fhPhi(0),fhEta(0),  fh2Pt(0),fh2Phi(0),fh2Eta(0)
  53 {
  54   //Default Ctor
  55
  56   //Initialize parameters
  57   InitParameters();
  58 }
  59
  60 //________________________________________________________________________
  61 TList *  AliAnaExample::GetCreateOutputObjects()
  62 {
  63   // Create histograms to be saved in output file and
  64   // store them in fOutputContainer
  65
  66   TList * outputContainer = new TList() ;
  67   outputContainer->SetName("ExampleHistos") ;
  68
  69   Int_t nptbins  = GetHistoPtBins();
  70   Int_t nphibins = GetHistoPhiBins();
  71   Int_t netabins = GetHistoEtaBins();
  72   Float_t ptmax  = GetHistoPtMax();
  73   Float_t phimax = GetHistoPhiMax();
  74   Float_t etamax = GetHistoEtaMax();
  75   Float_t ptmin  = GetHistoPtMin();
  76   Float_t phimin = GetHistoPhiMin();
  77   Float_t etamin = GetHistoEtaMin();
  78
  79   fhPt  = new TH1F ("hPt","p_T distribution", nptbins,ptmin,ptmax);
  80   fhPt->SetXTitle("p_{T} (GeV/c)");
  81   outputContainer->Add(fhPt);
  82
  83   fhPhi  = new TH1F ("hPhi","#phi distribution",nphibins,phimin,phimax);
  84   fhPhi->SetXTitle("#phi (rad)");
  85   outputContainer->Add(fhPhi);
  86
  87   fhEta  = new TH1F ("hEta","#eta distribution",netabins,etamin,etamax);
  88   fhEta->SetXTitle("#eta ");
  89   outputContainer->Add(fhEta);
  90
  91   if(IsDataMC()){
  92     fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
  93     fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
  94     fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
  95     outputContainer->Add(fh2Pt);
  96
  97     fh2Phi  = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", nphibins,phimin,phimax, nphibins,phimin,phimax);
  98     fh2Phi->SetXTitle("#phi_{rec} (rad)");
  99     fh2Phi->SetYTitle("#phi_{gen} (rad)");
 100     outputContainer->Add(fh2Phi);
 101
 102     fh2Eta  = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", netabins,etamin,etamax,netabins,etamin,etamax);
 103     fh2Eta->SetXTitle("#eta_{rec} ");
 104     fh2Eta->SetYTitle("#eta_{gen} ");
 105     outputContainer->Add(fh2Eta);
 106
 107   }
 108
 109   return outputContainer;
 110 }
 111
 112 //__________________________________________________
 113 void AliAnaExample::InitParameters()
 114 {
 115   //Initialize the parameters of the analysis.
 116   SetOutputAODClassName("AliAODPWG4Particle");
 117   SetOutputAODName("PWG4Particle");
 118   AddToHistogramsName("AnaExample_");
 119
 120   fPdg = 22; //Keep photons
 121   fDetector = "PHOS";
 122
 123 }
 124
 125 //__________________________________________________________________
 126 void AliAnaExample::Print(const Option_t * opt) const
 127 {
 128   //Print some relevant parameters set for the analysis
 129   if(! opt)
 130     return;
 131
 132   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 133   AliAnaPartCorrBaseClass::Print(" ");
 134
 135   printf("Select clusters with pdg %d \n",fPdg);
 136   printf("Select detector %s \n",fDetector.Data());
 137
 138 }
 139
 140 //__________________________________________________________________
 141 void  AliAnaExample::MakeAnalysisFillAOD()
 142 {
 143   //Do analysis and fill aods
 144
 145   //Some prints
 146   if(GetDebug() > 0){
 147     if(fDetector == "EMCAL" && GetEMCALClusters())printf("AliAnaExample::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetEMCALClusters()->GetEntriesFast());
 148     if(fDetector == "CTS"   && GetCTSTracks())  printf("AliAnaExample::MakeAnalysisFillAOD() - In CTS aod entries %d\n",   GetCTSTracks()  ->GetEntriesFast());
 149     if(fDetector == "PHOS"  && GetPHOSClusters()) printf("AliAnaExample::MakeAnalysisFillAOD() - In PHOS aod entries %d\n",  GetPHOSClusters() ->GetEntriesFast());
 150   }
 151
 152   //Get List with tracks or clusters
 153   TObjArray * partList = 0x0;
 154   if(fDetector == "CTS")        partList = GetCTSTracks();
 155   else if(fDetector == "EMCAL") partList = GetEMCALClusters();
 156   else if(fDetector == "PHOS")  partList = GetPHOSClusters();
 157
 158   if(!partList || partList->GetEntriesFast() == 0) return ;
 159
 160   //Fill AODParticle with PHOS/EMCAL aods
 161   if(fDetector == "EMCAL" || fDetector == "PHOS"){
 162
 163     //Get vertex for photon momentum calculation
 164
 165
 166     TLorentzVector mom ;
 167     for(Int_t i = 0; i < partList->GetEntriesFast(); i++){
 168
 169       AliVCluster * calo =  (AliVCluster*) (partList->At(i));
 170
 171       //Get the index where the cluster comes, to retrieve the corresponding vertex
 172       Int_t evtIndex = 0 ;
 173       if (GetMixedEvent()) {
 174         evtIndex=GetMixedEvent()->EventIndexForCaloCluster(calo->GetID()) ;
 175       }
 176
 177       //Fill AODParticle after some selection
 178       if(GetReader()->GetDataType() != AliCaloTrackReader::kMC){
 179         calo->GetMomentum(mom,GetVertex(evtIndex)) ;}//Assume that come from vertex in straight line
 180       else{
 181         Double_t vertex[]={0,0,0};
 182         calo->GetMomentum(mom,vertex) ;
 183       }
 184
 185       Int_t pdg = fPdg;
 186
 187       if(IsCaloPIDOn()){
 188         const Double_t *pid = calo->GetPID();
 189         pdg = GetCaloPID()->GetPdg(fDetector,pid,mom.E());
 190         //pdg = GetCaloPID()->GetPdg(fDetector,mom,
 191         //                calo->GetM02(), calo->GetM02(),
 192         //                calo->GetDispersion(), 0, 0);
 193       }
 194
 195       //Acceptance selection
 196       Bool_t in = kTRUE;
 197       if(IsFiducialCutOn())
 198         in =  GetFiducialCut()->IsInFiducialCut(mom,fDetector) ;
 199
 200       if(GetDebug() > 1) printf("AliAnaExample::MakeAnalysisFillAOD() - Cluster pt %2.2f, phi %2.2f, pdg %d, in fiducial cut %d \n",mom.Pt(), mom.Phi(), pdg, in);
 201
 202       //Select cluster if momentum, pdg and acceptance are good
 203       if(mom.Pt() > GetMinPt() && pdg ==fPdg && in) {
 204         AliAODPWG4Particle ph = AliAODPWG4Particle(mom);
 205         //AddAODParticleCorrelation(AliAODPWG4ParticleCorrelation(mom));
 206         ph.SetLabel(calo->GetLabel());
 207         ph.SetPdg(pdg);
 208         ph.SetDetector(fDetector);
 209         AddAODParticle(ph);
 210       }//selection
 211     }//loop
 212   }//Calorimeters
 213   else if(fDetector == "CTS"){ //Track analysis
 214     //Fill AODParticle with CTS aods
 215     TVector3 p3;
 216     for(Int_t i = 0; i < GetCTSTracks()->GetEntriesFast(); i++){
 217
 218       AliVTrack * track =  (AliVTrack*) (GetCTSTracks()->At(i));
 219
 220       //Fill AODParticle after some selection
 221       Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 222       p3.SetXYZ(mom[0],mom[1],mom[2]);
 223
 224       //Acceptance selection
 225       Bool_t in =  GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
 226       if(GetDebug() > 1) printf("AliAnaExample::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, in fiducial cut %d\n",p3.Pt(), p3.Phi(), in);
 227       if(p3.Pt() > GetMinPt() && in) {
 228         AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
 229         tr.SetDetector("CTS");
 230         AddAODParticle(tr);
 231       }//selection
 232     }//loop
 233   }//CTS analysis
 234
 235   if(GetDebug() > 0) {
 236     if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
 237       printf("AliAnaExample::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
 238     //  if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
 239     //printf("Example: final aod cell  entries %d\n", GetAODCaloCells()->GetNumberOfCells());
 240   }
 241 }
 242
 243 //__________________________________________________________________
 244 void  AliAnaExample::MakeAnalysisFillHistograms()
 245 {
 246   //Do analysis and fill histograms
 247
 248   //Loop on stored AODParticles
 249   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 250   if(GetDebug() > 0) printf("AliAnaExample::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
 251   for(Int_t iaod = 0; iaod < naod ; iaod++){
 252     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 253
 254     fhPt->Fill(ph->Pt());
 255     fhPhi->Fill(ph->Phi());
 256     fhEta->Fill(ph->Eta());
 257
 258     if(IsDataMC()){
 259     //Play with the MC stack if available
 260           Float_t ptprim  = 0;
 261           Float_t phiprim = 0;
 262           Float_t etaprim = 0;
 263           if(GetReader()->ReadStack()){
 264                   AliStack * stack =  GetMCStack() ;
 265
 266                   if(ph->GetLabel() < 0 || !stack) {
 267                           printf("AliAnaExample::MakeAnalysisFillHistograms() *** bad label or no stack ***:  label %d \n", ph->GetLabel());
 268                           continue;
 269                   }
 270
 271                   if(ph->GetLabel() >=  stack->GetNtrack()) {
 272                           printf("AliAnaExample::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
 273                           continue ;
 274                   }
 275
 276                   TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
 277       ptprim  = mom->Pt();
 278                   phiprim = mom->Phi();
 279                   etaprim = mom->Eta();
 280           }
 281           else if(GetReader()->ReadAODMCParticles()){
 282                   printf("AliAnaExample::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
 283           }
 284           fh2Pt->Fill(ph->Pt(), ptprim);
 285           fh2Phi->Fill(ph->Phi(), phiprim);
 286           fh2Eta->Fill(ph->Eta(), etaprim);
 287     }//Work with stack also
 288   }// aod branch loop
 289
 290  }
 291
 292 //________________________________________________________________________
 293 void AliAnaExample::ReadHistograms(TList* outputList)
 294 {
 295         // Needed when Terminate is executed in distributed environment
 296         // Refill analysis histograms of this class with corresponding histograms in output list.
 297
 298         // Histograms of this analsys are kept in the same list as other analysis, recover the position of
 299         // the first one and then add the next
 300         Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hPt"));
 301
 302         //Read histograms, must be in the same order as in GetCreateOutputObject.
 303         fhPt   = (TH1F *) outputList->At(index++);
 304         fhPhi  = (TH1F *) outputList->At(index++);
 305         fhEta  = (TH1F *) outputList->At(index++);
 306
 307         if(IsDataMC()){
 308           fh2Pt  = (TH2F *) outputList->At(index++);
 309           fh2Phi = (TH2F *) outputList->At(index++);
 310           fh2Eta = (TH2F *) outputList->At(index);
 311         }
 312
 313 }
 314
 315 //__________________________________________________________________
 316 void  AliAnaExample::Terminate(TList* outputList)
 317 {
 318
 319   //Do some plots to end
 320
 321   //Recover histograms from output histograms list, needed for distributed analysis.
 322   ReadHistograms(outputList);
 323
 324   printf(" AliAnaExample::Terminate()  *** %s Report:", GetName()) ;
 325   printf(" AliAnaExample::Terminate()        pt         : %5.3f , RMS : %5.3f \n", fhPt->GetMean(),   fhPt->GetRMS() ) ;
 326
 327   TCanvas  * c = new TCanvas("c", "PHOS ESD Test", 400, 10, 600, 700) ;
 328   c->Divide(1, 3);
 329
 330   c->cd(1) ;
 331   gPad->SetLogy();
 332   fhPt->SetLineColor(2);
 333   fhPt->Draw();
 334
 335   c->cd(2) ;
 336   fhPhi->SetLineColor(2);
 337   fhPhi->Draw();
 338
 339   c->cd(3) ;
 340   fhEta->SetLineColor(2);
 341   fhEta->Draw();
 342
 343   c->Print("Example.eps");
 344
 345   const Int_t buffersize = 1024;
 346   char line[buffersize] ;
 347   snprintf(line,buffersize, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
 348   gROOT->ProcessLine(line);
 349   snprintf(line,buffersize, ".!rm -fR *.eps");
 350   gROOT->ProcessLine(line);
 351
 352   printf("AliAnaExample::Terminate() - !! All the eps files are in %s.tar.gz !!!", GetName());
 353
 354 }