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