PWGJE/EMCALJetTasks/UserTasks/AliAnalysisTaskSOH.cxx

   1 // $Id$
   2 //
   3 // Simulation EMCal task.
   4 //
   5 // Author: Saehanseul Oh
   6
   7 #include "AliAnalysisTaskSOH.h"
   8
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TH3F.h>
  12 #include <THnSparse.h>
  13
  14 #include "TChain.h"
  15 #include "AliAnalysisManager.h"
  16 #include "AliAnalysisTask.h"
  17 #include "AliEMCALRecoUtils.h"
  18 #include "AliEMCALTrack.h"
  19 #include "AliESDCaloCluster.h"
  20 #include "AliESDEvent.h"
  21 #include "AliESDInputHandler.h"
  22 #include "AliESDtrack.h"
  23 #include "AliESDtrackCuts.h"
  24 #include "AliExternalTrackParam.h"
  25 #include "AliLog.h"
  26 #include "AliMCEvent.h"
  27 #include "AliHeader.h"
  28 #include "AliGenCocktailEventHeader.h"
  29
  30
  31 ClassImp(AliAnalysisTaskSOH)
  32
  33 //________________________________________________________________________
  34 AliAnalysisTaskSOH::AliAnalysisTaskSOH() :
  35   AliAnalysisTaskSE(),
  36   fESD(0),
  37   fMC(0),
  38   fZVtxMax(10),
  39   fEsdTrackCuts(0x0),
  40   fHybridTrackCuts1(0x0),
  41   fHybridTrackCuts2(0x0),
  42   fTrackIndices(0x0),
  43   fClusterIndices(0x0),
  44   fClusterArray(0x0),
  45   fMcProcess(kTRUE),
  46   fTrackProcess(kTRUE),
  47   fSFProcess(kFALSE),
  48   fClusterProcess(kFALSE),
  49   fOutputList(0x0),
  50   fHEventStat(0),
  51   fHScaleFactor(0),
  52   fHScaleFactor100HC(0),
  53   fHEOverPVsPt(0x0),
  54   fHEMCalResponsePion(0x0),
  55   fHEMCalResponseElec(0x0),
  56   fHEMCalResponseProton(0x0),
  57   fHEMCalRecdPhidEta(0x0),
  58   fHEMCalRecdPhidEtaP(0x0),
  59   fHEMCalRecdPhidEtaM(0x0),
  60   fHEMCalRecdPhidEta_Truth(0x0),
  61   fHEMCalRecdPhidEtaP_Truth(0x0),
  62   fHEMCalRecdPhidEtaM_Truth(0x0),
  63   fHEMCalRecdPhidEtaposEta(0x0),
  64   fHEMCalRecdPhidEtanegEta(0x0),
  65   fHPhotonEdiff100HC(0x0),
  66   fHPhotonEdiff70HC(0),
  67   fHPhotonEdiff30HC(0),
  68   fHPhotonEdiff0HC(0x0),
  69   fHPhotonEVsClsE(0x0),
  70   fHistEsub1Pch(0x0),
  71   fHistEsub2Pch(0x0),
  72   fHistEsub1PchRat(0x0),
  73   fHistEsub2PchRat(0x0),
  74   fHClsEoverMcE_All(0x0),
  75   fHClsEoverMcE_Photon(0x0),
  76   fHClsEoverMcE_Elec(0x0),
  77   fHClsEoverMcE_Pion(0x0),
  78   fHParGenPion_p(0x0),
  79   fHParGenPion_m(0x0),
  80   fHParGenPion_rmInj_p(0x0),
  81   fHParGenPion_rmInj_m(0x0),
  82   fHDetGenFakePion(0x0),
  83   fHDetRecFakePion(0x0),
  84   fHDetGenSecPion(0x0),
  85   fHDetRecSecPion(0x0)
  86 {
  87   for(Int_t i=0; i<3; i++)
  88   {
  89     fHDetGenPion_p[i]        = 0x0;
  90     fHDetRecPion_p[i]        = 0x0;
  91     fHDetGenPion_m[i]        = 0x0;
  92     fHDetRecPion_m[i]        = 0x0;
  93     fHDetGenPion_rmInj_p[i]  = 0x0;
  94     fHDetRecPion_rmInj_p[i]  = 0x0;
  95     fHDetGenPion_rmInj_m[i]  = 0x0;
  96     fHDetRecPion_rmInj_m[i]  = 0x0;
  97   }
  98   DefineInput (0, TChain::Class());
  99   DefineOutput(1, TList::Class());
 100
 101 }
 102
 103
 104 //________________________________________________________________________
 105 AliAnalysisTaskSOH::AliAnalysisTaskSOH(const char *name) :
 106   AliAnalysisTaskSE(name),
 107   fESD(0),
 108   fMC(0),
 109   fZVtxMax(10),
 110   fEsdTrackCuts(0x0),
 111   fHybridTrackCuts1(0x0),
 112   fHybridTrackCuts2(0x0),
 113   fTrackIndices(0x0),
 114   fClusterIndices(0x0),
 115   fClusterArray(0x0),
 116   fMcProcess(kTRUE),
 117   fTrackProcess(kTRUE),
 118   fSFProcess(kFALSE),
 119   fClusterProcess(kFALSE),
 120   fOutputList(0x0),
 121   fHEventStat(0),
 122   fHScaleFactor(0),
 123   fHScaleFactor100HC(0),
 124   fHEOverPVsPt(0x0),
 125   fHEMCalResponsePion(0x0),
 126   fHEMCalResponseElec(0x0),
 127   fHEMCalResponseProton(0x0),
 128   fHEMCalRecdPhidEta(0x0),
 129   fHEMCalRecdPhidEtaP(0x0),
 130   fHEMCalRecdPhidEtaM(0x0),
 131   fHEMCalRecdPhidEta_Truth(0x0),
 132   fHEMCalRecdPhidEtaP_Truth(0x0),
 133   fHEMCalRecdPhidEtaM_Truth(0x0),
 134   fHEMCalRecdPhidEtaposEta(0x0),
 135   fHEMCalRecdPhidEtanegEta(0x0),
 136   fHPhotonEdiff100HC(0x0),
 137   fHPhotonEdiff70HC(0),
 138   fHPhotonEdiff30HC(0),
 139   fHPhotonEdiff0HC(0x0),
 140   fHPhotonEVsClsE(0x0),
 141   fHistEsub1Pch(0x0),
 142   fHistEsub2Pch(0x0),
 143   fHistEsub1PchRat(0x0),
 144   fHistEsub2PchRat(0x0),
 145   fHClsEoverMcE_All(0x0),
 146   fHClsEoverMcE_Photon(0x0),
 147   fHClsEoverMcE_Elec(0x0),
 148   fHClsEoverMcE_Pion(0x0),
 149   fHParGenPion_p(0x0),
 150   fHParGenPion_m(0x0),
 151   fHParGenPion_rmInj_p(0x0),
 152   fHParGenPion_rmInj_m(0x0),
 153   fHDetGenFakePion(0x0),
 154   fHDetRecFakePion(0x0),
 155   fHDetGenSecPion(0x0),
 156   fHDetRecSecPion(0x0)
 157 {
 158   for(Int_t i=0; i<3; i++)
 159   {
 160     fHDetGenPion_p[i]        = 0x0;
 161     fHDetRecPion_p[i]        = 0x0;
 162     fHDetGenPion_m[i]        = 0x0;
 163     fHDetRecPion_m[i]        = 0x0;
 164     fHDetGenPion_rmInj_p[i]  = 0x0;
 165     fHDetRecPion_rmInj_p[i]  = 0x0;
 166     fHDetGenPion_rmInj_m[i]  = 0x0;
 167     fHDetRecPion_rmInj_m[i]  = 0x0;
 168   }
 169
 170   // Constructor
 171   // Output slot #1 writes into a TH1 container
 172   DefineOutput(1, TList::Class());
 173 }
 174
 175
 176 //________________________________________________________________________
 177 AliAnalysisTaskSOH::~AliAnalysisTaskSOH()
 178 {
 179   // Destructor.
 180
 181   if(fEsdTrackCuts) delete fEsdTrackCuts;
 182   if(fHybridTrackCuts1) delete fHybridTrackCuts1;
 183   if(fHybridTrackCuts2) delete fHybridTrackCuts2;
 184   if(fTrackIndices) delete fTrackIndices;
 185   if(fClusterIndices) delete fClusterIndices;
 186   if(fClusterArray) delete fClusterArray;
 187 }
 188
 189
 190 //________________________________________________________________________
 191 void AliAnalysisTaskSOH::UserCreateOutputObjects()
 192 {
 193   // Create histograms, called once.
 194
 195   OpenFile(1);
 196
 197   fOutputList = new TList();
 198   fOutputList->SetOwner(1);
 199
 200   fHEventStat = new TH1F("fHEventStat","Event statistics for analysis",8,0,8);
 201   fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
 202   fHEventStat->GetXaxis()->SetBinLabel(2,"cluster");
 203   fHEventStat->GetXaxis()->SetBinLabel(3,"good cluster");
 204   fHEventStat->GetXaxis()->SetBinLabel(4,"cls/0-truth");
 205   fHEventStat->GetXaxis()->SetBinLabel(5,"cls/1-truth");
 206   fHEventStat->GetXaxis()->SetBinLabel(6,"cls/2-truth");
 207   fHEventStat->GetXaxis()->SetBinLabel(7,"cls/2-goodtruth");
 208   fHEventStat->GetXaxis()->SetBinLabel(8,"cls/>3-truth");
 209   fOutputList->Add(fHEventStat);
 210
 211
 212   if(fSFProcess)
 213   {
 214     fHScaleFactor = new TH1F("fHScaleFactor", "Scale factor distribution without hadronic correction;Scale factor",100,0,10);
 215     fOutputList->Add(fHScaleFactor);
 216
 217     fHScaleFactor100HC = new TH1F("fHScaleFactor100HC", "Scale factor distribution with 100% hadronic correction;Scale factor",100,0,10);
 218     fOutputList->Add(fHScaleFactor100HC);
 219   }
 220
 221   if(fClusterProcess)
 222   {
 223     fHEOverPVsPt = new TH2F("fHEOverPVsPt", "E/P vs track p_{T}; p_{T} (GeV/c); E/P", 200 , 0, 4, 200, 0, 3.2);
 224     fOutputList->Add(fHEOverPVsPt);
 225
 226     fHEMCalResponsePion = new TH2F("fHEMCalResponsePion", "Pion E/P vs track p_{T}; p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);
 227     fOutputList->Add(fHEMCalResponsePion);
 228
 229     fHEMCalResponseElec = new TH2F("fHEMCalResponseElec", "Electron E/P vs track p_{T};  p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);
 230     fOutputList->Add(fHEMCalResponseElec);
 231
 232     fHEMCalResponseProton = new TH2F("fHEMCalResponseProton", "Proton E/P vs track p_{T};  p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);
 233     fOutputList->Add(fHEMCalResponseProton);
 234
 235     fHEMCalRecdPhidEta = new TH2F("fHEMCalRecdPhidEta","EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
 236     fOutputList->Add(fHEMCalRecdPhidEta);
 237
 238     fHEMCalRecdPhidEtaP = new TH2F("fHEMCalRecdPhidEtaP","EMCAL Charge+ Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
 239     fOutputList->Add(fHEMCalRecdPhidEtaP);
 240
 241     fHEMCalRecdPhidEtaM = new TH2F("fHEMCalRecdPhidEtaM","EMCAL Charge- Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
 242     fOutputList->Add(fHEMCalRecdPhidEtaM);
 243
 244     fHEMCalRecdPhidEta_Truth = new TH2F("fHEMCalRecdPhidEta_Truth","EMCAL Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
 245     fOutputList->Add(fHEMCalRecdPhidEta_Truth);
 246
 247     fHEMCalRecdPhidEtaP_Truth = new TH2F("fHEMCalRecdPhidEtaP_Truth","EMCAL Charge+ Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
 248     fOutputList->Add(fHEMCalRecdPhidEtaP_Truth);
 249
 250     fHEMCalRecdPhidEtaM_Truth = new TH2F("fHEMCalRecdPhidEtaM_Truth","EMCAL Charge- Cluster-Track(Truth matched) #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
 251     fOutputList->Add(fHEMCalRecdPhidEtaM_Truth);
 252
 253     fHEMCalRecdPhidEtaposEta = new TH2F("fHEMCalRecdPhidEtaposEta","(+eta track) EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
 254     fOutputList->Add(fHEMCalRecdPhidEtaposEta);
 255
 256     fHEMCalRecdPhidEtanegEta = new TH2F("fHEMCalRecdPhidEtanegEta","(-eta track) EMCAL Cluster-Track #Delta#phi-#Delta#eta; #Delta#eta; #Delta#phi",1000,-0.1,0.1,1000,-0.5,0.5);
 257     fOutputList->Add(fHEMCalRecdPhidEtanegEta);
 258
 259     fHPhotonEdiff100HC = new TH2F("fHPhotonEdiff100HC","Photon (E_{Truth}- E_{calc,100% HC})/E_{Truth} vs. E_{Truth}; E_{Truth} (GeV); (E_{Truth}- E_{calc,100% HC})/E_{Truth}",1000,0,10,600,-4.9,1.1);
 260     fOutputList->Add(fHPhotonEdiff100HC);
 261
 262     fHPhotonEdiff70HC = new TH2F("fHPhotonEdiff70HC","Photon (E_{Truth}- E_{calc,70% HC})/E_{Truth} vs. E_{Truth}; E_{Truth} (GeV); (E_{Truth}- E_{calc,30% HC})/E_{Truth}",1000,0,10,600,-4.9,1.1);
 263     fOutputList->Add(fHPhotonEdiff70HC);
 264
 265     fHPhotonEdiff30HC = new TH2F("fHPhotonEdiff30HC","Photon (E_{Truth}- E_{calc,30% HC})/E_{Truth} vs. E_{Truth}; E_{Truth} (GeV); (E_{Truth}- E_{calc,30% HC})/E_{Truth}",1000,0,10,600,-4.9,1.1);
 266     fOutputList->Add(fHPhotonEdiff30HC);
 267
 268     fHPhotonEdiff0HC = new TH2F("fHPhotonEdiff0HC","Photon (E_{Truth}- E_{calc,0% HC})/E_{Truth} vs. E_{Truth}; E_{Truth} (GeV); (E_{Truth}- E_{cls})/E_{Truth}",1000,0,10,600,-4.9,1.1);
 269     fOutputList->Add(fHPhotonEdiff0HC);
 270
 271     fHPhotonEVsClsE = new TH2F("fHPhotonEVsClsE","Cluster E vs. photon E_{Truth}; photon E_{Truth} (GeV); Cluster E (GeV)",500,0,5,500,0,5);
 272     fOutputList->Add(fHPhotonEVsClsE);
 273
 274     fHistEsub1Pch =new  TH2F("fHistEsub1Pch", "(subtracted E in 100% HC) vs. total track P, clusters with 1 matching track; total track P (GeV/c); E_{sub}(GeV)" , 1000, 0., 10, 1000, 0., 10.);
 275     fOutputList->Add(fHistEsub1Pch);
 276
 277     fHistEsub2Pch =new  TH2F("fHistEsub2Pch", "(subtracted E in 100% HC) vs. total track P, clusters with 2 matching tracks; total track P (GeV/c); E_{sub}(GeV)" , 1000, 0., 10, 1000, 0., 10.);
 278     fOutputList->Add(fHistEsub2Pch);
 279
 280     fHistEsub1PchRat =new  TH2F("fHistEsub1PchRat", "(subtracted E in 100% HC)/total track P vs. total track P, clusters with 1 matching track; total track P (GeV/c); E_{sub}/P_{tot}" , 1000, 0., 10, 1100, 0., 1.1);
 281     fOutputList->Add(fHistEsub1PchRat);
 282
 283     fHistEsub2PchRat =new  TH2F("fHistEsub2PchRat", "(subtracted E in 100% HC)/total track P vs. total track P, clusters with 2 matching tracks; total track P (GeV/c); E_{sub}/P_{tot}" , 1000, 0., 10, 1100, 0., 1.1);
 284     fOutputList->Add(fHistEsub2PchRat);
 285
 286     Int_t bins[4] = {150, 150, 100, 200};
 287     Double_t xmin[4] = {1.3, -0.8, 0, 0};
 288     Double_t xmax[4] = {3.2, 0.8, 10, 2};
 289
 290     fHClsEoverMcE_All = new THnSparseF("fHClsEoverMcE_All", "Cluster E/MC E, clusters with 1 matching particle; #phi; #eta; E (GeV); ClsE/McE", 4, bins, xmin, xmax);
 291     fOutputList->Add(fHClsEoverMcE_All);
 292
 293     fHClsEoverMcE_Photon = new THnSparseF("fHClsEoverMcE_Photon", "Cluster E/MC E, clusters with 1 matching photon; #phi; #eta; E (GeV); ClsE/McE", 4, bins, xmin, xmax);
 294     fOutputList->Add(fHClsEoverMcE_Photon);
 295
 296     fHClsEoverMcE_Elec = new THnSparseF("fHClsEoverMcE_Elec", "Cluster E/MC E, clusters with 1 matching electron; #phi; #eta; E (GeV); ClsE/McE", 4, bins, xmin, xmax);
 297     fOutputList->Add(fHClsEoverMcE_Elec);
 298
 299     fHClsEoverMcE_Pion = new THnSparseF("fHClsEoverMcE_Pion", "Cluster E/MC E, clusters with 1 matching pion; #phi; #eta; E (GeV); ClsE/McE",4, bins, xmin, xmax);
 300     fOutputList->Add(fHClsEoverMcE_Pion);
 301   }
 302
 303   fHParGenPion_p = new TH3F("fHParGenPion_p","Particle level truth Phi-Eta-p_{T} distribution of #pi+",  500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 304   fHParGenPion_p->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 305   fHParGenPion_p->GetYaxis()->SetTitle("#eta");
 306   fHParGenPion_p->GetZaxis()->SetTitle("#phi");
 307   fOutputList->Add(fHParGenPion_p);
 308
 309   fHParGenPion_m = new TH3F("fHParGenPion_m", "Particle level truth Phi-Eta-p_{T} distribution of all #pi-", 500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 310   fHParGenPion_m->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 311   fHParGenPion_m->GetYaxis()->SetTitle("#eta");
 312   fHParGenPion_m->GetZaxis()->SetTitle("#phi");
 313   fOutputList->Add(fHParGenPion_m);
 314
 315   fHParGenPion_rmInj_p = new TH3F("fHParGenPion_rmInj_p","Particle level truth Phi-Eta-p_{T} distribution of all #pi+ without injected signal",  500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 316   fHParGenPion_rmInj_p->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 317   fHParGenPion_rmInj_p->GetYaxis()->SetTitle("#eta");
 318   fHParGenPion_rmInj_p->GetZaxis()->SetTitle("#phi");
 319   fOutputList->Add(fHParGenPion_rmInj_p);
 320
 321   fHParGenPion_rmInj_m = new TH3F("fHParGenPion_rmInj_m","Particle level truth Phi-Eta-p_{T} distribution of #pi- without injected signal", 500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 322   fHParGenPion_rmInj_m->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 323   fHParGenPion_rmInj_m->GetYaxis()->SetTitle("#eta");
 324   fHParGenPion_rmInj_m->GetZaxis()->SetTitle("#phi");
 325   fOutputList->Add(fHParGenPion_rmInj_m);
 326
 327
 328
 329   const char* trackCut[3] = {"cut1","cut2", "cut3"};
 330   if(fMcProcess && fTrackProcess)
 331   {
 332     //Fake
 333     fHDetGenFakePion = new TH3F("fHDetGenFakePion", "fake charged pion track Phi-Eta-p_{T} distribution",500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 334     fHDetGenFakePion->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 335     fHDetGenFakePion->GetYaxis()->SetTitle("#eta");
 336     fHDetGenFakePion->GetZaxis()->SetTitle("#phi");
 337     fOutputList->Add(fHDetGenFakePion);
 338
 339     fHDetRecFakePion = new TH3F("fHDetRecFakePion", "fake charged pion track Phi-Eta-p_{T} distribution", 500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 340     fHDetRecFakePion->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 341     fHDetRecFakePion->GetYaxis()->SetTitle("#eta");
 342     fHDetRecFakePion->GetZaxis()->SetTitle("#phi");
 343     fOutputList->Add(fHDetRecFakePion);
 344
 345     //Secondary
 346     fHDetGenSecPion = new TH3F("fHDetGenSecPion", "secondary charged pion charged pion track Phi-Eta-p_{T} distribution", 500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 347     fHDetGenSecPion->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 348     fHDetGenSecPion->GetYaxis()->SetTitle("#eta");
 349     fHDetGenSecPion->GetZaxis()->SetTitle("#phi");
 350     fOutputList->Add(fHDetGenSecPion);
 351
 352     fHDetRecSecPion = new TH3F("fHDetRecSecPion", "secondary charged pion charged pion track Phi-Eta-p_{T} distribution", 500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 353     fHDetRecSecPion->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 354     fHDetRecSecPion->GetYaxis()->SetTitle("#eta");
 355     fHDetRecSecPion->GetZaxis()->SetTitle("#phi");
 356     fOutputList->Add(fHDetRecSecPion);
 357
 358     for(Int_t i=0; i<3; i++)
 359     {
 360       // pi+
 361       fHDetGenPion_p[i] = new TH3F(Form("fHDetGenPion_p_%s", trackCut[i]), Form("%s: Detector level truth Phi-Eta-p_{T} distribution of #pi+", trackCut[i]),  500, 0, 100, 60, -1.2, 1.2, 128, 0, 6.4);
 362       fHDetGenPion_p[i]->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 363       fHDetGenPion_p[i]->GetYaxis()->SetTitle("#eta");
 364       fHDetGenPion_p[i]->GetZaxis()->SetTitle("#phi");
 365       fOutputList->Add(fHDetGenPion_p[i]);
 366
 367       fHDetRecPion_p[i] = new TH3F(Form("fHDetRecPion_p_%s", trackCut[i]), Form("%s: Reconstructed track Phi-Eta-p_{T} distribution of all #pi+", trackCut[i]),  500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 368       fHDetRecPion_p[i]->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 369       fHDetRecPion_p[i]->GetYaxis()->SetTitle("#eta");
 370       fHDetRecPion_p[i]->GetZaxis()->SetTitle("#phi");
 371       fOutputList->Add(fHDetRecPion_p[i]);
 372
 373       // pi-
 374       fHDetGenPion_m[i] = new TH3F(Form("fHDetGenPion_m_%s", trackCut[i]), Form("%s: Detector level truth Phi-Eta-p_{T} distribution of #pi-", trackCut[i]),  500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 375       fHDetGenPion_m[i]->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 376       fHDetGenPion_m[i]->GetYaxis()->SetTitle("#eta");
 377       fHDetGenPion_m[i]->GetZaxis()->SetTitle("#phi");
 378       fOutputList->Add(fHDetGenPion_m[i]);
 379
 380       fHDetRecPion_m[i] = new TH3F(Form("fHDetRecPion_m_%s", trackCut[i]), Form("%s: Reconstructed track Phi-Eta-p_{T} distribution of #pi-", trackCut[i]),  500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 381       fHDetRecPion_m[i]->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 382       fHDetRecPion_m[i]->GetYaxis()->SetTitle("#eta");
 383       fHDetRecPion_m[i]->GetZaxis()->SetTitle("#phi");
 384       fOutputList->Add(fHDetRecPion_m[i]);
 385
 386         //pi+ without injected signal
 387       fHDetGenPion_rmInj_p[i] = new TH3F(Form("fHDetGenPion_rmInj_p_%s", trackCut[i]), Form("%s: Detector level truth Phi-Eta-p_{T} distribution of #pi+ without injected signal", trackCut[i]), 500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 388       fHDetGenPion_rmInj_p[i]->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 389       fHDetGenPion_rmInj_p[i]->GetYaxis()->SetTitle("#eta");
 390       fHDetGenPion_rmInj_p[i]->GetZaxis()->SetTitle("#phi");
 391       fOutputList->Add(fHDetGenPion_rmInj_p[i]);
 392
 393       fHDetRecPion_rmInj_p[i] = new TH3F(Form("fHDetRecPion_rmInj_p_%s", trackCut[i]), Form("%s: Reconstructed track Phi-Eta-p_{T} distribution of #pi+ without injected signal", trackCut[i]), 500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 394       fHDetRecPion_rmInj_p[i]->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 395       fHDetRecPion_rmInj_p[i]->GetYaxis()->SetTitle("#eta");
 396       fHDetRecPion_rmInj_p[i]->GetZaxis()->SetTitle("#phi");
 397       fOutputList->Add(fHDetRecPion_rmInj_p[i]);
 398
 399       //pi- charged particle without injected signal
 400       fHDetGenPion_rmInj_m[i] = new TH3F(Form("fHDetGenPion_rmInj_m_%s", trackCut[i]), Form("%s: Detector level truth Phi-Eta-p_{T} distribution of #pi- without injected signal", trackCut[i]), 500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 401       fHDetGenPion_rmInj_m[i]->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 402       fHDetGenPion_rmInj_m[i]->GetYaxis()->SetTitle("#eta");
 403       fHDetGenPion_rmInj_m[i]->GetZaxis()->SetTitle("#phi");
 404       fOutputList->Add(fHDetGenPion_rmInj_m[i]);
 405
 406       fHDetRecPion_rmInj_m[i] = new TH3F(Form("fHDetRecPion_rmInj_m_%s", trackCut[i]), Form("%s: Reconstructed track Phi-Eta-p_{T} distribution of #pi- without injected signal", trackCut[i]), 500, 0, 100, 100, -1.0, 1.0, 120, 0.0,2.*TMath::Pi());
 407       fHDetRecPion_rmInj_m[i]->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 408       fHDetRecPion_rmInj_m[i]->GetYaxis()->SetTitle("#eta");
 409       fHDetRecPion_rmInj_m[i]->GetZaxis()->SetTitle("#phi");
 410       fOutputList->Add(fHDetRecPion_rmInj_m[i]);
 411     }
 412   }
 413
 414   fTrackIndices = new TArrayI();
 415   fClusterIndices = new TArrayI();
 416
 417   fClusterArray = new TObjArray();
 418   fClusterArray->SetOwner(1);
 419
 420   PostData(1, fOutputList);
 421 }
 422
 423 //________________________________________________________________________
 424 void AliAnalysisTaskSOH::UserExec(Option_t *)
 425 {
 426   // Main loop, called for each event.
 427
 428   // Post output data.
 429   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 430   if (!fESD) {
 431     printf("ERROR: fESD not available\n");
 432     return;
 433   }
 434   if (!EsdVertexOk())  return;  // Vetex cut
 435
 436   fMC = MCEvent();
 437   if (!fMC) {
 438     printf("ERROR: fMC not available\n");
 439     return;
 440   }
 441
 442   fHEventStat->Fill(0.5);
 443
 444   if(fTrackIndices)
 445     fTrackIndices->Reset();
 446   if(fClusterIndices)
 447     fClusterIndices->Reset();
 448   if(fClusterArray)
 449     fClusterArray->Delete();
 450
 451   if(fTrackProcess)
 452     ProcessTrack();
 453   if(fClusterProcess)
 454     ProcessCluster();
 455   if(fMcProcess)
 456     ProcessMc();
 457   if(fSFProcess)
 458     ProcessScaleFactor();
 459
 460   PostData(1, fOutputList);
 461 }
 462
 463 //________________________________________________________________________
 464 void  AliAnalysisTaskSOH::ProcessTrack()
 465 {
 466   // Process track.
 467
 468   fTrackIndices->Set(fESD->GetNumberOfTracks());
 469   AliDebug(3,Form("%s:%d Selecting tracks",(char*)__FILE__,__LINE__));
 470
 471   Int_t isMth = 0;
 472   Int_t nTracks = 0;
 473
 474   Float_t ClsPos[3] = {-999,-999,-999};
 475   Double_t emcTrkpos[3] = {-999,-999,-999};
 476
 477   for(Int_t itr=0; itr<fESD->GetNumberOfTracks(); itr++)
 478   {
 479     AliESDtrack *esdtrack = fESD->GetTrack(itr);
 480     if(!esdtrack)continue;
 481     AliESDtrack *newTrack = GetAcceptTrack(esdtrack);
 482     if(!newTrack) continue;
 483     if(newTrack->Pt()<0.15 || TMath::Abs(newTrack->Eta())>0.9) {delete newTrack; continue;}
 484
 485     if(fClusterProcess)
 486     {
 487       Double_t clsE = -1;
 488       Int_t clsIndex = newTrack->GetEMCALcluster();
 489       if(newTrack->GetEMCALcluster()>-1)
 490       {
 491         AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsIndex);
 492         if(IsGoodCluster(cluster))
 493         {
 494           isMth=1;
 495
 496           cluster->GetPosition(ClsPos);
 497           TVector3 VClsPos(ClsPos[0], ClsPos[1], ClsPos[2]);
 498
 499           AliEMCALTrack EMCTrk(*newTrack);
 500           if(!EMCTrk.PropagateToGlobal(ClsPos[0], ClsPos[1], ClsPos[2], 0.0, 0.0)) {continue;}
 501           EMCTrk.GetXYZ(emcTrkpos);
 502           TVector3 VemcTrkPos(emcTrkpos[0],emcTrkpos[1],emcTrkpos[2]);
 503
 504           Double_t dPhi = VClsPos.Phi() - VemcTrkPos.Phi();
 505           if (dPhi < -1*TMath::Pi()) dPhi += (2*TMath::Pi());
 506           else if (dPhi > TMath::Pi()) dPhi -= (2*TMath::Pi());
 507
 508           Double_t dEta = VClsPos.Eta() - VemcTrkPos.Eta();
 509
 510           fHEMCalRecdPhidEta->Fill(dEta, dPhi);
 511
 512           if((newTrack->GetLabel())>-1 && (newTrack->GetLabel()) < fMC->GetNumberOfTracks())
 513           {
 514             AliVParticle *vParticle = fMC->GetTrack(newTrack->GetLabel());
 515             if(IsGoodMcParticle(vParticle, newTrack->GetLabel()))
 516             {
 517               fHEMCalRecdPhidEta_Truth->Fill(dEta, dPhi);
 518               if(vParticle->Charge() > 0) fHEMCalRecdPhidEtaP_Truth->Fill(dEta, dPhi);
 519               if(vParticle->Charge() < 0) fHEMCalRecdPhidEtaM_Truth->Fill(dEta, dPhi);
 520             }
 521           }
 522
 523           if(esdtrack->Charge() > 0) {fHEMCalRecdPhidEtaP->Fill(dEta, dPhi);}
 524           if(esdtrack->Charge() < 0) {fHEMCalRecdPhidEtaM->Fill(dEta, dPhi);}
 525
 526           if(VemcTrkPos.Eta() > 0) fHEMCalRecdPhidEtaposEta->Fill(dEta, dPhi);
 527           if(VemcTrkPos.Eta() < 0) fHEMCalRecdPhidEtanegEta->Fill(dEta, dPhi);
 528
 529           clsE = cluster->E();
 530           if(newTrack->P()>0) fHEOverPVsPt->Fill(newTrack->Pt(),clsE/newTrack->P());
 531         }
 532
 533         Int_t ipart = newTrack->GetLabel();
 534         if(ipart>-1 && ipart<fMC->GetNumberOfTracks())
 535         {
 536           AliVParticle* vParticle = fMC->GetTrack(ipart);
 537           if(isMth && vParticle)
 538           {
 539             if(TMath::Abs(vParticle->PdgCode())==211)
 540             {
 541               fHEMCalResponsePion->Fill(newTrack->Pt(),clsE/newTrack->P());
 542             }
 543             if(TMath::Abs(vParticle->PdgCode())==11)
 544             {
 545               fHEMCalResponseElec->Fill(newTrack->Pt(),clsE/newTrack->P());
 546             }
 547             if(TMath::Abs(vParticle->PdgCode())==2212)
 548             {
 549               fHEMCalResponseProton->Fill(newTrack->Pt(),clsE/newTrack->P());
 550             }
 551           }
 552         }
 553       }
 554     }
 555     if(newTrack) delete newTrack;
 556
 557     // Track Indices
 558     fTrackIndices->AddAt(itr,nTracks);
 559     nTracks++;
 560   }
 561
 562   fTrackIndices->Set(nTracks);
 563 }
 564
 565 //________________________________________________________________________
 566 void AliAnalysisTaskSOH::ProcessCluster()
 567 {
 568   // Process cluster.
 569
 570   Int_t nCluster = 0;
 571   TLorentzVector gamma;
 572   Double_t vertex[3] = {0, 0, 0};
 573   fESD->GetVertex()->GetXYZ(vertex);
 574   const Int_t nCaloClusters = fESD->GetNumberOfCaloClusters();
 575   fClusterIndices->Set(nCaloClusters);
 576   Float_t ClsPos[3] = {-999,-999,-999};
 577
 578   for(Int_t itr=0; itr<nCaloClusters; itr++)
 579   {
 580     fHEventStat->Fill(1.5);
 581     AliESDCaloCluster *cluster = fESD->GetCaloCluster(itr);
 582     if(!IsGoodCluster(cluster)) continue;
 583     cluster->GetMomentum(gamma, vertex);
 584     if (gamma.Pt() < 0.15) continue;
 585     fHEventStat->Fill(2.5);
 586
 587     cluster->GetPosition(ClsPos);
 588     TVector3 VClsPos(ClsPos[0], ClsPos[1], ClsPos[2]);
 589
 590     TArrayI *TrackLabels = cluster->GetTracksMatched();
 591
 592     if(TrackLabels->GetSize() == 1)
 593     {
 594       AliESDtrack *esdtrack = fESD->GetTrack(TrackLabels->operator[](0));
 595       AliESDtrack *newTrack = GetAcceptTrack(esdtrack);
 596       if(newTrack && TMath::Abs(newTrack->Eta())<0.7)
 597       {
 598         Double_t Esub = newTrack->P();
 599         if (Esub > cluster->E()) Esub = cluster->E();
 600         fHistEsub1Pch->Fill(newTrack->P(), Esub);
 601         fHistEsub1PchRat->Fill(newTrack->P(), Esub/newTrack->P());
 602       }
 603     }
 604
 605     if(TrackLabels->GetSize() == 2)
 606     {
 607       AliESDtrack *esdtrack1 = fESD->GetTrack(TrackLabels->operator[](0));
 608       AliESDtrack *esdtrack2 = fESD->GetTrack(TrackLabels->operator[](1));
 609       AliESDtrack *newTrack1 = GetAcceptTrack(esdtrack1);
 610       AliESDtrack *newTrack2 = GetAcceptTrack(esdtrack2);
 611       if(newTrack1 && newTrack2 && TMath::Abs(newTrack1->Eta())<0.7  && TMath::Abs(newTrack2->Eta())<0.7)
 612       {
 613         Double_t Esub = newTrack1->P() + newTrack2->P();
 614         if (Esub > cluster->E()) Esub = cluster->E();
 615         fHistEsub2Pch->Fill(newTrack1->P() + newTrack2->P(), Esub);
 616         fHistEsub2PchRat->Fill(newTrack1->P() + newTrack2->P(), Esub/(newTrack1->P() + newTrack2->P()));
 617       }
 618       else if(newTrack1 && !(newTrack2) && TMath::Abs(newTrack1->Eta())<0.7)
 619       {
 620         Double_t Esub = newTrack1->P();
 621         if (Esub > cluster->E()) Esub = cluster->E();
 622         fHistEsub1Pch->Fill(newTrack1->P(), Esub);
 623         fHistEsub1PchRat->Fill(newTrack1->P(), Esub/newTrack1->P());
 624       }
 625       else if (!(newTrack1) && newTrack2 && TMath::Abs(newTrack2->Eta())<0.7)
 626       {
 627         Double_t Esub = newTrack2->P();
 628         if (Esub > cluster->E()) Esub = cluster->E();
 629         fHistEsub1Pch->Fill(newTrack2->P(), Esub);
 630         fHistEsub1PchRat->Fill(newTrack2->P(), Esub/newTrack2->P());
 631       }
 632       else {;}
 633     }
 634
 635     TArrayI *MCLabels = cluster->GetLabelsArray();
 636
 637     if(MCLabels->GetSize() == 0) fHEventStat->Fill(3.5);
 638     if(MCLabels->GetSize() == 1)
 639     {
 640       fHEventStat->Fill(4.5);
 641       AliVParticle* vParticle1 = fMC->GetTrack(MCLabels->operator[](0));
 642       if(IsGoodMcParticle(vParticle1, MCLabels->operator[](0)))
 643       {
 644         Double_t Entries[4] = {VClsPos.Phi(), VClsPos.Eta(), vParticle1->E(), cluster->E()/vParticle1->E()};
 645         fHClsEoverMcE_All->Fill(Entries);
 646         if(vParticle1->PdgCode() == 22)
 647         {
 648           fHClsEoverMcE_Photon->Fill(Entries);
 649         }
 650         if(TMath::Abs(vParticle1->PdgCode()) == 11)
 651         {
 652           fHClsEoverMcE_Elec->Fill(Entries);
 653         }
 654         if(TMath::Abs(vParticle1->PdgCode()) == 211)
 655         {
 656           fHClsEoverMcE_Pion->Fill(Entries);
 657         }
 658       }
 659     }
 660     if(MCLabels->GetSize() == 2)
 661     {
 662       fHEventStat->Fill(5.5);
 663       AliVParticle* vParticle1 = fMC->GetTrack(MCLabels->operator[](0));
 664       AliVParticle* vParticle2 = fMC->GetTrack(MCLabels->operator[](1));
 665       if(IsGoodMcParticle(vParticle1, MCLabels->operator[](0)) && IsGoodMcParticle(vParticle2, MCLabels->operator[](1)))
 666       {
 667         fHEventStat->Fill(6.5);
 668         if((vParticle1->PdgCode()==22) && (vParticle2->PdgCode()==22)) {;}
 669         else if((vParticle1->PdgCode()!=22) && (vParticle2->PdgCode()!=22)) {;}
 670         else
 671         {
 672           fClusterIndices->AddAt(itr,nCluster);
 673           nCluster++;
 674         }
 675       }
 676     }
 677     if(MCLabels->GetSize() > 2) fHEventStat->Fill(7.5);
 678
 679     AliESDCaloCluster *newCluster = new AliESDCaloCluster(*cluster);
 680
 681     Double_t subE = 0;
 682     TArrayI arrayTrackMatched(fTrackIndices->GetSize());
 683     Int_t nGoodMatch = 0;
 684
 685     for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
 686     {
 687       AliESDtrack *trk = fESD->GetTrack(fTrackIndices->At(j));
 688       if(itr==trk->GetEMCALcluster())
 689       {
 690         arrayTrackMatched[nGoodMatch] = j;
 691         nGoodMatch ++;
 692         subE += trk->P();
 693       }
 694     }
 695
 696     arrayTrackMatched.Set(nGoodMatch);
 697     newCluster->AddTracksMatched(arrayTrackMatched);
 698
 699     Double_t clsE = newCluster->E();
 700     Double_t newE = clsE-subE;
 701     if(newE<0) newE = 0;
 702     newCluster->SetDispersion(newE);
 703     fClusterArray->Add(newCluster);
 704   }
 705
 706   fClusterIndices->Set(nCluster);
 707 }
 708 //________________________________________________________________________
 709 void AliAnalysisTaskSOH::ProcessMc()
 710 {
 711   // Process MC.
 712    for(Int_t i=0; i<fTrackIndices->GetSize(); i++)
 713   {
 714     AliESDtrack *esdtrack = fESD->GetTrack(fTrackIndices->At(i));
 715     if(!esdtrack)continue;
 716     AliESDtrack *newTrack = GetAcceptTrack(esdtrack);
 717     if(!newTrack) continue;
 718     if(newTrack->Pt()<0.15 || TMath::Abs(newTrack->Eta())>0.9) {delete newTrack; continue;}
 719
 720     Int_t index = esdtrack->GetLabel();
 721     if(index < 0)
 722     {
 723       AliVParticle *vParticle1  = (AliVParticle*)fMC->GetTrack(-1*index);
 724       if((TMath::Abs(vParticle1->PdgCode())==211) && IsGoodMcParticle(vParticle1, -1*index))
 725       {
 726         fHDetRecFakePion->Fill(newTrack->Pt(), newTrack->Eta(), newTrack->Phi());
 727         fHDetGenFakePion->Fill(vParticle1->Pt(), vParticle1->Eta(), vParticle1->Phi());
 728       }
 729     }
 730
 731     AliVParticle* vParticle2 = fMC->GetTrack(TMath::Abs(index));
 732     if(!IsGoodMcParticle(vParticle2, TMath::Abs(index)) && (TMath::Abs(vParticle2->PdgCode())==211))
 733     {
 734       fHDetRecSecPion->Fill(newTrack->Pt(), newTrack->Eta(), newTrack->Phi());
 735       fHDetGenSecPion->Fill(vParticle2->Pt(), vParticle2->Eta(), vParticle2->Phi());
 736     }
 737
 738     if(newTrack) delete newTrack;
 739   }
 740
 741   //tracking effciency
 742   AliHeader* header = (AliHeader*) fMC->Header();
 743   if (!header) AliFatal("fInjectedSignals set but no MC header found");
 744
 745   AliGenCocktailEventHeader* cocktailHeader = dynamic_cast<AliGenCocktailEventHeader*> (header->GenEventHeader());
 746   if (!cocktailHeader)
 747   {
 748     header->Dump();
 749     AliFatal("fInjectedSignals set but no MC cocktail header found");
 750   }
 751
 752   AliGenEventHeader* eventHeader = 0;
 753   eventHeader = dynamic_cast<AliGenEventHeader*> (cocktailHeader->GetHeaders()->First());
 754   if (!eventHeader) AliFatal("First event header not found");
 755
 756   for(Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
 757   {
 758     AliVParticle* vParticle = fMC->GetTrack(ipart);
 759     Int_t pdgCode = vParticle->PdgCode();
 760     AliMCParticle *McParticle  = (AliMCParticle*)fMC->GetTrack(ipart);
 761     if(!IsGoodMcParticle(vParticle, ipart)) continue;
 762
 763     if(TMath::Abs(vParticle->Eta())<0.9)
 764     {
 765       if(pdgCode==211)
 766       {
 767         fHParGenPion_p->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
 768         if(McParticle->GetMother() < eventHeader->NProduced()) fHParGenPion_rmInj_p->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
 769       }
 770
 771       else if(pdgCode==-211)
 772       {
 773         fHParGenPion_m->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
 774         if(McParticle->GetMother() < eventHeader->NProduced()) fHParGenPion_rmInj_m->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
 775       }
 776     }
 777
 778     for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
 779     {
 780       AliESDtrack *esdtrack = fESD->GetTrack(fTrackIndices->At(j));
 781       if(!esdtrack)continue;
 782       AliESDtrack *newTrack = GetAcceptTrack(esdtrack);
 783       if(!newTrack) continue;
 784       if(newTrack->Pt()<0.15 || TMath::Abs(newTrack->Eta())>0.9) {delete newTrack; continue;}
 785
 786       Int_t cutType = (Int_t)newTrack->GetTRDQuality();
 787
 788       if(newTrack->GetLabel()==ipart && TMath::Abs(vParticle->Eta())<0.9)
 789       {
 790         if(pdgCode==211)
 791         {
 792           fHDetGenPion_p[cutType]->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
 793           fHDetRecPion_p[cutType]->Fill(newTrack->Pt(), newTrack->Eta(), newTrack->Phi());
 794           if(McParticle->GetMother() < eventHeader->NProduced())
 795           {
 796             fHDetGenPion_rmInj_p[cutType]->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
 797             fHDetRecPion_rmInj_p[cutType]->Fill(newTrack->Pt(), newTrack->Eta(), newTrack->Phi());
 798           }
 799         }
 800         else if(pdgCode==-211)
 801         {
 802           fHDetGenPion_m[cutType]->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
 803           fHDetRecPion_m[cutType]->Fill(newTrack->Pt(), newTrack->Eta(), newTrack->Phi());
 804           if(McParticle->GetMother() < eventHeader->NProduced())
 805           {
 806             fHDetGenPion_rmInj_m[cutType]->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
 807             fHDetRecPion_rmInj_m[cutType]->Fill(newTrack->Pt(), newTrack->Eta(), newTrack->Phi());
 808           }
 809         }
 810
 811
 812     //cluster E vs. truth photon energy
 813         if(fClusterProcess)
 814          {
 815            for(Int_t k=0; k<fClusterIndices->GetSize(); k++)
 816            {
 817              AliESDCaloCluster *cluster = fESD->GetCaloCluster(fClusterIndices->At(k));
 818              Double_t clsE = cluster->E();
 819              TArrayI *MCLabels = cluster->GetLabelsArray();
 820              AliVParticle* vParticle1 = fMC->GetTrack(MCLabels->operator[](0));
 821              AliVParticle* vParticle2 = fMC->GetTrack(MCLabels->operator[](1));
 822
 823              if(vParticle1->PdgCode()==22 && vParticle2 == vParticle)
 824              {
 825                fHPhotonEdiff0HC->Fill(vParticle1->E(), (vParticle1->E() - clsE)/vParticle1->E());
 826                fHPhotonEVsClsE->Fill(vParticle1->E(), clsE);
 827
 828                if((clsE - 0.3*esdtrack->E())<0) fHPhotonEdiff30HC->Fill(vParticle1->E(), 1);
 829                else  fHPhotonEdiff30HC->Fill(vParticle1->E(), (vParticle1->E() + 0.3*esdtrack->E() - clsE)/vParticle1->E());
 830
 831                if((clsE - 0.7*esdtrack->E())<0) fHPhotonEdiff70HC->Fill(vParticle1->E(), 1);
 832                else  fHPhotonEdiff70HC->Fill(vParticle1->E(), (vParticle1->E() + 0.7*esdtrack->E() - clsE)/vParticle1->E());
 833
 834                if((clsE - esdtrack->E())<0) fHPhotonEdiff100HC->Fill(vParticle1->E(), 1);
 835                else  fHPhotonEdiff100HC->Fill(vParticle1->E(), (vParticle1->E() + esdtrack->E() - clsE)/vParticle1->E());
 836                continue;
 837              }
 838              if(vParticle2->PdgCode()==22 && vParticle1 == vParticle)
 839              {
 840                fHPhotonEdiff0HC->Fill(vParticle2->E(), (vParticle2->E() - clsE)/vParticle2->E());
 841                fHPhotonEVsClsE->Fill(vParticle2->E(), clsE);
 842
 843                if((clsE - 0.3*esdtrack->E())<0) fHPhotonEdiff30HC->Fill(vParticle2->E(), 1);
 844                else fHPhotonEdiff30HC->Fill(vParticle2->E(), (vParticle2->E() + 0.3*esdtrack->E() - clsE)/vParticle2->E());
 845
 846                if((clsE - 0.7*esdtrack->E())<0) fHPhotonEdiff70HC->Fill(vParticle2->E(), 1);
 847                else fHPhotonEdiff70HC->Fill(vParticle2->E(), (vParticle2->E() + 0.7*esdtrack->E() - clsE)/vParticle2->E());
 848
 849                if((clsE-esdtrack->E())<0) fHPhotonEdiff100HC->Fill(vParticle2->E(), 1);
 850                else fHPhotonEdiff100HC->Fill(vParticle2->E(), (vParticle2->E() + esdtrack->E() - clsE)/vParticle2->E());
 851              }
 852            }
 853          }
 854         if(newTrack) delete newTrack;
 855         break;
 856       }
 857       if(newTrack) delete newTrack;
 858     }
 859   }
 860 }
 861
 862
 863 //________________________________________________________________________
 864 void AliAnalysisTaskSOH::ProcessScaleFactor()
 865 {
 866   // Scale factor.
 867
 868   const Double_t phiMax = 180 * TMath::DegToRad();
 869   const Double_t phiMin = 80 * TMath::DegToRad();
 870   const Double_t TPCArea= 2*TMath::Pi()*1.8;
 871   const Double_t EMCArea = (phiMax-phiMin)*1.4;
 872
 873   Double_t PtEMC = 0;
 874   Double_t PtTPC = 0;
 875
 876   for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
 877   {
 878     AliESDtrack *trk = fESD->GetTrack(fTrackIndices->At(j));
 879     Double_t eta = trk->Eta();
 880     Double_t phi = trk->Phi();
 881     if(TMath::Abs(eta)<0.9) PtTPC += trk->Pt();
 882     if(TMath::Abs(eta)<0.7  && phi > phiMin && phi < phiMax ) PtEMC += trk->Pt();
 883   }
 884
 885   Double_t EtWithHadCorr = 0;
 886   Double_t EtWithoutHadCorr = 0;
 887   Double_t vertex[3] = {0, 0, 0};
 888   fESD->GetVertex()->GetXYZ(vertex);
 889   TLorentzVector gamma;
 890
 891   for(Int_t i=0; i<fClusterArray->GetEntriesFast(); i++)
 892   {
 893     AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(i);
 894     cluster->GetMomentum(gamma, vertex);
 895     Double_t sinTheta = TMath::Sqrt(1-TMath::Power(gamma.CosTheta(),2));
 896     EtWithoutHadCorr +=  cluster->E() * sinTheta;
 897     EtWithHadCorr += cluster->GetDispersion() * sinTheta;
 898   }
 899
 900   if(PtTPC>0)
 901   {
 902     fHScaleFactor->Fill((PtEMC+EtWithoutHadCorr)/EMCArea * TPCArea/PtTPC);
 903     fHScaleFactor100HC->Fill((PtEMC+EtWithHadCorr)/EMCArea * TPCArea/PtTPC);
 904   }
 905 }
 906
 907 //________________________________________________________________________
 908 AliESDtrack *AliAnalysisTaskSOH::GetAcceptTrack(AliESDtrack *esdtrack)
 909 {
 910   // Get accepted track.
 911   AliESDtrack *newTrack = 0x0;
 912   if(fEsdTrackCuts->AcceptTrack(esdtrack))
 913   {
 914     newTrack = new AliESDtrack(*esdtrack);
 915     newTrack->SetTRDQuality(0);
 916   }
 917   else if(fHybridTrackCuts1->AcceptTrack(esdtrack))
 918   {
 919     if(esdtrack->GetConstrainedParam())
 920     {
 921       newTrack = new AliESDtrack(*esdtrack);
 922       const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
 923       newTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
 924       newTrack->SetTRDQuality(1);
 925     }
 926     else
 927       return 0x0;
 928   }
 929   else if(fHybridTrackCuts2->AcceptTrack(esdtrack))
 930   {
 931     if(esdtrack->GetConstrainedParam())
 932     {
 933       newTrack = new AliESDtrack(*esdtrack);
 934       const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
 935       newTrack->Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
 936       newTrack->SetTRDQuality(2);
 937     }
 938     else
 939       return 0x0;
 940   }
 941   else
 942   {
 943     return 0x0;
 944   }
 945
 946   return newTrack;
 947 }
 948
 949 //________________________________________________________________________
 950 Bool_t AliAnalysisTaskSOH::IsGoodMcParticle(AliVParticle* vParticle, Int_t ipart)
 951 {
 952   // Return true if good MC particle.
 953
 954   if(!vParticle) return kFALSE;
 955   if(!fMC->IsPhysicalPrimary(ipart)) return kFALSE;
 956   if (TMath::Abs(vParticle->Eta())>2) return kFALSE;
 957   if(vParticle->Pt()<0.15) return kFALSE;
 958   return kTRUE;
 959 }
 960
 961 //________________________________________________________________________
 962 Bool_t AliAnalysisTaskSOH::IsGoodCluster(AliESDCaloCluster *cluster)
 963 {
 964   // Return true if good cluster.
 965
 966   if(!cluster) return kFALSE;
 967   if (!cluster->IsEMCAL()) return kFALSE;
 968   return kTRUE;
 969 }
 970
 971 //________________________________________________________________________
 972 void AliAnalysisTaskSOH::Terminate(Option_t *)
 973 {
 974   // Terminate analysis.
 975 }
 976
 977 //________________________________________________________________________
 978 Bool_t AliAnalysisTaskSOH::EsdVertexOk() const
 979 {
 980   // Modified from AliAnalyseLeadingTrackUE::VertexSelection()
 981
 982   const AliESDVertex* vtx = fESD->GetPrimaryVertex();
 983   if (!vtx) return kFALSE;
 984   Int_t nContributors = vtx->GetNContributors();
 985   Double_t zVertex    = vtx->GetZ();
 986   if( nContributors < 1 || TMath::Abs(zVertex) > fZVtxMax ) return kFALSE;
 987   return kTRUE;
 988 }