PWG4/PartCorrDep/AliAnaOmegaToPi0Gamma.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15 /* $Id: $ */
  16 //_________________________________________________________________________
  17 // class to extract omega(782)->pi0+gamma->3gamma
  18 //  Mar. 22, 2011: Additional method, espeically for EMCAL. A high E cluster is assumpted as pi0 (two photons are overlapped) without unfolding
  19 //
  20 //-- Author: Renzhuo Wan (IOPP-Wuhan, China)
  21 //_________________________________________________________________________
  22
  23 // --- ROOT system
  24 class TROOT;
  25
  26 // --- AliRoot system
  27 //class AliVEvent;
  28 // --- ROOT system ---
  29 #include "TH2F.h"
  30 #include "TLorentzVector.h"
  31 #include "TParticle.h"
  32 #include "TCanvas.h"
  33 #include "TFile.h"
  34 //---- AliRoot system ----
  35 #include "AliAnaOmegaToPi0Gamma.h"
  36 #include "AliCaloTrackReader.h"
  37 #include "AliCaloPID.h"
  38 #include "AliStack.h"
  39 #include "AliVEvent.h"
  40 #include "AliAODEvent.h"
  41 #include "AliAODMCParticle.h"
  42 ClassImp(AliAnaOmegaToPi0Gamma)
  43
  44 //______________________________________________________________________________
  45 AliAnaOmegaToPi0Gamma::AliAnaOmegaToPi0Gamma() : AliAnaPartCorrBaseClass(),
  46 fInputAODPi0(0), fInputAODGammaName(""),
  47 fEventsList(0x0),fNVtxZBin(0), fNCentBin(0), fNRpBin(0), fNBadChDistBin(0), fNpid(0),
  48 fNmaxMixEv(0), fVtxZCut(0), fCent(0), fRp(0),
  49 fPi0Mass(0),fPi0MassWindow(0),fPi0OverOmegaPtCut(0),
  50 fGammaOverOmegaPtCut(0), fEOverlapCluster(0),
  51 fhEtalon(0),
  52 fRealOmega0(0), fMixAOmega0(0),
  53 fMixBOmega0(0), fMixCOmega0(0),
  54 fRealOmega1(0), fMixAOmega1(0),
  55 fMixBOmega1(0), fMixCOmega1(0),
  56 fRealOmega2(0), fMixAOmega2(0),
  57 fMixBOmega2(0), fMixCOmega2(0),
  58 fhFakeOmega(0),
  59 fhOmegaPriPt(0)
  60 {
  61  //Default Ctor
  62  InitParameters();
  63 }
  64 /*
  65 //______________________________________________________________________________
  66 AliAnaOmegaToPi0Gamma::AliAnaOmegaToPi0Gamma(const AliAnaOmegaToPi0Gamma & ex) : AliAnaPartCorrBaseClass(ex),
  67 fInputAODPi0(new TClonesArray (*ex.fInputAODPi0)),
  68 fInputAODGammaName(ex.fInputAODGammaName),
  69 fEventsList(0x0),
  70 fNVtxZBin(ex.fNVtxZBin), fNCentBin(ex.fNCentBin), fNRpBin(ex.fNRpBin),
  71 fNBadChDistBin(ex.fNBadChDistBin),fNpid(ex.fNpid),
  72 fNmaxMixEv(ex.fNmaxMixEv),
  73 fVtxZCut(ex.fVtxZCut), fCent(ex.fCent), fRp(ex.fRp),
  74 fPi0Mass(ex.fPi0Mass),
  75 fPi0MassWindow(ex.fPi0MassWindow),
  76 fPi0OverOmegaPtCut(ex.fPi0OverOmegaPtCut),
  77 fGammaOverOmegaPtCut(ex.fGammaOverOmegaPtCut),
  78 fhEtalon(ex.fhEtalon),
  79 fRealOmega0(ex.fRealOmega0), fMixAOmega0(ex.fMixAOmega0),
  80 fMixBOmega0(ex.fMixBOmega0), fMixCOmega0(ex.fMixCOmega0),
  81 fRealOmega1(ex.fRealOmega1), fMixAOmega1(ex.fMixAOmega1),
  82 fMixBOmega1(ex.fMixBOmega1), fMixCOmega1(ex.fMixCOmega1),
  83 fRealOmega2(ex.fRealOmega2), fMixAOmega2(ex.fMixAOmega2),
  84 fMixBOmega2(ex.fMixBOmega2), fMixCOmega2(ex.fMixCOmega2),
  85 fhOmegaPriPt(ex.fhOmegaPriPt)
  86 {
  87  // cpy ctor
  88  //Do not need it
  89 }
  90 */
  91 /*
  92 //______________________________________________________________________________
  93 AliAnaOmegaToPi0Gamma & AliAnaOmegaToPi0Gamma::operator = (const AliAnaOmegaToPi0Gamma & ex)
  94 {
  95  // assignment operator
  96
  97  if(this == &ex)return *this;
  98    ((AliAnaPartCorrBaseClass *)this)->operator=(ex);
  99    fInputAODPi0 = new TClonesArray(*ex.fInputAODPi0);
 100    fInputAODGammaName = ex.fInputAODGammaName;
 101    fEventsList = ex.fEventsList;
 102
 103    fNVtxZBin=ex.fNVtxZBin;
 104    fNCentBin=ex.fNCentBin;
 105    fNRpBin=ex.fNRpBin;
 106    fNBadChDistBin=ex.fNBadChDistBin;
 107    fNpid=ex.fNpid;
 108    fNmaxMixEv =ex.fNmaxMixEv;
 109
 110    fVtxZCut=ex.fVtxZCut;
 111    fCent=ex.fCent;
 112    fRp=ex.fRp;
 113
 114    fPi0Mass=ex.fPi0Mass;
 115    fPi0MassWindow=ex.fPi0MassWindow;
 116    fPi0OverOmegaPtCut=ex.fPi0OverOmegaPtCut;
 117    fGammaOverOmegaPtCut=ex.fGammaOverOmegaPtCut;
 118
 119    fhEtalon=ex.fhEtalon;
 120    fRealOmega0=ex.fRealOmega0;
 121    fMixAOmega0=ex.fMixAOmega0;
 122    fMixBOmega0=ex.fMixBOmega0;
 123    fMixCOmega0=ex.fMixCOmega0;
 124    fRealOmega1=ex.fRealOmega1;
 125    fMixAOmega1=ex.fMixAOmega1;
 126    fMixBOmega1=ex.fMixBOmega1;
 127    fMixCOmega1=ex.fMixCOmega1;
 128    fRealOmega2=ex.fRealOmega2;
 129    fMixAOmega2=ex.fMixAOmega2;
 130    fMixBOmega2=ex.fMixBOmega2;
 131    fMixCOmega2=ex.fMixCOmega2;
 132    fhOmegaPriPt=ex.fhOmegaPriPt;
 133   return *this;
 134
 135 }
 136 */
 137
 138 //______________________________________________________________________________
 139 AliAnaOmegaToPi0Gamma::~AliAnaOmegaToPi0Gamma() {
 140
 141   //dtor
 142 //  Done by the maker
 143 //  if(fInputAODPi0){
 144 //    fInputAODPi0->Clear();
 145 //    delete fInputAODPi0;
 146 //  }
 147
 148   if(fEventsList){
 149      for(Int_t i=0;i<fNVtxZBin;i++){
 150         for(Int_t j=0;j<fNCentBin;j++){
 151            for(Int_t k=0;k<fNRpBin;k++){
 152                fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k]->Clear();
 153                delete fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k];
 154            }
 155         }
 156      }
 157   }
 158   delete [] fEventsList;
 159   fEventsList=0;
 160
 161  delete [] fVtxZCut;
 162  delete [] fCent;
 163  delete [] fRp;
 164
 165 }
 166
 167 //______________________________________________________________________________
 168 void AliAnaOmegaToPi0Gamma::InitParameters()
 169 {
 170 //Init parameters when first called the analysis
 171 //Set default parameters
 172  fInputAODGammaName="PhotonsDetector";
 173  fNVtxZBin=1;
 174  fNCentBin=1;
 175  fNRpBin=1;
 176  fNBadChDistBin=3;
 177  fNpid=1;
 178  fNmaxMixEv=8;
 179
 180  fPi0Mass=0.1348;
 181  fPi0MassWindow=0.015;
 182  fPi0OverOmegaPtCut=0.8;
 183  fGammaOverOmegaPtCut=0.2;
 184
 185  fEOverlapCluster=6;
 186 }
 187
 188
 189 //______________________________________________________________________________
 190 TList * AliAnaOmegaToPi0Gamma::GetCreateOutputObjects()
 191 {
 192   //
 193   fVtxZCut = new Double_t [fNVtxZBin];
 194   for(Int_t i=0;i<fNVtxZBin;i++) fVtxZCut[i]=10*(i+1);
 195
 196   fCent=new Double_t[fNCentBin];
 197   for(int i = 0;i<fNCentBin;i++)fCent[i]=0;
 198
 199   fRp=new Double_t[fNRpBin];
 200   for(int i = 0;i<fNRpBin;i++)fRp[i]=0;
 201   //
 202   Int_t nptbins   = GetHistoPtBins();
 203   Float_t ptmax   = GetHistoPtMax();
 204   Float_t ptmin   = GetHistoPtMin();
 205
 206   Int_t nmassbins = GetHistoMassBins();
 207   Float_t massmin = GetHistoMassMin();
 208   Float_t massmax = GetHistoMassMax();
 209
 210   fhEtalon = new TH2F("hEtalon","Histo with binning parameters", nptbins,ptmin,ptmax,nmassbins,massmin,massmax) ;
 211   fhEtalon->SetXTitle("P_{T} (GeV)") ;
 212   fhEtalon->SetYTitle("m_{inv} (GeV)") ;
 213
 214   // store them in fOutputContainer
 215   fEventsList = new TList*[fNVtxZBin*fNCentBin*fNRpBin];
 216   for(Int_t i=0;i<fNVtxZBin;i++){
 217     for(Int_t j=0;j<fNCentBin;j++){
 218       for(Int_t k=0;k<fNRpBin;k++){
 219         fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k] = new TList();
 220         fEventsList[i*fNCentBin*fNRpBin+j*fNRpBin+k]->SetOwner(kFALSE);
 221       }
 222     }
 223   }
 224
 225   TList * outputContainer = new TList() ;
 226   outputContainer->SetName(GetName());
 227   const Int_t buffersize = 255;
 228   char key[buffersize] ;
 229   char title[buffersize] ;
 230   const char * detector= fInputAODGammaName.Data();
 231   Int_t ndim=fNVtxZBin*fNCentBin*fNRpBin*fNBadChDistBin*fNpid;
 232
 233   fRealOmega0 =new TH2F*[ndim];
 234   fMixAOmega0 =new TH2F*[ndim];
 235   fMixBOmega0 =new TH2F*[ndim];
 236   fMixCOmega0 =new TH2F*[ndim];
 237
 238   fRealOmega1 =new TH2F*[ndim];
 239   fMixAOmega1 =new TH2F*[ndim];
 240   fMixBOmega1 =new TH2F*[ndim];
 241   fMixCOmega1 =new TH2F*[ndim];
 242
 243   fRealOmega2 =new TH2F*[ndim];
 244   fMixAOmega2 =new TH2F*[ndim];
 245   fMixBOmega2 =new TH2F*[ndim];
 246   fMixCOmega2 =new TH2F*[ndim];
 247
 248   fhFakeOmega = new TH2F*[fNCentBin];
 249
 250   for(Int_t i=0;i<fNVtxZBin;i++){
 251     for(Int_t j=0;j<fNCentBin;j++){
 252       for(Int_t k=0;k<fNRpBin;k++){ //at event level
 253         Int_t idim=i*fNCentBin*fNRpBin+j*fNRpBin+k;
 254         for(Int_t ipid=0;ipid<fNpid;ipid++){
 255           for(Int_t idist=0;idist<fNBadChDistBin;idist++){ //at particle level
 256
 257             Int_t index=idim*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
 258
 259             snprintf(key,buffersize,"RealToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 260             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);
 261             fRealOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
 262             fRealOmega0[index]->SetName(key) ;
 263             fRealOmega0[index]->SetTitle(title);
 264             outputContainer->Add(fRealOmega0[index]);
 265
 266             snprintf(key,buffersize,"MixAToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 267             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);
 268             fMixAOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
 269             fMixAOmega0[index]->SetName(key) ;
 270             fMixAOmega0[index]->SetTitle(title);
 271             outputContainer->Add(fMixAOmega0[index]);
 272
 273             snprintf(key,buffersize,"MixBToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 274             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);
 275             fMixBOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
 276             fMixBOmega0[index]->SetName(key) ;
 277             fMixBOmega0[index]->SetTitle(title);
 278             outputContainer->Add(fMixBOmega0[index]);
 279
 280             snprintf(key,buffersize,"MixCToPi0Gamma_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 281             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);
 282             fMixCOmega0[index]=(TH2F*)fhEtalon->Clone(key) ;
 283             fMixCOmega0[index]->SetName(key) ;
 284             fMixCOmega0[index]->SetTitle(title);
 285             outputContainer->Add(fMixCOmega0[index]);
 286
 287             snprintf(key,buffersize,"RealToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 288             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);
 289             fRealOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
 290             fRealOmega1[index]->SetName(key) ;
 291             fRealOmega1[index]->SetTitle(title);
 292             outputContainer->Add(fRealOmega1[index]);
 293
 294             snprintf(key,buffersize,"MixAToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 295             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);
 296             fMixAOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
 297             fMixAOmega1[index]->SetName(key) ;
 298             fMixAOmega1[index]->SetTitle(title);
 299             outputContainer->Add(fMixAOmega1[index]);
 300
 301             snprintf(key,buffersize,"MixBToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 302             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);
 303             fMixBOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
 304             fMixBOmega1[index]->SetName(key) ;
 305             fMixBOmega1[index]->SetTitle(title);
 306             outputContainer->Add(fMixBOmega1[index]);
 307
 308             snprintf(key,buffersize,"MixCToPi0Gamma1_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 309             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);
 310             fMixCOmega1[index]=(TH2F*)fhEtalon->Clone(key) ;
 311             fMixCOmega1[index]->SetName(key) ;
 312             fMixCOmega1[index]->SetTitle(title);
 313             outputContainer->Add(fMixCOmega1[index]);
 314
 315             snprintf(key,buffersize,"RealToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 316             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);
 317             fRealOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
 318             fRealOmega2[index]->SetName(key) ;
 319             fRealOmega2[index]->SetTitle(title);
 320             outputContainer->Add(fRealOmega2[index]);
 321
 322             snprintf(key,buffersize,"MixAToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 323             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);
 324             fMixAOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
 325             fMixAOmega2[index]->SetName(key) ;
 326             fMixAOmega2[index]->SetTitle(title);
 327             outputContainer->Add(fMixAOmega2[index]);
 328
 329             snprintf(key,buffersize,"MixBToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 330             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);
 331             fMixBOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
 332             fMixBOmega2[index]->SetName(key) ;
 333             fMixBOmega2[index]->SetTitle(title);
 334             outputContainer->Add(fMixBOmega2[index]);
 335
 336             snprintf(key,buffersize,"MixCToPi0Gamma2_Vz%dC%dRp%dPid%dDist%d",i,j,k,ipid,idist);
 337             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);
 338             fMixCOmega2[index]=(TH2F*)fhEtalon->Clone(key) ;
 339             fMixCOmega2[index]->SetName(key) ;
 340             fMixCOmega2[index]->SetTitle(title);
 341             outputContainer->Add(fMixCOmega2[index]);
 342           }
 343         }
 344       }
 345     }
 346   }
 347
 348   for(Int_t i=0;i<fNCentBin;i++){
 349      snprintf(key,buffersize, "fhFakeOmega%d",i);
 350      snprintf(title,buffersize,"FakePi0(high pt cluster)+Gamma with Centrality%d",i);
 351      fhFakeOmega[i] = (TH2F*)fhEtalon->Clone(key);
 352      fhFakeOmega[i]->SetTitle(title);
 353      outputContainer->Add(fhFakeOmega[i]);
 354   }
 355   if(IsDataMC()){
 356     snprintf(key,buffersize, "%sOmegaPri",detector);
 357     snprintf(title,buffersize,"primary #omega in %s",detector);
 358     fhOmegaPriPt=new TH1F(key, title,nptbins,ptmin,ptmax);
 359     fhOmegaPriPt->GetXaxis()->SetTitle("P_{T}");
 360     fhOmegaPriPt->GetYaxis()->SetTitle("dN/P_{T}");
 361     outputContainer->Add(fhOmegaPriPt);
 362   }
 363
 364   delete fhEtalon;
 365   return outputContainer;
 366 }
 367
 368 //______________________________________________________________________________
 369 void AliAnaOmegaToPi0Gamma::Print(const Option_t * /*opt*/) const
 370 {
 371   //Print some relevant parameters set in the analysis
 372   printf("**** Print %s %s ****\n", GetName(), GetTitle() ) ;
 373   AliAnaPartCorrBaseClass::Print(" ");
 374   printf("Omega->pi0+gamma->3gamma\n");
 375   printf("Cuts at event level:            \n");
 376   printf("Bins of vertex Z:                     %d \n", fNVtxZBin);
 377   printf("Bins of centrality:                   %d \n",fNCentBin);
 378   printf("Bins of Reaction plane:               %d\n",fNRpBin);
 379   printf("Cuts at AOD particle level:\n");
 380   printf("Number of PID:                        %d \n", fNpid);
 381   printf("Number of DistToBadChannel cuts:      %d\n", fNBadChDistBin);
 382   printf("number of events buffer to be mixed:  %d\n",fNmaxMixEv);
 383 }
 384
 385 //______________________________________________________________________________
 386 void AliAnaOmegaToPi0Gamma::MakeAnalysisFillHistograms()
 387 {
 388   //fill the MC AOD if needed first
 389   //-----------
 390   //need to be further implemented
 391   AliStack * stack = 0x0;
 392   // TParticle * primary = 0x0;
 393   TClonesArray * mcparticles0 = 0x0;
 394   //TClonesArray * mcparticles1 = 0x0;
 395   AliAODMCParticle * aodprimary = 0x0;
 396   Int_t pdg=0;
 397   Double_t pt=0;
 398   Double_t eta=0;
 399
 400   if(IsDataMC()){
 401     if(GetReader()->ReadStack()){
 402       stack =  GetMCStack() ;
 403       if(!stack){
 404         printf("AliAnaAcceptance::MakeAnalysisFillHistograms() - There is no stack!\n");
 405       }
 406       else{
 407         for(Int_t i=0 ; i<stack->GetNtrack(); i++){
 408           TParticle * prim = stack->Particle(i) ;
 409           pdg = prim->GetPdgCode() ;
 410           eta=prim->Eta();
 411           pt=prim->Pt();
 412           if(TMath::Abs(eta)<0.5) {
 413             if(pdg==223) fhOmegaPriPt->Fill(pt);
 414           }
 415         }
 416       }
 417     }
 418     else if(GetReader()->ReadAODMCParticles()){
 419       //Get the list of MC particles
 420       mcparticles0 = GetReader()->GetAODMCParticles(0);
 421       if(!mcparticles0 )     {
 422         if(GetDebug() > 0) printf("AliAnaAcceptance::MakeAnalysisFillHistograms() -  Standard MCParticles not available!\n");
 423       }
 424       //           if(GetReader()->GetSecondInputAODTree()){
 425       //               mcparticles1 = GetReader()->GetAODMCParticles(1);
 426       //               if(!mcparticles1 && GetDebug() > 0)     {
 427       //                   printf("AliAnaAcceptance::MakeAnalysisFillHistograms() -  Second input MCParticles not available!\n");
 428       //                }
 429       //           }
 430       else{
 431         for(Int_t i=0;i<mcparticles0->GetEntries();i++){
 432           aodprimary =(AliAODMCParticle*)mcparticles0->At(i);
 433           pdg = aodprimary->GetPdgCode() ;
 434           eta=aodprimary->Eta();
 435           pt=aodprimary->Pt();
 436           if(TMath::Abs(eta)<0.5) {
 437             if(pdg==223) fhOmegaPriPt->Fill(pt);
 438           }
 439
 440         }
 441       }//mcparticles0 exists
 442     }//AOD MC Particles
 443   }// is data and MC
 444
 445
 446   //process event from AOD brach
 447   //extract pi0, eta and omega analysis
 448   Int_t iRun=(GetReader()->GetInputEvent())->GetRunNumber() ;
 449   if(IsBadRun(iRun)) return ;
 450
 451   //vertex z
 452   Double_t vert[]={0,0,0} ;
 453   GetVertex(vert);
 454   Int_t curEventBin =0;
 455
 456   Int_t ivtxzbin=(Int_t)TMath::Abs(vert[2])/10;
 457   if(ivtxzbin>=fNVtxZBin)return;
 458
 459   //centrality
 460   Int_t currentCentrality = GetEventCentrality();
 461   if(currentCentrality == -1) return;
 462   Int_t optCent = GetReader()->GetCentralityOpt();
 463   Int_t icentbin=currentCentrality/(optCent/fNCentBin) ; //GetEventCentrality();
 464
 465   printf("-------------- %d  %d  %d  ",currentCentrality, optCent, icentbin);
 466   //reaction plane
 467   Int_t irpbin=0;
 468
 469   if(ivtxzbin==-1) return;
 470   curEventBin = ivtxzbin*fNCentBin*fNRpBin + icentbin*fNRpBin + irpbin;
 471   TClonesArray *aodGamma = (TClonesArray*) GetAODBranch(fInputAODGammaName); //photon array
 472   //  TClonesArray *aodGamma = (TClonesArray *) GetReader()->GetOutputEvent()->FindListObject(fInputAODGammaName); //photon array
 473   Int_t nphotons =0;
 474   if(aodGamma) nphotons= aodGamma->GetEntries();
 475   else return;
 476
 477   fInputAODPi0 = (TClonesArray*)GetInputAODBranch();  //pi0 array
 478   Int_t npi0s = 0;
 479   if(fInputAODPi0) npi0s= fInputAODPi0 ->GetEntries();
 480   else return;
 481
 482   if(nphotons<3 || npi0s<1)return; //for pi0, eta and omega->pi0+gamma->3gamma reconstruction
 483
 484   //reconstruction of omega(782)->pi0+gamma->3gamma
 485   //loop for pi0 and photon
 486   if(GetDebug() > 0) printf("omega->pi0+gamma->3gamma invariant mass analysis ! This event have %d photons and %d pi0 \n", nphotons, npi0s);
 487   for(Int_t i=0;i<npi0s;i++){
 488     AliAODPWG4Particle * pi0 = (AliAODPWG4Particle*) (fInputAODPi0->At(i)) ; //pi0
 489     TLorentzVector vpi0(pi0->Px(),pi0->Py(),pi0->Pz(),pi0->E());
 490     Int_t lab1=pi0->GetCaloLabel(0);  // photon1 from pi0 decay
 491     Int_t lab2=pi0->GetCaloLabel(1);  // photon2 from pi0 decay
 492     //for omega->pi0+gamma, it needs at least three photons per event
 493     //Get the two decay photons from pi0
 494     AliAODPWG4Particle * photon1 =0;
 495     AliAODPWG4Particle * photon2 =0;
 496     for(Int_t d1=0;d1<nphotons;d1++){
 497       for(Int_t d2=0;d2<nphotons;d2++){
 498         AliAODPWG4Particle * dp1 = (AliAODPWG4Particle*) (aodGamma->At(d1));
 499         AliAODPWG4Particle * dp2 = (AliAODPWG4Particle*) (aodGamma->At(d2));
 500         Int_t dlab1=dp1->GetCaloLabel(0);
 501         Int_t dlab2=dp2->GetCaloLabel(0);
 502         if(dlab1==lab1 && dlab2==lab2){
 503           photon1=dp1;
 504           photon2=dp2;
 505         }
 506         else continue;
 507       }
 508     }
 509     //caculate the asy and dist of the two photon from pi0 decay
 510     TLorentzVector dph1(photon1->Px(),photon1->Py(),photon1->Pz(),photon1->E());
 511     TLorentzVector dph2(photon2->Px(),photon2->Py(),photon2->Pz(),photon2->E());
 512
 513     Double_t pi0asy= TMath::Abs(dph1.E()-dph2.E())/(dph1.E()+dph2.E());
 514     //    Double_t phi1=dph1.Phi();
 515     //    Double_t phi2=dph2.Phi();
 516     //    Double_t eta1=dph1.Eta();
 517     //    Double_t eta2=dph2.Eta();
 518     //    Double_t pi0dist=TMath::Sqrt((phi1-phi2)*(phi1-phi2)+(eta1-eta2)*(eta1-eta2));
 519
 520     if(pi0->GetIdentifiedParticleType()==111  && nphotons>2 && npi0s
 521        && TMath::Abs(vpi0.M()-fPi0Mass)<fPi0MassWindow) { //pi0 candidates
 522
 523       //avoid the double counting
 524       Int_t * dc1= new Int_t[nphotons];
 525       Int_t * dc2= new Int_t[nphotons];
 526       Int_t index1=0;
 527       Int_t index2=0;
 528       for(Int_t k=0;k<i;k++){
 529         AliAODPWG4Particle * p3=(AliAODPWG4Particle*)(fInputAODPi0->At(k));
 530         Int_t lab4=p3->GetCaloLabel(0);
 531         Int_t lab5=p3->GetCaloLabel(1);
 532         if(lab1==lab4){ dc1[index1]=lab5;  index1++;  }
 533         if(lab2==lab5){ dc2[index2]=lab4;  index2++;  }
 534       }
 535
 536
 537       //loop the pi0 with third gamma
 538       for(Int_t j=0;j<nphotons;j++){
 539         AliAODPWG4Particle *photon3 = (AliAODPWG4Particle*) (aodGamma->At(j));
 540         TLorentzVector dph3(photon3->Px(),photon3->Py(),photon3->Pz(),photon3->E());
 541         Int_t lab3=photon3->GetCaloLabel(0);
 542         Double_t pi0gammapt=(vpi0+dph3).Pt();
 543         Double_t pi0gammamass=(vpi0+dph3).M();
 544         Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt;
 545         Double_t gammaOverOmegaPtRatio= dph3.Pt()/pi0gammapt;
 546
 547         //pi0, gamma pt cut
 548         if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
 549            gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
 550
 551         for(Int_t l=0;l<index1;l++) if(lab3==dc1[l]) lab3=-1;
 552         for(Int_t l=0;l<index2;l++) if(lab3==dc2[l]) lab3=-1;
 553
 554         if(lab3>0 && lab3!=lab1 && lab3!=lab2){
 555                 for(Int_t ipid=0;ipid<fNpid;ipid++){
 556             for(Int_t idist=0;idist<fNBadChDistBin;idist++){
 557               Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
 558               if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
 559                  photon2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
 560                  photon3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
 561                  photon1->DistToBad()>=idist &&
 562                  photon2->DistToBad()>=idist &&
 563                  photon3->DistToBad()>=idist ){
 564                 //fill the histograms
 565                 if(GetDebug() > 2) printf("Real: index  %d  pt  %2.3f  mass   %2.3f \n", index, pi0gammapt, pi0gammamass);
 566                 fRealOmega0[index]->Fill(pi0gammapt,pi0gammamass);
 567                 if(pi0asy<0.7) fRealOmega1[index]->Fill(pi0gammapt,pi0gammamass);
 568                 if(pi0asy<0.8) fRealOmega2[index]->Fill(pi0gammapt,pi0gammamass);
 569               }
 570             }
 571                 }
 572         }
 573       }
 574       delete []dc1;
 575       delete []dc2;
 576       if(GetDebug() > 0) printf("MixA: (r1_event1+r2_event1)+r3_event2 \n");
 577       //-------------------------
 578       //background analysis
 579       //three background
 580       // --A   (r1_event1+r2_event1)+r3_event2
 581       Int_t nMixed = fEventsList[curEventBin]->GetSize();
 582       for(Int_t im=0;im<nMixed;im++){
 583         TClonesArray* ev2= (TClonesArray*) (fEventsList[curEventBin]->At(im));
 584         for(Int_t mix1=0;mix1<ev2->GetEntries();mix1++){
 585           AliAODPWG4Particle *mix1ph = (AliAODPWG4Particle*) (ev2->At(mix1));
 586           TLorentzVector vmixph(mix1ph->Px(),mix1ph->Py(),mix1ph->Pz(),mix1ph->E());
 587           Double_t pi0gammapt=(vpi0+vmixph).Pt();
 588           Double_t pi0gammamass=(vpi0+vmixph).M();
 589           Double_t pi0OverOmegaPtRatio =vpi0.Pt()/pi0gammapt;
 590           Double_t gammaOverOmegaPtRatio= vmixph.Pt()/pi0gammapt;
 591
 592           //pi0, gamma pt cut
 593           if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
 594              gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
 595
 596                 for(Int_t ipid=0;ipid<fNpid;ipid++){
 597             for(Int_t idist=0;idist<fNBadChDistBin;idist++){
 598               Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
 599               if(photon1->IsPIDOK(ipid,AliCaloPID::kPhoton)&&
 600                  photon2->IsPIDOK(ipid,AliCaloPID::kPhoton)&&
 601                  mix1ph->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
 602                  photon1->DistToBad()>=idist &&
 603                  photon2->DistToBad()>=idist &&
 604                  mix1ph->DistToBad()>=idist ){
 605                 if(GetDebug() > 2) printf("MixA: index  %d   pt  %2.3f  mass   %2.3f \n",index, pi0gammapt, pi0gammamass);
 606                 //fill the histograms
 607                 fMixAOmega0[index]->Fill(pi0gammapt,pi0gammamass);
 608                 if(pi0asy<0.7)fMixAOmega1[index]->Fill(pi0gammapt,pi0gammamass);
 609                 if(pi0asy<0.8)fMixAOmega2[index]->Fill(pi0gammapt,pi0gammamass);
 610                 //printf("mix A  %d  %2.2f \n", index, pi0gammamass);
 611
 612               }
 613             }
 614           }
 615         }
 616       }
 617     }
 618   }
 619
 620   //
 621   // --B   (r1_event1+r2_event2)+r3_event2
 622   //
 623   if(GetDebug() >0)printf("MixB:  (r1_event1+r2_event2)+r3_event2 \n");
 624   for(Int_t i=0;i<nphotons;i++){
 625     AliAODPWG4Particle *ph1 = (AliAODPWG4Particle*) (aodGamma->At(i));
 626     TLorentzVector vph1(ph1->Px(),ph1->Py(),ph1->Pz(),ph1->E());
 627     //interrupt here...................
 628     //especially for EMCAL
 629     //we suppose the high pt clusters are overlapped pi0
 630
 631     for(Int_t j=i+1;j<nphotons;j++){
 632         AliAODPWG4Particle *ph2 = (AliAODPWG4Particle*) (aodGamma->At(j));
 633         TLorentzVector fakePi0, fakeOmega, vph;
 634
 635         if(ph1->E() > fEOverlapCluster && ph1->E() > ph2->E()) {
 636            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));
 637            vph.SetPxPyPzE(ph2->Px(),ph2->Py(),ph2->Pz(),ph2->E());
 638         }
 639         else if(ph2->E() > fEOverlapCluster && ph2->E() > ph1->E()) {
 640            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));
 641            vph.SetPxPyPzE(ph1->Px(), ph1->Py(),ph1->Pz(),ph1->E());
 642         }
 643         else continue;
 644
 645         fakeOmega=fakePi0+vph;
 646         for(Int_t ii=0;ii<fNCentBin;ii++){
 647            fhFakeOmega[icentbin]->Fill(fakeOmega.Pt(), fakeOmega.M());
 648         }
 649     }//j
 650
 651     //continue ................
 652     Int_t nMixed = fEventsList[curEventBin]->GetSize();
 653     for(Int_t ie=0;ie<nMixed;ie++){
 654       TClonesArray* ev2= (TClonesArray*) (fEventsList[curEventBin]->At(ie));
 655       for(Int_t mix1=0;mix1<ev2->GetEntries();mix1++){
 656         AliAODPWG4Particle *ph2 = (AliAODPWG4Particle*) (ev2->At(mix1));
 657         TLorentzVector vph2(ph2->Px(),ph2->Py(),ph2->Pz(),ph2->E());
 658         Double_t pi0asy = TMath::Abs(vph1.E()-vph2.E())/(vph1.E()+vph2.E());
 659         Double_t pi0mass=(vph1+vph2).M();
 660
 661         if(TMath::Abs(pi0mass-fPi0Mass)<fPi0MassWindow){//for pi0 selection
 662           for(Int_t mix2=(mix1+1);mix2<ev2->GetEntries();mix2++){
 663             AliAODPWG4Particle *ph3 = (AliAODPWG4Particle*) (ev2->At(mix2));
 664             TLorentzVector vph3(ph3->Px(),ph3->Py(),ph3->Pz(),ph3->E());
 665
 666             Double_t pi0gammapt=(vph1+vph2+vph3).Pt();
 667             Double_t pi0gammamass=(vph1+vph2+vph3).M();
 668             Double_t pi0OverOmegaPtRatio =(vph1+vph2).Pt()/pi0gammapt;
 669             Double_t gammaOverOmegaPtRatio= vph3.Pt()/pi0gammapt;
 670             //pi0, gamma pt cut
 671             if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
 672                gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
 673
 674             for(Int_t ipid=0;ipid<fNpid;ipid++){
 675               for(Int_t idist=0;idist<fNBadChDistBin;idist++){
 676                 Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
 677                 if(ph1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
 678                                ph2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
 679                    ph3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
 680                    ph1->DistToBad()>=idist &&
 681                    ph2->DistToBad()>=idist &&
 682                    ph3->DistToBad()>=idist ){
 683                   if(GetDebug() > 2) printf("MixB: index  %d   pt  %2.3f  mass   %2.3f \n", index, pi0gammapt, pi0gammamass);
 684                   //fill histograms
 685                   fMixBOmega0[index]->Fill(pi0gammapt,pi0gammamass);
 686                   if(pi0asy<0.7) fMixBOmega1[index]->Fill(pi0gammapt,pi0gammamass);
 687                   if(pi0asy<0.8) fMixBOmega2[index]->Fill(pi0gammapt,pi0gammamass);
 688                   //printf("mix B  %d  %2.2f \n", index, pi0gammamass);
 689                 }
 690               }
 691             }
 692           }
 693
 694           //
 695                 // --C   (r1_event1+r2_event2)+r3_event3
 696           //
 697           if(GetDebug() >0)printf("MixC: (r1_event1+r2_event2)+r3_event3\n");
 698           for(Int_t je=(ie+1);je<nMixed;je++){
 699             TClonesArray* ev3= (TClonesArray*) (fEventsList[curEventBin]->At(je));
 700             for(Int_t mix3=0;mix3<ev3->GetEntries();mix3++){
 701               AliAODPWG4Particle *ph3 = (AliAODPWG4Particle*) (ev3->At(mix3));
 702               TLorentzVector vph3(ph3->Px(),ph3->Py(),ph3->Pz(),ph3->E());
 703
 704               Double_t pi0gammapt=(vph1+vph2+vph3).Pt();
 705               Double_t pi0gammamass=(vph1+vph2+vph3).M();
 706               Double_t pi0OverOmegaPtRatio =(vph1+vph2).Pt()/pi0gammapt;
 707               Double_t gammaOverOmegaPtRatio= vph3.Pt()/pi0gammapt;
 708               //pi0, gamma pt cut
 709               if(pi0OverOmegaPtRatio>fPi0OverOmegaPtCut ||
 710                  gammaOverOmegaPtRatio<fGammaOverOmegaPtCut) continue;
 711
 712               for(Int_t ipid=0;ipid<fNpid;ipid++){
 713                 for(Int_t idist=0;idist<fNBadChDistBin;idist++){
 714                   Int_t index=curEventBin*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
 715                   if(ph1->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
 716                      ph2->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
 717                      ph3->IsPIDOK(ipid,AliCaloPID::kPhoton) &&
 718                      ph1->DistToBad()>=idist &&
 719                      ph2->DistToBad()>=idist &&
 720                      ph3->DistToBad()>=idist ){
 721                     if(GetDebug() > 2) printf("MixC: index  %d  pt  %2.3f  mass   %2.3f \n", index, pi0gammapt, pi0gammamass);
 722                     //fill histograms
 723                     fMixCOmega0[index]->Fill(pi0gammapt,pi0gammamass);
 724                     if(pi0asy<0.7) fMixCOmega1[index]->Fill(pi0gammapt,pi0gammamass);
 725                     if(pi0asy<0.8) fMixCOmega2[index]->Fill(pi0gammapt,pi0gammamass);
 726                     //printf("mix C  %d  %2.2f \n", index, pi0gammamass);
 727                   }
 728                 }
 729               }
 730             }
 731           }
 732         } //for pi0 selecton
 733       }
 734     }
 735   }
 736
 737
 738   //event buffer
 739   TClonesArray *currentEvent = new TClonesArray(*aodGamma);
 740   if(currentEvent->GetEntriesFast()>0){
 741     fEventsList[curEventBin]->AddFirst(currentEvent) ;
 742     currentEvent=0 ;
 743     if(fEventsList[curEventBin]->GetSize()>=fNmaxMixEv) {
 744       TClonesArray * tmp = (TClonesArray*) (fEventsList[curEventBin]->Last()) ;
 745       fEventsList[curEventBin]->RemoveLast() ;
 746       delete tmp ;
 747     }
 748   }
 749   else{
 750     delete currentEvent ;
 751     currentEvent=0 ;
 752   }
 753
 754 }
 755
 756 //______________________________________________________________________________
 757 void AliAnaOmegaToPi0Gamma::ReadHistograms(TList * outputList)
 758 {
 759  //read the histograms
 760  //for the finalization of the terminate analysis
 761
 762  Int_t index = outputList->IndexOf(outputList->FindObject(GetAddedHistogramsStringToName()+"RealToPi0Gamma_Vz0C0Rp0Pid0Dist0"));
 763
 764   Int_t ndim=fNVtxZBin*fNCentBin*fNRpBin*fNBadChDistBin*fNpid;
 765
 766  if(!fRealOmega0) fRealOmega0 =new TH2F*[ndim];
 767  if(!fMixAOmega0) fMixAOmega0 =new TH2F*[ndim];
 768  if(!fMixBOmega0) fMixBOmega0 =new TH2F*[ndim];
 769  if(!fMixCOmega0) fMixCOmega0 =new TH2F*[ndim];
 770
 771  if(!fRealOmega1) fRealOmega1 =new TH2F*[ndim];
 772  if(!fMixAOmega1) fMixAOmega1 =new TH2F*[ndim];
 773  if(!fMixBOmega1) fMixBOmega1 =new TH2F*[ndim];
 774  if(!fMixCOmega1) fMixCOmega1 =new TH2F*[ndim];
 775
 776  if(!fRealOmega2) fRealOmega2 =new TH2F*[ndim];
 777  if(!fMixAOmega2) fMixAOmega2 =new TH2F*[ndim];
 778  if(!fMixBOmega2) fMixBOmega2 =new TH2F*[ndim];
 779  if(!fMixCOmega2) fMixCOmega2 =new TH2F*[ndim];
 780
 781   for(Int_t i=0;i<fNVtxZBin;i++){
 782      for(Int_t j=0;j<fNCentBin;j++){
 783          for(Int_t k=0;k<fNRpBin;k++){ //at event level
 784              Int_t idim=i*fNCentBin*fNRpBin+j*fNRpBin+k;
 785              for(Int_t ipid=0;ipid<fNpid;ipid++){
 786                 for(Int_t idist=0;idist<fNBadChDistBin;idist++){ //at particle
 787                     Int_t ind=idim*fNpid*fNBadChDistBin+ipid*fNBadChDistBin+idist;
 788                     fRealOmega0[ind]= (TH2F*) outputList->At(index++);
 789                     fMixAOmega0[ind]= (TH2F*) outputList->At(index++);
 790                     fMixBOmega0[ind]= (TH2F*) outputList->At(index++);
 791                     fMixCOmega0[ind]= (TH2F*) outputList->At(index++);
 792
 793                     fRealOmega1[ind]= (TH2F*) outputList->At(index++);
 794                     fMixAOmega1[ind]= (TH2F*) outputList->At(index++);
 795                     fMixBOmega1[ind]= (TH2F*) outputList->At(index++);
 796                     fMixCOmega1[ind]= (TH2F*) outputList->At(index++);
 797
 798                     fRealOmega2[ind]= (TH2F*) outputList->At(index++);
 799                     fMixAOmega2[ind]= (TH2F*) outputList->At(index++);
 800                     fMixBOmega2[ind]= (TH2F*) outputList->At(index++);
 801                     fMixCOmega2[ind]= (TH2F*) outputList->At(index++);
 802
 803
 804                 }
 805               }
 806           }
 807       }
 808   }
 809
 810   if(IsDataMC()){
 811      fhOmegaPriPt  = (TH1F*)  outputList->At(index++);
 812   }
 813
 814 }
 815
 816 //______________________________________________________________________________
 817 void AliAnaOmegaToPi0Gamma::Terminate(TList * outputList)
 818 {
 819 // //Do some calculations and plots from the final histograms.
 820   if(GetDebug() >= 0) printf("AliAnaOmegaToPi0Gamma::Terminate() \n");
 821   ReadHistograms(outputList);
 822   const Int_t buffersize = 255;
 823   char cvs1[buffersize];
 824   snprintf(cvs1, buffersize, "Neutral_%s_IVM",fInputAODGammaName.Data());
 825
 826   TCanvas * cvsIVM = new TCanvas(cvs1, cvs1, 400, 10, 600, 700) ;
 827   cvsIVM->Divide(2, 2);
 828
 829   cvsIVM->cd(1);
 830   char dec[buffersize];
 831   snprintf(dec,buffersize,"h2Real_%s",fInputAODGammaName.Data());
 832   TH2F * h2Real= (TH2F*)fRealOmega0[0]->Clone(dec);
 833   h2Real->GetXaxis()->SetRangeUser(4,6);
 834   TH1F * hRealOmega = (TH1F*) h2Real->ProjectionY();
 835   hRealOmega->SetTitle("RealPi0Gamma 4<pt<6");
 836   hRealOmega->SetLineColor(2);
 837   hRealOmega->Draw();
 838
 839   cvsIVM->cd(2);
 840   snprintf(dec,buffersize,"hMixA_%s",fInputAODGammaName.Data());
 841   TH2F *h2MixA= (TH2F*)fMixAOmega0[0]->Clone(dec);
 842   h2MixA->GetXaxis()->SetRangeUser(4,6);
 843   TH1F * hMixAOmega = (TH1F*) h2MixA->ProjectionY();
 844   hMixAOmega->SetTitle("MixA 4<pt<6");
 845   hMixAOmega->SetLineColor(2);
 846   hMixAOmega->Draw();
 847
 848   cvsIVM->cd(3);
 849   snprintf(dec,buffersize,"hMixB_%s",fInputAODGammaName.Data());
 850   TH2F * h2MixB= (TH2F*)fMixBOmega0[0]->Clone(dec);
 851   h2MixB->GetXaxis()->SetRangeUser(4,6);
 852   TH1F * hMixBOmega = (TH1F*) h2MixB->ProjectionY();
 853   hMixBOmega->SetTitle("MixB 4<pt<6");
 854   hMixBOmega->SetLineColor(2);
 855   hMixBOmega->Draw();
 856
 857   cvsIVM->cd(4);
 858   snprintf(dec,buffersize,"hMixC_%s",fInputAODGammaName.Data());
 859   TH2F *h2MixC= (TH2F*)fMixCOmega0[0]->Clone(dec);
 860   h2MixC->GetXaxis()->SetRangeUser(4,6);
 861   TH1F * hMixCOmega = (TH1F*) h2MixC->ProjectionY();
 862   hMixCOmega->SetTitle("MixC 4<pt<6");
 863   hMixCOmega->SetLineColor(2);
 864   hMixCOmega->Draw();
 865
 866   char eps[buffersize];
 867   snprintf(eps,buffersize,"CVS_%s_IVM.eps",fInputAODGammaName.Data());
 868   cvsIVM->Print(eps);
 869   cvsIVM->Modified();
 870
 871 }