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