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