]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
changed methods in AliAODPWG4Particle (G. Conesa)
authorkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 25 Jul 2011 17:52:43 +0000 (17:52 +0000)
committerkleinb <kleinb@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 25 Jul 2011 17:52:43 +0000 (17:52 +0000)
PWG4/PartCorrBase/AliCaloPID.cxx
PWG4/PartCorrDep/AliAnaExample.cxx
PWG4/PartCorrDep/AliAnaOmegaToPi0Gamma.cxx
PWG4/PartCorrDep/AliAnaPhoton.cxx
PWG4/PartCorrDep/AliAnaPi0EbE.cxx
STEER/AOD/AliAODPWG4Particle.h

index f768e7d2f774df05a7b607797e04d217fc37f7ac..b5dd16458a763be36ec163c19a4d39545972a98e 100755 (executable)
 
 //_________________________________________________________________________
 // Class for PID selection with calorimeters
-// The Output of the 2 main methods GetPdg is a PDG number identifying the cluster, 
+// The Output of the 2 main methods GetIdentifiedParticleType is a PDG number identifying the cluster, 
 // being kPhoton, kElectron, kPi0 ... as defined in the header file
-//   - GetPdg(const TString calo, const Double_t * pid, const Float_t energy)
+//   - GetIdentifiedParticleType(const TString calo, const Double_t * pid, const Float_t energy)
 //      Reads the PID weights array of the ESDs and depending on its magnitude identifies the particle
-//   - GetPdg(const TString calo,const TLorentzVector mom, const AliVCluster * cluster)
+//   - GetIdentifiedParticleType(const TString calo,const TLorentzVector mom, const AliVCluster * cluster)
 //      Recalcultes PID, the bayesian or any new one to be implemented in the future
 //      Right now only the possibility to recalculate EMCAL with bayesian and simple PID.
 //      In order to recalculate Bayesian, it is necessary to load the EMCALUtils library
@@ -212,11 +212,11 @@ void AliCaloPID::InitParameters()
 }
 
 //_______________________________________________________________
-Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t energy) {
+Int_t AliCaloPID::GetIdentifiedParticleType(const TString calo, const Double_t * pid, const Float_t energy) {
   //Return most probable identity of the particle.
   
   if(!pid){ 
-    printf("AliCaloPID::GetPdg() - pid pointer not initialized!!!\n");
+    printf("AliCaloPID::GetIdentifiedParticleType() - pid pointer not initialized!!!\n");
     abort();
   }
   
@@ -241,7 +241,7 @@ Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t
     
   }
   
-  if(fDebug > 0)  printf("AliCaloPID::GetPdg: calo %s, ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, hadrons: pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
+  if(fDebug > 0)  printf("AliCaloPID::GetIdentifiedParticleType: calo %s, ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, hadrons: pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
                         calo.Data(),pid[AliVCluster::kPhoton], pid[AliVCluster::kPi0],
                         pid[AliVCluster::kElectron], pid[AliVCluster::kEleCon],
                         pid[AliVCluster::kPion], pid[AliVCluster::kKaon], pid[AliVCluster::kProton],
@@ -278,20 +278,20 @@ Int_t AliCaloPID::GetPdg(const TString calo, const Double_t * pid, const Float_t
     else                                                     pdg = kNeutralUnknown ;
   }
   
-  if(fDebug > 0)printf("AliCaloPID::GetPdg:Final Pdg: %d, cluster energy %2.2f \n", pdg,energy);
+  if(fDebug > 0)printf("AliCaloPID::GetIdentifiedParticleType:Final Pdg: %d, cluster energy %2.2f \n", pdg,energy);
    //printf("PDG1\n");
   return pdg ;
   
 }
 
 //_______________________________________________________________
-Int_t AliCaloPID::GetPdg(const TString calo,const TLorentzVector mom, const AliVCluster * cluster) {
+Int_t AliCaloPID::GetIdentifiedParticleType(const TString calo,const TLorentzVector mom, const AliVCluster * cluster) {
   //Recalculated PID with all parameters
   
   Float_t lambda0 = cluster->GetM02();
   Float_t energy  = mom.E();   
   
-  if(fDebug > 0) printf("AliCaloPID::GetPdg: Calorimeter %s, E %3.2f, l0 %3.2f, l1 %3.2f, disp %3.2f, tof %1.11f, distCPV %3.2f, distToBC %1.1f, NMax %d\n",
+  if(fDebug > 0) printf("AliCaloPID::GetIdentifiedParticleType: Calorimeter %s, E %3.2f, l0 %3.2f, l1 %3.2f, disp %3.2f, tof %1.11f, distCPV %3.2f, distToBC %1.1f, NMax %d\n",
                         calo.Data(),energy,lambda0,cluster->GetM20(),cluster->GetDispersion(),cluster->GetTOF(), 
                         cluster->GetEmcCpvDistance(), cluster->GetDistanceToBadChannel(),cluster->GetNExMax());
   
@@ -300,7 +300,7 @@ Int_t AliCaloPID::GetPdg(const TString calo,const TLorentzVector mom, const AliV
          if(fRecalculateBayesian){       
                  if(fDebug > 0)  {
                          const Double_t  *pid0 = cluster->GetPID();
-                         printf("AliCaloPID::GetPdg: BEFORE calo %s, ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, hadrons: pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
+                         printf("AliCaloPID::GetIdentifiedParticleType: BEFORE calo %s, ph %0.2f, pi0 %0.2f, el %0.2f, conv el %0.2f, hadrons: pion %0.2f, kaon %0.2f, proton %0.2f , neutron %0.2f, kaon %0.2f \n",
                calo.Data(),pid0[AliVCluster::kPhoton], pid0[AliVCluster::kPi0],
                pid0[AliVCluster::kElectron], pid0[AliVCluster::kEleCon],
                pid0[AliVCluster::kPion], pid0[AliVCluster::kKaon], pid0[AliVCluster::kProton],
@@ -310,18 +310,15 @@ Int_t AliCaloPID::GetPdg(const TString calo,const TLorentzVector mom, const AliV
       fEMCALPIDUtils->ComputePID(energy, lambda0);
       Double_t pid[AliPID::kSPECIESN];
       for(Int_t i = 0; i < AliPID::kSPECIESN; i++) pid[i] = fEMCALPIDUtils->GetPIDFinal(i);
-      return GetPdg(calo, pid, energy);
+      return GetIdentifiedParticleType(calo, pid, energy);
                  
     }
   }
   
   // If no use of bayesian, simplest PID  
   if(lambda0 < fDispCut) return kPhoton ;
-  //else return  kNeutralHadron ; 
-  else return  kPi0 ;
-  
-  return  kNeutralHadron ; 
-  
+  else return  kPi0 ; // Wild guess, it can be hadron but not photon
+    
 }
 
 //__________________________________________________
@@ -445,12 +442,12 @@ void AliCaloPID::SetPIDBits(const TString calo, const AliVCluster * cluster, Ali
   ph->SetChargedBit(isNeutral);
   
   //Set PID pdg
-  ph->SetPdg(GetPdg(calo,cluster->GetPID(),ph->E()));
+  ph->SetIdentifiedParticleType(GetIdentifiedParticleType(calo,cluster->GetPID(),ph->E()));
   
   if(fDebug > 0){ 
     printf("AliCaloPID::SetPIDBits: TOF %e, Lambda0 %2.2f, Lambda1 %2.2f\n",tof , l0, l1);     
     printf("AliCaloPID::SetPIDBits: pdg %d, bits: TOF %d, Dispersion %d, Charge %d\n",
-            ph->GetPdg(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit()); 
+            ph->GetIdentifiedParticleType(), ph->GetTOFBit() , ph->GetDispBit() , ph->GetChargedBit()); 
   }
 }
 
index 4abbd5be1105c95204c3bda2d1a076382abd74f4..b0a2decca526a61270eccb8f3e817c79be8cbde9 100755 (executable)
@@ -186,8 +186,8 @@ void  AliAnaExample::MakeAnalysisFillAOD()
       
       if(IsCaloPIDOn()){
         const Double_t *pid = calo->GetPID();
-        pdg = GetCaloPID()->GetPdg(fDetector,pid,mom.E());
-        //pdg = GetCaloPID()->GetPdg(fDetector,mom,
+        pdg = GetCaloPID()->GetIdentifiedParticleType(fDetector,pid,mom.E());
+        //pdg = GetCaloPID()->GetIdentifiedParticleType(fDetector,mom,
         //               calo->GetM02(), calo->GetM02(),
         //               calo->GetDispersion(), 0, 0); 
       }
@@ -204,7 +204,7 @@ void  AliAnaExample::MakeAnalysisFillAOD()
         AliAODPWG4Particle ph = AliAODPWG4Particle(mom);
         //AddAODParticleCorrelation(AliAODPWG4ParticleCorrelation(mom));
         ph.SetLabel(calo->GetLabel());
-        ph.SetPdg(pdg);
+        ph.SetIdentifiedParticleType(pdg);
         ph.SetDetector(fDetector);
         AddAODParticle(ph);
       }//selection
index 63cdeac995a084b47464fa788a0cbe0bdb78455f..7c5bf69c72ab6ee848f1c833e1445d9fc03be907 100644 (file)
-/**************************************************************************\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();
+}
index 8580b30bd6fb849219a2c05c0b7aaf8047aa57d2..4cfd0bccfdf397afffc3adcc598fcc4dc22e1130 100755 (executable)
@@ -1029,10 +1029,10 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
     // MC
     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
       //Get most probable PID, check PID weights (in MC this option is mandatory)
-      aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
-      if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());        
+      aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
+      if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());     
       //If primary is not photon, skip it.
-      if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;
+      if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;
     }  
     //...............................................
     // Data, PID check on
@@ -1040,14 +1040,14 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
       //Get most probable PID, 2 options check PID weights 
       //or redo PID, recommended option for EMCal.             
       if(!IsCaloPIDRecalculationOn())
-        aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
+        aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
       else
-        aodph.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
+        aodph.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
       
-      if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetPdg());
+      if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PDG of identified particle %d\n",aodph.GetIdentifiedParticleType());
       
       //If cluster does not pass pid, not photon, skip it.
-      if(aodph.GetPdg() != AliCaloPID::kPhoton) continue ;                     
+      if(aodph.GetIdentifiedParticleType() != AliCaloPID::kPhoton) continue ;                  
       
     }
     //...............................................
@@ -1059,7 +1059,7 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
       if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - PID Bits set \n");              
     }
     
-    if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetPdg());
+    if(GetDebug() > 1) printf("AliAnaPhoton::MakeAnalysisFillAOD() - Photon selection cuts passed: pT %3.2f, pdg %d\n",aodph.Pt(), aodph.GetIdentifiedParticleType());
     
     //--------------------------------------------------------------------------------------
     //Play with the MC stack if available
@@ -1300,7 +1300,7 @@ void  AliAnaPhoton::MakeAnalysisFillAOD()
           //Set the indeces of the original caloclusters  
           aodpair.SetCaloLabel(calo->GetID(),id2);
           aodpair.SetDetector(fCalorimeter);
-          aodpair.SetPdg(aodph.GetPdg());
+          aodpair.SetIdentifiedParticleType(aodph.GetIdentifiedParticleType());
           aodpair.SetTag(aodph.GetTag());
           aodpair.SetTagged(kTRUE);
           //Add AOD with pair object to aod branch
@@ -1406,10 +1406,10 @@ void  AliAnaPhoton::MakeAnalysisFillHistograms()
        
        for(Int_t iaod = 0; iaod < naod ; iaod++){
          AliAODPWG4Particle* ph =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
-         Int_t pdg = ph->GetPdg();
+         Int_t pdg = ph->GetIdentifiedParticleType();
          
          if(GetDebug() > 3) 
-           printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetPdg(),ph->GetTag(), (ph->GetDetector()).Data()) ;
+           printf("AliAnaPhoton::MakeAnalysisFillHistograms() - PDG %d, MC TAG %d, Calorimeter %s\n", ph->GetIdentifiedParticleType(),ph->GetTag(), (ph->GetDetector()).Data()) ;
          
          //If PID used, fill histos with photons in Calorimeter fCalorimeter
          if(IsCaloPIDOn() && pdg != AliCaloPID::kPhoton) continue; 
index 14e7fec8490254326e59b1a40d75c7ffde2dce98..c8166c050d47b0338016a38d80171393b93c6446 100755 (executable)
@@ -5,7 +5,7 @@
  * 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 iGetEntriesFast(s hereby granted   *
+ * 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     *
@@ -19,7 +19,7 @@
 // Pi0 identified by one of the following:
 //  -Invariant mass of 2 cluster in calorimeter
 //  -Shower shape analysis in calorimeter
-//  -Invariant mass of one cluster in calorimeter and one photon reconstructed in TPC (in near future)
+//  -Invariant mass of one cluster in calorimeter and one photon reconstructed in CTS
 //
 // -- Author: Gustavo Conesa (LNF-INFN) &  Raphaelle Ichou (SUBATECH)
 //////////////////////////////////////////////////////////////////////////////
@@ -588,7 +588,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeter()
         mom = mom1+mom2;
         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
         //pi0.SetLabel(calo->GetLabel());
-        pi0.SetPdg(AliCaloPID::kPi0);
+        pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
         pi0.SetDetector(photon1->GetDetector());
         pi0.SetTag(tag);  
         //Set the indeces of the original caloclusters  
@@ -722,7 +722,7 @@ void  AliAnaPi0EbE::MakeInvMassInCalorimeterAndCTS()
         mom = mom1+mom2;
         AliAODPWG4Particle pi0 = AliAODPWG4Particle(mom);
         //pi0.SetLabel(calo->GetLabel());
-        pi0.SetPdg(AliCaloPID::kPi0);
+        pi0.SetIdentifiedParticleType(AliCaloPID::kPi0);
         pi0.SetDetector(photon1->GetDetector());
         pi0.SetTag(tag);
         //Set the indeces of the original tracks or caloclusters  
@@ -825,11 +825,11 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
     //PID selection or bit setting
     if(GetReader()->GetDataType() == AliCaloTrackReader::kMC){
       //Get most probable PID, check PID weights (in MC this option is mandatory)
-      aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
+      aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
       if(GetDebug() > 1) 
-        printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: PDG of identified particle %d\n",aodpi0.GetPdg());
+        printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - FillAOD: PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
       //If primary is not pi0, skip it.
-      if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;
+      if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;
     }                                  
     else if(IsCaloPIDOn()){
       //Skip matched clusters with tracks
@@ -838,14 +838,14 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
       //Get most probable PID, 2 options check PID weights 
       //or redo PID, recommended option for EMCal.             
       if(!IsCaloPIDRecalculationOn())
-        aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
+        aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,calo->GetPID(),mom.E()));//PID with weights
       else
-        aodpi0.SetPdg(GetCaloPID()->GetPdg(fCalorimeter,mom,calo));//PID recalculated
+        aodpi0.SetIdentifiedParticleType(GetCaloPID()->GetIdentifiedParticleType(fCalorimeter,mom,calo));//PID recalculated
       
-      if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetPdg());
+      if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PDG of identified particle %d\n",aodpi0.GetIdentifiedParticleType());
       
       //If cluster does not pass pid, not pi0, skip it.
-      if(aodpi0.GetPdg() != AliCaloPID::kPi0) continue ;                       
+      if(aodpi0.GetIdentifiedParticleType() != AliCaloPID::kPi0) continue ;                    
       
     }
     else{
@@ -855,7 +855,7 @@ void  AliAnaPi0EbE::MakeShowerShapeIdentification()
       if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - PID Bits set \n");            
     }
     
-    if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetPdg());
+    if(GetDebug() > 1) printf("AliAnaPi0EbE::MakeShowerShapeIdentification() - Pi0 selection cuts passed: pT %3.2f, pdg %d\n",aodpi0.Pt(), aodpi0.GetIdentifiedParticleType());
     
     //Play with the MC stack if available
     //Check origin of the candidates
@@ -895,7 +895,7 @@ void  AliAnaPi0EbE::MakeAnalysisFillHistograms()
   for(Int_t iaod = 0; iaod < naod ; iaod++){
     
     AliAODPWG4Particle* pi0 =  (AliAODPWG4Particle*) (GetOutputAODBranch()->At(iaod));
-    Int_t pdg = pi0->GetPdg();
+    Int_t pdg = pi0->GetIdentifiedParticleType();
          
     if(IsCaloPIDOn() && pdg != AliCaloPID::kPi0) continue;              
     
index 0fdbe511bc714f445f8ac74a005e1e7ee38ef335..a3b0c2e8a442015dcf8c9daeac511903b1b1a05d 100755 (executable)
-#ifndef ALIAODPWG4PARTICLE_H\r
-#define ALIAODPWG4PARTICLE_H\r
-/* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *\r
- * See cxx source for full Copyright notice                               */\r
-\r
-/* $Id: AliAODPWG4Particle.h  $ */\r
-\r
-//-------------------------------------------------------------------------\r
-//     Copy of AOD photon class, adapted for particle identification\r
-//     and correlations analysis\r
-//     Author: Yves Schutz, CERN, Gustavo Conesa, INFN\r
-//-------------------------------------------------------------------------\r
-\r
-//-- ROOT system --\r
-#include <TLorentzVector.h>\r
-class TString;\r
-\r
-//-- Analysis system\r
-#include "AliVParticle.h"\r
-\r
-class AliAODPWG4Particle : public AliVParticle {\r
-  \r
- public:\r
-  AliAODPWG4Particle();\r
-  AliAODPWG4Particle(Double_t px, Double_t py, Double_t pz, Double_t e);\r
-  AliAODPWG4Particle(TLorentzVector & p);\r
-  \r
-  virtual ~AliAODPWG4Particle();\r
-  virtual void Clear(const Option_t* /*opt*/);\r
-\r
-  AliAODPWG4Particle(const AliAODPWG4Particle& photon); \r
-  AliAODPWG4Particle& operator=(const AliAODPWG4Particle& photon);\r
-\r
-  //enumerated type for various b-tags of electrons\r
-  enum btagTypes {kDVMTag0, kDVMTag1, kDVMTag2, kTransverseIPTag, kUnknownTag};\r
-\r
-  void SetBTagBit(Int_t &tag, const UInt_t set) const {\r
-    // Set bit of type set (btagTypes) in tag                                             \r
-    tag |= (1<<set) ;\r
-  }\r
-\r
-  Bool_t CheckBTagBit(const Int_t tag, const UInt_t test) const {\r
-    // Check if in fBtag the bit test (btagTypes) is set.                                   \r
-    if (tag & (1<<test) ) return  kTRUE ;\r
-    else return kFALSE ;\r
-  }\r
-\r
-  // AliVParticle methods\r
-  virtual Double_t Px()         const { return fMomentum->Px();      }\r
-  virtual Double_t Py()         const { return fMomentum->Py();      }\r
-  virtual Double_t Pz()         const { return fMomentum->Pz();      }\r
-  virtual Double_t Pt()         const { return fMomentum->Pt();      }\r
-  virtual Double_t P()          const { return fMomentum->P();       }\r
-  virtual Bool_t   PxPyPz(Double_t p[3]) const { p[0] = Px(); p[1] = Py(); p[2] = Pz(); return kTRUE; }\r
-  virtual Double_t OneOverPt()  const { return 1. / fMomentum->Pt(); }\r
-  virtual Double_t Phi()        const;\r
-  virtual Double_t Theta()      const { return fMomentum->Theta();   }\r
-  virtual Double_t E()          const { return fMomentum->E();       }\r
-  virtual Double_t M()          const { return fMomentum->M();       }\r
-  virtual Double_t Eta()        const { return fMomentum->Eta();     }\r
-  virtual Double_t Y()          const { return fMomentum->Rapidity();}\r
-  virtual Double_t Xv()         const {return -999.;} // put reasonable values here\r
-  virtual Double_t Yv()         const {return -999.;} //\r
-  virtual Double_t Zv()         const {return -999.;} //\r
-  virtual Bool_t   XvYvZv(Double_t x[3]) const { x[0] = Xv(); x[1] = Yv(); x[2] = Zv(); return kTRUE; }  \r
-  virtual void     Print(Option_t* /*option*/) const;\r
-  \r
-  //\r
-  //Dummy\r
-  virtual Short_t Charge()      const { return 0;}\r
-  virtual const Double_t* PID() const { return NULL;}\r
-  //\r
-  \r
-  virtual Int_t   GetPdg()               const {return fPdg ; }\r
-  virtual Int_t   GetTag()               const {return fTag ; }\r
-  virtual Int_t   GetBtag()              const {return fBtag ; }\r
-  virtual Int_t   GetLabel()             const {return fLabel ; }\r
-  virtual Int_t   GetCaloLabel (Int_t i) const {return fCaloLabel[i] ; }\r
-  virtual Int_t   GetTrackLabel(Int_t i) const {return fTrackLabel[i] ; }\r
-  virtual TString GetDetector()          const {return fDetector ; }\r
-  virtual Bool_t  GetDispBit(void)       const {return fDisp ; }\r
-  virtual Bool_t  GetTOFBit(void)        const {return fTof ; }\r
-  virtual Bool_t  GetChargedBit(void)    const {return fCharged ; }\r
-  virtual Int_t   DistToBad()            const {return fBadDist ; } \r
-  virtual Int_t   GetInputFileIndex()    const {return fInputFileIndex ; }\r
-  virtual Int_t   GetFiducialArea()      const {return fFidArea;}\r
-\r
-  virtual Bool_t   IsTagged()    const {return fTagged ; }\r
\r
-  virtual void SetLabel(Int_t l)         { fLabel = l ; }\r
-  virtual void SetCaloLabel (Int_t a, Int_t b) { fCaloLabel [0] = a; fCaloLabel [1] = b  ; }\r
-  virtual void SetTrackLabel(Int_t a, Int_t b) { fTrackLabel[0] = a; fTrackLabel[1] = b  ; }\r
-  virtual void SetTrackLabel(Int_t a, Int_t b, Int_t c, Int_t d) \r
-  { fTrackLabel[0] = a; fTrackLabel[1] = b  ; fTrackLabel[2] = c; fTrackLabel[3] = d; }\r
-  \r
-  virtual void SetPdg(Int_t pdg)         { fPdg = pdg ; }\r
-  virtual void SetTag(Int_t tag)         { fTag = tag ; }\r
-  virtual void SetTagged(Bool_t tag)     { fTagged = tag ; }\r
-  virtual void SetBtag(Int_t tag)        { fBtag = tag ; }\r
-  virtual void SetDetector(TString d)    { fDetector = d ; }\r
-  virtual void SetDispBit(Bool_t disp)   { fDisp = disp ; } \r
-  virtual void SetTOFBit(Bool_t tof)     { fTof = tof ;} \r
-  virtual void SetChargedBit(Bool_t ch)  { fCharged = ch ; }\r
-  virtual void SetDistToBad(Int_t dist)  { fBadDist=dist ;} \r
-  virtual void SetInputFileIndex(Int_t i){ fInputFileIndex = i;}\r
-  virtual void SetFiducialArea(Int_t a)  {fFidArea=a;}\r
-       \r
-  TLorentzVector * Momentum() const                { return fMomentum ; }\r
-  virtual void     SetMomentum(TLorentzVector *lv) { fMomentum = lv ; }\r
-  \r
-  Bool_t IsPIDOK(const Int_t ipid, const Int_t pdgwanted) const;\r
-  Double_t GetPairMass(AliAODPWG4Particle * p)const{ return (*(p->fMomentum)+*fMomentum).M(); }\r
-  // Dummy\r
-  Int_t    PdgCode() const {return 0;}\r
-  \r
- private:\r
-  TLorentzVector* fMomentum;  // Photon 4-momentum vector\r
-  Int_t      fPdg ;          // id of particle\r
-  Int_t      fTag ;          // tag of particle (decay, fragment, prompt photon), MC\r
-  Int_t      fBtag;          // tag particle from B.\r
-  Int_t      fLabel ;        // MC label\r
-  Int_t      fCaloLabel[2];  // CaloCluster index, 1 for photons, 2 for pi0.\r
-  Int_t      fTrackLabel[4]; // Track lable, 1 for pions, 2 for conversion photons \r
-  TString    fDetector ;     // Detector where particle was measured.\r
-  Bool_t     fDisp ;         // Dispersion bit\r
-  Bool_t     fTof ;          // TOF bit\r
-  Bool_t     fCharged ;      // Charged bit\r
-  Bool_t     fTagged ;       // If photon tagged (pi0 decay)\r
-  Int_t      fBadDist ;      // Distance to bad module in module units\r
-  Int_t      fFidArea ;      // Type of fiducial area hit by this photon\r
-  Int_t      fInputFileIndex;// 0, standard input, 1 first input added. \r
-                                // Only possible one for now, more in future?\r
-       \r
-  ClassDef(AliAODPWG4Particle, 4);\r
-};\r
-\r
-inline Double_t AliAODPWG4Particle::Phi() const\r
-{\r
-  // Return phi\r
-  Double_t phi = fMomentum->Phi();\r
-  if (phi < 0.) phi += TMath::TwoPi();\r
-  return phi;\r
-}\r
-\r
-#endif //ALIAODPWG4PARTICLE_H\r
+#ifndef ALIAODPWG4PARTICLE_H
+#define ALIAODPWG4PARTICLE_H
+/* Copyright(c) 1998-2007, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id: AliAODPWG4Particle.h  $ */
+
+//-------------------------------------------------------------------------
+//     Copy of AOD photon class, adapted for particle identification
+//     and correlations analysis
+//     Author: Yves Schutz, CERN, Gustavo Conesa, INFN
+//-------------------------------------------------------------------------
+
+//-- ROOT system --
+#include <TLorentzVector.h>
+class TString;
+
+//-- Analysis system
+#include "AliVParticle.h"
+
+class AliAODPWG4Particle : public AliVParticle {
+  
+ public:
+  AliAODPWG4Particle();
+  AliAODPWG4Particle(Double_t px, Double_t py, Double_t pz, Double_t e);
+  AliAODPWG4Particle(TLorentzVector & p);
+  
+  virtual ~AliAODPWG4Particle();
+  virtual void Clear(const Option_t* /*opt*/);
+
+  AliAODPWG4Particle(const AliAODPWG4Particle& photon); 
+  AliAODPWG4Particle& operator=(const AliAODPWG4Particle& photon);
+
+  //enumerated type for various b-tags of electrons
+  enum btagTypes {kDVMTag0, kDVMTag1, kDVMTag2, kTransverseIPTag, kUnknownTag};
+
+  void SetBTagBit(Int_t &tag, const UInt_t set) const {
+    // Set bit of type set (btagTypes) in tag                                             
+    tag |= (1<<set) ;
+  }
+
+  Bool_t CheckBTagBit(const Int_t tag, const UInt_t test) const {
+    // Check if in fBtag the bit test (btagTypes) is set.                                   
+    if (tag & (1<<test) ) return  kTRUE ;
+    else return kFALSE ;
+  }
+
+  // AliVParticle methods
+  virtual Double_t Px()         const { return fMomentum->Px();      }
+  virtual Double_t Py()         const { return fMomentum->Py();      }
+  virtual Double_t Pz()         const { return fMomentum->Pz();      }
+  virtual Double_t Pt()         const { return fMomentum->Pt();      }
+  virtual Double_t P()          const { return fMomentum->P();       }
+  virtual Bool_t   PxPyPz(Double_t p[3]) const { p[0] = Px(); p[1] = Py(); p[2] = Pz(); return kTRUE; }
+  virtual Double_t OneOverPt()  const { return 1. / fMomentum->Pt(); }
+  virtual Double_t Phi()        const;
+  virtual Double_t Theta()      const { return fMomentum->Theta();   }
+  virtual Double_t E()          const { return fMomentum->E();       }
+  virtual Double_t M()          const { return fMomentum->M();       }
+  virtual Double_t Eta()        const { return fMomentum->Eta();     }
+  virtual Double_t Y()          const { return fMomentum->Rapidity();}
+  virtual Double_t Xv()         const {return -999.;} // put reasonable values here
+  virtual Double_t Yv()         const {return -999.;} //
+  virtual Double_t Zv()         const {return -999.;} //
+  virtual Bool_t   XvYvZv(Double_t x[3]) const { x[0] = Xv(); x[1] = Yv(); x[2] = Zv(); return kTRUE; }  
+  virtual void     Print(Option_t* /*option*/) const;
+  
+  //
+  //Dummy
+  virtual Short_t Charge()      const { return 0;}
+  virtual const Double_t* PID() const { return NULL;}
+  //
+  
+  virtual Int_t   GetIdentifiedParticleType() const {return fPdg ; }
+  virtual Int_t   GetTag()               const {return fTag ; }
+  virtual Int_t   GetBtag()              const {return fBtag ; }
+  virtual Int_t   GetLabel()             const {return fLabel ; }
+  virtual Int_t   GetCaloLabel (Int_t i) const {return fCaloLabel[i] ; }
+  virtual Int_t   GetTrackLabel(Int_t i) const {return fTrackLabel[i] ; }
+  virtual TString GetDetector()          const {return fDetector ; }
+  virtual Bool_t  GetDispBit(void)       const {return fDisp ; }
+  virtual Bool_t  GetTOFBit(void)        const {return fTof ; }
+  virtual Bool_t  GetChargedBit(void)    const {return fCharged ; }
+  virtual Int_t   DistToBad()            const {return fBadDist ; } 
+  virtual Int_t   GetInputFileIndex()    const {return fInputFileIndex ; }
+  virtual Int_t   GetFiducialArea()      const {return fFidArea;}
+
+  virtual Bool_t   IsTagged()    const {return fTagged ; }
+  virtual void SetLabel(Int_t l)         { fLabel = l ; }
+  virtual void SetCaloLabel (Int_t a, Int_t b) { fCaloLabel [0] = a; fCaloLabel [1] = b  ; }
+  virtual void SetTrackLabel(Int_t a, Int_t b) { fTrackLabel[0] = a; fTrackLabel[1] = b  ; }
+  virtual void SetTrackLabel(Int_t a, Int_t b, Int_t c, Int_t d) 
+  { fTrackLabel[0] = a; fTrackLabel[1] = b  ; fTrackLabel[2] = c; fTrackLabel[3] = d; }
+  
+  virtual void SetIdentifiedParticleType(Int_t pdg) { fPdg = pdg ; }
+  virtual void SetTag(Int_t tag)         { fTag = tag ; }
+  virtual void SetTagged(Bool_t tag)     { fTagged = tag ; }
+  virtual void SetBtag(Int_t tag)        { fBtag = tag ; }
+  virtual void SetDetector(TString d)    { fDetector = d ; }
+  virtual void SetDispBit(Bool_t disp)   { fDisp = disp ; } 
+  virtual void SetTOFBit(Bool_t tof)     { fTof = tof ;} 
+  virtual void SetChargedBit(Bool_t ch)  { fCharged = ch ; }
+  virtual void SetDistToBad(Int_t dist)  { fBadDist=dist ;} 
+  virtual void SetInputFileIndex(Int_t i){ fInputFileIndex = i;}
+  virtual void SetFiducialArea(Int_t a)  {fFidArea=a;}
+       
+  TLorentzVector * Momentum() const                { return fMomentum ; }
+  virtual void     SetMomentum(TLorentzVector *lv) { fMomentum = lv ; }
+  
+  Bool_t IsPIDOK(const Int_t ipid, const Int_t pdgwanted) const;
+  Double_t GetPairMass(AliAODPWG4Particle * p)const{ return (*(p->fMomentum)+*fMomentum).M(); }
+  // Dummy
+  Int_t    PdgCode() const {return 0;}
+  
+ private:
+  TLorentzVector* fMomentum;  // Photon 4-momentum vector
+  Int_t      fPdg ;          // type of identified particle, same code as PDG, but this is not a MonteCarlo particle 
+  Int_t      fTag ;          // tag of particle (decay, fragment, prompt photon), MC
+  Int_t      fBtag;          // tag particle from B.
+  Int_t      fLabel ;        // MC label
+  Int_t      fCaloLabel[2];  // CaloCluster index, 1 for photons, 2 for pi0.
+  Int_t      fTrackLabel[4]; // Track lable, 1 for pions, 2 for conversion photons 
+  TString    fDetector ;     // Detector where particle was measured.
+  Bool_t     fDisp ;         // Dispersion bit
+  Bool_t     fTof ;          // TOF bit
+  Bool_t     fCharged ;      // Charged bit
+  Bool_t     fTagged ;       // If photon tagged (pi0 decay)
+  Int_t      fBadDist ;      // Distance to bad module in module units
+  Int_t      fFidArea ;      // Type of fiducial area hit by this photon
+  Int_t      fInputFileIndex;// 0, standard input, 1 first input added. 
+                                  // Only possible one for now, more in future?
+       
+  ClassDef(AliAODPWG4Particle, 4);
+};
+
+inline Double_t AliAODPWG4Particle::Phi() const
+{
+  // Return phi
+  Double_t phi = fMomentum->Phi();
+  if (phi < 0.) phi += TMath::TwoPi();
+  return phi;
+}
+
+#endif //ALIAODPWG4PARTICLE_H