This version is bug fixed to work on analysis train (D. Blau)
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 5 Mar 2010 09:34:30 +0000 (09:34 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 5 Mar 2010 09:34:30 +0000 (09:34 +0000)
PWG4/PartCorrDep/AliAnaCalorimeterQA.cxx
PWG4/PartCorrDep/AliAnalysisTaskTaggedPhotons.cxx
PWG4/PartCorrDep/AliAnalysisTaskTaggedPhotons.h

index 72eab2a..d409941 100755 (executable)
@@ -408,7 +408,8 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
        fhEtaPhiE->SetYTitle("#phi (rad)");
        fhEtaPhiE->SetZTitle("E (GeV) ");
        outputContainer->Add(fhEtaPhiE);
-       
+
+
        //Track Matching
 
        fhECharged  = new TH1F ("hECharged","E reconstructed clusters, matched with track", nptbins,ptmin,ptmax); 
@@ -434,6 +435,7 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
        fhEtaPhiECharged->SetZTitle("E (GeV) ");
        outputContainer->Add(fhEtaPhiECharged); 
 
+
        fhEChargedNoOut  = new TH1F ("hEChargedNoOut","E reconstructed clusters, matched with track, no output track params", nptbins,ptmin,ptmax); 
        fhEChargedNoOut->SetXTitle("E (GeV)");
        outputContainer->Add(fhEChargedNoOut);
@@ -500,6 +502,7 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
        fhNCellsPerCluster->SetZTitle("#eta");
        outputContainer->Add(fhNCellsPerCluster);
        
+       
        fhNCellsPerClusterMIP  = new TH3F ("hNCellsPerClusterMIP","# cells per cluster vs energy vs #eta, smaller bin for MIP search", 
                                                                                100,0.,1., 6,0,5,netabins,etamin,etamax); 
        fhNCellsPerClusterMIP->SetXTitle("E (GeV)");
@@ -507,6 +510,7 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
        fhNCellsPerClusterMIP->SetZTitle("#eta");
        outputContainer->Add(fhNCellsPerClusterMIP);
        
+       
        fhNCellsPerClusterMIPCharged  = new TH3F ("hNCellsPerClusterMIPCharged","# cells per track-matched cluster vs energy vs #eta, smaller bin for MIP search", 
                                                                                                100,0.,1., 6,0,5,netabins,etamin,etamax); 
        fhNCellsPerClusterMIPCharged->SetXTitle("E (GeV)");
@@ -514,6 +518,7 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
        fhNCellsPerClusterMIPCharged->SetZTitle("#eta");
        outputContainer->Add(fhNCellsPerClusterMIPCharged);
        
+       
        fhNClusters  = new TH1F ("hNClusters","# clusters", nbins,nmin,nmax); 
        fhNClusters->SetXTitle("number of clusters");
        outputContainer->Add(fhNClusters);
@@ -563,7 +568,6 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
        fhXYZ->SetYTitle("y (cm)");
        fhXYZ->SetZTitle("z (cm) ");
        outputContainer->Add(fhXYZ);
-       
                
        fhRCellE  = new TH2F ("hRCellE","Cell R position vs cell energy",rbins,rmin,rmax,nptbins,ptmin,ptmax); 
        fhRCellE->SetXTitle("r = #sqrt{x^{2}+y^{2}} (cm)");
@@ -588,9 +592,10 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
        fhXYZCell  = new TH3F ("hXYZCell","Cell : x vs y vs z",xbins,xmin,xmax,ybins,ymin,ymax,zbins,zmin,zmax); 
        fhXYZCell->SetXTitle("x (cm)");
        fhXYZCell->SetYTitle("y (cm)");
-       fhXYZCell->SetZTitle("z (cm) ");
+       fhXYZCell->SetZTitle("z (cm)");
        outputContainer->Add(fhXYZCell);
-       
+
+
        Float_t dx = TMath::Abs(xmin)+TMath::Abs(xmax);
        Float_t dy = TMath::Abs(ymin)+TMath::Abs(ymax);
        Float_t dz = TMath::Abs(zmin)+TMath::Abs(zmax);
@@ -641,7 +646,8 @@ TList *  AliAnaCalorimeterQA::GetCreateOutputObjects()
        fhEtaPhiAmp->SetYTitle("#phi (rad)");
        fhEtaPhiAmp->SetZTitle("E (GeV) ");
        outputContainer->Add(fhEtaPhiAmp);              
-       
+
+
        //Calo cells
        fhNCells  = new TH1F ("hNCells","# cells", colmax*rowmax*fNModules,0,colmax*rowmax*fNModules); 
        fhNCells->SetXTitle("n cells");
@@ -1435,7 +1441,6 @@ void AliAnaCalorimeterQA::Print(const Option_t * opt) const
 void  AliAnaCalorimeterQA::MakeAnalysisFillHistograms() 
 {
        //Fill Calorimeter QA histograms
-       
        TLorentzVector mom ;
        TLorentzVector mom2 ;
        TRefArray * caloClusters = new TRefArray();
index ba5ebae..ac12481 100644 (file)
@@ -1,60 +1,60 @@
-/**************************************************************************\r
- * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- *                                                                        *\r
- * Author: The ALICE Off-line Project.                                    *\r
- * Contributors are mentioned in the code where appropriate.              *\r
- *                                                                        *\r
- * Permission to use, copy, modify and distribute this software and its   *\r
- * documentation strictly for non-commercial purposes is hereby granted   *\r
- * without fee, provided that the above copyright notice appears in all   *\r
- * copies and that both the copyright notice and this permission notice   *\r
- * appear in the supporting documentation. The authors make no claims     *\r
- * about the suitability of this software for any purpose. It is          *\r
- * provided "as is" without express or implied warranty.                  *\r
- **************************************************************************/\r
-\r
-/* $Id: */\r
-\r
-//_________________________________________________________________________\r
-// Analysis for Tagged Photons\r
-// Prepares all necessary histograms for calculation of \r
-// the yield of pi0 decay photon in calorimeter:\r
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/* $Id: */
+
+//_________________________________________________________________________
+// Analysis for Tagged Photons
+// Prepares all necessary histograms for calculation of 
+// the yield of pi0 decay photon in calorimeter:
 // Marks photons which makes pi0 with some other and
 // fill invariant mass distributions for estimate background below pi0
 // peak so that estimate proportion of fake pairs. 
 // Fills as well controll MC histograms with clasification of the photon origin 
 // and check of the ptoportion of truly tagged photons.
-// \r
-//\r
-//*-- Dmitry Blau \r
-//////////////////////////////////////////////////////////////////////////////\r
-\r
-#include <TH1.h>\r
-#include <TH2.h>\r
-#include <TH3.h>\r
-#include <TFile.h>\r
-#include <TROOT.h>\r
-\r
-#include "AliAnalysisTaskTaggedPhotons.h" \r
-#include "AliAnalysisManager.h"\r
-#include "AliESDVertex.h"\r
-#include "AliESDEvent.h" \r
-#include "AliAODEvent.h" \r
-#include "AliESDCaloCluster.h" \r
-#include "AliAODPWG4Particle.h"\r
-#include "AliAnalysisManager.h"\r
-#include "AliLog.h"\r
-#include "TGeoManager.h"\r
+// 
+//
+//*-- Dmitry Blau 
+//////////////////////////////////////////////////////////////////////////////
+
+#include <TH1.h>
+#include <TH2.h>
+#include <TH3.h>
+#include <TFile.h>
+#include <TROOT.h>
+
+#include "AliAnalysisTaskTaggedPhotons.h" 
+#include "AliAnalysisManager.h"
+#include "AliESDVertex.h"
+#include "AliESDEvent.h" 
+#include "AliAODEvent.h" 
+#include "AliESDCaloCluster.h" 
+#include "AliAODPWG4Particle.h"
+#include "AliAnalysisManager.h"
+#include "AliLog.h"
+#include "TGeoManager.h"
 #include "AliMCAnalysisUtils.h"
-#include "AliMCEventHandler.h"\r
-#include "AliMCEvent.h"\r
-#include "AliStack.h"\r
-#include "AliPHOSGeoUtils.h"\r
-#include "AliEMCALGeoUtils.h"\r
+#include "AliMCEventHandler.h"
+#include "AliMCEvent.h"
+#include "AliStack.h"
+#include "AliPHOSGeoUtils.h"
+#include "AliEMCALGeoUtils.h"
+
 
 
-\r
-//______________________________________________________________________________\r
+//______________________________________________________________________________
 AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons() : AliAnalysisTaskSE(),
   fPHOSgeom(0x0),
   fEMCALgeom(0x0),
@@ -77,12 +77,12 @@ AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons() : AliAnalysisTaskSE
   fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
   fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
   fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
-{\r
-  //Default constructor\r
-}\r
-//______________________________________________________________________________\r
-AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const char *name) : \r
-  AliAnalysisTaskSE(name),\r
+{
+  //Default constructor
+}
+//______________________________________________________________________________
+AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const char *name) : 
+  AliAnalysisTaskSE(name),
   fPHOSgeom(0x0),
   fEMCALgeom(0x0),
   fStack(0x0),fPHOS(1),
@@ -104,16 +104,16 @@ AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const char *name) :
   fhTaggedMCTrue(0x0),fhMCMissedTagging(0x0),fhMCFakeTagged(0x0),
   fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
   fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
-{\r
-  // Constructor.\r
-\r
-  // Output slots \r
-  DefineOutput(1,  TList::Class()) ; \r
-}\r
-\r
-//____________________________________________________________________________\r
-AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) :\r
-  AliAnalysisTaskSE(ap.GetName()),  \r
+{
+  // Constructor.
+
+  // Output slots 
+  DefineOutput(1,  TList::Class()) ; 
+}
+
+//____________________________________________________________________________
+AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) :
+  AliAnalysisTaskSE(ap.GetName()),  
   fPHOSgeom(0x0),
   fEMCALgeom(0x0),
   fStack(0x0),fPHOS(ap.fPHOS),
@@ -137,40 +137,40 @@ AliAnalysisTaskTaggedPhotons::AliAnalysisTaskTaggedPhotons(const AliAnalysisTask
   fhInvMassReal(0x0),fhInvMassMixed(0x0),fhMCMissedTaggingMass(0x0),
   fhConversionRadius(0x0),fhInteractionRadius(0x0),fhEvents(0x0)
 {
-  // cpy ctor\r
-}\r
-\r
-//_____________________________________________________________________________\r
-AliAnalysisTaskTaggedPhotons& AliAnalysisTaskTaggedPhotons::operator = (const AliAnalysisTaskTaggedPhotons& ap)\r
-{\r
-// assignment operator\r
-\r
-  this->~AliAnalysisTaskTaggedPhotons();\r
-  new(this) AliAnalysisTaskTaggedPhotons(ap);\r
-  return *this;\r
-}\r
-\r
-//______________________________________________________________________________\r
-AliAnalysisTaskTaggedPhotons::~AliAnalysisTaskTaggedPhotons()\r
-{\r
-  // dtor\r
-  if(fOutputList) {\r
-    fOutputList->Clear() ; \r
-    delete fOutputList ;\r
-  }\r
-}\r
-\r
-\r
-//________________________________________________________________________\r
-void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()\r
-{ \r
-\r
-
-  //Load geometry\r
-  //if file "geometry.root" exists, load it\r
-  //otherwise use misaligned matrixes stored in ESD\r
-  TFile *geoFile = new TFile("geometry.root","read");\r
-  if(geoFile->IsZombie()){ //no file, load geo matrixes from ESD\r
+  // cpy ctor
+}
+
+//_____________________________________________________________________________
+AliAnalysisTaskTaggedPhotons& AliAnalysisTaskTaggedPhotons::operator = (const AliAnalysisTaskTaggedPhotons& ap)
+{
+// assignment operator
+
+  this->~AliAnalysisTaskTaggedPhotons();
+  new(this) AliAnalysisTaskTaggedPhotons(ap);
+  return *this;
+}
+
+//______________________________________________________________________________
+AliAnalysisTaskTaggedPhotons::~AliAnalysisTaskTaggedPhotons()
+{
+  // dtor
+  if(fOutputList) {
+    fOutputList->Clear() ; 
+    delete fOutputList ;
+  }
+}
+
+
+//________________________________________________________________________
+void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
+{ 
+
+
+  //Load geometry
+  //if file "geometry.root" exists, load it
+  //otherwise use misaligned matrixes stored in ESD
+  TFile *geoFile = new TFile("geometry.root","read");
+  if(geoFile->IsZombie()){ //no file, load geo matrixes from ESD
     AliInfo("Can not find file geometry.root, reading misalignment matrixes from ESD/AOD") ;
     AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent()) ;
     AliAODEvent * aod = 0x0 ;
@@ -179,7 +179,7 @@ void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
     if(!esd && !aod)
       AliFatal("Can not read geometry even from ESD/AOD.") ;
     if(fPHOS){//reading PHOS matrixes
-      fPHOSgeom = new AliPHOSGeoUtils("IHEP","");\r
+      fPHOSgeom = new AliPHOSGeoUtils("IHEP","");
       for(Int_t mod=0; mod<5; mod++){
         if(esd){
           const TGeoHMatrix* m=esd->GetPHOSMatrix(mod) ;
@@ -204,26 +204,26 @@ void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
         }
       }
     }
-  }\r
-  else{\r
-    gGeoManager = (TGeoManager*)geoFile->Get("Geometry");\r
+  }
+  else{
+    gGeoManager = (TGeoManager*)geoFile->Get("Geometry");
     //Geometry will be misaligned from GeoManager
     if(fPHOS){
-      fPHOSgeom = new AliPHOSGeoUtils("IHEP","");\r
+      fPHOSgeom = new AliPHOSGeoUtils("IHEP","");
     }
     else{
-      fEMCALgeom = new AliEMCALGeoUtils("EMCAL_COMPLETE");\r
+      fEMCALgeom = new AliEMCALGeoUtils("EMCAL_COMPLETE");
     }
-  }\r
-\r
-\r
-  if(fPHOSgeom==NULL && fEMCALgeom==NULL){\r
-    AliError("Error loading Geometry\n");\r
-  }\r
-  else\r
-    AliInfo("Geometry loaded... OK\n");\r
-\r
-  //Evaluate active PHOS/EMCAL area\r
+  }
+
+
+  if(fPHOSgeom==NULL && fEMCALgeom==NULL){
+    AliError("Error loading Geometry\n");
+  }
+  else
+    AliInfo("Geometry loaded... OK\n");
+
+  //Evaluate active PHOS/EMCAL area
   if(fZmax<=fZmin){ //not set yet
     if(fPHOS){
       //Check active modules in the current configuration
@@ -257,10 +257,10 @@ void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
       }
       else{
         //Use approximate range
-        fZmax= 2.25*56/2. ;\r
-        fZmin=-2.25*56/2. ;\r
-        fPhimax=220./180.*TMath::Pi() ;\r
-        fPhimin=320./180.*TMath::Pi() ;\r
+        fZmax= 2.25*56/2. ;
+        fZmin=-2.25*56/2. ;
+        fPhimax=220./180.*TMath::Pi() ;
+        fPhimin=320./180.*TMath::Pi() ;
       }
     }
     else{ //Similar for EMCAL <--Gustavo, Could you please have a look?
@@ -300,220 +300,224 @@ void AliAnalysisTaskTaggedPhotons::UserCreateOutputObjects()
       }
     }
   }
-\r
-\r
-  // Create the outputs containers\r
-\r
-  OpenFile(1) ; \r
-  const Int_t nPtBins=52 ;\r
-  Double_t ptBins[nPtBins+1] ;\r
-  for(Int_t i=0;i<=20;i++)ptBins[i]=0.1*i ; //0-2 GeV:  0.1 GeV/bin\r
-  for(Int_t i=21;i<=30;i++)ptBins[i]=2.+0.2*(i-20) ; //2-4 GeV:  0.2 GeV/bin                   \r
-  for(Int_t i=31;i<=38;i++)ptBins[i]=4.+0.5*(i-30) ; //4-8 GeV:  0.5 GeV/bin                   \r
-  for(Int_t i=39;i<=42;i++)ptBins[i]=8.+1.0*(i-38) ; //8-12 GeV: 1. GeV/bin                   \r
-  for(Int_t i=43;i<=52;i++)ptBins[i]=12.+2.0*(i-42) ; //12-30 GeV: 2. GeV/bin                   \r
-\r
-\r
-  //Reconstructed spectra\r
-    fhRecAll      = new TH1D("fhRecAll","Spectrum of all reconstructed particles", nPtBins, ptBins ) ;\r
-    fhRecAllArea1 = new TH1D("fhRecAllArea1","Spectrum of rec particles in Fid. Area 1", nPtBins, ptBins ) ;\r
-    fhRecAllArea2 = new TH1D("fhRecAllArea2","Spectrum of rec particles in Fid. Area 2", nPtBins, ptBins ) ;\r
-    fhRecAllArea3 = new TH1D("fhRecAllArea3","Spectrum of rec particles in Fid. Area 3", nPtBins, ptBins ) ;\r
-\r
-    //Sort registered particles spectra according MC information\r
-    fhRecPhoton   = new TH1D("fhRecPhoton","Spectrum of rec. with primary==22 and no PID criteria", nPtBins, ptBins ) ;\r
-    fhRecOther    = new TH1D("fhRecOther"," Spectrum of rec. with primary!=22 and no PID criteria", nPtBins, ptBins ) ;\r
-    fhRecPhotonPID[0]= new TH1D("fhRecPhotonPID0","Spectrum of rec. with primary==22 and PID=0", nPtBins, ptBins ) ;\r
-    fhRecPhotonPID[1]= new TH1D("fhRecPhotonPID1","Spectrum of rec. with primary==22 and PID=1", nPtBins, ptBins ) ;\r
-    fhRecPhotonPID[2]= new TH1D("fhRecPhotonPID2","Spectrum of rec. with primary==22 and PID=2", nPtBins, ptBins ) ;\r
-    fhRecPhotonPID[3]= new TH1D("fhRecPhotonPID3","Spectrum of rec. with primary==22 and PID=3", nPtBins, ptBins ) ;\r
-    fhRecOtherPID[0]= new TH1D("fhRecOtherPID0","Spectrum of rec. with primary!=22 and PID=0", nPtBins, ptBins ) ;\r
-    fhRecOtherPID[1]= new TH1D("fhRecOtherPID1","Spectrum of rec. with primary!=22 and PID=1", nPtBins, ptBins ) ;\r
-    fhRecOtherPID[2]= new TH1D("fhRecOtherPID2","Spectrum of rec. with primary!=22 and PID=2", nPtBins, ptBins ) ;\r
-    fhRecOtherPID[3]= new TH1D("fhRecOtherPID3","Spectrum of rec. with primary!=22 and PID=3", nPtBins, ptBins ) ;\r
-    fhRecPhotPi0    = new TH1D("fhRecPhotPi0","Spectrum of rec. photons from pi0 decays", nPtBins, ptBins ) ;\r
-    fhRecPhotEta    = new TH1D("fhRecPhotEta","Spectrum of rec. photons from eta decays", nPtBins, ptBins ) ;\r
-    fhRecPhotOmega  = new TH1D("fhRecPhotOmega","Spectrum of rec. photons from omega decays", nPtBins, ptBins ) ;\r
-    fhRecPhotEtapr  = new TH1D("fhRecPhotEtapr","Spectrum of rec. photons from eta prime decays", nPtBins, ptBins ) ;\r
-    fhRecPhotConv   = new TH1D("fhRecPhotConv"," Spectrum of rec. photons from conversion", nPtBins, ptBins ) ;\r
-    fhRecPhotHadron = new TH1D("fhRecPhotHadron","Spectrum of rec. photons from hadron-matter interactions", nPtBins, ptBins ) ;\r
-    fhRecPhotDirect = new TH1D("fhRecPhotDirect","Spectrum of rec. photons direct or no primary", nPtBins, ptBins ) ;\r
-    fhRecPhotOther  = new TH1D("fhRecPhotOther","Spectrum of rec. photons from other hadron decays", nPtBins, ptBins ) ;\r
-\r
-    //MC tagging: reasons of partner loss etc.\r
-    fhDecWMCPartner        = new TH1D("fhDecWMCPartner","pi0 decay photon which partner should be registered according to MC", nPtBins, ptBins ) ;\r
-    fhDecWMissedPartnerNotPhoton = new TH1D("fhDecWMissedPartnerNotPhoton","Decay photon with missed non-photon partner", nPtBins, ptBins ) ;\r
-    fhDecWMissedPartnerAll  = new TH1D("fhDecWMissedPartnerAll","Decay photons with partner missed due to some reason", nPtBins, ptBins ) ;\r
-    fhDecWMissedPartnerEmin = new TH1D("fhDecWMissedPartnerEmin","Decay photons with partner missed due to low energy", nPtBins, ptBins ) ;\r
-    fhDecWMissedPartnerConv = new TH1D("fhDecWMissedPartnerConv","Decay photons with partner missed due to conversion", nPtBins, ptBins ) ;\r
-    fhDecWMissedPartnerStack= new TH1D("fhDecWMissedPartnerStack","Decay photons with partner not in Stack", nPtBins, ptBins ) ;\r
-    fhDecWMissedPartnerGeom0 = new TH1D("fhDecWMissedPartnerGeom0","Decay photons with partner missed due geometry", nPtBins, ptBins ) ;\r
-    fhDecWMissedPartnerGeom1 = new TH1D("fhDecWMissedPartnerGeom1","Decay photons with partner missed due geometry Fid. area. 1", nPtBins, ptBins ) ;\r
-    fhDecWMissedPartnerGeom2 = new TH1D("fhDecWMissedPartnerGeom2","Decay photons with partner missed due geometry Fid. area. 2", nPtBins, ptBins ) ;\r
-    fhDecWMissedPartnerGeom3 = new TH1D("fhDecWMissedPartnerGeom3","Decay photons with partner missed due geometry Fid. area. 3", nPtBins, ptBins ) ;\r
-\r
-    //MC tagging: Decay partners spectra\r
-    fhPartnerMCReg      = new TH1D("fhPartnerMCReg","Spectrum of decay partners which should be registered (MC)", nPtBins, ptBins ) ;\r
-    fhPartnerMissedEmin = new TH1D("fhPartnerMissedEmin","Spectrum of decay partners lost due to Emin cut", nPtBins, ptBins ) ;\r
-    fhPartnerMissedConv = new TH1D("fhPartnerMissedConv","Spectrum of decay partners lost due to conversion", nPtBins, ptBins ) ;\r
-    fhPartnerMissedGeo  = new TH1D("fhPartnerMissedGeo","Spectrum of decay partners lost due to acceptance", nPtBins, ptBins ) ;\r
-\r
-    //Tagging\r
-    fhTaggedAll    = new TH1D("fhTaggedAll","Spectrum of all tagged photons", nPtBins, ptBins ) ;\r
-    fhTaggedArea1  = new TH1D("fhTaggedArea1","Spectrum of all tagged photons Fid. area1", nPtBins, ptBins ) ;\r
-    fhTaggedArea2  = new TH1D("fhTaggedArea2","Spectrum of all tagged photons Fid. area2", nPtBins, ptBins ) ;\r
-    fhTaggedArea3  = new TH1D("fhTaggedArea3","Spectrum of all tagged photons Fid. area3", nPtBins, ptBins ) ;\r
-    fhTaggedPID[0] = new TH1D("fhTaggedPID0","Spectrum of tagged photons for PID=0", nPtBins, ptBins ) ;\r
-    fhTaggedPID[1] = new TH1D("fhTaggedPID1","Spectrum of tagged photons for PID=1", nPtBins, ptBins ) ;\r
-    fhTaggedPID[2] = new TH1D("fhTaggedPID2","Spectrum of tagged photons for PID=2", nPtBins, ptBins ) ;\r
-    fhTaggedPID[3] = new TH1D("fhTaggedPID3","Spectrum of tagged photons for PID=3", nPtBins, ptBins ) ;\r
-    fhTaggedMult   = new TH1D("fhTaggedMult","Spectrum of multiply tagged photons", nPtBins, ptBins ) ;\r
-\r
-    //Tagging: use MC information if available\r
-    fhTaggedMCTrue    = new TH1D("fhTaggedMCTrue","Spectrum of correctly tagged pi0 decay photons", nPtBins, ptBins ) ;\r
-    fhMCMissedTagging = new TH1D("fhMCMissedTagging","Spectrum of pi0 decay photons missed tagging due to wrong pair mass", nPtBins, ptBins ) ;\r
-    fhMCFakeTagged    = new TH1D("fhMCFakeTagged","Spectrum of photons wrongly tagged according to MC", nPtBins, ptBins ) ;\r
-\r
-    //Invariant mass distributions for fake corrections\r
-    const Int_t nmass=1000 ;\r
-    Double_t masses[nmass+1] ;\r
-    for(Int_t i=0;i<=nmass;i++)masses[i]=0.001*i ;\r
-    fhInvMassReal = new TH2D("fhInvMassReal","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;\r
-    fhInvMassMixed  = new TH2D("fhInvMassMixed","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;\r
-    fhMCMissedTaggingMass= new TH2D("fhMCMissedTaggingMass","Inv mass of pairs missed tagging",nmass,masses,nPtBins, ptBins ) ;\r
-\r
-    //Conversion and annihilation radius distributions\r
-    fhConversionRadius  = new TH1D("fhConversionRadius","Radis of photon production (conversion)",100,0.,500.) ;\r
-    fhInteractionRadius = new TH1D("fhInteractionRadius","Radis of photon production (hadron interaction)",100,0.,500.);\r
-\r
-    fhEvents               = new TH1D("hEvents", "Number of Events processed", 1, 0, 1);\r
-\r
-// create output container\r
-  \r
-  fOutputList = new TList() ;\r
-  fEventList = new TList() ; \r
-  fOutputList->SetName(GetName()) ; \r
-\r
-  fOutputList->AddAt(fhRecAll,                               0) ;\r
-  fOutputList->AddAt(fhRecAllArea1,                          1) ;\r
-  fOutputList->AddAt(fhRecAllArea2,                          2) ;\r
-  fOutputList->AddAt(fhRecAllArea3,                          3) ;\r
-\r
-  fOutputList->AddAt(fhRecPhoton,                            4) ;\r
-  fOutputList->AddAt(fhRecOther,                             5) ;\r
-  fOutputList->AddAt(fhRecPhotonPID[0],                      6) ;\r
-  fOutputList->AddAt(fhRecPhotonPID[1],                      7) ;\r
-  fOutputList->AddAt(fhRecPhotonPID[2],                      8) ;\r
-  fOutputList->AddAt(fhRecPhotonPID[3],                      9) ;\r
-  fOutputList->AddAt(fhRecOtherPID[0],                      10) ;\r
-  fOutputList->AddAt(fhRecOtherPID[1],                      11) ;\r
-  fOutputList->AddAt(fhRecOtherPID[2],                      12) ;\r
-  fOutputList->AddAt(fhRecOtherPID[3],                      13) ;\r
-  fOutputList->AddAt(fhRecPhotPi0,                          14) ;\r
-  fOutputList->AddAt(fhRecPhotEta,                          15) ;\r
-  fOutputList->AddAt(fhRecPhotOmega,                        16) ;\r
-  fOutputList->AddAt(fhRecPhotEtapr,                        17) ;\r
-  fOutputList->AddAt(fhRecPhotConv,                         18) ;\r
-  fOutputList->AddAt(fhRecPhotHadron,                       19) ;\r
-  fOutputList->AddAt(fhRecPhotDirect,                       20) ;\r
-  fOutputList->AddAt(fhRecPhotOther,                        21) ;\r
-\r
-  fOutputList->AddAt(fhDecWMCPartner,                       22) ;\r
-  fOutputList->AddAt(fhDecWMissedPartnerNotPhoton,          23) ;\r
-  fOutputList->AddAt(fhDecWMissedPartnerAll,                24) ;\r
-  fOutputList->AddAt(fhDecWMissedPartnerEmin,               25) ;\r
-  fOutputList->AddAt(fhDecWMissedPartnerConv,               26) ;\r
-  fOutputList->AddAt(fhDecWMissedPartnerStack,              27) ;\r
-  fOutputList->AddAt(fhDecWMissedPartnerGeom0,              28) ;\r
-  fOutputList->AddAt(fhDecWMissedPartnerGeom1,              29) ;\r
-  fOutputList->AddAt(fhDecWMissedPartnerGeom2,              30) ;\r
-  fOutputList->AddAt(fhDecWMissedPartnerGeom3,              31) ;\r
-\r
-  fOutputList->AddAt(fhPartnerMCReg,                        32) ;\r
-  fOutputList->AddAt(fhPartnerMissedEmin,                   33) ;\r
-  fOutputList->AddAt(fhPartnerMissedConv,                   34) ;\r
-  fOutputList->AddAt(fhPartnerMissedGeo,                    35) ;\r
-\r
-  fOutputList->AddAt(fhTaggedAll,                           36) ;\r
-  fOutputList->AddAt(fhTaggedArea1,                         37) ;\r
-  fOutputList->AddAt(fhTaggedArea2,                         38) ;\r
-  fOutputList->AddAt(fhTaggedArea3,                         39) ;\r
-  fOutputList->AddAt(fhTaggedPID[0],                        40) ;\r
-  fOutputList->AddAt(fhTaggedPID[1],                        41) ;\r
-  fOutputList->AddAt(fhTaggedPID[2],                        42) ;\r
-  fOutputList->AddAt(fhTaggedPID[3],                        43) ;\r
-  fOutputList->AddAt(fhTaggedMult,                          44) ;\r
-\r
-  fOutputList->AddAt(fhTaggedMCTrue,                        45) ;\r
-  fOutputList->AddAt(fhMCMissedTagging,                     46) ;\r
-  fOutputList->AddAt(fhMCFakeTagged,                        47) ;\r
-\r
-  fOutputList->AddAt(fhInvMassReal,                         48) ;\r
-  fOutputList->AddAt(fhInvMassMixed,                        49) ;\r
-  fOutputList->AddAt(fhMCMissedTaggingMass,                 50) ;\r
-\r
-  fOutputList->AddAt(fhConversionRadius,                    51) ;\r
-  fOutputList->AddAt(fhInteractionRadius,                   52) ;\r
-\r
-  fOutputList->AddAt(fhEvents,                              53) ;\r
-\r
-}\r
-\r
-//______________________________________________________________________________\r
-void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *) \r
-{\r
+
+
+  // Create the outputs containers
+
+  OpenFile(1) ; 
+  const Int_t nPtBins=52 ;
+  Double_t ptBins[nPtBins+1] ;
+  for(Int_t i=0;i<=20;i++)ptBins[i]=0.1*i ; //0-2 GeV:  0.1 GeV/bin
+  for(Int_t i=21;i<=30;i++)ptBins[i]=2.+0.2*(i-20) ; //2-4 GeV:  0.2 GeV/bin                   
+  for(Int_t i=31;i<=38;i++)ptBins[i]=4.+0.5*(i-30) ; //4-8 GeV:  0.5 GeV/bin                   
+  for(Int_t i=39;i<=42;i++)ptBins[i]=8.+1.0*(i-38) ; //8-12 GeV: 1. GeV/bin                   
+  for(Int_t i=43;i<=52;i++)ptBins[i]=12.+2.0*(i-42) ; //12-30 GeV: 2. GeV/bin                   
+
+
+  //Reconstructed spectra
+    fhRecAll      = new TH1D("fhRecAll","Spectrum of all reconstructed particles", nPtBins, ptBins ) ;
+    fhRecAllArea1 = new TH1D("fhRecAllArea1","Spectrum of rec particles in Fid. Area 1", nPtBins, ptBins ) ;
+    fhRecAllArea2 = new TH1D("fhRecAllArea2","Spectrum of rec particles in Fid. Area 2", nPtBins, ptBins ) ;
+    fhRecAllArea3 = new TH1D("fhRecAllArea3","Spectrum of rec particles in Fid. Area 3", nPtBins, ptBins ) ;
+
+    //Sort registered particles spectra according MC information
+    fhRecPhoton   = new TH1D("fhRecPhoton","Spectrum of rec. with primary==22 and no PID criteria", nPtBins, ptBins ) ;
+    fhRecOther    = new TH1D("fhRecOther"," Spectrum of rec. with primary!=22 and no PID criteria", nPtBins, ptBins ) ;
+    fhRecPhotonPID[0]= new TH1D("fhRecPhotonPID0","Spectrum of rec. with primary==22 and PID=0", nPtBins, ptBins ) ;
+    fhRecPhotonPID[1]= new TH1D("fhRecPhotonPID1","Spectrum of rec. with primary==22 and PID=1", nPtBins, ptBins ) ;
+    fhRecPhotonPID[2]= new TH1D("fhRecPhotonPID2","Spectrum of rec. with primary==22 and PID=2", nPtBins, ptBins ) ;
+    fhRecPhotonPID[3]= new TH1D("fhRecPhotonPID3","Spectrum of rec. with primary==22 and PID=3", nPtBins, ptBins ) ;
+    fhRecOtherPID[0]= new TH1D("fhRecOtherPID0","Spectrum of rec. with primary!=22 and PID=0", nPtBins, ptBins ) ;
+    fhRecOtherPID[1]= new TH1D("fhRecOtherPID1","Spectrum of rec. with primary!=22 and PID=1", nPtBins, ptBins ) ;
+    fhRecOtherPID[2]= new TH1D("fhRecOtherPID2","Spectrum of rec. with primary!=22 and PID=2", nPtBins, ptBins ) ;
+    fhRecOtherPID[3]= new TH1D("fhRecOtherPID3","Spectrum of rec. with primary!=22 and PID=3", nPtBins, ptBins ) ;
+    fhRecPhotPi0    = new TH1D("fhRecPhotPi0","Spectrum of rec. photons from pi0 decays", nPtBins, ptBins ) ;
+    fhRecPhotEta    = new TH1D("fhRecPhotEta","Spectrum of rec. photons from eta decays", nPtBins, ptBins ) ;
+    fhRecPhotOmega  = new TH1D("fhRecPhotOmega","Spectrum of rec. photons from omega decays", nPtBins, ptBins ) ;
+    fhRecPhotEtapr  = new TH1D("fhRecPhotEtapr","Spectrum of rec. photons from eta prime decays", nPtBins, ptBins ) ;
+    fhRecPhotConv   = new TH1D("fhRecPhotConv"," Spectrum of rec. photons from conversion", nPtBins, ptBins ) ;
+    fhRecPhotHadron = new TH1D("fhRecPhotHadron","Spectrum of rec. photons from hadron-matter interactions", nPtBins, ptBins ) ;
+    fhRecPhotDirect = new TH1D("fhRecPhotDirect","Spectrum of rec. photons direct or no primary", nPtBins, ptBins ) ;
+    fhRecPhotOther  = new TH1D("fhRecPhotOther","Spectrum of rec. photons from other hadron decays", nPtBins, ptBins ) ;
+
+    //MC tagging: reasons of partner loss etc.
+    fhDecWMCPartner        = new TH1D("fhDecWMCPartner","pi0 decay photon which partner should be registered according to MC", nPtBins, ptBins ) ;
+    fhDecWMissedPartnerNotPhoton = new TH1D("fhDecWMissedPartnerNotPhoton","Decay photon with missed non-photon partner", nPtBins, ptBins ) ;
+    fhDecWMissedPartnerAll  = new TH1D("fhDecWMissedPartnerAll","Decay photons with partner missed due to some reason", nPtBins, ptBins ) ;
+    fhDecWMissedPartnerEmin = new TH1D("fhDecWMissedPartnerEmin","Decay photons with partner missed due to low energy", nPtBins, ptBins ) ;
+    fhDecWMissedPartnerConv = new TH1D("fhDecWMissedPartnerConv","Decay photons with partner missed due to conversion", nPtBins, ptBins ) ;
+    fhDecWMissedPartnerStack= new TH1D("fhDecWMissedPartnerStack","Decay photons with partner not in Stack", nPtBins, ptBins ) ;
+    fhDecWMissedPartnerGeom0 = new TH1D("fhDecWMissedPartnerGeom0","Decay photons with partner missed due geometry", nPtBins, ptBins ) ;
+    fhDecWMissedPartnerGeom1 = new TH1D("fhDecWMissedPartnerGeom1","Decay photons with partner missed due geometry Fid. area. 1", nPtBins, ptBins ) ;
+    fhDecWMissedPartnerGeom2 = new TH1D("fhDecWMissedPartnerGeom2","Decay photons with partner missed due geometry Fid. area. 2", nPtBins, ptBins ) ;
+    fhDecWMissedPartnerGeom3 = new TH1D("fhDecWMissedPartnerGeom3","Decay photons with partner missed due geometry Fid. area. 3", nPtBins, ptBins ) ;
+
+    //MC tagging: Decay partners spectra
+    fhPartnerMCReg      = new TH1D("fhPartnerMCReg","Spectrum of decay partners which should be registered (MC)", nPtBins, ptBins ) ;
+    fhPartnerMissedEmin = new TH1D("fhPartnerMissedEmin","Spectrum of decay partners lost due to Emin cut", nPtBins, ptBins ) ;
+    fhPartnerMissedConv = new TH1D("fhPartnerMissedConv","Spectrum of decay partners lost due to conversion", nPtBins, ptBins ) ;
+    fhPartnerMissedGeo  = new TH1D("fhPartnerMissedGeo","Spectrum of decay partners lost due to acceptance", nPtBins, ptBins ) ;
+
+    //Tagging
+    fhTaggedAll    = new TH1D("fhTaggedAll","Spectrum of all tagged photons", nPtBins, ptBins ) ;
+    fhTaggedArea1  = new TH1D("fhTaggedArea1","Spectrum of all tagged photons Fid. area1", nPtBins, ptBins ) ;
+    fhTaggedArea2  = new TH1D("fhTaggedArea2","Spectrum of all tagged photons Fid. area2", nPtBins, ptBins ) ;
+    fhTaggedArea3  = new TH1D("fhTaggedArea3","Spectrum of all tagged photons Fid. area3", nPtBins, ptBins ) ;
+    fhTaggedPID[0] = new TH1D("fhTaggedPID0","Spectrum of tagged photons for PID=0", nPtBins, ptBins ) ;
+    fhTaggedPID[1] = new TH1D("fhTaggedPID1","Spectrum of tagged photons for PID=1", nPtBins, ptBins ) ;
+    fhTaggedPID[2] = new TH1D("fhTaggedPID2","Spectrum of tagged photons for PID=2", nPtBins, ptBins ) ;
+    fhTaggedPID[3] = new TH1D("fhTaggedPID3","Spectrum of tagged photons for PID=3", nPtBins, ptBins ) ;
+    fhTaggedMult   = new TH1D("fhTaggedMult","Spectrum of multiply tagged photons", nPtBins, ptBins ) ;
+
+    //Tagging: use MC information if available
+    fhTaggedMCTrue    = new TH1D("fhTaggedMCTrue","Spectrum of correctly tagged pi0 decay photons", nPtBins, ptBins ) ;
+    fhMCMissedTagging = new TH1D("fhMCMissedTagging","Spectrum of pi0 decay photons missed tagging due to wrong pair mass", nPtBins, ptBins ) ;
+    fhMCFakeTagged    = new TH1D("fhMCFakeTagged","Spectrum of photons wrongly tagged according to MC", nPtBins, ptBins ) ;
+
+    //Invariant mass distributions for fake corrections
+    const Int_t nmass=1000 ;
+    Double_t masses[nmass+1] ;
+    for(Int_t i=0;i<=nmass;i++)masses[i]=0.001*i ;
+    fhInvMassReal = new TH2D("fhInvMassReal","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
+    fhInvMassMixed  = new TH2D("fhInvMassMixed","Two-photon inv. mass vs first photon pt",nmass,masses,nPtBins, ptBins ) ;
+    fhMCMissedTaggingMass= new TH2D("fhMCMissedTaggingMass","Inv mass of pairs missed tagging",nmass,masses,nPtBins, ptBins ) ;
+
+    //Conversion and annihilation radius distributions
+    fhConversionRadius  = new TH1D("fhConversionRadius","Radis of photon production (conversion)",100,0.,500.) ;
+    fhInteractionRadius = new TH1D("fhInteractionRadius","Radis of photon production (hadron interaction)",100,0.,500.);
+
+    fhEvents               = new TH1D("hEvents", "Number of Events processed", 1, 0, 1);
+
+// create output container
+  
+  fOutputList = new TList() ;
+  fEventList = new TList() ; 
+  fOutputList->SetName(GetName()) ; 
+
+  fOutputList->AddAt(fhRecAll,                               0) ;
+  fOutputList->AddAt(fhRecAllArea1,                          1) ;
+  fOutputList->AddAt(fhRecAllArea2,                          2) ;
+  fOutputList->AddAt(fhRecAllArea3,                          3) ;
+
+  fOutputList->AddAt(fhRecPhoton,                            4) ;
+  fOutputList->AddAt(fhRecOther,                             5) ;
+  fOutputList->AddAt(fhRecPhotonPID[0],                      6) ;
+  fOutputList->AddAt(fhRecPhotonPID[1],                      7) ;
+  fOutputList->AddAt(fhRecPhotonPID[2],                      8) ;
+  fOutputList->AddAt(fhRecPhotonPID[3],                      9) ;
+  fOutputList->AddAt(fhRecOtherPID[0],                      10) ;
+  fOutputList->AddAt(fhRecOtherPID[1],                      11) ;
+  fOutputList->AddAt(fhRecOtherPID[2],                      12) ;
+  fOutputList->AddAt(fhRecOtherPID[3],                      13) ;
+  fOutputList->AddAt(fhRecPhotPi0,                          14) ;
+  fOutputList->AddAt(fhRecPhotEta,                          15) ;
+  fOutputList->AddAt(fhRecPhotOmega,                        16) ;
+  fOutputList->AddAt(fhRecPhotEtapr,                        17) ;
+  fOutputList->AddAt(fhRecPhotConv,                         18) ;
+  fOutputList->AddAt(fhRecPhotHadron,                       19) ;
+  fOutputList->AddAt(fhRecPhotDirect,                       20) ;
+  fOutputList->AddAt(fhRecPhotOther,                        21) ;
+
+  fOutputList->AddAt(fhDecWMCPartner,                       22) ;
+  fOutputList->AddAt(fhDecWMissedPartnerNotPhoton,          23) ;
+  fOutputList->AddAt(fhDecWMissedPartnerAll,                24) ;
+  fOutputList->AddAt(fhDecWMissedPartnerEmin,               25) ;
+  fOutputList->AddAt(fhDecWMissedPartnerConv,               26) ;
+  fOutputList->AddAt(fhDecWMissedPartnerStack,              27) ;
+  fOutputList->AddAt(fhDecWMissedPartnerGeom0,              28) ;
+  fOutputList->AddAt(fhDecWMissedPartnerGeom1,              29) ;
+  fOutputList->AddAt(fhDecWMissedPartnerGeom2,              30) ;
+  fOutputList->AddAt(fhDecWMissedPartnerGeom3,              31) ;
+
+  fOutputList->AddAt(fhPartnerMCReg,                        32) ;
+  fOutputList->AddAt(fhPartnerMissedEmin,                   33) ;
+  fOutputList->AddAt(fhPartnerMissedConv,                   34) ;
+  fOutputList->AddAt(fhPartnerMissedGeo,                    35) ;
+
+  fOutputList->AddAt(fhTaggedAll,                           36) ;
+  fOutputList->AddAt(fhTaggedArea1,                         37) ;
+  fOutputList->AddAt(fhTaggedArea2,                         38) ;
+  fOutputList->AddAt(fhTaggedArea3,                         39) ;
+  fOutputList->AddAt(fhTaggedPID[0],                        40) ;
+  fOutputList->AddAt(fhTaggedPID[1],                        41) ;
+  fOutputList->AddAt(fhTaggedPID[2],                        42) ;
+  fOutputList->AddAt(fhTaggedPID[3],                        43) ;
+  fOutputList->AddAt(fhTaggedMult,                          44) ;
+
+  fOutputList->AddAt(fhTaggedMCTrue,                        45) ;
+  fOutputList->AddAt(fhMCMissedTagging,                     46) ;
+  fOutputList->AddAt(fhMCFakeTagged,                        47) ;
+
+  fOutputList->AddAt(fhInvMassReal,                         48) ;
+  fOutputList->AddAt(fhInvMassMixed,                        49) ;
+  fOutputList->AddAt(fhMCMissedTaggingMass,                 50) ;
+
+  fOutputList->AddAt(fhConversionRadius,                    51) ;
+  fOutputList->AddAt(fhInteractionRadius,                   52) ;
+
+  fOutputList->AddAt(fhEvents,                              53) ;
+
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *) 
+{
   //Fill all histograms
 
-  fhEvents->Fill(0.);\r
-\r
-  // Processing of one event\r
-  if(fDebug>1)\r
-    AliInfo(Form("\n\n Processing event # %lld",  Entry())) ; \r
-  AliESDEvent* esd = (AliESDEvent*)InputEvent();\r
-\r
-  //MC stack init\r
-  AliMCEventHandler* mctruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());\r
-  fStack = mctruth->MCEvent()->Stack();\r
-\r
-  if(!fStack && gDebug>1)\r
-    AliInfo("No stack! \n");\r
-\r
-  //************************  PHOS/EMCAL *************************************\r
-  TRefArray * caloClustersArr  = new TRefArray();  \r
+  fhEvents->Fill(0.);
+
+  // Processing of one event
+  if(fDebug>1)
+    AliInfo(Form("\n\n Processing event # %lld",  Entry())) ; 
+  AliESDEvent* esd = (AliESDEvent*)InputEvent();
+
+  //MC stack init
+  fStack=0x0 ;
+  if(AliAnalysisManager::GetAnalysisManager()){
+    AliMCEventHandler* mctruth = (AliMCEventHandler*)((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
+    if(mctruth)
+      fStack = mctruth->MCEvent()->Stack();
+  }
+
+  if(!fStack && gDebug>1)
+    AliInfo("No stack! \n");
+
+  //************************  PHOS/EMCAL *************************************
+  TRefArray * caloClustersArr  = new TRefArray();  
   if(fPHOS)
-    esd->GetPHOSClusters(caloClustersArr);\r
+    esd->GetPHOSClusters(caloClustersArr);
   else
-    esd->GetEMCALClusters(caloClustersArr);\r
-  const Int_t kNumberOfClusters = caloClustersArr->GetEntries() ;  \r
-\r
-  TClonesArray * fCaloPhotonsArr   = new TClonesArray("AliAODPWG4Particle",kNumberOfClusters);\r
-  Int_t inList = 0; //counter of caloClusters\r
-\r
-  Int_t cluster ; \r
-\r
-  // loop over Clusters\r
-  for(cluster = 0 ; cluster < kNumberOfClusters ; cluster++) {\r
-    AliESDCaloCluster * caloCluster = static_cast<AliESDCaloCluster*>(caloClustersArr->At(cluster)) ;\r
-  \r
-    if((fPHOS && !caloCluster->IsPHOS()) ||\r
-       (!fPHOS && caloCluster->IsPHOS()))\r
-      continue ; \r
-\r
-    Double_t v[3] ; //vertex ;\r
-    esd->GetVertex()->GetXYZ(v) ;\r
-\r
-    TLorentzVector momentum ;\r
-    caloCluster->GetMomentum(momentum, v);\r
-    new ((*fCaloPhotonsArr)[inList]) AliAODPWG4Particle(momentum.Px(),momentum.Py(),momentum.Pz(),caloCluster->E() );\r
-    AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(inList));\r
-    inList++;\r
-\r
+    esd->GetEMCALClusters(caloClustersArr);
+  const Int_t kNumberOfClusters = caloClustersArr->GetEntries() ;  
+
+  TClonesArray * fCaloPhotonsArr   = new TClonesArray("AliAODPWG4Particle",kNumberOfClusters);
+  Int_t inList = 0; //counter of caloClusters
+
+  Int_t cluster ; 
+
+  // loop over Clusters
+  for(cluster = 0 ; cluster < kNumberOfClusters ; cluster++) {
+    AliESDCaloCluster * caloCluster = static_cast<AliESDCaloCluster*>(caloClustersArr->At(cluster)) ;
+  
+    if((fPHOS && !caloCluster->IsPHOS()) ||
+       (!fPHOS && caloCluster->IsPHOS()))
+      continue ; 
+
+    Double_t v[3] ; //vertex ;
+    esd->GetVertex()->GetXYZ(v) ;
+
+    TLorentzVector momentum ;
+    caloCluster->GetMomentum(momentum, v);
+    new ((*fCaloPhotonsArr)[inList]) AliAODPWG4Particle(momentum.Px(),momentum.Py(),momentum.Pz(),caloCluster->E() );
+    AliAODPWG4Particle *p = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(inList));
+    inList++;
+
     p->SetCaloLabel(cluster,-1); //This and partner cluster
-    p->SetDistToBad(Int_t(caloCluster->GetDistanceToBadChannel()));\r
+    p->SetDistToBad(Int_t(caloCluster->GetDistanceToBadChannel()));
 
-    p->SetTag(AliMCAnalysisUtils::kMCUnknown);\r
-    p->SetTagged(kFALSE);   //Reconstructed pairs found\r
-    p->SetLabel(caloCluster->GetLabel());\r
+    p->SetTag(AliMCAnalysisUtils::kMCUnknown);
+    p->SetTagged(kFALSE);   //Reconstructed pairs found
+    p->SetLabel(caloCluster->GetLabel());
     Float_t pos[3] ;
     caloCluster->GetPosition(pos) ;
     p->SetFiducialArea(GetFiducialArea(pos)) ;
@@ -522,486 +526,434 @@ void AliAnalysisTaskTaggedPhotons::UserExec(Option_t *)
     p->SetDispBit(TestDisp(caloCluster->GetM02(),caloCluster->GetM20(),caloCluster->E())) ;
     p->SetTOFBit(TestTOF(caloCluster->GetTOF(),caloCluster->E())) ;
     p->SetChargedBit(TestCharged(caloCluster->GetEmcCpvDistance(),caloCluster->E())) ;
-\r
-    fhRecAll->Fill( p->Pt() ) ; //All recontructed particles\r
-    Int_t iFidArea = p->GetFiducialArea(); \r
-    if(iFidArea>0){\r
-      fhRecAllArea1->Fill(p->Pt() ) ; \r
-      if(iFidArea>1){\r
-        fhRecAllArea2->Fill(p->Pt() ) ; \r
-        if(iFidArea>2){\r
-          fhRecAllArea3->Fill(p->Pt() ) ; \r
-        }\r
-      }\r
-    }\r
-\r
-    if(fStack){\r
-      TParticle * prim = fStack->Particle(caloCluster->GetLabel()) ;\r
-      if(fDebug>2)\r
-        printf("Pdgcode = %d\n",prim->GetPdgCode());\r
-\r
-      if(prim->GetPdgCode()!=22){ //not photon\r
-//<--DP        p->SetPhoton(kFALSE);\r
-        fhRecOther->Fill(p->Pt()); //not photons real spectra\r
-        for(Int_t iPID=0; iPID<4; iPID++){\r
-          if(p->IsPIDOK(iPID,22))\r
-            fhRecOtherPID[iPID]->Fill(p->Pt());\r
-        }\r
-      }\r
-      else{ //Primary photon (as in MC)\r
-//<--DP        p->SetPhoton(kTRUE);\r
-        fhRecPhoton->Fill(p->Pt()); //Reconstructed with primary photon\r
-        for(Int_t iPID=0; iPID<4; iPID++){\r
-          if(p->IsPIDOK(iPID,22))\r
-            fhRecPhotonPID[iPID]->Fill(p->Pt());\r
-        }\r
-        Int_t pi0i=prim->GetFirstMother();\r
-        Int_t grandpaPDG=-1 ;\r
-        TParticle * pi0p = 0;\r
-        if(pi0i>=0){\r
-          pi0p=fStack->Particle(pi0i);\r
-          grandpaPDG=pi0p->GetPdgCode() ;\r
-        }\r
-        switch(grandpaPDG){\r
-        case 111: //Pi0 decay\r
-                //Primary decay photon (as in MC)\r
-                fhRecPhotPi0->Fill(p->Pt());\r
-                break ;\r
-        case  11:\r
-        case -11: //electron/positron conversion\r
-//<--DP                p->SetConverted(1);\r
-                fhRecPhotConv->Fill(p->Pt());  //Reconstructed with photon from conversion primary\r
-                fhConversionRadius->Fill(prim->R());\r
-                break ;\r
-        case -2212:\r
-        case -2112: //antineutron & antiproton conversion\r
-                fhRecPhotHadron->Fill(p->Pt());  //Reconstructed with photon from antibaryon annihilation\r
-                fhInteractionRadius->Fill(prim->R());\r
-                break ;\r
-          \r
-        case 221: //eta decay\r
-                fhRecPhotEta->Fill(p->Pt());\r
-                break ;  \r
-          \r
-        case 223: //omega meson decay\r
-                fhRecPhotOmega->Fill(p->Pt());\r
-                break ;\r
-            \r
-        case 331: //eta' decay\r
-                fhRecPhotEtapr->Fill(p->Pt());\r
-                break ;\r
-              \r
-        case -1: //direct photon or no primary\r
-                fhRecPhotDirect->Fill(p->Pt());\r
-                break ;\r
-              \r
-        default:  \r
-                fhRecPhotOther->Fill(p->Pt());\r
-                break ;\r
-        }  \r
-\r
-        //Now classify pi0 decay photon\r
-        if(grandpaPDG==111){\r
-//<--DP          p->Pi0Decay(kTRUE); //Mark this photon as primary decayed\r
-//<--DP          p->Pi0Id(pi0i); //remember id of the parent pi0\r
-\r
-          //Now check if second (partner) photon from pi0 decay hits PHOS or not\r
-          //i.e. both photons can be tagged or it's the systematic error\r
-//<--DP          p->SetPartnerPt(0.); \r
-          Int_t indexdecay1,indexdecay2;\r
-\r
-          indexdecay1=pi0p->GetFirstDaughter();\r
-          indexdecay2=pi0p->GetLastDaughter();\r
-          Int_t indexdecay=-1;\r
-          if(fDebug>2)\r
-                printf("checking pi0 decay...index1=%d, index2=%d, index_pi0=%d, index_ph_prim=%d\n", indexdecay1,indexdecay2,pi0i,caloCluster->GetLabel());\r
-                \r
-          if(indexdecay1!=caloCluster->GetLabel()) \r
-            indexdecay=indexdecay1;\r
-          if(indexdecay2!=caloCluster->GetLabel()) \r
-            indexdecay=indexdecay2;\r
-          if(indexdecay==-1){\r
-            if(fDebug>2){\r
-              printf("Probably the other photon is not in the stack!\n");\r
-              printf("Number of daughters: %d\n",pi0p->GetNDaughters());\r
-            }\r
+
+    fhRecAll->Fill( p->Pt() ) ; //All recontructed particles
+    Int_t iFidArea = p->GetFiducialArea(); 
+    if(iFidArea>0){
+      fhRecAllArea1->Fill(p->Pt() ) ; 
+      if(iFidArea>1){
+        fhRecAllArea2->Fill(p->Pt() ) ; 
+        if(iFidArea>2){
+          fhRecAllArea3->Fill(p->Pt() ) ; 
+        }
+      }
+    }
+
+    if(fStack){
+      TParticle * prim = fStack->Particle(caloCluster->GetLabel()) ;
+      if(fDebug>2)
+        printf("Pdgcode = %d\n",prim->GetPdgCode());
+
+      if(prim->GetPdgCode()!=22){ //not photon
+        fhRecOther->Fill(p->Pt()); //not photons real spectra
+        for(Int_t iPID=0; iPID<4; iPID++){
+          if(p->IsPIDOK(iPID,22))
+            fhRecOtherPID[iPID]->Fill(p->Pt());
+        }
+      }
+      else{ //Primary photon (as in MC)
+        fhRecPhoton->Fill(p->Pt()); //Reconstructed with primary photon
+        for(Int_t iPID=0; iPID<4; iPID++){
+          if(p->IsPIDOK(iPID,22))
+            fhRecPhotonPID[iPID]->Fill(p->Pt());
+        }
+        Int_t pi0i=prim->GetFirstMother();
+        Int_t grandpaPDG=-1 ;
+        TParticle * pi0p = 0;
+        if(pi0i>=0){
+          pi0p=fStack->Particle(pi0i);
+          grandpaPDG=pi0p->GetPdgCode() ;
+        }
+        switch(grandpaPDG){
+        case 111: //Pi0 decay
+                //Primary decay photon (as in MC)
+                fhRecPhotPi0->Fill(p->Pt());
+                break ;
+        case  11:
+        case -11: //electron/positron conversion
+                fhRecPhotConv->Fill(p->Pt());  //Reconstructed with photon from conversion primary
+                fhConversionRadius->Fill(prim->R());
+                break ;
+        case -2212:
+        case -2112: //antineutron & antiproton conversion
+                fhRecPhotHadron->Fill(p->Pt());  //Reconstructed with photon from antibaryon annihilation
+                fhInteractionRadius->Fill(prim->R());
+                break ;
+          
+        case 221: //eta decay
+                fhRecPhotEta->Fill(p->Pt());
+                break ;  
+          
+        case 223: //omega meson decay
+                fhRecPhotOmega->Fill(p->Pt());
+                break ;
+            
+        case 331: //eta' decay
+                fhRecPhotEtapr->Fill(p->Pt());
+                break ;
+              
+        case -1: //direct photon or no primary
+                fhRecPhotDirect->Fill(p->Pt());
+                break ;
+              
+        default:  
+                fhRecPhotOther->Fill(p->Pt());
+                break ;
+        }  
+
+        //Now classify pi0 decay photon
+        if(grandpaPDG==111){
+//<--DP          p->Pi0Decay(kTRUE); //Mark this photon as primary decayed
+//<--DP          p->Pi0Id(pi0i); //remember id of the parent pi0
+
+          //Now check if second (partner) photon from pi0 decay hits PHOS or not
+          //i.e. both photons can be tagged or it's the systematic error
+//<--DP          p->SetPartnerPt(0.); 
+          Int_t indexdecay1,indexdecay2;
+
+          indexdecay1=pi0p->GetFirstDaughter();
+          indexdecay2=pi0p->GetLastDaughter();
+          Int_t indexdecay=-1;
+          if(fDebug>2)
+                printf("checking pi0 decay...index1=%d, index2=%d, index_pi0=%d, index_ph_prim=%d\n", indexdecay1,indexdecay2,pi0i,caloCluster->GetLabel());
+                
+          if(indexdecay1!=caloCluster->GetLabel()) 
+            indexdecay=indexdecay1;
+          if(indexdecay2!=caloCluster->GetLabel()) 
+            indexdecay=indexdecay2;
+          if(indexdecay==-1){
+            if(fDebug>2){
+              printf("Probably the other photon is not in the stack!\n");
+              printf("Number of daughters: %d\n",pi0p->GetNDaughters());
+            }
             fhDecWMissedPartnerStack->Fill(p->Pt()) ;
-          }  \r
-          else{\r
-            TParticle *partner = fStack->Particle(indexdecay);\r
-//<--DP            p->SetPartnerPt(partner->Pt());\r
-            if(partner->GetPdgCode()==22){ \r
-              Bool_t isPartnerLost=kFALSE; //If partner is lost for some reason\r
-              if(partner->GetNDaughters()!=0){ //this photon is converted before it is registered by some detector\r
-                if(fDebug>2)\r
-                  printf("P_Conv, daughters=%d\n",partner->GetNDaughters());\r
-//<--DP                p->SetConvertedPartner(1);\r
-                fhPartnerMissedConv->Fill(partner->Pt());\r
-                fhDecWMissedPartnerConv->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
-                isPartnerLost=kTRUE;\r
-              }\r
+          }  
+          else{
+            TParticle *partner = fStack->Particle(indexdecay);
+//<--DP            p->SetPartnerPt(partner->Pt());
+            if(partner->GetPdgCode()==22){ 
+              Bool_t isPartnerLost=kFALSE; //If partner is lost for some reason
+              if(partner->GetNDaughters()!=0){ //this photon is converted before it is registered by some detector
+                if(fDebug>2)
+                  printf("P_Conv, daughters=%d\n",partner->GetNDaughters());
+//<--DP                p->SetConvertedPartner(1);
+                fhPartnerMissedConv->Fill(partner->Pt());
+                fhDecWMissedPartnerConv->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
+                isPartnerLost=kTRUE;
+              }
               Bool_t impact = kFALSE ;
               if(fPHOS){
                 Int_t modulenum ;
                 Double_t ztmp,xtmp ;
                 impact=fPHOSgeom->ImpactOnEmc(partner,modulenum,ztmp,xtmp) ;
-                if(fDebug>2){\r
-                  printf("Impact on PHOS: module: %d, x tower: %f, z tower: %f\n", modulenum,xtmp,ztmp);\r
+                if(fDebug>2){
+                  printf("Impact on PHOS: module: %d, x tower: %f, z tower: %f\n", modulenum,xtmp,ztmp);
                 }
               }
               else{
                 impact = fEMCALgeom->Impact(partner) ;
               }
-              if(!impact){ //this photon cannot hit PHOS\r
-                if(fDebug>2)\r
-                  printf("P_Geo\n");\r
-                fhPartnerMissedGeo->Fill(partner->Pt());\r
-                fhDecWMissedPartnerGeom0->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
-                if(iFidArea>0){\r
-                  fhDecWMissedPartnerGeom1->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
-                  if(iFidArea>1){\r
-                    fhDecWMissedPartnerGeom2->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
-                    if(iFidArea>2){\r
-                      fhDecWMissedPartnerGeom3->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
-                    }\r
-                  }\r
-                }\r
-                isPartnerLost=kTRUE;\r
-              }\r
-              if(!isPartnerLost && partner->Energy()<fMinEnergyCut){ //energy is not enough to be registered by PHOS\r
-                if(fDebug>2)\r
-                  printf("P_Reg, E=%f\n",partner->Energy());\r
-                fhPartnerMissedEmin->Fill(partner->Pt());  //Spectrum of missed partners\r
-                fhDecWMissedPartnerEmin->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner\r
-                isPartnerLost=kTRUE;\r
-              }\r
-              if(!isPartnerLost){\r
-//                p->SetMCTagged(1); //set this photon as primary tagged\r
-                fhDecWMCPartner->Fill(p->Pt());\r
-                fhPartnerMCReg->Fill(partner->Pt());\r
-                if(fDebug>2){\r
-                  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
-                }\r
-              }\r
-              else{\r
-                fhDecWMissedPartnerAll->Fill(p->Pt());\r
-              }\r
-            }//Partner - photon\r
-            else{//partner not photon\r
-              fhDecWMissedPartnerNotPhoton->Fill(p->Pt());                \r
-            }\r
-          }\r
-        }\r
-      }\r
-    }\r
-  } //PHOS/EMCAL clusters\r
-    \r
-  if(fDebug>1)   \r
-    printf("number of clusters: %d\n",inList);\r
-\r
-  //Invariant Mass analysis\r
-  for(Int_t phosPhoton1 = 0 ; phosPhoton1 < inList-1 ; phosPhoton1++) {\r
-    AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton1));        \r
-    for(Int_t phosPhoton2 = phosPhoton1 + 1 ; phosPhoton2 < inList ; phosPhoton2++) {\r
-      AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton2));\r
-\r
-      Double_t invMass = p1->GetPairMass(p2);\r
-      fhInvMassReal->Fill(invMass,p1->Pt());\r
-      fhInvMassReal->Fill(invMass,p2->Pt());\r
-      if(fDebug>2)\r
-          printf("Pair i=%d,j=%d, M=%f\n",phosPhoton1,phosPhoton2,invMass);\r
-\r
-      Bool_t makePi01=IsInPi0Band(invMass,p1->Pt());\r
-      Bool_t makePi02=IsInPi0Band(invMass,p2->Pt());\r
-\r
-\r
-      if(makePi01 && p1->IsTagged()){//Multiple tagging\r
-        fhTaggedMult->Fill(p1->Pt());\r
-      }  \r
-      if(makePi01 && !p1->IsTagged()){//Each photon should enter histogram once, even if tagged several times\r
-        fhTaggedAll->Fill(p1->Pt());\r
-        Int_t iFidArea = p1->GetFiducialArea(); \r
-        if(iFidArea>0){\r
-          fhTaggedArea1->Fill(p1->Pt() ) ; \r
-          if(iFidArea>1){\r
-            fhTaggedArea2->Fill(p1->Pt() ) ; \r
-            if(iFidArea>2){\r
-              fhTaggedArea3->Fill(p1->Pt() ) ; \r
-            }\r
-          }\r
-        }\r
-\r
-        for(Int_t iPID=0; iPID<4; iPID++){\r
-          if(p1->IsPIDOK(iPID,22))\r
-            fhTaggedPID[iPID]->Fill(p1->Pt());\r
-        }\r
-\r
-        p1->SetTagged(kTRUE) ;\r
-      }  \r
-      if(makePi02 && p2->IsTagged()){//Multiple tagging\r
-        fhTaggedMult->Fill(p2->Pt());\r
-      }  \r
-      if(makePi02 && !p2->IsTagged()){//How should be account for multiply tagged photons?\r
-        fhTaggedAll->Fill(p2->Pt());\r
-        p2->SetTagged(kTRUE) ;\r
-      }\r
-      \r
-      //Now get use MC information\r
-      //First chesk if this is true pi0 pair\r
-      if(IsSamePi0(p1,p2)){ //Correctly tagged - from the same pi0\r
-//        p1->SetTrueTagged(1);\r
-//        p2->SetTrueTagged(1);\r
-        if(makePi01)//Correctly tagged photons\r
-          fhTaggedMCTrue->Fill(p1->Pt());\r
-        else{ //Decay pair missed tagging      \r
-          fhMCMissedTagging->Fill(p1->Pt());\r
-          fhMCMissedTaggingMass->Fill(invMass,p1->Pt()) ;\r
-          //Clussify why missed tagging (todo)\r
-          //Converted\r
-          //Partner not a photon\r
-          //Tagged not a photon\r
-          //Just wrong inv.mass          \r
-        }  \r
-        if(makePi02)\r
-          fhTaggedMCTrue->Fill(p2->Pt());\r
-        else{      \r
-          fhMCMissedTagging->Fill(p2->Pt());\r
-          fhMCMissedTaggingMass->Fill(invMass,p2->Pt()) ;\r
-          //Clussify why missed tagging (todo)\r
-          //Converted\r
-          //Partner not a photon\r
-          //Tagged not a photon\r
-          //Just wrong inv.mass          \r
-        }  \r
-      }\r
-      else{//Fake tagged - not from the same pi0\r
-        if(makePi01)//Fake pair\r
-          fhMCFakeTagged->Fill(p1->Pt());\r
-        if(makePi02)\r
-          fhMCFakeTagged->Fill(p2->Pt());\r
-      }\r
-    }\r
-  }\r
-\r
-  //Fill Mixed InvMass distributions:\r
-  TIter nextEv(fEventList) ;\r
-  while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){\r
-    Int_t nPhotons2 = event2->GetEntriesFast() ;\r
-    for(Int_t i=0; i < inList ; i++){\r
-      AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(i)) ;\r
-      for(Int_t j=0; j < nPhotons2 ; j++){\r
-        AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(event2->At(j)) ;\r
-        Double_t invMass = p1->GetPairMass(p2);\r
-        fhInvMassMixed->Fill(invMass,p1->Pt());\r
-        fhInvMassMixed->Fill(invMass,p2->Pt());\r
-      }\r
-    }\r
-  }\r
-\r
-  //Remove old events\r
-  fEventList->AddFirst(fCaloPhotonsArr);\r
-  if(fEventList->GetSize() > 10){\r
-    TClonesArray *tmp = static_cast <TClonesArray*> (fEventList->Last());\r
-    fEventList->Remove(tmp);\r
-    delete tmp;\r
-  }\r
-\r
-  PostData(1, fOutputList);\r
-}\r
-\r
-\r
-//______________________________________________________________________________\r
-void AliAnalysisTaskTaggedPhotons::Init()\r
-{\r
-  // Intialisation of parameters\r
-  AliInfo("Doing initialisation") ; \r
-  SetPhotonId(0.9) ; \r
-  SetMinEnergyCut(0.4);\r
-  SetPi0MeanParameters(0.136,0.,0.0,0.0);\r
-//  SetPi0MeanParameters(0.1377,-0.002566,0.001216,-0.0001256);\r
-  SetPi0SigmaParameters(0.004508,0.005497,0.00000006);\r
-}\r
-\r
-//______________________________________________________________________________\r
-void AliAnalysisTaskTaggedPhotons::Terminate(Option_t *)\r
-{\r
-  // Processing when the event loop is ended\r
-\r
-  //Write everything to the file\r
-  char outname[55];\r
-  if(fPHOS)\r
-    sprintf(outname,"Tagging_PHOS.root") ;\r
-  else  \r
-    sprintf(outname,"Tagging_EMCAL.root") ;\r
-  TFile *outfile = new TFile (outname,"recreate");\r
-\r
-fhRecAll->Write();\r
-fhRecAllArea1->Write();\r
-fhRecAllArea2->Write();\r
-fhRecAllArea3->Write();\r
-fhRecPhoton->Write();\r
-fhRecOther->Write();\r
-fhRecPhotonPID[0]->Write();\r
-fhRecPhotonPID[1]->Write();\r
-fhRecPhotonPID[2]->Write();\r
-fhRecPhotonPID[3]->Write();\r
-fhRecOtherPID[0]->Write();\r
-fhRecOtherPID[1]->Write();\r
-fhRecOtherPID[2]->Write();\r
-fhRecOtherPID[3]->Write();\r
-fhRecPhotPi0->Write();\r
-fhRecPhotEta->Write();\r
-fhRecPhotOmega->Write();\r
-fhRecPhotEtapr->Write();\r
-fhRecPhotConv->Write();\r
-fhRecPhotHadron->Write();\r
-fhRecPhotDirect->Write();\r
-fhRecPhotOther->Write();\r
-fhDecWMCPartner->Write();\r
-fhDecWMissedPartnerNotPhoton->Write();\r
-fhDecWMissedPartnerAll->Write();\r
-fhDecWMissedPartnerEmin->Write();\r
-fhDecWMissedPartnerConv->Write();\r
-fhDecWMissedPartnerStack->Write();\r
-fhDecWMissedPartnerGeom0->Write();\r
-fhDecWMissedPartnerGeom1->Write();\r
-fhDecWMissedPartnerGeom2->Write();\r
-fhDecWMissedPartnerGeom3->Write();\r
-fhPartnerMCReg->Write();\r
-fhPartnerMissedEmin->Write();\r
-fhPartnerMissedConv->Write();\r
-fhPartnerMissedGeo->Write();\r
-fhTaggedAll->Write();\r
-fhTaggedArea1->Write();\r
-fhTaggedArea2->Write();\r
-fhTaggedArea3->Write();\r
-\r
-fhTaggedPID[0]->Write();\r
-fhTaggedPID[1]->Write();\r
-fhTaggedPID[2]->Write();\r
-fhTaggedPID[3]->Write();\r
-fhTaggedMult->Write();\r
-fhTaggedMCTrue->Write();\r
-fhMCMissedTagging->Write();\r
-fhMCFakeTagged->Write();\r
-fhInvMassReal->Write();\r
-fhInvMassMixed->Write();\r
-fhMCMissedTaggingMass->Write();\r
-fhConversionRadius->Write();\r
-fhInteractionRadius->Write();\r
-fhEvents->Write();\r
-\r
-/*\r
-  fhAllPhotons->Write() ;\r
-  fhNotPhotons->Write() ;\r
-  fhAllPhotonsPrimary->Write() ;\r
-  fhNotPhotonsPrimary->Write() ;\r
-  fhfakeNotPhotons->Write() ;\r
-\r
-  fhTaggedPhotons->Write();\r
-  fhfakeTaggedPhotons->Write();\r
-  fhDecayNotTaggedPhotons->Write();\r
-  fhstrangeNotTaggedPhotons->Write();\r
-  fhstrangeNotTaggedPhotonsPair->Write();\r
-  fhstrangeNotTaggedPhotonsRegCut->Write();\r
-  fhstrangeNotTaggedPhotonsPairRegCut->Write();\r
-\r
-  fhPi0DecayPhotonsPrimary->Write();\r
-  fhEtaDecayPhotonsPrimary->Write();\r
-  fhOmegaDecayPhotonsPrimary->Write();\r
-  fhEtaSDecayPhotonsPrimary->Write();\r
-  fhOtherDecayPhotonsPrimary->Write();\r
-  fhDecayPhotonsPrimary->Write();\r
-  fhConvertedPhotonsPrimary->Write();\r
-  fhConvertedPhotonsPrimaryHadronsDecays->Write();\r
-  fhCoordsConvertion->Write();\r
-  fhCoordsConvertion2->Write();\r
-\r
-  fhPHOSInvariantMassReal->Write();\r
-  fhPHOSInvariantMassMixed->Write();\r
-\r
-  fhPHOSPi0->Write();\r
-\r
-  fhPi0DecayPhotonsGeomfake->Write();\r
-  fhPi0DecayPhotonsTaggedPrimary->Write();\r
-  fhPi0DecayPhotonsTaggedPrimaryPair->Write();\r
-  fhPi0DecayPhotonsTaggedPrimaryRegCut->Write();\r
-  fhPi0DecayPhotonsTaggedPrimaryPairRegCut->Write();\r
-\r
-  fhPi0DecayPhotonsBigDecay->Write();\r
-  fhPi0DecayPhotonsPConv->Write();\r
-  fhPi0DecayPhotonsPGeo->Write();\r
-  fhPi0DecayPhotonsPReg->Write();\r
-\r
-  fhfakeTaggedPhotonsConv->Write();\r
-  fhfakeTaggedPhotonsPID->Write();\r
-\r
-  fhTrackRefCoords->Write();\r
-  fhEvents->Write();\r
-*/\r
-\r
-outfile->Close();\r
-\r
-}\r
-//______________________________________________________________________________\r
-Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)const\r
-{\r
-  //Parameterization of the pi0 peak region\r
-  Double_t mpi0mean = fPi0MeanP0 + fPi0MeanP1 * pt + fPi0MeanP2 * pt*pt + fPi0MeanP3 * pt*pt*pt;\r
-  Double_t mpi0sigma = TMath::Sqrt(fPi0SigmaP0 * fPi0SigmaP0 / pt + fPi0SigmaP1 * fPi0SigmaP1 + fPi0SigmaP2 * fPi0SigmaP2 / pt / pt);\r
\r
-  return (m>mpi0mean-2*mpi0sigma && m<mpi0mean+2*mpi0sigma) ;\r
-}\r
-//______________________________________________________________________________\r
-Bool_t AliAnalysisTaskTaggedPhotons::IsSamePi0(const AliAODPWG4Particle *p1, const AliAODPWG4Particle *p2)const{\r
-  //Looks through parents and finds if there was commont pi0 among ancestors\r
-\r
-  if(!fStack)\r
-    return kFALSE ; //can not say anything\r
-\r
-  Int_t prim1 = p1->GetLabel();\r
-  while(prim1!=-1){ \r
-    Int_t prim2 = p2->GetLabel();\r
-    while(prim2!=-1){ \r
-      if(prim1==prim2){\r
-        if(fStack->Particle(prim1)->GetPdgCode()==111)\r
-          return kTRUE ;\r
-        else\r
-          return kFALSE ;\r
-      }\r
-      prim2=fStack->Particle(prim2)->GetFirstMother() ;\r
-    }\r
-    prim1=fStack->Particle(prim1)->GetFirstMother() ;\r
-  }\r
-  return kFALSE ;\r
-}\r
-//______________________________________________________________________________\r
-Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(Float_t * pos)const{\r
-  //calculates in which kind of fiducial area photon hit\r
-  Double_t phi=TMath::ATan2(pos[1],pos[0]) ;\r
+              if(!impact){ //this photon cannot hit PHOS
+                if(fDebug>2)
+                  printf("P_Geo\n");
+                fhPartnerMissedGeo->Fill(partner->Pt());
+                fhDecWMissedPartnerGeom0->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
+                if(iFidArea>0){
+                  fhDecWMissedPartnerGeom1->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
+                  if(iFidArea>1){
+                    fhDecWMissedPartnerGeom2->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
+                    if(iFidArea>2){
+                      fhDecWMissedPartnerGeom3->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
+                    }
+                  }
+                }
+                isPartnerLost=kTRUE;
+              }
+              if(!isPartnerLost && partner->Energy()<fMinEnergyCut){ //energy is not enough to be registered by PHOS
+                if(fDebug>2)
+                  printf("P_Reg, E=%f\n",partner->Energy());
+                fhPartnerMissedEmin->Fill(partner->Pt());  //Spectrum of missed partners
+                fhDecWMissedPartnerEmin->Fill(p->Pt()) ;  //Spectrum of tagged with missed partner
+                isPartnerLost=kTRUE;
+              }
+              if(!isPartnerLost){
+//                p->SetMCTagged(1); //set this photon as primary tagged
+                fhDecWMCPartner->Fill(p->Pt());
+                fhPartnerMCReg->Fill(partner->Pt());
+                if(fDebug>2){
+                  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());
+                }
+              }
+              else{
+                fhDecWMissedPartnerAll->Fill(p->Pt());
+              }
+            }//Partner - photon
+            else{//partner not photon
+              fhDecWMissedPartnerNotPhoton->Fill(p->Pt());                
+            }
+          }
+        }
+      }
+    }
+  } //PHOS/EMCAL clusters
+    
+  if(fDebug>1)   
+    printf("number of clusters: %d\n",inList);
+
+  //Invariant Mass analysis
+  for(Int_t phosPhoton1 = 0 ; phosPhoton1 < inList-1 ; phosPhoton1++) {
+    AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton1));        
+    for(Int_t phosPhoton2 = phosPhoton1 + 1 ; phosPhoton2 < inList ; phosPhoton2++) {
+      AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(phosPhoton2));
+
+      Double_t invMass = p1->GetPairMass(p2);
+      fhInvMassReal->Fill(invMass,p1->Pt());
+      fhInvMassReal->Fill(invMass,p2->Pt());
+      if(fDebug>2)
+          printf("Pair i=%d,j=%d, M=%f\n",phosPhoton1,phosPhoton2,invMass);
+
+      Bool_t makePi01=IsInPi0Band(invMass,p1->Pt());
+      Bool_t makePi02=IsInPi0Band(invMass,p2->Pt());
+
+
+      if(makePi01 && p1->IsTagged()){//Multiple tagging
+        fhTaggedMult->Fill(p1->Pt());
+      }  
+      if(makePi01 && !p1->IsTagged()){//Each photon should enter histogram once, even if tagged several times
+        fhTaggedAll->Fill(p1->Pt());
+        Int_t iFidArea = p1->GetFiducialArea(); 
+        if(iFidArea>0){
+          fhTaggedArea1->Fill(p1->Pt() ) ; 
+          if(iFidArea>1){
+            fhTaggedArea2->Fill(p1->Pt() ) ; 
+            if(iFidArea>2){
+              fhTaggedArea3->Fill(p1->Pt() ) ; 
+            }
+          }
+        }
+
+        for(Int_t iPID=0; iPID<4; iPID++){
+          if(p1->IsPIDOK(iPID,22))
+            fhTaggedPID[iPID]->Fill(p1->Pt());
+        }
+
+        p1->SetTagged(kTRUE) ;
+      }  
+      if(makePi02 && p2->IsTagged()){//Multiple tagging
+        fhTaggedMult->Fill(p2->Pt());
+      }  
+      if(makePi02 && !p2->IsTagged()){//How should be account for multiply tagged photons?
+        fhTaggedAll->Fill(p2->Pt());
+        p2->SetTagged(kTRUE) ;
+      }
+      
+      //Now get use MC information
+      //First chesk if this is true pi0 pair
+      if(IsSamePi0(p1,p2)){ //Correctly tagged - from the same pi0
+//        p1->SetTrueTagged(1);
+//        p2->SetTrueTagged(1);
+        if(makePi01)//Correctly tagged photons
+          fhTaggedMCTrue->Fill(p1->Pt());
+        else{ //Decay pair missed tagging      
+          fhMCMissedTagging->Fill(p1->Pt());
+          fhMCMissedTaggingMass->Fill(invMass,p1->Pt()) ;
+          //Clussify why missed tagging (todo)
+          //Converted
+          //Partner not a photon
+          //Tagged not a photon
+          //Just wrong inv.mass          
+        }  
+        if(makePi02)
+          fhTaggedMCTrue->Fill(p2->Pt());
+        else{      
+          fhMCMissedTagging->Fill(p2->Pt());
+          fhMCMissedTaggingMass->Fill(invMass,p2->Pt()) ;
+          //Clussify why missed tagging (todo)
+          //Converted
+          //Partner not a photon
+          //Tagged not a photon
+          //Just wrong inv.mass          
+        }  
+      }
+      else{//Fake tagged - not from the same pi0
+        if(makePi01)//Fake pair
+          fhMCFakeTagged->Fill(p1->Pt());
+        if(makePi02)
+          fhMCFakeTagged->Fill(p2->Pt());
+      }
+    }
+  }
+
+  //Fill Mixed InvMass distributions:
+  TIter nextEv(fEventList) ;
+  while(TClonesArray * event2 = static_cast<TClonesArray*>(nextEv())){
+    Int_t nPhotons2 = event2->GetEntriesFast() ;
+    for(Int_t i=0; i < inList ; i++){
+      AliAODPWG4Particle * p1 = static_cast<AliAODPWG4Particle*>(fCaloPhotonsArr->At(i)) ;
+      for(Int_t j=0; j < nPhotons2 ; j++){
+        AliAODPWG4Particle * p2 = static_cast<AliAODPWG4Particle*>(event2->At(j)) ;
+        Double_t invMass = p1->GetPairMass(p2);
+        fhInvMassMixed->Fill(invMass,p1->Pt());
+        fhInvMassMixed->Fill(invMass,p2->Pt());
+      }
+    }
+  }
+
+  //Remove old events
+  fEventList->AddFirst(fCaloPhotonsArr);
+  if(fEventList->GetSize() > 10){
+    TClonesArray *tmp = static_cast <TClonesArray*> (fEventList->Last());
+    fEventList->Remove(tmp);
+    delete tmp;
+  }
+
+  PostData(1, fOutputList);
+}
+
+
+//______________________________________________________________________________
+void AliAnalysisTaskTaggedPhotons::Init()
+{
+  // Intialisation of parameters
+  AliInfo("Doing initialisation") ; 
+  SetPhotonId(0.9) ; 
+  SetMinEnergyCut(0.4);
+  SetPi0MeanParameters(0.136,0.,0.0,0.0);
+//  SetPi0MeanParameters(0.1377,-0.002566,0.001216,-0.0001256);
+  SetPi0SigmaParameters(0.004508,0.005497,0.00000006);
+}
+
+//______________________________________________________________________________
+void AliAnalysisTaskTaggedPhotons::Terminate(Option_t *)
+{
+  // Processing when the event loop is ended
+
+  //Write everything to the file
+  char outname[55];
+  if(fPHOS)
+    sprintf(outname,"Tagging_PHOS.root") ;
+  else  
+    sprintf(outname,"Tagging_EMCAL.root") ;
+  TFile *outfile = new TFile (outname,"recreate");
+
+fhRecAll->Write();
+fhRecAllArea1->Write();
+fhRecAllArea2->Write();
+fhRecAllArea3->Write();
+fhRecPhoton->Write();
+fhRecOther->Write();
+fhRecPhotonPID[0]->Write();
+fhRecPhotonPID[1]->Write();
+fhRecPhotonPID[2]->Write();
+fhRecPhotonPID[3]->Write();
+fhRecOtherPID[0]->Write();
+fhRecOtherPID[1]->Write();
+fhRecOtherPID[2]->Write();
+fhRecOtherPID[3]->Write();
+fhRecPhotPi0->Write();
+fhRecPhotEta->Write();
+fhRecPhotOmega->Write();
+fhRecPhotEtapr->Write();
+fhRecPhotConv->Write();
+fhRecPhotHadron->Write();
+fhRecPhotDirect->Write();
+fhRecPhotOther->Write();
+fhDecWMCPartner->Write();
+fhDecWMissedPartnerNotPhoton->Write();
+fhDecWMissedPartnerAll->Write();
+fhDecWMissedPartnerEmin->Write();
+fhDecWMissedPartnerConv->Write();
+fhDecWMissedPartnerStack->Write();
+fhDecWMissedPartnerGeom0->Write();
+fhDecWMissedPartnerGeom1->Write();
+fhDecWMissedPartnerGeom2->Write();
+fhDecWMissedPartnerGeom3->Write();
+fhPartnerMCReg->Write();
+fhPartnerMissedEmin->Write();
+fhPartnerMissedConv->Write();
+fhPartnerMissedGeo->Write();
+fhTaggedAll->Write();
+fhTaggedArea1->Write();
+fhTaggedArea2->Write();
+fhTaggedArea3->Write();
+
+fhTaggedPID[0]->Write();
+fhTaggedPID[1]->Write();
+fhTaggedPID[2]->Write();
+fhTaggedPID[3]->Write();
+fhTaggedMult->Write();
+fhTaggedMCTrue->Write();
+fhMCMissedTagging->Write();
+fhMCFakeTagged->Write();
+fhInvMassReal->Write();
+fhInvMassMixed->Write();
+fhMCMissedTaggingMass->Write();
+fhConversionRadius->Write();
+fhInteractionRadius->Write();
+fhEvents->Write();
+
+outfile->Close();
+
+}
+//______________________________________________________________________________
+Bool_t AliAnalysisTaskTaggedPhotons::IsInPi0Band(Double_t m, Double_t pt)const
+{
+  //Parameterization of the pi0 peak region
+  Double_t mpi0mean = fPi0MeanP0 + fPi0MeanP1 * pt + fPi0MeanP2 * pt*pt + fPi0MeanP3 * pt*pt*pt;
+  Double_t mpi0sigma = TMath::Sqrt(fPi0SigmaP0 * fPi0SigmaP0 / pt + fPi0SigmaP1 * fPi0SigmaP1 + fPi0SigmaP2 * fPi0SigmaP2 / pt / pt);
+  return (m>mpi0mean-2*mpi0sigma && m<mpi0mean+2*mpi0sigma) ;
+}
+//______________________________________________________________________________
+Bool_t AliAnalysisTaskTaggedPhotons::IsSamePi0(const AliAODPWG4Particle *p1, const AliAODPWG4Particle *p2)const{
+  //Looks through parents and finds if there was commont pi0 among ancestors
+
+  if(!fStack)
+    return kFALSE ; //can not say anything
+
+  Int_t prim1 = p1->GetLabel();
+  while(prim1!=-1){ 
+    Int_t prim2 = p2->GetLabel();
+    while(prim2!=-1){ 
+      if(prim1==prim2){
+        if(fStack->Particle(prim1)->GetPdgCode()==111)
+          return kTRUE ;
+        else
+          return kFALSE ;
+      }
+      prim2=fStack->Particle(prim2)->GetFirstMother() ;
+    }
+    prim1=fStack->Particle(prim1)->GetFirstMother() ;
+  }
+  return kFALSE ;
+}
+//______________________________________________________________________________
+Int_t AliAnalysisTaskTaggedPhotons::GetFiducialArea(Float_t * pos)const{
+  //calculates in which kind of fiducial area photon hit
+  Double_t phi=TMath::ATan2(pos[1],pos[0]) ;
   Double_t z=pos[2] ;
-  while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;\r
-  while(phi<0.)phi+=TMath::TwoPi() ;\r
-  if(fPHOS){\r
-    //From active PHOS area remove bands in 10 cm\r
-    const Double_t kphi=TMath::ATan(10./460.) ; //angular band width\r
-    Double_t dzMax=TMath::Ceil((fZmax-z)/10.) ;\r
-    Double_t dzMin=TMath::Ceil((z-fZmin)/10.) ;\r
-    Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);\r
-    Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);\r
-    return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin)); \r
-  }\r
-  else{//EMCAL\r
-    //From active EMCAL area remove bands in 20 cm\r
-    const Double_t kphi=TMath::ATan(20./428.) ; //angular band width\r
-    Double_t dzMax=TMath::Ceil((fZmax-z)/20.) ;\r
-    Double_t dzMin=TMath::Ceil((z-fZmin)/20.) ;\r
-    Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);\r
-    Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);\r
-    return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin)); \r
-  }\r
-}\r
-//______________________________________________________________________________\r
+  while(phi>TMath::TwoPi())phi-=TMath::TwoPi() ;
+  while(phi<0.)phi+=TMath::TwoPi() ;
+  if(fPHOS){
+    //From active PHOS area remove bands in 10 cm
+    const Double_t kphi=TMath::ATan(10./460.) ; //angular band width
+    Double_t dzMax=TMath::Ceil((fZmax-z)/10.) ;
+    Double_t dzMin=TMath::Ceil((z-fZmin)/10.) ;
+    Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);
+    Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);
+    return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin)); 
+  }
+  else{//EMCAL
+    //From active EMCAL area remove bands in 20 cm
+    const Double_t kphi=TMath::ATan(20./428.) ; //angular band width
+    Double_t dzMax=TMath::Ceil((fZmax-z)/20.) ;
+    Double_t dzMin=TMath::Ceil((z-fZmin)/20.) ;
+    Double_t dphiMax=TMath::Ceil((fPhimax-phi)/kphi);
+    Double_t dphiMin=TMath::Ceil((phi-fPhimin)/kphi);
+    return (Int_t)TMath::Min(TMath::Min(dzMax,dzMin),TMath::Min(dphiMax,dphiMin)); 
+  }
+}
+//______________________________________________________________________________
 Bool_t  AliAnalysisTaskTaggedPhotons::TestDisp(Double_t l0, Double_t l1, Double_t e)const{
   //test if dispersion corresponds to those of photon
   if(fPHOS){
index 12830f4..0deabd4 100644 (file)
-#ifndef ALIANALYSISTASKTAGGEDPHOTONS_H\r
-#define ALIANALYSISTASKTAGGEDPHOTONS_H\r
-/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *\r
- * See cxx source for full Copyright notice     */\r
-//______________________________________________________________________________\r
-// Analysis for PHOS Tagged Photons \r
+#ifndef ALIANALYSISTASKTAGGEDPHOTONS_H
+#define ALIANALYSISTASKTAGGEDPHOTONS_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice     */
+//______________________________________________________________________________
+// Analysis for PHOS Tagged Photons 
 // marks photons making pi0 with any other photon
 // and calculates necessary corrections for fake pairs and
 // decay partners escaped acceptance. If MC info is present 
 // fills set of controll histograms.
-//\r
-//*-- Dmitry Blau \r
-//////////////////////////////////////////////////////////////////////////////\r
-\r
-#include "AliAnalysisTaskSE.h"  \r
-\r
-class AliStack ; \r
-class AliESDEvent ; \r
-//class AliAODEvent ; \r
-class TH1D ; \r
-class TH1I ; \r
-class TH2D ;\r
-class TH2F ;\r
-class TH3D ;\r
+//
+//*-- Dmitry Blau 
+//////////////////////////////////////////////////////////////////////////////
+
+#include "AliAnalysisTaskSE.h"  
+
+class AliStack ; 
+class AliESDEvent ; 
+//class AliAODEvent ; 
+class TH1D ; 
+class TH1I ; 
+class TH2D ;
+class TH2F ;
+class TH3D ;
 class AliPHOSGeoUtils;
 class AliEMCALGeoUtils;
-class AliAODPWG4Particle;\r
-\r
-class AliAnalysisTaskTaggedPhotons : public AliAnalysisTaskSE {\r
-\r
-public:\r
-  AliAnalysisTaskTaggedPhotons() ;\r
-  AliAnalysisTaskTaggedPhotons(const char *name) ;\r
-  AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) ;   \r
-  AliAnalysisTaskTaggedPhotons& operator = (const AliAnalysisTaskTaggedPhotons& ap) ;\r
-  virtual ~AliAnalysisTaskTaggedPhotons() ;\r
-   \r
-  virtual void UserCreateOutputObjects(); \r
-  virtual void Init() ; \r
-  virtual void LocalInit() { Init() ; }\r
-  virtual void UserExec(Option_t * opt = "") ;\r
-  virtual void Terminate(Option_t * opt = "") ;\r
-\r
-  void SetDebugLevel(Int_t level) { fDebug = level ; }\r
-  void SetPHOS(Bool_t isPHOS=kTRUE){fPHOS=isPHOS ;} //Which calirimeter to analyze\r
-  \r
-  void SetPhotonId(Float_t threshold) { fPhotonId = threshold ; }\r
-  void SetMinEnergyCut(Float_t threshold) { fMinEnergyCut = threshold; }\r
-\r
-  //Parameterization of pi0 peak position\r
-  void SetPi0MeanParameters(Float_t p0, Float_t p1, Float_t p2, Float_t p3)\r
-      { fPi0MeanP0 = p0; fPi0MeanP1 = p1; fPi0MeanP2 = p2; fPi0MeanP3 = p3;}\r
-  //Parameterization of pi0 peak width    \r
-  void SetPi0SigmaParameters(Float_t p0, Float_t p1, Float_t p2){ fPi0SigmaP0 = p0; fPi0SigmaP1 = p1;fPi0SigmaP2 = p2; }\r
-\r
-protected:\r
-  Float_t GetPhotonId() const { return fPhotonId ; }\r
-  Int_t   GetFiducialArea(Float_t * pos)const ; //what kind of fiducial area hit the photon\r
-  Bool_t  IsSamePi0(const AliAODPWG4Particle *p1, const AliAODPWG4Particle *p2) const; //Check MC genealogy\r
-  Bool_t  IsInPi0Band(Double_t m, Double_t pt)const; //Check if invariant mass is within pi0 peak\r
+class AliAODPWG4Particle;
+
+class AliAnalysisTaskTaggedPhotons : public AliAnalysisTaskSE {
+
+public:
+  AliAnalysisTaskTaggedPhotons() ;
+  AliAnalysisTaskTaggedPhotons(const char *name) ;
+  AliAnalysisTaskTaggedPhotons(const AliAnalysisTaskTaggedPhotons& ap) ;   
+  AliAnalysisTaskTaggedPhotons& operator = (const AliAnalysisTaskTaggedPhotons& ap) ;
+  virtual ~AliAnalysisTaskTaggedPhotons() ;
+   
+  virtual void UserCreateOutputObjects(); 
+  virtual void Init() ; 
+  virtual void LocalInit() { Init() ; }
+  virtual void UserExec(Option_t * opt = "") ;
+  virtual void Terminate(Option_t * opt = "") ;
+
+  void SetDebugLevel(Int_t level) { fDebug = level ; }
+  void SetPHOS(Bool_t isPHOS=kTRUE){fPHOS=isPHOS ;} //Which calirimeter to analyze
+  
+  void SetPhotonId(Float_t threshold) { fPhotonId = threshold ; }
+  void SetMinEnergyCut(Float_t threshold) { fMinEnergyCut = threshold; }
+
+  //Parameterization of pi0 peak position
+  void SetPi0MeanParameters(Float_t p0, Float_t p1, Float_t p2, Float_t p3)
+      { fPi0MeanP0 = p0; fPi0MeanP1 = p1; fPi0MeanP2 = p2; fPi0MeanP3 = p3;}
+  //Parameterization of pi0 peak width    
+  void SetPi0SigmaParameters(Float_t p0, Float_t p1, Float_t p2){ fPi0SigmaP0 = p0; fPi0SigmaP1 = p1;fPi0SigmaP2 = p2; }
+
+protected:
+  Float_t GetPhotonId() const { return fPhotonId ; }
+  Int_t   GetFiducialArea(Float_t * pos)const ; //what kind of fiducial area hit the photon
+  Bool_t  IsSamePi0(const AliAODPWG4Particle *p1, const AliAODPWG4Particle *p2) const; //Check MC genealogy
+  Bool_t  IsInPi0Band(Double_t m, Double_t pt)const; //Check if invariant mass is within pi0 peak
   Bool_t  TestDisp(Double_t l0, Double_t l1, Double_t e)const  ;
   Bool_t  TestTOF(Double_t /*t*/,Double_t /*en*/)const{return kTRUE;} 
   Bool_t  TestCharged(Double_t /*dr*/,Double_t /*en*/)const{return kTRUE;} 
-\r
-private:\r
-\r
-  AliPHOSGeoUtils *fPHOSgeom;    //!PHOS geometry\r
-  AliEMCALGeoUtils *fEMCALgeom;  //!EMCAL geometry\r
-
-  AliStack        *fStack ;      //!Pointer to MC stack\r
-  Bool_t           fPHOS ;       //Choose Calorimeter: PHOS/EMCAL\r
-\r
-  // task parameters\r
-  Float_t   fPhotonId ;          // threshold for photon identification (Bayesian)\r
-  Float_t   fMinEnergyCut;       // min energy of partner photon\r
-  Float_t   fPi0MeanP0;          // Parameterization of pi0 mass:\r
-  Float_t   fPi0MeanP1;          // m_mean_pi0 = p[0] + p[1]*pt + p[2]*pt^2 + p[3]*pt^3\r
-  Float_t   fPi0MeanP2;          // m_mean_pi0 = p[0] + p[1]*pt + p[2]*pt^2 + p[3]*pt^3\r
-  Float_t   fPi0MeanP3;          // m_mean_pi0 = p[0] + p[1]*pt + p[2]*pt^2 + p[3]*pt^3\r
-\r
-  Float_t   fPi0SigmaP0;        // sigma_m_pi0 = sqrt ( p0*p0/x + p1*p1 + p2*p2/x/x)\r
-  Float_t   fPi0SigmaP1;\r      // sigma_m_pi0 = sqrt ( p0*p0/x + p1*p1 + p2*p2/x/x)
-  Float_t   fPi0SigmaP2;\r      // sigma_m_pi0 = sqrt ( p0*p0/x + p1*p1 + p2*p2/x/x)
-\r
-  //Fiducial area parameters\r
+
+private:
+
+  AliPHOSGeoUtils *fPHOSgeom;    //PHOS geometry
+  AliEMCALGeoUtils *fEMCALgeom;  //EMCAL geometry
+
+  AliStack        *fStack ;      //Pointer to MC stack
+  Bool_t           fPHOS ;       //Choose Calorimeter: PHOS/EMCAL
+
+  // task parameters
+  Float_t   fPhotonId ;          // threshold for photon identification (Bayesian)
+  Float_t   fMinEnergyCut;       // min energy of partner photon
+  Float_t   fPi0MeanP0;          // Parameterization of pi0 mass:
+  Float_t   fPi0MeanP1;          // m_mean_pi0 = p[0] + p[1]*pt + p[2]*pt^2 + p[3]*pt^3
+  Float_t   fPi0MeanP2;          // m_mean_pi0 = p[0] + p[1]*pt + p[2]*pt^2 + p[3]*pt^3
+  Float_t   fPi0MeanP3;          // m_mean_pi0 = p[0] + p[1]*pt + p[2]*pt^2 + p[3]*pt^3
+
+  Float_t   fPi0SigmaP0;        // sigma_m_pi0 = sqrt ( p0*p0/x + p1*p1 + p2*p2/x/x)
+  Float_t   fPi0SigmaP1;      // sigma_m_pi0 = sqrt ( p0*p0/x + p1*p1 + p2*p2/x/x)
+  Float_t   fPi0SigmaP2;      // sigma_m_pi0 = sqrt ( p0*p0/x + p1*p1 + p2*p2/x/x)
+
+  //Fiducial area parameters
   Float_t fZmax ;               //Rectangular
   Float_t fZmin ;               //area
   Float_t fPhimax ;             //covered by
   Float_t fPhimin ;             //full calorimeter
-\r
-  TList   * fOutputList ;        //! output data list\r
-  TList   * fEventList ;         //!  event list for mixed InvMass\r
-\r
-  // Histograms\r
-  //Reconstructed spectra\r
-  TH1D    * fhRecAll;                  // Spectrum of all reconstructed particles\r
-  TH1D    * fhRecAllArea1;             // Spectrum of rec particles in Fid. Area 1\r
-  TH1D    * fhRecAllArea2;             // Spectrum of rec particles in Fid. Area 2\r
-  TH1D    * fhRecAllArea3;             // Spectrum of rec particles in Fid. Area 3\r
-\r
-  //Sort registered particles spectra according MC information\r
-  TH1D    * fhRecPhoton;               // Spectrum of rec. with primary==22 and no PID criteria\r
-  TH1D    * fhRecOther;                // Spectrum of rec. with primary!=22 and no PID criteria\r
-  TH1D    * fhRecPhotonPID[4];         // Spectrum of rec. with primary==22 and different PID criteria\r
-  TH1D    * fhRecOtherPID[4];          // Spectrum of rec. with primary!=22 and different PID criteria\r
-  TH1D    * fhRecPhotPi0 ;             // Spectrum of rec. photons from pi0 decays\r
-  TH1D    * fhRecPhotEta ;             // Spectrum of rec. photons from eta decays\r
-  TH1D    * fhRecPhotOmega ;           // Spectrum of rec. photons from omega decays\r
-  TH1D    * fhRecPhotEtapr ;           // Spectrum of rec. photons from eta prime decays\r
-  TH1D    * fhRecPhotConv ;            // Spectrum of rec. photons from conversion\r
-  TH1D    * fhRecPhotHadron ;          // Spectrum of rec. photons from hadron-matter interactions\r
-  TH1D    * fhRecPhotDirect ;          // Spectrum of rec. photons direct or no primary\r
-  TH1D    * fhRecPhotOther ;           // Spectrum of rec. photons from other hadron decays\r
-\r
-  //MC tagging: reasons of partner loss etc.\r
-  TH1D  * fhDecWMCPartner ;         //pi0 decay photon which partner should be registered according to MC\r
-  TH1D  * fhDecWMissedPartnerNotPhoton ; //Decay photon with missed non-photon partner\r
-  TH1D  * fhDecWMissedPartnerAll ;  //Decay photons with partner missed due to some reason (sum of below)\r
-  TH1D  * fhDecWMissedPartnerEmin;  //Decay photons with partner missed due to low energy\r
-  TH1D  * fhDecWMissedPartnerConv;  //Decay photons with partner missed due to conversion\r
-  TH1D  * fhDecWMissedPartnerStack;  //Decay photons with partner not present in Stack\r
-  TH1D  * fhDecWMissedPartnerGeom0; //Decay photons with partner missed due geometry\r
-  TH1D  * fhDecWMissedPartnerGeom1; //Decay photons with partner missed due geometry Fid. area. 1\r
-  TH1D  * fhDecWMissedPartnerGeom2; //Decay photons with partner missed due geometry Fid. area. 2\r
-  TH1D  * fhDecWMissedPartnerGeom3; //Decay photons with partner missed due geometry Fid. area. 3\r
-\r
-  //MC tagging: Decay partners spectra\r
-  TH1D  * fhPartnerMCReg ;       //Spectrum of decay partners which should be registered (MC)\r
-  TH1D  * fhPartnerMissedEmin ;  //Spectrum of decay partners lost due to Emin cut  \r
-  TH1D  * fhPartnerMissedConv ;  //Spectrum of decay partners lost due to conversion\r
-  TH1D  * fhPartnerMissedGeo  ;  //Spectrum of decay partners lost due to acceptance\r
-\r
-  //Tagging\r
-  TH1D  * fhTaggedAll ;      //Spectrum of all tagged photons\r
-  TH1D  * fhTaggedArea1 ;    //Spectrum of all tagged photons Fid. area1\r
-  TH1D  * fhTaggedArea2 ;    //Spectrum of all tagged photons Fid. area2\r
-  TH1D  * fhTaggedArea3 ;    //Spectrum of all tagged photons Fid. area3\r
-  TH1D  * fhTaggedPID[4] ;   //Spectrum of tagged photons for different PID criteria\r
-  TH1D  * fhTaggedMult ;     //Spectrum of multiply tagged photons\r
-\r
-  //Tagging: use MC information if available\r
-  TH1D *  fhTaggedMCTrue ;     //Spectrum of correctly tagged pi0 decay photons\r
-  TH1D *  fhMCMissedTagging ;  //Spectrum of pi0 decay photons missed tagging due to wrong pair mass\r
-  TH1D *  fhMCFakeTagged ;     //Spectrum of photons wrongly tagged according to MC\r
-\r
-  //Invariant mass distributions for fake corrections\r
-  TH2D * fhInvMassReal ;    //Two-photon inv. mass vs first photon pt\r
-  TH2D * fhInvMassMixed ;   //Two-photon inv. mass vs first photon pt\r
-  TH2D * fhMCMissedTaggingMass ; //Inv mass of pairs missed tagging\r
-  \r
-  //Conversion and annihilation radius distributions\r
-  TH1D * fhConversionRadius ;      // Radis of photon production (conversion)\r
-  TH1D * fhInteractionRadius ;     // Radis of photon production (hadron interaction)\r
-\r
-  TH1D    * fhEvents;                  //  number of processed Events\r
\r
-  ClassDef(AliAnalysisTaskTaggedPhotons, 1);   // a PHOS photon analysis task \r
-};\r
-#endif // ALIANALYSISTASKTAGGEDPHOTONS_H\r
+
+  TList   * fOutputList ;        // output data list
+  TList   * fEventList ;         //  event list for mixed InvMass
+
+  // Histograms
+  //Reconstructed spectra
+  TH1D    * fhRecAll;                  // Spectrum of all reconstructed particles
+  TH1D    * fhRecAllArea1;             // Spectrum of rec particles in Fid. Area 1
+  TH1D    * fhRecAllArea2;             // Spectrum of rec particles in Fid. Area 2
+  TH1D    * fhRecAllArea3;             // Spectrum of rec particles in Fid. Area 3
+
+  //Sort registered particles spectra according MC information
+  TH1D    * fhRecPhoton;               // Spectrum of rec. with primary==22 and no PID criteria
+  TH1D    * fhRecOther;                // Spectrum of rec. with primary!=22 and no PID criteria
+  TH1D    * fhRecPhotonPID[4];         // Spectrum of rec. with primary==22 and different PID criteria
+  TH1D    * fhRecOtherPID[4];          // Spectrum of rec. with primary!=22 and different PID criteria
+  TH1D    * fhRecPhotPi0 ;             // Spectrum of rec. photons from pi0 decays
+  TH1D    * fhRecPhotEta ;             // Spectrum of rec. photons from eta decays
+  TH1D    * fhRecPhotOmega ;           // Spectrum of rec. photons from omega decays
+  TH1D    * fhRecPhotEtapr ;           // Spectrum of rec. photons from eta prime decays
+  TH1D    * fhRecPhotConv ;            // Spectrum of rec. photons from conversion
+  TH1D    * fhRecPhotHadron ;          // Spectrum of rec. photons from hadron-matter interactions
+  TH1D    * fhRecPhotDirect ;          // Spectrum of rec. photons direct or no primary
+  TH1D    * fhRecPhotOther ;           // Spectrum of rec. photons from other hadron decays
+
+  //MC tagging: reasons of partner loss etc.
+  TH1D  * fhDecWMCPartner ;         //pi0 decay photon which partner should be registered according to MC
+  TH1D  * fhDecWMissedPartnerNotPhoton ; //Decay photon with missed non-photon partner
+  TH1D  * fhDecWMissedPartnerAll ;  //Decay photons with partner missed due to some reason (sum of below)
+  TH1D  * fhDecWMissedPartnerEmin;  //Decay photons with partner missed due to low energy
+  TH1D  * fhDecWMissedPartnerConv;  //Decay photons with partner missed due to conversion
+  TH1D  * fhDecWMissedPartnerStack;  //Decay photons with partner not present in Stack
+  TH1D  * fhDecWMissedPartnerGeom0; //Decay photons with partner missed due geometry
+  TH1D  * fhDecWMissedPartnerGeom1; //Decay photons with partner missed due geometry Fid. area. 1
+  TH1D  * fhDecWMissedPartnerGeom2; //Decay photons with partner missed due geometry Fid. area. 2
+  TH1D  * fhDecWMissedPartnerGeom3; //Decay photons with partner missed due geometry Fid. area. 3
+
+  //MC tagging: Decay partners spectra
+  TH1D  * fhPartnerMCReg ;       //Spectrum of decay partners which should be registered (MC)
+  TH1D  * fhPartnerMissedEmin ;  //Spectrum of decay partners lost due to Emin cut  
+  TH1D  * fhPartnerMissedConv ;  //Spectrum of decay partners lost due to conversion
+  TH1D  * fhPartnerMissedGeo  ;  //Spectrum of decay partners lost due to acceptance
+
+  //Tagging
+  TH1D  * fhTaggedAll ;      //Spectrum of all tagged photons
+  TH1D  * fhTaggedArea1 ;    //Spectrum of all tagged photons Fid. area1
+  TH1D  * fhTaggedArea2 ;    //Spectrum of all tagged photons Fid. area2
+  TH1D  * fhTaggedArea3 ;    //Spectrum of all tagged photons Fid. area3
+  TH1D  * fhTaggedPID[4] ;   //Spectrum of tagged photons for different PID criteria
+  TH1D  * fhTaggedMult ;     //Spectrum of multiply tagged photons
+
+  //Tagging: use MC information if available
+  TH1D *  fhTaggedMCTrue ;     //Spectrum of correctly tagged pi0 decay photons
+  TH1D *  fhMCMissedTagging ;  //Spectrum of pi0 decay photons missed tagging due to wrong pair mass
+  TH1D *  fhMCFakeTagged ;     //Spectrum of photons wrongly tagged according to MC
+
+  //Invariant mass distributions for fake corrections
+  TH2D * fhInvMassReal ;    //Two-photon inv. mass vs first photon pt
+  TH2D * fhInvMassMixed ;   //Two-photon inv. mass vs first photon pt
+  TH2D * fhMCMissedTaggingMass ; //Inv mass of pairs missed tagging
+  
+  //Conversion and annihilation radius distributions
+  TH1D * fhConversionRadius ;      // Radis of photon production (conversion)
+  TH1D * fhInteractionRadius ;     // Radis of photon production (hadron interaction)
+
+  TH1D    * fhEvents;                  //  number of processed Events
+  ClassDef(AliAnalysisTaskTaggedPhotons, 1);   // a PHOS photon analysis task 
+};
+#endif // ALIANALYSISTASKTAGGEDPHOTONS_H