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