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