PWG4/PartCorrDep/AliAnaPi0.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 // Class to collect two-photon invariant mass distributions for
  19 // extractin raw pi0 yield.
  20 //
  21 //-- Author: Dmitri Peressounko (RRC "KI")
  22 //-- Adapted to PartCorr frame by Lamia Benhabib (SUBATECH)
  23 //-- and Gustavo Conesa (INFN-Frascati)
  24 //_________________________________________________________________________
  25
  26
  27 // --- ROOT system ---
  28 #include "TH3.h"
  29 //#include "Riostream.h"
  30 #include "TCanvas.h"
  31 #include "TPad.h"
  32 #include "TROOT.h"
  33 #include "TClonesArray.h"
  34 #include "TObjString.h"
  35
  36 //---- AliRoot system ----
  37 #include "AliAnaPi0.h"
  38 #include "AliCaloTrackReader.h"
  39 #include "AliCaloPID.h"
  40 #include "AliStack.h"
  41 #include "AliFidutialCut.h"
  42 #include "TParticle.h"
  43 #include "AliAODCaloCluster.h"
  44 #include "AliVEvent.h"
  45
  46 #ifdef __PHOSUTIL__
  47         #include "AliPHOSGeoUtils.h"
  48 #endif
  49
  50 #ifdef __EMCALUTIL__
  51         #include "AliEMCALGeoUtils.h"
  52 #endif
  53
  54 ClassImp(AliAnaPi0)
  55
  56 //________________________________________________________________________________________________________________________________________________
  57 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
  58 fNCentrBin(0),fNZvertBin(0),fNrpBin(0),
  59 fNPID(0),fNmaxMixEv(0), fZvtxCut(0.),fCalorimeter(""),
  60 fEventsList(0x0), fhEtalon(0x0),
  61 fhRe1(0x0),fhMi1(0x0),fhRe2(0x0),fhMi2(0x0),fhRe3(0x0),fhMi3(0x0),fhEvents(0x0),
  62 fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0)
  63 #ifdef __PHOSUTIL__
  64 ,fPHOSGeo(0x0)
  65 #endif
  66 #ifdef __EMCALUTIL__
  67 ,fEMCALGeo(0x0)
  68 #endif
  69
  70 {
  71 //Default Ctor
  72  InitParameters();
  73
  74 }
  75
  76 //________________________________________________________________________________________________________________________________________________
  77 AliAnaPi0::AliAnaPi0(const AliAnaPi0 & ex) : AliAnaPartCorrBaseClass(ex),
  78 fNCentrBin(ex.fNCentrBin),fNZvertBin(ex.fNZvertBin),fNrpBin(ex.fNrpBin),
  79 fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv),fZvtxCut(ex.fZvtxCut), fCalorimeter(ex.fCalorimeter),
  80 fEventsList(ex.fEventsList), fhEtalon(ex.fhEtalon),
  81 fhRe1(ex.fhRe1),fhMi1(ex.fhMi1),fhRe2(ex.fhRe2),fhMi2(ex.fhMi2),fhRe3(ex.fhRe3),fhMi3(ex.fhMi3),fhEvents(ex.fhEvents),
  82 fhPrimPt(ex.fhPrimPt), fhPrimAccPt(ex.fhPrimAccPt), fhPrimY(ex.fhPrimY),
  83 fhPrimAccY(ex.fhPrimAccY), fhPrimPhi(ex.fhPrimPhi), fhPrimAccPhi(ex.fhPrimAccPhi)
  84 #ifdef __PHOSUTIL__
  85 ,fPHOSGeo(ex.fPHOSGeo)
  86 #endif
  87 #ifdef __EMCALUTIL__
  88 ,fEMCALGeo(ex.fEMCALGeo)
  89 #endif
  90 {
  91   // cpy ctor
  92   //Do not need it
  93 }
  94
  95 //________________________________________________________________________________________________________________________________________________
  96 AliAnaPi0 & AliAnaPi0::operator = (const AliAnaPi0 & ex)
  97 {
  98   // assignment operator
  99
 100   if(this == &ex)return *this;
 101   ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
 102
 103   fNCentrBin = ex.fNCentrBin  ; fNZvertBin = ex.fNZvertBin  ; fNrpBin = ex.fNrpBin  ;
 104   fNPID = ex.fNPID  ; fNmaxMixEv = ex.fNmaxMixEv  ; fZvtxCut = ex.fZvtxCut  ;  fCalorimeter = ex.fCalorimeter  ;
 105   fEventsList = ex.fEventsList  ;  fhEtalon = ex.fhEtalon  ;
 106   fhRe1 = ex.fhRe1  ; fhMi1 = ex.fhMi1  ; fhRe2 = ex.fhRe2  ; fhMi2 = ex.fhMi2  ;
 107   fhRe3 = ex.fhRe3  ; fhMi3 = ex.fhMi3  ; fhEvents = ex.fhEvents  ;
 108   fhPrimPt = ex.fhPrimPt  ;  fhPrimAccPt = ex.fhPrimAccPt  ;  fhPrimY = ex.fhPrimY  ;
 109   fhPrimAccY = ex.fhPrimAccY  ;  fhPrimPhi = ex.fhPrimPhi  ;  fhPrimAccPhi = ex.fhPrimAccPhi ;
 110
 111   return *this;
 112
 113 }
 114
 115 //________________________________________________________________________________________________________________________________________________
 116 AliAnaPi0::~AliAnaPi0() {
 117   // Remove event containers
 118   if(fEventsList){
 119     for(Int_t ic=0; ic<fNCentrBin; ic++){
 120       for(Int_t iz=0; iz<fNZvertBin; iz++){
 121         for(Int_t irp=0; irp<fNrpBin; irp++){
 122           fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp]->Delete() ;
 123           delete fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] ;
 124         }
 125       }
 126     }
 127     delete[] fEventsList;
 128     fEventsList=0 ;
 129   }
 130
 131 #ifdef __PHOSUTIL__
 132   if(fPHOSGeo) delete fPHOSGeo ;
 133 #endif
 134
 135 #ifdef __EMCALUTIL__
 136   if(fEMCALGeo) delete fEMCALGeo ;
 137 #endif
 138
 139 }
 140
 141 //________________________________________________________________________________________________________________________________________________
 142 void AliAnaPi0::InitParameters()
 143 {
 144 //Init parameters when first called the analysis
 145 //Set default parameters
 146   SetInputAODName("PWG4Particle");
 147
 148   AddToHistogramsName("AnaPi0_");
 149
 150   fNCentrBin = 1;
 151   fNZvertBin = 1;
 152   fNrpBin    = 1;
 153   fNPID      = 9;
 154   fNmaxMixEv = 10;
 155   fZvtxCut   = 40;
 156   fCalorimeter  = "PHOS";
 157 }
 158 //________________________________________________________________________________________________________________________________________________
 159 void AliAnaPi0::Init()
 160 {
 161   //Init some data members needed in analysis
 162
 163   //Histograms binning and range
 164   if(!fhEtalon){                                                   //  p_T      alpha   d m_gg
 165     fhEtalon = new TH3D("hEtalon","Histo with binning parameters",50,0.,25.,10,0.,1.,200,0.,1.) ;
 166     fhEtalon->SetXTitle("P_{T} (GeV)") ;
 167     fhEtalon->SetYTitle("#alpha") ;
 168     fhEtalon->SetZTitle("m_{#gamma#gamma} (GeV)") ;
 169   }
 170
 171 }
 172
 173 //________________________________________________________________________________________________________________________________________________
 174 TList * AliAnaPi0::GetCreateOutputObjects()
 175 {
 176   // Create histograms to be saved in output file and
 177   // store them in fOutputContainer
 178
 179   //create event containers
 180   fEventsList = new TList*[fNCentrBin*fNZvertBin*fNrpBin] ;
 181
 182   for(Int_t ic=0; ic<fNCentrBin; ic++){
 183     for(Int_t iz=0; iz<fNZvertBin; iz++){
 184       for(Int_t irp=0; irp<fNrpBin; irp++){
 185         fEventsList[ic*fNZvertBin*fNrpBin+iz*fNrpBin+irp] = new TList() ;
 186       }
 187     }
 188   }
 189
 190   //If Geometry library loaded, do geometry selection during analysis.
 191 #ifdef __PHOSUTIL__
 192   printf("AliAnaPi0::GetCreateOutputObjects() - PHOS geometry initialized!\n");
 193   fPHOSGeo = new AliPHOSGeoUtils("PHOSgeo") ;
 194 #endif
 195
 196 #ifdef __EMCALUTIL__
 197   printf("AliAnaPi0::GetCreateOutputObjects() - EMCAL geometry initialized!\n");
 198   fEMCALGeo = new AliEMCALGeoUtils("EMCALgeo") ;
 199 #endif
 200
 201   TList * outputContainer = new TList() ;
 202   outputContainer->SetName(GetName());
 203
 204   fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
 205   fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
 206   fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
 207   fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
 208   fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
 209   fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
 210
 211   char key[255] ;
 212   char title[255] ;
 213
 214   for(Int_t ic=0; ic<fNCentrBin; ic++){
 215     for(Int_t ipid=0; ipid<fNPID; ipid++){
 216
 217       //Distance to bad module 1
 218       sprintf(key,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
 219       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 220
 221       fhEtalon->Clone(key);
 222       fhRe1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 223       fhRe1[ic*fNPID+ipid]->SetName(key) ;
 224       fhRe1[ic*fNPID+ipid]->SetTitle(title) ;
 225       outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
 226
 227       sprintf(key,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
 228       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 229       fhMi1[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 230       fhMi1[ic*fNPID+ipid]->SetName(key) ;
 231       fhMi1[ic*fNPID+ipid]->SetTitle(title) ;
 232       outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
 233
 234       //Distance to bad module 2
 235       sprintf(key,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
 236       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 237       fhRe2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 238       fhRe2[ic*fNPID+ipid]->SetName(key) ;
 239       fhRe2[ic*fNPID+ipid]->SetTitle(title) ;
 240       outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
 241
 242       sprintf(key,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
 243       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 244       fhMi2[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 245       fhMi2[ic*fNPID+ipid]->SetName(key) ;
 246       fhMi2[ic*fNPID+ipid]->SetTitle(title) ;
 247       outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
 248
 249       //Distance to bad module 3
 250       sprintf(key,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
 251       sprintf(title,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 252       fhRe3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 253       fhRe3[ic*fNPID+ipid]->SetName(key) ;
 254       fhRe3[ic*fNPID+ipid]->SetTitle(title) ;
 255       outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
 256
 257       sprintf(key,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
 258       sprintf(title,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 259       fhMi3[ic*fNPID+ipid]=(TH3D*)fhEtalon->Clone(key) ;
 260       fhMi3[ic*fNPID+ipid]->SetName(key) ;
 261       fhMi3[ic*fNPID+ipid]->SetTitle(title) ;
 262       outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
 263     }
 264   }
 265
 266
 267   fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
 268                     fNZvertBin,0.,1.*fNZvertBin,fNrpBin,0.,1.*fNrpBin) ;
 269   outputContainer->Add(fhEvents) ;
 270
 271   //Histograms filled only if MC data is requested
 272   if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
 273     if(fhEtalon->GetXaxis()->GetXbins() && fhEtalon->GetXaxis()->GetXbins()->GetSize()){ //Variable bin size
 274       fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
 275       fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",fhEtalon->GetXaxis()->GetNbins(),
 276                              fhEtalon->GetXaxis()->GetXbins()->GetArray()) ;
 277     }
 278     else{
 279       fhPrimPt = new TH1D("hPrimPt","Primary pi0 pt",fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
 280       fhPrimAccPt = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",
 281                              fhEtalon->GetXaxis()->GetNbins(),fhEtalon->GetXaxis()->GetXmin(),fhEtalon->GetXaxis()->GetXmax()) ;
 282     }
 283     outputContainer->Add(fhPrimPt) ;
 284     outputContainer->Add(fhPrimAccPt) ;
 285
 286     fhPrimY = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",100,-5.,5.) ;
 287     outputContainer->Add(fhPrimY) ;
 288
 289     fhPrimAccY = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",100,-5.,5.) ;
 290     outputContainer->Add(fhPrimAccY) ;
 291
 292     fhPrimPhi = new TH1D("hPrimaryPhi","Azimithal of primary pi0",180,0.,360.) ;
 293     outputContainer->Add(fhPrimPhi) ;
 294
 295     fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",180,-0.,360.) ;
 296     outputContainer->Add(fhPrimAccPhi) ;
 297   }
 298
 299   //Save parameters used for analysis
 300   TString parList ; //this will be list of parameters used for this analysis.
 301   char onePar[255] ;
 302   sprintf(onePar,"--- AliAnaPi0 ---\n") ;
 303   parList+=onePar ;
 304   sprintf(onePar,"Number of bins in Centrality:  %d \n",fNCentrBin) ;
 305   parList+=onePar ;
 306   sprintf(onePar,"Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
 307   parList+=onePar ;
 308   sprintf(onePar,"Number of bins in Reac. Plain: %d \n",fNrpBin) ;
 309   parList+=onePar ;
 310   sprintf(onePar,"Depth of event buffer: %d \n",fNmaxMixEv) ;
 311   parList+=onePar ;
 312   sprintf(onePar,"Number of different PID used:  %d \n",fNPID) ;
 313   parList+=onePar ;
 314   sprintf(onePar,"Cuts: \n") ;
 315   parList+=onePar ;
 316   sprintf(onePar,"Z vertex position: -%f < z < %f \n",fZvtxCut,fZvtxCut) ;
 317   parList+=onePar ;
 318   sprintf(onePar,"Calorimeter: %s \n",fCalorimeter.Data()) ;
 319   parList+=onePar ;
 320
 321   TObjString *oString= new TObjString(parList) ;
 322   outputContainer->Add(oString);
 323
 324   return outputContainer;
 325 }
 326
 327 //_________________________________________________________________________________________________________________________________________________
 328 void AliAnaPi0::Print(const Option_t * /*opt*/) const
 329 {
 330   //Print some relevant parameters set for the analysis
 331   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 332   AliAnaPartCorrBaseClass::Print(" ");
 333
 334   printf("Number of bins in Centrality:  %d \n",fNCentrBin) ;
 335   printf("Number of bins in Z vert. pos: %d \n",fNZvertBin) ;
 336   printf("Number of bins in Reac. Plain: %d \n",fNrpBin) ;
 337   printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
 338   printf("Number of different PID used:  %d \n",fNPID) ;
 339   printf("Cuts: \n") ;
 340   printf("Z vertex position: -%2.3f < z < %2.3f \n",fZvtxCut,fZvtxCut) ;
 341   printf("------------------------------------------------------\n") ;
 342 }
 343
 344
 345 //____________________________________________________________________________________________________________________________________________________
 346 void AliAnaPi0::MakeAnalysisFillHistograms()
 347 {
 348   //Process one event and extract photons from AOD branch
 349   // filled with AliAnaPhoton and fill histos with invariant mass
 350
 351   //Apply some cuts on event: vertex position and centrality range
 352   Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
 353   if(IsBadRun(iRun)) return ;
 354
 355   Double_t vert[]={0,0,0} ; //vertex ;
 356   GetReader()->GetVertex(vert);
 357   if(vert[2]<-fZvtxCut || vert[2]> fZvtxCut) return ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
 358
 359   //Get Centrality and calculate centrality bin
 360   //Does not exist in ESD yet???????
 361   Int_t curCentrBin=0 ;
 362
 363   //Get Reaction Plain position and calculate RP bin
 364   //does not exist in ESD yet????
 365   Int_t curRPBin=0 ;
 366
 367   Int_t curZvertBin=(Int_t)(0.5*fNZvertBin*(vert[2]+fZvtxCut)/fZvtxCut) ;
 368
 369   fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
 370
 371   Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
 372   if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
 373
 374   for(Int_t i1=0; i1<nPhot-1; i1++){
 375     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
 376     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
 377     for(Int_t i2=i1+1; i2<nPhot; i2++){
 378       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
 379       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
 380       Double_t m  = (photon1 + photon2).M() ;
 381       Double_t pt = (photon1 + photon2).Pt();
 382       Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
 383       if(GetDebug() > 2)
 384         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
 385                p1->Pt(), p2->Pt(), pt,m,a);
 386       for(Int_t ipid=0; ipid<fNPID; ipid++)
 387         {
 388           if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
 389             fhRe1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 390             if(p1->DistToBad()>0 && p2->DistToBad()>0){
 391               fhRe2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 392               if(p1->DistToBad()>1 && p2->DistToBad()>1){
 393                 fhRe3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 394               }
 395             }
 396           }
 397         }
 398     }
 399   }
 400
 401   //Fill mixed
 402
 403   TList * evMixList=fEventsList[curCentrBin*fNZvertBin*fNrpBin+curZvertBin*fNrpBin+curRPBin] ;
 404   Int_t nMixed = evMixList->GetSize() ;
 405   for(Int_t ii=0; ii<nMixed; ii++){
 406     TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
 407     Int_t nPhot2=ev2->GetEntriesFast() ;
 408     Double_t m = -999;
 409     if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d\n", ii, nPhot);
 410
 411     for(Int_t i1=0; i1<nPhot; i1++){
 412       AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
 413       TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
 414       for(Int_t i2=0; i2<nPhot2; i2++){
 415         AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
 416
 417         TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
 418         m =           (photon1+photon2).M() ;
 419         Double_t pt = (photon1 + photon2).Pt();
 420         Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
 421         if(GetDebug() > 2)
 422           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
 423                  p1->Pt(), p2->Pt(), pt,m,a);
 424         for(Int_t ipid=0; ipid<fNPID; ipid++){
 425           if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
 426             fhMi1[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 427             if(p1->DistToBad()>0 && p2->DistToBad()>0){
 428               fhMi2[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 429               if(p1->DistToBad()>1 && p2->DistToBad()>1){
 430                 fhMi3[curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 431               }
 432
 433             }
 434           }
 435         }
 436       }
 437     }
 438   }
 439
 440   TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
 441   //Add current event to buffer and Remove redandant events
 442   if(currentEvent->GetEntriesFast()>0){
 443     evMixList->AddFirst(currentEvent) ;
 444     currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
 445     if(evMixList->GetSize()>=fNmaxMixEv)
 446       {
 447         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
 448         evMixList->RemoveLast() ;
 449         delete tmp ;
 450       }
 451   }
 452   else{ //empty event
 453     delete currentEvent ;
 454     currentEvent=0 ;
 455   }
 456
 457   //Acceptance
 458   if(GetReader()->ReadStack()){
 459           AliStack * stack = GetMCStack();
 460           if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
 461              for(Int_t i=0 ; i<stack->GetNprimary(); i++){
 462                  TParticle * prim = stack->Particle(i) ;
 463                          if(prim->GetPdgCode() == 111){
 464                      Double_t pi0Pt = prim->Pt() ;
 465                      //printf("pi0, pt %2.2f\n",pi0Pt);
 466                     if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
 467                     Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
 468                     Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
 469                 if(TMath::Abs(pi0Y) < 0.5){
 470                    fhPrimPt->Fill(pi0Pt) ;
 471                 }
 472                 fhPrimY  ->Fill(pi0Y) ;
 473                     fhPrimPhi->Fill(phi) ;
 474
 475                 //Check if both photons hit Calorimeter
 476                 Int_t iphot1=prim->GetFirstDaughter() ;
 477                 Int_t iphot2=prim->GetLastDaughter() ;
 478                 if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
 479                    TParticle * phot1 = stack->Particle(iphot1) ;
 480                    TParticle * phot2 = stack->Particle(iphot2) ;
 481                    if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
 482                    //printf("2 photons: photon 1: pt %2.2f, phi %3.2f, eta %1.2f; photon 2: pt %2.2f, phi %3.2f, eta %1.2f\n",
 483                            //   phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
 484                    Bool_t inacceptance = kFALSE;
 485
 486                            if(fCalorimeter == "PHOS"){
 487
 488 #ifdef __PHOSUTIL__
 489             Int_t mod ;
 490             Double_t x,z ;
 491             if(fPHOSGeo->ImpactOnEmc(phot1,mod,z,x) && fPHOSGeo->ImpactOnEmc(phot2,mod,z,x))
 492               inacceptance = kTRUE;
 493             //printf("In REAL PHOS acceptance? %d\n",inacceptance);
 494 #else
 495             TLorentzVector lv1, lv2;
 496             phot1->Momentum(lv1);
 497             phot2->Momentum(lv2);
 498             if(GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) && GetFidutialCut()->IsInFidutialCut(lv2,fCalorimeter))
 499               inacceptance = kTRUE ;
 500              if(GetDebug() > 2) printf("In %s fidutial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 501 #endif
 502
 503 }
 504
 505                            else if(fCalorimeter == "EMCAL"){
 506
 507 #ifdef __EMCALUTIL__
 508             if(fEMCALGeo->Impact(phot1) && fEMCALGeo->Impact(phot2))
 509               inacceptance = kTRUE;
 510             if(GetDebug() > 2) printf("In REAL EMCAL acceptance? %d\n",inacceptance);
 511 #else
 512             TLorentzVector lv1, lv2;
 513             phot1->Momentum(lv1);
 514             phot2->Momentum(lv2);
 515             if(GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) && GetFidutialCut()->IsInFidutialCut(lv2,fCalorimeter))
 516               inacceptance = kTRUE ;
 517             //printf("In %s fidutial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 518 #endif
 519                    }
 520
 521                    if(inacceptance){
 522                       fhPrimAccPt->Fill(pi0Pt) ;
 523                       fhPrimAccPhi->Fill(phi) ;
 524                       fhPrimAccY->Fill(pi0Y) ;
 525                    }//Accepted
 526                }// 2 photons
 527              }//Check daughters exist
 528        }// Primary pi0
 529      }//loop on primaries
 530    }//stack exists and data is MC
 531   }//read stack
 532   else if(GetReader()->ReadAODMCParticles()){
 533           printf("AliAnaPi0::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
 534   }
 535
 536 }
 537
 538 //________________________________________________________________________
 539 void AliAnaPi0::ReadHistograms(TList* outputList)
 540 {
 541         // Needed when Terminate is executed in distributed environment
 542         // Refill analysis histograms of this class with corresponding histograms in output list.
 543
 544         // Histograms of this analsys are kept in the same list as other analysis, recover the position of
 545         // the first one and then add the next.
 546         Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
 547
 548         if(!fhRe1) fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
 549         if(!fhRe2) fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
 550         if(!fhRe3) fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
 551         if(!fhMi1) fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
 552         if(!fhMi2) fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
 553         if(!fhMi3) fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
 554
 555     for(Int_t ic=0; ic<fNCentrBin; ic++){
 556         for(Int_t ipid=0; ipid<fNPID; ipid++){
 557                         fhRe1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 558                         fhMi1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 559                         fhRe2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 560                         fhMi2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 561                         fhRe3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 562                         fhMi3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 563                 }
 564         }
 565
 566         fhEvents = (TH3D *) outputList->At(index++);
 567
 568         //Histograms filled only if MC data is requested
 569         if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
 570                  fhPrimPt     = (TH1D*)  outputList->At(index++);
 571                  fhPrimAccPt  = (TH1D*)  outputList->At(index++);
 572                  fhPrimY      = (TH1D*)  outputList->At(index++);
 573                  fhPrimAccY   = (TH1D*)  outputList->At(index++);
 574                  fhPrimPhi    = (TH1D*)  outputList->At(index++);
 575                  fhPrimAccPhi = (TH1D*)  outputList->At(index);
 576         }
 577
 578 }
 579
 580
 581 //____________________________________________________________________________________________________________________________________________________
 582 void AliAnaPi0::Terminate(TList* outputList)
 583 {
 584   //Do some calculations and plots from the final histograms.
 585
 586   printf(" *** %s Terminate:\n", GetName()) ;
 587
 588   //Recover histograms from output histograms list, needed for distributed analysis.
 589   ReadHistograms(outputList);
 590
 591   if (!fhRe1) {
 592      printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
 593      return;
 594   }
 595
 596   printf("AliAnaPi0::Terminate()         Mgg Real        : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(),   fhRe1[0]->GetRMS() ) ;
 597
 598   char nameIM[128];
 599   sprintf(nameIM,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
 600   TCanvas  * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
 601   cIM->Divide(2, 2);
 602
 603   cIM->cd(1) ;
 604   //gPad->SetLogy();
 605   TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ();
 606   hIMAllPt->SetLineColor(2);
 607   hIMAllPt->SetTitle("No cut on  p_{T, #gamma#gamma} ");
 608   hIMAllPt->Draw();
 609
 610   cIM->cd(2) ;
 611   TH3F * hRe1Pt5 = (TH3F*)fhRe1[0]->Clone("IMPt5");
 612   hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
 613   TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D("z");
 614   hIMPt5->SetLineColor(2);
 615   hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
 616   hIMPt5->Draw();
 617
 618   cIM->cd(3) ;
 619   TH3F * hRe1Pt10 =  (TH3F*)fhRe1[0]->Clone("IMPt10");
 620   hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
 621   TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D("z");
 622   hIMPt10->SetLineColor(2);
 623   hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
 624   hIMPt10->Draw();
 625
 626   cIM->cd(4) ;
 627   TH3F * hRe1Pt20 =  (TH3F*)fhRe1[0]->Clone("IMPt20");
 628   hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
 629   TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D("z");
 630   hIMPt20->SetLineColor(2);
 631   hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
 632   hIMPt20->Draw();
 633
 634   char nameIMF[128];
 635   sprintf(nameIMF,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
 636   cIM->Print(nameIMF);
 637
 638   char namePt[128];
 639   sprintf(namePt,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
 640   TCanvas  * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
 641   cPt->Divide(2, 2);
 642
 643   cPt->cd(1) ;
 644   //gPad->SetLogy();
 645   TH1D * hPt = (TH1D*) fhRe1[0]->Project3D("x");
 646   hPt->SetLineColor(2);
 647   hPt->SetTitle("No cut on  M_{#gamma#gamma} ");
 648   hPt->Draw();
 649
 650   cPt->cd(2) ;
 651   TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone("PtIM1");
 652   hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
 653   TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
 654   hPtIM1->SetLineColor(2);
 655   hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
 656   hPtIM1->Draw();
 657
 658   cPt->cd(3) ;
 659   TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone("PtIM2");
 660   hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
 661   TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
 662   hPtIM2->SetLineColor(2);
 663   hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
 664   hPtIM2->Draw();
 665
 666   cPt->cd(4) ;
 667   TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone("PtIM3");
 668   hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
 669   TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
 670   hPtIM3->SetLineColor(2);
 671   hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
 672   hPtIM3->Draw();
 673
 674   char namePtF[128];
 675   sprintf(namePtF,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
 676   cPt->Print(namePtF);
 677
 678
 679   char line[1024] ;
 680   sprintf(line, ".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
 681   gROOT->ProcessLine(line);
 682   sprintf(line, ".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
 683   gROOT->ProcessLine(line);
 684
 685   printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
 686
 687 }
 688
 689
 690
 691
 692