PWGGA/PHOSTasks/PHOS_PbPb/AliAnalysisTaskPi0FlowMC.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 // Extension to Pi0FLOw, mimicing AliPHOSHijingEfficiency
  17 // by Dmitri Peressounko, 05.02.2013
  18 // Authors: Henrik Qvigstad, Dmitri Peressounko
  19 // Date   : 05.04.2013
  20 /* $Id$ */
  21
  22 #include "TChain.h"
  23 #include "TTree.h"
  24 #include "TObjArray.h"
  25 #include "TF1.h"
  26 #include "TFile.h"
  27 #include "TH1F.h"
  28 #include "TH2F.h"
  29 #include "TH2I.h"
  30 #include "TH3F.h"
  31 #include "TParticle.h"
  32 #include "TCanvas.h"
  33 #include "TStyle.h"
  34 #include "TRandom.h"
  35 #include "TROOT.h"
  36 #include "THashList.h"
  37 #include "TArray.h"
  38 #include "TArrayD.h"
  39
  40
  41 #include "AliAnalysisManager.h"
  42 #include "AliMCEventHandler.h"
  43 #include "AliMCEvent.h"
  44 #include "AliStack.h"
  45 #include "AliAnalysisTaskSE.h"
  46 #include "AliPHOSHijingEfficiency.h"
  47 #include "AliCaloPhoton.h"
  48 #include "AliPHOSGeometry.h"
  49 #include "AliPHOSEsdCluster.h"
  50 #include "AliPHOSCalibData.h"
  51 #include "AliESDEvent.h"
  52 #include "AliESDCaloCells.h"
  53 #include "AliESDVertex.h"
  54 #include "AliESDtrackCuts.h"
  55 #include "AliLog.h"
  56 #include "AliPID.h"
  57 #include "AliCDBManager.h"
  58 #include "AliCentrality.h"
  59 #include "AliESDtrackCuts.h"
  60 #include "AliEventplane.h"
  61 #include "TProfile.h"
  62 #include <TPDGCode.h>
  63 #include "AliOADBContainer.h"
  64
  65
  66 #include "AliAnalysisTaskPi0Flow.h"
  67 #include "AliAnalysisTaskPi0FlowMC.h"
  68
  69 ClassImp(AliAnalysisTaskPi0FlowMC);
  70
  71 //TODO: rnlin?
  72
  73 //TODO: Geometry IHEP?
  74 //TODO: PHOS matrix?
  75 //TODO: Centrality.?
  76
  77 const Double_t AliAnalysisTaskPi0FlowMC::kRCut = 1.;
  78
  79 AliAnalysisTaskPi0FlowMC::AliAnalysisTaskPi0FlowMC(const char* name, AliAnalysisTaskPi0Flow::Period period)
  80 : AliAnalysisTaskPi0Flow(name, period),
  81   fStack(0x0)
  82 {
  83 }
  84
  85 AliAnalysisTaskPi0FlowMC::~AliAnalysisTaskPi0FlowMC()
  86 {
  87 }
  88
  89 void AliAnalysisTaskPi0FlowMC::UserCreateOutputObjects()
  90 {
  91   // Do standard Pi0Flow CreateOuputObjects
  92   AliAnalysisTaskPi0Flow::UserCreateOutputObjects();
  93
  94   // MC Generated histograms
  95   char key[55];
  96   for(Int_t cent=0; cent < fCentEdges.GetSize()-1; cent++){
  97     snprintf(key,55,"hMC_rap_gamma_cen%d",cent) ;
  98     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
  99     snprintf(key,55,"hMC_rap_pi0_cen%d",cent) ;
 100     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",200,-1.,1.)) ;
 101     snprintf(key,55,"hMC_rap_eta_cen%d",cent) ;
 102     fOutputContainer->Add(new TH1F(key,"Rapidity eta",200,-1.,1.)) ;
 103     snprintf(key,55,"hMC_phi_gamma_cen%d",cent) ;
 104     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 105     snprintf(key,55,"hMC_phi_pi0_cen%d",cent) ;
 106     fOutputContainer->Add(new TH1F(key,"Phi pi0",200,0.,TMath::TwoPi())) ;
 107     snprintf(key,55,"hMC_phi_eta_cen%d",cent) ;
 108     fOutputContainer->Add(new TH1F(key,"Phi eta",200,0.,TMath::TwoPi())) ;
 109     snprintf(key,55,"hMC_all_gamma_cen%d",cent) ;
 110     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
 111     snprintf(key,55,"hMC_all_pi0_cen%d",cent) ;
 112     fOutputContainer->Add(new TH1F(key,"Rapidity pi0",250,0.,25.)) ;
 113     snprintf(key,55,"hMC_all_eta_cen%d",cent) ;
 114     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 115     snprintf(key,55,"hMC_unitEta_gamma_cen%d",cent) ;
 116     fOutputContainer->Add(new TH1F(key,"Pt photon",250,0.,25.)) ;
 117     snprintf(key,55,"hMC_unitEta_pi0_cen%d",cent) ;
 118     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 119     snprintf(key,55,"hMC_unitEta_eta_cen%d",cent) ;
 120     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 121
 122     snprintf(key,55,"hMC_all_K0S_cen%d",cent) ;
 123     fOutputContainer->Add(new TH1F(key,"Rapidity photon",250,0.,25.)) ;
 124
 125     snprintf(key,55,"hMC_unitEta_K0S_cen%d",cent) ;
 126     fOutputContainer->Add(new TH1F(key,"Rapidity eta",250,0.,25.)) ;
 127
 128   }
 129   fOutputContainer->Add(new TH2F("hMC_gamma_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
 130   fOutputContainer->Add(new TH2F("hMC_pi0_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
 131   fOutputContainer->Add(new TH2F("hMC_eta_vertex","Creation vertex",25,0.,25.,1000,0.,500.)) ;
 132
 133
 134   Int_t nPt      = 200;
 135   Double_t ptMin = 0;
 136   Double_t ptMax = 20;
 137   fOutputContainer->Add(new TH2F("Vertex","Pi0 creation vertex",nPt,ptMin,ptMax,5000,0.,500.));
 138   fOutputContainer->Add(new TH3F("hSecondPi0RphiZ","Secondary pi0 vertex",450,0.,450.,100,0.,TMath::TwoPi(),200,-100.,100.));
 139   fOutputContainer->Add(new TH2F("hSecondPi0RE","Secondary pi0 vertex",450,0.,450.,200,0.,20.));
 140   fOutputContainer->Add(new TH3F("hMass_R","Mass vs radius any parent",50,0.,0.25,100,0.,10.,300,0.,600.));
 141   fOutputContainer->Add(new TH3F("Real_pi_R","All clusters",50,0.,0.25,100,0.,10.,250,0.,500.));
 142   fOutputContainer->Add(new TH3F("Real_pi_Z","All clusters",50,0.,0.25,100,0.,10.,100,-100.,100.));
 143 //  fOutputContainer->Add(new TH2F(Form("Real_npnp_RZ"),"All clusters",250,0.,500.,100,-100.,100.));
 144 //  fOutputContainer->Add(new TH3F(Form("Real_mass_R"),"All clusters",50,0.,0.25,100,0.,10.,300,0.,600.));
 145
 146   const Int_t nM       = 500;
 147   const Double_t mMin  = 0.0;
 148   const Double_t mMax  = 1.0;
 149
 150   for(Int_t cen=0; cen < fCentEdges.GetSize()-1; cen++){
 151     fOutputContainer->Add(new TH1F(Form("hPrimPhot_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 152     fOutputContainer->Add(new TH1F(Form("hPrimEl_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 153     fOutputContainer->Add(new TH1F(Form("hPrimPi0_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 154     fOutputContainer->Add(new TH1F(Form("hPrimEta_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 155     fOutputContainer->Add(new TH1F(Form("hPrimPipm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 156     fOutputContainer->Add(new TH1F(Form("hPrimP_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 157     fOutputContainer->Add(new TH1F(Form("hPrimPbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 158     fOutputContainer->Add(new TH1F(Form("hPrimN_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 159     fOutputContainer->Add(new TH1F(Form("hPrimNbar_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 160     fOutputContainer->Add(new TH1F(Form("hPrimK0S_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 161     fOutputContainer->Add(new TH1F(Form("hPrimK0L_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 162     fOutputContainer->Add(new TH1F(Form("hPrimKpm_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 163     fOutputContainer->Add(new TH1F(Form("hPrimOther_cen%d",cen),"Primary spetrum",nPt,ptMin,ptMax));
 164
 165     //pairs from common parents
 166     fOutputContainer->Add(new TH2F(Form("hParentAll_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 167     fOutputContainer->Add(new TH2F(Form("hParentK0s_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 168     fOutputContainer->Add(new TH2F(Form("hParentGamma_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 169     fOutputContainer->Add(new TH2F(Form("hParentEl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 170     fOutputContainer->Add(new TH2F(Form("hParentOther_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 171     fOutputContainer->Add(new TH2F(Form("hParentPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 172     fOutputContainer->Add(new TH2F(Form("hParentDirPi0_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 173
 174     //common parent - pi0
 175     fOutputContainer->Add(new TH2F(Form("hParentPi0NoPrim_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 176     fOutputContainer->Add(new TH2F(Form("hParentPi0Eta_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 177     fOutputContainer->Add(new TH2F(Form("hParentPi0Omega_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 178     fOutputContainer->Add(new TH2F(Form("hParentPi0Pipm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 179     fOutputContainer->Add(new TH2F(Form("hParentPi0Kpm_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 180     fOutputContainer->Add(new TH2F(Form("hParentPi0Ks_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 181     fOutputContainer->Add(new TH2F(Form("hParentPi0Kl_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 182     fOutputContainer->Add(new TH2F(Form("hParentPi0pn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 183     fOutputContainer->Add(new TH2F(Form("hParentPi0antipn_cen%d",cen),"All clusters",nM,mMin,mMax,nPt,ptMin,ptMax));
 184
 185   }
 186
 187
 188   //Photon contaminations
 189   fOutputContainer->Add(new TH2F("hPipmGammaConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 190   fOutputContainer->Add(new TH2F("hPipmElConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 191   fOutputContainer->Add(new TH2F("hPipmNConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 192   fOutputContainer->Add(new TH2F("hPipmOtherConvR","Conversion radius" ,200,0.,20.,1000,0.,500.));
 193   fOutputContainer->Add(new TH2F("hPipmGammaConvRZ","Conversion radius" ,400,-200.,200.,1000,0.,500.));
 194
 195    const Int_t nTypes=24 ;
 196    char partTypes[nTypes][55] ;
 197    snprintf(partTypes[0],55,"hGammaNoPrim") ; //
 198    snprintf(partTypes[1],55,"hGammaPhot") ; //
 199    snprintf(partTypes[2],55,"hGammaEl") ; //
 200    snprintf(partTypes[3],55,"hGammaPi0") ; //
 201    snprintf(partTypes[4],55,"hGammaEta") ; //
 202    snprintf(partTypes[5],55,"hhGammaOmega") ; //
 203    snprintf(partTypes[6],55,"hGammaPipm") ; //
 204    snprintf(partTypes[7],55,"hGammaP") ; //
 205    snprintf(partTypes[8],55,"hGammaPbar") ; //
 206    snprintf(partTypes[9],55,"hGammaN") ; //
 207    snprintf(partTypes[10],55,"hGammaNbar") ; //
 208    snprintf(partTypes[11],55,"hGammaK0S") ; //
 209    snprintf(partTypes[12],55,"hGammaK0L") ; //
 210    snprintf(partTypes[13],55,"hGammaKpm") ; //
 211    snprintf(partTypes[14],55,"hGammaKstar") ; //
 212    snprintf(partTypes[15],55,"hGammaDelta") ; //
 213    snprintf(partTypes[16],55,"hGammaOtherCharged") ; //
 214    snprintf(partTypes[17],55,"hGammaOtherNeutral") ; //
 215    snprintf(partTypes[18],55,"hGammaPipmGamma") ; //
 216    snprintf(partTypes[19],55,"hGammaPipmEl") ; //
 217    snprintf(partTypes[20],55,"hGammaPipmOther") ; //
 218    snprintf(partTypes[21],55,"hGammaPipmDirect") ; //
 219    snprintf(partTypes[22],55,"hGammaPipmp") ; //
 220    snprintf(partTypes[23],55,"hGammaPipmn") ; //
 221
 222    const Int_t nPID=12 ;
 223    char cPID[12][25] ;
 224    snprintf(cPID[0],25,"All") ;
 225    snprintf(cPID[1],25,"Allcore") ;
 226    snprintf(cPID[2],25,"CPV") ;
 227    snprintf(cPID[3],25,"CPVcore") ;
 228    snprintf(cPID[4],25,"CPV2") ;
 229    snprintf(cPID[5],25,"CPV2core") ;
 230    snprintf(cPID[6],25,"Disp") ;
 231    snprintf(cPID[7],25,"Dispcore") ;
 232    snprintf(cPID[8],25,"Disp2") ;
 233    snprintf(cPID[9],25,"Disp2core") ;
 234    snprintf(cPID[10],25,"Both") ;
 235    snprintf(cPID[11],25,"Bothcore") ;
 236
 237    for(Int_t itype=0; itype<nTypes; itype++){
 238      for(Int_t iPID=0; iPID<nPID; iPID++){
 239        for(Int_t cen=0; cen<5; cen++){
 240          fOutputContainer->Add(new TH1F(Form("%s_%s_cen%d",partTypes[itype],cPID[iPID],cen),"Cluster parents",nPt,ptMin,ptMax));
 241        }
 242      }
 243    }
 244
 245   PostData(1, fOutputContainer);
 246 }
 247
 248 void AliAnalysisTaskPi0FlowMC::UserExec(Option_t* option)
 249 {
 250   fStack = GetMCStack();
 251
 252   AliAnalysisTaskPi0Flow::UserExec(option);
 253 }
 254
 255
 256 void AliAnalysisTaskPi0FlowMC::SelectPhotonClusters()
 257 {
 258   AliAnalysisTaskPi0Flow::SelectPhotonClusters();
 259
 260   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
 261     AliCaloPhoton * photon = (AliCaloPhoton*)fCaloPhotonsPHOS->At(i1);
 262     AliVCluster* cluster = photon->GetCluster();
 263     Bool_t sure=  kTRUE;
 264     Int_t primary=FindPrimary(cluster,sure) ;
 265     photon->SetPrimary(primary);
 266     photon->SetWeight(PrimaryWeight(primary)) ;
 267
 268     TParticle* p = fStack->Particle(primary);
 269
 270     if(p->R() >kRCut) {
 271       if(p->GetPdgCode()==11 || p->GetPdgCode()==-11) continue;
 272       else { fCaloPhotonsPHOS->Remove(photon); fCaloPhotonsPHOS->Compress(); }
 273     }
 274   }
 275 }
 276
 277 void AliAnalysisTaskPi0FlowMC::FillSelectedClusterHistograms()
 278 {
 279   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
 280     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
 281
 282     Double_t dphiA=ph1->Phi()-fRPV0A ;
 283     while(dphiA<0)dphiA+=TMath::Pi() ;
 284     while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
 285
 286     Double_t dphiC=ph1->Phi()-fRPV0C ;
 287     while(dphiC<0)dphiC+=TMath::Pi() ;
 288     while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
 289
 290     Double_t dphiT=ph1->Phi()-fRP ;
 291     while(dphiT<0)dphiT+=TMath::Pi() ;
 292     while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
 293
 294     Double_t pt = ph1->Pt() ;
 295     Double_t ptcore = ph1->GetMomV2()->Pt() ;
 296     Double_t w = ph1->GetWeight();
 297
 298     FillHistogram(Form("hPhotPhiV0AAll_cen%d",fCentBin),pt,dphiA, w) ;
 299     FillHistogram(Form("hPhotPhiV0CAll_cen%d",fCentBin),pt,dphiC, w) ;
 300     if(fHaveTPCRP)
 301       FillHistogram(Form("hPhotPhiTPCAll_cen%d",fCentBin),pt,dphiT, w) ;
 302     FillHistogram(Form("hPhotPhiV0AAllcore_cen%d",fCentBin),ptcore,dphiA, w) ;
 303     FillHistogram(Form("hPhotPhiV0CAllcore_cen%d",fCentBin),ptcore,dphiC, w) ;
 304     if(fHaveTPCRP)
 305       FillHistogram(Form("hPhotPhiTPCAllcore_cen%d",fCentBin),ptcore,dphiT, w) ;
 306
 307     FillHistogram(Form("hPhotAll_cen%d",fCentBin),pt, w) ;
 308     FillHistogram(Form("hPhotAllcore_cen%d",fCentBin),ptcore, w) ;
 309     if(ph1->IsntUnfolded()){
 310       FillHistogram(Form("hPhotAllwou_cen%d",fCentBin),pt, w) ;
 311       FillHistogram(Form("hPhotPhiV0AAllwou_cen%d",fCentBin),pt,dphiA, w) ;
 312       FillHistogram(Form("hPhotPhiV0CAllwou_cen%d",fCentBin),pt,dphiC, w) ;
 313       if(fHaveTPCRP)
 314         FillHistogram(Form("hPhotPhiTPCAllwou_cen%d",fCentBin),pt,dphiT, w) ;
 315     }
 316     if(ph1->IsCPVOK()){
 317       FillHistogram(Form("hPhotPhiV0ACPV_cen%d",fCentBin),pt,dphiA, w) ;
 318       FillHistogram(Form("hPhotPhiV0CCPV_cen%d",fCentBin),pt,dphiC, w) ;
 319       if(fHaveTPCRP)
 320         FillHistogram(Form("hPhotPhiTPCCPV_cen%d",fCentBin),pt,dphiT, w) ;
 321
 322       FillHistogram(Form("hPhotPhiV0ACPVcore_cen%d",fCentBin),ptcore,dphiA, w) ;
 323       FillHistogram(Form("hPhotPhiV0CCPVcore_cen%d",fCentBin),ptcore,dphiC, w) ;
 324       if(fHaveTPCRP)
 325         FillHistogram(Form("hPhotPhiTPCCPVcore_cen%d",fCentBin),ptcore,dphiT, w) ;
 326
 327       FillHistogram(Form("hPhotCPV_cen%d",fCentBin),pt, w) ;
 328       FillHistogram(Form("hPhotCPVcore_cen%d",fCentBin),ptcore, w) ;
 329     }
 330     if(ph1->IsCPV2OK()){
 331       FillHistogram(Form("hPhotPhiV0ACPV2_cen%d",fCentBin),pt,dphiA, w) ;
 332       FillHistogram(Form("hPhotPhiV0CCPV2_cen%d",fCentBin),pt,dphiC, w) ;
 333       if(fHaveTPCRP)
 334         FillHistogram(Form("hPhotPhiTPCCPV2_cen%d",fCentBin),pt,dphiT, w) ;
 335
 336       FillHistogram(Form("hPhotPhiV0ACPV2core_cen%d",fCentBin),ptcore,dphiA, w) ;
 337       FillHistogram(Form("hPhotPhiV0CCPV2core_cen%d",fCentBin),ptcore,dphiC, w) ;
 338       if(fHaveTPCRP)
 339         FillHistogram(Form("hPhotPhiTPCCPV2core_cen%d",fCentBin),ptcore,dphiT, w) ;
 340       FillHistogram(Form("hPhotCPV2_cen%d",fCentBin),pt, w) ;
 341       FillHistogram(Form("hPhotCPV2core_cen%d",fCentBin),ptcore, w) ;
 342     }
 343     if(ph1->IsDispOK()){
 344       FillHistogram(Form("hPhotPhiV0ADisp_cen%d",fCentBin),pt,dphiA, w) ;
 345       FillHistogram(Form("hPhotPhiV0CDisp_cen%d",fCentBin),pt,dphiC, w) ;
 346       if(fHaveTPCRP)
 347         FillHistogram(Form("hPhotPhiTPCDisp_cen%d",fCentBin),pt,dphiT, w) ;
 348
 349       FillHistogram(Form("hPhotPhiV0ADispcore_cen%d",fCentBin),ptcore,dphiA, w) ;
 350       FillHistogram(Form("hPhotPhiV0CDispcore_cen%d",fCentBin),ptcore,dphiC, w) ;
 351       if(fHaveTPCRP)
 352         FillHistogram(Form("hPhotPhiTPCDispcore_cen%d",fCentBin),ptcore,dphiT, w) ;
 353
 354       if(ph1->IsntUnfolded()){
 355         FillHistogram(Form("hPhotPhiV0ADispwou_cen%d",fCentBin),pt,dphiA, w) ;
 356         FillHistogram(Form("hPhotPhiV0CDispwou_cen%d",fCentBin),pt,dphiC, w) ;
 357         if(fHaveTPCRP)
 358           FillHistogram(Form("hPhotPhiTPCDispwou_cen%d",fCentBin),pt,dphiT, w) ;
 359
 360       }
 361       FillHistogram(Form("hPhotDisp_cen%d",fCentBin),pt, w) ;
 362       FillHistogram(Form("hPhotDispcore_cen%d",fCentBin),ptcore, w) ;
 363       if(ph1->IsntUnfolded()){
 364         FillHistogram(Form("hPhotDispwou_cen%d",fCentBin),pt, w) ;
 365       }
 366       if(ph1->IsCPVOK()){
 367         FillHistogram(Form("hPhotPhiV0ABoth_cen%d",fCentBin),pt,dphiA, w) ;
 368         FillHistogram(Form("hPhotPhiV0CBoth_cen%d",fCentBin),pt,dphiC, w) ;
 369         if(fHaveTPCRP)
 370           FillHistogram(Form("hPhotPhiTPCBoth_cen%d",fCentBin),pt,dphiT, w) ;
 371
 372         FillHistogram(Form("hPhotPhiV0ABothcore_cen%d",fCentBin),ptcore,dphiA, w) ;
 373         FillHistogram(Form("hPhotPhiV0CBothcore_cen%d",fCentBin),ptcore,dphiC, w) ;
 374         if(fHaveTPCRP)
 375           FillHistogram(Form("hPhotPhiTPCBothcore_cen%d",fCentBin),ptcore,dphiT, w) ;
 376
 377         FillHistogram(Form("hPhotBoth_cen%d",fCentBin),pt, w) ;
 378         FillHistogram(Form("hPhotBothcore_cen%d",fCentBin),ptcore, w) ;
 379       }
 380     }
 381     if(ph1->IsDisp2OK()){
 382       FillHistogram(Form("hPhotPhiV0ADisp2_cen%d",fCentBin),pt,dphiA, w) ;
 383       FillHistogram(Form("hPhotPhiV0CDisp2_cen%d",fCentBin),pt,dphiC, w) ;
 384       if(fHaveTPCRP)
 385         FillHistogram(Form("hPhotPhiTPCDisp2_cen%d",fCentBin),pt,dphiT, w) ;
 386       FillHistogram(Form("hPhotPhiV0ADisp2core_cen%d",fCentBin),ptcore,dphiA, w) ;
 387       FillHistogram(Form("hPhotPhiV0CDisp2core_cen%d",fCentBin),ptcore,dphiC, w) ;
 388       if(fHaveTPCRP)
 389         FillHistogram(Form("hPhotPhiTPCDisp2core_cen%d",fCentBin),ptcore,dphiT, w) ;
 390
 391       FillHistogram(Form("hPhotDisp2_cen%d",fCentBin),pt, w) ;
 392       FillHistogram(Form("hPhotDisp2core_cen%d",fCentBin),ptcore, w) ;
 393       if(ph1->IsCPVOK()){
 394         FillHistogram(Form("hPhotPhiV0ABoth2_cen%d",fCentBin),pt,dphiA, w) ;
 395         FillHistogram(Form("hPhotPhiV0CBoth2_cen%d",fCentBin),pt,dphiC, w) ;
 396         if(fHaveTPCRP)
 397           FillHistogram(Form("hPhotPhiTPCBoth2_cen%d",fCentBin),pt,dphiT, w) ;
 398
 399         FillHistogram(Form("hPhotPhiV0ABoth2core_cen%d",fCentBin),ptcore,dphiA, w) ;
 400         FillHistogram(Form("hPhotPhiV0CBoth2core_cen%d",fCentBin),ptcore,dphiC, w) ;
 401         if(fHaveTPCRP)
 402           FillHistogram(Form("hPhotPhiTPCBoth2core_cen%d",fCentBin),ptcore,dphiT, w) ;
 403
 404         FillHistogram(Form("hPhotBoth2_cen%d",fCentBin),pt, w) ;
 405         FillHistogram(Form("hPhotBoth2core_cen%d",fCentBin),ptcore, w) ;
 406       }
 407     }
 408   }
 409 }
 410
 411 void AliAnalysisTaskPi0FlowMC::ConsiderPi0s()
 412 {
 413   char key[55];
 414   for (Int_t i1=0; i1 < fCaloPhotonsPHOS->GetEntriesFast()-1; i1++) {
 415     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
 416     const Double_t w1 = ph1->GetWeight();
 417     for (Int_t i2=i1+1; i2<fCaloPhotonsPHOS->GetEntriesFast(); i2++) {
 418       AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
 419       TLorentzVector p12  = *ph1  + *ph2;
 420       TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 421
 422       const Double_t w2 = ph2->GetWeight();
 423       Double_t w = TMath::Sqrt(w1*w2);
 424
 425       FillHistogram("hPHOSphi",fCentrality,p12.Pt(),p12.Phi(), w) ;
 426       Double_t dphiA=p12.Phi()-fRPV0A ;
 427       while(dphiA<0)dphiA+=TMath::Pi() ;
 428       while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
 429
 430       Double_t dphiC=p12.Phi()-fRPV0C ;
 431       while(dphiC<0)dphiC+=TMath::Pi() ;
 432       while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
 433
 434       Double_t dphiT=p12.Phi()-fRP ;
 435       while(dphiT<0)dphiT+=TMath::Pi() ;
 436       while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
 437
 438       Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 439       Double_t m=p12.M() ;
 440       Double_t mcore=pv12.M() ;
 441       Double_t pt=p12.Pt() ;
 442       Double_t ptcore=pv12.Pt() ;
 443       Double_t pt1=ph1->Pt() ;
 444       Double_t pt2=ph2->Pt() ;
 445       Double_t ptcore1=ph1->GetMomV2()->Pt() ;
 446       Double_t ptcore2=ph2->GetMomV2()->Pt() ;
 447
 448       FillHistogram(Form("hMassPtV0AAll_cen%d",fCentBin),m,pt,dphiA, w) ;
 449       FillHistogram(Form("hMassPtV0CAll_cen%d",fCentBin),m,pt,dphiC, w) ;
 450       if(fHaveTPCRP)
 451         FillHistogram(Form("hMassPtTPCAll_cen%d",fCentBin),m,pt,dphiT, w) ;
 452
 453       FillHistogram(Form("hMassPtV0AAllcore_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
 454       FillHistogram(Form("hMassPtV0CAllcore_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
 455       if(fHaveTPCRP)
 456         FillHistogram(Form("hMassPtTPCAllcore_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
 457
 458
 459       FillHistogram(Form("hPi0All_cen%d",fCentBin),m,pt, w) ;
 460       FillHistogram(Form("hPi0Allcore_cen%d",fCentBin),mcore,ptcore, w) ;
 461       if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
 462         FillHistogram(Form("hPi0Allwou_cen%d",fCentBin),m,pt, w) ;
 463         FillHistogram(Form("hMassPtV0AAllwou_cen%d",fCentBin),m,pt,dphiA, w) ;
 464         FillHistogram(Form("hMassPtV0CAllwou_cen%d",fCentBin),m,pt,dphiC, w) ;
 465         if(fHaveTPCRP)
 466           FillHistogram(Form("hMassPtTPCAllwou_cen%d",fCentBin),m,pt,dphiT, w) ;
 467       }
 468
 469       FillHistogram(Form("hSingleAll_cen%d",fCentBin),m,pt1, w) ;
 470       FillHistogram(Form("hSingleAll_cen%d",fCentBin),m,pt2, w) ;
 471       FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),mcore,ptcore1, w) ;
 472       FillHistogram(Form("hSingleAllcore_cen%d",fCentBin),mcore,ptcore2, w) ;
 473       if(ph1->IsntUnfolded())
 474         FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),m,pt1, w) ;
 475       if(ph2->IsntUnfolded())
 476         FillHistogram(Form("hSingleAllwou_cen%d",fCentBin),m,pt2, w) ;
 477       if(ph1->IsCPVOK()){
 478         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),m,pt1, w) ;
 479         FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),mcore,ptcore1, w) ;
 480       }
 481       if(ph2->IsCPVOK()){
 482         FillHistogram(Form("hSingleCPV_cen%d",fCentBin),m,pt2, w) ;
 483         FillHistogram(Form("hSingleCPVcore_cen%d",fCentBin),mcore,ptcore2, w) ;
 484       }
 485       if(ph1->IsCPV2OK()){
 486         FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),m,pt1, w) ;
 487         FillHistogram(Form("hSingleCPV2core_cen%d",fCentBin),mcore,ptcore2, w) ;
 488       }
 489       if(ph2->IsCPV2OK()){
 490         FillHistogram(Form("hSingleCPV2_cen%d",fCentBin),m,pt2, w) ;
 491         FillHistogram(Form("hSingleCPV2core_cen%d",fCentBin),mcore,ptcore2, w) ;
 492       }
 493       if(ph1->IsDispOK()){
 494         FillHistogram(Form("hSingleDisp_cen%d",fCentBin),m,pt1, w) ;
 495         if(ph1->IsntUnfolded()){
 496           FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),m,pt1, w) ;
 497         }
 498         FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),mcore,ptcore1, w) ;
 499       }
 500       if(ph2->IsDispOK()){
 501         FillHistogram(Form("hSingleDisp_cen%d",fCentBin),m,pt2, w) ;
 502         if(ph1->IsntUnfolded()){
 503           FillHistogram(Form("hSingleDispwou_cen%d",fCentBin),m,pt2, w) ;
 504         }
 505         FillHistogram(Form("hSingleDispcore_cen%d",fCentBin),mcore,ptcore2, w) ;
 506       }
 507       if(ph1->IsDisp2OK()){
 508         FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),m,pt1, w) ;
 509         FillHistogram(Form("hSingleDisp2core_cen%d",fCentBin),mcore,ptcore1, w) ;
 510       }
 511       if(ph2->IsDisp2OK()){
 512         FillHistogram(Form("hSingleDisp2_cen%d",fCentBin),m,pt2, w) ;
 513         FillHistogram(Form("hSingleDisp2core_cen%d",fCentBin),mcore,ptcore1, w) ;
 514       }
 515       if(ph1->IsDispOK() && ph1->IsCPVOK()){
 516         FillHistogram(Form("hSingleBoth_cen%d",fCentBin),m,pt1, w) ;
 517         FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),mcore,ptcore1, w) ;
 518       }
 519       if(ph2->IsDispOK() && ph2->IsCPVOK()){
 520         FillHistogram(Form("hSingleBoth_cen%d",fCentBin),m,pt2, w) ;
 521         FillHistogram(Form("hSingleBothcore_cen%d",fCentBin),mcore,ptcore2, w) ;
 522       }
 523       if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
 524         FillHistogram(Form("hSingleBoth2_cen%d",fCentBin),m,pt1, w) ;
 525         FillHistogram(Form("hSingleBoth2core_cen%d",fCentBin),mcore,ptcore1, w) ;
 526       }
 527       if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
 528         FillHistogram(Form("hSingleBoth2_cen%d",fCentBin),m,pt2, w) ;
 529         FillHistogram(Form("hSingleBoth2core_cen%d",fCentBin),mcore,ptcore2, w) ;
 530       }
 531
 532
 533       if(a<kAlphaCut){
 534         FillHistogram(Form("hPi0All_a07_cen%d",fCentBin),m,pt, w) ;
 535       }
 536
 537       if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 538         snprintf(key,55,"hMassPtCPV_cen%d",fCentBin) ;
 539         FillHistogram(Form("hMassPtV0ACPV_cen%d",fCentBin),m,pt,dphiA, w) ;
 540         FillHistogram(Form("hMassPtV0CCPV_cen%d",fCentBin),m,pt,dphiC, w) ;
 541         if(fHaveTPCRP)
 542           FillHistogram(Form("hMassPtTPCCPV_cen%d",fCentBin),m,pt,dphiT, w) ;
 543
 544         FillHistogram(Form("hMassPtV0ACPVcore_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
 545         FillHistogram(Form("hMassPtV0CCPVcore_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
 546         if(fHaveTPCRP)
 547           FillHistogram(Form("hMassPtTPCCPVcore_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
 548
 549         FillHistogram(Form("hPi0CPV_cen%d",fCentBin),m,pt, w) ;
 550         FillHistogram(Form("hPi0CPVcore_cen%d",fCentBin),mcore, ptcore, w) ;
 551
 552         if(a<kAlphaCut){
 553           FillHistogram(Form("hPi0CPV_a07_cen%d",fCentBin),m,pt, w) ;
 554         }
 555       }
 556       if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
 557         FillHistogram(Form("hMassPtV0ACPV2_cen%d",fCentBin),m,pt,dphiA, w) ;
 558         FillHistogram(Form("hMassPtV0CCPV2_cen%d",fCentBin),m,pt,dphiC, w) ;
 559         if(fHaveTPCRP)
 560           FillHistogram(Form("hMassPtTPCCPV2_cen%d",fCentBin),m,pt,dphiT, w) ;
 561         FillHistogram(Form("hMassPtV0ACPV2core_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
 562         FillHistogram(Form("hMassPtV0CCPV2core_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
 563         if(fHaveTPCRP)
 564           FillHistogram(Form("hMassPtTPCCPV2core_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
 565
 566         FillHistogram(Form("hPi0CPV2_cen%d",fCentBin),m,pt, w) ;
 567         FillHistogram(Form("hPi0CPV2core_cen%d",fCentBin),mcore, ptcore, w) ;
 568         if(a<kAlphaCut){
 569           FillHistogram(Form("hPi0CPV2_a07_cen%d",fCentBin),m,pt, w) ;
 570         }
 571       }
 572       if(ph1->IsDispOK() && ph2->IsDispOK()){
 573         snprintf(key,55,"hMassPtDisp_cen%d",fCentBin) ;
 574         FillHistogram(Form("hMassPtV0ADisp_cen%d",fCentBin),m,pt,dphiA, w) ;
 575         FillHistogram(Form("hMassPtV0CDisp_cen%d",fCentBin),m,pt,dphiC, w) ;
 576         if(fHaveTPCRP)
 577           FillHistogram(Form("hMassPtTPCDisp_cen%d",fCentBin),m,pt,dphiT, w) ;
 578
 579         FillHistogram(Form("hMassPtV0ADispcore_cen%d",fCentBin),mcore, ptcore,dphiA, w) ;
 580         FillHistogram(Form("hMassPtV0CDispcore_cen%d",fCentBin),mcore, ptcore,dphiC, w) ;
 581         if(fHaveTPCRP)
 582           FillHistogram(Form("hMassPtTPCDispcore_cen%d",fCentBin),mcore, ptcore,dphiT, w) ;
 583
 584         FillHistogram(Form("hPi0Disp_cen%d",fCentBin),m,pt, w) ;
 585         FillHistogram(Form("hPi0Dispcore_cen%d",fCentBin),mcore, ptcore, w) ;
 586
 587         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
 588           FillHistogram(Form("hPi0Dispwou_cen%d",fCentBin),m,pt, w) ;
 589
 590           FillHistogram(Form("hMassPtV0ADispwou_cen%d",fCentBin),m,pt,dphiA, w) ;
 591           FillHistogram(Form("hMassPtV0CDispwou_cen%d",fCentBin),m,pt,dphiC, w) ;
 592           if(fHaveTPCRP)
 593             FillHistogram(Form("hMassPtTPCDispwou_cen%d",fCentBin),m,pt,dphiT, w) ;
 594         }
 595
 596         if(a<kAlphaCut){
 597           FillHistogram(Form("hPi0Disp_a07_cen%d",fCentBin),m,pt, w) ;
 598         }
 599         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 600           FillHistogram(Form("hMassPtV0ABoth_cen%d",fCentBin),m,pt,dphiA, w) ;
 601           FillHistogram(Form("hMassPtV0CBoth_cen%d",fCentBin),m,pt,dphiC, w) ;
 602           if(fHaveTPCRP)
 603             FillHistogram(Form("hMassPtTPCBoth_cen%d",fCentBin),m,pt,dphiT, w) ;
 604
 605           FillHistogram(Form("hMassPtV0ABothcore_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
 606           FillHistogram(Form("hMassPtV0CBothcore_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
 607           if(fHaveTPCRP)
 608             FillHistogram(Form("hMassPtTPCBothcore_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
 609
 610           FillHistogram(Form("hPi0Both_cen%d",fCentBin),m,pt, w) ;
 611           FillHistogram(Form("hPi0Bothcore_cen%d",fCentBin),mcore,ptcore, w) ;
 612
 613           if(a<kAlphaCut){
 614             snprintf(key,55,"hPi0Both_a07_cen%d",fCentBin) ;
 615             FillHistogram(Form("hPi0Both_a07_cen%d",fCentBin),m,pt, w) ;
 616           }
 617           if(ph1->Module()==1 && ph2->Module()==1)
 618             FillHistogram("hPi0M11", m, pt, w) ;
 619           else if(ph1->Module()==2 && ph2->Module()==2)
 620             FillHistogram("hPi0M22", m, pt, w) ;
 621           else if(ph1->Module()==3 && ph2->Module()==3)
 622             FillHistogram("hPi0M33", m, pt, w) ;
 623           else if(ph1->Module()==1 && ph2->Module()==2)
 624             FillHistogram("hPi0M12", m, pt, w) ;
 625           else if(ph1->Module()==1 && ph2->Module()==3)
 626             FillHistogram("hPi0M13", m, pt, w) ;
 627           else if(ph1->Module()==2 && ph2->Module()==3)
 628             FillHistogram("hPi0M23", m, pt, w) ;
 629
 630         }
 631
 632       }
 633
 634
 635       if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
 636         FillHistogram(Form("hPi0Disp2_cen%d",fCentBin),m,pt,w) ;
 637         FillHistogram(Form("hPi0Disp2core_cen%d",fCentBin),mcore, ptcore, w) ;
 638
 639         FillHistogram(Form("hMassPtV0ADisp2_cen%d",fCentBin),m,pt,dphiA, w) ;
 640         FillHistogram(Form("hMassPtV0CDisp2_cen%d",fCentBin),m,pt,dphiC, w) ;
 641         if(fHaveTPCRP)
 642           FillHistogram(Form("hMassPtTPCDisp2_cen%d",fCentBin),m,pt,dphiT, w) ;
 643
 644         FillHistogram(Form("hMassPtV0ADisp2core_cen%d",fCentBin),mcore, ptcore,dphiA, w) ;
 645         FillHistogram(Form("hMassPtV0CDisp2core_cen%d",fCentBin),mcore, ptcore,dphiC, w) ;
 646         if(fHaveTPCRP)
 647           FillHistogram(Form("hMassPtTPCDisp2core_cen%d",fCentBin),mcore, ptcore,dphiT, w) ;
 648
 649         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 650           FillHistogram(Form("hMassPtV0ABoth2_cen%d",fCentBin),m,pt,dphiA, w) ;
 651           FillHistogram(Form("hMassPtV0CBoth2_cen%d",fCentBin),m,pt,dphiC, w) ;
 652           if(fHaveTPCRP)
 653             FillHistogram(Form("hMassPtTPCBoth2_cen%d",fCentBin),m,pt,dphiT, w) ;
 654
 655           FillHistogram(Form("hMassPtV0ABoth2core_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
 656           FillHistogram(Form("hMassPtV0CBoth2core_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
 657           if(fHaveTPCRP)
 658             FillHistogram(Form("hMassPtTPCBoth2core_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
 659
 660           FillHistogram(Form("hPi0Both2_cen%d",fCentBin),m,pt, w) ;
 661           FillHistogram(Form("hPi0Both2core_cen%d",fCentBin),mcore,ptcore, w) ;
 662         }
 663
 664       }
 665     } // end of loop i2
 666   } // end of loop i1
 667 }
 668
 669 //________________________________________________________________________
 670 void AliAnalysisTaskPi0FlowMC::ConsiderPi0sMix()
 671 {
 672   char key[55];
 673
 674   TList * arrayList = GetCaloPhotonsPHOSList(fVtxBin, fCentBin, fEMRPBin);
 675
 676   for (Int_t i1=0; i1<fCaloPhotonsPHOS->GetEntriesFast(); i1++) {
 677     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
 678     const Double_t w1 = ph1->GetWeight();
 679     for(Int_t evi=0; evi<arrayList->GetEntries();evi++){
 680       TObjArray * mixPHOS = static_cast<TObjArray*>(arrayList->At(evi));
 681       for(Int_t i2=0; i2<mixPHOS->GetEntriesFast();i2++){
 682         AliCaloPhoton * ph2=(AliCaloPhoton*)mixPHOS->At(i2) ;
 683         TLorentzVector p12  = *ph1  + *ph2;
 684         TLorentzVector pv12 = *(ph1->GetMomV2()) + *(ph2->GetMomV2());
 685
 686         const Double_t w2 = ph2->GetWeight();
 687         Double_t w = TMath::Sqrt(w1*w2);
 688
 689         Double_t dphiA=p12.Phi()-fRPV0A ;
 690         while(dphiA<0)dphiA+=TMath::Pi() ;
 691         while(dphiA>TMath::Pi())dphiA-=TMath::Pi() ;
 692
 693         Double_t dphiC=p12.Phi()-fRPV0C ;
 694         while(dphiC<0)dphiC+=TMath::Pi() ;
 695         while(dphiC>TMath::Pi())dphiC-=TMath::Pi() ;
 696
 697         Double_t dphiT=p12.Phi()-fRP ;
 698         while(dphiT<0)dphiT+=TMath::Pi() ;
 699         while(dphiT>TMath::Pi())dphiT-=TMath::Pi() ;
 700
 701
 702         Double_t a=TMath::Abs((ph1->E()-ph2->E())/(ph1->E()+ph2->E())) ;
 703         Double_t m=p12.M() ;
 704         Double_t mcore=pv12.M() ;
 705         Double_t pt=p12.Pt() ;
 706         Double_t ptcore=pv12.Pt() ;
 707         Double_t pt1=ph1->Pt() ;
 708         Double_t pt2=ph2->Pt() ;
 709         Double_t ptcore1=ph1->GetMomV2()->Pt() ;
 710         Double_t ptcore2=ph2->GetMomV2()->Pt() ;
 711
 712
 713         snprintf(key,55,"hMiMassPtAll_cen%d",fCentBin) ;
 714         FillHistogram(Form("hMiMassPtV0AAll_cen%d",fCentBin),m,pt,dphiA, w) ;
 715         FillHistogram(Form("hMiMassPtV0CAll_cen%d",fCentBin),m,pt,dphiC, w) ;
 716         if(fHaveTPCRP)
 717           FillHistogram(Form("hMiMassPtTPCAll_cen%d",fCentBin),m,pt,dphiT, w) ;
 718
 719         FillHistogram(Form("hMiMassPtV0AAllcore_cen%d",fCentBin),mcore, ptcore, dphiA, w) ;
 720         FillHistogram(Form("hMiMassPtV0CAllcore_cen%d",fCentBin),mcore, ptcore, dphiC, w) ;
 721         if(fHaveTPCRP)
 722           FillHistogram(Form("hMiMassPtTPCAllcore_cen%d",fCentBin),mcore, ptcore, dphiT, w) ;
 723
 724         FillHistogram(Form("hMiPi0All_cen%d",fCentBin),m,pt, w) ;
 725         FillHistogram(Form("hMiPi0Allcore_cen%d",fCentBin),mcore,ptcore, w) ;
 726         if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
 727           FillHistogram(Form("hMiPi0Allwou_cen%d",fCentBin),m,pt, w) ;
 728           FillHistogram(Form("hMiMassPtV0AAllwou_cen%d",fCentBin),m,pt,dphiA, w) ;
 729           FillHistogram(Form("hMiMassPtV0CAllwou_cen%d",fCentBin),m,pt,dphiC, w) ;
 730           if(fHaveTPCRP)
 731             FillHistogram(Form("hMiMassPtTPCAllwou_cen%d",fCentBin),m,pt,dphiT, w) ;
 732         }
 733
 734         FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),m,pt1, w) ;
 735         FillHistogram(Form("hMiSingleAll_cen%d",fCentBin),m,pt2, w) ;
 736         FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),mcore,ptcore1, w) ;
 737         FillHistogram(Form("hMiSingleAllcore_cen%d",fCentBin),mcore,ptcore2, w) ;
 738         if(ph1->IsntUnfolded())
 739           FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),m,pt1, w) ;
 740         if(ph2->IsntUnfolded())
 741           FillHistogram(Form("hMiSingleAllwou_cen%d",fCentBin),m,pt2, w) ;
 742         if(ph1->IsCPVOK()){
 743           FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),m,pt1, w) ;
 744           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),mcore,ptcore1, w) ;
 745         }
 746         if(ph2->IsCPVOK()){
 747           FillHistogram(Form("hMiSingleCPV_cen%d",fCentBin),m,pt2, w) ;
 748           FillHistogram(Form("hMiSingleCPVcore_cen%d",fCentBin),mcore,ptcore2, w) ;
 749         }
 750         if(ph1->IsCPV2OK()){
 751           FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),m,pt1, w) ;
 752           FillHistogram(Form("hMiSingleCPV2core_cen%d",fCentBin),mcore,ptcore1, w) ;
 753         }
 754         if(ph2->IsCPV2OK()){
 755           FillHistogram(Form("hMiSingleCPV2_cen%d",fCentBin),m,pt2, w) ;
 756           FillHistogram(Form("hMiSingleCPV2core_cen%d",fCentBin),mcore,ptcore2, w) ;
 757         }
 758         if(ph1->IsDispOK()){
 759           FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),m,pt1, w) ;
 760           if(ph1->IsntUnfolded()){
 761             FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),m,pt1, w) ;
 762           }
 763           FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),mcore,ptcore1, w) ;
 764         }
 765         if(ph2->IsDispOK()){
 766           FillHistogram(Form("hMiSingleDisp_cen%d",fCentBin),m,pt2, w) ;
 767           if(ph1->IsntUnfolded()){
 768             FillHistogram(Form("hMiSingleDispwou_cen%d",fCentBin),m,pt2, w) ;
 769           }
 770           FillHistogram(Form("hMiSingleDispcore_cen%d",fCentBin),mcore,ptcore2, w) ;
 771         }
 772         if(ph1->IsDisp2OK()){
 773           FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),m,pt1, w) ;
 774           FillHistogram(Form("hMiSingleDisp2core_cen%d",fCentBin),mcore,ptcore1, w) ;
 775         }
 776         if(ph2->IsDisp2OK()){
 777           FillHistogram(Form("hMiSingleDisp2_cen%d",fCentBin),m,pt2, w) ;
 778           FillHistogram(Form("hMiSingleDisp2core_cen%d",fCentBin),mcore,ptcore2, w) ;
 779         }
 780         if(ph1->IsDispOK() && ph1->IsCPVOK()){
 781           snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin) ;
 782           FillHistogram(key,m,pt1, w) ;
 783           snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin);
 784           FillHistogram(key,mcore,ptcore1, w) ;
 785         }
 786         if(ph2->IsDispOK() && ph2->IsCPVOK()){
 787           snprintf(key,55,"hMiSingleBoth_cen%d",fCentBin);
 788           FillHistogram(key,m,pt2, w) ;
 789           snprintf(key,55,"hMiSingleBothcore_cen%d",fCentBin);
 790           FillHistogram(key,mcore,ptcore2, w) ;
 791         }
 792         if(ph1->IsDisp2OK() && ph1->IsCPVOK()){
 793           FillHistogram(Form("hMiSingleBoth2_cen%d",fCentBin),m,pt1, w) ;
 794           FillHistogram(Form("hMiSingleBoth2core_cen%d",fCentBin),mcore,ptcore1, w) ;
 795         }
 796         if(ph2->IsDisp2OK() && ph2->IsCPVOK()){
 797           FillHistogram(Form("hMiSingleBoth2_cen%d",fCentBin),m,pt2, w) ;
 798           FillHistogram(Form("hMiSingleBoth2core_cen%d",fCentBin),mcore,ptcore2, w) ;
 799         }
 800
 801
 802
 803         if(a<kAlphaCut){
 804           FillHistogram(Form("hMiPi0All_a07_cen%d",fCentBin),m,pt, w) ;
 805         }
 806         if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 807           FillHistogram(Form("hMiMassPtV0ACPV_cen%d",fCentBin),m,pt,dphiA, w) ;
 808           FillHistogram(Form("hMiMassPtV0CCPV_cen%d",fCentBin),m,pt,dphiC, w) ;
 809           if(fHaveTPCRP)
 810             FillHistogram(Form("hMiMassPtTPCCPV_cen%d",fCentBin),m,pt,dphiT, w) ;
 811
 812           FillHistogram(Form("hMiMassPtV0ACPVcore_cen%d",fCentBin),mcore, ptcore,dphiA, w) ;
 813           FillHistogram(Form("hMiMassPtV0CCPVcore_cen%d",fCentBin),mcore, ptcore,dphiC, w) ;
 814           if(fHaveTPCRP)
 815             FillHistogram(Form("hMiMassPtTPCCPVcore_cen%d",fCentBin),mcore, ptcore,dphiT, w) ;
 816
 817           FillHistogram(Form("hMiPi0CPV_cen%d",fCentBin),m,pt, w) ;
 818           FillHistogram(Form("hMiPi0CPVcore_cen%d",fCentBin),mcore, ptcore, w) ;
 819
 820           if(a<kAlphaCut){
 821             FillHistogram(Form("hMiPi0CPV_a07_cen%d",fCentBin),m,pt, w) ;
 822           }
 823         }
 824         if(ph1->IsCPV2OK() && ph2->IsCPV2OK()){
 825           FillHistogram(Form("hMiPi0CPV2_cen%d",fCentBin),m,pt, w) ;
 826           FillHistogram(Form("hMiPi0CPV2core_cen%d",fCentBin),mcore, ptcore, w) ;
 827
 828           FillHistogram(Form("hMiMassPtV0ACPV2_cen%d",fCentBin),m,pt,dphiA, w) ;
 829           FillHistogram(Form("hMiMassPtV0CCPV2_cen%d",fCentBin),m,pt,dphiC, w) ;
 830           if(fHaveTPCRP)
 831             FillHistogram(Form("hMiMassPtTPCCPV2_cen%d",fCentBin),m,pt,dphiT, w) ;
 832           FillHistogram(Form("hMiMassPtV0ACPV2core_cen%d",fCentBin),mcore,ptcore,dphiA, w) ;
 833           FillHistogram(Form("hMiMassPtV0CCPV2core_cen%d",fCentBin),mcore,ptcore,dphiC, w) ;
 834           if(fHaveTPCRP)
 835             FillHistogram(Form("hMiMassPtTPCCPV2core_cen%d",fCentBin),mcore,ptcore,dphiT, w) ;
 836
 837           if(a<kAlphaCut){
 838             FillHistogram(Form("hMiPi0CPV2_a07_cen%d",fCentBin),m,pt, w) ;
 839           }
 840         }
 841         if(ph1->IsDispOK() && ph2->IsDispOK()){
 842           FillHistogram(Form("hMiMassPtV0ADisp_cen%d",fCentBin),m,pt,dphiA, w) ;
 843           FillHistogram(Form("hMiMassPtV0CDisp_cen%d",fCentBin),m,pt,dphiC, w) ;
 844           if(fHaveTPCRP)
 845             FillHistogram(Form("hMiMassPtTPCDisp_cen%d",fCentBin),m,pt,dphiT, w) ;
 846
 847           FillHistogram(Form("hMiMassPtV0ADispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA, w) ;
 848           FillHistogram(Form("hMiMassPtV0CDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC, w) ;
 849           if(fHaveTPCRP)
 850             FillHistogram(Form("hMiMassPtTPCDispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT, w) ;
 851
 852
 853           FillHistogram(Form("hMiPi0Disp_cen%d",fCentBin),m,pt, w) ;
 854           FillHistogram(Form("hMiPi0Dispcore_cen%d",fCentBin),pv12.M(),pv12.Pt(), w) ;
 855           if(ph1->IsntUnfolded() && ph2->IsntUnfolded()){
 856             FillHistogram(Form("hMiPi0Dispwou_cen%d",fCentBin),m,pt, w) ;
 857             FillHistogram(Form("hMiMassPtV0ADispwou_cen%d",fCentBin),m,pt,dphiA, w) ;
 858             FillHistogram(Form("hMiMassPtV0CDispwou_cen%d",fCentBin),m,pt,dphiC, w) ;
 859             if(fHaveTPCRP)
 860               FillHistogram(Form("hMiMassPtTPCDispwou_cen%d",fCentBin),m,pt,dphiT, w) ;
 861           }
 862
 863           if(a<kAlphaCut){
 864             FillHistogram(Form("hMiPi0Disp_a07_cen%d",fCentBin),m,pt, w) ;
 865           }
 866           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 867             FillHistogram(Form("hMiMassPtV0ABoth_cen%d",fCentBin),m,pt,dphiA, w) ;
 868             FillHistogram(Form("hMiMassPtV0CBoth_cen%d",fCentBin),m,pt,dphiC, w) ;
 869             if(fHaveTPCRP)
 870               FillHistogram(Form("hMiMassPtTPCBoth_cen%d",fCentBin),m,pt,dphiT, w) ;
 871
 872             FillHistogram(Form("hMiMassPtV0ABothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA, w) ;
 873             FillHistogram(Form("hMiMassPtV0CBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC, w) ;
 874             if(fHaveTPCRP)
 875               FillHistogram(Form("hMiMassPtTPCBothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT, w) ;
 876
 877             FillHistogram(Form("hMiPi0Both_cen%d",fCentBin),m,pt, w) ;
 878             FillHistogram(Form("hMiPi0Bothcore_cen%d",fCentBin),pv12.M(),pv12.Pt(), w) ;
 879
 880             if(a<kAlphaCut){
 881               FillHistogram(Form("hMiPi0Both_a07_cen%d",fCentBin),m,pt, w) ;
 882             }
 883           }
 884         }
 885
 886         if(ph1->IsDisp2OK() && ph2->IsDisp2OK()){
 887           FillHistogram(Form("hMiMassPtV0ADisp2_cen%d",fCentBin),m,pt,dphiA, w) ;
 888           FillHistogram(Form("hMiMassPtV0CDisp2_cen%d",fCentBin),m,pt,dphiC, w) ;
 889           if(fHaveTPCRP)
 890             FillHistogram(Form("hMiMassPtTPCDisp2_cen%d",fCentBin),m,pt,dphiT, w) ;
 891
 892           FillHistogram(Form("hMiMassPtV0ADisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA, w) ;
 893           FillHistogram(Form("hMiMassPtV0CDisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC, w) ;
 894           if(fHaveTPCRP)
 895             FillHistogram(Form("hMiMassPtTPCDisp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT, w) ;
 896
 897
 898           FillHistogram(Form("hMiPi0Disp2_cen%d",fCentBin),m,pt, w) ;
 899           FillHistogram(Form("hMiPi0Disp2core_cen%d",fCentBin),pv12.M(),pv12.Pt(), w) ;
 900
 901           if(ph1->IsCPVOK() && ph2->IsCPVOK()){
 902             FillHistogram(Form("hMiMassPtV0ABoth2_cen%d",fCentBin),m,pt,dphiA, w) ;
 903             FillHistogram(Form("hMiMassPtV0CBoth2_cen%d",fCentBin),m,pt,dphiC, w) ;
 904             if(fHaveTPCRP)
 905               FillHistogram(Form("hMiMassPtTPCBoth2_cen%d",fCentBin),m,pt,dphiT, w) ;
 906
 907             FillHistogram(Form("hMiMassPtV0ABoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiA, w) ;
 908             FillHistogram(Form("hMiMassPtV0CBoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiC, w) ;
 909             if(fHaveTPCRP)
 910               FillHistogram(Form("hMiMassPtTPCBoth2core_cen%d",fCentBin),pv12.M(),pv12.Pt(),dphiT, w) ;
 911
 912             FillHistogram(Form("hMiPi0Both2_cen%d",fCentBin),m,pt, w) ;
 913             FillHistogram(Form("hMiPi0Both2core_cen%d",fCentBin),pv12.M(),pv12.Pt(), w) ;
 914
 915           }
 916         }
 917       } // end of loop i2
 918     }
 919   } // end of loop i1
 920 }
 921
 922
 923 void AliAnalysisTaskPi0FlowMC::ProcessMC()
 924 {
 925   FillMCHist();
 926   FillSecondaries() ;
 927 }
 928
 929
 930
 931 //___________________________________________________________________________
 932 void AliAnalysisTaskPi0FlowMC::FillMCHist(){
 933   //fill histograms for efficiensy etc. calculation
 934
 935   //---------First pi0/eta-----------------------------
 936   char partName[10] ;
 937   char hkey[55] ;
 938
 939   if(!fStack) return ;
 940   for(Int_t i=0;i<fStack->GetNtrack();i++){
 941      TParticle* particle =  fStack->Particle(i);
 942     if(particle->GetPdgCode() == kPi0)
 943       snprintf(partName,10,"pi0") ;
 944     else
 945       if(particle->GetPdgCode() == kEta)
 946         snprintf(partName,10,"eta") ;
 947       else
 948         if(particle->GetPdgCode() == kGamma)
 949            snprintf(partName,10,"gamma") ;
 950         else
 951           if(particle->GetPdgCode() == 310)
 952             snprintf(partName,10,"K0S") ;
 953           else
 954             continue ;
 955
 956     //Primary particle
 957     Double_t r=particle->R() ;
 958     Double_t pt = particle->Pt() ;
 959     //Distribution over vertex
 960     FillHistogram(Form("hMC_%s_vertex",partName),pt,r) ;
 961
 962     if(r >kRCut)
 963       continue ;
 964
 965     Double_t phi=particle->Phi() ;
 966     while(phi<0.)phi+=TMath::TwoPi() ;
 967     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 968
 969     Double_t phig = 180./TMath::Pi()*phi; // phi in deg
 970
 971     //Total number of pi0 with creation radius <1 cm
 972     Double_t weight = PrimaryParticleWeight(particle) ;
 973     snprintf(hkey,55,"hMC_all_%s_cen%d",partName,fCentBin) ;
 974
 975     FillHistogram(hkey,pt,weight) ;
 976
 977     if(TMath::Abs(particle->Y())<0.135 && phig>260. && phig<320.){
 978       snprintf(hkey,55,"hMC_unitEta_%s_cen%d",partName,fCentBin) ;
 979       FillHistogram(hkey,pt,weight) ;
 980
 981       snprintf(hkey,55,"hMC_rap_%s_cen%d",partName,fCentBin) ;
 982       FillHistogram(hkey,particle->Y(),weight) ;
 983
 984       snprintf(hkey,55,"hMC_phi_%s_cen%d",partName,fCentBin) ;
 985       FillHistogram(hkey,phi,weight) ;
 986     }
 987
 988   }
 989 }
 990
 991 //________________________________________________________________________
 992 void AliAnalysisTaskPi0FlowMC::FillSecondaries(){
 993   //Sort secondaires
 994
 995   //Fill spectra of primary particles
 996   //with proper weight
 997   if( fDebug )
 998     AliInfo("start");
 999   for(Int_t i=0; i<fStack->GetNtrack(); i++){
1000     TParticle * p = fStack->Particle(i) ;
1001     if(p->R()>kRCut)
1002       continue ;
1003     if( TMath::Abs(p->Pt())<1.e-6 )
1004       continue;
1005     if(TMath::Abs(p->Y())>0.5)
1006       continue ;
1007     Double_t w = PrimaryParticleWeight(p) ;
1008     Int_t primPdgCode=p->GetPdgCode() ;
1009       switch(primPdgCode){
1010         case  kGamma: FillHistogram(Form("hPrimPhot_cen%d",fCentBin),p->Pt(),w);
1011                   break ;
1012         case  kElectron:
1013         case -kElectron:
1014                   FillHistogram(Form("hPrimEl_cen%d",fCentBin),p->Pt(),w);
1015                   break ;
1016         case  kPi0:
1017                   FillHistogram(Form("hPrimPi0_cen%d",fCentBin),p->Pt(),w);
1018                   break ;
1019         case  kEta:
1020                   FillHistogram(Form("hPrimEta_cen%d",fCentBin),p->Pt(),w);
1021                   break ;
1022         case  kPiPlus:
1023         case  kPiMinus:
1024                   FillHistogram(Form("hPrimPipm_cen%d",fCentBin),p->Pt(),w);
1025                   break ;
1026         case  kProton:  //p
1027                   FillHistogram(Form("hPrimP_cen%d",fCentBin),p->Pt(),w);
1028                   break ;
1029         case kProtonBar:  //pbar
1030                   FillHistogram(Form("hPrimPbar_cen%d",fCentBin),p->Pt(),w);
1031                   break ;
1032         case  kNeutron:  //n
1033                   FillHistogram(Form("hPrimN_cen%d",fCentBin),p->Pt(),w);
1034                   break ;
1035         case  kNeutronBar:  //nbar
1036                   FillHistogram(Form("hPrimNbar_cen%d",fCentBin),p->Pt(),w);
1037                   break ;
1038         case  310:  //nbar
1039                   FillHistogram(Form("hPrimK0S_cen%d",fCentBin),p->Pt(),w);
1040                   break ;
1041         case  130:  //nbar
1042                   FillHistogram(Form("hPrimK0L_cen%d",fCentBin),p->Pt(),w);
1043                   break ;
1044         case  321:  //K+
1045         case -321:  //K-
1046                   FillHistogram(Form("hPrimKpm_cen%d",fCentBin),p->Pt(),w);
1047                   break ;
1048         default:           //other
1049                   FillHistogram(Form("hPrimOther_cen%d",fCentBin),p->Pt(),w);
1050       }
1051   }
1052   if(fDebug)
1053     AliInfo("Origins of secondary pi0s");
1054   //Origins of secondary pi0s
1055   for(Int_t i=0; i<fStack->GetNtrack(); i++){
1056     TParticle * p = fStack->Particle(i) ;
1057     if(p->GetPdgCode()!=111)
1058       continue ;
1059     FillHistogram("Vertex",p->Pt(),p->R());
1060     if(p->R()<kRCut)
1061       continue ;
1062     Double_t phi=p->Phi() ;
1063     while(phi<0.)phi+=TMath::TwoPi() ;
1064     while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
1065     FillHistogram("hSecondPi0RphiZ",p->R(),phi,p->Vz()) ;
1066     Double_t w = PrimaryParticleWeight(p) ;
1067     FillHistogram("hSecondPi0RE",p->R(),p->Pt(),w) ;
1068   }
1069
1070   TLorentzVector p1;
1071
1072   Int_t inPHOS=fCaloPhotonsPHOS->GetEntries() ;
1073   for (Int_t i1=0; i1<inPHOS-1; i1++) {
1074     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1075     Double_t w1=ph1->GetWeight() ;
1076     for (Int_t i2=i1+1; i2<inPHOS; i2++) {
1077       AliCaloPhoton * ph2=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i2) ;
1078       TLorentzVector p12  = *ph1  + *ph2;
1079       Double_t w2=ph2->GetWeight() ;
1080       Double_t w = TMath::Sqrt(w1*w2) ;
1081       FillHistogram(Form("hParentAll_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1082       Int_t prim=FindCommonParent(ph1->GetPrimary(),ph2->GetPrimary()) ;
1083       if(prim>-1){
1084         TParticle * particle = (TParticle *)fStack->Particle(prim);
1085         FillHistogram("hMass_R",p12.M(),p12.Pt(),TMath::Sqrt(particle->R()*particle->R()+particle->Vz()*particle->Vz())) ;
1086
1087
1088         Int_t pdgCode=particle->GetPdgCode() ;
1089         if(pdgCode!=111){ //common parent - not pi0
1090           if(pdgCode==22)
1091             FillHistogram(Form("hParentGamma_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1092           else{
1093             if(pdgCode==11 || pdgCode==-11){
1094               FillHistogram(Form("hParentEl_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1095             }
1096             else{
1097               if(InPi0mass(p12.M() ,p12.Pt())){
1098                 if(fDebug >1) AliInfo(Form("Common parent: %d",pdgCode));
1099               }
1100               FillHistogram(Form("hParentOther_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1101             }
1102           }//Not photons
1103         }//Common parent not pi0
1104         else{ //common parent - pi0
1105           FillHistogram(Form("hParentPi0_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1106           FillHistogram(Form("Real_pi_R"),p12.M(),p12.Pt(),particle->R(),w) ;
1107           FillHistogram(Form("Real_pi_Z"),p12.M(),p12.Pt(),particle->Vz(),w) ;
1108           if(particle->R()<kRCut && TMath::Abs(particle->Vz())<fMaxAbsVertexZ){
1109             FillHistogram(Form("hParentDirPi0_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1110             continue ;
1111           }
1112           //Common particle pi0, created off-vertex
1113           Int_t primPi0=particle->GetFirstMother();
1114           if(primPi0==-1){
1115             FillHistogram(Form("hParentPi0NoPrim_cen%d",fCentBin),p12.M(),p12.Pt(),w) ;
1116           }
1117           else{
1118             Int_t primPdgCode=fStack->Particle(primPi0)->GetPdgCode() ;
1119             switch(primPdgCode){
1120             case 221: FillHistogram(Form("hParentPi0Eta_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //eta
1121                       break ;
1122             case 223: FillHistogram(Form("hParentPi0Omega_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //omega
1123                       break ;
1124             case  211:  //pi+-
1125             case -211: FillHistogram(Form("hParentPi0Pipm_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
1126                       break ;
1127             case  321:  //K+-
1128             case -321: FillHistogram(Form("hParentPi0Kpm_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
1129                       break ;
1130             case 310: FillHistogram(Form("hParentPi0Ks_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // K0S
1131                       break ;
1132             case 130: FillHistogram(Form("hParentPi0Kl_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // K0L
1133                       break ;
1134             case  2212:  //p
1135             case  2112:  //n
1136                       FillHistogram(Form("hParentPi0pn_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // pn
1137                       break ;
1138             case -2212:  //pbar
1139             case -2112:  //nbar
1140                       FillHistogram(Form("hParentPi0antipn_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; // pn
1141                       break ;
1142             default:       //other
1143                       FillHistogram(Form("hParentPi0Other_cen%d",fCentBin),p12.M(),p12.Pt(),w) ; //
1144             }//switch
1145           }//pi0 with primary
1146         }//common parent - pi0
1147       }//there is common primary
1148     }//seond photon loop
1149   }//first photon loop
1150
1151
1152   //Now look at photon contaiminations
1153   for (Int_t i1=0; i1<inPHOS-1; i1++) {
1154     AliCaloPhoton * ph1=(AliCaloPhoton*)fCaloPhotonsPHOS->At(i1) ;
1155     Int_t iprim = ph1->GetPrimary() ;
1156     if(iprim<0)
1157       FillAllHistograms(Form("hGammaNoPrim_cen%d",fCentBin),ph1) ; //
1158     else{
1159       //Find primary at vertex
1160       TParticle * primPHOS = fStack->Particle(iprim) ;
1161       Int_t iprimV=primPHOS->GetFirstMother();
1162       TParticle * primVtx = primPHOS ;
1163       while((iprimV>-1) && primVtx->R()>kRCut){
1164         primVtx = fStack->Particle(iprimV) ;
1165         iprimV=primVtx->GetFirstMother();
1166       }
1167
1168       //photon
1169       Int_t primPdgCode=primVtx->GetPdgCode() ;
1170       switch(primPdgCode){
1171         case  22: FillAllHistograms("hGammaPhot",ph1);
1172                   break ;
1173         case  11:
1174         case -11:
1175                   FillAllHistograms("hGammaEl",ph1);
1176                   break ;
1177         case  111:
1178                   FillAllHistograms("hGammaPi0",ph1);
1179                   break ;
1180         case  221:
1181                   FillAllHistograms("hGammaEta",ph1);
1182                   break ;
1183         case 223: FillAllHistograms("hGammaOmega",ph1) ; //omega
1184                   break ;
1185         case  211:
1186         case -211:
1187                   FillAllHistograms("hGammaPipm",ph1);
1188                   //Find particle entered PHOS
1189                   if(primVtx == primPHOS)
1190                     FillAllHistograms("hGammaPipmDirect",ph1);
1191                   else{
1192                     Int_t primPdgPHOS=primPHOS->GetPdgCode() ;
1193                     if(primPdgPHOS==22){
1194                        FillAllHistograms("hGammaPipmGamma",ph1);
1195                        FillHistogram("hPipmGammaConvR",ph1->Pt(),primPHOS->R());
1196                        FillHistogram("hPipmGammaConvRZ",primPHOS->Vz(),primPHOS->R());
1197                        break ;
1198                     }
1199                     if(TMath::Abs(primPdgPHOS)==11){
1200                        FillAllHistograms("hGammaPipmEl",ph1);
1201                        FillHistogram("hPipmElConvR",ph1->Pt(),primPHOS->R());
1202                        break ;
1203                     }
1204                     if(TMath::Abs(primPdgPHOS)==2212){
1205                        FillAllHistograms("hGammaPipmp",ph1);
1206                        FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
1207                        break ;
1208                     }
1209                     if(TMath::Abs(primPdgPHOS)==2112){
1210                        FillAllHistograms("hGammaPipmn",ph1);
1211                        FillHistogram("hPipmNConvR",ph1->Pt(),primPHOS->R());
1212                        break ;
1213                     }
1214                     FillAllHistograms("hGammaPipmOther",ph1);
1215                     FillHistogram("hPipmOtherConvR",ph1->Pt(),primPHOS->R());
1216                   }
1217                   break ;
1218         case  2212:  //p
1219                   FillAllHistograms("hGammaP",ph1);
1220                   break ;
1221         case -2212:  //pbar
1222                   FillAllHistograms("hGammaPbar",ph1);
1223                   break ;
1224         case  2112:  //n
1225                   FillAllHistograms("hGammaN",ph1);
1226                   break ;
1227         case -2112:  //nbar
1228                   FillAllHistograms("hGammaNbar",ph1) ; // pn
1229                   break ;
1230         case  310:  //nbar
1231                   FillAllHistograms("hGammaK0S",ph1) ; // pn
1232                   break ;
1233         case  130:  //nbar
1234                   FillAllHistograms("hGammaK0L",ph1) ; // pn
1235                   break ;
1236         case  321:  //K+
1237         case -321:  //K-
1238                   FillAllHistograms("hGammaKpm",ph1) ; // pn
1239                   break ;
1240         case -323:
1241         case  323:
1242         case -313:
1243         case  313: FillAllHistograms("hGammaKstar",ph1) ; // K*(892)
1244                   break ;
1245
1246         case -2224 : //Deltas
1247         case  2224 : //Deltas
1248         case -2214 : //Deltas
1249         case  2214 : //Deltas
1250         case -2114 : //Deltas
1251         case  2114 : //Deltas
1252         case -1114 : //Deltas
1253         case  1114 : //Deltas
1254                   FillAllHistograms("hGammaDelta",ph1) ; // pn
1255                   break ;
1256         default:           //other
1257             if(primVtx->GetPDG()->Charge())
1258               FillAllHistograms("hGammaOtherCharged",ph1) ; //
1259             else
1260               FillAllHistograms("hGammaOtherNeutral",ph1) ; //
1261       }
1262     }
1263
1264   }//single photons
1265 }
1266
1267 //_____________________________________________________________________________
1268 void AliAnalysisTaskPi0FlowMC::FillAllHistograms(const char * particleName,AliCaloPhoton * ph)
1269 {
1270   //Fill All PID histograms
1271
1272   Double_t w=ph->GetWeight() ;
1273   Double_t pt = ph->Pt() ;
1274   Double_t ptC=ph->GetMomV2()->Pt() ;
1275   FillHistogram(Form("%s_All_cen%d",particleName,fCentBin),pt,w) ;
1276   FillHistogram(Form("%s_Allcore_cen%d",particleName,fCentBin),ptC,w) ;
1277   if(ph->IsCPVOK()){
1278     FillHistogram(Form("%s_CPV_cen%d",particleName,fCentBin),pt,w) ;
1279     FillHistogram(Form("%s_CPVcore_cen%d",particleName,fCentBin),ptC,w) ;
1280   }
1281   if(ph->IsCPV2OK()){
1282     FillHistogram(Form("%s_CPV2_cen%d",particleName,fCentBin),pt,w) ;
1283     FillHistogram(Form("%s_CPV2core_cen%d",particleName,fCentBin),ptC,w) ;
1284   }
1285   if(ph->IsDispOK()){
1286     FillHistogram(Form("%s_Disp_cen%d",particleName,fCentBin),pt,w) ;
1287     FillHistogram(Form("%s_Dispcore_cen%d",particleName,fCentBin),ptC,w) ;
1288     if(ph->IsDisp2OK()){
1289       FillHistogram(Form("%s_Disp2_cen%d",particleName,fCentBin),pt,w) ;
1290       FillHistogram(Form("%s_Disp2core_cen%d",particleName,fCentBin),ptC,w) ;
1291     }
1292     if(ph->IsCPVOK()){
1293       FillHistogram(Form("%s_Both_cen%d",particleName,fCentBin),pt,w) ;
1294       FillHistogram(Form("%s_Bothcore_cen%d",particleName,fCentBin),ptC,w) ;
1295     }
1296   }
1297 }
1298
1299
1300 //___________________________________________________________________________
1301 Double_t AliAnalysisTaskPi0FlowMC::PrimaryWeight(Int_t /* primary */){
1302   return 1.;
1303 }
1304 //________________________________________________________________________
1305 Double_t AliAnalysisTaskPi0FlowMC::PrimaryParticleWeight(TParticle * /* particle */){
1306   return 1.;
1307 }
1308
1309 //___________________________________________________________________________
1310 AliStack* AliAnalysisTaskPi0FlowMC::GetMCStack()
1311 {
1312   fStack = 0;
1313   AliVEventHandler* eventHandler = AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler();
1314   if(eventHandler){
1315     AliMCEventHandler* mcEventHandler = dynamic_cast<AliMCEventHandler*> (eventHandler);
1316     if( mcEventHandler)
1317       fStack = static_cast<AliMCEventHandler*>(AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler())->MCEvent()->Stack();
1318   }
1319
1320   if( ! fStack )
1321     AliError("Could not get MC Stack!");
1322   return fStack;
1323 }
1324
1325 Int_t AliAnalysisTaskPi0FlowMC::FindPrimary(AliVCluster*clu,  Bool_t&sure){
1326   //Finds primary and estimates if it unique one?
1327   //First check can it be photon/electron
1328   const Double_t emFraction=0.9; //part of energy of cluster to be assigned to EM particle
1329   Int_t n=clu->GetNLabels() ;
1330   for(Int_t i=0;  i<n;  i++){
1331     TParticle*  p=  fStack->Particle(clu->GetLabelAt(i)) ;
1332     Int_t pdg = p->GetPdgCode() ;
1333     if(pdg==22  ||  pdg==11 || pdg == -11){
1334       if(p->Energy()>emFraction*clu->E()){
1335         sure=kTRUE ;
1336         return clu->GetLabelAt(i);
1337       }
1338     }
1339   }
1340
1341   Double_t*  Ekin=  new  Double_t[n] ;
1342   for(Int_t i=0;  i<n;  i++){
1343     TParticle*  p=  fStack->Particle(clu->GetLabelAt(i)) ;
1344     Ekin[i]=p->P() ;  // estimate of kinetic energy
1345     if(p->GetPdgCode()==-2212  ||  p->GetPdgCode()==-2112){
1346       Ekin[i]+=1.8  ;  //due to annihilation
1347     }
1348   }
1349   Int_t iMax=0;
1350   Double_t eMax=0.,eSubMax=0. ;
1351   for(Int_t i=0;  i<n;  i++){
1352     if(Ekin[i]>eMax){
1353       eSubMax=eMax;
1354       eMax=Ekin[i];
1355       iMax=i;
1356     }
1357   }
1358   if(eSubMax>0.8*eMax)//not obvious primary
1359     sure=kFALSE;
1360   else
1361     sure=kTRUE;
1362   delete[]  Ekin;
1363   return  clu->GetLabelAt(iMax) ;
1364 }
1365
1366 //________________________________________________________________________
1367 Int_t AliAnalysisTaskPi0FlowMC::FindCommonParent(Int_t iPart, Int_t jPart){
1368   //check if there is a common parent for particles i and j
1369   // -1: no common parent or wrong iPart/jPart
1370
1371   if(iPart==-1 || iPart>=fStack->GetNtrack() ||
1372      jPart==-1 || jPart>=fStack->GetNtrack()) return -1;
1373
1374   Int_t iprim1=iPart;
1375   while(iprim1>-1){
1376      Int_t iprim2=jPart;
1377      while(iprim2>-1){
1378        if(iprim1==iprim2)
1379          return iprim1 ;
1380        iprim2=((TParticle *)fStack->Particle(iprim2))->GetFirstMother();
1381      }
1382      iprim1=((TParticle *)fStack->Particle(iprim1))->GetFirstMother();
1383   }
1384   return -1;
1385 }
1386 //________________________________________________________________________
1387 Bool_t AliAnalysisTaskPi0FlowMC::HaveParent(Int_t iPart, Int_t pdgParent){
1388   //check if there is a common parent for particles i and j
1389   // -1: no common parent or wrong iPart/jPart
1390
1391   if(iPart==-1 || iPart>=fStack->GetNtrack()) return -1;
1392
1393   Int_t iprim1=iPart;
1394   while(iprim1>-1){
1395     TParticle * tmp = fStack->Particle(iprim1) ;
1396     if(tmp->GetPdgCode()==pdgParent)
1397       return kTRUE ;
1398     iprim1=tmp->GetFirstMother();
1399   }
1400   return kFALSE;
1401 }
1402 //________________________________________________________________________
1403 Bool_t AliAnalysisTaskPi0FlowMC::InPi0mass(Double_t m, Double_t /*pt*/){
1404
1405  return TMath::Abs(m-0.135)<0.007*2.5 ;
1406 }
1407