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 "TClonesArray.h"
  32 #include "TParticle.h"
  33
  34 //---- AliRoot system ----
  35 #include "AliAnaExample.h"
  36 #include "AliCaloTrackReader.h"
  37 #include "AliLog.h"
  38 #include "AliAODPWG4Particle.h"
  39 #include "AliESDCaloCells.h"
  40 #include "AliStack.h"
  41 #include "AliCaloPID.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         Int_t nptbins  = GetHistoNPtBins();
 109         Int_t nphibins = GetHistoNPhiBins();
 110         Int_t netabins = GetHistoNEtaBins();
 111         Float_t ptmax  = GetHistoPtMax();
 112         Float_t phimax = GetHistoPhiMax();
 113         Float_t etamax = GetHistoEtaMax();
 114         Float_t ptmin  = GetHistoPtMin();
 115         Float_t phimin = GetHistoPhiMin();
 116         Float_t etamin = GetHistoEtaMin();
 117
 118         fhPt  = new TH1F ("hPt","p_T distribution", nptbins,ptmin,ptmax);
 119         fhPt->SetXTitle("p_{T} (GeV/c)");
 120         outputContainer->Add(fhPt);
 121
 122         fhPhi  = new TH1F ("hPhi","#phi distribution",nphibins,phimin,phimax);
 123         fhPhi->SetXTitle("#phi (rad)");
 124         outputContainer->Add(fhPhi);
 125
 126         fhEta  = new TH1F ("hEta","#eta distribution",netabins,etamin,etamax);
 127         fhEta->SetXTitle("#eta ");
 128         outputContainer->Add(fhEta);
 129
 130         if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
 131
 132                 //Calo cells
 133                 fhNCells  = new TH1F ("hNCells","# cells per event", 100,0,1000);
 134                 fhNCells->SetXTitle("n cells");
 135                 outputContainer->Add(fhNCells);
 136
 137                 fhAmplitude  = new TH1F ("hAmplitude","#eta distribution", 100,0,1000);
 138                 fhAmplitude->SetXTitle("Amplitude ");
 139                 outputContainer->Add(fhAmplitude);
 140         }
 141
 142         if(IsDataMC()){
 143                 fh2Pt  = new TH2F ("h2Pt","p_T distribution, reconstructed vs generated", nptbins,ptmin,ptmax,nptbins,ptmin,ptmax);
 144                 fh2Pt->SetXTitle("p_{T,rec} (GeV/c)");
 145                 fh2Pt->SetYTitle("p_{T,gen} (GeV/c)");
 146                 outputContainer->Add(fh2Pt);
 147
 148                 fh2Phi  = new TH2F ("h2Phi","#phi distribution, reconstructed vs generated", nphibins,phimin,phimax, nphibins,phimin,phimax);
 149                 fh2Phi->SetXTitle("#phi_{rec} (rad)");
 150                 fh2Phi->SetYTitle("#phi_{gen} (rad)");
 151                 outputContainer->Add(fh2Phi);
 152
 153                 fh2Eta  = new TH2F ("h2Eta","#eta distribution, reconstructed vs generated", netabins,etamin,etamax,netabins,etamin,etamax);
 154                 fh2Eta->SetXTitle("#eta_{rec} ");
 155                 fh2Eta->SetYTitle("#eta_{gen} ");
 156                 outputContainer->Add(fh2Eta);
 157
 158         }
 159         return outputContainer;
 160 }
 161
 162  //__________________________________________________
 163 void AliAnaExample::InitParameters()
 164 {
 165   //Initialize the parameters of the analysis.
 166   SetOutputAODClassName("AliAODPWG4Particle");
 167   SetOutputAODName("example");
 168   fPdg = 22; //Keep photons
 169   fDetector = "PHOS";
 170
 171 }
 172
 173 //__________________________________________________________________
 174 void AliAnaExample::Print(const Option_t * opt) const
 175 {
 176   //Print some relevant parameters set for the analysis
 177   if(! opt)
 178     return;
 179
 180   printf("Min Pt = %3.2f\n", GetMinPt());
 181   printf("Max Pt = %3.2f\n", GetMaxPt());
 182   printf("Select clusters with pdg %d \n",fPdg);
 183   printf("Select detector %s \n",fDetector.Data());
 184
 185 }
 186
 187 //__________________________________________________________________
 188 void  AliAnaExample::MakeAnalysisFillAOD()
 189 {
 190         //Do analysis and fill aods
 191
 192         //Some prints
 193         if(GetDebug() > 0){
 194                 if(fDetector == "EMCAL" && GetAODEMCAL())printf("Example : in EMCAL aod entries %d\n", GetAODEMCAL()->GetEntriesFast());
 195                 if(fDetector == "CTS" && GetAODCTS())printf("Example : in CTS aod entries %d\n", GetAODCTS()->GetEntriesFast());
 196                 if(fDetector == "PHOS" && GetAODPHOS())printf("Example : in PHOS aod entries %d\n", GetAODPHOS()->GetEntriesFast());
 197         }
 198
 199         //Get List with tracks or clusters
 200         TClonesArray * partList = new TClonesArray;
 201         if(fDetector == "CTS") partList = GetAODCTS();
 202         else if(fDetector == "EMCAL") partList = GetAODEMCAL();
 203         else if(fDetector == "PHOS") partList = GetAODPHOS();
 204
 205         if(!partList || partList->GetEntriesFast() == 0) return ;
 206
 207         //Fill AODCaloClusters and AODParticle with PHOS/EMCAL aods
 208         if(fDetector == "EMCAL" || fDetector == "PHOS"){
 209
 210                 //WORK WITH CALOCLUSTERS
 211                 if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC){
 212                         ConnectAODCaloClusters(); //Do Only when filling AODCaloClusters
 213                         if(GetDebug() > 0) printf("Example: in calo clusters aod entries %d\n", GetAODCaloClusters()->GetEntriesFast());
 214                 }
 215                 //Get vertex for photon momentum calculation
 216                 Double_t v[3] ; //vertex ;
 217                 GetReader()->GetVertex(v);
 218
 219                 TLorentzVector mom ;
 220                 for(Int_t i = 0; i < partList->GetEntriesFast(); i++){
 221
 222                         AliAODCaloCluster * calo =  (AliAODCaloCluster*) (partList->At(i));
 223
 224                         //Fill AODCaloClusters
 225                         if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
 226                                 AddAODCaloCluster(AliAODCaloCluster(*(calo)));
 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(IsFidutialCutOn())
 244                                 in =  GetFidutialCut()->IsInFidutialCut(mom,fDetector) ;
 245
 246                         if(GetDebug() > 1) printf("cluster pt %2.2f, phi %2.2f, pdg %d, in fidutial 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) AliFatal("No CELLS available for analysis");
 274
 275                         //Some prints
 276                         if(GetDebug() > 0 && esdCell )
 277                                 printf("Example : in ESD %s cell entries %d\n", fDetector.Data(), esdCell->GetNumberOfCells());
 278
 279                         //Fill AODCells in file
 280                         Int_t ncells = esdCell->GetNumberOfCells() ;
 281                         GetAODCaloCells()->CreateContainer(ncells);
 282
 283                         GetAODCaloCells()->SetType((AliAODCaloCells::AODCells_t) esdCell->GetType());
 284
 285                         for (Int_t iCell = 0; iCell < ncells; iCell++) {
 286                                 if(GetDebug() > 2)  printf("cell : amp %f, absId %d,  time %f\n", esdCell->GetAmplitude(iCell), esdCell->GetCellNumber(iCell), esdCell->GetTime(iCell));
 287
 288                                 GetAODCaloCells()->SetCell(iCell,esdCell->GetCellNumber(iCell),esdCell->GetAmplitude(iCell));
 289                         }
 290                         GetAODCaloCells()->Sort();
 291                 }
 292         }//cluster-cell analysis
 293         else if(fDetector == "CTS"){ //Track analysis
 294                 //Fill AODParticle with CTS aods
 295                 TVector3 p3;
 296                 for(Int_t i = 0; i < GetAODCTS()->GetEntriesFast(); i++){
 297
 298                         AliAODTrack * track =  (AliAODTrack*) (GetAODCTS()->At(i));
 299
 300                         //Fill AODParticle after some selection
 301                         Double_t mom[3] = {track->Px(),track->Py(),track->Pz()};
 302                         p3.SetXYZ(mom[0],mom[1],mom[2]);
 303
 304                         //Acceptance selection
 305                         Bool_t in =  GetFidutialCut()->IsInFidutialCut(mom,"CTS") ;
 306                         if(GetDebug() > 1) printf("track pt %2.2f, phi %2.2f, in fidutial cut %d\n",p3.Pt(), p3.Phi(), in);
 307                         if(p3.Pt() > GetMinPt() && in) {
 308                                 AliAODPWG4Particle tr = AliAODPWG4Particle(mom[0],mom[1],mom[2],0);
 309                                 tr.SetDetector("CTS");
 310                                 AddAODParticle(tr);
 311                         }//selection
 312                 }//loop
 313         }//CTS analysis
 314
 315         if(GetDebug() > 0) {
 316                 if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
 317                         printf("Example: final aod calocluster entries %d\n", GetAODCaloClusters()->GetEntriesFast());
 318                 printf("Example: final aod branch entries %d\n", GetOutputAODBranch()->GetEntriesFast());
 319                 if(fDetector!="CTS" && GetReader()->GetDataType()!= AliCaloTrackReader::kMC)
 320                         printf("Example: final aod cell  entries %d\n", GetAODCaloCells()->GetNumberOfCells());
 321         }
 322 }
 323
 324 //__________________________________________________________________
 325 void  AliAnaExample::MakeAnalysisFillHistograms()
 326 {
 327         //Do analysis and fill histograms
 328
 329         //Loop on stored AODParticles
 330         Int_t naod = GetOutputAODBranch()->GetEntriesFast();
 331         if(GetDebug() > 0) printf("histo aod branch entries %d\n", naod);
 332         for(Int_t iaod = 0; iaod < naod ; iaod++){
 333                 AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
 334
 335                 fhPt->Fill(ph->Pt());
 336                 fhPhi->Fill(ph->Phi());
 337                 fhEta->Fill(ph->Eta());
 338
 339                 if(IsDataMC()){
 340                         //Play with the MC stack if available
 341                         AliStack * stack =  GetMCStack() ;
 342
 343                         if(ph->GetLabel() < 0 || !stack) {
 344                                 printf("*** bad label or no stack ***:  label %d \n", ph->GetLabel());
 345                                 continue;
 346                         }
 347
 348                         if(ph->GetLabel() >=  stack->GetNtrack()) {
 349                                 printf("*** large label ***:  label %d, n tracks %d \n", ph->GetLabel(), stack->GetNtrack());
 350                                 continue ;
 351                         }
 352
 353                         TParticle * mom = GetMCStack()->Particle(ph->GetLabel());
 354
 355                         fh2Pt->Fill(ph->Pt(), mom->Pt());
 356                         fh2Phi->Fill(ph->Phi(), mom->Phi());
 357                         fh2Eta->Fill(ph->Eta(), mom->Eta());
 358                 }//Work with stack also
 359         }// aod branch loop
 360
 361         // CaloCells histograms
 362         if(GetReader()->GetDataType()!= AliCaloTrackReader::kMC) {
 363                 if(GetAODCaloCells()){
 364
 365                         Int_t ncells = GetAODCaloCells()->GetNumberOfCells();
 366                         fhNCells->Fill(ncells) ;
 367
 368                         for (Int_t iCell = 0; iCell < ncells; iCell++) {
 369                                 if(GetDebug() > 2)  printf("cell : amp %f, absId %d \n", GetAODCaloCells()->GetAmplitude(iCell), GetAODCaloCells()->GetCellNumber(iCell));
 370                                 fhAmplitude->Fill(GetAODCaloCells()->GetAmplitude(iCell));
 371                         }
 372                 }//calo cells container exist
 373         }
 374 }