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