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