PWG4/PartCorrDep/AliAnalysisTaskTaggedPhotons.cxx

   1 /**************************************************************************\r
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
   3  *                                                                        *\r
   4  * Author: The ALICE Off-line Project.                                    *\r
   5  * Contributors are mentioned in the code where appropriate.              *\r
   6  *                                                                        *\r
   7  * Permission to use, copy, modify and distribute this software and its   *\r
   8  * documentation strictly for non-commercial purposes is hereby granted   *\r
   9  * without fee, provided that the above copyright notice appears in all   *\r
  10  * copies and that both the copyright notice and this permission notice   *\r
  11  * appear in the supporting documentation. The authors make no claims     *\r
  12  * about the suitability of this software for any purpose. It is          *\r
  13  * provided "as is" without express or implied warranty.                  *\r
  14  **************************************************************************/\r
  15 \r
  16 /* $Id: */\r
  17 \r
  18 //_________________________________________________________________________\r
  19 // Analysis for Tagged Photons\r
  20 // Prepares all necessary histograms for calculation of \r
  21 // the yield of pi0 decay photon in calorimeter:\r
  22 // Marks photons which makes pi0 with some other and
  23 // fill invariant mass distributions for estimate background below pi0
  24 // peak so that estimate proportion of fake pairs.
  25 // Fills as well controll MC histograms with clasification of the photon origin
  26 // and check of the ptoportion of truly tagged photons.
  27 // \r
  28 //\r
  29 //*-- Dmitry Blau \r
  30 //////////////////////////////////////////////////////////////////////////////\r
  31 \r
  32 #include <TH1.h>\r
  33 #include <TH2.h>\r
  34 #include <TH3.h>\r
  35 #include <TFile.h>\r
  36 #include <TROOT.h>\r
  37 \r
  38 #include "AliAnalysisTaskTaggedPhotons.h" \r
  39 #include "AliAnalysisManager.h"\r
  40 #include "AliESDVertex.h"\r
  41 #include "AliESDEvent.h" \r
  42 #include "AliAODEvent.h" \r
  43 #include "AliESDCaloCluster.h" \r
  44 #include "AliAODPWG4Particle.h"\r
  45 #include "AliAnalysisManager.h"\r
  46 #include "AliLog.h"\r
  47 #include "TGeoManager.h"\r
  48 #include "AliMCAnalysisUtils.h"
  49 #include "AliMCEventHandler.h"\r
  50 #include "AliMCEvent.h"\r
  51 #include "AliStack.h"\r
  52 #include "AliPHOSGeoUtils.h"\r
  53 #include "AliEMCALGeoUtils.h"\r
  54
  55
  56 \r
  57 //______________________________________________________________________________\r
  58 AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons() : AliAnalysisTaskSE(),
  59   fPHOSgeom(0x0),
  60   fEMCALgeom(0x0),
  61   fStack(0x0),fPHOS(1),
  62   fPhotonId(1.0),fMinEnergyCut(0.4),
  63   fPi0MeanP0(0.136),fPi0MeanP1(0.0),fPi0MeanP2(0.0),fPi0MeanP3(0.0),
  64   fPi0SigmaP0(0.004508),fPi0SigmaP1(0.005497),fPi0SigmaP2(0.00000006),
  65   fZmax(0.),fZmin(0.),fPhimax(0.),fPhimin(0.),
  66   fOutputList(0x0),fEventList(0x0),
  67   fhRecAll(0x0),fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
  68   fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
  69   fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
  70   fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
  71   fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
  72   fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerStack(0x0),
  73   fhDecWMissedPartnerGeom0(0x0),
  74   fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
  75   fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
  76   fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
  77   fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
  78   fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
  79   fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
  80 {\r
  81   //Default constructor\r
  82 }\r
  83 //______________________________________________________________________________\r
  84 AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const char *name) : \r
  85   AliAnalysisTaskSE(name),\r
  86   fPHOSgeom(0x0),
  87   fEMCALgeom(0x0),
  88   fStack(0x0),fPHOS(1),
  89   fPhotonId(1.0),fMinEnergyCut(0.4),
  90   fPi0MeanP0(0.136),fPi0MeanP1(0.0),fPi0MeanP2(0.0),fPi0MeanP3(0.0),
  91   fPi0SigmaP0(0.004508),fPi0SigmaP1(0.005497),fPi0SigmaP2(0.00000006),
  92   fZmax(0.),fZmin(0.),fPhimax(0.),fPhimin(0.),
  93   fOutputList(0x0),fEventList(0x0),
  94   fhRecAll(0x0),fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
  95   fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
  96   fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
  97   fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
  98   fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
  99   fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerStack(0x0),
 100   fhDecWMissedPartnerGeom0(0x0),
 101   fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
 102   fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),fhPartnerMissedGeo(0x0),
 103   fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
 104   fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
 105   fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
 106   fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
 107 {\r
 108   // Constructor.\r
 109 \r
 110   // Output slots \r
 111   DefineOutput(1,  TList::Class()) ; \r
 112 }\r
 113 \r
 114 //____________________________________________________________________________\r
 115 AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) :\r
 116   AliAnalysisTaskSE(ap.GetName()),  \r
 117   fPHOSgeom(0x0),
 118   fEMCALgeom(0x0),
 119   fStack(0x0),fPHOS(ap.fPHOS),
 120   fPhotonId(ap.fPhotonId),fMinEnergyCut(ap.fMinEnergyCut),
 121   fPi0MeanP0(ap.fPi0MeanP0),fPi0MeanP1(ap.fPi0MeanP1),fPi0MeanP2(ap.fPi0MeanP2),fPi0MeanP3(ap.fPi0MeanP3),
 122   fPi0SigmaP0(ap.fPi0SigmaP0),fPi0SigmaP1(ap.fPi0SigmaP1),fPi0SigmaP2(ap.fPi0SigmaP2),
 123   fZmax(ap.fZmax),fZmin(ap.fZmin),fPhimax(ap.fPhimax),fPhimin(ap.fPhimin),
 124   fOutputList(0x0),fEventList(0x0),
 125   fhRecAll(0x0),fhRecAllArea1(0x0),fhRecAllArea2(0x0),fhRecAllArea3(0x0),
 126   fhRecPhoton(0x0),fhRecOther(0x0),fhRecPhotPi0(0x0),fhRecPhotEta(0x0),
 127   fhRecPhotOmega(0x0),fhRecPhotEtapr(0x0),fhRecPhotConv(0x0),fhRecPhotHadron(0x0),
 128   fhRecPhotDirect(0x0),fhRecPhotOther(0x0),
 129   fhDecWMCPartner(0x0),fhDecWMissedPartnerNotPhoton(0x0),fhDecWMissedPartnerAll(0x0),
 130   fhDecWMissedPartnerEmin(0x0),fhDecWMissedPartnerConv(0x0),fhDecWMissedPartnerStack(0x0),
 131   fhDecWMissedPartnerGeom0(0x0),
 132   fhDecWMissedPartnerGeom1(0x0),fhDecWMissedPartnerGeom2(0x0),fhDecWMissedPartnerGeom3(0x0),
 133   fhPartnerMCReg(0x0),fhPartnerMissedEmin(0x0),fhPartnerMissedConv(0x0),
 134   fhPartnerMissedGeo(0x0),
 135   fhTaggedAll(0x0),fhTaggedArea1(0x0),fhTaggedArea2(0x0),fhTaggedArea3(0x0),fhTaggedMult(0x0),
 136   fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
 137   fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
 138   fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
 139 {
 140   // cpy ctor\r
 141 }\r
 142 \r
 143 //_____________________________________________________________________________\r
 144 AliAnalysisTaskTaggedPhotons& AliAnalysisTaskTaggedPhotons::operator = (const AliAnalysisTaskTaggedPhotons& ap)\r
 145 {\r
 146 // assignment operator\r
 147 \r
 148   this->~AliAnalysisTaskTaggedPhotons();\r
 149   new(this) AliAnalysisTaskTaggedPhotons(ap);\r
 150   return *this;\r
 151 }\r
 152 \r
 153 //______________________________________________________________________________\r
 154 AliAnalysisTaskTaggedPhotons::~AliAnalysisTaskTaggedPhotons()\r
 155 {\r
 156   // dtor\r
 157   if(fOutputList) {\r
 158     fOutputList->Clear() ; \r
 159     delete fOutputList ;\r
 160   }\r
 161 }\r
 162 \r
 163 \r
 164 //________________________________________________________________________\r
 165 void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()\r
 166 { \r
 167 \r
 168
 169   //Load geometry\r
 170   //if file "geometry.root" exists, load it\r
 171   //otherwise use misaligned matrixes stored in ESD\r
 172   TFile *geoFile = new TFile("geometry.root","read");\r
 173   if(geoFile->IsZombie()){ //no file, load geo matrixes from ESD\r
 174     AliInfo("Can not find file geometry.root, reading misalignment matrixes from ESD/AOD") ;
 175     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
 176     AliAODEvent * aod = 0x0 ;
 177     if(!esd)
 178       aod=dynamic_cast<AliAODEvent*>(InputEvent()) ;
 179     if(!esd && !aod)
 180       AliFatal("Can not read geometry even from ESD/AOD.") ;
 181     if(fPHOS){//reading PHOS matrixes
 182       fPHOSgeom = new AliPHOSGeoUtils("IHEP","");\r
 183       for(Int_t mod=0; mod<5; mod++){
 184         if(esd){
 185           const TGeoHMatrix* m=esd->GetPHOSMatrix(mod) ;
 186           fPHOSgeom->SetMisalMatrix(m, mod) ;
 187         }
 188         else{
 189           const TGeoHMatrix* m=aod->GetHeader()->GetPHOSMatrix(mod) ;
 190           fPHOSgeom->SetMisalMatrix(m, mod) ;
 191         }
 192       }
 193     }
 194     else{ //EMCAL
 195       fEMCALgeom = new AliEMCALGeoUtils("");
 196       for(Int_t mod=0; mod < 12; mod++){ //<---Gustavo, could you check???
 197         if(esd){
 198           const TGeoHMatrix* m=esd->GetEMCALMatrix(mod) ;
 199           fEMCALgeom->SetMisalMatrix(m, mod) ;
 200         }
 201         else{
 202           const TGeoHMatrix* m=aod->GetHeader()->GetEMCALMatrix(mod) ;
 203           fEMCALgeom->SetMisalMatrix(m, mod) ;
 204         }
 205       }
 206     }
 207   }\r
 208   else{\r
 209     gGeoManager = (TGeoManager*)geoFile->Get("Geometry");\r
 210     //Geometry will be misaligned from GeoManager
 211     if(fPHOS){
 212       fPHOSgeom = new AliPHOSGeoUtils("IHEP","");\r
 213     }
 214     else{
 215       fEMCALgeom = new AliEMCALGeoUtils("EMCAL_COMPLETE");\r
 216     }
 217   }\r
 218 \r
 219 \r
 220   if(fPHOSgeom==NULL && fEMCALgeom==NULL){\r
 221     AliError("Error loading Geometry\n");\r
 222   }\r
 223   else\r
 224     AliInfo("Geometry loaded... OK\n");\r
 225 \r
 226   //Evaluate active PHOS/EMCAL area\r
 227   if(fZmax<=fZmin){ //not set yet
 228     if(fPHOS){
 229       //Check active modules in the current configuration
 230       if(fPHOSgeom){
 231         fZmax= 999. ;
 232         fZmin=-999. ;
 233         fPhimax=-999. ;
 234         fPhimin= 999. ;
 235         for(Int_t imod=1; imod<=5; imod++){
 236           //Find exact coordinates of PHOS corners
 237           Int_t relId[4]={imod,0,1,1} ;
 238           Int_t absId ;
 239           fPHOSgeom->RelToAbsNumbering(relId,absId) ;
 240           TVector3 pos ;
 241           fPHOSgeom->RelPosInAlice(absId,pos) ;
 242           Double_t phi=pos.Phi() ;
 243           while(phi<0.)phi+=TMath::TwoPi() ;
 244           while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 245           fPhimin=TMath::Min(fPhimin,float(phi)) ;
 246           fZmin=TMath::Max(fZmin,float(pos.Z())) ;
 247           relId[2]=64 ;
 248           relId[3]=56 ;
 249           fPHOSgeom->RelToAbsNumbering(relId,absId) ;
 250           fPHOSgeom->RelPosInAlice(absId,pos) ;
 251           phi=pos.Phi() ;
 252           while(phi<0.)phi+=TMath::TwoPi() ;
 253           while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 254           fPhimax=TMath::Max(fPhimax,float(phi)) ;
 255           fZmax=TMath::Min(fZmax,float(pos.Z())) ;
 256         }
 257       }
 258       else{
 259         //Use approximate range
 260         fZmax= 2.25*56/2. ;\r
 261         fZmin=-2.25*56/2. ;\r
 262         fPhimax=220./180.*TMath::Pi() ;\r
 263         fPhimin=320./180.*TMath::Pi() ;\r
 264       }
 265     }
 266     else{ //Similar for EMCAL <--Gustavo, Could you please have a look?
 267       if(fEMCALgeom){
 268         fZmax= 999. ;
 269         fZmin=-999. ;
 270         fPhimax=-999. ;
 271         fPhimin= 999. ;
 272         for(Int_t imod=0; imod<12; imod++){
 273
 274           //Find exact coordinates of SM corners
 275           Int_t absId = fEMCALgeom->GetAbsCellIdFromCellIndexes(imod, 0, 0);
 276           TVector3 pos ;
 277           //Get the position of this tower.
 278           fEMCALgeom->RelPosCellInSModule(absId,pos);
 279           Double_t phi=pos.Phi() ;
 280           while(phi<0.)phi+=TMath::TwoPi() ;
 281           while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 282           fPhimin=TMath::Min(fPhimin,float(phi)) ;
 283           fZmin=TMath::Max(fZmin,float(pos.Z())) ;
 284           absId = fEMCALgeom->GetAbsCellIdFromCellIndexes(imod, 24, 48);
 285           fEMCALgeom->RelPosCellInSModule(absId,pos);
 286           phi=pos.Phi() ;
 287           while(phi<0.)phi+=TMath::TwoPi() ;
 288           while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
 289           fPhimax=TMath::Max(fPhimax,float(phi)) ;
 290           fZmax=TMath::Min(fZmax,float(pos.Z())) ;
 291
 292         }
 293       }
 294       else{
 295         //Use approximate range
 296         fZmax= 325. ;
 297         fZmin=-325. ;
 298         fPhimax=80./180.*TMath::Pi() ;
 299         fPhimin=180./180.*TMath::Pi() ;
 300       }
 301     }
 302   }
 303 \r
 304 \r
 305   // Create the outputs containers\r
 306 \r
 307   OpenFile(1) ; \r
 308   const Int_t nPtBins=52 ;\r
 309   Double_t ptBins[nPtBins+1] ;\r
 310   for(Int_t i=0;i<=20;i++)ptBins[i]=0.1*i ; //0-2 GeV:  0.1 GeV/bin\r
 311   for(Int_t i=21;i<=30;i++)ptBins[i]=2.+0.2*(i-20) ; //2-4 GeV:  0.2 GeV/bin                   \r
 312   for(Int_t i=31;i<=38;i++)ptBins[i]=4.+0.5*(i-30) ; //4-8 GeV:  0.5 GeV/bin                   \r
 313   for(Int_t i=39;i<=42;i++)ptBins[i]=8.+1.0*(i-38) ; //8-12 GeV: 1. GeV/bin                   \r
 314   for(Int_t i=43;i<=52;i++)ptBins[i]=12.+2.0*(i-42) ; //12-30 GeV: 2. GeV/bin                   \r
 315 \r
 316 \r
 317   //Reconstructed spectra\r
 318     fhRecAll      = new TH1D("fhRecAll","Spectrum of all reconstructed particles", nPtBins, ptBins ) ;\r
 319     fhRecAllArea1 = new TH1D("fhRecAllArea1","Spectrum of rec particles in Fid. Area 1", nPtBins, ptBins ) ;\r
 320     fhRecAllArea2 = new TH1D("fhRecAllArea2","Spectrum of rec particles in Fid. Area 2", nPtBins, ptBins ) ;\r
 321     fhRecAllArea3 = new TH1D("fhRecAllArea3","Spectrum of rec particles in Fid. Area 3", nPtBins, ptBins ) ;\r
 322 \r
 323     //Sort registered particles spectra according MC information\r
 324     fhRecPhoton   = new TH1D("fhRecPhoton","Spectrum of rec. with primary==22 and no PID criteria", nPtBins, ptBins ) ;\r
 325     fhRecOther    = new TH1D("fhRecOther"," Spectrum of rec. with primary!=22 and no PID criteria", nPtBins, ptBins ) ;\r
 326     fhRecPhotonPID[0]= new TH1D("fhRecPhotonPID0","Spectrum of rec. with primary==22 and PID=0", nPtBins, ptBins ) ;\r
 327     fhRecPhotonPID[1]= new TH1D("fhRecPhotonPID1","Spectrum of rec. with primary==22 and PID=1", nPtBins, ptBins ) ;\r
 328     fhRecPhotonPID[2]= new TH1D("fhRecPhotonPID2","Spectrum of rec. with primary==22 and PID=2", nPtBins, ptBins ) ;\r
 329     fhRecPhotonPID[3]= new TH1D("fhRecPhotonPID3","Spectrum of rec. with primary==22 and PID=3", nPtBins, ptBins ) ;\r
 330     fhRecOtherPID[0]= new TH1D("fhRecOtherPID0","Spectrum of rec. with primary!=22 and PID=0", nPtBins, ptBins ) ;\r
 331     fhRecOtherPID[1]= new TH1D("fhRecOtherPID1","Spectrum of rec. with primary!=22 and PID=1", nPtBins, ptBins ) ;\r
 332     fhRecOtherPID[2]= new TH1D("fhRecOtherPID2","Spectrum of rec. with primary!=22 and PID=2", nPtBins, ptBins ) ;\r
 333     fhRecOtherPID[3]= new TH1D("fhRecOtherPID3","Spectrum of rec. with primary!=22 and PID=3", nPtBins, ptBins ) ;\r
 334     fhRecPhotPi0    = new TH1D("fhRecPhotPi0","Spectrum of rec. photons from pi0 decays", nPtBins, ptBins ) ;\r
 335     fhRecPhotEta    = new TH1D("fhRecPhotEta","Spectrum of rec. photons from eta decays", nPtBins, ptBins ) ;\r
 336     fhRecPhotOmega  = new TH1D("fhRecPhotOmega","Spectrum of rec. photons from omega decays", nPtBins, ptBins ) ;\r
 337     fhRecPhotEtapr  = new TH1D("fhRecPhotEtapr","Spectrum of rec. photons from eta prime decays", nPtBins, ptBins ) ;\r
 338     fhRecPhotConv   = new TH1D("fhRecPhotConv"," Spectrum of rec. photons from conversion", nPtBins, ptBins ) ;\r
 339     fhRecPhotHadron = new TH1D("fhRecPhotHadron","Spectrum of rec. photons from hadron-matter interactions", nPtBins, ptBins ) ;\r
 340     fhRecPhotDirect = new TH1D("fhRecPhotDirect","Spectrum of rec. photons direct or no primary", nPtBins, ptBins ) ;\r
 341     fhRecPhotOther  = new TH1D("fhRecPhotOther","Spectrum of rec. photons from other hadron decays", nPtBins, ptBins ) ;\r
 342 \r
 343     //MC tagging: reasons of partner loss etc.\r
 344     fhDecWMCPartner        = new TH1D("fhDecWMCPartner","pi0 decay photon which partner should be registered according to MC", nPtBins, ptBins ) ;\r
 345     fhDecWMissedPartnerNotPhoton = new TH1D("fhDecWMissedPartnerNotPhoton","Decay photon with missed non-photon partner", nPtBins, ptBins ) ;\r
 346     fhDecWMissedPartnerAll  = new TH1D("fhDecWMissedPartnerAll","Decay photons with partner missed due to some reason", nPtBins, ptBins ) ;\r
 347     fhDecWMissedPartnerEmin = new TH1D("fhDecWMissedPartnerEmin","Decay photons with partner missed due to low energy", nPtBins, ptBins ) ;\r
 348     fhDecWMissedPartnerConv = new TH1D("fhDecWMissedPartnerConv","Decay photons with partner missed due to conversion", nPtBins, ptBins ) ;\r
 349     fhDecWMissedPartnerStack= new TH1D("fhDecWMissedPartnerStack","Decay photons with partner not in Stack", nPtBins, ptBins ) ;\r
 350     fhDecWMissedPartnerGeom0 = new TH1D("fhDecWMissedPartnerGeom0","Decay photons with partner missed due geometry", nPtBins, ptBins ) ;\r
 351     fhDecWMissedPartnerGeom1 = new TH1D("fhDecWMissedPartnerGeom1","Decay photons with partner missed due geometry Fid. area. 1", nPtBins, ptBins ) ;\r
 352     fhDecWMissedPartnerGeom2 = new TH1D("fhDecWMissedPartnerGeom2","Decay photons with partner missed due geometry Fid. area. 2", nPtBins, ptBins ) ;\r
 353     fhDecWMissedPartnerGeom3 = new TH1D("fhDecWMissedPartnerGeom3","Decay photons with partner missed due geometry Fid. area. 3", nPtBins, ptBins ) ;\r
 354 \r
 355     //MC tagging: Decay partners spectra\r
 356     fhPartnerMCReg      = new TH1D("fhPartnerMCReg","Spectrum of decay partners which should be registered (MC)", nPtBins, ptBins ) ;\r
 357     fhPartnerMissedEmin = new TH1D("fhPartnerMissedEmin","Spectrum of decay partners lost due to Emin cut", nPtBins, ptBins ) ;\r
 358     fhPartnerMissedConv = new TH1D("fhPartnerMissedConv","Spectrum of decay partners lost due to conversion", nPtBins, ptBins ) ;\r
 359     fhPartnerMissedGeo  = new TH1D("fhPartnerMissedGeo","Spectrum of decay partners lost due to acceptance", nPtBins, ptBins ) ;\r
 360 \r
 361     //Tagging\r
 362     fhTaggedAll    = new TH1D("fhTaggedAll","Spectrum of all tagged photons", nPtBins, ptBins ) ;\r
 363     fhTaggedArea1  = new TH1D("fhTaggedArea1","Spectrum of all tagged photons Fid. area1", nPtBins, ptBins ) ;\r
 364     fhTaggedArea2  = new TH1D("fhTaggedArea2","Spectrum of all tagged photons Fid. area2", nPtBins, ptBins ) ;\r
 365     fhTaggedArea3  = new TH1D("fhTaggedArea3","Spectrum of all tagged photons Fid. area3", nPtBins, ptBins ) ;\r
 366     fhTaggedPID[0] = new TH1D("fhTaggedPID0","Spectrum of tagged photons for PID=0", nPtBins, ptBins ) ;\r
 367     fhTaggedPID[1] = new TH1D("fhTaggedPID1","Spectrum of tagged photons for PID=1", nPtBins, ptBins ) ;\r
 368     fhTaggedPID[2] = new TH1D("fhTaggedPID2","Spectrum of tagged photons for PID=2", nPtBins, ptBins ) ;\r
 369     fhTaggedPID[3] = new TH1D("fhTaggedPID3","Spectrum of tagged photons for PID=3", nPtBins, ptBins ) ;\r
 370     fhTaggedMult   = new TH1D("fhTaggedMult","Spectrum of multiply tagged photons", nPtBins, ptBins ) ;\r
 371 \r
 372     //Tagging: use MC information if available\r
 373     fhTaggedMCTrue    = new TH1D("fhTaggedMCTrue","Spectrum of correctly tagged pi0 decay photons", nPtBins, ptBins ) ;\r
 374     fhMCMissedTagging = new TH1D("fhMCMissedTagging","Spectrum of pi0 decay photons missed tagging due to wrong pair mass", nPtBins, ptBins ) ;\r
 375     fhMCFakeTagged    = new TH1D("fhMCFakeTagged","Spectrum of photons wrongly tagged according to MC", nPtBins, ptBins ) ;\r
 376 \r
 377     //Invariant mass distributions for fake corrections\r
 378     const Int_t nmass=1000 ;\r
 379     Double_t masses[nmass+1] ;\r
 380     for(Int_t i=0;i<=nmass;i++)masses[i]=0.001*i ;\r
 381     fhInvMassReal = new TH2D("fhInvMassReal","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;\r
 382     fhInvMassMixed  = new TH2D("fhInvMassMixed","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;\r
 383     fhMCMissedTaggingMass= new TH2D("fhMCMissedTaggingMass","Inv mass of pairs missed tagging",nmass,masses,nPtBins, ptBins ) ;\r
 384 \r
 385     //Conversion and annihilation radius distributions\r
 386     fhConversionRadius  = new TH1D("fhConversionRadius","Radis of photon production (conversion)",100,0.,500.) ;\r
 387     fhInteractionRadius = new TH1D("fhInteractionRadius","Radis of photon production (hadron interaction)",100,0.,500.);\r
 388 \r
 389     fhEvents               = new TH1D("hEvents", "Number of Events processed", 1, 0, 1);\r
 390 \r
 391 // create output container\r
 392   \r
 393   fOutputList = new TList() ;\r
 394   fEventList = new TList() ; \r
 395   fOutputList->SetName(GetName()) ; \r
 396 \r
 397   fOutputList->AddAt(fhRecAll,                               0) ;\r
 398   fOutputList->AddAt(fhRecAllArea1,                          1) ;\r
 399   fOutputList->AddAt(fhRecAllArea2,                          2) ;\r
 400   fOutputList->AddAt(fhRecAllArea3,                          3) ;\r
 401 \r
 402   fOutputList->AddAt(fhRecPhoton,                            4) ;\r
 403   fOutputList->AddAt(fhRecOther,                             5) ;\r
 404   fOutputList->AddAt(fhRecPhotonPID[0],                      6) ;\r
 405   fOutputList->AddAt(fhRecPhotonPID[1],                      7) ;\r
 406   fOutputList->AddAt(fhRecPhotonPID[2],                      8) ;\r
 407   fOutputList->AddAt(fhRecPhotonPID[3],                      9) ;\r
 408   fOutputList->AddAt(fhRecOtherPID[0],                      10) ;\r
 409   fOutputList->AddAt(fhRecOtherPID[1],                      11) ;\r
 410   fOutputList->AddAt(fhRecOtherPID[2],                      12) ;\r
 411   fOutputList->AddAt(fhRecOtherPID[3],                      13) ;\r
 412   fOutputList->AddAt(fhRecPhotPi0,                          14) ;\r
 413   fOutputList->AddAt(fhRecPhotEta,                          15) ;\r
 414   fOutputList->AddAt(fhRecPhotOmega,                        16) ;\r
 415   fOutputList->AddAt(fhRecPhotEtapr,                        17) ;\r
 416   fOutputList->AddAt(fhRecPhotConv,                         18) ;\r
 417   fOutputList->AddAt(fhRecPhotHadron,                       19) ;\r
 418   fOutputList->AddAt(fhRecPhotDirect,                       20) ;\r
 419   fOutputList->AddAt(fhRecPhotOther,                        21) ;\r
 420 \r
 421   fOutputList->AddAt(fhDecWMCPartner,                       22) ;\r
 422   fOutputList->AddAt(fhDecWMissedPartnerNotPhoton,          23) ;\r
 423   fOutputList->AddAt(fhDecWMissedPartnerAll,                24) ;\r
 424   fOutputList->AddAt(fhDecWMissedPartnerEmin,               25) ;\r
 425   fOutputList->AddAt(fhDecWMissedPartnerConv,               26) ;\r
 426   fOutputList->AddAt(fhDecWMissedPartnerStack,              27) ;\r
 427   fOutputList->AddAt(fhDecWMissedPartnerGeom0,              28) ;\r
 428   fOutputList->AddAt(fhDecWMissedPartnerGeom1,              29) ;\r
 429   fOutputList->AddAt(fhDecWMissedPartnerGeom2,              30) ;\r
 430   fOutputList->AddAt(fhDecWMissedPartnerGeom3,              31) ;\r
 431 \r
 432   fOutputList->AddAt(fhPartnerMCReg,                        32) ;\r
 433   fOutputList->AddAt(fhPartnerMissedEmin,                   33) ;\r
 434   fOutputList->AddAt(fhPartnerMissedConv,                   34) ;\r
 435   fOutputList->AddAt(fhPartnerMissedGeo,                    35) ;\r
 436 \r
 437   fOutputList->AddAt(fhTaggedAll,                           36) ;\r
 438   fOutputList->AddAt(fhTaggedArea1,                         37) ;\r
 439   fOutputList->AddAt(fhTaggedArea2,                         38) ;\r
 440   fOutputList->AddAt(fhTaggedArea3,                         39) ;\r
 441   fOutputList->AddAt(fhTaggedPID[0],                        40) ;\r
 442   fOutputList->AddAt(fhTaggedPID[1],                        41) ;\r
 443   fOutputList->AddAt(fhTaggedPID[2],                        42) ;\r
 444   fOutputList->AddAt(fhTaggedPID[3],                        43) ;\r
 445   fOutputList->AddAt(fhTaggedMult,                          44) ;\r
 446 \r
 447   fOutputList->AddAt(fhTaggedMCTrue,                        45) ;\r
 448   fOutputList->AddAt(fhMCMissedTagging,                     46) ;\r
 449   fOutputList->AddAt(fhMCFakeTagged,                        47) ;\r
 450 \r
 451   fOutputList->AddAt(fhInvMassReal,                         48) ;\r
 452   fOutputList->AddAt(fhInvMassMixed,                        49) ;\r
 453   fOutputList->AddAt(fhMCMissedTaggingMass,                 50) ;\r
 454 \r
 455   fOutputList->AddAt(fhConversionRadius,                    51) ;\r
 456   fOutputList->AddAt(fhInteractionRadius,                   52) ;\r
 457 \r
 458   fOutputList->AddAt(fhEvents,                              53) ;\r
 459 \r
 460 }\r
 461 \r
 462 //______________________________________________________________________________\r
 463 void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *) \r
 464 {\r
 465   //Fill all histograms
 466
 467   fhEvents->Fill(0.);\r
 468 \r
 469   // Processing of one event\r
 470   if(fDebug>1)\r
 471     AliInfo(Form("\n\n Processing event # %lld",  Entry())) ; \r
 472   AliESDEvent* esd = (AliESDEvent*)InputEvent();\r
 473 \r
 474   //MC stack init\r
 475   AliMCEventHandler* mctruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());\r
 476   fStack = mctruth->MCEvent()->Stack();\r
 477 \r
 478   if(!fStack && gDebug>1)\r
 479     AliInfo("No stack! \n");\r
 480 \r
 481   //************************  PHOS/EMCAL *************************************\r
 482   TRefArray * caloClustersArr  = new TRefArray();  \r
 483   if(fPHOS)
 484     esd->GetPHOSClusters(caloClustersArr);\r
 485   else
 486     esd->GetEMCALClusters(caloClustersArr);\r
 487   const Int_t kNumberOfClusters = caloClustersArr->GetEntries() ;  \r
 488 \r
 489   TClonesArray * fCaloPhotonsArr   = new TClonesArray("AliAODPWG4Particle",kNumberOfClusters);\r
 490   Int_t inList = 0; //counter of caloClusters\r
 491 \r
 492   Int_t cluster ; \r
 493 \r
 494   // loop over Clusters\r
 495   for(cluster = 0 ; cluster < kNumberOfClusters ; cluster++) {\r
 496     AliESDCaloCluster * caloCluster = static_cast<AliESDCaloCluster*>(caloClustersArr->At(cluster)) ;\r
 497   \r
 498     if((fPHOS && !caloCluster->IsPHOS()) ||\r
 499        (!fPHOS && caloCluster->IsPHOS()))\r
 500       continue ; \r
 501 \r
 502     Double_t v[3] ; //vertex ;\r
 503     esd->GetVertex()->GetXYZ(v) ;\r
 504 \r
 505     TLorentzVector momentum ;\r
 506     caloCluster->GetMomentum(momentum, v);\r
 507     new ((*fCaloPhotonsArr)[inList]) AliAODPWG4Particle(momentum.Px(),momentum.Py(),momentum.Pz(),caloCluster->E() );\r
 508     AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(inList));\r
 509     inList++;\r
 510 \r
 511     p->SetCaloLabel(cluster,-1); //This and partner cluster
 512     p->SetDistToBad(Int_t(caloCluster->GetDistanceToBadChannel()));\r
 513
 514     p->SetTag(AliMCAnalysisUtils::kMCUnknown);\r
 515     p->SetTagged(kFALSE);   //Reconstructed pairs found\r
 516     p->SetLabel(caloCluster->GetLabel());\r
 517     Float_t pos[3] ;
 518     caloCluster->GetPosition(pos) ;
 519     p->SetFiducialArea(GetFiducialArea(pos)) ;
 520
 521     //PID criteria
 522     p->SetDispBit(TestDisp(caloCluster->GetM02(),caloCluster->GetM20(),caloCluster->E())) ;
 523     p->SetTOFBit(TestTOF(caloCluster->GetTOF(),caloCluster->E())) ;
 524     p->SetChargedBit(TestCharged(caloCluster->GetEmcCpvDistance(),caloCluster->E())) ;
 525 \r
 526     fhRecAll->Fill( p->Pt() ) ; //All recontructed particles\r
 527     Int_t iFidArea = p->GetFiducialArea(); \r
 528     if(iFidArea>0){\r
 529       fhRecAllArea1->Fill(p->Pt() ) ; \r
 530       if(iFidArea>1){\r
 531         fhRecAllArea2->Fill(p->Pt() ) ; \r
 532         if(iFidArea>2){\r
 533           fhRecAllArea3->Fill(p->Pt() ) ; \r
 534         }\r
 535       }\r
 536     }\r
 537 \r
 538     if(fStack){\r
 539       TParticle * prim = fStack->Particle(caloCluster->GetLabel()) ;\r
 540       if(fDebug>2)\r
 541         printf("Pdgcode = %d\n",prim->GetPdgCode());\r
 542 \r
 543       if(prim->GetPdgCode()!=22){ //not photon\r
 544 //<--DP        p->SetPhoton(kFALSE);\r
 545         fhRecOther->Fill(p->Pt()); //not photons real spectra\r
 546         for(Int_t iPID=0; iPID<4; iPID++){\r
 547           if(p->IsPIDOK(iPID,22))\r
 548             fhRecOtherPID[iPID]->Fill(p->Pt());\r
 549         }\r
 550       }\r
 551       else{ //Primary photon (as in MC)\r
 552 //<--DP        p->SetPhoton(kTRUE);\r
 553         fhRecPhoton->Fill(p->Pt()); //Reconstructed with primary photon\r
 554         for(Int_t iPID=0; iPID<4; iPID++){\r
 555           if(p->IsPIDOK(iPID,22))\r
 556             fhRecPhotonPID[iPID]->Fill(p->Pt());\r
 557         }\r
 558         Int_t pi0i=prim->GetFirstMother();\r
 559         Int_t grandpaPDG=-1 ;\r
 560         TParticle * pi0p = 0;\r
 561         if(pi0i>=0){\r
 562           pi0p=fStack->Particle(pi0i);\r
 563           grandpaPDG=pi0p->GetPdgCode() ;\r
 564         }\r
 565         switch(grandpaPDG){\r
 566         case 111: //Pi0 decay\r
 567                 //Primary decay photon (as in MC)\r
 568                 fhRecPhotPi0->Fill(p->Pt());\r
 569                 break ;\r
 570         case  11:\r
 571         case -11: //electron/positron conversion\r
 572 //<--DP                p->SetConverted(1);\r
 573                 fhRecPhotConv->Fill(p->Pt());  //Reconstructed with photon from conversion primary\r
 574                 fhConversionRadius->Fill(prim->R());\r
 575                 break ;\r
 576         case -2212:\r
 577         case -2112: //antineutron & antiproton conversion\r
 578                 fhRecPhotHadron->Fill(p->Pt());  //Reconstructed with photon from antibaryon annihilation\r
 579                 fhInteractionRadius->Fill(prim->R());\r
 580                 break ;\r
 581           \r
 582         case 221: //eta decay\r
 583                 fhRecPhotEta->Fill(p->Pt());\r
 584                 break ;  \r
 585           \r
 586         case 223: //omega meson decay\r
 587                 fhRecPhotOmega->Fill(p->Pt());\r
 588                 break ;\r
 589             \r
 590         case 331: //eta' decay\r
 591                 fhRecPhotEtapr->Fill(p->Pt());\r
 592                 break ;\r
 593               \r
 594         case -1: //direct photon or no primary\r
 595                 fhRecPhotDirect->Fill(p->Pt());\r
 596                 break ;\r
 597               \r
 598         default:  \r
 599                 fhRecPhotOther->Fill(p->Pt());\r
 600                 break ;\r
 601         }  \r
 602 \r
 603         //Now classify pi0 decay photon\r
 604         if(grandpaPDG==111){\r
 605 //<--DP          p->Pi0Decay(kTRUE); //Mark this photon as primary decayed\r
 606 //<--DP          p->Pi0Id(pi0i); //remember id of the parent pi0\r
 607 \r
 608           //Now check if second (partner) photon from pi0 decay hits PHOS or not\r
 609           //i.e. both photons can be tagged or it's the systematic error\r
 610 //<--DP          p->SetPartnerPt(0.); \r
 611           Int_t indexdecay1,indexdecay2;\r
 612 \r
 613           indexdecay1=pi0p->GetFirstDaughter();\r
 614           indexdecay2=pi0p->GetLastDaughter();\r
 615           Int_t indexdecay=-1;\r
 616           if(fDebug>2)\r
 617                 printf("checking pi0 decay...index1=%d, index2=%d, index_pi0=%d, index_ph_prim=%d\n", indexdecay1,indexdecay2,pi0i,caloCluster->GetLabel());\r
 618                 \r
 619           if(indexdecay1!=caloCluster->GetLabel()) \r
 620             indexdecay=indexdecay1;\r
 621           if(indexdecay2!=caloCluster->GetLabel()) \r
 622             indexdecay=indexdecay2;\r
 623           if(indexdecay==-1){\r
 624             if(fDebug>2){\r
 625               printf("Probably the other photon is not in the stack!\n");\r
 626               printf("Number of daughters: %d\n",pi0p->GetNDaughters());\r
 627             }\r
 628             fhDecWMissedPartnerStack->Fill(p->Pt()) ;
 629           }  \r
 630           else{\r
 631             TParticle *partner = fStack->Particle(indexdecay);\r
 632 //<--DP            p->SetPartnerPt(partner->Pt());\r
 633             if(partner->GetPdgCode()==22){ \r
 634               Bool_t isPartnerLost=kFALSE; //If partner is lost for some reason\r
 635               if(partner->GetNDaughters()!=0){ //this photon is converted before it is registered by some detector\r
 636                 if(fDebug>2)\r
 637                   printf("P_Conv, daughters=%d\n",partner->GetNDaughters());\r
 638 //<--DP                p->SetConvertedPartner(1);\r
 639                 fhPartnerMissedConv->Fill(partner->Pt());\r
 640                 fhDecWMissedPartnerConv->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
 641                 isPartnerLost=kTRUE;\r
 642               }\r
 643               Bool_t impact = kFALSE ;
 644               if(fPHOS){
 645                 Int_t modulenum ;
 646                 Double_t ztmp,xtmp ;
 647                 impact=fPHOSgeom->ImpactOnEmc(partner,modulenum,ztmp,xtmp) ;
 648                 if(fDebug>2){\r
 649                   printf("Impact on PHOS: module: %d, x tower: %f, z tower: %f\n", modulenum,xtmp,ztmp);\r
 650                 }
 651               }
 652               else{
 653                 impact = fEMCALgeom->Impact(partner) ;
 654               }
 655               if(!impact){ //this photon cannot hit PHOS\r
 656                 if(fDebug>2)\r
 657                   printf("P_Geo\n");\r
 658                 fhPartnerMissedGeo->Fill(partner->Pt());\r
 659                 fhDecWMissedPartnerGeom0->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
 660                 if(iFidArea>0){\r
 661                   fhDecWMissedPartnerGeom1->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
 662                   if(iFidArea>1){\r
 663                     fhDecWMissedPartnerGeom2->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
 664                     if(iFidArea>2){\r
 665                       fhDecWMissedPartnerGeom3->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
 666                     }\r
 667                   }\r
 668                 }\r
 669                 isPartnerLost=kTRUE;\r
 670               }\r
 671               if(!isPartnerLost && partner->Energy()<fMinEnergyCut){ //energy is not enough to be registered by PHOS\r
 672                 if(fDebug>2)\r
 673                   printf("P_Reg, E=%f\n",partner->Energy());\r
 674                 fhPartnerMissedEmin->Fill(partner->Pt());  //Spectrum of missed partners\r
 675                 fhDecWMissedPartnerEmin->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
 676                 isPartnerLost=kTRUE;\r
 677               }\r
 678               if(!isPartnerLost){\r
 679 //                p->SetMCTagged(1); //set this photon as primary tagged\r
 680                 fhDecWMCPartner->Fill(p->Pt());\r
 681                 fhPartnerMCReg->Fill(partner->Pt());\r
 682                 if(fDebug>2){\r
 683                   printf("both photons are inside PHOS. Energy: %f, Pt of pair photon: %f, E of pair photon: %f, Px: %f Py: %f Pz: %f, num of daughters: %d \n", caloCluster->E(),partner->Pt(),partner->Energy(),partner->Px(),partner->Py(),partner->Pz(),partner->GetNDaughters());\r
 684                 }\r
 685               }\r
 686               else{\r
 687                 fhDecWMissedPartnerAll->Fill(p->Pt());\r
 688               }\r
 689             }//Partner - photon\r
 690             else{//partner not photon\r
 691               fhDecWMissedPartnerNotPhoton->Fill(p->Pt());                \r
 692             }\r
 693           }\r
 694         }\r
 695       }\r
 696     }\r
 697   } //PHOS/EMCAL clusters\r
 698     \r
 699   if(fDebug>1)   \r
 700     printf("number of clusters: %d\n",inList);\r
 701 \r
 702   //Invariant Mass analysis\r
 703   for(Int_t phosPhoton1 = 0 ; phosPhoton1 < inList-1 ; phosPhoton1++) {\r
 704     AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton1));        \r
 705     for(Int_t phosPhoton2 = phosPhoton1 + 1 ; phosPhoton2 < inList ; phosPhoton2++) {\r
 706       AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton2));\r
 707 \r
 708       Double_t invMass = p1->GetPairMass(p2);\r
 709       fhInvMassReal->Fill(invMass,p1->Pt());\r
 710       fhInvMassReal->Fill(invMass,p2->Pt());\r
 711       if(fDebug>2)\r
 712           printf("Pair i=%d,j=%d, M=%f\n",phosPhoton1,phosPhoton2,invMass);\r
 713 \r
 714       Bool_t makePi01=IsInPi0Band(invMass,p1->Pt());\r
 715       Bool_t makePi02=IsInPi0Band(invMass,p2->Pt());\r
 716 \r
 717 \r
 718       if(makePi01 && p1->IsTagged()){//Multiple tagging\r
 719         fhTaggedMult->Fill(p1->Pt());\r
 720       }  \r
 721       if(makePi01 && !p1->IsTagged()){//Each photon should enter histogram once, even if tagged several times\r
 722         fhTaggedAll->Fill(p1->Pt());\r
 723         Int_t iFidArea = p1->GetFiducialArea(); \r
 724         if(iFidArea>0){\r
 725           fhTaggedArea1->Fill(p1->Pt() ) ; \r
 726           if(iFidArea>1){\r
 727             fhTaggedArea2->Fill(p1->Pt() ) ; \r
 728             if(iFidArea>2){\r
 729               fhTaggedArea3->Fill(p1->Pt() ) ; \r
 730             }\r
 731           }\r
 732         }\r
 733 \r
 734         for(Int_t iPID=0; iPID<4; iPID++){\r
 735           if(p1->IsPIDOK(iPID,22))\r
 736             fhTaggedPID[iPID]->Fill(p1->Pt());\r
 737         }\r
 738 \r
 739         p1->SetTagged(kTRUE) ;\r
 740       }  \r
 741       if(makePi02 && p2->IsTagged()){//Multiple tagging\r
 742         fhTaggedMult->Fill(p2->Pt());\r
 743       }  \r
 744       if(makePi02 && !p2->IsTagged()){//How should be account for multiply tagged photons?\r
 745         fhTaggedAll->Fill(p2->Pt());\r
 746         p2->SetTagged(kTRUE) ;\r
 747       }\r
 748       \r
 749       //Now get use MC information\r
 750       //First chesk if this is true pi0 pair\r
 751       if(IsSamePi0(p1,p2)){ //Correctly tagged - from the same pi0\r
 752 //        p1->SetTrueTagged(1);\r
 753 //        p2->SetTrueTagged(1);\r
 754         if(makePi01)//Correctly tagged photons\r
 755           fhTaggedMCTrue->Fill(p1->Pt());\r
 756         else{ //Decay pair missed tagging      \r
 757           fhMCMissedTagging->Fill(p1->Pt());\r
 758           fhMCMissedTaggingMass->Fill(invMass,p1->Pt()) ;\r
 759           //Clussify why missed tagging (todo)\r
 760           //Converted\r
 761           //Partner not a photon\r
 762           //Tagged not a photon\r
 763           //Just wrong inv.mass          \r
 764         }  \r
 765         if(makePi02)\r
 766           fhTaggedMCTrue->Fill(p2->Pt());\r
 767         else{      \r
 768           fhMCMissedTagging->Fill(p2->Pt());\r
 769           fhMCMissedTaggingMass->Fill(invMass,p2->Pt()) ;\r
 770           //Clussify why missed tagging (todo)\r
 771           //Converted\r
 772           //Partner not a photon\r
 773           //Tagged not a photon\r
 774           //Just wrong inv.mass          \r
 775         }  \r
 776       }\r
 777       else{//Fake tagged - not from the same pi0\r
 778         if(makePi01)//Fake pair\r
 779           fhMCFakeTagged->Fill(p1->Pt());\r
 780         if(makePi02)\r
 781           fhMCFakeTagged->Fill(p2->Pt());\r
 782       }\r
 783     }\r
 784   }\r
 785 \r
 786   //Fill Mixed InvMass distributions:\r
 787   TIter nextEv(fEventList) ;\r
 788   while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){\r
 789     Int_t nPhotons2 = event2->GetEntriesFast() ;\r
 790     for(Int_t i=0; i < inList ; i++){\r
 791       AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(i)) ;\r
 792       for(Int_t j=0; j < nPhotons2 ; j++){\r
 793         AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(event2->At(j)) ;\r
 794         Double_t invMass = p1->GetPairMass(p2);\r
 795         fhInvMassMixed->Fill(invMass,p1->Pt());\r
 796         fhInvMassMixed->Fill(invMass,p2->Pt());\r
 797       }\r
 798     }\r
 799   }\r
 800 \r
 801   //Remove old events\r
 802   fEventList->AddFirst(fCaloPhotonsArr);\r
 803   if(fEventList->GetSize() > 10){\r
 804     TClonesArray *tmp = static_cast <TClonesArray*> (fEventList->Last());\r
 805     fEventList->Remove(tmp);\r
 806     delete tmp;\r
 807   }\r
 808 \r
 809   PostData(1, fOutputList);\r
 810 }\r
 811 \r
 812 \r
 813 //______________________________________________________________________________\r
 814 void AliAnalysisTaskTaggedPhotons::Init()\r
 815 {\r
 816   // Intialisation of parameters\r
 817   AliInfo("Doing initialisation") ; \r
 818   SetPhotonId(0.9) ; \r
 819   SetMinEnergyCut(0.4);\r
 820   SetPi0MeanParameters(0.136,0.,0.0,0.0);\r
 821 //  SetPi0MeanParameters(0.1377,-0.002566,0.001216,-0.0001256);\r
 822   SetPi0SigmaParameters(0.004508,0.005497,0.00000006);\r
 823 }\r
 824 \r
 825 //______________________________________________________________________________\r
 826 void AliAnalysisTaskTaggedPhotons::Terminate(Option_t *)\r
 827 {\r
 828   // Processing when the event loop is ended\r
 829 \r
 830   //Write everything to the file\r
 831   char outname[55];\r
 832   if(fPHOS)\r
 833     sprintf(outname,"Tagging_PHOS.root") ;\r
 834   else  \r
 835     sprintf(outname,"Tagging_EMCAL.root") ;\r
 836   TFile *outfile = new TFile (outname,"recreate");\r
 837 \r
 838 fhRecAll->Write();\r
 839 fhRecAllArea1->Write();\r
 840 fhRecAllArea2->Write();\r
 841 fhRecAllArea3->Write();\r
 842 fhRecPhoton->Write();\r
 843 fhRecOther->Write();\r
 844 fhRecPhotonPID[0]->Write();\r
 845 fhRecPhotonPID[1]->Write();\r
 846 fhRecPhotonPID[2]->Write();\r
 847 fhRecPhotonPID[3]->Write();\r
 848 fhRecOtherPID[0]->Write();\r
 849 fhRecOtherPID[1]->Write();\r
 850 fhRecOtherPID[2]->Write();\r
 851 fhRecOtherPID[3]->Write();\r
 852 fhRecPhotPi0->Write();\r
 853 fhRecPhotEta->Write();\r
 854 fhRecPhotOmega->Write();\r
 855 fhRecPhotEtapr->Write();\r
 856 fhRecPhotConv->Write();\r
 857 fhRecPhotHadron->Write();\r
 858 fhRecPhotDirect->Write();\r
 859 fhRecPhotOther->Write();\r
 860 fhDecWMCPartner->Write();\r
 861 fhDecWMissedPartnerNotPhoton->Write();\r
 862 fhDecWMissedPartnerAll->Write();\r
 863 fhDecWMissedPartnerEmin->Write();\r
 864 fhDecWMissedPartnerConv->Write();\r
 865 fhDecWMissedPartnerStack->Write();\r
 866 fhDecWMissedPartnerGeom0->Write();\r
 867 fhDecWMissedPartnerGeom1->Write();\r
 868 fhDecWMissedPartnerGeom2->Write();\r
 869 fhDecWMissedPartnerGeom3->Write();\r
 870 fhPartnerMCReg->Write();\r
 871 fhPartnerMissedEmin->Write();\r
 872 fhPartnerMissedConv->Write();\r
 873 fhPartnerMissedGeo->Write();\r
 874 fhTaggedAll->Write();\r
 875 fhTaggedArea1->Write();\r
 876 fhTaggedArea2->Write();\r
 877 fhTaggedArea3->Write();\r
 878 \r
 879 fhTaggedPID[0]->Write();\r
 880 fhTaggedPID[1]->Write();\r
 881 fhTaggedPID[2]->Write();\r
 882 fhTaggedPID[3]->Write();\r
 883 fhTaggedMult->Write();\r
 884 fhTaggedMCTrue->Write();\r
 885 fhMCMissedTagging->Write();\r
 886 fhMCFakeTagged->Write();\r
 887 fhInvMassReal->Write();\r
 888 fhInvMassMixed->Write();\r
 889 fhMCMissedTaggingMass->Write();\r
 890 fhConversionRadius->Write();\r
 891 fhInteractionRadius->Write();\r
 892 fhEvents->Write();\r
 893 \r
 894 /*\r
 895   fhAllPhotons->Write() ;\r
 896   fhNotPhotons->Write() ;\r
 897   fhAllPhotonsPrimary->Write() ;\r
 898   fhNotPhotonsPrimary->Write() ;\r
 899   fhfakeNotPhotons->Write() ;\r
 900 \r
 901   fhTaggedPhotons->Write();\r
 902   fhfakeTaggedPhotons->Write();\r
 903   fhDecayNotTaggedPhotons->Write();\r
 904   fhstrangeNotTaggedPhotons->Write();\r
 905   fhstrangeNotTaggedPhotonsPair->Write();\r
 906   fhstrangeNotTaggedPhotonsRegCut->Write();\r
 907   fhstrangeNotTaggedPhotonsPairRegCut->Write();\r
 908 \r
 909   fhPi0DecayPhotonsPrimary->Write();\r
 910   fhEtaDecayPhotonsPrimary->Write();\r
 911   fhOmegaDecayPhotonsPrimary->Write();\r
 912   fhEtaSDecayPhotonsPrimary->Write();\r
 913   fhOtherDecayPhotonsPrimary->Write();\r
 914   fhDecayPhotonsPrimary->Write();\r
 915   fhConvertedPhotonsPrimary->Write();\r
 916   fhConvertedPhotonsPrimaryHadronsDecays->Write();\r
 917   fhCoordsConvertion->Write();\r
 918   fhCoordsConvertion2->Write();\r
 919 \r
 920   fhPHOSInvariantMassReal->Write();\r
 921   fhPHOSInvariantMassMixed->Write();\r
 922 \r
 923   fhPHOSPi0->Write();\r
 924 \r
 925   fhPi0DecayPhotonsGeomfake->Write();\r
 926   fhPi0DecayPhotonsTaggedPrimary->Write();\r
 927   fhPi0DecayPhotonsTaggedPrimaryPair->Write();\r
 928   fhPi0DecayPhotonsTaggedPrimaryRegCut->Write();\r
 929   fhPi0DecayPhotonsTaggedPrimaryPairRegCut->Write();\r
 930 \r
 931   fhPi0DecayPhotonsBigDecay->Write();\r
 932   fhPi0DecayPhotonsPConv->Write();\r
 933   fhPi0DecayPhotonsPGeo->Write();\r
 934   fhPi0DecayPhotonsPReg->Write();\r
 935 \r
 936   fhfakeTaggedPhotonsConv->Write();\r
 937   fhfakeTaggedPhotonsPID->Write();\r
 938 \r
 939   fhTrackRefCoords->Write();\r
 940   fhEvents->Write();\r
 941 */\r
 942 \r
 943 outfile->Close();\r
 944 \r
 945 }\r
 946 //______________________________________________________________________________\r
 947 Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)const\r
 948 {\r
 949   //Parameterization of the pi0 peak region\r
 950   Double_t mpi0mean = fPi0MeanP0 + fPi0MeanP1 * pt + fPi0MeanP2 * pt*pt + fPi0MeanP3 * pt*pt*pt;\r
 951   Double_t mpi0sigma = TMath::Sqrt(fPi0SigmaP0 * fPi0SigmaP0 / pt + fPi0SigmaP1 * fPi0SigmaP1 + fPi0SigmaP2 * fPi0SigmaP2 / pt / pt);\r
 952  \r
 953   return (m>mpi0mean-2*mpi0sigma && m<mpi0mean+2*mpi0sigma) ;\r
 954 }\r
 955 //______________________________________________________________________________\r
 956 Bool_t AliAnalysisTaskTaggedPhotons::IsSamePi0(const AliAODPWG4Particle *p1, const AliAODPWG4Particle *p2)const{\r
 957   //Looks through parents and finds if there was commont pi0 among ancestors\r
 958 \r
 959   if(!fStack)\r
 960     return kFALSE ; //can not say anything\r
 961 \r
 962   Int_t prim1 = p1->GetLabel();\r
 963   while(prim1!=-1){ \r
 964     Int_t prim2 = p2->GetLabel();\r
 965     while(prim2!=-1){ \r
 966       if(prim1==prim2){\r
 967         if(fStack->Particle(prim1)->GetPdgCode()==111)\r
 968           return kTRUE ;\r
 969         else\r
 970           return kFALSE ;\r
 971       }\r
 972       prim2=fStack->Particle(prim2)->GetFirstMother() ;\r
 973     }\r
 974     prim1=fStack->Particle(prim1)->GetFirstMother() ;\r
 975   }\r
 976   return kFALSE ;\r
 977 }\r
 978 //______________________________________________________________________________\r
 979 Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(Float_t * pos)const{\r
 980   //calculates in which kind of fiducial area photon hit\r
 981   Double_t phi=TMath::ATan2(pos[1],pos[0]) ;\r
 982   Double_t z=pos[2] ;
 983   while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;\r
 984   while(phi<0.)phi+=TMath::TwoPi() ;\r
 985   if(fPHOS){\r
 986     //From active PHOS area remove bands in 10 cm\r
 987     const Double_t kphi=TMath::ATan(10./460.) ; //angular band width\r
 988     Double_t dzMax=TMath::Ceil((fZmax-z)/10.) ;\r
 989     Double_t dzMin=TMath::Ceil((z-fZmin)/10.) ;\r
 990     Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);\r
 991     Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);\r
 992     return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin)); \r
 993   }\r
 994   else{//EMCAL\r
 995     //From active EMCAL area remove bands in 20 cm\r
 996     const Double_t kphi=TMath::ATan(20./428.) ; //angular band width\r
 997     Double_t dzMax=TMath::Ceil((fZmax-z)/20.) ;\r
 998     Double_t dzMin=TMath::Ceil((z-fZmin)/20.) ;\r
 999     Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);\r
1000     Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);\r
1001     return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin)); \r
1002   }\r
1003 }\r
1004 //______________________________________________________________________________\r
1005 Bool_t  AliAnalysisTaskTaggedPhotons::TestDisp(Double_t l0, Double_t l1, Double_t e)const{
1006   //test if dispersion corresponds to those of photon
1007   if(fPHOS){
1008     Double_t l0mean=1.38736+0.490405*TMath::Exp(-e*0.286170) ;
1009     Double_t l1mean=1.09786-0.323469*TMath::Exp(-e*0.918719) ;
1010     Double_t l0sigma=0.159905+0.829831/e-0.158067/e/e ;
1011     Double_t l1sigma=0.133170+0.404387/e-0.0426302/e/e ;
1012     Double_t c =-0.382233 ;
1013     return ((l0-l0mean)*(l0-l0mean)/l0sigma/l0sigma + (l1-l1mean)*(l1-l1mean)/l1sigma/l1sigma+c*(l0-l0mean)*(l1-l1mean)/l0sigma/l1sigma)<1. ;
1014   }
1015   else{ //EMCAL: not ready yet
1016    return kTRUE ;
1017
1018   }
1019
1020 }
1021