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