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 AODCaloClusters, ESDCaloCells 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 "AliESDCaloCells.h"
  43 #include "AliStack.h"
  44 #include "AliCaloPID.h"
  45 #include "AliFiducialCut.h"
  46 #include "AliAODCaloCells.h"
  47 #include "AliAODCaloCluster.h"
  48 #include "AliAODTrack.h"
  49
  50 ClassImp(AliAnaExample)
  51
  52 //____________________________________________________________________________
  53   AliAnaExample::AliAnaExample() :
  54     AliAnaPartCorrBaseClass(),fPdg(0),  fDetector(""), fhPt(0),fhPhi(0),fhEta(0),  fh2Pt(0),fh2Phi(0),fh2Eta(0),
  55     fhNCells(0), fhAmplitude(0)
  56 {
  57   //Default Ctor
  58
  59   //Initialize parameters
  60   InitParameters();
  61 }
  62
  63 //____________________________________________________________________________
  64 AliAnaExample::AliAnaExample(const AliAnaExample & ex) :
  65   AliAnaPartCorrBaseClass(ex), fPdg(ex.fPdg), fDetector(ex.fDetector), fhPt(ex.fhPt),  fhPhi(ex.fhPhi),fhEta(ex.fhEta),
  66   fh2Pt(ex.fh2Pt),  fh2Phi(ex.fh2Phi),fh2Eta(ex.fh2Eta), fhNCells(ex.fhNCells), fhAmplitude(ex.fhAmplitude)
  67 {
  68   // cpy ctor
  69
  70 }
  71
  72 //_________________________________________________________________________
  73 AliAnaExample & AliAnaExample::operator = (const AliAnaExample & ex)
  74 {
  75   // assignment operator
  76
  77   if(this == &ex)return *this;
  78   ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
  79
  80   fPdg = ex.fPdg;
  81   fDetector = ex.fDetector;
  82   fhPt = ex.fhPt;
  83   fhPhi = ex.fhPhi;
  84   fhEta = ex.fhEta;
  85   fh2Pt = ex.fh2Pt;
  86   fh2Phi = ex.fh2Phi;
  87   fh2Eta = ex.fh2Eta;
  88   fhNCells = ex.fhNCells;
  89   fhAmplitude = ex.fhAmplitude;
  90
  91   return *this;
  92
  93 }
  94
  95 // //____________________________________________________________________________
  96 // AliAnaExample::~AliAnaExample()
  97 // {
  98 //   // Remove all pointers except analysis output pointers.
  99
 100 //   ;
 101 // }
 102
 103
 104 //________________________________________________________________________
 105 TList *  AliAnaExample::GetCreateOutputObjects()
 106 {
 107   // Create histograms to be saved in output file and
 108   // store them in fOutputContainer
 109
 110   TList * outputContainer = new TList() ;
 111   outputContainer->SetName("ExampleHistos") ;
 112
 113   Int_t nptbins  = GetHistoNPtBins();
 114   Int_t nphibins = GetHistoNPhiBins();
 115   Int_t netabins = GetHistoNEtaBins();
 116   Float_t ptmax  = GetHistoPtMax();
 117   Float_t phimax = GetHistoPhiMax();
 118   Float_t etamax = GetHistoEtaMax();
 119   Float_t ptmin  = GetHistoPtMin();
 120   Float_t phimin = GetHistoPhiMin();
 121   Float_t etamin = GetHistoEtaMin();
 122
 123   fhPt  = new TH1F ("hPt","p_T distribution", nptbins,ptmin,ptmax);
 124   fhPt->SetXTitle("p_{T} (GeV/c)");
 125   outputContainer->Add(fhPt);
 126
 127   fhPhi  = new TH1F ("hPhi","#phi distribution",nphibins,phimin,phimax);
 128   fhPhi->SetXTitle("#phi (rad)");
 129   outputContainer->Add(fhPhi);
 130
 131   fhEta  = new TH1F ("hEta","#eta distribution",netabins,etamin,etamax);
 132   fhEta->SetXTitle("#eta ");
 133   outputContainer->Add(fhEta);
 134
 135   if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
 136
 137     //Calo cells
 138     fhNCells  = new TH1F ("hNCells","# cells per event", 100,0,1000);
 139     fhNCells->SetXTitle("n cells");
 140     outputContainer->Add(fhNCells);
 141
 142     fhAmplitude  = new TH1F ("hAmplitude","#eta distribution", 100,0,1000);
 143     fhAmplitude->SetXTitle("Amplitude ");
 144     outputContainer->Add(fhAmplitude);
 145   }
 146
 147   if(IsDataMC()){
 148     fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 149     fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
 150     fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
 151     outputContainer->Add(fh2Pt);
 152
 153     fh2Phi  = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", nphibins,phimin,phimax, nphibins,phimin,phimax);
 154     fh2Phi->SetXTitle("#phi_{rec} (rad)");
 155     fh2Phi->SetYTitle("#phi_{gen} (rad)");
 156     outputContainer->Add(fh2Phi);
 157
 158     fh2Eta  = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", netabins,etamin,etamax,netabins,etamin,etamax);
 159     fh2Eta->SetXTitle("#eta_{rec} ");
 160     fh2Eta->SetYTitle("#eta_{gen} ");
 161     outputContainer->Add(fh2Eta);
 162
 163   }
 164
 165   return outputContainer;
 166 }
 167
 168 //__________________________________________________
 169 void AliAnaExample::InitParameters()
 170 {
 171   //Initialize the parameters of the analysis.
 172   SetOutputAODClassName("AliAODPWG4Particle");
 173   SetOutputAODName("PWG4Particle");
 174   AddToHistogramsName("AnaExample_");
 175
 176   fPdg = 22; //Keep photons
 177   fDetector = "PHOS";
 178
 179 }
 180
 181 //__________________________________________________________________
 182 void AliAnaExample::Print(const Option_t * opt) const
 183 {
 184   //Print some relevant parameters set for the analysis
 185   if(! opt)
 186     return;
 187
 188   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 189   AliAnaPartCorrBaseClass::Print(" ");
 190
 191   printf("Select clusters with pdg %d \n",fPdg);
 192   printf("Select detector %s \n",fDetector.Data());
 193
 194 }
 195
 196 //__________________________________________________________________
 197 void  AliAnaExample::MakeAnalysisFillAOD()
 198 {
 199   //Do analysis and fill aods
 200
 201   //Some prints
 202   if(GetDebug() > 0){
 203     if(fDetector == "EMCAL" && GetAODEMCAL())printf("AliAnaExample::MakeAnalysisFillAOD() - In EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
 204     if(fDetector == "CTS" && GetAODCTS())printf("AliAnaExample::MakeAnalysisFillAOD() - In CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
 205     if(fDetector == "PHOS" && GetAODPHOS())printf("AliAnaExample::MakeAnalysisFillAOD() - In PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
 206   }
 207
 208   //Get List with tracks or clusters
 209   TObjArray * partList = new TObjArray();
 210   if(fDetector == "CTS") partList = GetAODCTS();
 211   else if(fDetector == "EMCAL") partList = GetAODEMCAL();
 212   else if(fDetector == "PHOS") partList = GetAODPHOS();
 213
 214   if(!partList || partList->GetEntriesFast() == 0) return ;
 215
 216   //Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods
 217   if(fDetector == "EMCAL" || fDetector == "PHOS"){
 218
 219     //Get vertex for photon momentum calculation
 220     Double_t v[3] ; //vertex ;
 221     GetReader()->GetVertex(v);
 222
 223     TLorentzVector mom ;
 224     for(Int_t i = 0; i < partList->GetEntriesFast(); i++){
 225
 226       AliAODCaloCluster * calo =  (AliAODCaloCluster*) (partList->At(i));
 227
 228       //Fill AODParticle after some selection
 229       calo->GetMomentum(mom,v);
 230       Int_t pdg = fPdg;
 231
 232       if(IsCaloPIDOn()){
 233         Double_t pid[13];
 234         calo->GetPID(pid);
 235         pdg = GetCaloPID()->GetPdg(fDetector,pid,mom.E());
 236         //pdg = GetCaloPID()->GetPdg(fDetector,mom,
 237         //                calo->GetM02(), calo->GetM02(),
 238         //                calo->GetDispersion(), 0, 0);
 239       }
 240
 241       //Acceptance selection
 242       Bool_t in = kTRUE;
 243       if(IsFiducialCutOn())
 244         in =  GetFiducialCut()->IsInFiducialCut(mom,fDetector) ;
 245
 246       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);
 247
 248       //Select cluster if momentum, pdg and acceptance are good
 249       if(mom.Pt() > GetMinPt() && pdg ==fPdg && in) {
 250         AliAODPWG4Particle ph = AliAODPWG4Particle(mom);
 251         //AddAODParticleCorrelation(AliAODPWG4ParticleCorrelation(mom));
 252         ph.SetLabel(calo->GetLabel(0));
 253         ph.SetPdg(pdg);
 254         ph.SetDetector(fDetector);
 255         AddAODParticle(ph);
 256       }//selection
 257     }//loop
 258
 259     if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
 260       //WORK WITH ESDCALOCELLS
 261       //Don't connect in the same analysis PHOS and EMCAL cells.
 262
 263       AliESDCaloCells * esdCell = new AliESDCaloCells ;
 264       if(fDetector == "PHOS") {
 265         ConnectAODPHOSCells(); //Do Only when filling AODCaloCells
 266         esdCell = (AliESDCaloCells *) GetPHOSCells();
 267       }
 268       else  {
 269         ConnectAODEMCALCells(); //Do Only when filling AODCaloCells
 270         esdCell = (AliESDCaloCells *) GetEMCALCells();
 271       }
 272
 273       if(!esdCell) {
 274         printf("AliAnaExample::MakeAnalysisFillAOD() - STOP: No CELLS available for analysis");
 275         abort();
 276       }
 277       //Some prints
 278       if(GetDebug() > 0 && esdCell )
 279         printf("AliAnaExample::MakeAnalysisFillAOD() - In ESD %s cell entries %d\n", fDetector.Data(), esdCell->GetNumberOfCells());
 280
 281       //Fill AODCells in file
 282       Int_t ncells = esdCell->GetNumberOfCells() ;
 283       GetAODCaloCells()->CreateContainer(ncells);
 284
 285       GetAODCaloCells()->SetType((AliAODCaloCells::AODCells_t) esdCell->GetType());
 286
 287       for (Int_t iCell = 0; iCell < ncells; iCell++) {
 288         if(GetDebug() > 2)  printf("AliAnaExample::MakeAnalysisFillAOD() - Cell : amp %f, absId %d,  time %f\n", esdCell->GetAmplitude(iCell), esdCell->GetCellNumber(iCell), esdCell->GetTime(iCell));
 289
 290         GetAODCaloCells()->SetCell(iCell,esdCell->GetCellNumber(iCell),esdCell->GetAmplitude(iCell));
 291       }
 292       GetAODCaloCells()->Sort();
 293     }
 294   }//cluster-cell analysis
 295   else if(fDetector == "CTS"){ //Track analysis
 296     //Fill AODParticle with CTS aods
 297     TVector3 p3;
 298     for(Int_t i = 0; i < GetAODCTS()->GetEntriesFast(); i++){
 299
 300       AliAODTrack * track =  (AliAODTrack*) (GetAODCTS()->At(i));
 301
 302       //Fill AODParticle after some selection
 303       Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 304       p3.SetXYZ(mom[0],mom[1],mom[2]);
 305
 306       //Acceptance selection
 307       Bool_t in =  GetFiducialCut()->IsInFiducialCut(mom,"CTS") ;
 308       if(GetDebug() > 1) printf("AliAnaExample::MakeAnalysisFillAOD() - Track pt %2.2f, phi %2.2f, in fiducial cut %d\n",p3.Pt(), p3.Phi(), in);
 309       if(p3.Pt() > GetMinPt() && in) {
 310         AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
 311         tr.SetDetector("CTS");
 312         AddAODParticle(tr);
 313       }//selection
 314     }//loop
 315   }//CTS analysis
 316
 317   if(GetDebug() > 0) {
 318     if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
 319       //printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntriesFast());
 320       printf("AliAnaExample::MakeAnalysisFillAOD() - Final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
 321     //  if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
 322     //printf("Example: final aod cell  entries %d\n", GetAODCaloCells()->GetNumberOfCells());
 323   }
 324 }
 325
 326 //__________________________________________________________________
 327 void  AliAnaExample::MakeAnalysisFillHistograms()
 328 {
 329   //Do analysis and fill histograms
 330
 331   //Loop on stored AODParticles
 332   Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 333   if(GetDebug() > 0) printf("AliAnaExample::MakeAnalysisFillHistograms() - Histo aod branch entries %d\n", naod);
 334   for(Int_t iaod = 0; iaod < naod ; iaod++){
 335     AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 336
 337     fhPt->Fill(ph->Pt());
 338     fhPhi->Fill(ph->Phi());
 339     fhEta->Fill(ph->Eta());
 340
 341     if(IsDataMC()){
 342       //Play with the MC stack if available
 343           Float_t ptprim  = 0;
 344           Float_t phiprim = 0;
 345           Float_t etaprim = 0;
 346           if(GetReader()->ReadStack()){
 347                   AliStack * stack =  GetMCStack() ;
 348
 349                   if(ph->GetLabel() < 0 || !stack) {
 350                           printf("AliAnaExample::MakeAnalysisFillHistograms() *** bad label or no stack ***:  label %d \n", ph->GetLabel());
 351                           continue;
 352                   }
 353
 354                   if(ph->GetLabel() >=  stack->GetNtrack()) {
 355                           printf("AliAnaExample::MakeAnalysisFillHistograms() *** large label ***:  label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
 356                           continue ;
 357                   }
 358
 359                   TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
 360           ptprim  = mom->Pt();
 361                   phiprim = mom->Phi();
 362                   etaprim = mom->Eta();
 363           }
 364           else if(GetReader()->ReadAODMCParticles()){
 365                   printf("AliAnaExample::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
 366           }
 367           fh2Pt->Fill(ph->Pt(), ptprim);
 368           fh2Phi->Fill(ph->Phi(), phiprim);
 369           fh2Eta->Fill(ph->Eta(), etaprim);
 370     }//Work with stack also
 371   }// aod branch loop
 372
 373   // CaloCells histograms
 374   if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
 375     if(GetAODCaloCells()){
 376
 377       Int_t ncells = GetAODCaloCells()->GetNumberOfCells();
 378       fhNCells->Fill(ncells) ;
 379
 380       for (Int_t iCell = 0; iCell < ncells; iCell++) {
 381         if(GetDebug() > 2)  printf("AliAnaExample::MakeAnalysisFillHistograms() - Cell : amp %f, absId %d \n", GetAODCaloCells()->GetAmplitude(iCell), GetAODCaloCells()->GetCellNumber(iCell));
 382         fhAmplitude->Fill(GetAODCaloCells()->GetAmplitude(iCell));
 383       }
 384     }//calo cells container exist
 385   }
 386 }
 387
 388 //________________________________________________________________________
 389 void AliAnaExample::ReadHistograms(TList* outputList)
 390 {
 391         // Needed when Terminate is executed in distributed environment
 392         // Refill analysis histograms of this class with corresponding histograms in output list.
 393
 394         // Histograms of this analsys are kept in the same list as other analysis, recover the position of
 395         // the first one and then add the next
 396         Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hPt"));
 397
 398         //Read histograms, must be in the same order as in GetCreateOutputObject.
 399         fhPt   = (TH1F *) outputList->At(index++);
 400         fhPhi  = (TH1F *) outputList->At(index++);
 401         fhEta  = (TH1F *) outputList->At(index++);
 402
 403         if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
 404          fhNCells     = (TH1F *) outputList->At(index++);
 405          fhAmplitude  = (TH1F *) outputList->At(index++);
 406         }
 407
 408         if(IsDataMC()){
 409           fh2Pt  = (TH2F *) outputList->At(index++);
 410           fh2Phi = (TH2F *) outputList->At(index++);
 411           fh2Eta = (TH2F *) outputList->At(index);
 412         }
 413
 414 }
 415
 416 //__________________________________________________________________
 417 void  AliAnaExample::Terminate(TList* outputList)
 418 {
 419
 420   //Do some plots to end
 421
 422   //Recover histograms from output histograms list, needed for distributed analysis.
 423   ReadHistograms(outputList);
 424
 425   printf(" AliAnaExample::Terminate()  *** %s Report:", GetName()) ;
 426   printf(" AliAnaExample::Terminate()        pt         : %5.3f , RMS : %5.3f \n", fhPt->GetMean(),   fhPt->GetRMS() ) ;
 427
 428   TCanvas  * c = new TCanvas("c", "PHOS ESD Test", 400, 10, 600, 700) ;
 429   c->Divide(1, 3);
 430
 431   c->cd(1) ;
 432   gPad->SetLogy();
 433   fhPt->SetLineColor(2);
 434   fhPt->Draw();
 435
 436   c->cd(2) ;
 437   fhPhi->SetLineColor(2);
 438   fhPhi->Draw();
 439
 440   c->cd(3) ;
 441   fhEta->SetLineColor(2);
 442   fhEta->Draw();
 443
 444   c->Print("Example.eps");
 445
 446   char line[1024] ;
 447   sprintf(line, ".!tar -zcf %s.tar.gz *.eps", GetName()) ;
 448   gROOT->ProcessLine(line);
 449   sprintf(line, ".!rm -fR *.eps");
 450   gROOT->ProcessLine(line);
 451
 452   printf("AliAnaExample::Terminate() - !! All the eps files are in %s.tar.gz !!!", GetName());
 453
 454 }