PWG4/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 /* History of cvs commits:
  18  *
  19  * $Log$
  20  *
  21  */
  22
  23 //_________________________________________________________________________
  24 // Example class on how to read AODCaloClusters, ESDCaloCells and AODTracks and how
  25 // fill AODs with PWG4PartCorr analysis frame
  26 // Select the type of detector information that you want to analyze, CTS (tracking), PHOS or EMCAL
  27 // Select the PID custer type of the calorimeters
  28 // Set min momentum of the cluster/tracks
  29 // Fill few histograms
  30 //
  31 //-- Author: Gustavo Conesa (INFN-LNF)
  32 //-- Author: Gustavo Conesa (LNF-INFN)
  33 //_________________________________________________________________________
  34
  35
  36 // --- ROOT system ---
  37 #include "Riostream.h"
  38 #include "TClonesArray.h"
  39 #include "TH2F.h"
  40 #include "TParticle.h"
  41
  42 //---- AliRoot system ----
  43 #include "AliAnaExample.h"
  44 #include "AliCaloTrackReader.h"
  45 #include "AliLog.h"
  46 #include "AliAODParticleCorrelation.h"
  47 #include "AliESDCaloCells.h"
  48
  49 ClassImp(AliAnaExample)
  50
  51 //____________________________________________________________________________
  52   AliAnaExample::AliAnaExample() :
  53     AliAnaBaseClass(),fPdg(0),  fDetector(""), fhPt(0),fhPhi(0),fhEta(0),  fh2Pt(0),fh2Phi(0),fh2Eta(0),
  54     fhNCells(0), fhAmplitude(0)
  55 {
  56   //Default Ctor
  57
  58   //Initialize parameters
  59   InitParameters();
  60 }
  61
  62 //____________________________________________________________________________
  63 AliAnaExample::AliAnaExample(const AliAnaExample & ex) :
  64   AliAnaBaseClass(ex), fPdg(ex.fPdg), fDetector(ex.fDetector), fhPt(ex.fhPt),  fhPhi(ex.fhPhi),fhEta(ex.fhEta),
  65   fh2Pt(ex.fh2Pt),  fh2Phi(ex.fh2Phi),fh2Eta(ex.fh2Eta), fhNCells(ex.fhNCells), fhAmplitude(ex.fhAmplitude)
  66 {
  67   // cpy ctor
  68
  69 }
  70
  71 //_________________________________________________________________________
  72 AliAnaExample & AliAnaExample::operator = (const AliAnaExample & ex)
  73 {
  74   // assignment operator
  75
  76   if(this == &ex)return *this;
  77   ((AliAnaBaseClass *)this)->operator=(ex);
  78
  79   fPdg = ex.fPdg;
  80   fDetector = ex.fDetector;
  81   fhPt = ex.fhPt;
  82   fhPhi = ex.fhPhi;
  83   fhEta = ex.fhEta;
  84   fh2Pt = ex.fh2Pt;
  85   fh2Phi = ex.fh2Phi;
  86   fh2Eta = ex.fh2Eta;
  87   fhNCells = ex.fhNCells;
  88   fhAmplitude = ex.fhAmplitude;
  89
  90   return *this;
  91
  92 }
  93
  94 // //____________________________________________________________________________
  95 // AliAnaExample::~AliAnaExample()
  96 // {
  97 //   // Remove all pointers except analysis output pointers.
  98
  99 //   ;
 100 // }
 101
 102
 103 //________________________________________________________________________
 104 TList *  AliAnaExample::GetCreateOutputObjects()
 105 {
 106   // Create histograms to be saved in output file and
 107   // store them in fOutputContainer
 108
 109   AliDebug(1,"Init parton histograms");
 110
 111   TList * outputContainer = new TList() ;
 112   outputContainer->SetName("ExampleHistos") ;
 113
 114   fhPt  = new TH1F ("hPt","p_T distribution", 100,0,100);
 115   fhPt->SetXTitle("p_{T} (GeV/c)");
 116   outputContainer->Add(fhPt);
 117
 118   fhPhi  = new TH1F ("hPhi","#phi distribution", 100,0,TMath::TwoPi());
 119   fhPhi->SetXTitle("#phi (rad)");
 120   outputContainer->Add(fhPhi);
 121
 122   fhEta  = new TH1F ("hEta","#eta distribution", 100,-2,2);
 123   fhEta->SetXTitle("#eta ");
 124   outputContainer->Add(fhEta);
 125
 126   //Calo cells
 127   fhNCells  = new TH1F ("hNCells","# cells per event", 100,0,1000);
 128   fhNCells->SetXTitle("n cells");
 129   outputContainer->Add(fhNCells);
 130
 131   fhAmplitude  = new TH1F ("hAmplitude","#eta distribution", 100,0,1000);
 132   fhAmplitude->SetXTitle("Amplitude ");
 133   outputContainer->Add(fhAmplitude);
 134
 135  if(IsDataMC()){
 136    fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", 100,0,100,100,0,100);
 137    fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
 138    fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
 139    outputContainer->Add(fh2Pt);
 140
 141    fh2Phi  = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", 100,0,TMath::TwoPi(), 100,0,TMath::TwoPi());
 142    fh2Phi->SetXTitle("#phi_{rec} (rad)");
 143    fh2Phi->SetYTitle("#phi_{gen} (rad)");
 144    outputContainer->Add(fh2Phi);
 145
 146    fh2Eta  = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", 100,-1,1,100,-1,1);
 147    fh2Eta->SetXTitle("#eta_{rec} ");
 148    fh2Eta->SetYTitle("#eta_{gen} ");
 149    outputContainer->Add(fh2Eta);
 150  }
 151  return outputContainer;
 152 }
 153
 154  //__________________________________________________
 155 void AliAnaExample::InitParameters()
 156 {
 157   //Initialize the parameters of the analysis.
 158
 159   fPdg = 22; //Keep photons
 160   fDetector = "PHOS";
 161
 162 }
 163
 164 //__________________________________________________________________
 165 void AliAnaExample::Print(const Option_t * opt) const
 166 {
 167   //Print some relevant parameters set for the analysis
 168   if(! opt)
 169     return;
 170
 171   printf("Min Pt = %3.2f\n", GetMinPt());
 172   printf("Max Pt = %3.2f\n", GetMaxPt());
 173   printf("Select clusters with pdg %d \n",fPdg);
 174   printf("Select detector %s \n",fDetector.Data());
 175
 176 }
 177
 178 //__________________________________________________________________
 179 void  AliAnaExample::MakeAnalysisFillAOD()
 180 {
 181   //Do analysis and fill aods
 182
 183   //Some prints
 184   if(GetDebug() > 0){
 185     if(fDetector == "EMCAL" && GetAODEMCAL())printf("Example : in emcal aod entries %d\n", GetAODEMCAL()->GetEntries());
 186     if(fDetector == "CTS" && GetAODCTS())printf("Example : in CTS aod entries %d\n", GetAODCTS()->GetEntries());
 187     if(fDetector == "PHOS" && GetAODPHOS())printf("Example : in PHOS aod entries %d\n", GetAODPHOS()->GetEntries());
 188   }
 189
 190   //Get List with tracks or clusters
 191   TClonesArray * partList = new TClonesArray;
 192   if(fDetector == "CTS") partList = GetAODCTS();
 193   else if(fDetector == "EMCAL") partList = GetAODEMCAL();
 194   else if(fDetector == "PHOS") partList = GetAODPHOS();
 195
 196   if(!partList || partList->GetEntries() == 0) return ;
 197
 198   //Fill AODCaloClusters and AODParticleCorrelation with PHOS/EMCAL aods
 199   if(fDetector == "EMCAL" || fDetector == "PHOS"){
 200
 201     //WORK WITH CALOCLUSTERS
 202     ConnectAODCaloClusters(); //Do Only when filling AODCaloClusters
 203     if(GetDebug() > 0) printf("Example: in calo clusters aod entries %d\n", GetAODCaloClusters()->GetEntries());
 204
 205     //Get vertex for photon momentum calculation
 206     Double_t v[3] ; //vertex ;
 207     GetReader()->GetVertex(v);
 208
 209     TLorentzVector mom ;
 210     for(Int_t i = 0; i < partList->GetEntries(); i++){
 211
 212       AliAODCaloCluster * calo =  dynamic_cast<AliAODCaloCluster*> (partList->At(i));
 213
 214       //Fill AODCaloClusters
 215       AddAODCaloCluster(AliAODCaloCluster(*(calo)));
 216
 217       //Fill AODParticleCorrelation after some selection
 218       calo->GetMomentum(mom,v);
 219       Int_t pdg = 0;
 220       if(IsCaloPIDOn()){
 221         Double_t pid[13];
 222         calo->GetPID(pid);
 223         pdg = GetCaloPID()->GetPdg(fDetector,pid,mom.E());
 224         //pdg = GetCaloPID()->GetPdg(fDetector,mom,
 225         //                calo->GetM02(), calo->GetM02(),
 226         //                calo->GetDispersion(), 0, 0);
 227       }
 228
 229       //Acceptance selection
 230       Bool_t in = kTRUE;
 231       if(IsFidutialCutOn())
 232         in =  GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
 233
 234       if(GetDebug() > 1) printf("cluster pt %2.2f, phi %2.2f, pdg %d, in fidutial cut %d \n",mom.Pt(), mom.Phi(), pdg, in);
 235
 236       //Select cluster if momentum, pdg and acceptance are good
 237       if(mom.Pt() > GetMinPt() && pdg ==fPdg && in) {
 238         AliAODParticleCorrelation ph = AliAODParticleCorrelation(mom);
 239         //AddAODParticleCorrelation(AliAODParticleCorrelation(mom));
 240         ph.SetLabel(calo->GetLabel(0));
 241         ph.SetPdg(pdg);
 242         ph.SetDetector(fDetector);
 243         AddAODParticleCorrelation(ph);
 244       }//selection
 245     }//loop
 246
 247     //WORK WITH ESDCALOCELLS
 248     //Don't connect in the same analysis PHOS and EMCAL cells.
 249     AliESDCaloCells * esdCell = new AliESDCaloCells ;
 250     if(fDetector == "PHOS") {
 251       ConnectAODPHOSCells(); //Do Only when filling AODCaloCells
 252       esdCell = (AliESDCaloCells *) GetPHOSCells();
 253     }
 254     else  {
 255       ConnectAODEMCALCells(); //Do Only when filling AODCaloCells
 256       esdCell = (AliESDCaloCells *) GetEMCALCells();
 257     }
 258     //Some prints
 259     if(GetDebug() > 0 && esdCell )
 260       printf("Example : in ESD %s cell entries %d\n", fDetector.Data(), esdCell->GetNumberOfCells());
 261
 262     //Fill AODCells in file
 263     Int_t ncells = esdCell->GetNumberOfCells() ;
 264     GetAODCaloCells()->CreateContainer(ncells);
 265
 266     GetAODCaloCells()->SetType((AliAODCaloCells::AODCells_t) esdCell->GetType());
 267
 268     for (Int_t iCell = 0; iCell < ncells; iCell++) {
 269       if(GetDebug() > 2)  printf("cell : amp %f, absId %d,  time %f\n", esdCell->GetAmplitude(iCell), esdCell->GetCellNumber(iCell), esdCell->GetTime(iCell));
 270
 271       GetAODCaloCells()->SetCell(iCell,esdCell->GetCellNumber(iCell),esdCell->GetAmplitude(iCell));
 272     }
 273     GetAODCaloCells()->Sort();
 274
 275   }//cluster-cell analysis
 276   else if(fDetector == "CTS"){ //Track analysis
 277     //Fill AODParticleCorrelation with CTS aods
 278     TVector3 p3;
 279     for(Int_t i = 0; i < GetAODCTS()->GetEntries(); i++){
 280
 281       AliAODTrack * track =  dynamic_cast<AliAODTrack*> (GetAODCTS()->At(i));
 282
 283       //Fill AODParticleCorrelation after some selection
 284       Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 285       p3.SetXYZ(mom[0],mom[1],mom[2]);
 286
 287       //Acceptance selection
 288       Bool_t in =  GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
 289       if(GetDebug() > 1) printf("track pt %2.2f, phi %2.2f, in fidutial cut %d\n",p3.Pt(), p3.Phi(), in);
 290       if(p3.Pt() > GetMinPt() && in) {
 291         AliAODParticleCorrelation tr = AliAODParticleCorrelation(mom[0],mom[1],mom[2],0);
 292         //AddAODParticleCorrelation(AliAODParticleCorrelation(mom));
 293         tr.SetDetector("CTS");
 294         AddAODParticleCorrelation(tr);
 295
 296       }//selection
 297     }//loop
 298   }//CTS analysis
 299
 300   if(GetDebug() > 0) {
 301     printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntries());
 302     printf("Example: final aod branch entries %d\n", GetAODBranch()->GetEntries());
 303     printf("Example: final aod cell  entries %d\n", GetAODCaloCells()->GetNumberOfCells());
 304   }
 305 }
 306
 307 //__________________________________________________________________
 308 void  AliAnaExample::MakeAnalysisFillHistograms()
 309 {
 310   //Do analysis and fill histograms
 311
 312   //Loop on stored AODParticles
 313   Int_t naod = GetAODBranch()->GetEntriesFast();
 314   if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
 315   for(Int_t iaod = 0; iaod < naod ; iaod++){
 316    AliAODParticleCorrelation* ph =  dynamic_cast<AliAODParticleCorrelation*> (GetAODBranch()->At(iaod));
 317
 318    fhPt->Fill(ph->Pt());
 319    fhPhi->Fill(ph->Phi());
 320    fhEta->Fill(ph->Eta());
 321
 322    if(IsDataMC()){
 323      //Play with the MC stack if available
 324      AliStack * stack =  GetMCStack() ;
 325
 326      if(ph->GetLabel() < 0 || !stack) {
 327        printf("*** bad label or no stack ***:  label %d \n", ph->GetLabel());
 328        continue;
 329      }
 330
 331      if(ph->GetLabel() >=  stack->GetNtrack()) {
 332        printf("*** large label ***:  label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
 333        continue ;
 334      }
 335
 336      TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
 337
 338      fh2Pt->Fill(ph->Pt(), mom->Pt());
 339      fh2Phi->Fill(ph->Phi(), mom->Phi());
 340      fh2Eta->Fill(ph->Eta(), mom->Eta());
 341    }//Work with stack also
 342   }// aod branch loop
 343
 344   // CaloCells histograms
 345   if(GetAODCaloCells()){
 346
 347     Int_t ncells = GetAODCaloCells()->GetNumberOfCells();
 348     fhNCells->Fill(ncells) ;
 349
 350     for (Int_t iCell = 0; iCell < ncells; iCell++) {
 351       if(GetDebug() > 2)  printf("cell : amp %f, absId %d \n", GetAODCaloCells()->GetAmplitude(iCell), GetAODCaloCells()->GetCellNumber(iCell));
 352       fhAmplitude->Fill(GetAODCaloCells()->GetAmplitude(iCell));
 353     }
 354   }//calo cells container exist
 355 }