-/**************************************************************************\r
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- * *\r
- * Author: The ALICE Off-line Project. *\r
- * Contributors are mentioned in the code where appropriate. *\r
- * *\r
- * Permission to use, copy, modify and distribute this software and its *\r
- * documentation strictly for non-commercial purposes is hereby granted *\r
- * without fee, provided that the above copyright notice appears in all *\r
- * copies and that both the copyright notice and this permission notice *\r
- * appear in the supporting documentation. The authors make no claims *\r
- * about the suitability of this software for any purpose. It is *\r
- * provided "as is" without express or implied warranty. *\r
- **************************************************************************/\r
-/* $Id: $ */\r
-//_________________________________________________________________________\r
-// class to extract omega(782)->pi0+gamma->3gamma\r
-// Mar. 22, 2011: Additional method, espeically for EMCAL. A high E cluster is assumpted as pi0 (two photons are overlapped) without unfolding\r
-//\r
-//-- Author: Renzhuo Wan (IOPP-Wuhan, China)\r
-//_________________________________________________________________________\r
-\r
-// --- ROOT system\r
-class TROOT;\r
-\r
-// --- AliRoot system\r
-//class AliVEvent;\r
-// --- ROOT system ---\r
-#include "TH2F.h"\r
-#include "TLorentzVector.h"\r
-#include "TParticle.h"\r
-#include "TCanvas.h"\r
-#include "TFile.h"\r
-//---- AliRoot system ----\r
-#include "AliAnaOmegaToPi0Gamma.h"\r
-#include "AliCaloTrackReader.h"\r
-#include "AliCaloPID.h"\r
-#include "AliStack.h"\r
-#include "AliVEvent.h"\r
-#include "AliAODEvent.h"\r
-#include "AliAODMCParticle.h"\r
-ClassImp(AliAnaOmegaToPi0Gamma)\r
-\r
-//______________________________________________________________________________\r
-AliAnaOmegaToPi0Gamma::AliAnaOmegaToPi0Gamma() : AliAnaPartCorrBaseClass(),\r
-fInputAODPi0(0), fInputAODGammaName(""),\r
-fEventsList(0x0),fNVtxZBin(0), fNCentBin(0), fNRpBin(0), fNBadChDistBin(0), fNpid(0),\r
-fNmaxMixEv(0), fVtxZCut(0), fCent(0), fRp(0), \r
-fPi0Mass(0),fPi0MassWindow(0),fPi0OverOmegaPtCut(0),\r
-fGammaOverOmegaPtCut(0), fEOverlapCluster(0),\r
-fhEtalon(0),\r
-fRealOmega0(0), fMixAOmega0(0),\r
-fMixBOmega0(0), fMixCOmega0(0),\r
-fRealOmega1(0), fMixAOmega1(0),\r
-fMixBOmega1(0), fMixCOmega1(0),\r
-fRealOmega2(0), fMixAOmega2(0),\r
-fMixBOmega2(0), fMixCOmega2(0),\r
-fhFakeOmega(0),\r
-fhOmegaPriPt(0)\r
-{\r
- //Default Ctor\r
- InitParameters();\r
-}\r
-/*\r
-//______________________________________________________________________________\r
-AliAnaOmegaToPi0Gamma::AliAnaOmegaToPi0Gamma(const AliAnaOmegaToPi0Gamma & ex) : AliAnaPartCorrBaseClass(ex),\r
-fInputAODPi0(new TClonesArray (*ex.fInputAODPi0)),\r
-fInputAODGammaName(ex.fInputAODGammaName),\r
-fEventsList(0x0), \r
-fNVtxZBin(ex.fNVtxZBin), fNCentBin(ex.fNCentBin), fNRpBin(ex.fNRpBin),\r
-fNBadChDistBin(ex.fNBadChDistBin),fNpid(ex.fNpid),\r
-fNmaxMixEv(ex.fNmaxMixEv),\r
-fVtxZCut(ex.fVtxZCut), fCent(ex.fCent), fRp(ex.fRp), \r
-fPi0Mass(ex.fPi0Mass),\r
-fPi0MassWindow(ex.fPi0MassWindow),\r
-fPi0OverOmegaPtCut(ex.fPi0OverOmegaPtCut),\r
-fGammaOverOmegaPtCut(ex.fGammaOverOmegaPtCut),\r
-fhEtalon(ex.fhEtalon),\r
-fRealOmega0(ex.fRealOmega0), fMixAOmega0(ex.fMixAOmega0),\r
-fMixBOmega0(ex.fMixBOmega0), fMixCOmega0(ex.fMixCOmega0),\r
-fRealOmega1(ex.fRealOmega1), fMixAOmega1(ex.fMixAOmega1),\r
-fMixBOmega1(ex.fMixBOmega1), fMixCOmega1(ex.fMixCOmega1),\r
-fRealOmega2(ex.fRealOmega2), fMixAOmega2(ex.fMixAOmega2),\r
-fMixBOmega2(ex.fMixBOmega2), fMixCOmega2(ex.fMixCOmega2),\r
-fhOmegaPriPt(ex.fhOmegaPriPt)\r
-{\r
- // cpy ctor\r
- //Do not need it\r
-}\r
-*/\r
-/*\r
-//______________________________________________________________________________\r
-AliAnaOmegaToPi0Gamma & AliAnaOmegaToPi0Gamma::operator = (const AliAnaOmegaToPi0Gamma & ex)\r
-{\r
- // assignment operator\r
- \r
- if(this == &ex)return *this;\r
- ((AliAnaPartCorrBaseClass *)this)->operator=(ex);\r
- fInputAODPi0 = new TClonesArray(*ex.fInputAODPi0);\r
- fInputAODGammaName = ex.fInputAODGammaName;\r
- fEventsList = ex.fEventsList;\r
-\r
- fNVtxZBin=ex.fNVtxZBin;\r
- fNCentBin=ex.fNCentBin;\r
- fNRpBin=ex.fNRpBin;\r
- fNBadChDistBin=ex.fNBadChDistBin;\r
- fNpid=ex.fNpid;\r
- fNmaxMixEv =ex.fNmaxMixEv;\r
-\r
- fVtxZCut=ex.fVtxZCut;\r
- fCent=ex.fCent;\r
- fRp=ex.fRp;\r
-\r
- fPi0Mass=ex.fPi0Mass;\r
- fPi0MassWindow=ex.fPi0MassWindow;\r
- fPi0OverOmegaPtCut=ex.fPi0OverOmegaPtCut;\r
- fGammaOverOmegaPtCut=ex.fGammaOverOmegaPtCut;\r
-\r
- fhEtalon=ex.fhEtalon;\r
- fRealOmega0=ex.fRealOmega0;\r
- fMixAOmega0=ex.fMixAOmega0;\r
- fMixBOmega0=ex.fMixBOmega0;\r
- fMixCOmega0=ex.fMixCOmega0;\r
- fRealOmega1=ex.fRealOmega1;\r
- fMixAOmega1=ex.fMixAOmega1;\r
- fMixBOmega1=ex.fMixBOmega1;\r
- fMixCOmega1=ex.fMixCOmega1;\r
- fRealOmega2=ex.fRealOmega2;\r
- fMixAOmega2=ex.fMixAOmega2;\r
- fMixBOmega2=ex.fMixBOmega2;\r
- fMixCOmega2=ex.fMixCOmega2;\r
- fhOmegaPriPt=ex.fhOmegaPriPt;\r
- return *this;\r
- \r
-}\r
-*/\r
-\r
-//______________________________________________________________________________\r
-AliAnaOmegaToPi0Gamma::~AliAnaOmegaToPi0Gamma() {\r
-\r
- //dtor\r
-// Done by the maker \r
-// if(fInputAODPi0){\r
-// fInputAODPi0->Clear();\r
-// delete fInputAODPi0;\r
-// } \r
-\r
- if(fEventsList){\r
- for(Int_t i=0;i<fNVtxZBin;i++){\r
- for(Int_t j=0;j<fNCentBin;j++){\r
- for(Int_t k=0;k<fNRpBin;k++){\r
- fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k]->Clear();\r
- delete fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k];\r
- }\r
- }\r
- }\r
- }\r
- delete [] fEventsList;\r
- fEventsList=0;\r
-\r
- delete [] fVtxZCut;\r
- delete [] fCent;\r
- delete [] fRp;\r
-\r
-}\r
-\r
-//______________________________________________________________________________\r
-void AliAnaOmegaToPi0Gamma::InitParameters()\r
-{\r
-//Init parameters when first called the analysis\r
-//Set default parameters\r
- fInputAODGammaName="PhotonsDetector"; \r
- fNVtxZBin=1; \r
- fNCentBin=1; \r
- fNRpBin=1; \r
- fNBadChDistBin=3; \r
- fNpid=1; \r
- fNmaxMixEv=8; \r
- \r
- fPi0Mass=0.1348; \r
- fPi0MassWindow=0.015; \r
- fPi0OverOmegaPtCut=0.8; \r
- fGammaOverOmegaPtCut=0.2; \r
- \r
- fEOverlapCluster=6;\r
-}\r
-\r
-\r
-//______________________________________________________________________________\r
-TList * AliAnaOmegaToPi0Gamma::GetCreateOutputObjects()\r
-{\r
- //\r
- fVtxZCut = new Double_t [fNVtxZBin];\r
- for(Int_t i=0;i<fNVtxZBin;i++) fVtxZCut[i]=10*(i+1);\r
- \r
- fCent=new Double_t[fNCentBin];\r
- for(int i = 0;i<fNCentBin;i++)fCent[i]=0;\r
- \r
- fRp=new Double_t[fNRpBin];\r
- for(int i = 0;i<fNRpBin;i++)fRp[i]=0;\r
- //\r
- Int_t nptbins = GetHistoPtBins();\r
- Float_t ptmax = GetHistoPtMax();\r
- Float_t ptmin = GetHistoPtMin();\r
- \r
- Int_t nmassbins = GetHistoMassBins();\r
- Float_t massmin = GetHistoMassMin();\r
- Float_t massmax = GetHistoMassMax();\r
- \r
- fhEtalon = new TH2F("hEtalon","Histo with binning parameters", nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;\r
- fhEtalon->SetXTitle("P_{T} (GeV)") ;\r
- fhEtalon->SetYTitle("m_{inv} (GeV)") ;\r
- \r
- // store them in fOutputContainer\r
- fEventsList = new TList*[fNVtxZBin*fNCentBin*fNRpBin];\r
- for(Int_t i=0;i<fNVtxZBin;i++){\r
- for(Int_t j=0;j<fNCentBin;j++){\r
- for(Int_t k=0;k<fNRpBin;k++){\r
- fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k] = new TList();\r
- fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k]->SetOwner(kFALSE); \r
- }\r
- }\r
- }\r
- \r
- TList * outputContainer = new TList() ; \r
- outputContainer->SetName(GetName());\r
- const Int_t buffersize = 255;\r
- char key[buffersize] ;\r
- char title[buffersize] ;\r
- const char * detector= fInputAODGammaName.Data();\r
- Int_t ndim=fNVtxZBin*fNCentBin*fNRpBin*fNBadChDistBin*fNpid;\r
- \r
- fRealOmega0 =new TH2F*[ndim];\r
- fMixAOmega0 =new TH2F*[ndim];\r
- fMixBOmega0 =new TH2F*[ndim];\r
- fMixCOmega0 =new TH2F*[ndim];\r
- \r
- fRealOmega1 =new TH2F*[ndim];\r
- fMixAOmega1 =new TH2F*[ndim];\r
- fMixBOmega1 =new TH2F*[ndim];\r
- fMixCOmega1 =new TH2F*[ndim];\r
- \r
- fRealOmega2 =new TH2F*[ndim];\r
- fMixAOmega2 =new TH2F*[ndim];\r
- fMixBOmega2 =new TH2F*[ndim];\r
- fMixCOmega2 =new TH2F*[ndim];\r
- \r
- fhFakeOmega = new TH2F*[fNCentBin];\r
- \r
- for(Int_t i=0;i<fNVtxZBin;i++){\r
- for(Int_t j=0;j<fNCentBin;j++){\r
- for(Int_t k=0;k<fNRpBin;k++){ //at event level\r
- Int_t idim=i*fNCentBin*fNRpBin+j*fNRpBin+k;\r
- for(Int_t ipid=0;ipid<fNpid;ipid++){\r
- for(Int_t idist=0;idist<fNBadChDistBin;idist++){ //at particle level\r
- \r
- Int_t index=idim*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;\r
- \r
- snprintf(key,buffersize,"RealToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s Real Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fRealOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fRealOmega0[index]->SetName(key) ;\r
- fRealOmega0[index]->SetTitle(title);\r
- outputContainer->Add(fRealOmega0[index]);\r
- \r
- snprintf(key,buffersize,"MixAToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s MixA Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fMixAOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fMixAOmega0[index]->SetName(key) ;\r
- fMixAOmega0[index]->SetTitle(title);\r
- outputContainer->Add(fMixAOmega0[index]);\r
- \r
- snprintf(key,buffersize,"MixBToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s MixB Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fMixBOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fMixBOmega0[index]->SetName(key) ;\r
- fMixBOmega0[index]->SetTitle(title);\r
- outputContainer->Add(fMixBOmega0[index]);\r
- \r
- snprintf(key,buffersize,"MixCToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s MixC Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fMixCOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fMixCOmega0[index]->SetName(key) ;\r
- fMixCOmega0[index]->SetTitle(title);\r
- outputContainer->Add(fMixCOmega0[index]);\r
- \r
- snprintf(key,buffersize,"RealToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s Real Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fRealOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fRealOmega1[index]->SetName(key) ;\r
- fRealOmega1[index]->SetTitle(title);\r
- outputContainer->Add(fRealOmega1[index]);\r
- \r
- snprintf(key,buffersize,"MixAToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s MixA Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fMixAOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fMixAOmega1[index]->SetName(key) ;\r
- fMixAOmega1[index]->SetTitle(title);\r
- outputContainer->Add(fMixAOmega1[index]);\r
- \r
- snprintf(key,buffersize,"MixBToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s MixB Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fMixBOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fMixBOmega1[index]->SetName(key) ;\r
- fMixBOmega1[index]->SetTitle(title);\r
- outputContainer->Add(fMixBOmega1[index]);\r
- \r
- snprintf(key,buffersize,"MixCToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s MixC Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fMixCOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fMixCOmega1[index]->SetName(key) ;\r
- fMixCOmega1[index]->SetTitle(title);\r
- outputContainer->Add(fMixCOmega1[index]);\r
- \r
- snprintf(key,buffersize,"RealToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s Real Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fRealOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fRealOmega2[index]->SetName(key) ;\r
- fRealOmega2[index]->SetTitle(title);\r
- outputContainer->Add(fRealOmega2[index]);\r
- \r
- snprintf(key,buffersize,"MixAToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s MixA Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fMixAOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fMixAOmega2[index]->SetName(key) ;\r
- fMixAOmega2[index]->SetTitle(title);\r
- outputContainer->Add(fMixAOmega2[index]);\r
- \r
- snprintf(key,buffersize,"MixBToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s MixB Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fMixBOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fMixBOmega2[index]->SetName(key) ;\r
- fMixBOmega2[index]->SetTitle(title);\r
- outputContainer->Add(fMixBOmega2[index]);\r
- \r
- snprintf(key,buffersize,"MixCToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);\r
- snprintf(title,buffersize, "%s MixC Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);\r
- fMixCOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;\r
- fMixCOmega2[index]->SetName(key) ;\r
- fMixCOmega2[index]->SetTitle(title);\r
- outputContainer->Add(fMixCOmega2[index]);\r
- }\r
- }\r
- }\r
- } \r
- }\r
-\r
- for(Int_t i=0;i<fNCentBin;i++){\r
- snprintf(key,buffersize, "fhFakeOmega%d",i);\r
- snprintf(title,buffersize,"FakePi0(high pt cluster)+Gamma with Centrality%d",i);\r
- fhFakeOmega[i] = (TH2F*)fhEtalon->Clone(key);\r
- fhFakeOmega[i]->SetTitle(title);\r
- outputContainer->Add(fhFakeOmega[i]); \r
- }\r
- if(IsDataMC()){\r
- snprintf(key,buffersize, "%sOmegaPri",detector);\r
- snprintf(title,buffersize,"primary #omega in %s",detector);\r
- fhOmegaPriPt=new TH1F(key, title,nptbins,ptmin,ptmax);\r
- fhOmegaPriPt->GetXaxis()->SetTitle("P_{T}");\r
- fhOmegaPriPt->GetYaxis()->SetTitle("dN/P_{T}");\r
- outputContainer->Add(fhOmegaPriPt);\r
- }\r
- \r
- delete fhEtalon;\r
- return outputContainer;\r
-}\r
-\r
-//______________________________________________________________________________\r
-void AliAnaOmegaToPi0Gamma::Print(const Option_t * /*opt*/) const\r
-{\r
- //Print some relevant parameters set in the analysis\r
- printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;\r
- AliAnaPartCorrBaseClass::Print(" ");\r
- printf("Omega->pi0+gamma->3gamma\n");\r
- printf("Cuts at event level: \n");\r
- printf("Bins of vertex Z: %d \n", fNVtxZBin);\r
- printf("Bins of centrality: %d \n",fNCentBin);\r
- printf("Bins of Reaction plane: %d\n",fNRpBin);\r
- printf("Cuts at AOD particle level:\n");\r
- printf("Number of PID: %d \n", fNpid);\r
- printf("Number of DistToBadChannel cuts: %d\n", fNBadChDistBin);\r
- printf("number of events buffer to be mixed: %d\n",fNmaxMixEv);\r
-} \r
-\r
-//______________________________________________________________________________\r
-void AliAnaOmegaToPi0Gamma::MakeAnalysisFillHistograms() \r
-{\r
- //fill the MC AOD if needed first\r
- //-----------\r
- //need to be further implemented\r
- AliStack * stack = 0x0;\r
- // TParticle * primary = 0x0;\r
- TClonesArray * mcparticles0 = 0x0;\r
- //TClonesArray * mcparticles1 = 0x0;\r
- AliAODMCParticle * aodprimary = 0x0;\r
- Int_t pdg=0;\r
- Double_t pt=0;\r
- Double_t eta=0;\r
- \r
- if(IsDataMC()){\r
- if(GetReader()->ReadStack()){\r
- stack = GetMCStack() ;\r
- if(!stack){\r
- printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - There is no stack!\n");\r
- }\r
- else{\r
- for(Int_t i=0 ; i<stack->GetNtrack(); i++){\r
- TParticle * prim = stack->Particle(i) ;\r
- pdg = prim->GetPdgCode() ;\r
- eta=prim->Eta();\r
- pt=prim->Pt();\r
- if(TMath::Abs(eta)<0.5) {\r
- if(pdg==223) fhOmegaPriPt->Fill(pt);\r
- }\r
- }\r
- }\r
- }\r
- else if(GetReader()->ReadAODMCParticles()){\r
- //Get the list of MC particles\r
- mcparticles0 = GetReader()->GetAODMCParticles(0);\r
- if(!mcparticles0 ) {\r
- if(GetDebug() > 0) printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");\r
- }\r
- // if(GetReader()->GetSecondInputAODTree()){\r
- // mcparticles1 = GetReader()->GetAODMCParticles(1);\r
- // if(!mcparticles1 && GetDebug() > 0) {\r
- // printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");\r
- // }\r
- // }\r
- else{\r
- for(Int_t i=0;i<mcparticles0->GetEntries();i++){\r
- aodprimary =(AliAODMCParticle*)mcparticles0->At(i);\r
- pdg = aodprimary->GetPdgCode() ;\r
- eta=aodprimary->Eta();\r
- pt=aodprimary->Pt();\r
- if(TMath::Abs(eta)<0.5) {\r
- if(pdg==223) fhOmegaPriPt->Fill(pt);\r
- }\r
- \r
- }\r
- }//mcparticles0 exists\r
- }//AOD MC Particles\r
- }// is data and MC\r
- \r
- \r
- //process event from AOD brach \r
- //extract pi0, eta and omega analysis\r
- Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;\r
- if(IsBadRun(iRun)) return ; \r
- \r
- //vertex z\r
- Double_t vert[]={0,0,0} ;\r
- GetVertex(vert);\r
- Int_t curEventBin =0;\r
- \r
- Int_t ivtxzbin=(Int_t)TMath::Abs(vert[2])/10;\r
- if(ivtxzbin>=fNVtxZBin)return;\r
- \r
- //centrality\r
- Int_t currentCentrality = GetEventCentrality();\r
- if(currentCentrality == -1) return;\r
- Int_t optCent = GetReader()->GetCentralityOpt();\r
- Int_t icentbin=currentCentrality/(optCent/fNCentBin) ; //GetEventCentrality();\r
- \r
- printf("-------------- %d %d %d ",currentCentrality, optCent, icentbin);\r
- //reaction plane\r
- Int_t irpbin=0;\r
- \r
- if(ivtxzbin==-1) return; \r
- curEventBin = ivtxzbin*fNCentBin*fNRpBin + icentbin*fNRpBin + irpbin;\r
- TClonesArray *aodGamma = (TClonesArray*) GetAODBranch(fInputAODGammaName); //photon array\r
- // TClonesArray *aodGamma = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaName); //photon array\r
- Int_t nphotons =0;\r
- if(aodGamma) nphotons= aodGamma->GetEntries(); \r
- else return;\r
- \r
- fInputAODPi0 = (TClonesArray*)GetInputAODBranch(); //pi0 array\r
- Int_t npi0s = 0;\r
- if(fInputAODPi0) npi0s= fInputAODPi0 ->GetEntries();\r
- else return;\r
- \r
- if(nphotons<3 || npi0s<1)return; //for pi0, eta and omega->pi0+gamma->3gamma reconstruction\r
- \r
- //reconstruction of omega(782)->pi0+gamma->3gamma\r
- //loop for pi0 and photon\r
- if(GetDebug() > 0) printf("omega->pi0+gamma->3gamma invariant mass analysis ! This event have %d photons and %d pi0 \n", nphotons, npi0s);\r
- for(Int_t i=0;i<npi0s;i++){\r
- AliAODPWG4Particle * pi0 = (AliAODPWG4Particle*) (fInputAODPi0->At(i)) ; //pi0\r
- TLorentzVector vpi0(pi0->Px(),pi0->Py(),pi0->Pz(),pi0->E());\r
- Int_t lab1=pi0->GetCaloLabel(0); // photon1 from pi0 decay\r
- Int_t lab2=pi0->GetCaloLabel(1); // photon2 from pi0 decay\r
- //for omega->pi0+gamma, it needs at least three photons per event\r
- //Get the two decay photons from pi0\r
- AliAODPWG4Particle * photon1 =0;\r
- AliAODPWG4Particle * photon2 =0;\r
- for(Int_t d1=0;d1<nphotons;d1++){\r
- for(Int_t d2=0;d2<nphotons;d2++){\r
- AliAODPWG4Particle * dp1 = (AliAODPWG4Particle*) (aodGamma->At(d1));\r
- AliAODPWG4Particle * dp2 = (AliAODPWG4Particle*) (aodGamma->At(d2));\r
- Int_t dlab1=dp1->GetCaloLabel(0);\r
- Int_t dlab2=dp2->GetCaloLabel(0);\r
- if(dlab1==lab1 && dlab2==lab2){\r
- photon1=dp1;\r
- photon2=dp2;\r
- }\r
- else continue;\r
- }\r
- }\r
- //caculate the asy and dist of the two photon from pi0 decay\r
- TLorentzVector dph1(photon1->Px(),photon1->Py(),photon1->Pz(),photon1->E());\r
- TLorentzVector dph2(photon2->Px(),photon2->Py(),photon2->Pz(),photon2->E());\r
- \r
- Double_t pi0asy= TMath::Abs(dph1.E()-dph2.E())/(dph1.E()+dph2.E());\r
- // Double_t phi1=dph1.Phi();\r
- // Double_t phi2=dph2.Phi();\r
- // Double_t eta1=dph1.Eta();\r
- // Double_t eta2=dph2.Eta();\r
- // Double_t pi0dist=TMath::Sqrt((phi1-phi2)*(phi1-phi2)+(eta1-eta2)*(eta1-eta2));\r
- \r
- if(pi0->GetPdg()==111 && nphotons>2 && npi0s\r
- && TMath::Abs(vpi0.M()-fPi0Mass)<fPi0MassWindow) { //pi0 candidates\r
- \r
- //avoid the double counting\r
- Int_t * dc1= new Int_t[nphotons];\r
- Int_t * dc2= new Int_t[nphotons];\r
- Int_t index1=0;\r
- Int_t index2=0;\r
- for(Int_t k=0;k<i;k++){\r
- AliAODPWG4Particle * p3=(AliAODPWG4Particle*)(fInputAODPi0->At(k));\r
- Int_t lab4=p3->GetCaloLabel(0);\r
- Int_t lab5=p3->GetCaloLabel(1);\r
- if(lab1==lab4){ dc1[index1]=lab5; index1++; }\r
- if(lab2==lab5){ dc2[index2]=lab4; index2++; }\r
- }\r
- \r
- \r
- //loop the pi0 with third gamma\r
- for(Int_t j=0;j<nphotons;j++){\r
- AliAODPWG4Particle *photon3 = (AliAODPWG4Particle*) (aodGamma->At(j));\r
- TLorentzVector dph3(photon3->Px(),photon3->Py(),photon3->Pz(),photon3->E());\r
- Int_t lab3=photon3->GetCaloLabel(0);\r
- Double_t pi0gammapt=(vpi0+dph3).Pt();\r
- Double_t pi0gammamass=(vpi0+dph3).M();\r
- Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt; \r
- Double_t gammaOverOmegaPtRatio= dph3.Pt()/pi0gammapt;\r
- \r
- //pi0, gamma pt cut \r
- if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut || \r
- gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;\r
- \r
- for(Int_t l=0;l<index1;l++) if(lab3==dc1[l]) lab3=-1;\r
- for(Int_t l=0;l<index2;l++) if(lab3==dc2[l]) lab3=-1;\r
- \r
- if(lab3>0 && lab3!=lab1 && lab3!=lab2){\r
- for(Int_t ipid=0;ipid<fNpid;ipid++){\r
- for(Int_t idist=0;idist<fNBadChDistBin;idist++){\r
- Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;\r
- if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
- photon2->IsPIDOK(ipid,AliCaloPID::kPhoton) && \r
- photon3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
- photon1->DistToBad()>=idist &&\r
- photon2->DistToBad()>=idist &&\r
- photon3->DistToBad()>=idist ){\r
- //fill the histograms\r
- if(GetDebug() > 2) printf("Real: index %d pt %2.3f mass %2.3f \n", index, pi0gammapt, pi0gammamass);\r
- fRealOmega0[index]->Fill(pi0gammapt,pi0gammamass);\r
- if(pi0asy<0.7) fRealOmega1[index]->Fill(pi0gammapt,pi0gammamass);\r
- if(pi0asy<0.8) fRealOmega2[index]->Fill(pi0gammapt,pi0gammamass);\r
- }\r
- }\r
- }\r
- }\r
- } \r
- delete []dc1;\r
- delete []dc2;\r
- if(GetDebug() > 0) printf("MixA: (r1_event1+r2_event1)+r3_event2 \n");\r
- //-------------------------\r
- //background analysis\r
- //three background\r
- // --A (r1_event1+r2_event1)+r3_event2\r
- Int_t nMixed = fEventsList[curEventBin]->GetSize();\r
- for(Int_t im=0;im<nMixed;im++){\r
- TClonesArray* ev2= (TClonesArray*) (fEventsList[curEventBin]->At(im));\r
- for(Int_t mix1=0;mix1<ev2->GetEntries();mix1++){\r
- AliAODPWG4Particle *mix1ph = (AliAODPWG4Particle*) (ev2->At(mix1)); \r
- TLorentzVector vmixph(mix1ph->Px(),mix1ph->Py(),mix1ph->Pz(),mix1ph->E());\r
- Double_t pi0gammapt=(vpi0+vmixph).Pt();\r
- Double_t pi0gammamass=(vpi0+vmixph).M();\r
- Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt;\r
- Double_t gammaOverOmegaPtRatio= vmixph.Pt()/pi0gammapt;\r
- \r
- //pi0, gamma pt cut \r
- if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut || \r
- gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;\r
- \r
- for(Int_t ipid=0;ipid<fNpid;ipid++){\r
- for(Int_t idist=0;idist<fNBadChDistBin;idist++){\r
- Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;\r
- if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton)&&\r
- photon2->IsPIDOK(ipid,AliCaloPID::kPhoton)&&\r
- mix1ph->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
- photon1->DistToBad()>=idist &&\r
- photon2->DistToBad()>=idist &&\r
- mix1ph->DistToBad()>=idist ){\r
- if(GetDebug() > 2) printf("MixA: index %d pt %2.3f mass %2.3f \n",index, pi0gammapt, pi0gammamass);\r
- //fill the histograms\r
- fMixAOmega0[index]->Fill(pi0gammapt,pi0gammamass);\r
- if(pi0asy<0.7)fMixAOmega1[index]->Fill(pi0gammapt,pi0gammamass);\r
- if(pi0asy<0.8)fMixAOmega2[index]->Fill(pi0gammapt,pi0gammamass);\r
- //printf("mix A %d %2.2f \n", index, pi0gammamass);\r
- \r
- }\r
- }\r
- }\r
- }\r
- }\r
- }\r
- }\r
- \r
- //\r
- // --B (r1_event1+r2_event2)+r3_event2\r
- //\r
- if(GetDebug() >0)printf("MixB: (r1_event1+r2_event2)+r3_event2 \n");\r
- for(Int_t i=0;i<nphotons;i++){\r
- AliAODPWG4Particle *ph1 = (AliAODPWG4Particle*) (aodGamma->At(i)); \r
- TLorentzVector vph1(ph1->Px(),ph1->Py(),ph1->Pz(),ph1->E());\r
- //interrupt here...................\r
- //especially for EMCAL\r
- //we suppose the high pt clusters are overlapped pi0 \r
-\r
- for(Int_t j=i+1;j<nphotons;j++){\r
- AliAODPWG4Particle *ph2 = (AliAODPWG4Particle*) (aodGamma->At(j));\r
- TLorentzVector fakePi0, fakeOmega, vph;\r
-\r
- if(ph1->E() > fEOverlapCluster && ph1->E() > ph2->E()) {\r
- fakePi0.SetPxPyPzE(ph1->Px(),ph1->Py(),ph1->Pz(),TMath::Sqrt(ph1->Px()*ph1->Px()+ph1->Py()*ph1->Py()+ph1->Pz()*ph1->Pz()+0.135*0.135));\r
- vph.SetPxPyPzE(ph2->Px(),ph2->Py(),ph2->Pz(),ph2->E());\r
- }\r
- else if(ph2->E() > fEOverlapCluster && ph2->E() > ph1->E()) {\r
- fakePi0.SetPxPyPzE(ph2->Px(),ph2->Py(),ph2->Pz(),TMath::Sqrt(ph2->Px()*ph2->Px()+ph2->Py()*ph2->Py()+ph2->Pz()*ph2->Pz()+0.135*0.135));\r
- vph.SetPxPyPzE(ph1->Px(), ph1->Py(),ph1->Pz(),ph1->E());\r
- }\r
- else continue;\r
-\r
- fakeOmega=fakePi0+vph;\r
- for(Int_t ii=0;ii<fNCentBin;ii++){ \r
- fhFakeOmega[icentbin]->Fill(fakeOmega.Pt(), fakeOmega.M());\r
- }\r
- }//j\r
-\r
- //continue ................\r
- Int_t nMixed = fEventsList[curEventBin]->GetSize();\r
- for(Int_t ie=0;ie<nMixed;ie++){\r
- TClonesArray* ev2= (TClonesArray*) (fEventsList[curEventBin]->At(ie));\r
- for(Int_t mix1=0;mix1<ev2->GetEntries();mix1++){\r
- AliAODPWG4Particle *ph2 = (AliAODPWG4Particle*) (ev2->At(mix1));\r
- TLorentzVector vph2(ph2->Px(),ph2->Py(),ph2->Pz(),ph2->E());\r
- Double_t pi0asy = TMath::Abs(vph1.E()-vph2.E())/(vph1.E()+vph2.E()); \r
- Double_t pi0mass=(vph1+vph2).M();\r
- \r
- if(TMath::Abs(pi0mass-fPi0Mass)<fPi0MassWindow){//for pi0 selection\r
- for(Int_t mix2=(mix1+1);mix2<ev2->GetEntries();mix2++){\r
- AliAODPWG4Particle *ph3 = (AliAODPWG4Particle*) (ev2->At(mix2));\r
- TLorentzVector vph3(ph3->Px(),ph3->Py(),ph3->Pz(),ph3->E());\r
- \r
- Double_t pi0gammapt=(vph1+vph2+vph3).Pt();\r
- Double_t pi0gammamass=(vph1+vph2+vph3).M(); \r
- Double_t pi0OverOmegaPtRatio =(vph1+vph2).Pt()/pi0gammapt;\r
- Double_t gammaOverOmegaPtRatio= vph3.Pt()/pi0gammapt;\r
- //pi0, gamma pt cut \r
- if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||\r
- gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;\r
- \r
- for(Int_t ipid=0;ipid<fNpid;ipid++){\r
- for(Int_t idist=0;idist<fNBadChDistBin;idist++){\r
- Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;\r
- if(ph1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
- ph2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
- ph3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
- ph1->DistToBad()>=idist &&\r
- ph2->DistToBad()>=idist &&\r
- ph3->DistToBad()>=idist ){\r
- if(GetDebug() > 2) printf("MixB: index %d pt %2.3f mass %2.3f \n", index, pi0gammapt, pi0gammamass);\r
- //fill histograms\r
- fMixBOmega0[index]->Fill(pi0gammapt,pi0gammamass);\r
- if(pi0asy<0.7) fMixBOmega1[index]->Fill(pi0gammapt,pi0gammamass);\r
- if(pi0asy<0.8) fMixBOmega2[index]->Fill(pi0gammapt,pi0gammamass);\r
- //printf("mix B %d %2.2f \n", index, pi0gammamass);\r
- }\r
- } \r
- }\r
- }\r
- \r
- //\r
- // --C (r1_event1+r2_event2)+r3_event3\r
- //\r
- if(GetDebug() >0)printf("MixC: (r1_event1+r2_event2)+r3_event3\n");\r
- for(Int_t je=(ie+1);je<nMixed;je++){\r
- TClonesArray* ev3= (TClonesArray*) (fEventsList[curEventBin]->At(je));\r
- for(Int_t mix3=0;mix3<ev3->GetEntries();mix3++){\r
- AliAODPWG4Particle *ph3 = (AliAODPWG4Particle*) (ev3->At(mix3));\r
- TLorentzVector vph3(ph3->Px(),ph3->Py(),ph3->Pz(),ph3->E());\r
- \r
- Double_t pi0gammapt=(vph1+vph2+vph3).Pt();\r
- Double_t pi0gammamass=(vph1+vph2+vph3).M();\r
- Double_t pi0OverOmegaPtRatio =(vph1+vph2).Pt()/pi0gammapt;\r
- Double_t gammaOverOmegaPtRatio= vph3.Pt()/pi0gammapt;\r
- //pi0, gamma pt cut \r
- if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||\r
- gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;\r
- \r
- for(Int_t ipid=0;ipid<fNpid;ipid++){\r
- for(Int_t idist=0;idist<fNBadChDistBin;idist++){\r
- Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;\r
- if(ph1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
- ph2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
- ph3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&\r
- ph1->DistToBad()>=idist &&\r
- ph2->DistToBad()>=idist &&\r
- ph3->DistToBad()>=idist ){\r
- if(GetDebug() > 2) printf("MixC: index %d pt %2.3f mass %2.3f \n", index, pi0gammapt, pi0gammamass);\r
- //fill histograms\r
- fMixCOmega0[index]->Fill(pi0gammapt,pi0gammamass);\r
- if(pi0asy<0.7) fMixCOmega1[index]->Fill(pi0gammapt,pi0gammamass);\r
- if(pi0asy<0.8) fMixCOmega2[index]->Fill(pi0gammapt,pi0gammamass);\r
- //printf("mix C %d %2.2f \n", index, pi0gammamass);\r
- }\r
- }\r
- }\r
- }\r
- }\r
- } //for pi0 selecton \r
- }\r
- }\r
- }\r
- \r
- \r
- //event buffer \r
- TClonesArray *currentEvent = new TClonesArray(*aodGamma);\r
- if(currentEvent->GetEntriesFast()>0){\r
- fEventsList[curEventBin]->AddFirst(currentEvent) ;\r
- currentEvent=0 ; \r
- if(fEventsList[curEventBin]->GetSize()>=fNmaxMixEv) {\r
- TClonesArray * tmp = (TClonesArray*) (fEventsList[curEventBin]->Last()) ;\r
- fEventsList[curEventBin]->RemoveLast() ;\r
- delete tmp ;\r
- }\r
- }\r
- else{ \r
- delete currentEvent ;\r
- currentEvent=0 ;\r
- }\r
- \r
-}\r
-\r
-//______________________________________________________________________________\r
-void AliAnaOmegaToPi0Gamma::ReadHistograms(TList * outputList)\r
-{\r
- //read the histograms \r
- //for the finalization of the terminate analysis\r
-\r
- Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"RealToPi0Gamma_Vz0C0Rp0Pid0Dist0"));\r
-\r
- Int_t ndim=fNVtxZBin*fNCentBin*fNRpBin*fNBadChDistBin*fNpid;\r
-\r
- if(!fRealOmega0) fRealOmega0 =new TH2F*[ndim];\r
- if(!fMixAOmega0) fMixAOmega0 =new TH2F*[ndim];\r
- if(!fMixBOmega0) fMixBOmega0 =new TH2F*[ndim];\r
- if(!fMixCOmega0) fMixCOmega0 =new TH2F*[ndim];\r
-\r
- if(!fRealOmega1) fRealOmega1 =new TH2F*[ndim];\r
- if(!fMixAOmega1) fMixAOmega1 =new TH2F*[ndim];\r
- if(!fMixBOmega1) fMixBOmega1 =new TH2F*[ndim];\r
- if(!fMixCOmega1) fMixCOmega1 =new TH2F*[ndim];\r
-\r
- if(!fRealOmega2) fRealOmega2 =new TH2F*[ndim];\r
- if(!fMixAOmega2) fMixAOmega2 =new TH2F*[ndim];\r
- if(!fMixBOmega2) fMixBOmega2 =new TH2F*[ndim];\r
- if(!fMixCOmega2) fMixCOmega2 =new TH2F*[ndim];\r
-\r
- for(Int_t i=0;i<fNVtxZBin;i++){\r
- for(Int_t j=0;j<fNCentBin;j++){\r
- for(Int_t k=0;k<fNRpBin;k++){ //at event level\r
- Int_t idim=i*fNCentBin*fNRpBin+j*fNRpBin+k;\r
- for(Int_t ipid=0;ipid<fNpid;ipid++){ \r
- for(Int_t idist=0;idist<fNBadChDistBin;idist++){ //at particle\r
- Int_t ind=idim*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;\r
- fRealOmega0[ind]= (TH2F*) outputList->At(index++);\r
- fMixAOmega0[ind]= (TH2F*) outputList->At(index++);\r
- fMixBOmega0[ind]= (TH2F*) outputList->At(index++);\r
- fMixCOmega0[ind]= (TH2F*) outputList->At(index++);\r
-\r
- fRealOmega1[ind]= (TH2F*) outputList->At(index++);\r
- fMixAOmega1[ind]= (TH2F*) outputList->At(index++);\r
- fMixBOmega1[ind]= (TH2F*) outputList->At(index++);\r
- fMixCOmega1[ind]= (TH2F*) outputList->At(index++);\r
-\r
- fRealOmega2[ind]= (TH2F*) outputList->At(index++);\r
- fMixAOmega2[ind]= (TH2F*) outputList->At(index++);\r
- fMixBOmega2[ind]= (TH2F*) outputList->At(index++);\r
- fMixCOmega2[ind]= (TH2F*) outputList->At(index++);\r
- \r
- \r
- }\r
- }\r
- }\r
- }\r
- }\r
- \r
- if(IsDataMC()){\r
- fhOmegaPriPt = (TH1F*) outputList->At(index++);\r
- }\r
-\r
-}\r
-\r
-//______________________________________________________________________________\r
-void AliAnaOmegaToPi0Gamma::Terminate(TList * outputList) \r
-{\r
-// //Do some calculations and plots from the final histograms.\r
- if(GetDebug() >= 0) printf("AliAnaOmegaToPi0Gamma::Terminate() \n");\r
- ReadHistograms(outputList);\r
- const Int_t buffersize = 255;\r
- char cvs1[buffersize]; \r
- snprintf(cvs1, buffersize, "Neutral_%s_IVM",fInputAODGammaName.Data());\r
-\r
- TCanvas * cvsIVM = new TCanvas(cvs1, cvs1, 400, 10, 600, 700) ;\r
- cvsIVM->Divide(2, 2);\r
-\r
- cvsIVM->cd(1);\r
- char dec[buffersize];\r
- snprintf(dec,buffersize,"h2Real_%s",fInputAODGammaName.Data());\r
- TH2F * h2Real= (TH2F*)fRealOmega0[0]->Clone(dec);\r
- h2Real->GetXaxis()->SetRangeUser(4,6);\r
- TH1F * hRealOmega = (TH1F*) h2Real->ProjectionY();\r
- hRealOmega->SetTitle("RealPi0Gamma 4<pt<6");\r
- hRealOmega->SetLineColor(2);\r
- hRealOmega->Draw();\r
-\r
- cvsIVM->cd(2);\r
- snprintf(dec,buffersize,"hMixA_%s",fInputAODGammaName.Data());\r
- TH2F *h2MixA= (TH2F*)fMixAOmega0[0]->Clone(dec);\r
- h2MixA->GetXaxis()->SetRangeUser(4,6);\r
- TH1F * hMixAOmega = (TH1F*) h2MixA->ProjectionY();\r
- hMixAOmega->SetTitle("MixA 4<pt<6");\r
- hMixAOmega->SetLineColor(2);\r
- hMixAOmega->Draw();\r
-\r
- cvsIVM->cd(3);\r
- snprintf(dec,buffersize,"hMixB_%s",fInputAODGammaName.Data());\r
- TH2F * h2MixB= (TH2F*)fMixBOmega0[0]->Clone(dec);\r
- h2MixB->GetXaxis()->SetRangeUser(4,6);\r
- TH1F * hMixBOmega = (TH1F*) h2MixB->ProjectionY();\r
- hMixBOmega->SetTitle("MixB 4<pt<6");\r
- hMixBOmega->SetLineColor(2);\r
- hMixBOmega->Draw();\r
-\r
- cvsIVM->cd(4);\r
- snprintf(dec,buffersize,"hMixC_%s",fInputAODGammaName.Data());\r
- TH2F *h2MixC= (TH2F*)fMixCOmega0[0]->Clone(dec);\r
- h2MixC->GetXaxis()->SetRangeUser(4,6);\r
- TH1F * hMixCOmega = (TH1F*) h2MixC->ProjectionY();\r
- hMixCOmega->SetTitle("MixC 4<pt<6");\r
- hMixCOmega->SetLineColor(2);\r
- hMixCOmega->Draw();\r
-\r
- char eps[buffersize];\r
- snprintf(eps,buffersize,"CVS_%s_IVM.eps",fInputAODGammaName.Data());\r
- cvsIVM->Print(eps);\r
- cvsIVM->Modified();\r
- \r
-}\r
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+/* $Id: $ */
+//_________________________________________________________________________
+// class to extract omega(782)->pi0+gamma->3gamma
+// Mar. 22, 2011: Additional method, espeically for EMCAL. A high E cluster is assumpted as pi0 (two photons are overlapped) without unfolding
+//
+//-- Author: Renzhuo Wan (IOPP-Wuhan, China)
+//_________________________________________________________________________
+
+// --- ROOT system
+class TROOT;
+
+// --- AliRoot system
+//class AliVEvent;
+// --- ROOT system ---
+#include "TH2F.h"
+#include "TLorentzVector.h"
+#include "TParticle.h"
+#include "TCanvas.h"
+#include "TFile.h"
+//---- AliRoot system ----
+#include "AliAnaOmegaToPi0Gamma.h"
+#include "AliCaloTrackReader.h"
+#include "AliCaloPID.h"
+#include "AliStack.h"
+#include "AliVEvent.h"
+#include "AliAODEvent.h"
+#include "AliAODMCParticle.h"
+ClassImp(AliAnaOmegaToPi0Gamma)
+
+//______________________________________________________________________________
+AliAnaOmegaToPi0Gamma::AliAnaOmegaToPi0Gamma() : AliAnaPartCorrBaseClass(),
+fInputAODPi0(0), fInputAODGammaName(""),
+fEventsList(0x0),fNVtxZBin(0), fNCentBin(0), fNRpBin(0), fNBadChDistBin(0), fNpid(0),
+fNmaxMixEv(0), fVtxZCut(0), fCent(0), fRp(0),
+fPi0Mass(0),fPi0MassWindow(0),fPi0OverOmegaPtCut(0),
+fGammaOverOmegaPtCut(0), fEOverlapCluster(0),
+fhEtalon(0),
+fRealOmega0(0), fMixAOmega0(0),
+fMixBOmega0(0), fMixCOmega0(0),
+fRealOmega1(0), fMixAOmega1(0),
+fMixBOmega1(0), fMixCOmega1(0),
+fRealOmega2(0), fMixAOmega2(0),
+fMixBOmega2(0), fMixCOmega2(0),
+fhFakeOmega(0),
+fhOmegaPriPt(0)
+{
+ //Default Ctor
+ InitParameters();
+}
+/*
+//______________________________________________________________________________
+AliAnaOmegaToPi0Gamma::AliAnaOmegaToPi0Gamma(const AliAnaOmegaToPi0Gamma & ex) : AliAnaPartCorrBaseClass(ex),
+fInputAODPi0(new TClonesArray (*ex.fInputAODPi0)),
+fInputAODGammaName(ex.fInputAODGammaName),
+fEventsList(0x0),
+fNVtxZBin(ex.fNVtxZBin), fNCentBin(ex.fNCentBin), fNRpBin(ex.fNRpBin),
+fNBadChDistBin(ex.fNBadChDistBin),fNpid(ex.fNpid),
+fNmaxMixEv(ex.fNmaxMixEv),
+fVtxZCut(ex.fVtxZCut), fCent(ex.fCent), fRp(ex.fRp),
+fPi0Mass(ex.fPi0Mass),
+fPi0MassWindow(ex.fPi0MassWindow),
+fPi0OverOmegaPtCut(ex.fPi0OverOmegaPtCut),
+fGammaOverOmegaPtCut(ex.fGammaOverOmegaPtCut),
+fhEtalon(ex.fhEtalon),
+fRealOmega0(ex.fRealOmega0), fMixAOmega0(ex.fMixAOmega0),
+fMixBOmega0(ex.fMixBOmega0), fMixCOmega0(ex.fMixCOmega0),
+fRealOmega1(ex.fRealOmega1), fMixAOmega1(ex.fMixAOmega1),
+fMixBOmega1(ex.fMixBOmega1), fMixCOmega1(ex.fMixCOmega1),
+fRealOmega2(ex.fRealOmega2), fMixAOmega2(ex.fMixAOmega2),
+fMixBOmega2(ex.fMixBOmega2), fMixCOmega2(ex.fMixCOmega2),
+fhOmegaPriPt(ex.fhOmegaPriPt)
+{
+ // cpy ctor
+ //Do not need it
+}
+*/
+/*
+//______________________________________________________________________________
+AliAnaOmegaToPi0Gamma & AliAnaOmegaToPi0Gamma::operator = (const AliAnaOmegaToPi0Gamma & ex)
+{
+ // assignment operator
+
+ if(this == &ex)return *this;
+ ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
+ fInputAODPi0 = new TClonesArray(*ex.fInputAODPi0);
+ fInputAODGammaName = ex.fInputAODGammaName;
+ fEventsList = ex.fEventsList;
+
+ fNVtxZBin=ex.fNVtxZBin;
+ fNCentBin=ex.fNCentBin;
+ fNRpBin=ex.fNRpBin;
+ fNBadChDistBin=ex.fNBadChDistBin;
+ fNpid=ex.fNpid;
+ fNmaxMixEv =ex.fNmaxMixEv;
+
+ fVtxZCut=ex.fVtxZCut;
+ fCent=ex.fCent;
+ fRp=ex.fRp;
+
+ fPi0Mass=ex.fPi0Mass;
+ fPi0MassWindow=ex.fPi0MassWindow;
+ fPi0OverOmegaPtCut=ex.fPi0OverOmegaPtCut;
+ fGammaOverOmegaPtCut=ex.fGammaOverOmegaPtCut;
+
+ fhEtalon=ex.fhEtalon;
+ fRealOmega0=ex.fRealOmega0;
+ fMixAOmega0=ex.fMixAOmega0;
+ fMixBOmega0=ex.fMixBOmega0;
+ fMixCOmega0=ex.fMixCOmega0;
+ fRealOmega1=ex.fRealOmega1;
+ fMixAOmega1=ex.fMixAOmega1;
+ fMixBOmega1=ex.fMixBOmega1;
+ fMixCOmega1=ex.fMixCOmega1;
+ fRealOmega2=ex.fRealOmega2;
+ fMixAOmega2=ex.fMixAOmega2;
+ fMixBOmega2=ex.fMixBOmega2;
+ fMixCOmega2=ex.fMixCOmega2;
+ fhOmegaPriPt=ex.fhOmegaPriPt;
+ return *this;
+
+}
+*/
+
+//______________________________________________________________________________
+AliAnaOmegaToPi0Gamma::~AliAnaOmegaToPi0Gamma() {
+
+ //dtor
+// Done by the maker
+// if(fInputAODPi0){
+// fInputAODPi0->Clear();
+// delete fInputAODPi0;
+// }
+
+ if(fEventsList){
+ for(Int_t i=0;i<fNVtxZBin;i++){
+ for(Int_t j=0;j<fNCentBin;j++){
+ for(Int_t k=0;k<fNRpBin;k++){
+ fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k]->Clear();
+ delete fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k];
+ }
+ }
+ }
+ }
+ delete [] fEventsList;
+ fEventsList=0;
+
+ delete [] fVtxZCut;
+ delete [] fCent;
+ delete [] fRp;
+
+}
+
+//______________________________________________________________________________
+void AliAnaOmegaToPi0Gamma::InitParameters()
+{
+//Init parameters when first called the analysis
+//Set default parameters
+ fInputAODGammaName="PhotonsDetector";
+ fNVtxZBin=1;
+ fNCentBin=1;
+ fNRpBin=1;
+ fNBadChDistBin=3;
+ fNpid=1;
+ fNmaxMixEv=8;
+
+ fPi0Mass=0.1348;
+ fPi0MassWindow=0.015;
+ fPi0OverOmegaPtCut=0.8;
+ fGammaOverOmegaPtCut=0.2;
+
+ fEOverlapCluster=6;
+}
+
+
+//______________________________________________________________________________
+TList * AliAnaOmegaToPi0Gamma::GetCreateOutputObjects()
+{
+ //
+ fVtxZCut = new Double_t [fNVtxZBin];
+ for(Int_t i=0;i<fNVtxZBin;i++) fVtxZCut[i]=10*(i+1);
+
+ fCent=new Double_t[fNCentBin];
+ for(int i = 0;i<fNCentBin;i++)fCent[i]=0;
+
+ fRp=new Double_t[fNRpBin];
+ for(int i = 0;i<fNRpBin;i++)fRp[i]=0;
+ //
+ Int_t nptbins = GetHistoPtBins();
+ Float_t ptmax = GetHistoPtMax();
+ Float_t ptmin = GetHistoPtMin();
+
+ Int_t nmassbins = GetHistoMassBins();
+ Float_t massmin = GetHistoMassMin();
+ Float_t massmax = GetHistoMassMax();
+
+ fhEtalon = new TH2F("hEtalon","Histo with binning parameters", nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
+ fhEtalon->SetXTitle("P_{T} (GeV)") ;
+ fhEtalon->SetYTitle("m_{inv} (GeV)") ;
+
+ // store them in fOutputContainer
+ fEventsList = new TList*[fNVtxZBin*fNCentBin*fNRpBin];
+ for(Int_t i=0;i<fNVtxZBin;i++){
+ for(Int_t j=0;j<fNCentBin;j++){
+ for(Int_t k=0;k<fNRpBin;k++){
+ fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k] = new TList();
+ fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k]->SetOwner(kFALSE);
+ }
+ }
+ }
+
+ TList * outputContainer = new TList() ;
+ outputContainer->SetName(GetName());
+ const Int_t buffersize = 255;
+ char key[buffersize] ;
+ char title[buffersize] ;
+ const char * detector= fInputAODGammaName.Data();
+ Int_t ndim=fNVtxZBin*fNCentBin*fNRpBin*fNBadChDistBin*fNpid;
+
+ fRealOmega0 =new TH2F*[ndim];
+ fMixAOmega0 =new TH2F*[ndim];
+ fMixBOmega0 =new TH2F*[ndim];
+ fMixCOmega0 =new TH2F*[ndim];
+
+ fRealOmega1 =new TH2F*[ndim];
+ fMixAOmega1 =new TH2F*[ndim];
+ fMixBOmega1 =new TH2F*[ndim];
+ fMixCOmega1 =new TH2F*[ndim];
+
+ fRealOmega2 =new TH2F*[ndim];
+ fMixAOmega2 =new TH2F*[ndim];
+ fMixBOmega2 =new TH2F*[ndim];
+ fMixCOmega2 =new TH2F*[ndim];
+
+ fhFakeOmega = new TH2F*[fNCentBin];
+
+ for(Int_t i=0;i<fNVtxZBin;i++){
+ for(Int_t j=0;j<fNCentBin;j++){
+ for(Int_t k=0;k<fNRpBin;k++){ //at event level
+ Int_t idim=i*fNCentBin*fNRpBin+j*fNRpBin+k;
+ for(Int_t ipid=0;ipid<fNpid;ipid++){
+ for(Int_t idist=0;idist<fNBadChDistBin;idist++){ //at particle level
+
+ Int_t index=idim*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
+
+ snprintf(key,buffersize,"RealToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s Real Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fRealOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fRealOmega0[index]->SetName(key) ;
+ fRealOmega0[index]->SetTitle(title);
+ outputContainer->Add(fRealOmega0[index]);
+
+ snprintf(key,buffersize,"MixAToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s MixA Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fMixAOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fMixAOmega0[index]->SetName(key) ;
+ fMixAOmega0[index]->SetTitle(title);
+ outputContainer->Add(fMixAOmega0[index]);
+
+ snprintf(key,buffersize,"MixBToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s MixB Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fMixBOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fMixBOmega0[index]->SetName(key) ;
+ fMixBOmega0[index]->SetTitle(title);
+ outputContainer->Add(fMixBOmega0[index]);
+
+ snprintf(key,buffersize,"MixCToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s MixC Pi0GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fMixCOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fMixCOmega0[index]->SetName(key) ;
+ fMixCOmega0[index]->SetTitle(title);
+ outputContainer->Add(fMixCOmega0[index]);
+
+ snprintf(key,buffersize,"RealToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s Real Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fRealOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fRealOmega1[index]->SetName(key) ;
+ fRealOmega1[index]->SetTitle(title);
+ outputContainer->Add(fRealOmega1[index]);
+
+ snprintf(key,buffersize,"MixAToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s MixA Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fMixAOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fMixAOmega1[index]->SetName(key) ;
+ fMixAOmega1[index]->SetTitle(title);
+ outputContainer->Add(fMixAOmega1[index]);
+
+ snprintf(key,buffersize,"MixBToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s MixB Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fMixBOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fMixBOmega1[index]->SetName(key) ;
+ fMixBOmega1[index]->SetTitle(title);
+ outputContainer->Add(fMixBOmega1[index]);
+
+ snprintf(key,buffersize,"MixCToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s MixC Pi0(A<0.7)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fMixCOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fMixCOmega1[index]->SetName(key) ;
+ fMixCOmega1[index]->SetTitle(title);
+ outputContainer->Add(fMixCOmega1[index]);
+
+ snprintf(key,buffersize,"RealToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s Real Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fRealOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fRealOmega2[index]->SetName(key) ;
+ fRealOmega2[index]->SetTitle(title);
+ outputContainer->Add(fRealOmega2[index]);
+
+ snprintf(key,buffersize,"MixAToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s MixA Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fMixAOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fMixAOmega2[index]->SetName(key) ;
+ fMixAOmega2[index]->SetTitle(title);
+ outputContainer->Add(fMixAOmega2[index]);
+
+ snprintf(key,buffersize,"MixBToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s MixB Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fMixBOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fMixBOmega2[index]->SetName(key) ;
+ fMixBOmega2[index]->SetTitle(title);
+ outputContainer->Add(fMixBOmega2[index]);
+
+ snprintf(key,buffersize,"MixCToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
+ snprintf(title,buffersize, "%s MixC Pi0(A<0.8)GammaIVM vz_%2.1f_ct_%2.1f_Rp_%2.1f_pid_%d_dist_%d",detector,fVtxZCut[i],fCent[j],fRp[k],ipid,idist);
+ fMixCOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
+ fMixCOmega2[index]->SetName(key) ;
+ fMixCOmega2[index]->SetTitle(title);
+ outputContainer->Add(fMixCOmega2[index]);
+ }
+ }
+ }
+ }
+ }
+
+ for(Int_t i=0;i<fNCentBin;i++){
+ snprintf(key,buffersize, "fhFakeOmega%d",i);
+ snprintf(title,buffersize,"FakePi0(high pt cluster)+Gamma with Centrality%d",i);
+ fhFakeOmega[i] = (TH2F*)fhEtalon->Clone(key);
+ fhFakeOmega[i]->SetTitle(title);
+ outputContainer->Add(fhFakeOmega[i]);
+ }
+ if(IsDataMC()){
+ snprintf(key,buffersize, "%sOmegaPri",detector);
+ snprintf(title,buffersize,"primary #omega in %s",detector);
+ fhOmegaPriPt=new TH1F(key, title,nptbins,ptmin,ptmax);
+ fhOmegaPriPt->GetXaxis()->SetTitle("P_{T}");
+ fhOmegaPriPt->GetYaxis()->SetTitle("dN/P_{T}");
+ outputContainer->Add(fhOmegaPriPt);
+ }
+
+ delete fhEtalon;
+ return outputContainer;
+}
+
+//______________________________________________________________________________
+void AliAnaOmegaToPi0Gamma::Print(const Option_t * /*opt*/) const
+{
+ //Print some relevant parameters set in the analysis
+ printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
+ AliAnaPartCorrBaseClass::Print(" ");
+ printf("Omega->pi0+gamma->3gamma\n");
+ printf("Cuts at event level: \n");
+ printf("Bins of vertex Z: %d \n", fNVtxZBin);
+ printf("Bins of centrality: %d \n",fNCentBin);
+ printf("Bins of Reaction plane: %d\n",fNRpBin);
+ printf("Cuts at AOD particle level:\n");
+ printf("Number of PID: %d \n", fNpid);
+ printf("Number of DistToBadChannel cuts: %d\n", fNBadChDistBin);
+ printf("number of events buffer to be mixed: %d\n",fNmaxMixEv);
+}
+
+//______________________________________________________________________________
+void AliAnaOmegaToPi0Gamma::MakeAnalysisFillHistograms()
+{
+ //fill the MC AOD if needed first
+ //-----------
+ //need to be further implemented
+ AliStack * stack = 0x0;
+ // TParticle * primary = 0x0;
+ TClonesArray * mcparticles0 = 0x0;
+ //TClonesArray * mcparticles1 = 0x0;
+ AliAODMCParticle * aodprimary = 0x0;
+ Int_t pdg=0;
+ Double_t pt=0;
+ Double_t eta=0;
+
+ if(IsDataMC()){
+ if(GetReader()->ReadStack()){
+ stack = GetMCStack() ;
+ if(!stack){
+ printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - There is no stack!\n");
+ }
+ else{
+ for(Int_t i=0 ; i<stack->GetNtrack(); i++){
+ TParticle * prim = stack->Particle(i) ;
+ pdg = prim->GetPdgCode() ;
+ eta=prim->Eta();
+ pt=prim->Pt();
+ if(TMath::Abs(eta)<0.5) {
+ if(pdg==223) fhOmegaPriPt->Fill(pt);
+ }
+ }
+ }
+ }
+ else if(GetReader()->ReadAODMCParticles()){
+ //Get the list of MC particles
+ mcparticles0 = GetReader()->GetAODMCParticles(0);
+ if(!mcparticles0 ) {
+ if(GetDebug() > 0) printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - Standard MCParticles not available!\n");
+ }
+ // if(GetReader()->GetSecondInputAODTree()){
+ // mcparticles1 = GetReader()->GetAODMCParticles(1);
+ // if(!mcparticles1 && GetDebug() > 0) {
+ // printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - Second input MCParticles not available!\n");
+ // }
+ // }
+ else{
+ for(Int_t i=0;i<mcparticles0->GetEntries();i++){
+ aodprimary =(AliAODMCParticle*)mcparticles0->At(i);
+ pdg = aodprimary->GetPdgCode() ;
+ eta=aodprimary->Eta();
+ pt=aodprimary->Pt();
+ if(TMath::Abs(eta)<0.5) {
+ if(pdg==223) fhOmegaPriPt->Fill(pt);
+ }
+
+ }
+ }//mcparticles0 exists
+ }//AOD MC Particles
+ }// is data and MC
+
+
+ //process event from AOD brach
+ //extract pi0, eta and omega analysis
+ Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
+ if(IsBadRun(iRun)) return ;
+
+ //vertex z
+ Double_t vert[]={0,0,0} ;
+ GetVertex(vert);
+ Int_t curEventBin =0;
+
+ Int_t ivtxzbin=(Int_t)TMath::Abs(vert[2])/10;
+ if(ivtxzbin>=fNVtxZBin)return;
+
+ //centrality
+ Int_t currentCentrality = GetEventCentrality();
+ if(currentCentrality == -1) return;
+ Int_t optCent = GetReader()->GetCentralityOpt();
+ Int_t icentbin=currentCentrality/(optCent/fNCentBin) ; //GetEventCentrality();
+
+ printf("-------------- %d %d %d ",currentCentrality, optCent, icentbin);
+ //reaction plane
+ Int_t irpbin=0;
+
+ if(ivtxzbin==-1) return;
+ curEventBin = ivtxzbin*fNCentBin*fNRpBin + icentbin*fNRpBin + irpbin;
+ TClonesArray *aodGamma = (TClonesArray*) GetAODBranch(fInputAODGammaName); //photon array
+ // TClonesArray *aodGamma = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaName); //photon array
+ Int_t nphotons =0;
+ if(aodGamma) nphotons= aodGamma->GetEntries();
+ else return;
+
+ fInputAODPi0 = (TClonesArray*)GetInputAODBranch(); //pi0 array
+ Int_t npi0s = 0;
+ if(fInputAODPi0) npi0s= fInputAODPi0 ->GetEntries();
+ else return;
+
+ if(nphotons<3 || npi0s<1)return; //for pi0, eta and omega->pi0+gamma->3gamma reconstruction
+
+ //reconstruction of omega(782)->pi0+gamma->3gamma
+ //loop for pi0 and photon
+ if(GetDebug() > 0) printf("omega->pi0+gamma->3gamma invariant mass analysis ! This event have %d photons and %d pi0 \n", nphotons, npi0s);
+ for(Int_t i=0;i<npi0s;i++){
+ AliAODPWG4Particle * pi0 = (AliAODPWG4Particle*) (fInputAODPi0->At(i)) ; //pi0
+ TLorentzVector vpi0(pi0->Px(),pi0->Py(),pi0->Pz(),pi0->E());
+ Int_t lab1=pi0->GetCaloLabel(0); // photon1 from pi0 decay
+ Int_t lab2=pi0->GetCaloLabel(1); // photon2 from pi0 decay
+ //for omega->pi0+gamma, it needs at least three photons per event
+ //Get the two decay photons from pi0
+ AliAODPWG4Particle * photon1 =0;
+ AliAODPWG4Particle * photon2 =0;
+ for(Int_t d1=0;d1<nphotons;d1++){
+ for(Int_t d2=0;d2<nphotons;d2++){
+ AliAODPWG4Particle * dp1 = (AliAODPWG4Particle*) (aodGamma->At(d1));
+ AliAODPWG4Particle * dp2 = (AliAODPWG4Particle*) (aodGamma->At(d2));
+ Int_t dlab1=dp1->GetCaloLabel(0);
+ Int_t dlab2=dp2->GetCaloLabel(0);
+ if(dlab1==lab1 && dlab2==lab2){
+ photon1=dp1;
+ photon2=dp2;
+ }
+ else continue;
+ }
+ }
+ //caculate the asy and dist of the two photon from pi0 decay
+ TLorentzVector dph1(photon1->Px(),photon1->Py(),photon1->Pz(),photon1->E());
+ TLorentzVector dph2(photon2->Px(),photon2->Py(),photon2->Pz(),photon2->E());
+
+ Double_t pi0asy= TMath::Abs(dph1.E()-dph2.E())/(dph1.E()+dph2.E());
+ // Double_t phi1=dph1.Phi();
+ // Double_t phi2=dph2.Phi();
+ // Double_t eta1=dph1.Eta();
+ // Double_t eta2=dph2.Eta();
+ // Double_t pi0dist=TMath::Sqrt((phi1-phi2)*(phi1-phi2)+(eta1-eta2)*(eta1-eta2));
+
+ if(pi0->GetIdentifiedParticleType()==111 && nphotons>2 && npi0s
+ && TMath::Abs(vpi0.M()-fPi0Mass)<fPi0MassWindow) { //pi0 candidates
+
+ //avoid the double counting
+ Int_t * dc1= new Int_t[nphotons];
+ Int_t * dc2= new Int_t[nphotons];
+ Int_t index1=0;
+ Int_t index2=0;
+ for(Int_t k=0;k<i;k++){
+ AliAODPWG4Particle * p3=(AliAODPWG4Particle*)(fInputAODPi0->At(k));
+ Int_t lab4=p3->GetCaloLabel(0);
+ Int_t lab5=p3->GetCaloLabel(1);
+ if(lab1==lab4){ dc1[index1]=lab5; index1++; }
+ if(lab2==lab5){ dc2[index2]=lab4; index2++; }
+ }
+
+
+ //loop the pi0 with third gamma
+ for(Int_t j=0;j<nphotons;j++){
+ AliAODPWG4Particle *photon3 = (AliAODPWG4Particle*) (aodGamma->At(j));
+ TLorentzVector dph3(photon3->Px(),photon3->Py(),photon3->Pz(),photon3->E());
+ Int_t lab3=photon3->GetCaloLabel(0);
+ Double_t pi0gammapt=(vpi0+dph3).Pt();
+ Double_t pi0gammamass=(vpi0+dph3).M();
+ Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt;
+ Double_t gammaOverOmegaPtRatio= dph3.Pt()/pi0gammapt;
+
+ //pi0, gamma pt cut
+ if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
+ gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
+
+ for(Int_t l=0;l<index1;l++) if(lab3==dc1[l]) lab3=-1;
+ for(Int_t l=0;l<index2;l++) if(lab3==dc2[l]) lab3=-1;
+
+ if(lab3>0 && lab3!=lab1 && lab3!=lab2){
+ for(Int_t ipid=0;ipid<fNpid;ipid++){
+ for(Int_t idist=0;idist<fNBadChDistBin;idist++){
+ Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
+ if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
+ photon2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
+ photon3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
+ photon1->DistToBad()>=idist &&
+ photon2->DistToBad()>=idist &&
+ photon3->DistToBad()>=idist ){
+ //fill the histograms
+ if(GetDebug() > 2) printf("Real: index %d pt %2.3f mass %2.3f \n", index, pi0gammapt, pi0gammamass);
+ fRealOmega0[index]->Fill(pi0gammapt,pi0gammamass);
+ if(pi0asy<0.7) fRealOmega1[index]->Fill(pi0gammapt,pi0gammamass);
+ if(pi0asy<0.8) fRealOmega2[index]->Fill(pi0gammapt,pi0gammamass);
+ }
+ }
+ }
+ }
+ }
+ delete []dc1;
+ delete []dc2;
+ if(GetDebug() > 0) printf("MixA: (r1_event1+r2_event1)+r3_event2 \n");
+ //-------------------------
+ //background analysis
+ //three background
+ // --A (r1_event1+r2_event1)+r3_event2
+ Int_t nMixed = fEventsList[curEventBin]->GetSize();
+ for(Int_t im=0;im<nMixed;im++){
+ TClonesArray* ev2= (TClonesArray*) (fEventsList[curEventBin]->At(im));
+ for(Int_t mix1=0;mix1<ev2->GetEntries();mix1++){
+ AliAODPWG4Particle *mix1ph = (AliAODPWG4Particle*) (ev2->At(mix1));
+ TLorentzVector vmixph(mix1ph->Px(),mix1ph->Py(),mix1ph->Pz(),mix1ph->E());
+ Double_t pi0gammapt=(vpi0+vmixph).Pt();
+ Double_t pi0gammamass=(vpi0+vmixph).M();
+ Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt;
+ Double_t gammaOverOmegaPtRatio= vmixph.Pt()/pi0gammapt;
+
+ //pi0, gamma pt cut
+ if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
+ gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
+
+ for(Int_t ipid=0;ipid<fNpid;ipid++){
+ for(Int_t idist=0;idist<fNBadChDistBin;idist++){
+ Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
+ if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton)&&
+ photon2->IsPIDOK(ipid,AliCaloPID::kPhoton)&&
+ mix1ph->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
+ photon1->DistToBad()>=idist &&
+ photon2->DistToBad()>=idist &&
+ mix1ph->DistToBad()>=idist ){
+ if(GetDebug() > 2) printf("MixA: index %d pt %2.3f mass %2.3f \n",index, pi0gammapt, pi0gammamass);
+ //fill the histograms
+ fMixAOmega0[index]->Fill(pi0gammapt,pi0gammamass);
+ if(pi0asy<0.7)fMixAOmega1[index]->Fill(pi0gammapt,pi0gammamass);
+ if(pi0asy<0.8)fMixAOmega2[index]->Fill(pi0gammapt,pi0gammamass);
+ //printf("mix A %d %2.2f \n", index, pi0gammamass);
+
+ }
+ }
+ }
+ }
+ }
+ }
+ }
+
+ //
+ // --B (r1_event1+r2_event2)+r3_event2
+ //
+ if(GetDebug() >0)printf("MixB: (r1_event1+r2_event2)+r3_event2 \n");
+ for(Int_t i=0;i<nphotons;i++){
+ AliAODPWG4Particle *ph1 = (AliAODPWG4Particle*) (aodGamma->At(i));
+ TLorentzVector vph1(ph1->Px(),ph1->Py(),ph1->Pz(),ph1->E());
+ //interrupt here...................
+ //especially for EMCAL
+ //we suppose the high pt clusters are overlapped pi0
+
+ for(Int_t j=i+1;j<nphotons;j++){
+ AliAODPWG4Particle *ph2 = (AliAODPWG4Particle*) (aodGamma->At(j));
+ TLorentzVector fakePi0, fakeOmega, vph;
+
+ if(ph1->E() > fEOverlapCluster && ph1->E() > ph2->E()) {
+ fakePi0.SetPxPyPzE(ph1->Px(),ph1->Py(),ph1->Pz(),TMath::Sqrt(ph1->Px()*ph1->Px()+ph1->Py()*ph1->Py()+ph1->Pz()*ph1->Pz()+0.135*0.135));
+ vph.SetPxPyPzE(ph2->Px(),ph2->Py(),ph2->Pz(),ph2->E());
+ }
+ else if(ph2->E() > fEOverlapCluster && ph2->E() > ph1->E()) {
+ fakePi0.SetPxPyPzE(ph2->Px(),ph2->Py(),ph2->Pz(),TMath::Sqrt(ph2->Px()*ph2->Px()+ph2->Py()*ph2->Py()+ph2->Pz()*ph2->Pz()+0.135*0.135));
+ vph.SetPxPyPzE(ph1->Px(), ph1->Py(),ph1->Pz(),ph1->E());
+ }
+ else continue;
+
+ fakeOmega=fakePi0+vph;
+ for(Int_t ii=0;ii<fNCentBin;ii++){
+ fhFakeOmega[icentbin]->Fill(fakeOmega.Pt(), fakeOmega.M());
+ }
+ }//j
+
+ //continue ................
+ Int_t nMixed = fEventsList[curEventBin]->GetSize();
+ for(Int_t ie=0;ie<nMixed;ie++){
+ TClonesArray* ev2= (TClonesArray*) (fEventsList[curEventBin]->At(ie));
+ for(Int_t mix1=0;mix1<ev2->GetEntries();mix1++){
+ AliAODPWG4Particle *ph2 = (AliAODPWG4Particle*) (ev2->At(mix1));
+ TLorentzVector vph2(ph2->Px(),ph2->Py(),ph2->Pz(),ph2->E());
+ Double_t pi0asy = TMath::Abs(vph1.E()-vph2.E())/(vph1.E()+vph2.E());
+ Double_t pi0mass=(vph1+vph2).M();
+
+ if(TMath::Abs(pi0mass-fPi0Mass)<fPi0MassWindow){//for pi0 selection
+ for(Int_t mix2=(mix1+1);mix2<ev2->GetEntries();mix2++){
+ AliAODPWG4Particle *ph3 = (AliAODPWG4Particle*) (ev2->At(mix2));
+ TLorentzVector vph3(ph3->Px(),ph3->Py(),ph3->Pz(),ph3->E());
+
+ Double_t pi0gammapt=(vph1+vph2+vph3).Pt();
+ Double_t pi0gammamass=(vph1+vph2+vph3).M();
+ Double_t pi0OverOmegaPtRatio =(vph1+vph2).Pt()/pi0gammapt;
+ Double_t gammaOverOmegaPtRatio= vph3.Pt()/pi0gammapt;
+ //pi0, gamma pt cut
+ if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
+ gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
+
+ for(Int_t ipid=0;ipid<fNpid;ipid++){
+ for(Int_t idist=0;idist<fNBadChDistBin;idist++){
+ Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
+ if(ph1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
+ ph2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
+ ph3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
+ ph1->DistToBad()>=idist &&
+ ph2->DistToBad()>=idist &&
+ ph3->DistToBad()>=idist ){
+ if(GetDebug() > 2) printf("MixB: index %d pt %2.3f mass %2.3f \n", index, pi0gammapt, pi0gammamass);
+ //fill histograms
+ fMixBOmega0[index]->Fill(pi0gammapt,pi0gammamass);
+ if(pi0asy<0.7) fMixBOmega1[index]->Fill(pi0gammapt,pi0gammamass);
+ if(pi0asy<0.8) fMixBOmega2[index]->Fill(pi0gammapt,pi0gammamass);
+ //printf("mix B %d %2.2f \n", index, pi0gammamass);
+ }
+ }
+ }
+ }
+
+ //
+ // --C (r1_event1+r2_event2)+r3_event3
+ //
+ if(GetDebug() >0)printf("MixC: (r1_event1+r2_event2)+r3_event3\n");
+ for(Int_t je=(ie+1);je<nMixed;je++){
+ TClonesArray* ev3= (TClonesArray*) (fEventsList[curEventBin]->At(je));
+ for(Int_t mix3=0;mix3<ev3->GetEntries();mix3++){
+ AliAODPWG4Particle *ph3 = (AliAODPWG4Particle*) (ev3->At(mix3));
+ TLorentzVector vph3(ph3->Px(),ph3->Py(),ph3->Pz(),ph3->E());
+
+ Double_t pi0gammapt=(vph1+vph2+vph3).Pt();
+ Double_t pi0gammamass=(vph1+vph2+vph3).M();
+ Double_t pi0OverOmegaPtRatio =(vph1+vph2).Pt()/pi0gammapt;
+ Double_t gammaOverOmegaPtRatio= vph3.Pt()/pi0gammapt;
+ //pi0, gamma pt cut
+ if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
+ gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
+
+ for(Int_t ipid=0;ipid<fNpid;ipid++){
+ for(Int_t idist=0;idist<fNBadChDistBin;idist++){
+ Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
+ if(ph1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
+ ph2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
+ ph3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
+ ph1->DistToBad()>=idist &&
+ ph2->DistToBad()>=idist &&
+ ph3->DistToBad()>=idist ){
+ if(GetDebug() > 2) printf("MixC: index %d pt %2.3f mass %2.3f \n", index, pi0gammapt, pi0gammamass);
+ //fill histograms
+ fMixCOmega0[index]->Fill(pi0gammapt,pi0gammamass);
+ if(pi0asy<0.7) fMixCOmega1[index]->Fill(pi0gammapt,pi0gammamass);
+ if(pi0asy<0.8) fMixCOmega2[index]->Fill(pi0gammapt,pi0gammamass);
+ //printf("mix C %d %2.2f \n", index, pi0gammamass);
+ }
+ }
+ }
+ }
+ }
+ } //for pi0 selecton
+ }
+ }
+ }
+
+
+ //event buffer
+ TClonesArray *currentEvent = new TClonesArray(*aodGamma);
+ if(currentEvent->GetEntriesFast()>0){
+ fEventsList[curEventBin]->AddFirst(currentEvent) ;
+ currentEvent=0 ;
+ if(fEventsList[curEventBin]->GetSize()>=fNmaxMixEv) {
+ TClonesArray * tmp = (TClonesArray*) (fEventsList[curEventBin]->Last()) ;
+ fEventsList[curEventBin]->RemoveLast() ;
+ delete tmp ;
+ }
+ }
+ else{
+ delete currentEvent ;
+ currentEvent=0 ;
+ }
+
+}
+
+//______________________________________________________________________________
+void AliAnaOmegaToPi0Gamma::ReadHistograms(TList * outputList)
+{
+ //read the histograms
+ //for the finalization of the terminate analysis
+
+ Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"RealToPi0Gamma_Vz0C0Rp0Pid0Dist0"));
+
+ Int_t ndim=fNVtxZBin*fNCentBin*fNRpBin*fNBadChDistBin*fNpid;
+
+ if(!fRealOmega0) fRealOmega0 =new TH2F*[ndim];
+ if(!fMixAOmega0) fMixAOmega0 =new TH2F*[ndim];
+ if(!fMixBOmega0) fMixBOmega0 =new TH2F*[ndim];
+ if(!fMixCOmega0) fMixCOmega0 =new TH2F*[ndim];
+
+ if(!fRealOmega1) fRealOmega1 =new TH2F*[ndim];
+ if(!fMixAOmega1) fMixAOmega1 =new TH2F*[ndim];
+ if(!fMixBOmega1) fMixBOmega1 =new TH2F*[ndim];
+ if(!fMixCOmega1) fMixCOmega1 =new TH2F*[ndim];
+
+ if(!fRealOmega2) fRealOmega2 =new TH2F*[ndim];
+ if(!fMixAOmega2) fMixAOmega2 =new TH2F*[ndim];
+ if(!fMixBOmega2) fMixBOmega2 =new TH2F*[ndim];
+ if(!fMixCOmega2) fMixCOmega2 =new TH2F*[ndim];
+
+ for(Int_t i=0;i<fNVtxZBin;i++){
+ for(Int_t j=0;j<fNCentBin;j++){
+ for(Int_t k=0;k<fNRpBin;k++){ //at event level
+ Int_t idim=i*fNCentBin*fNRpBin+j*fNRpBin+k;
+ for(Int_t ipid=0;ipid<fNpid;ipid++){
+ for(Int_t idist=0;idist<fNBadChDistBin;idist++){ //at particle
+ Int_t ind=idim*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
+ fRealOmega0[ind]= (TH2F*) outputList->At(index++);
+ fMixAOmega0[ind]= (TH2F*) outputList->At(index++);
+ fMixBOmega0[ind]= (TH2F*) outputList->At(index++);
+ fMixCOmega0[ind]= (TH2F*) outputList->At(index++);
+
+ fRealOmega1[ind]= (TH2F*) outputList->At(index++);
+ fMixAOmega1[ind]= (TH2F*) outputList->At(index++);
+ fMixBOmega1[ind]= (TH2F*) outputList->At(index++);
+ fMixCOmega1[ind]= (TH2F*) outputList->At(index++);
+
+ fRealOmega2[ind]= (TH2F*) outputList->At(index++);
+ fMixAOmega2[ind]= (TH2F*) outputList->At(index++);
+ fMixBOmega2[ind]= (TH2F*) outputList->At(index++);
+ fMixCOmega2[ind]= (TH2F*) outputList->At(index++);
+
+
+ }
+ }
+ }
+ }
+ }
+
+ if(IsDataMC()){
+ fhOmegaPriPt = (TH1F*) outputList->At(index++);
+ }
+
+}
+
+//______________________________________________________________________________
+void AliAnaOmegaToPi0Gamma::Terminate(TList * outputList)
+{
+// //Do some calculations and plots from the final histograms.
+ if(GetDebug() >= 0) printf("AliAnaOmegaToPi0Gamma::Terminate() \n");
+ ReadHistograms(outputList);
+ const Int_t buffersize = 255;
+ char cvs1[buffersize];
+ snprintf(cvs1, buffersize, "Neutral_%s_IVM",fInputAODGammaName.Data());
+
+ TCanvas * cvsIVM = new TCanvas(cvs1, cvs1, 400, 10, 600, 700) ;
+ cvsIVM->Divide(2, 2);
+
+ cvsIVM->cd(1);
+ char dec[buffersize];
+ snprintf(dec,buffersize,"h2Real_%s",fInputAODGammaName.Data());
+ TH2F * h2Real= (TH2F*)fRealOmega0[0]->Clone(dec);
+ h2Real->GetXaxis()->SetRangeUser(4,6);
+ TH1F * hRealOmega = (TH1F*) h2Real->ProjectionY();
+ hRealOmega->SetTitle("RealPi0Gamma 4<pt<6");
+ hRealOmega->SetLineColor(2);
+ hRealOmega->Draw();
+
+ cvsIVM->cd(2);
+ snprintf(dec,buffersize,"hMixA_%s",fInputAODGammaName.Data());
+ TH2F *h2MixA= (TH2F*)fMixAOmega0[0]->Clone(dec);
+ h2MixA->GetXaxis()->SetRangeUser(4,6);
+ TH1F * hMixAOmega = (TH1F*) h2MixA->ProjectionY();
+ hMixAOmega->SetTitle("MixA 4<pt<6");
+ hMixAOmega->SetLineColor(2);
+ hMixAOmega->Draw();
+
+ cvsIVM->cd(3);
+ snprintf(dec,buffersize,"hMixB_%s",fInputAODGammaName.Data());
+ TH2F * h2MixB= (TH2F*)fMixBOmega0[0]->Clone(dec);
+ h2MixB->GetXaxis()->SetRangeUser(4,6);
+ TH1F * hMixBOmega = (TH1F*) h2MixB->ProjectionY();
+ hMixBOmega->SetTitle("MixB 4<pt<6");
+ hMixBOmega->SetLineColor(2);
+ hMixBOmega->Draw();
+
+ cvsIVM->cd(4);
+ snprintf(dec,buffersize,"hMixC_%s",fInputAODGammaName.Data());
+ TH2F *h2MixC= (TH2F*)fMixCOmega0[0]->Clone(dec);
+ h2MixC->GetXaxis()->SetRangeUser(4,6);
+ TH1F * hMixCOmega = (TH1F*) h2MixC->ProjectionY();
+ hMixCOmega->SetTitle("MixC 4<pt<6");
+ hMixCOmega->SetLineColor(2);
+ hMixCOmega->Draw();
+
+ char eps[buffersize];
+ snprintf(eps,buffersize,"CVS_%s_IVM.eps",fInputAODGammaName.Data());
+ cvsIVM->Print(eps);
+ cvsIVM->Modified();
+
+}