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