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