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 "TH2D.h"
  30 //#include "Riostream.h"
  31 #include "TCanvas.h"
  32 #include "TPad.h"
  33 #include "TROOT.h"
  34 #include "TClonesArray.h"
  35 #include "TObjString.h"
  36
  37 //---- AliRoot system ----
  38 #include "AliAnaPi0.h"
  39 #include "AliCaloTrackReader.h"
  40 #include "AliCaloPID.h"
  41 #include "AliStack.h"
  42 #include "AliFiducialCut.h"
  43 #include "TParticle.h"
  44 #include "AliVEvent.h"
  45 #include "AliESDCaloCluster.h"
  46 #include "AliESDEvent.h"
  47 #include "AliAODEvent.h"
  48 #include "AliNeutralMesonSelection.h"
  49 #include "AliMixedEvent.h"
  50
  51
  52 ClassImp(AliAnaPi0)
  53
  54 //________________________________________________________________________________________________________________________________________________
  55 AliAnaPi0::AliAnaPi0() : AliAnaPartCorrBaseClass(),
  56 fDoOwnMix(kFALSE),fNCentrBin(0),//fNZvertBin(0),fNrpBin(0),
  57 fNPID(0),fNmaxMixEv(0), fCalorimeter(""),
  58 fNModules(12), fUseAngleCut(kFALSE), fEventsList(0x0), fMultiCutAna(kFALSE),
  59 fNPtCuts(0),fNAsymCuts(0), fNCellNCuts(0),fNPIDBits(0),
  60 fhReMod(0x0),fhReDiffMod(0x0),
  61 fhRe1(0x0),      fhMi1(0x0),      fhRe2(0x0),      fhMi2(0x0),      fhRe3(0x0),      fhMi3(0x0),
  62 fhReInvPt1(0x0), fhMiInvPt1(0x0), fhReInvPt2(0x0), fhMiInvPt2(0x0), fhReInvPt3(0x0), fhMiInvPt3(0x0),
  63 fhRePtNCellAsymCuts(0x0), fhRePIDBits(0x0),
  64 fhRePtMult(0x0), fhRePtMultAsy07(0x0)  ,
  65 fhEvents(0x0), fhRealOpeningAngle(0x0),fhRealCosOpeningAngle(0x0),
  66 fhPrimPt(0x0), fhPrimAccPt(0x0), fhPrimY(0x0), fhPrimAccY(0x0), fhPrimPhi(0x0), fhPrimAccPhi(0x0),
  67 fhPrimOpeningAngle(0x0),fhPrimCosOpeningAngle(0x0)
  68 {
  69 //Default Ctor
  70  InitParameters();
  71
  72 }
  73
  74 //________________________________________________________________________________________________________________________________________________
  75 AliAnaPi0::~AliAnaPi0() {
  76   // Remove event containers
  77
  78   if(fDoOwnMix && fEventsList){
  79     for(Int_t ic=0; ic<fNCentrBin; ic++){
  80       for(Int_t iz=0; iz<GetNZvertBin(); iz++){
  81         for(Int_t irp=0; irp<GetNRPBin(); irp++){
  82           fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp]->Delete() ;
  83           delete fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] ;
  84         }
  85       }
  86     }
  87     delete[] fEventsList;
  88     fEventsList=0 ;
  89   }
  90
  91 }
  92
  93 //________________________________________________________________________________________________________________________________________________
  94 void AliAnaPi0::InitParameters()
  95 {
  96 //Init parameters when first called the analysis
  97 //Set default parameters
  98   SetInputAODName("PWG4Particle");
  99
 100   AddToHistogramsName("AnaPi0_");
 101   fNModules = 12; // set maximum to maximum number of EMCAL modules
 102   fNCentrBin = 1;
 103 //  fNZvertBin = 1;
 104 //  fNrpBin    = 1;
 105   fNPID      = 9;
 106   fNmaxMixEv = 10;
 107
 108   fCalorimeter  = "PHOS";
 109   fUseAngleCut = kFALSE;
 110
 111   fMultiCutAna = kFALSE;
 112
 113   fNPtCuts = 3;
 114   fPtCuts[0] = 0.; fPtCuts[1] = 0.2;   fPtCuts[2] = 0.3;fPtCuts[3] = 0.;   fPtCuts[4] = 0.;
 115
 116   fNAsymCuts = 3;
 117   fAsymCuts[0] = 0.7;  fAsymCuts[1] = 0.8;   fAsymCuts[2] = 1.;    fAsymCuts[3] = 0.;   fAsymCuts[4] = 0.;
 118
 119   fNCellNCuts = 3;
 120   fCellNCuts[0] = 1; fCellNCuts[1] = 2;   fCellNCuts[2] = 3;  fCellNCuts[3] = 0;   fCellNCuts[4] = 0;
 121
 122   fNPIDBits = 3;
 123   fPIDBits[0] = 2;   fPIDBits[1] = 4;   fPIDBits[2] = 6; // check dispersion, neutral, dispersion&&neutral
 124   fPIDBits[3] = 0;   fPIDBits[4] = 0;
 125 }
 126
 127
 128 //________________________________________________________________________________________________________________________________________________
 129 TObjString * AliAnaPi0::GetAnalysisCuts()
 130 {
 131          //Save parameters used for analysis
 132          TString parList ; //this will be list of parameters used for this analysis.
 133    const Int_t buffersize = 255;
 134          char onePar[buffersize] ;
 135          snprintf(onePar,buffersize,"--- AliAnaPi0 ---\n") ;
 136          parList+=onePar ;
 137          snprintf(onePar,buffersize,"Number of bins in Centrality:  %d \n",fNCentrBin) ;
 138          parList+=onePar ;
 139          snprintf(onePar,buffersize,"Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
 140          parList+=onePar ;
 141          snprintf(onePar,buffersize,"Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
 142          parList+=onePar ;
 143          snprintf(onePar,buffersize,"Depth of event buffer: %d \n",fNmaxMixEv) ;
 144          parList+=onePar ;
 145          snprintf(onePar,buffersize,"Number of different PID used:  %d \n",fNPID) ;
 146          parList+=onePar ;
 147          snprintf(onePar,buffersize,"Cuts: \n") ;
 148          parList+=onePar ;
 149          snprintf(onePar,buffersize,"Z vertex position: -%f < z < %f \n",GetZvertexCut(),GetZvertexCut()) ;
 150          parList+=onePar ;
 151          snprintf(onePar,buffersize,"Calorimeter: %s \n",fCalorimeter.Data()) ;
 152          parList+=onePar ;
 153          snprintf(onePar,buffersize,"Number of modules: %d \n",fNModules) ;
 154          parList+=onePar ;
 155   if(fMultiCutAna){
 156     snprintf(onePar, buffersize," pT cuts: n = %d, pt > ",fNPtCuts) ;
 157     for(Int_t i = 0; i < fNPtCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fPtCuts[i]);
 158     parList+=onePar ;
 159
 160     snprintf(onePar,buffersize, " N cell in cluster cuts: n = %d, nCell > ",fNCellNCuts) ;
 161     for(Int_t i = 0; i < fNCellNCuts; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fCellNCuts[i]);
 162     parList+=onePar ;
 163
 164     snprintf(onePar,buffersize," Asymmetry cuts: n = %d, asymmetry < ",fNAsymCuts) ;
 165     for(Int_t i = 0; i < fNAsymCuts; i++) snprintf(onePar,buffersize,"%s %2.2f;",onePar,fAsymCuts[i]);
 166     parList+=onePar ;
 167
 168     snprintf(onePar,buffersize," PID selection bits: n = %d, PID bit =\n",fNPIDBits) ;
 169     for(Int_t i = 0; i < fNPIDBits; i++) snprintf(onePar,buffersize,"%s %d;",onePar,fPIDBits[i]);
 170     parList+=onePar ;
 171   }
 172
 173          return new TObjString(parList) ;
 174 }
 175
 176 //________________________________________________________________________________________________________________________________________________
 177 TList * AliAnaPi0::GetCreateOutputObjects()
 178 {
 179   // Create histograms to be saved in output file and
 180   // store them in fOutputContainer
 181
 182   //create event containers
 183   fEventsList = new TList*[fNCentrBin*GetNZvertBin()*GetNRPBin()] ;
 184
 185   for(Int_t ic=0; ic<fNCentrBin; ic++){
 186     for(Int_t iz=0; iz<GetNZvertBin(); iz++){
 187       for(Int_t irp=0; irp<GetNRPBin(); irp++){
 188         fEventsList[ic*GetNZvertBin()*GetNRPBin()+iz*GetNRPBin()+irp] = new TList() ;
 189       }
 190     }
 191   }
 192
 193   TList * outputContainer = new TList() ;
 194   outputContainer->SetName(GetName());
 195
 196   fhReMod = new TH3D*[fNModules] ;
 197   fhReDiffMod = new TH3D*[fNModules+1] ;
 198
 199   fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
 200   fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
 201   fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
 202   fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
 203   fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
 204   fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
 205
 206   fhReInvPt1 = new TH3D*[fNCentrBin*fNPID] ;
 207   fhReInvPt2 = new TH3D*[fNCentrBin*fNPID] ;
 208   fhReInvPt3 = new TH3D*[fNCentrBin*fNPID] ;
 209   fhMiInvPt1 = new TH3D*[fNCentrBin*fNPID] ;
 210   fhMiInvPt2 = new TH3D*[fNCentrBin*fNPID] ;
 211   fhMiInvPt3 = new TH3D*[fNCentrBin*fNPID] ;
 212
 213   const Int_t buffersize = 255;
 214   char key[buffersize] ;
 215   char title[buffersize] ;
 216
 217   Int_t nptbins   = GetHistoPtBins();
 218   Int_t nphibins  = GetHistoPhiBins();
 219   Int_t netabins  = GetHistoEtaBins();
 220   Float_t ptmax   = GetHistoPtMax();
 221   Float_t phimax  = GetHistoPhiMax();
 222   Float_t etamax  = GetHistoEtaMax();
 223   Float_t ptmin   = GetHistoPtMin();
 224   Float_t phimin  = GetHistoPhiMin();
 225   Float_t etamin  = GetHistoEtaMin();
 226
 227   Int_t nmassbins = GetHistoMassBins();
 228   Int_t nasymbins = GetHistoAsymmetryBins();
 229   Float_t massmax = GetHistoMassMax();
 230   Float_t asymmax = GetHistoAsymmetryMax();
 231   Float_t massmin = GetHistoMassMin();
 232   Float_t asymmin = GetHistoAsymmetryMin();
 233
 234   for(Int_t ic=0; ic<fNCentrBin; ic++){
 235     for(Int_t ipid=0; ipid<fNPID; ipid++){
 236
 237       //Distance to bad module 1
 238       snprintf(key, buffersize,"hRe_cen%d_pid%d_dist1",ic,ipid) ;
 239       snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 240       fhRe1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 241       outputContainer->Add(fhRe1[ic*fNPID+ipid]) ;
 242
 243       //Distance to bad module 2
 244       snprintf(key, buffersize,"hRe_cen%d_pid%d_dist2",ic,ipid) ;
 245       snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 246       fhRe2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 247       outputContainer->Add(fhRe2[ic*fNPID+ipid]) ;
 248
 249       //Distance to bad module 3
 250       snprintf(key, buffersize,"hRe_cen%d_pid%d_dist3",ic,ipid) ;
 251       snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 252       fhRe3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 253       outputContainer->Add(fhRe3[ic*fNPID+ipid]) ;
 254
 255       //Inverse pT
 256       //Distance to bad module 1
 257       snprintf(key, buffersize,"hReInvPt_cen%d_pid%d_dist1",ic,ipid) ;
 258       snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 259       fhReInvPt1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 260       outputContainer->Add(fhReInvPt1[ic*fNPID+ipid]) ;
 261
 262       //Distance to bad module 2
 263       snprintf(key, buffersize,"hReInvPt_cen%d_pid%d_dist2",ic,ipid) ;
 264       snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 265       fhReInvPt2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 266       outputContainer->Add(fhReInvPt2[ic*fNPID+ipid]) ;
 267
 268       //Distance to bad module 3
 269       snprintf(key, buffersize,"hReInvPt_cen%d_pid%d_dist3",ic,ipid) ;
 270       snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 271       fhReInvPt3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 272       outputContainer->Add(fhReInvPt3[ic*fNPID+ipid]) ;
 273
 274       if(fDoOwnMix){
 275         //Distance to bad module 1
 276         snprintf(key, buffersize,"hMi_cen%d_pid%d_dist1",ic,ipid) ;
 277         snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 278         fhMi1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 279         outputContainer->Add(fhMi1[ic*fNPID+ipid]) ;
 280
 281         //Distance to bad module 2
 282         snprintf(key, buffersize,"hMi_cen%d_pid%d_dist2",ic,ipid) ;
 283         snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 284         fhMi2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 285         outputContainer->Add(fhMi2[ic*fNPID+ipid]) ;
 286
 287         //Distance to bad module 3
 288         snprintf(key, buffersize,"hMi_cen%d_pid%d_dist3",ic,ipid) ;
 289         snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 290         fhMi3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 291         outputContainer->Add(fhMi3[ic*fNPID+ipid]) ;
 292
 293         //Inverse pT
 294         //Distance to bad module 1
 295         snprintf(key, buffersize,"hMiInvPt_cen%d_pid%d_dist1",ic,ipid) ;
 296         snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 297         fhMiInvPt1[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 298         outputContainer->Add(fhMiInvPt1[ic*fNPID+ipid]) ;
 299
 300         //Distance to bad module 2
 301         snprintf(key, buffersize,"hMiInvPt_cen%d_pid%d_dist2",ic,ipid) ;
 302         snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 303         fhMiInvPt2[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 304         outputContainer->Add(fhMiInvPt2[ic*fNPID+ipid]) ;
 305
 306         //Distance to bad module 3
 307         snprintf(key, buffersize,"hMiInvPt_cen%d_pid%d_dist3",ic,ipid) ;
 308         snprintf(title, buffersize,"Mixed m_{#gamma#gamma} distr. for centrality=%d and PID=%d",ic,ipid) ;
 309         fhMiInvPt3[ic*fNPID+ipid] = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 310         outputContainer->Add(fhMiInvPt3[ic*fNPID+ipid]) ;
 311
 312
 313       }
 314     }
 315   }
 316
 317   if(fMultiCutAna){
 318
 319     fhRePIDBits         = new TH2D*[fNPIDBits];
 320     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
 321       snprintf(key,   buffersize,"hRe_pidbit%d",ipid) ;
 322       snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for PIDBit=%d",fPIDBits[ipid]) ;
 323       fhRePIDBits[ipid] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 324       outputContainer->Add(fhRePIDBits[ipid]) ;
 325     }// pid bit loop
 326
 327     fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 328     for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
 329       for(Int_t icell=0; icell<fNCellNCuts; icell++){
 330         for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
 331           snprintf(key,   buffersize,"hRe_pt%d_cell%d_asym%d",ipt,icell,iasym) ;
 332           snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for pt >%2.2f, ncell>%d and asym >%2.2f ",fPtCuts[ipt],fCellNCuts[icell], fAsymCuts[iasym]) ;
 333           Int_t index = ((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym;
 334           //printf("ipt %d, icell %d, iassym %d, index %d\n",ipt, icell, iasym, index);
 335           fhRePtNCellAsymCuts[index] = new TH2D(key,title,nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 336           outputContainer->Add(fhRePtNCellAsymCuts[index]) ;
 337         }
 338       }
 339     }
 340
 341     fhRePtMult      = new TH3D("hRePtMult","(p_{T},C,M)_{#gamma#gamma}, 0<A<1.0",      nptbins,ptmin,ptmax,40,0.,160.,nmassbins,massmin,massmax);
 342     fhRePtMultAsy07 = new TH3D("hRePtMultAsy07","(p_{T},C,M)_{#gamma#gamma}, 0<A<0.7",nptbins,ptmin,ptmax,40,0.,160.,nmassbins,massmin,massmax);
 343     outputContainer->Add(fhRePtMult) ;
 344     outputContainer->Add(fhRePtMultAsy07) ;
 345
 346   }// multi cuts analysis
 347
 348   fhEvents=new TH3D("hEvents","Number of events",fNCentrBin,0.,1.*fNCentrBin,
 349                     GetNZvertBin(),0.,1.*GetNZvertBin(),GetNRPBin(),0.,1.*GetNRPBin()) ;
 350   outputContainer->Add(fhEvents) ;
 351
 352   fhRealOpeningAngle  = new TH2D
 353   ("hRealOpeningAngle","Angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,0,0.5);
 354   fhRealOpeningAngle->SetYTitle("#theta(rad)");
 355   fhRealOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 356   outputContainer->Add(fhRealOpeningAngle) ;
 357
 358   fhRealCosOpeningAngle  = new TH2D
 359   ("hRealCosOpeningAngle","Cosinus of angle between all #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,200,-1,1);
 360   fhRealCosOpeningAngle->SetYTitle("cos (#theta) ");
 361   fhRealCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 362   outputContainer->Add(fhRealCosOpeningAngle) ;
 363
 364
 365   //Histograms filled only if MC data is requested
 366   if(IsDataMC()){
 367
 368     fhPrimPt     = new TH1D("hPrimPt","Primary pi0 pt",nptbins,ptmin,ptmax) ;
 369     fhPrimAccPt  = new TH1D("hPrimAccPt","Primary pi0 pt with both photons in acceptance",nptbins,ptmin,ptmax) ;
 370     outputContainer->Add(fhPrimPt) ;
 371     outputContainer->Add(fhPrimAccPt) ;
 372
 373     fhPrimY      = new TH1D("hPrimaryRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
 374     outputContainer->Add(fhPrimY) ;
 375
 376     fhPrimAccY   = new TH1D("hPrimAccRapidity","Rapidity of primary pi0",netabins,etamin,etamax) ;
 377     outputContainer->Add(fhPrimAccY) ;
 378
 379     fhPrimPhi    = new TH1D("hPrimaryPhi","Azimithal of primary pi0",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 380     outputContainer->Add(fhPrimPhi) ;
 381
 382     fhPrimAccPhi = new TH1D("hPrimAccPhi","Azimithal of primary pi0 with accepted daughters",nphibins,phimin*TMath::RadToDeg(),phimax*TMath::RadToDeg()) ;
 383     outputContainer->Add(fhPrimAccPhi) ;
 384
 385
 386     fhPrimOpeningAngle  = new TH2D
 387     ("hPrimOpeningAngle","Angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,0,0.5);
 388     fhPrimOpeningAngle->SetYTitle("#theta(rad)");
 389     fhPrimOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 390     outputContainer->Add(fhPrimOpeningAngle) ;
 391
 392     fhPrimCosOpeningAngle  = new TH2D
 393     ("hPrimCosOpeningAngle","Cosinus of angle between all primary #gamma pair vs E_{#pi^{0}}",nptbins,ptmin,ptmax,100,-1,1);
 394     fhPrimCosOpeningAngle->SetYTitle("cos (#theta) ");
 395     fhPrimCosOpeningAngle->SetXTitle("E_{ #pi^{0}} (GeV)");
 396     outputContainer->Add(fhPrimCosOpeningAngle) ;
 397
 398   }
 399
 400   TString * pairname = new TString[fNModules];
 401   if(fCalorimeter=="EMCAL"){
 402     pairname[0]="A side (0-2)";
 403     pairname[1]="C side (1-3)";
 404     pairname[2]="Sector 0 (0-1)";
 405     pairname[3]="Sector 1 (2-3)";
 406     for(Int_t i = 4 ; i < fNModules ; i++) pairname[i]="";}
 407   if(fCalorimeter=="PHOS") {
 408     pairname[0]="(0-1)";
 409     pairname[1]="(0-2)";
 410     pairname[2]="(1-3)";
 411     for(Int_t i = 3 ; i < fNModules ; i++) pairname[i]="";}
 412
 413   for(Int_t imod=0; imod<fNModules; imod++){
 414     //Module dependent invariant mass
 415     snprintf(key, buffersize,"hReMod_%d",imod) ;
 416     snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Module %d",imod) ;
 417     fhReMod[imod]  = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 418     outputContainer->Add(fhReMod[imod]) ;
 419
 420     snprintf(key, buffersize,"hReDiffMod_%d",imod) ;
 421     snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for Different Modules: %s",(pairname[imod]).Data()) ;
 422     fhReDiffMod[imod]  = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 423     outputContainer->Add(fhReDiffMod[imod]) ;
 424   }
 425
 426   delete [] pairname;
 427
 428   snprintf(key, buffersize,"hReDiffMod_%d",fNModules) ;
 429   snprintf(title, buffersize,"Real m_{#gamma#gamma} distr. for all Modules Combination") ;
 430   fhReDiffMod[fNModules]  = new TH3D(key,title,nptbins,ptmin,ptmax,nasymbins,asymmin,asymmax,nmassbins,massmin,massmax) ;
 431   outputContainer->Add(fhReDiffMod[fNModules]) ;
 432
 433
 434 //  for(Int_t i = 0; i < outputContainer->GetEntries() ; i++){
 435 //
 436 //    printf("Histogram %d, name: %s\n ",i, outputContainer->At(i)->GetName());
 437 //
 438 //  }
 439
 440   return outputContainer;
 441 }
 442
 443 //_________________________________________________________________________________________________________________________________________________
 444 void AliAnaPi0::Print(const Option_t * /*opt*/) const
 445 {
 446   //Print some relevant parameters set for the analysis
 447   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 448   AliAnaPartCorrBaseClass::Print(" ");
 449
 450   printf("Number of bins in Centrality:  %d \n",fNCentrBin) ;
 451   printf("Number of bins in Z vert. pos: %d \n",GetNZvertBin()) ;
 452   printf("Number of bins in Reac. Plain: %d \n",GetNRPBin()) ;
 453   printf("Depth of event buffer: %d \n",fNmaxMixEv) ;
 454   printf("Number of different PID used:  %d \n",fNPID) ;
 455   printf("Cuts: \n") ;
 456   printf("Z vertex position: -%2.3f < z < %2.3f \n",GetZvertexCut(),GetZvertexCut()) ;
 457   printf("Number of modules:             %d \n",fNModules) ;
 458   printf("Select pairs with their angle: %d \n",fUseAngleCut) ;
 459   if(fMultiCutAna){
 460     printf("pT cuts: n = %d, \n",fNPtCuts) ;
 461     printf("\tpT > ");
 462     for(Int_t i = 0; i < fNPtCuts; i++) printf("%2.2f ",fPtCuts[i]);
 463     printf("GeV/c\n");
 464
 465     printf("N cell in cluster cuts: n = %d, \n",fNCellNCuts) ;
 466     printf("\tnCell > ");
 467     for(Int_t i = 0; i < fNCellNCuts; i++) printf("%d ",fCellNCuts[i]);
 468     printf("\n");
 469
 470     printf("Asymmetry cuts: n = %d, \n",fNAsymCuts) ;
 471     printf("\tasymmetry < ");
 472     for(Int_t i = 0; i < fNAsymCuts; i++) printf("%2.2f ",fAsymCuts[i]);
 473     printf("\n");
 474
 475     printf("PID selection bits: n = %d, \n",fNPIDBits) ;
 476     printf("\tPID bit = ");
 477     for(Int_t i = 0; i < fNPIDBits; i++) printf("%d ",fPIDBits[i]);
 478     printf("\n");
 479
 480   }
 481   printf("------------------------------------------------------\n") ;
 482 }
 483
 484 //_____________________________________________________________
 485 void AliAnaPi0::FillAcceptanceHistograms(){
 486   //Fill acceptance histograms if MC data is available
 487
 488   if(IsDataMC() && GetReader()->ReadStack()){
 489     AliStack * stack = GetMCStack();
 490     if(stack && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){
 491       for(Int_t i=0 ; i<stack->GetNprimary(); i++){
 492         TParticle * prim = stack->Particle(i) ;
 493         if(prim->GetPdgCode() == 111){
 494           Double_t pi0Pt = prim->Pt() ;
 495           //printf("pi0, pt %2.2f\n",pi0Pt);
 496           if(prim->Energy() == TMath::Abs(prim->Pz()))  continue ; //Protection against floating point exception
 497           Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;
 498           Double_t phi   = TMath::RadToDeg()*prim->Phi() ;
 499           if(TMath::Abs(pi0Y) < 0.5){
 500             fhPrimPt->Fill(pi0Pt) ;
 501           }
 502           fhPrimY  ->Fill(pi0Y) ;
 503           fhPrimPhi->Fill(phi) ;
 504
 505           //Check if both photons hit Calorimeter
 506           Int_t iphot1=prim->GetFirstDaughter() ;
 507           Int_t iphot2=prim->GetLastDaughter() ;
 508           if(iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){
 509             TParticle * phot1 = stack->Particle(iphot1) ;
 510             TParticle * phot2 = stack->Particle(iphot2) ;
 511             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){
 512               //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",
 513               //        phot1->Pt(), phot1->Phi()*180./3.1415, phot1->Eta(), phot2->Pt(), phot2->Phi()*180./3.1415, phot2->Eta());
 514
 515               TLorentzVector lv1, lv2;
 516               phot1->Momentum(lv1);
 517               phot2->Momentum(lv2);
 518
 519               Bool_t inacceptance = kFALSE;
 520               if(fCalorimeter == "PHOS"){
 521                 if(GetPHOSGeometry() && GetCaloUtils()->IsPHOSGeoMatrixSet()){
 522                   Int_t mod ;
 523                   Double_t x,z ;
 524                   if(GetPHOSGeometry()->ImpactOnEmc(phot1,mod,z,x) && GetPHOSGeometry()->ImpactOnEmc(phot2,mod,z,x))
 525                     inacceptance = kTRUE;
 526                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 527                 }
 528                 else{
 529
 530                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
 531                     inacceptance = kTRUE ;
 532                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 533                 }
 534
 535               }
 536               else if(fCalorimeter == "EMCAL" && GetCaloUtils()->IsEMCALGeoMatrixSet()){
 537                 if(GetEMCALGeometry()){
 538                   if(GetEMCALGeometry()->Impact(phot1) && GetEMCALGeometry()->Impact(phot2))
 539                     inacceptance = kTRUE;
 540                   if(GetDebug() > 2) printf("In %s Real acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 541                 }
 542                 else{
 543                   if(GetFiducialCut()->IsInFiducialCut(lv1,fCalorimeter) && GetFiducialCut()->IsInFiducialCut(lv2,fCalorimeter))
 544                     inacceptance = kTRUE ;
 545                   if(GetDebug() > 2) printf("In %s fiducial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);
 546                 }
 547               }
 548
 549               if(inacceptance){
 550
 551                 fhPrimAccPt->Fill(pi0Pt) ;
 552                 fhPrimAccPhi->Fill(phi) ;
 553                 fhPrimAccY->Fill(pi0Y) ;
 554                 Double_t angle  = lv1.Angle(lv2.Vect());
 555                 fhPrimOpeningAngle   ->Fill(pi0Pt,angle);
 556                 fhPrimCosOpeningAngle->Fill(pi0Pt,TMath::Cos(angle));
 557
 558               }//Accepted
 559             }// 2 photons
 560           }//Check daughters exist
 561         }// Primary pi0
 562       }//loop on primaries
 563     }//stack exists and data is MC
 564   }//read stack
 565   else if(GetReader()->ReadAODMCParticles()){
 566     if(GetDebug() >= 0)  printf("AliAnaPi0::FillAcceptanceHistograms() - Acceptance calculation with MCParticles not implemented yet\n");
 567   }
 568 }
 569
 570 //____________________________________________________________________________________________________________________________________________________
 571 void AliAnaPi0::MakeAnalysisFillHistograms()
 572 {
 573   //Process one event and extract photons from AOD branch
 574   // filled with AliAnaPhoton and fill histos with invariant mass
 575
 576   //In case of MC data, fill acceptance histograms
 577   FillAcceptanceHistograms();
 578
 579   //Apply some cuts on event: vertex position and centrality range
 580   Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
 581   if(IsBadRun(iRun)) return ;
 582
 583   Int_t nPhot = GetInputAODBranch()->GetEntriesFast() ;
 584   if(GetDebug() > 1)
 585     printf("AliAnaPi0::MakeAnalysisFillHistograms() - Photon entries %d\n", nPhot);
 586   if(nPhot < 2 )
 587     return ;
 588   Int_t module1 = -1;
 589   Int_t module2 = -1;
 590   Double_t vert[] = {0.0, 0.0, 0.0} ; //vertex
 591   Int_t evtIndex1 = 0 ;
 592   Int_t currentEvtIndex = -1 ;
 593   Int_t curCentrBin = 0 ;
 594   Int_t curRPBin    = 0 ;
 595   Int_t curZvertBin = 0 ;
 596
 597   for(Int_t i1=0; i1<nPhot-1; i1++){
 598     AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
 599     // get the event index in the mixed buffer where the photon comes from
 600     // in case of mixing with analysis frame, not own mixing
 601     evtIndex1 = GetEventIndex(p1, vert) ;
 602     //printf("charge = %d\n", track->Charge());
 603     if ( evtIndex1 == -1 )
 604       return ;
 605     if ( evtIndex1 == -2 )
 606       continue ;
 607     if (evtIndex1 != currentEvtIndex) {
 608       //Get Reaction Plan position and calculate RP bin
 609       //does not exist in ESD yet????
 610       curCentrBin = 0 ;
 611       curRPBin    = 0 ;
 612       curZvertBin = (Int_t)(0.5*GetNZvertBin()*(vert[2]+GetZvertexCut())/GetZvertexCut()) ;
 613       fhEvents->Fill(curCentrBin+0.5,curZvertBin+0.5,curRPBin+0.5) ;
 614       currentEvtIndex = evtIndex1 ;
 615     }
 616
 617     //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 1 Evt %d  Vertex : %f,%f,%f\n",evtIndex1, GetVertex(evtIndex1)[0] ,GetVertex(evtIndex1)[1],GetVertex(evtIndex1)[2]);
 618
 619     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
 620     //Get Module number
 621     module1 = GetModuleNumber(p1);
 622     for(Int_t i2=i1+1; i2<nPhot; i2++){
 623       AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i2)) ;
 624       Int_t evtIndex2 = GetEventIndex(p2, vert) ;
 625       if ( evtIndex2 == -1 )
 626         return ;
 627       if ( evtIndex2 == -2 )
 628         continue ;
 629       if (GetMixedEvent() && (evtIndex1 == evtIndex2))
 630         continue ;
 631       //printf("AliAnaPi0::MakeAnalysisFillHistograms(): Photon 2 Evt %d  Vertex : %f,%f,%f\n",evtIndex2, GetVertex(evtIndex2)[0] ,GetVertex(evtIndex2)[1],GetVertex(evtIndex2)[2]);
 632       TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
 633       //Get module number
 634       module2 = GetModuleNumber(p2);
 635       Double_t m  = (photon1 + photon2).M() ;
 636       Double_t pt = (photon1 + photon2).Pt();
 637       Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
 638       if(GetDebug() > 2)
 639         printf("AliAnaPi0::MakeAnalysisFillHistograms() - Current Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
 640                p1->Pt(), p2->Pt(), pt,m,a);
 641       //Check if opening angle is too large or too small compared to what is expected
 642       Double_t angle   = photon1.Angle(photon2.Vect());
 643       //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
 644       //printf("angle %f\n",angle);
 645       if(fUseAngleCut && angle < 0.1)
 646         continue;
 647       fhRealOpeningAngle   ->Fill(pt,angle);
 648       fhRealCosOpeningAngle->Fill(pt,TMath::Cos(angle));
 649
 650       //Fill module dependent histograms
 651       //if(module1==module2) printf("mod1 %d\n",module1);
 652       if(module1==module2 && module1 >=0 && module1<fNModules)
 653         fhReMod[module1]->Fill(pt,a,m) ;
 654       else
 655         fhReDiffMod[fNModules]->Fill(pt,a,m) ;
 656
 657       if(fCalorimeter=="EMCAL"){
 658         if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[0]->Fill(pt,a,m) ;
 659         if((module1==1 && module2==3) || (module1==3 && module2==1)) fhReDiffMod[1]->Fill(pt,a,m) ;
 660         if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[2]->Fill(pt,a,m) ;
 661         if((module1==2 && module2==3) || (module1==3 && module2==2)) fhReDiffMod[3]->Fill(pt,a,m) ;
 662       }
 663       else {
 664         if((module1==0 && module2==1) || (module1==1 && module2==0)) fhReDiffMod[0]->Fill(pt,a,m) ;
 665         if((module1==0 && module2==2) || (module1==2 && module2==0)) fhReDiffMod[1]->Fill(pt,a,m) ;
 666         if((module1==1 && module2==2) || (module1==2 && module2==1)) fhReDiffMod[2]->Fill(pt,a,m) ;
 667       }
 668
 669       for(Int_t ipid=0; ipid<fNPID; ipid++){
 670         if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
 671           fhRe1     [curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 672           fhReInvPt1[curCentrBin*fNPID+ipid]->Fill(pt,a,m,1./pt) ;
 673           if(p1->DistToBad()>0 && p2->DistToBad()>0){
 674             fhRe2     [curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 675             fhReInvPt2[curCentrBin*fNPID+ipid]->Fill(pt,a,m,1./pt) ;
 676             if(p1->DistToBad()>1 && p2->DistToBad()>1){
 677               fhRe3     [curCentrBin*fNPID+ipid]->Fill(pt,a,m) ;
 678               fhReInvPt3[curCentrBin*fNPID+ipid]->Fill(pt,a,m,1./pt) ;
 679             }// bad 3
 680           }// bad2
 681         }// bad 1
 682       }// pid loop
 683
 684       //Multi cuts analysis
 685       if(fMultiCutAna){
 686         //Histograms for different PID bits selection
 687         for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
 688
 689           if(p1->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton)    &&
 690              p2->IsPIDOK(fPIDBits[ipid],AliCaloPID::kPhoton))   fhRePIDBits[ipid]->Fill(pt,m) ;
 691
 692           //printf("ipt %d, ipid%d, name %s\n",ipt, ipid, fhRePtPIDCuts[ipt*fNPIDBits+ipid]->GetName());
 693         } // pid bit cut loop
 694
 695         //Several pt,ncell and asymetry cuts
 696         //Get the number of cells
 697         Int_t ncell1 = 0;
 698         Int_t ncell2 = 0;
 699         AliVEvent * event = GetReader()->GetInputEvent();
 700         if(event){
 701           for(Int_t iclus = 0; iclus < event->GetNumberOfCaloClusters(); iclus++){
 702             AliVCluster *cluster = event->GetCaloCluster(iclus);
 703
 704             Bool_t is = kFALSE;
 705             if     (fCalorimeter == "EMCAL" && GetReader()->IsEMCALCluster(cluster)) is = kTRUE;
 706             else if(fCalorimeter == "PHOS"  && GetReader()->IsPHOSCluster (cluster)) is = kTRUE;
 707
 708             if(is){
 709               if      (p1->GetCaloLabel(0) == cluster->GetID()) ncell1 = cluster->GetNCells();
 710               else if (p2->GetCaloLabel(0) == cluster->GetID()) ncell2 = cluster->GetNCells();
 711             } // PHOS or EMCAL cluster as requested in analysis
 712
 713             if(ncell2 > 0 &&  ncell1 > 0) break; // No need to continue the iteration
 714
 715           }
 716           //printf("e 1: %2.2f, e 2: %2.2f, ncells: n1 %d, n2 %d\n", p1->E(), p2->E(),ncell1,ncell2);
 717         }
 718         for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
 719           for(Int_t icell=0; icell<fNCellNCuts; icell++){
 720             for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
 721
 722               if(p1->Pt() >   fPtCuts[ipt]      && p2->Pt() > fPtCuts[ipt]        &&
 723                  a        <   fAsymCuts[iasym]                                    &&
 724                  ncell1   >=  fCellNCuts[icell] && ncell2   >= fCellNCuts[icell]) fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->Fill(pt,m) ;
 725
 726               //printf("ipt %d, icell%d, iasym %d, name %s\n",ipt, icell, iasym,  fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym]->GetName());
 727             }// pid bit cut loop
 728           }// icell loop
 729         }// pt cut loop
 730
 731         fhRePtMult->Fill(pt,GetTrackMultiplicity(),m) ;
 732         if(a < 0.7)
 733           fhRePtMultAsy07->Fill(pt,GetTrackMultiplicity(),m) ;
 734
 735       }// multiple cuts analysis
 736     }// second same event particle
 737   }// first cluster
 738
 739   if(fDoOwnMix){
 740     //Fill mixed
 741     TList * evMixList=fEventsList[curCentrBin*GetNZvertBin()*GetNRPBin()+curZvertBin*GetNRPBin()+curRPBin] ;
 742     Int_t nMixed = evMixList->GetSize() ;
 743     for(Int_t ii=0; ii<nMixed; ii++){
 744       TClonesArray* ev2= (TClonesArray*) (evMixList->At(ii));
 745       Int_t nPhot2=ev2->GetEntriesFast() ;
 746       Double_t m = -999;
 747       if(GetDebug() > 1) printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed event %d photon entries %d\n", ii, nPhot);
 748
 749       for(Int_t i1=0; i1<nPhot; i1++){
 750         AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i1)) ;
 751         TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());
 752         for(Int_t i2=0; i2<nPhot2; i2++){
 753           AliAODPWG4Particle * p2 = (AliAODPWG4Particle*) (ev2->At(i2)) ;
 754
 755           TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());
 756           m =           (photon1+photon2).M() ;
 757           Double_t pt = (photon1 + photon2).Pt();
 758           Double_t a  = TMath::Abs(p1->E()-p2->E())/(p1->E()+p2->E()) ;
 759
 760           //Check if opening angle is too large or too small compared to what is expected
 761           Double_t angle   = photon1.Angle(photon2.Vect());
 762           //if(fUseAngleCut && !GetNeutralMesonSelection()->IsAngleInWindow((photon1+photon2).E(),angle)) continue;
 763           if(fUseAngleCut && angle < 0.1) continue;
 764
 765           if(GetDebug() > 2)
 766             printf("AliAnaPi0::MakeAnalysisFillHistograms() - Mixed Event: pT: photon1 %2.2f, photon2 %2.2f; Pair: pT %2.2f, mass %2.3f, a %f2.3\n",
 767                    p1->Pt(), p2->Pt(), pt,m,a);
 768           for(Int_t ipid=0; ipid<fNPID; ipid++){
 769             if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){
 770               fhMi1     [curCentrBin*fNPID+ipid]->Fill(pt,   a,m) ;
 771               fhMiInvPt1[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
 772               if(p1->DistToBad()>0 && p2->DistToBad()>0){
 773                 fhMi2     [curCentrBin*fNPID+ipid]->Fill(pt,   a,m) ;
 774                 fhMiInvPt2[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
 775                 if(p1->DistToBad()>1 && p2->DistToBad()>1){
 776                   fhMi3     [curCentrBin*fNPID+ipid]->Fill(pt,   a,m) ;
 777                   fhMiInvPt3[curCentrBin*fNPID+ipid]->Fill(1./pt,a,m) ;
 778                 }
 779
 780               }
 781             }
 782           }//loop for histograms
 783         }// second cluster loop
 784       }//first cluster loop
 785     }//loop on mixed events
 786
 787     TClonesArray *currentEvent = new TClonesArray(*GetInputAODBranch());
 788     //Ad  d current event to buffer and Remove redundant events
 789     if(currentEvent->GetEntriesFast()>0){
 790       evMixList->AddFirst(currentEvent) ;
 791       currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer
 792       if(evMixList->GetSize()>=fNmaxMixEv)
 793       {
 794         TClonesArray * tmp = (TClonesArray*) (evMixList->Last()) ;
 795         evMixList->RemoveLast() ;
 796         delete tmp ;
 797       }
 798     }
 799     else{ //empty event
 800       delete currentEvent ;
 801       currentEvent=0 ;
 802     }
 803   }// DoOwnMix
 804
 805 }
 806
 807 //________________________________________________________________________
 808 void AliAnaPi0::ReadHistograms(TList* outputList)
 809 {
 810   // Needed when Terminate is executed in distributed environment
 811   // Refill analysis histograms of this class with corresponding histograms in output list.
 812
 813   // Histograms of this analsys are kept in the same list as other analysis, recover the position of
 814   // the first one and then add the next.
 815   Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"hRe_cen0_pid0_dist1"));
 816
 817   if(!fhRe1) fhRe1 = new TH3D*[fNCentrBin*fNPID] ;
 818   if(!fhRe2) fhRe2 = new TH3D*[fNCentrBin*fNPID] ;
 819   if(!fhRe3) fhRe3 = new TH3D*[fNCentrBin*fNPID] ;
 820   if(!fhMi1) fhMi1 = new TH3D*[fNCentrBin*fNPID] ;
 821   if(!fhMi2) fhMi2 = new TH3D*[fNCentrBin*fNPID] ;
 822   if(!fhMi3) fhMi3 = new TH3D*[fNCentrBin*fNPID] ;
 823   if(!fhReInvPt1) fhReInvPt1  = new TH3D*[fNCentrBin*fNPID] ;
 824   if(!fhReInvPt2) fhReInvPt2  = new TH3D*[fNCentrBin*fNPID] ;
 825   if(!fhReInvPt3) fhReInvPt3  = new TH3D*[fNCentrBin*fNPID] ;
 826   if(!fhMiInvPt1) fhMiInvPt1  = new TH3D*[fNCentrBin*fNPID] ;
 827   if(!fhMiInvPt2) fhMiInvPt2  = new TH3D*[fNCentrBin*fNPID] ;
 828   if(!fhMiInvPt3) fhMiInvPt3  = new TH3D*[fNCentrBin*fNPID] ;
 829   if(!fhReMod)    fhReMod     = new TH3D*[fNModules] ;
 830   if(!fhReDiffMod)fhReDiffMod = new TH3D*[fNModules+1] ;
 831
 832
 833   for(Int_t ic=0; ic<fNCentrBin; ic++){
 834     for(Int_t ipid=0; ipid<fNPID; ipid++){
 835       fhRe1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 836       fhRe2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 837       fhRe3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 838
 839       fhReInvPt1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 840       fhReInvPt2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 841       fhReInvPt3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 842
 843       fhMi1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 844       fhMi2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 845       fhMi3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 846
 847       fhMiInvPt1[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 848       fhMiInvPt2[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 849       fhMiInvPt3[ic*fNPID+ipid] = (TH3D*) outputList->At(index++);
 850     }
 851   }
 852
 853   if(fMultiCutAna){
 854
 855     if(!fhRePtNCellAsymCuts) fhRePtNCellAsymCuts = new TH2D*[fNPtCuts*fNAsymCuts*fNCellNCuts];
 856     if(!fhRePIDBits)         fhRePIDBits         = new TH2D*[fNPIDBits];
 857
 858     for(Int_t ipid=0; ipid<fNPIDBits; ipid++){
 859       fhRePIDBits[ipid] = (TH2D*) outputList->At(index++);
 860     }// ipid loop
 861
 862     for(Int_t ipt=0; ipt<fNPtCuts; ipt++){
 863       for(Int_t icell=0; icell<fNCellNCuts; icell++){
 864         for(Int_t iasym=0; iasym<fNAsymCuts; iasym++){
 865           fhRePtNCellAsymCuts[((ipt*fNCellNCuts)+icell)*fNAsymCuts + iasym] = (TH2D*) outputList->At(index++);
 866         }// iasym
 867       }// icell loop
 868     }// ipt loop
 869     fhRePtMult      = (TH3D*) outputList->At(index++);
 870     fhRePtMultAsy07 = (TH3D*) outputList->At(index++);
 871   }// multi cut analysis
 872
 873   fhEvents = (TH3D *) outputList->At(index++);
 874
 875   //Histograms filled only if MC data is requested
 876   if(IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC) ){
 877     fhPrimPt     = (TH1D*)  outputList->At(index++);
 878     fhPrimAccPt  = (TH1D*)  outputList->At(index++);
 879     fhPrimY      = (TH1D*)  outputList->At(index++);
 880     fhPrimAccY   = (TH1D*)  outputList->At(index++);
 881     fhPrimPhi    = (TH1D*)  outputList->At(index++);
 882     fhPrimAccPhi = (TH1D*)  outputList->At(index++);
 883   }
 884
 885   for(Int_t imod=0; imod < fNModules; imod++)
 886     fhReMod[imod] = (TH3D*) outputList->At(index++);
 887
 888
 889 }
 890
 891
 892 //____________________________________________________________________________________________________________________________________________________
 893 void AliAnaPi0::Terminate(TList* outputList)
 894 {
 895   //Do some calculations and plots from the final histograms.
 896
 897   printf(" *** %s Terminate:\n", GetName()) ;
 898
 899   //Recover histograms from output histograms list, needed for distributed analysis.
 900   ReadHistograms(outputList);
 901
 902   if (!fhRe1) {
 903     printf("AliAnaPi0::Terminate() - Error: Remote output histograms not imported in AliAnaPi0 object");
 904     return;
 905   }
 906
 907   printf("AliAnaPi0::Terminate()         Mgg Real        : %5.3f , RMS : %5.3f \n", fhRe1[0]->GetMean(),   fhRe1[0]->GetRMS() ) ;
 908
 909   const Int_t buffersize = 255;
 910
 911   char nameIM[buffersize];
 912   snprintf(nameIM, buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
 913   TCanvas  * cIM = new TCanvas(nameIM, "", 400, 10, 600, 700) ;
 914   cIM->Divide(2, 2);
 915
 916   cIM->cd(1) ;
 917   //gPad->SetLogy();
 918   TH1D * hIMAllPt = (TH1D*) fhRe1[0]->ProjectionZ(Form("IMPtAll_%s",fCalorimeter.Data()));
 919   hIMAllPt->SetLineColor(2);
 920   hIMAllPt->SetTitle("No cut on  p_{T, #gamma#gamma} ");
 921   hIMAllPt->Draw();
 922
 923   cIM->cd(2) ;
 924   TH3F * hRe1Pt5 = (TH3F*)fhRe1[0]->Clone(Form("IMPt5_%s",fCalorimeter.Data()));
 925   hRe1Pt5->GetXaxis()->SetRangeUser(0,5);
 926   TH1D * hIMPt5 = (TH1D*) hRe1Pt5->Project3D(Form("IMPt5_%s_pz",fCalorimeter.Data()));
 927   hIMPt5->SetLineColor(2);
 928   hIMPt5->SetTitle("0 < p_{T, #gamma#gamma} < 5 GeV/c");
 929   hIMPt5->Draw();
 930
 931   cIM->cd(3) ;
 932   TH3F * hRe1Pt10 =  (TH3F*)fhRe1[0]->Clone(Form("IMPt10_%s",fCalorimeter.Data()));
 933   hRe1Pt10->GetXaxis()->SetRangeUser(5,10);
 934   TH1D * hIMPt10 = (TH1D*) hRe1Pt10->Project3D(Form("IMPt10_%s_pz",fCalorimeter.Data()));
 935   hIMPt10->SetLineColor(2);
 936   hIMPt10->SetTitle("5 < p_{T, #gamma#gamma} < 10 GeV/c");
 937   hIMPt10->Draw();
 938
 939   cIM->cd(4) ;
 940   TH3F * hRe1Pt20 =  (TH3F*)fhRe1[0]->Clone(Form("IMPt20_%s",fCalorimeter.Data()));
 941   hRe1Pt20->GetXaxis()->SetRangeUser(10,20);
 942   TH1D * hIMPt20 = (TH1D*) hRe1Pt20->Project3D(Form("IMPt20_%s_pz",fCalorimeter.Data()));
 943   hIMPt20->SetLineColor(2);
 944   hIMPt20->SetTitle("10 < p_{T, #gamma#gamma} < 20 GeV/c");
 945   hIMPt20->Draw();
 946
 947   char nameIMF[buffersize];
 948   snprintf(nameIMF,buffersize,"AliAnaPi0_%s_Mgg.eps",fCalorimeter.Data());
 949   cIM->Print(nameIMF);
 950
 951   char namePt[buffersize];
 952   snprintf(namePt,buffersize,"AliAnaPi0_%s_cPt",fCalorimeter.Data());
 953   TCanvas  * cPt = new TCanvas(namePt, "", 400, 10, 600, 700) ;
 954   cPt->Divide(2, 2);
 955
 956   cPt->cd(1) ;
 957   //gPad->SetLogy();
 958   TH1D * hPt = (TH1D*) fhRe1[0]->Project3D("x");
 959   hPt->SetLineColor(2);
 960   hPt->SetTitle("No cut on  M_{#gamma#gamma} ");
 961   hPt->Draw();
 962
 963   cPt->cd(2) ;
 964   TH3F * hRe1IM1 = (TH3F*)fhRe1[0]->Clone(Form("Pt1_%s",fCalorimeter.Data()));
 965   hRe1IM1->GetZaxis()->SetRangeUser(0.05,0.21);
 966   TH1D * hPtIM1 = (TH1D*) hRe1IM1->Project3D("x");
 967   hPtIM1->SetLineColor(2);
 968   hPtIM1->SetTitle("0.05 < M_{#gamma#gamma} < 0.21 GeV/c^{2}");
 969   hPtIM1->Draw();
 970
 971   cPt->cd(3) ;
 972   TH3F * hRe1IM2 = (TH3F*)fhRe1[0]->Clone(Form("Pt2_%s",fCalorimeter.Data()));
 973   hRe1IM2->GetZaxis()->SetRangeUser(0.09,0.17);
 974   TH1D * hPtIM2 = (TH1D*) hRe1IM2->Project3D("x");
 975   hPtIM2->SetLineColor(2);
 976   hPtIM2->SetTitle("0.09 < M_{#gamma#gamma} < 0.17 GeV/c^{2}");
 977   hPtIM2->Draw();
 978
 979   cPt->cd(4) ;
 980   TH3F * hRe1IM3 = (TH3F*)fhRe1[0]->Clone(Form("Pt3_%s",fCalorimeter.Data()));
 981   hRe1IM3->GetZaxis()->SetRangeUser(0.11,0.15);
 982   TH1D * hPtIM3 = (TH1D*) hRe1IM1->Project3D("x");
 983   hPtIM3->SetLineColor(2);
 984   hPtIM3->SetTitle("0.11 < M_{#gamma#gamma} < 0.15 GeV/c^{2}");
 985   hPtIM3->Draw();
 986
 987   char namePtF[buffersize];
 988   snprintf(namePtF,buffersize,"AliAnaPi0_%s_Pt.eps",fCalorimeter.Data());
 989   cPt->Print(namePtF);
 990
 991   char line[buffersize] ;
 992   snprintf(line,buffersize,".!tar -zcf %s_%s.tar.gz *.eps", GetName(),fCalorimeter.Data()) ;
 993   gROOT->ProcessLine(line);
 994   snprintf(line, buffersize,".!rm -fR AliAnaPi0_%s*.eps",fCalorimeter.Data());
 995   gROOT->ProcessLine(line);
 996
 997   printf(" AliAnaPi0::Terminate() - !! All the eps files are in %s_%s.tar.gz !!!\n", GetName(), fCalorimeter.Data());
 998
 999 }
1000   //____________________________________________________________________________________________________________________________________________________
1001 Int_t AliAnaPi0::GetEventIndex(AliAODPWG4Particle * part, Double_t * vert)
1002 {
1003   // retieves the event index and checks the vertex
1004   //    in the mixed buffer returns -2 if vertex NOK
1005   //    for normal events   returns 0 if vertex OK and -1 if vertex NOK
1006
1007   Int_t evtIndex = -1 ;
1008   if(GetReader()->GetDataType()!=AliCaloTrackReader::kMC){
1009
1010     if (GetMixedEvent()){
1011
1012       evtIndex = GetMixedEvent()->EventIndexForCaloCluster(part->GetCaloLabel(0)) ;
1013       GetVertex(vert,evtIndex);
1014
1015       if(TMath::Abs(vert[2])> GetZvertexCut())
1016         evtIndex = -2 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
1017     } else {// Single event
1018
1019       GetVertex(vert);
1020
1021       if(TMath::Abs(vert[2])> GetZvertexCut())
1022         evtIndex = -1 ; //Event can not be used (vertex, centrality,... cuts not fulfilled)
1023       else
1024         evtIndex = 0 ;
1025     }
1026   }//No MC reader
1027   else {
1028     evtIndex = 0;
1029     vert[0] = 0. ;
1030     vert[1] = 0. ;
1031     vert[2] = 0. ;
1032   }
1033
1034   return evtIndex ;
1035 }
1036