PWG4/PartCorrDep/AliAnaNeutralMeson.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 /* $Id: $ */\r
  16 \r
  17 //_________________________________________________________________________\r
  18 // a general class to fill two-photon and pi0+gamma invariant mass hisograms \r
  19 // to be used to extract pi0, eta and omega raw yield.\r
  20 // also for PHOS acceptance\r
  21 //--\r
  22 //-- by Renzhuo Wan (Iopp-wuhan) May 13,2009\r
  23 //_________________________________________________________________________\r
  24 \r
  25 // --- ROOT system ---\r
  26 #include "TH3F.h"\r
  27 #include "TH2F.h"\r
  28 //#include "Riostream.h"\r
  29 #include "TCanvas.h"\r
  30 #include "TPad.h"\r
  31 #include "TROOT.h"\r
  32 #include "TLorentzVector.h"\r
  33 #include "TString.h"\r
  34 #include "TParticle.h"\r
  35 //---- AliRoot system ----\r
  36 #include "AliAnaNeutralMeson.h"\r
  37 #include "AliCaloTrackReader.h"\r
  38 #include "AliCaloPID.h"\r
  39 #include "AliStack.h"\r
  40 #include "AliFidutialCut.h"\r
  41 #include "AliVEvent.h"\r
  42 #ifdef __PHOSUTIL__\r
  43    #include "AliPHOSGeoUtils.h"\r
  44 #endif\r
  45 #ifdef __EMCALUTIL__\r
  46    #include "AliEMCALGeoUtils.h"\r
  47 #endif\r
  48 \r
  49 \r
  50 ClassImp(AliAnaNeutralMeson)\r
  51 \r
  52 //______________________________________________________________________________\r
  53 AliAnaNeutralMeson::AliAnaNeutralMeson() : AliAnaPartCorrBaseClass(),\r
  54 fAnaPi0Eta(0), fAnaOmega(0), fNPID(0),fNmaxMixEv(0), fNAsy(0),fNMod(0), \r
  55 fNHistos(0), fPerPtBin(0),\r
  56 fAsyCut(0), fModCut(0),\r
  57 fDist(0), fAsy(0), fPt(0), fMass(0),\r
  58 fInvMassCut(0),\r
  59 fNbinsAsy(0), fMinAsy(0), fMaxAsy(0),\r
  60 fNbinsPt(0), fMinPt(0), fMaxPt(0),\r
  61 fNbinsM(0), fMinM(0), fMaxM(0),\r
  62 fEventsListPi0Eta(0), fEventsListOmega(0),\r
  63 fPi0Mass(0), fPi0MassPeakWidthCut(0),  fhNClusters(0),\r
  64 fhRecPhoton(0), fhRecPhotonEtaPhi(0),\r
  65 fReal2Gamma(0), fMix2Gamma(0),fRealOmega(0),\r
  66 fMixOmegaA(0),fMixOmegaB(0), fMixOmegaC(0),\r
  67 fRealTwoGammaAsyPtM(0), fMixTwoGammaAsyPtM(0),\r
  68 fRealPi0GammaAsyPtM(0), fMixAPi0GammaAsyPtM(0), fMixBPi0GammaAsyPtM(0), fMixCPi0GammaAsyPtM(0),\r
  69 fhPrimPhotonPt(0), fhPrimPhotonAccPt(0), fhPrimPhotonY(0), fhPrimPhotonAccY(0), \r
  70 fhPrimPhotonPhi(0), fhPrimPhotonAccPhi(0),\r
  71 fhPrimPi0Pt(0), fhPrimPi0AccPt(0), fhPrimPi0Y(0), fhPrimPi0AccY(0), fhPrimPi0Phi(0), fhPrimPi0AccPhi(0), \r
  72 fhPrimEtaPt(0), fhPrimEtaAccPt(0), fhPrimEtaY(0), fhPrimEtaAccY(0), fhPrimEtaPhi(0), fhPrimEtaAccPhi(0),\r
  73 fhPrimOmegaPt(0), fhPrimOmegaAccPt(0), fhPrimOmegaY(0), fhPrimOmegaAccY(0), \r
  74 fhPrimOmegaPhi(0), fhPrimOmegaAccPhi(0),\r
  75 fCalorimeter("")\r
  76 #ifdef __PHOSUTIL__\r
  77 ,fPHOSGeo(0x0)\r
  78 #endif\r
  79 #ifdef __EMCALUTIL__\r
  80 ,fEMCALGeo(0x0)\r
  81 #endif\r
  82 {\r
  83  //Default Ctor\r
  84  InitParameters();\r
  85 }\r
  86 \r
  87 //______________________________________________________________________________\r
  88 AliAnaNeutralMeson::AliAnaNeutralMeson(const AliAnaNeutralMeson & ex) : AliAnaPartCorrBaseClass(ex),  \r
  89 fAnaPi0Eta(ex.fAnaPi0Eta), fAnaOmega(ex.fAnaOmega), fNPID(ex.fNPID),fNmaxMixEv(ex.fNmaxMixEv), fNAsy(ex.fNAsy),fNMod(ex.fNMod),\r
  90 fNHistos(ex.fNHistos), fPerPtBin(ex.fPerPtBin),\r
  91 fAsyCut(ex.fAsyCut), fModCut(ex.fModCut),\r
  92 fDist(ex.fDist), fAsy(ex.fAsy), fPt(ex.fPt), fMass(ex.fMass),\r
  93 fInvMassCut(ex.fInvMassCut), \r
  94 fNbinsAsy(ex.fNbinsAsy), fMinAsy(ex.fMinAsy), fMaxAsy(ex.fMaxAsy),\r
  95 fNbinsPt(ex.fNbinsPt), fMinPt(ex.fMinPt), fMaxPt(ex.fMaxPt),\r
  96 fNbinsM(ex.fNbinsM), fMinM(ex.fMinM), fMaxM(ex.fMaxM),\r
  97 fEventsListPi0Eta(ex.fEventsListPi0Eta), fEventsListOmega(ex.fEventsListOmega),\r
  98 fPi0Mass(ex.fPi0Mass), fPi0MassPeakWidthCut(ex.fPi0MassPeakWidthCut), fhNClusters(ex.fhNClusters),\r
  99 fhRecPhoton(ex.fhRecPhoton), fhRecPhotonEtaPhi(ex.fhRecPhotonEtaPhi),\r
 100 fReal2Gamma(ex.fReal2Gamma), fMix2Gamma(ex.fMix2Gamma),fRealOmega(ex.fRealOmega),\r
 101 fMixOmegaA(ex.fMixOmegaA),fMixOmegaB(ex.fMixOmegaB), fMixOmegaC(ex.fMixOmegaC),\r
 102 fRealTwoGammaAsyPtM(ex.fRealTwoGammaAsyPtM), fMixTwoGammaAsyPtM(ex.fMixTwoGammaAsyPtM),\r
 103 fRealPi0GammaAsyPtM(ex.fRealPi0GammaAsyPtM), fMixAPi0GammaAsyPtM(ex.fMixAPi0GammaAsyPtM),\r
 104 fMixBPi0GammaAsyPtM(ex.fMixBPi0GammaAsyPtM), fMixCPi0GammaAsyPtM(ex.fMixCPi0GammaAsyPtM),\r
 105 fhPrimPhotonPt(ex.fhPrimPhotonPt), fhPrimPhotonAccPt(ex.fhPrimPhotonAccPt), fhPrimPhotonY(ex.fhPrimPhotonY), \r
 106 fhPrimPhotonAccY(ex.fhPrimPhotonAccY), fhPrimPhotonPhi(ex.fhPrimPhotonPhi), \r
 107 fhPrimPhotonAccPhi(ex.fhPrimPhotonAccPhi),\r
 108 fhPrimPi0Pt(ex.fhPrimPi0Pt), fhPrimPi0AccPt(ex.fhPrimPi0AccPt), fhPrimPi0Y(ex.fhPrimPi0Y),\r
 109 fhPrimPi0AccY(ex.fhPrimPi0AccY), fhPrimPi0Phi(ex.fhPrimPi0Phi), fhPrimPi0AccPhi(ex.fhPrimPi0AccPhi), \r
 110 fhPrimEtaPt(ex.fhPrimEtaPt), fhPrimEtaAccPt(ex.fhPrimEtaAccPt), fhPrimEtaY(ex.fhPrimEtaY), \r
 111 fhPrimEtaAccY(ex.fhPrimEtaAccY), fhPrimEtaPhi(ex.fhPrimEtaPhi), fhPrimEtaAccPhi(ex.fhPrimEtaAccPhi),\r
 112 fhPrimOmegaPt(ex.fhPrimOmegaPt), fhPrimOmegaAccPt(ex.fhPrimOmegaAccPt), fhPrimOmegaY(ex.fhPrimOmegaY),\r
 113 fhPrimOmegaAccY(ex.fhPrimOmegaAccY), fhPrimOmegaPhi(ex.fhPrimOmegaPhi), fhPrimOmegaAccPhi(ex.fhPrimOmegaAccPhi),\r
 114 fCalorimeter("")\r
 115 #ifdef __PHOSUTIL__\r
 116 ,fPHOSGeo(ex.fPHOSGeo)\r
 117 #endif\r
 118 #ifdef __EMCALUTIL__\r
 119 ,fEMCALGeo(ex.fEMCALGeo)\r
 120 #endif\r
 121 {\r
 122  // cpy ctor\r
 123  //Do not need it\r
 124 }\r
 125 \r
 126 //______________________________________________________________________________\r
 127 AliAnaNeutralMeson & AliAnaNeutralMeson::operator = (const AliAnaNeutralMeson & ex)\r
 128 {\r
 129  // assignment operator\r
 130 \r
 131  if(this == &ex)return *this;\r
 132    ((AliAnaPartCorrBaseClass *)this)->operator=(ex);\r
 133         \r
 134   fAnaPi0Eta = ex.fAnaPi0Eta;  fAnaOmega = ex.fAnaOmega;\r
 135   fNPID=ex.fNPID; fNmaxMixEv=ex.fNmaxMixEv; fNAsy=ex.fNAsy;fNMod=ex.fNMod;\r
 136   fNHistos=ex.fNHistos; fPerPtBin=ex.fPerPtBin;\r
 137   fAsyCut=ex.fAsyCut; fModCut=ex.fModCut;\r
 138   fDist=ex.fDist; fAsy=ex.fAsy; fPt=ex.fPt; fMass=ex.fMass;\r
 139   fInvMassCut=ex.fInvMassCut;\r
 140   fNbinsAsy=ex.fNbinsAsy; fMinAsy=ex.fMinAsy; fMaxAsy=ex.fMaxAsy;\r
 141   fNbinsPt=ex.fNbinsPt; fMinPt=ex.fMinPt; fMaxPt=ex.fMaxPt;\r
 142   fNbinsM=ex.fNbinsM; fMinM=ex.fMinM; fMaxM=ex.fMaxM;\r
 143   fEventsListPi0Eta=ex.fEventsListPi0Eta; fEventsListOmega=ex.fEventsListOmega;\r
 144   fPi0Mass=ex.fPi0Mass; fPi0MassPeakWidthCut=ex.fPi0MassPeakWidthCut; fhNClusters =ex.fhNClusters;\r
 145   fhRecPhoton=ex.fhRecPhoton; fhRecPhotonEtaPhi =ex.fhRecPhotonEtaPhi;\r
 146   fReal2Gamma=ex.fReal2Gamma; fMix2Gamma=ex.fMix2Gamma;fRealOmega=ex.fRealOmega;\r
 147   fMixOmegaA=ex.fMixOmegaA; fMixOmegaB=ex.fMixOmegaB; fMixOmegaC=ex.fMixOmegaC;\r
 148   fRealTwoGammaAsyPtM=ex.fRealTwoGammaAsyPtM; fMixTwoGammaAsyPtM =ex.fMixTwoGammaAsyPtM;\r
 149   fRealPi0GammaAsyPtM=ex.fRealPi0GammaAsyPtM; fMixAPi0GammaAsyPtM=ex.fMixAPi0GammaAsyPtM;\r
 150   fMixBPi0GammaAsyPtM=ex.fMixBPi0GammaAsyPtM; fMixCPi0GammaAsyPtM=ex.fMixCPi0GammaAsyPtM;\r
 151   fhPrimPhotonPt=ex.fhPrimPhotonPt; fhPrimPhotonAccPt=ex.fhPrimPhotonAccPt; fhPrimPhotonY=ex.fhPrimPhotonY;\r
 152   fhPrimPhotonAccY=ex.fhPrimPhotonAccY; fhPrimPhotonPhi=ex.fhPrimPhotonPhi;\r
 153   fhPrimPhotonAccPhi=ex.fhPrimPhotonAccPhi;\r
 154   fhPrimPi0Pt=ex.fhPrimPi0Pt; fhPrimPi0AccPt=ex.fhPrimPi0AccPt; fhPrimPi0Y=ex.fhPrimPi0Y;\r
 155   fhPrimPi0AccY=ex.fhPrimPi0AccY; fhPrimPi0Phi=ex.fhPrimPi0Phi; fhPrimPi0AccPhi=ex.fhPrimPi0AccPhi;\r
 156   fhPrimEtaPt=ex.fhPrimEtaPt; fhPrimEtaAccPt=ex.fhPrimEtaAccPt; fhPrimEtaY=ex.fhPrimEtaY;\r
 157   fhPrimEtaAccY=ex.fhPrimEtaAccY; fhPrimEtaPhi=ex.fhPrimEtaPhi; fhPrimEtaAccPhi=ex.fhPrimEtaAccPhi;\r
 158   fhPrimOmegaPt=ex.fhPrimOmegaPt; fhPrimOmegaAccPt=ex.fhPrimOmegaAccPt; fhPrimOmegaY=ex.fhPrimOmegaY;\r
 159   fhPrimOmegaAccY=ex.fhPrimOmegaAccY; fhPrimOmegaPhi=ex.fhPrimOmegaPhi; fhPrimOmegaAccPhi=ex.fhPrimOmegaAccPhi;\r
 160         \r
 161   return *this;\r
 162         \r
 163 }\r
 164 \r
 165 //______________________________________________________________________________\r
 166 AliAnaNeutralMeson::~AliAnaNeutralMeson() {\r
 167 \r
 168   //dtor\r
 169         \r
 170  if(fEventsListPi0Eta){\r
 171      delete[] fEventsListPi0Eta;\r
 172      fEventsListPi0Eta=0 ;\r
 173   }\r
 174  if(fEventsListOmega){\r
 175      delete[] fEventsListOmega;\r
 176      fEventsListOmega=0 ;\r
 177  }      \r
 178 \r
 179 #ifdef __PHOSUTIL__\r
 180     if(fPHOSGeo) delete fPHOSGeo ;\r
 181 #endif  \r
 182 #ifdef __EMCALUTIL__\r
 183     if(fEMCALGeo) delete fEMCALGeo ;\r
 184 #endif  \r
 185 }\r
 186 \r
 187 //______________________________________________________________________________\r
 188 void AliAnaNeutralMeson::SetCalorimeter(TString det)\r
 189 {\r
 190  //set the detector\r
 191  //for PHOS,  we consider three cases with different PHOS number \r
 192  //for EMCAL, we use the default one, I don't know, needs to be confirmed\r
 193  Int_t n=0;\r
 194  Int_t * nmod=0;\r
 195  if(det == "PHOS"){  \r
 196      n      = 3;\r
 197      nmod= new Int_t [n];\r
 198      nmod[0]=1; nmod[1]=3; nmod[2]=5;\r
 199   }\r
 200   \r
 201   \r
 202   if(det == "EMCAL") { //the default number of EMCAL modules\r
 203      n = 1;\r
 204      nmod= new Int_t [n]; \r
 205      nmod[0]=0; \r
 206   }  \r
 207   fCalorimeter = det;\r
 208   fNMod = n;\r
 209   fModCut = nmod;\r
 210 }\r
 211 \r
 212 //______________________________________________________________________________\r
 213 void AliAnaNeutralMeson::InitParameters()\r
 214 {\r
 215 //Init parameters when first called the analysis\r
 216 //Set default parameters\r
 217   SetInputAODName("photons");\r
 218   AddToHistogramsName("AnaNeutralMesons_");\r
 219         \r
 220   fInvMassCut = 1.;\r
 221  \r
 222   fAnaPi0Eta=kTRUE;\r
 223   fAnaOmega=kTRUE;\r
 224  \r
 225   fNPID      = 9;  \r
 226   fNAsy      = 4;\r
 227   fAsyCut = new Double_t[fNAsy] ;\r
 228   fAsyCut[0]=0.7; fAsyCut[1]=0.8; fAsyCut[2]=0.9; fAsyCut[3]=1;\r
 229  \r
 230   fNHistos=20; // set the max pt range from 0 to 20 GeV\r
 231   fPerPtBin=1.;//set the invariant mass spectrum in pt bin\r
 232 \r
 233   fNmaxMixEv = 6; //event buffer size\r
 234   fPi0Mass=0.135; //nominal pi0 mass in GeV/c^2\r
 235   fPi0MassPeakWidthCut=0.015; //the Pi0MassCut should dependent on pt, here we assume it 0.015 \r
 236  \r
 237   fNbinsAsy =200; \r
 238   fMinAsy =0;\r
 239   fMaxAsy =1;\r
 240   fNbinsPt =200;\r
 241   fMinPt = 0;\r
 242   fMaxPt=20;\r
 243   fNbinsM = 200;\r
 244   fMinM=0;\r
 245   fMaxM=1.;\r
 246 }\r
 247 \r
 248 //______________________________________________________________________________\r
 249 void AliAnaNeutralMeson::Init()\r
 250 {  \r
 251 //Init some data members needed in analysis\r
 252  \r
 253   fEventsListPi0Eta = new TList ;\r
 254   fEventsListOmega = new TList ;\r
 255         \r
 256 #ifdef __PHOSUTIL__\r
 257   printf("PHOS geometry initialized!\n");\r
 258   fPHOSGeo = new AliPHOSGeoUtils("PHOSgeo") ;\r
 259 #endif  \r
 260         \r
 261 #ifdef __EMCALUTIL__\r
 262   printf("EMCAL geometry initialized!\n");\r
 263   fEMCALGeo = new AliEMCALGeoUtils("EMCALgeo") ;\r
 264 #endif  \r
 265 \r
 266 }\r
 267 \r
 268 //______________________________________________________________________________\r
 269 TList * AliAnaNeutralMeson::GetCreateOutputObjects()\r
 270 {  \r
 271  // Create histograms to be saved in output file and \r
 272  // store them in fOutputContainer\r
 273         \r
 274  TList * outputContainer = new TList() ; \r
 275  outputContainer->SetName(GetName()); \r
 276  \r
 277  char key[255] ;\r
 278  char title[255] ;\r
 279  const char * detector= fCalorimeter;\r
 280 \r
 281  fhNClusters = new TH1I("fhNClusters","N clusters per event",100,0,100); //\r
 282  outputContainer->Add(fhNClusters);\r
 283 \r
 284  fhRecPhoton  = new TH1D*[fNPID*fNMod];\r
 285  fhRecPhotonEtaPhi= new TH2F*[fNPID*fNMod];\r
 286  // plot the photon cluster pt, eta and phi distribution by using the PID and module cut \r
 287  for(Int_t ipid=0;ipid<fNPID;ipid++){\r
 288     for(Int_t jmod=0;jmod<fNMod;jmod++){\r
 289        sprintf(key,"fhRecPhoton_%s_pid_%d_nmod_%d",detector,ipid,fModCut[jmod]);\r
 290        sprintf(title,"Rec. photon with pid#%d in %d modules with %s",ipid, fModCut[jmod],detector );\r
 291        fhRecPhoton[ipid*fNMod+jmod] = new TH1D(key,title,fNbinsPt,fMinPt, fMaxPt);\r
 292        fhRecPhoton[ipid*fNMod+jmod]->GetXaxis()->SetTitle("pt (GeV/c)");\r
 293        fhRecPhoton[ipid*fNMod+jmod]->GetYaxis()->SetTitle("Entries");\r
 294        outputContainer->Add(fhRecPhoton[ipid*fNMod+jmod]);\r
 295  \r
 296        sprintf(key,"fhRecPhotonEtaPhi_%s_pid_%d_nmod_%d",detector, ipid,fModCut[jmod]);\r
 297        sprintf(title,"Rec. photon eta vs phi with pid#%d in %d modules with %s",ipid, fModCut[jmod], detector);\r
 298        fhRecPhotonEtaPhi[ipid*fNMod+jmod] = new TH2F(key,title,200,-1,1,200,0,7);\r
 299        fhRecPhotonEtaPhi[ipid*fNMod+jmod]->GetXaxis()->SetTitle("#eta");\r
 300        fhRecPhotonEtaPhi[ipid*fNMod+jmod]->GetYaxis()->SetTitle("#phi (rad.)");\r
 301        outputContainer->Add(fhRecPhotonEtaPhi[ipid*fNMod+jmod]);\r
 302     }\r
 303  } \r
 304 \r
 305 \r
 306  if(fAnaPi0Eta){\r
 307     fReal2Gamma = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
 308     fMix2Gamma  = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
 309     \r
 310     fRealTwoGammaAsyPtM = new TH3F*[fNPID];\r
 311     fMixTwoGammaAsyPtM  = new TH3F*[fNPID];\r
 312     for(Int_t ipid=0;ipid<fNPID;ipid++){\r
 313        sprintf(key,"RealTwoGammaAsyPtM_%s_pid_%d",detector,ipid);\r
 314        sprintf(title, "%s Real 2#gamma asy:pt:ivm with pid#%d",detector,ipid);\r
 315        fRealTwoGammaAsyPtM[ipid] = new TH3F(key, title,\r
 316                                fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
 317        fRealTwoGammaAsyPtM[ipid]->GetYaxis()->SetTitle("2#gamma pt");\r
 318        fRealTwoGammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#gamma1}-E_{#gamma2}|/(E_{#gamma1}+E_{#gamma2})");\r
 319        fRealTwoGammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{2#gamma}");\r
 320        outputContainer->Add(fRealTwoGammaAsyPtM[ipid]);\r
 321  \r
 322        sprintf(key,"MixTwoGammaAsyPtM_%s_pid_%d",detector,ipid);\r
 323        sprintf(title, "%s Mix 2#gamma asy:pt:ivm with pid#%d",detector,ipid);\r
 324        fMixTwoGammaAsyPtM[ipid] = new TH3F(key, title,\r
 325                                fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
 326        fMixTwoGammaAsyPtM[ipid]->GetYaxis()->SetTitle("2#gamma pt");\r
 327        fMixTwoGammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#gamma1}-E_{#gamma2}|/(E_{#gamma1}+E_{#gamma2})");\r
 328        fMixTwoGammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{2#gamma}");\r
 329        outputContainer->Add(fMixTwoGammaAsyPtM[ipid]);   \r
 330  \r
 331       for(Int_t jmod=0; jmod<fNMod; jmod++){\r
 332          for(Int_t kasy=0;kasy<fNAsy;kasy++){\r
 333             for(Int_t mhist=0;mhist<fNHistos;mhist++){\r
 334                Double_t pt1=mhist*fPerPtBin;\r
 335                Double_t pt2=(mhist+1)*fPerPtBin;\r
 336                Int_t index = ipid*fNMod*fNAsy*fNHistos+jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
 337  \r
 338                sprintf(key,"fReal2Gamma_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector, ipid,fModCut[jmod],kasy,mhist);\r
 339                sprintf(title,"real 2#gamma IVM with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
 340                                                    ipid,fAsyCut[kasy], fModCut[jmod],detector, pt1, pt2);\r
 341                fReal2Gamma[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
 342                fReal2Gamma[index]->GetYaxis()->SetTitle("dN/dM_{#gamma#gamma}");\r
 343                fReal2Gamma[index]->GetXaxis()->SetTitle("M_{#gamma#gamma} (GeV/c^{2})");\r
 344                outputContainer->Add(fReal2Gamma[index]);\r
 345  \r
 346                sprintf(key,"fMix2Gamma_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector,ipid,fModCut[jmod],kasy,mhist);\r
 347                sprintf(title,"Mix 2#gamma IVM with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
 348                                                   ipid,fAsyCut[kasy],fModCut[jmod],detector,pt1, pt2);\r
 349                fMix2Gamma[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
 350                fMix2Gamma[index]->GetYaxis()->SetTitle("dN/dM_{#gamma#gamma}");\r
 351                fMix2Gamma[index]->GetXaxis()->SetTitle("M_{#gamma#gamma} (GeV/c^{2})");\r
 352                outputContainer->Add(fMix2Gamma[index]);\r
 353             } \r
 354          }\r
 355       }\r
 356    }\r
 357  }\r
 358  \r
 359  if(fAnaOmega){\r
 360     fRealOmega  = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
 361     fMixOmegaA  = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
 362     fMixOmegaB  = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
 363     fMixOmegaC  = new TH1F*[fNPID*fNAsy*fNMod*fNHistos];\r
 364     fRealPi0GammaAsyPtM=new TH3F*[fNPID];\r
 365     fMixAPi0GammaAsyPtM=new TH3F*[fNPID];\r
 366     fMixBPi0GammaAsyPtM=new TH3F*[fNPID];\r
 367     fMixCPi0GammaAsyPtM=new TH3F*[fNPID];\r
 368    \r
 369     for(Int_t ipid=0;ipid<fNPID;ipid++){\r
 370        sprintf(key,"RealPi0GammaAsyPtM_%s_pid_%d",detector,ipid);\r
 371        sprintf(title,"%s Real #omega->#pi^{0}#gamma asy:pt:ivm with pid#%d",detector,ipid);\r
 372        fRealPi0GammaAsyPtM[ipid] = new TH3F(key, title,\r
 373                                fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
 374        fRealPi0GammaAsyPtM[ipid]->GetYaxis()->SetTitle("#pi^{0}#gamma pt");\r
 375        fRealPi0GammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#pi^{0}}-E_{#gamma}|/(E_{#pi^{0}}+E_{#gamma})");\r
 376        fRealPi0GammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{#pi^{0}#gamma}");\r
 377        outputContainer->Add(fRealPi0GammaAsyPtM[ipid]);\r
 378  \r
 379        sprintf(key,"MixAPi0GammaAsyPtM_%s_pid_%d",detector, ipid);\r
 380        sprintf(title,"%s MixA #omega->#pi^{0}#gamma asy:pt:ivm with pid#%d",detector,ipid);\r
 381        fMixAPi0GammaAsyPtM[ipid] = new TH3F(key,title,\r
 382                               fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
 383        fMixAPi0GammaAsyPtM[ipid]->GetYaxis()->SetTitle("#pi^{0}#gamma pt");\r
 384        fMixAPi0GammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#pi^{0}}-E_{#gamma}|/(E_{#pi^{0}}+E_{#gamma})");\r
 385        fMixAPi0GammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{#pi^{0}#gamma}");\r
 386        outputContainer->Add(fMixAPi0GammaAsyPtM[ipid]);\r
 387 \r
 388        sprintf(key,"MixBPi0GammaAsyPtM_%s_pid_%d",detector,ipid);\r
 389        sprintf(title,"%s MixB #omega->#pi^{0}#gamma asy:pt:ivm with pid#%d",detector,ipid); \r
 390        fMixBPi0GammaAsyPtM[ipid] = new TH3F(key, title,\r
 391                                fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
 392        fMixBPi0GammaAsyPtM[ipid]->GetYaxis()->SetTitle("#pi^{0}#gamma pt");\r
 393        fMixBPi0GammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#pi^{0}}-E_{#gamma}|/(E_{#pi^{0}}+E_{#gamma})");\r
 394        fMixBPi0GammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{#pi^{0}#gamma}");\r
 395        outputContainer->Add(fMixBPi0GammaAsyPtM[ipid]);\r
 396  \r
 397        sprintf(key,"MixCPi0GammaAsyPtM_%s_pid_%d",detector,ipid);\r
 398        sprintf(title,"%s MixC #omega->#pi^{0}#gamma asy:pt:ivm with pid#%d",detector,ipid);\r
 399        fMixCPi0GammaAsyPtM[ipid] = new TH3F(key, title,\r
 400                                fNbinsAsy,fMinAsy, fMaxAsy, fNbinsPt,fMinPt, fMaxPt,fNbinsM,fMinM, fMaxM);\r
 401        fMixCPi0GammaAsyPtM[ipid]->GetYaxis()->SetTitle("#pi^{0}#gamma pt");\r
 402        fMixCPi0GammaAsyPtM[ipid]->GetXaxis()->SetTitle("A=|E_{#pi^{0}}-E_{#gamma}|/(E_{#pi^{0}}+E_{#gamma})");\r
 403        fMixCPi0GammaAsyPtM[ipid]->GetZaxis()->SetTitle("M_{#pi^{0}#gamma}");\r
 404        outputContainer->Add(fMixCPi0GammaAsyPtM[ipid]);\r
 405 \r
 406        for(Int_t jmod=0; jmod<fNMod; jmod++){\r
 407          for(Int_t kasy=0;kasy<fNAsy;kasy++){\r
 408             for(Int_t mhist=0;mhist<fNHistos;mhist++){\r
 409                Double_t pt1=mhist*fPerPtBin;\r
 410                Double_t pt2=(mhist+1)*fPerPtBin;\r
 411                Int_t index = ipid*fNMod*fNAsy*fNHistos+jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
 412                sprintf(key,"fRealOmega_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector,ipid,fModCut[jmod],kasy,mhist);\r
 413                sprintf(title,"Real #omega with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
 414                                               ipid,fAsyCut[kasy],fModCut[jmod], detector, pt1, pt2);\r
 415                fRealOmega[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
 416                fRealOmega[index]->GetYaxis()->SetTitle("dN/dM_{#pi^{0}#gamma} ");       \r
 417                fRealOmega[index]->GetXaxis()->SetTitle("M_{#pi^{0}#gamma} (GeV/c^{2})");\r
 418                outputContainer->Add(fRealOmega[index]);\r
 419               \r
 420                sprintf(key,"fMixOmegaA_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector, ipid,fModCut[jmod],kasy,mhist);\r
 421                sprintf(title,"Mix #omega A with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
 422                                                ipid,fAsyCut[kasy],fModCut[jmod],detector,pt1, pt2);\r
 423                fMixOmegaA[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
 424                fMixOmegaA[index]->GetYaxis()->SetTitle("dN/dM_{#pi^{0}#gamma}");\r
 425                fMixOmegaA[index]->GetXaxis()->SetTitle("M_{#pi^{0}#gamma} (GeV/c^{2})");\r
 426                outputContainer->Add(fMixOmegaA[index]);\r
 427               \r
 428                sprintf(key,"fMixOmegaB_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector,ipid,fModCut[jmod],kasy,mhist);\r
 429                sprintf(title,"Mix #omega B with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
 430                                                ipid,fAsyCut[kasy],fModCut[jmod],detector,pt1, pt2);\r
 431                fMixOmegaB[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
 432                fMixOmegaB[index]->GetYaxis()->SetTitle("dN/dM_{#pi^{0}#gamma}");\r
 433                fMixOmegaB[index]->GetXaxis()->SetTitle("M_{#pi^{0}#gamma} (GeV/c^{2})");\r
 434                outputContainer->Add(fMixOmegaB[index]);\r
 435               \r
 436                sprintf(key,"fMixOmegaC_%s_pid_%d_mod_%d_asy_%d_pt_%d",detector,ipid,fModCut[jmod],kasy,mhist);\r
 437                sprintf(title,"Mix #omega C with pid#%d asy<%2.1f in %d modules with %s at %2.2f<pt<%2.2f",\r
 438                                                ipid,fAsyCut[kasy],fModCut[jmod],detector, pt1, pt2);\r
 439                fMixOmegaC[index] = new TH1F(key, title,fNbinsM,fMinM,fMaxM);\r
 440                fMixOmegaC[index]->GetYaxis()->SetTitle("dN/dM_{#pi^{0}#gamma}");\r
 441                fMixOmegaC[index]->GetXaxis()->SetTitle("M_{#pi^{0}#gamma} (GeV/c^{2})");\r
 442                outputContainer->Add(fMixOmegaC[index]);\r
 443             }\r
 444          }\r
 445       }\r
 446    }\r
 447  }\r
 448 \r
 449  //Histograms filled only if MC data is requested       \r
 450  //pt distribution\r
 451  //currently only PHOS included                                                                 \r
 452  if(fCalorimeter=="PHOS" && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){\r
 453       fhPrimPhotonPt = new TH1D("hPrimPhotonPt","Primary phton pt at |y|=1",fNbinsPt,fMinPt, fMaxPt) ;\r
 454       outputContainer->Add(fhPrimPhotonPt);\r
 455 \r
 456       fhPrimPhotonAccPt = new TH1D("hPrimPhotonAccPt","Primary photon pt in acceptance",fNbinsPt,fMinPt, fMaxPt);\r
 457       outputContainer->Add(fhPrimPhotonAccPt);\r
 458  \r
 459       fhPrimPi0Pt = new TH1D("hPrimPi0Pt","Primary pi0 pt at |y|=1",fNbinsPt,fMinPt, fMaxPt) ;\r
 460       outputContainer->Add(fhPrimPi0Pt);\r
 461 \r
 462       fhPrimPi0AccPt = new TH1D("hPrimPi0AccPt","Primary pi0 pt with both photons in acceptance",fNbinsPt,fMinPt, fMaxPt);\r
 463       outputContainer->Add(fhPrimPi0AccPt);\r
 464  \r
 465       fhPrimEtaPt = new TH1D("hPrimEtaPt","Primary eta pt at |y|=1",fNbinsPt,fMinPt, fMaxPt) ;\r
 466       outputContainer->Add(fhPrimEtaPt);\r
 467 \r
 468       fhPrimEtaAccPt = new TH1D("hPrimEtaAccPt","Primary eta pt with both photons in acceptance",fNbinsPt,fMinPt, fMaxPt) ;\r
 469       outputContainer->Add(fhPrimEtaAccPt);\r
 470  \r
 471       fhPrimOmegaPt = new TH1D("hPrimOmegaPt","Primary Omega pt at |y|=1",fNbinsPt,fMinPt, fMaxPt) ;\r
 472       outputContainer->Add(fhPrimOmegaPt);\r
 473 \r
 474       fhPrimOmegaAccPt = new TH1D("hPrimOmegaAccPt","Primary Omega pt with omega->pi0+gamma in acceptance",fNbinsPt,fMinPt, fMaxPt);\r
 475       outputContainer->Add(fhPrimOmegaAccPt);\r
 476  \r
 477       //for photon\r
 478       fhPrimPhotonY = new TH1D("hPrimPhotonY","Rapidity of primary pi0",100,-5.,5.) ;\r
 479       outputContainer->Add(fhPrimPhotonY) ;\r
 480  \r
 481       fhPrimPhotonAccY = new TH1D("hPrimPhotonAccY","Rapidity of primary pi0",100,-5.,5.) ;\r
 482       outputContainer->Add(fhPrimPhotonAccY) ;\r
 483  \r
 484       fhPrimPhotonPhi = new TH1D("hPrimPhotonPhi","Azimithal of primary photon",180,0.,360.) ;\r
 485       outputContainer->Add(fhPrimPhotonPhi) ;\r
 486  \r
 487       fhPrimPhotonAccPhi = new TH1D("hPrimPhotonAccPhi","Azimithal of primary photon in Acceptance",180,-0.,360.) ;\r
 488       outputContainer->Add(fhPrimPhotonAccPhi) ;\r
 489  \r
 490       //for pi0\r
 491       fhPrimPi0Y = new TH1D("hPrimPi0Y","Rapidity of primary pi0",100,-5.,5.) ;\r
 492       outputContainer->Add(fhPrimPi0Y) ;\r
 493  \r
 494       fhPrimPi0AccY = new TH1D("hPrimPi0AccY","Rapidity of primary pi0",100,-5.,5.) ;\r
 495       outputContainer->Add(fhPrimPi0AccY) ;\r
 496  \r
 497       fhPrimPi0Phi = new TH1D("hPrimPi0Phi","Azimithal of primary pi0",180,0.,360.) ;\r
 498       outputContainer->Add(fhPrimPi0Phi) ;\r
 499  \r
 500       fhPrimPi0AccPhi = new TH1D("hPrimPi0AccPhi","Azimithal of primary pi0 with accepted daughters",180,-0.,360.) ;\r
 501       outputContainer->Add(fhPrimPi0AccPhi) ;\r
 502  \r
 503       //for eta\r
 504       fhPrimEtaY = new TH1D("hPrimEtaY","Rapidity of primary Eta",100,-5.,5.) ;\r
 505       outputContainer->Add(fhPrimEtaY) ;\r
 506  \r
 507       fhPrimEtaAccY = new TH1D("hPrimEtaAccY","Rapidity of primary eta",100,-5.,5.) ;\r
 508       outputContainer->Add(fhPrimEtaAccY) ;\r
 509  \r
 510       fhPrimEtaPhi = new TH1D("hPrimEtaPhi","Azimithal of primary eta",180,0.,360.) ;\r
 511       outputContainer->Add(fhPrimEtaPhi) ;\r
 512  \r
 513       fhPrimEtaAccPhi = new TH1D("hPrimEtaAccPhi","Azimithal of primary eta with accepted daughters",180,-0.,360.) ;\r
 514       outputContainer->Add(fhPrimEtaAccPhi) ;\r
 515  \r
 516       //for omega\r
 517       fhPrimOmegaY = new TH1D("hPrimOmegaY","Rapidity of primary Omega",100,-5.,5.) ;\r
 518       outputContainer->Add(fhPrimOmegaY) ;\r
 519  \r
 520       fhPrimOmegaAccY = new TH1D("hPrimOmegaAccY","Rapidity of primary omega->pi0+gamma",100,-5.,5.) ;\r
 521       outputContainer->Add(fhPrimOmegaAccY) ;\r
 522  \r
 523       fhPrimOmegaPhi = new TH1D("hPrimOmegaPhi","Azimithal of primary omega",180,0.,360.) ;\r
 524       outputContainer->Add(fhPrimOmegaPhi) ;\r
 525  \r
 526       fhPrimOmegaAccPhi = new TH1D("hPrimOmegaAccPhi","Azimithal of primary omega->pi0+gamma with accepted daughters",180,-0.,360.) ;\r
 527       outputContainer->Add(fhPrimOmegaAccPhi);\r
 528  }\r
 529  return outputContainer;\r
 530 }\r
 531 \r
 532 //______________________________________________________________________________\r
 533 void AliAnaNeutralMeson::Print(const Option_t * /*opt*/) const\r
 534 {\r
 535   //Print some relevant parameters set for the analysis\r
 536   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
 537   AliAnaPartCorrBaseClass::Print(" ");\r
 538   printf("performing the analysis of 2gamma IVM  %d\n", fAnaPi0Eta); //IVM: Invariant Mass\r
 539   printf("performing the analysis of pi0+gamma->3gamma IVM  %d\n", fAnaOmega);\r
 540   printf("Number of Events to be mixed: %d \n",fNmaxMixEv) ;\r
 541   printf("Number of different PID used:  %d \n",fNPID) ;\r
 542   printf("Number of different Asy cut:  %d\n",fNAsy);\r
 543   printf("Asy bins:  %d  Min:  %2.1f  Max:  %2.1f\n", fNbinsAsy, fMinAsy, fMaxAsy);\r
 544   printf("Pt bins:  %d  Min:  %2.1f  Max:  %2.1f  GeV/c\n", fNbinsPt, fMinPt, fMaxPt);\r
 545   printf("IVM bins:  %d  Min:  %2.1f  Max:  %2.1f  GeV/c^2\n", fNbinsM, fMinM, fMaxM);\r
 546   printf("produce the %d histograms per %2.1fGeV/c pt bin \n", fNHistos,fPerPtBin);\r
 547   printf("Vertec, centrality Cuts, RP are not used at present \n") ;\r
 548   printf("------------------------------------------------------\n") ;\r
 549 \r
 550 } \r
 551 \r
 552 //______________________________________________________________________________\r
 553 void AliAnaNeutralMeson::MakeAnalysisFillHistograms() \r
 554 {\r
 555  //process event from AOD brach\r
 556  //extract pi0, eta and omega analysis\r
 557  Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;\r
 558  if(IsBadRun(iRun)) return ;    \r
 559 \r
 560  //vertex cut not used yet\r
 561  //centrality not used yet\r
 562  //reaction plane not used yet\r
 563 \r
 564  TClonesArray *  array1 = (TClonesArray*)GetInputAODBranch();\r
 565  Int_t nPhot = array1 ->GetEntries();\r
 566  fhNClusters->Fill(nPhot);\r
 567 \r
 568  //fill photon clusters\r
 569  for(Int_t i=0;i<nPhot;i++){\r
 570      AliAODPWG4Particle * p1 = (AliAODPWG4Particle*) (GetInputAODBranch()->At(i)) ;\r
 571      TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
 572      for(Int_t ipid=0;ipid<fNPID;ipid++){\r
 573         for(Int_t jmod=0;jmod<fNMod;jmod++){\r
 574            if(IsPhotonSelected(p1, ipid, fModCut[jmod])){\r
 575                 fhRecPhoton[ipid*fNMod+jmod]->Fill(photon1.Pt());\r
 576                 Double_t phi = photon1.Phi() ;\r
 577                 if(phi<0) phi = photon1.Phi()+2*TMath::Pi();\r
 578                 fhRecPhotonEtaPhi[ipid*fNMod+jmod]->Fill(photon1.Eta(),phi);\r
 579            } \r
 580         }\r
 581      }\r
 582  }\r
 583 \r
 584 /////////////////////////////////////////////\r
 585 //real events\r
 586  if(nPhot>=2 && fAnaPi0Eta){\r
 587 //   printf("real nPhot:  %d \n",nPhot);\r
 588  \r
 589     RealPi0Eta(array1); //real events for pi0 and eta\r
 590  \r
 591     TList * evMixListPi0Eta=fEventsListPi0Eta ;        //with cluster larger than 1\r
 592     Int_t nMixed = evMixListPi0Eta->GetSize() ;\r
 593     ////////////////////////////////////////////\r
 594     //Mixing events for pi0 and eta to reconstruct the background\r
 595     for(Int_t i1=0; i1<nMixed; i1++){\r
 596         TClonesArray* array2= (TClonesArray*) (evMixListPi0Eta->At(i1));\r
 597         MixPi0Eta(array1, array2);\r
 598      }\r
 599     ////////////////////////////////////////////\r
 600  \r
 601    //event buffer for pi0 and eta\r
 602     TClonesArray *currentEvent = new TClonesArray(*array1);\r
 603  \r
 604     if(currentEvent->GetEntriesFast()>=2){\r
 605         evMixListPi0Eta->AddFirst(currentEvent) ;\r
 606         currentEvent=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer\r
 607         if(evMixListPi0Eta->GetSize()>=fNmaxMixEv) {\r
 608             TClonesArray * tmp = (TClonesArray*) (evMixListPi0Eta->Last()) ;\r
 609             evMixListPi0Eta->RemoveLast() ;\r
 610             delete tmp ;\r
 611         }\r
 612      }\r
 613      else{ //empty event\r
 614           delete currentEvent ;\r
 615           currentEvent=0 ;\r
 616      }\r
 617  \r
 618  }\r
 619  \r
 620  if(nPhot>=3 && fAnaOmega) {\r
 621     RealOmega(array1); //real events for omega\r
 622     TList * evMixListOmega = fEventsListOmega ; //with cluster larger than 2\r
 623     Int_t nMixedOmega = evMixListOmega->GetSize() ;\r
 624  \r
 625     ////////////////////////////////////////////\r
 626     //mix for omega\r
 627     for(Int_t i1=0; i1<nMixedOmega; i1++){\r
 628         TClonesArray* array2= (TClonesArray*) (evMixListOmega->At(i1));\r
 629         MixOmegaAB(array1, array2);\r
 630         //third event\r
 631         for(Int_t i2=i1+1;i2<nMixedOmega;i2++){\r
 632             TClonesArray * array3 = (TClonesArray*)(evMixListOmega->At(i2));\r
 633 //            printf("omega  nPhot:  %d   nPhot2:  %d   nPhot3:  %d \n", nPhot,nPhot2,nPhot3);\r
 634             MixOmegaC(array1,array2,array3);\r
 635         }\r
 636      }\r
 637  \r
 638     //fill event buffer with 2 more clusters events for omega extraction\r
 639     TClonesArray *currentEventOmega = new TClonesArray(*array1);\r
 640     if(currentEventOmega->GetEntriesFast()>=3){\r
 641          evMixListOmega->AddFirst(currentEventOmega) ;\r
 642          currentEventOmega=0 ; //Now list of particles belongs to buffer and it will be deleted with buffer\r
 643          if(evMixListOmega->GetSize()>=fNmaxMixEv) {\r
 644              TClonesArray * tmp = (TClonesArray*) (evMixListOmega->Last()) ;\r
 645              evMixListOmega->RemoveLast() ;\r
 646              delete tmp ;\r
 647           }\r
 648      }  \r
 649      else{ //empty event\r
 650            delete currentEventOmega ;\r
 651            currentEventOmega=0 ;\r
 652      }      \r
 653   \r
 654  }          \r
 655 \r
 656  if(fCalorimeter=="PHOS" && (IsDataMC() || (GetReader()->GetDataType() == AliCaloTrackReader::kMC)) ){\r
 657          if(GetReader()->ReadStack()){\r
 658                  AliStack * stack = GetMCStack();\r
 659                  //photon acceptance\r
 660                  PhotonAcceptance(stack);\r
 661                  //pi0 and eta acceptance\r
 662                  Pi0EtaAcceptance(stack);\r
 663                  //omega acceptance      \r
 664                  OmegaAcceptance(stack);\r
 665          }\r
 666         else if(GetReader()->ReadAODMCParticles()){\r
 667                   printf("AliAnaNeutralMeson::MakeAnalysisFillHistograms() - Acceptance calculation with MCParticles not implemented yet\n");\r
 668         }\r
 669  }\r
 670 \r
 671 }\r
 672 \r
 673 //______________________________________________________________________________\r
 674 Bool_t AliAnaNeutralMeson::IsPhotonSelected(AliAODPWG4Particle *p1, Int_t ipid, Int_t nmod)\r
 675 {\r
 676  //select the photon to be analyzed based on pid and which module \r
 677 \r
 678  TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
 679  Int_t mod = 20 ;\r
 680 \r
 681  if(fCalorimeter == "PHOS"){\r
 682 #ifdef __PHOSUTIL__\r
 683           Double_t vtx[3]={0,0,0};   \r
 684 //        if(GetReader()->GetDataType() != AliCaloTrackReader::kMC) GetReader()->GetVertex(vtx);\r
 685         Double_t x = 0 ,z = 0 ;\r
 686         fPHOSGeo->ImpactOnEmc(vtx,p1->Theta(),p1->Phi(), mod,z,x) ;\r
 687 #endif\r
 688  }\r
 689          \r
 690 \r
 691  else if(fCalorimeter == "EMCAL"){\r
 692   \r
 693 #ifdef __EMCALUTIL__\r
 694   Int_t AbsID=0;\r
 695   //TVector3 vtx(p1->Vx(),p1->Vy(),p1->Vz());//CHECK\r
 696   TVector3 vtx(0,0,0);\r
 697   TVector3 vimpact(0,0,0);\r
 698   fEMCALGeo->ImpactOnEmcal(vtx,p1->Theta(),p1->Phi(),AbsID,vimpact);\r
 699 \r
 700   //Get from the absid the supermodule, tower and eta/phi numbers\r
 701   Int_t tower = -1;\r
 702   Int_t phiT = -1;\r
 703   Int_t etaT = -1;\r
 704   fEMCALGeo->GetCellIndex(AbsID,mod,tower,phiT,etaT); \r
 705 \r
 706 #endif\r
 707   }\r
 708 \r
 709   if(mod<=nmod && p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) return kTRUE;\r
 710   else return kFALSE;\r
 711 \r
 712 }\r
 713 \r
 714 //______________________________________________________________________________\r
 715 void AliAnaNeutralMeson::Get2GammaAsyDistPtM(AliAODPWG4Particle *p1, AliAODPWG4Particle *p2,\r
 716                                             Double_t &asy, Double_t &dist, Double_t &pt, Double_t &mass)\r
 717 {\r
 718   //get the asy, dist, pair pt and pair inv mass from the 2gammas\r
 719   TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
 720   TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
 721   if(photon1.Pt()>0 && photon2.Pt()>0){\r
 722      asy=TMath::Abs(photon1.E()-photon2.E())/(photon1.E()+photon2.E());\r
 723      Double_t phi1 = photon1.Phi();\r
 724      Double_t phi2 = photon2.Phi();\r
 725      Double_t eta1 = photon1.Eta();\r
 726      Double_t eta2 = photon2.Eta();\r
 727      dist=TMath::Sqrt((phi1-phi2)*(phi1-phi2)+(eta1-eta2)*(eta1-eta2));\r
 728      pt=(photon1+photon2).Pt();\r
 729      mass=(photon1+photon2).M();\r
 730 //   printf("asy:  %2.3f   dist:   %2.3f   pt:  %2.3f   m:  %2.3f\n",asy,dist,pt,mass);\r
 731    }\r
 732    else { asy=0; dist =0; pt =0; mass=0;}\r
 733 }\r
 734 \r
 735 //______________________________________________________________________________\r
 736 void AliAnaNeutralMeson::GetAsyDistPtM(TLorentzVector photon1, TLorentzVector photon2,\r
 737                                             Double_t &asy, Double_t &dist, Double_t &pt, Double_t &mass)\r
 738 {\r
 739   //get the asy, dist, pair pt and pair inv mass from the 2gammas \r
 740   if(photon1.Pt()>0 && photon2.Pt()>0){\r
 741      asy=TMath::Abs(photon1.E()-photon2.E())/(photon1.E()+photon2.E());\r
 742      Double_t phi1 = photon1.Phi();\r
 743      Double_t phi2 = photon2.Phi();\r
 744      Double_t eta1 = photon1.Eta();\r
 745      Double_t eta2 = photon2.Eta();\r
 746      dist=TMath::Sqrt((phi1-phi2)*(phi1-phi2)+(eta1-eta2)*(eta1-eta2));\r
 747      pt=(photon1+photon2).Pt();\r
 748      mass=(photon1+photon2).M();\r
 749 //   printf("asy:  %2.3f   dist:   %2.3f   pt:  %2.3f   m:  %2.3f\n",asy,dist,pt,mass);\r
 750    }\r
 751    else { asy=0; dist =0; pt =0; mass=0;}\r
 752 }\r
 753 \r
 754 //______________________________________________________________________________ \r
 755 void AliAnaNeutralMeson::RealPi0Eta(TClonesArray * array)\r
 756 {\r
 757  //get the 2gamma invariant mass distribution from one event\r
 758  Int_t nclusters =  array->GetEntries();\r
 759  for(Int_t ie=0;ie<nclusters;ie++){\r
 760      for(Int_t je=ie+1;je<nclusters;je++){\r
 761          AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array->At(ie);\r
 762          AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array->At(je);\r
 763          //Asy vs. pt\r
 764          Get2GammaAsyDistPtM(p1,p2,fAsy,fDist,fPt,fMass);\r
 765          for(Int_t ipid=0;ipid<fNPID;ipid++){\r
 766 \r
 767              if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){ \r
 768                 fRealTwoGammaAsyPtM[ipid]->Fill(fAsy, fPt, fMass); //default module with pid\r
 769              }\r
 770 \r
 771              for(Int_t jmod=0;jmod<fNMod;jmod++){\r
 772                  for(Int_t kasy=0;kasy<fNAsy;kasy++){\r
 773                      if(IsPhotonSelected(p1,ipid,fModCut[jmod]) && IsPhotonSelected(p2,ipid,fModCut[jmod])){      \r
 774                            Get2GammaAsyDistPtM(p1,p2,fAsy,fDist,fPt,fMass);\r
 775                            //filling the histo\r
 776                            for(Int_t mhist=0;mhist<fNHistos;mhist++){\r
 777                                Double_t pt1 =mhist *fPerPtBin;\r
 778                                Double_t pt2 =(mhist+1)*fPerPtBin;\r
 779                                Int_t index = ipid*fNMod*fNAsy*fNHistos+jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
 780                                if(fMass<fInvMassCut && fAsy<fAsyCut[kasy] \r
 781                                   && fPt<pt2 && fPt>pt1 &&fAsy && fDist && fPt && fMass)\r
 782                                fReal2Gamma[index]->Fill(fMass);\r
 783                            }//hist\r
 784                       } //photon selection    \r
 785                    }//asy\r
 786                }//mod\r
 787           }//pid\r
 788      }//je\r
 789  }//ie\r
 790  \r
 791 }\r
 792 \r
 793 //______________________________________________________________________________ \r
 794 void AliAnaNeutralMeson::MixPi0Eta(TClonesArray * array1, TClonesArray * array2)\r
 795 {\r
 796  //mixing events for pi0 and eta invariant mass analysis\r
 797  // printf("mixed events for pi0 and eta invariant mass analysis!\n");\r
 798  Int_t nclusters1 =  array1->GetEntries();\r
 799  Int_t nclusters2 =  array2->GetEntries();\r
 800  for(Int_t ie=0;ie<nclusters1;ie++){\r
 801      AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array1->At(ie);\r
 802      for(Int_t je=0;je<nclusters2;je++){\r
 803          AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array2->At(je);\r
 804          TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
 805          TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
 806          Get2GammaAsyDistPtM(p1,p2, fAsy,fDist,fPt,fMass);\r
 807 \r
 808          for(Int_t ipid=0;ipid<fNPID;ipid++){\r
 809              if((p1->IsPIDOK(ipid,AliCaloPID::kPhoton)) && (p2->IsPIDOK(ipid,AliCaloPID::kPhoton))){\r
 810                 fMixTwoGammaAsyPtM[ipid]->Fill(fAsy,fPt,fMass); //default module with pid\r
 811              }\r
 812              for(Int_t jmod=0;jmod<fNMod;jmod++){ \r
 813                 for(Int_t kasy=0;kasy<fNAsy;kasy++){\r
 814                     if(IsPhotonSelected(p1,ipid,fModCut[jmod]) && IsPhotonSelected(p2,ipid,fModCut[jmod])){\r
 815                         for(Int_t mhist=0;mhist<fNHistos;mhist++){\r
 816                             Double_t pt1 =mhist *fPerPtBin;\r
 817                             Double_t pt2 =(mhist+1)*fPerPtBin;\r
 818                             Int_t index = ipid*fNMod*fNAsy*fNHistos+jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
 819                             if(fMass<fInvMassCut && fAsy<fAsyCut[kasy] \r
 820                                && fPt<pt2 && fPt>pt1 &&fAsy&& fDist && fPt && fMass)\r
 821                             fMix2Gamma[index]->Fill(fMass);\r
 822                         } //hist\r
 823                     } //photon selected\r
 824                 } //asy\r
 825               } //mod\r
 826          }//pid\r
 827     } //cluster2\r
 828  } //cluster1                 \r
 829 }                   \r
 830 \r
 831 //______________________________________________________________________________ \r
 832 void AliAnaNeutralMeson::RealOmega(TClonesArray * array)\r
 833 {\r
 834  //omega invariant mass extraction\r
 835  //\r
 836  // 1. select the pi0 candidate from 2gamma \r
 837  // 2. combine the third photon with the pi0 candidate to get the omega invariant mass\r
 838  //\r
 839  Int_t nclusters =  array->GetEntries();\r
 840  for(Int_t ie=0;ie<nclusters;ie++){//ptc1\r
 841      for(Int_t je=ie+1;je<nclusters;je++){ //ptc2\r
 842          AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array->At(ie);\r
 843          AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array->At(je);\r
 844          //pid, mod and asy cut\r
 845          Double_t f2gammaAsy;\r
 846          Get2GammaAsyDistPtM(p1,p2,f2gammaAsy,fDist,fPt,fMass);\r
 847 \r
 848          if(TMath::Abs(fMass-fPi0Mass)<fPi0MassPeakWidthCut){//select the pi0 candidate\r
 849              for(Int_t ke=0;ke<nclusters;ke++){ //ptc3\r
 850                  AliAODPWG4Particle * p3 = (AliAODPWG4Particle *)array->At(ke);              \r
 851                  if(p3 != p2 && p3 != p1 ){ ////p3!=p1 p3!=p2\r
 852                      TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
 853                      TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
 854                      Double_t px12 =(photon1+photon2).Px();\r
 855                      Double_t py12 =(photon1+photon2).Py();\r
 856                      Double_t pz12 =(photon1+photon2).Pz();\r
 857                      Double_t e12 =(photon1+photon2).E();\r
 858                      TLorentzVector photon12(px12, py12,pz12,e12);\r
 859  \r
 860                      TLorentzVector photon3(p3->Px(),p3->Py(),p3->Pz(),p3->E());\r
 861                      GetAsyDistPtM(photon12, photon3,fAsy,fDist,fPt,fMass);\r
 862                      for(Int_t ipid=0;ipid<fNPID;ipid++){ //pid\r
 863                          if(p1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&   \r
 864                                p2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
 865                                p3->IsPIDOK(ipid,AliCaloPID::kPhoton)){\r
 866                              fRealPi0GammaAsyPtM[ipid]->Fill(fAsy, fPt, fMass); //default module, without pid\r
 867                          }\r
 868 \r
 869                          for(Int_t jmod=0;jmod<fNMod;jmod++){ //mod\r
 870                              for(Int_t kasy=0;kasy<fNAsy;kasy++){ //asy\r
 871                                 if(IsPhotonSelected(p1,ipid,fModCut[jmod]) && \r
 872                                    IsPhotonSelected(p2,ipid,fModCut[jmod]) &&\r
 873                                    IsPhotonSelected(p3,ipid,fModCut[jmod])){ //photon selection\r
 874                                  \r
 875                                    //filling the histo\r
 876                                   for(Int_t mhist=0;mhist<fNHistos;mhist++){ //hist\r
 877                                       Double_t pt1 =mhist *fPerPtBin;\r
 878                                       Double_t pt2 =(mhist+1)*fPerPtBin;\r
 879                                       Int_t index = ipid*fNMod*fNAsy*fNHistos+\r
 880                                                     jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
 881                                       if(fMass<fInvMassCut && f2gammaAsy<fAsyCut[kasy] && \r
 882                                          fPt<pt2 && fPt>pt1&& fAsy<fAsyCut[kasy] &&fAsy&&fDist&&fPt&&fMass)\r
 883                                       fRealOmega[index]->Fill(fMass);\r
 884                                    }//hist\r
 885                                  } //photon selection\r
 886                                } //asy\r
 887                           }//mod\r
 888                       }//pid\r
 889                  } //p3!=p1 p3!=p2\r
 890              } //ptc3\r
 891          } //pi0 candidate\r
 892      }//je  ptc2   \r
 893  } //ie ptc1\r
 894 }\r
 895 \r
 896 \r
 897 //______________________________________________________________________________\r
 898 void AliAnaNeutralMeson::MixOmegaAB(TClonesArray * array1, TClonesArray * array2)\r
 899 {\r
 900  //omega background\r
 901  //three omega background, we classify it A, B and C\r
 902  // --A\r
 903  // (r1_event1+r2_event1)+r3_event2\r
 904  //\r
 905  // --B\r
 906  // (r1_event1+r2_event2)+r3_event2\r
 907  // \r
 908  // --C\r
 909  // (r1_event1+r2_event2)+r3_event3\r
 910  //\r
 911  Int_t nclusters1 =  array1->GetEntries();\r
 912  Int_t nclusters2 =  array2->GetEntries();\r
 913  for(Int_t ie=0;ie<nclusters1;ie++){\r
 914      //mix Omega A\r
 915      for(Int_t je=ie+1;je<nclusters1;je++){ //.............\r
 916          AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array1->At(ie);\r
 917          AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array1->At(je);\r
 918          Double_t f2gammaAsy ;\r
 919          Get2GammaAsyDistPtM(p1,p2,f2gammaAsy,fDist,fPt,fMass);\r
 920          if(TMath::Abs(fMass-fPi0Mass)<fPi0MassPeakWidthCut ){ //pi0 candidate \r
 921             for(Int_t ke=0;ke<nclusters2;ke++){ //third photon\r
 922                 AliAODPWG4Particle * p3 = (AliAODPWG4Particle *)array2->At(ke);\r
 923                 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
 924                 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
 925                 Double_t px12 =(photon1+photon2).Px();\r
 926                 Double_t py12 =(photon1+photon2).Py();\r
 927                 Double_t pz12 =(photon1+photon2).Pz();\r
 928                 Double_t e12 =(photon1+photon2).E();\r
 929                 TLorentzVector photon12(px12, py12,pz12,e12);\r
 930                 TLorentzVector photon3(p3->Px(),p3->Py(),p3->Pz(),p3->E());\r
 931                 GetAsyDistPtM(photon12, photon3,fAsy,fDist,fPt,fMass);\r
 932                 //pid, mod and asy cut\r
 933                 for(Int_t ipid=0;ipid<fNPID;ipid++){ //pid\r
 934                     if(p1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
 935                             p2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
 936                             p3->IsPIDOK(ipid,AliCaloPID::kPhoton)){\r
 937                          fMixAPi0GammaAsyPtM[ipid]->Fill(fAsy, fPt, fMass); //default module, with pid\r
 938                     }\r
 939                     for(Int_t jmod=0;jmod<fNMod;jmod++){ //mod\r
 940                         for(Int_t kasy=0;kasy<fNAsy;kasy++){  //asy\r
 941                            if(IsPhotonSelected(p1,ipid,fModCut[jmod]) && \r
 942                               IsPhotonSelected(p2,ipid,fModCut[jmod]) &&\r
 943                               IsPhotonSelected(p3,ipid,fModCut[jmod])){ //photon selection\r
 944                                //fill the hist\r
 945                                for(Int_t mhist=0;mhist<fNHistos;mhist++){ //hist\r
 946                                    Double_t pt1 =mhist *fPerPtBin;\r
 947                                    Double_t pt2 =(mhist+1)*fPerPtBin;\r
 948                                    Int_t index = ipid*fNMod*fNAsy*fNHistos+                                                                                        jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
 949                                    if(fMass<fInvMassCut && f2gammaAsy<fAsyCut[kasy] &&\r
 950                                       fPt<pt2 && fPt>pt1&& fAsy<fAsyCut[kasy] &&fAsy&&fDist&&fPt&&fMass)\r
 951                                    fMixOmegaA[index]->Fill(fMass);\r
 952                               }//hist\r
 953                             }//photon selection\r
 954                         } //asy\r
 955                     }//mod\r
 956                 } //pid\r
 957             } //ptc3\r
 958          } //pi0 candidate\r
 959       } //mixA je ptc2\r
 960 \r
 961      //mix omega B\r
 962      for(Int_t je=0;je<nclusters2;je++){ //mixb...............\r
 963          AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array1->At(ie);\r
 964          AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array2->At(je);\r
 965          Double_t f2gammaAsy;\r
 966          Get2GammaAsyDistPtM(p1,p2,f2gammaAsy,fDist,fPt,fMass);\r
 967          if(TMath::Abs(fMass-fPi0Mass)<fPi0MassPeakWidthCut){\r
 968             for(Int_t ke=0;ke<nclusters2;ke++){ //third photon\r
 969                 AliAODPWG4Particle * p3 = (AliAODPWG4Particle *)array2->At(ke);\r
 970                 if(p3 != p2){ //p3 != p2\r
 971                     TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
 972                     TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
 973                     Double_t px12 =(photon1+photon2).Px();\r
 974                     Double_t py12 =(photon1+photon2).Py();\r
 975                     Double_t pz12 =(photon1+photon2).Pz();\r
 976                     Double_t e12 =(photon1+photon2).E();\r
 977                     TLorentzVector photon12(px12, py12,pz12,e12);\r
 978                     TLorentzVector photon3(p3->Px(),p3->Py(),p3->Pz(),p3->E());\r
 979                     GetAsyDistPtM(photon12, photon3,fAsy,fDist,fPt,fMass);\r
 980 \r
 981                     for(Int_t ipid=0;ipid<fNPID;ipid++){ //pid\r
 982                         if(p1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
 983                               p2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
 984                               p3->IsPIDOK(ipid,AliCaloPID::kPhoton)){\r
 985                            fMixBPi0GammaAsyPtM[ipid]->Fill(fAsy, fPt, fMass); //default module, with pid\r
 986                         }\r
 987                         for(Int_t jmod=0;jmod<fNMod;jmod++){ //mod\r
 988                             for(Int_t kasy=0;kasy<fNAsy;kasy++){ //asy\r
 989                                 if(IsPhotonSelected(p1,ipid,fModCut[jmod]) &&\r
 990                                    IsPhotonSelected(p2,ipid,fModCut[jmod]) &&\r
 991                                    IsPhotonSelected(p3,ipid,fModCut[jmod])){ //photon selectin\r
 992                                    //fill the hist\r
 993                                    for(Int_t mhist=0;mhist<fNHistos;mhist++){ //hist\r
 994                                        Double_t pt1 =mhist *fPerPtBin;\r
 995                                        Double_t pt2 =(mhist+1)*fPerPtBin;\r
 996                                        Int_t index = ipid*fNMod*fNAsy*fNHistos+                                                                                        jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
 997                                        if(fMass<fInvMassCut && f2gammaAsy<fAsyCut[kasy] && \r
 998                                           fPt<pt2 && fPt>pt1&& fAsy<fAsyCut[kasy] &&fAsy&&fDist&&fPt&&fMass)\r
 999                                        fMixOmegaB[index]->Fill(fMass);\r
1000                                   }//hist \r
1001                                 } //photon selection\r
1002                             }//asy\r
1003                         }//mod\r
1004                     }//pid\r
1005                 }    \r
1006             }//ptc3          \r
1007          } //pi0 candidate\r
1008      } //mix B \r
1009 \r
1010  }//ie\r
1011 \r
1012 }\r
1013 \r
1014 //______________________________________________________________________________ \r
1015 void AliAnaNeutralMeson::MixOmegaC(TClonesArray * array1, TClonesArray * array2, TClonesArray * array3)\r
1016 {\r
1017   //omega background\r
1018  //three omega background, we classify it A, B and C\r
1019  // --A\r
1020  // (r1_event1+r2_event1)+r3_event2\r
1021  //\r
1022  // --B\r
1023  // (r1_event1+r2_event2)+r3_event2\r
1024  // \r
1025  // --C\r
1026  // (r1_event1+r2_event2)+r3_event3\r
1027  Int_t nclusters1 =  array1->GetEntries();\r
1028  Int_t nclusters2 =  array2->GetEntries();\r
1029  Int_t nclusters3 =  array3->GetEntries();\r
1030  for(Int_t ie=0;ie<nclusters1;ie++){\r
1031      for(Int_t je=0;je<nclusters2;je++){\r
1032          AliAODPWG4Particle * p1 = (AliAODPWG4Particle *)array1->At(ie);\r
1033          AliAODPWG4Particle * p2 = (AliAODPWG4Particle *)array2->At(je);\r
1034          Double_t f2gammaAsy;\r
1035          Get2GammaAsyDistPtM(p1,p2,f2gammaAsy,fDist,fPt,fMass);\r
1036          if(TMath::Abs(fMass-fPi0Mass)<fPi0MassPeakWidthCut){\r
1037             for(Int_t ke=0;ke<nclusters3;ke++){ \r
1038                 AliAODPWG4Particle * p3 = (AliAODPWG4Particle *)array3->At(ke); //third photon\r
1039                 TLorentzVector photon1(p1->Px(),p1->Py(),p1->Pz(),p1->E());\r
1040                 TLorentzVector photon2(p2->Px(),p2->Py(),p2->Pz(),p2->E());\r
1041                 Double_t px12 =(photon1+photon2).Px();\r
1042                 Double_t py12 =(photon1+photon2).Py();\r
1043                 Double_t pz12 =(photon1+photon2).Pz();\r
1044                 Double_t e12 =(photon1+photon2).E();\r
1045                 TLorentzVector photon12(px12, py12,pz12,e12);\r
1046                 TLorentzVector photon3(p3->Px(),p3->Py(),p3->Pz(),p3->E());\r
1047                 GetAsyDistPtM(photon12, photon3,fAsy,fDist,fPt,fMass);\r
1048 \r
1049                 for(Int_t ipid=0;ipid<fNPID;ipid++){\r
1050                     if(p1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
1051                             p2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
1052                             p3->IsPIDOK(ipid,AliCaloPID::kPhoton)){\r
1053                        fMixCPi0GammaAsyPtM[ipid]->Fill(fAsy, fPt, fMass); //default module, with pid\r
1054                     }\r
1055                     for(Int_t jmod=0;jmod<fNMod;jmod++){\r
1056                         for(Int_t kasy=0;kasy<fNAsy;kasy++){\r
1057                            if(IsPhotonSelected(p1,ipid,fModCut[jmod]) && \r
1058                               IsPhotonSelected(p2,ipid,fModCut[jmod])&&\r
1059                               IsPhotonSelected(p3,ipid,fModCut[jmod])){ \r
1060                                \r
1061                               //filling the histo\r
1062                                for(Int_t mhist=0;mhist<fNHistos;mhist++){\r
1063                                    Double_t pt1 =mhist *fPerPtBin;\r
1064                                    Double_t pt2 =(mhist+1)*fPerPtBin;\r
1065                                    Int_t index = ipid*fNMod*fNAsy*fNHistos+                                                                                        jmod*fNAsy*fNHistos+kasy*fNHistos+mhist;\r
1066                                    if(fMass<fInvMassCut && f2gammaAsy<fAsyCut[kasy] && \r
1067                                       fPt<pt2 && fPt>pt1&& fAsy<fAsyCut[kasy] &&fAsy&&fDist&&fPt&&fMass)\r
1068                                    fMixOmegaC[index]->Fill(fMass);\r
1069                               }//hist\r
1070                             }//photon selection\r
1071                         }//asy\r
1072                     }//mod\r
1073                 }//pid\r
1074             }//ke\r
1075         }//pi0 candidate\r
1076      } //je\r
1077  }//ie\r
1078 \r
1079 } \r
1080  \r
1081 //______________________________________________________________________________\r
1082 void AliAnaNeutralMeson::PhotonAcceptance(AliStack * stack)\r
1083 {\r
1084  //photon Acceptance\r
1085  for(Int_t i=0 ; i<stack->GetNtrack(); i++){\r
1086      TParticle * prim = stack->Particle(i) ;\r
1087      Int_t pdg = prim->GetPdgCode() ;\r
1088 //   Int_t ndau = prim->GetNDaughters();\r
1089 //   Int_t iphot1=prim->GetFirstDaughter() ;\r
1090 //   Int_t iphot2=prim->GetLastDaughter() ;\r
1091      if(pdg ==22 && stack->IsPhysicalPrimary(i)) {\r
1092           Double_t photonPt = prim->Pt() ;\r
1093           //printf("photon, pt %2.2f\n",photonPt);                           \r
1094           Double_t photonY  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;   \r
1095           Double_t phi   = TMath::RadToDeg()*prim->Phi() ;\r
1096           if(TMath::Abs(photonY) < 0.5){  fhPrimPhotonPt->Fill(photonPt);  }\r
1097               fhPrimPhotonY  ->Fill(photonY) ;\r
1098               fhPrimPhotonPhi->Fill(phi) ;\r
1099               //Check if both photons hit Calorimeter\r
1100               Bool_t inacceptance = kFALSE;\r
1101 \r
1102         if(fCalorimeter == "PHOS"){\r
1103 \r
1104 #ifdef __PHOSUTIL__\r
1105             Int_t mod ;\r
1106             Double_t x,z ;\r
1107             if(fPHOSGeo->ImpactOnEmc(prim,mod,z,x)) \r
1108               inacceptance = kTRUE;\r
1109             //printf("In REAL PHOS acceptance? %d\n",inacceptance);\r
1110 #else\r
1111             TLorentzVector lv1;\r
1112             prim->Momentum(lv1);\r
1113             if(GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter)) \r
1114               inacceptance = kTRUE ;\r
1115              if(GetDebug() > 2) printf("In %s fidutial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);\r
1116 #endif                                                                                                    \r
1117 \r
1118 }          \r
1119 \r
1120         else if(fCalorimeter == "EMCAL"){\r
1121 \r
1122 #ifdef __EMCALUTIL__\r
1123 \r
1124             if(fEMCALGeo->Impact(prim)) \r
1125               inacceptance = kTRUE;\r
1126             if(GetDebug() > 2) printf("In REAL EMCAL acceptance? %d\n",inacceptance);\r
1127 #else\r
1128             TLorentzVector lv1;\r
1129             prim->Momentum(lv1);\r
1130             if(GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter)) \r
1131               inacceptance = kTRUE ;\r
1132             //printf("In %s fidutial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);\r
1133 #endif                                                                                                    \r
1134         }         \r
1135 \r
1136               if(inacceptance){\r
1137                   fhPrimPhotonAccPt->Fill(photonPt) ;\r
1138                   fhPrimPhotonAccPhi->Fill(phi) ;\r
1139                   fhPrimPhotonAccY->Fill(photonY) ;\r
1140                }//Accepted\r
1141  \r
1142   }// Primary photon\r
1143  }//loop on primaries  \r
1144 }\r
1145  \r
1146 //______________________________________________________________________________\r
1147 void AliAnaNeutralMeson::Pi0EtaAcceptance(AliStack * stack)\r
1148 {\r
1149  //pi0 and eta Acceptance\r
1150  for(Int_t i=0 ; i<stack->GetNprimary(); i++){\r
1151      TParticle * prim = stack->Particle(i) ;\r
1152      Int_t pdg = prim->GetPdgCode() ;\r
1153      if(pdg==111 || pdg==221){\r
1154         Double_t pi0Pt = prim->Pt() ;\r
1155         //printf("  pi0, pt %2.2f\n",pi0Pt);\r
1156         Double_t pi0Y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;\r
1157         Double_t phi   = TMath::RadToDeg()*prim->Phi();\r
1158         if(pdg==111){\r
1159            if(TMath::Abs(pi0Y) < 0.5)   fhPrimPi0Pt->Fill(pi0Pt) ;\r
1160               fhPrimPi0Y  ->Fill(pi0Y) ;\r
1161               fhPrimPi0Phi->Fill(phi) ;\r
1162          } \r
1163         else{\r
1164            if(TMath::Abs(pi0Y) < 0.5)  fhPrimEtaPt->Fill(pi0Pt) ;\r
1165               fhPrimEtaY  ->Fill(pi0Y) ;\r
1166               fhPrimEtaPhi->Fill(phi) ;\r
1167         }\r
1168         //Check if both photons hit Calorimeter\r
1169         Int_t ndau = prim->GetNDaughters();\r
1170         Int_t iphot1=prim->GetFirstDaughter() ;\r
1171         Int_t iphot2=prim->GetLastDaughter() ;\r
1172         if(ndau==2 && ndau>1 && iphot1>-1 && iphot1<stack->GetNtrack() \r
1173                    && iphot2>-1 && iphot2<stack->GetNtrack()){   \r
1174             TParticle * phot1 = stack->Particle(iphot1) ;\r
1175             TParticle * phot2 = stack->Particle(iphot2) ;\r
1176             if(phot1 && phot2 && phot1->GetPdgCode()==22 && phot2->GetPdgCode()==22){\r
1177                 Bool_t inacceptance = kFALSE;              \r
1178 \r
1179                 if(fCalorimeter == "PHOS"){\r
1180   \r
1181 #ifdef __PHOSUTIL__\r
1182                  Int_t mod, mod1,mod2 ;        \r
1183                  Double_t x, z, x1,z1,x2,z2;\r
1184                  if(fPHOSGeo->ImpactOnEmc(prim, mod, z,x)  && fPHOSGeo->ImpactOnEmc(phot1,mod1,z1,x1) \r
1185                           && fPHOSGeo->ImpactOnEmc(phot2,mod2,z2,x2))  inacceptance = kTRUE;\r
1186                                     \r
1187                  //printf("In REAL PHOS acceptance? %d\n",inacceptance);\r
1188 #else \r
1189                  TLorentzVector lv0, lv1, lv3;\r
1190                  prim->Momentum(lv0);\r
1191                  phot1->Momentum(lv1);\r
1192                  phot2->Momentum(lv3);\r
1193                  if(GetFidutialCut()->IsInFidutialCut(lv0,fCalorimeter) \r
1194                     &&GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) \r
1195                     && GetFidutialCut()->IsInFidutialCut(lv3,fCalorimeter))  inacceptance = kTRUE ;\r
1196 #endif\r
1197         }\r
1198 \r
1199 \r
1200                 else if(fCalorimeter == "EMCAL"){\r
1201   \r
1202 #ifdef __EMCALUTIL__\r
1203                  if(fEMCALGeo->Impact(prim)  && fEMCALGeo->Impact(phot1) \r
1204                           && fEMCALGeo->Impact(phot2))  inacceptance = kTRUE;\r
1205                                     \r
1206                  //printf("In REAL EMCAL acceptance? %d\n",inacceptance);\r
1207 #else \r
1208                  TLorentzVector lv0, lv1, lv3;\r
1209                  prim->Momentum(lv0);\r
1210                  phot1->Momentum(lv1);\r
1211                  phot2->Momentum(lv3);\r
1212                  if(GetFidutialCut()->IsInFidutialCut(lv0,fCalorimeter) \r
1213                     &&GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) \r
1214                     && GetFidutialCut()->IsInFidutialCut(lv3,fCalorimeter))  inacceptance = kTRUE ;\r
1215 #endif\r
1216         }\r
1217 \r
1218 //                printf("In %s fFidutial cut acceptance? %d\n",fCalorimeter.Data(),inacceptance);\r
1219                  if(inacceptance){\r
1220                     if(pdg==111){\r
1221                         fhPrimPi0AccPt->Fill(pi0Pt) ;\r
1222                         fhPrimPi0AccPhi->Fill(phi) ;                                                                                      fhPrimPi0AccY->Fill(pi0Y) ;\r
1223                      }\r
1224                      else{\r
1225                         fhPrimEtaAccPt->Fill(pi0Pt) ;                                                                                     fhPrimEtaAccPhi->Fill(phi) ;\r
1226                         fhPrimEtaAccY->Fill(pi0Y) ;\r
1227                     }\r
1228                   }//Accepted\r
1229            }// 2 photons      \r
1230        }//Check daughters exist\r
1231     }// Primary pi0 \r
1232   }//loop on primaries    \r
1233 }\r
1234  \r
1235 //______________________________________________________________________________\r
1236 void AliAnaNeutralMeson::OmegaAcceptance(AliStack * stack)\r
1237 {\r
1238  //acceptance for omega->pi0+gamma       \r
1239  for(Int_t i=0 ; i<stack->GetNprimary(); i++){ //aaa\r
1240      TParticle * prim = stack->Particle(i) ;\r
1241      Int_t pdg = prim->GetPdgCode();\r
1242      if(pdg == 223){ //bbb\r
1243         Double_t pt = prim->Pt() ;\r
1244         Double_t y  = 0.5*TMath::Log((prim->Energy()-prim->Pz())/(prim->Energy()+prim->Pz())) ;\r
1245         Double_t phi   = TMath::RadToDeg()*prim->Phi() ;\r
1246         if(TMath::Abs(y) < 0.5){  fhPrimOmegaPt->Fill(pt) ;  }\r
1247         fhPrimOmegaY  ->Fill(y) ;  \r
1248         fhPrimOmegaPhi->Fill(phi) ;\r
1249         //check if omega->pi0+gamma->(gamma+gamma)+gamma hit Calorimeter\r
1250         Int_t ndau = prim->GetNDaughters();\r
1251         Int_t iphot1=prim->GetFirstDaughter() ;\r
1252         Int_t iphot2=prim->GetLastDaughter() ;\r
1253         if(ndau==2 && iphot1>-1 && iphot1<stack->GetNtrack() && iphot2>-1 && iphot2<stack->GetNtrack()){ //ccc\r
1254            TParticle * phot1 = stack->Particle(iphot1) ;\r
1255            TParticle * phot2 = stack->Particle(iphot2) ;\r
1256            Int_t pdg1 = phot1->GetPdgCode();\r
1257            Int_t pdg2 = phot2->GetPdgCode();        \r
1258            Int_t ndau2 = phot2->GetNDaughters();\r
1259            if(ndau2==2 && phot1 && pdg1==22 &&phot2 && pdg2==111){ //ddd\r
1260               Int_t jphot1 = phot2->GetFirstDaughter() ;\r
1261               Int_t jphot2 = phot2->GetLastDaughter() ;\r
1262               if(jphot1>-1 && jphot1<stack->GetNtrack() &&\r
1263                               jphot2>-1 && jphot2<stack->GetNtrack()){  //eee                \r
1264                   TParticle * phot3 = stack->Particle(jphot1) ;\r
1265                   TParticle * phot4 = stack->Particle(jphot2) ;\r
1266                   Int_t pdg3 = phot3->GetPdgCode();\r
1267                   Int_t pdg4 = phot4->GetPdgCode();\r
1268                   if(phot3 && pdg3==22 && phot4 && pdg4==22){ //fff\r
1269                      Bool_t inacceptance = kFALSE;\r
1270 \r
1271         if(fCalorimeter == "PHOS"){\r
1272 #ifdef __PHOSUTIL__\r
1273                      Int_t mod, mod1,mod2,mod3 ;\r
1274                      Double_t x,z, x1,z1,x2,z2, x3,z3 ;\r
1275                      //3 gammas hit Calorimeter\r
1276                      if(fPHOSGeo->ImpactOnEmc(prim,mod,z,x) \r
1277                                                && fPHOSGeo->ImpactOnEmc(phot1,mod1,z1,x1) \r
1278                                                && fPHOSGeo->ImpactOnEmc(phot3,mod2,z2,x2) \r
1279                                                && fPHOSGeo->ImpactOnEmc(phot4,mod3,z3,x3))\r
1280                       inacceptance = kTRUE;\r
1281  \r
1282 #else\r
1283                       TLorentzVector lv0,lv1, lv2, lv3;                                                         \r
1284                       prim->Momentum(lv0);\r
1285                       phot1->Momentum(lv1);\r
1286                       phot3->Momentum(lv2);\r
1287                       phot4->Momentum(lv3);\r
1288                       if(GetFidutialCut()->IsInFidutialCut(lv0,fCalorimeter) &&\r
1289                          GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) && \r
1290                          GetFidutialCut()->IsInFidutialCut(lv2,fCalorimeter) && \r
1291                          GetFidutialCut()->IsInFidutialCut(lv3,fCalorimeter))  \r
1292                       inacceptance = kTRUE ;\r
1293 #endif\r
1294         }\r
1295 \r
1296         else if(fCalorimeter == "EMCAL"){\r
1297 #ifdef __EMCALUTIL__\r
1298                      //3 gammas hit Calorimeter\r
1299                      if(fEMCALGeo->Impact(prim) \r
1300                                                && fEMCALGeo->Impact(phot1) \r
1301                                                && fEMCALGeo->Impact(phot3) \r
1302                                                && fEMCALGeo->Impact(phot4))\r
1303                       inacceptance = kTRUE;\r
1304  \r
1305 #else\r
1306                       TLorentzVector lv0,lv1, lv2, lv3;                                                         \r
1307                       prim->Momentum(lv0);\r
1308                       phot1->Momentum(lv1);\r
1309                       phot3->Momentum(lv2);\r
1310                       phot4->Momentum(lv3);\r
1311                       if(GetFidutialCut()->IsInFidutialCut(lv0,fCalorimeter) &&\r
1312                          GetFidutialCut()->IsInFidutialCut(lv1,fCalorimeter) && \r
1313                          GetFidutialCut()->IsInFidutialCut(lv2,fCalorimeter) && \r
1314                          GetFidutialCut()->IsInFidutialCut(lv3,fCalorimeter))  \r
1315                       inacceptance = kTRUE ;\r
1316 #endif\r
1317              }\r
1318 \r
1319                       if(inacceptance){ //ggg\r
1320                          fhPrimOmegaAccPt->Fill(pt) ;\r
1321                          fhPrimOmegaAccPhi->Fill(phi) ;\r
1322                          fhPrimOmegaAccY->Fill(y) ;\r
1323                        } //ggg\r
1324                    }//fff\r
1325                 }//eee\r
1326              }//ddd\r
1327          }//ccc\r
1328        } //bbb\r
1329   }//aaa                                  \r
1330 }                                                               \r
1331  \r
1332 void AliAnaNeutralMeson::GetMCDistAsy(TParticle *p1,TParticle *p2,Double_t &dist,Double_t &asy)\r
1333 {   \r
1334   //get dist and asy for MC particle                                                            \r
1335   Double_t eta1 = p1->Eta();        \r
1336   Double_t eta2 = p2->Eta(); \r
1337   Double_t phi1 = p1->Phi();\r
1338   Double_t phi2 = p2->Phi();                                            \r
1339   Double_t e1 = p1->Energy(); \r
1340   Double_t e2 = p2->Energy();                                                               \r
1341   dist = TMath::Sqrt( (eta1-eta2)*(eta1-eta2) + (phi1-phi2)*(phi1-phi2));\r
1342   if(e1+e2) asy = TMath::Abs(e1-e2)/(e1+e2);\r
1343 }                                   \r
1344 \r
1345 //______________________________________________________________________________\r
1346 //void AliAnaNeutralMeson::Terminate(TList * outList) \r
1347 //{\r
1348 // //Do some calculations and plots from the final histograms.\r
1349 //  printf("AliAnaNeutralMeson::Terminate() - Not implemented yet\n");\r
1350 //\r
1351 //}\r
1352 \r
1353 \r