PWGGA/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
  13 #include "AliAnalysisManager.h"
  14 #include "AliAnalysisTask.h"
  15 #include "AliEMCALRecoUtils.h"
  16 #include "AliEMCALTrack.h"
  17 #include "AliESDCaloCluster.h"
  18 #include "AliESDEvent.h"
  19 #include "AliESDInputHandler.h"
  20 #include "AliESDtrack.h"
  21 #include "AliESDtrackCuts.h"
  22 #include "AliExternalTrackParam.h"
  23 #include "AliLog.h"
  24 #include "AliMCEvent.h"
  25
  26 ClassImp(AliAnalysisTaskSOH)
  27
  28 //________________________________________________________________________
  29 AliAnalysisTaskSOH::AliAnalysisTaskSOH() :
  30   AliAnalysisTaskSE(),
  31   fESD(0),
  32   fMC(0),
  33   fEsdTrackCuts(0x0),
  34   fHybridTrackCuts1(0x0),
  35   fHybridTrackCuts2(0x0),
  36   fTrackIndices(0x0),
  37   fClusterIndices(0x0),
  38   fClusterArray(0x0),
  39   fOutputList(0x0),
  40   fHEventStat(0),
  41   fHTrkEffParGenPtEtaPhi(0),
  42   fHTrkEffDetGenPtEtaPhi(0),
  43   fHTrkEffDetRecPtEtaPhi(0),
  44   fHTrkEffDetRecFakePtEtaPhi(0),
  45   fHScaleFactor(0),
  46   fHScaleFactor100HC(0),
  47   fHEOverPVsPt(0x0),
  48   fHEMCalResponsePion(0x0),
  49   fHEMCalResponseElec(0x0),
  50   fHEMCalResponseProton(0x0),
  51   fHEMCalRecdPhidEta(0x0),
  52   fHEMCalRecdPhidEtaP(0x0),
  53   fHEMCalRecdPhidEtaM(0x0),
  54   fHEMCalRecdPhidEta_Truth(0x0),
  55   fHEMCalRecdPhidEtaP_Truth(0x0),
  56   fHEMCalRecdPhidEtaM_Truth(0x0),
  57   fHEMCalRecdPhidEtaposEta(0x0),
  58   fHEMCalRecdPhidEtanegEta(0x0),
  59   fHPhotonEdiff100HC(0x0),
  60   fHPhotonEdiff70HC(0),
  61   fHPhotonEdiff30HC(0),
  62   fHPhotonEdiff0HC(0x0),
  63   fHPhotonEVsClsE(0x0),
  64   fHistEsub1Pch(0x0),
  65   fHistEsub2Pch(0x0),
  66   fHistEsub1PchRat(0x0),
  67   fHistEsub2PchRat(0x0)
  68 {
  69   // Constructor
  70
  71 }
  72
  73
  74 //________________________________________________________________________
  75 AliAnalysisTaskSOH::AliAnalysisTaskSOH(const char *name) :
  76   AliAnalysisTaskSE(name),
  77   fESD(0),
  78   fMC(0),
  79   fEsdTrackCuts(0x0),
  80   fHybridTrackCuts1(0x0),
  81   fHybridTrackCuts2(0x0),
  82   fTrackIndices(0x0),
  83   fClusterIndices(0x0),
  84   fClusterArray(0x0),
  85   fOutputList(0x0),
  86   fHEventStat(0),
  87   fHTrkEffParGenPtEtaPhi(0),
  88   fHTrkEffDetGenPtEtaPhi(0),
  89   fHTrkEffDetRecPtEtaPhi(0),
  90   fHTrkEffDetRecFakePtEtaPhi(0),
  91   fHScaleFactor(0),
  92   fHScaleFactor100HC(0),
  93   fHEOverPVsPt(0x0),
  94   fHEMCalResponsePion(0x0),
  95   fHEMCalResponseElec(0x0),
  96   fHEMCalResponseProton(0x0),
  97   fHEMCalRecdPhidEta(0x0),
  98   fHEMCalRecdPhidEtaP(0x0),
  99   fHEMCalRecdPhidEtaM(0x0),
 100   fHEMCalRecdPhidEta_Truth(0x0),
 101   fHEMCalRecdPhidEtaP_Truth(0x0),
 102   fHEMCalRecdPhidEtaM_Truth(0x0),
 103   fHEMCalRecdPhidEtaposEta(0x0),
 104   fHEMCalRecdPhidEtanegEta(0x0),
 105   fHPhotonEdiff100HC(0x0),
 106   fHPhotonEdiff70HC(0),
 107   fHPhotonEdiff30HC(0),
 108   fHPhotonEdiff0HC(0x0),
 109   fHPhotonEVsClsE(0x0),
 110   fHistEsub1Pch(0x0),
 111   fHistEsub2Pch(0x0),
 112   fHistEsub1PchRat(0x0),
 113   fHistEsub2PchRat(0x0)
 114 {
 115   // Constructor
 116
 117   // Output slot #1 writes into a TH1 container
 118   DefineOutput(1, TList::Class());
 119 }
 120
 121
 122 //________________________________________________________________________
 123 AliAnalysisTaskSOH::~AliAnalysisTaskSOH()
 124 {
 125   // Destructor.
 126
 127   delete fEsdTrackCuts;
 128   delete fHybridTrackCuts1;
 129   delete fHybridTrackCuts2;
 130   delete fTrackIndices;
 131   delete fClusterIndices;
 132   delete fClusterArray;
 133 }
 134
 135
 136 //________________________________________________________________________
 137 void AliAnalysisTaskSOH::UserCreateOutputObjects()
 138 {
 139   // Create histograms, called once.
 140
 141   OpenFile(1);
 142
 143   fOutputList = new TList();
 144   fOutputList->SetOwner(1);
 145
 146   fHEventStat = new TH1F("fHEventStat","Event statistics for analysis",8,0,8);
 147   fHEventStat->GetXaxis()->SetBinLabel(1,"Event");
 148   fHEventStat->GetXaxis()->SetBinLabel(2,"cluster");
 149   fHEventStat->GetXaxis()->SetBinLabel(3,"good cluster");
 150   fHEventStat->GetXaxis()->SetBinLabel(4,"cls/0-truth");
 151   fHEventStat->GetXaxis()->SetBinLabel(5,"cls/1-truth");
 152   fHEventStat->GetXaxis()->SetBinLabel(6,"cls/2-truth");
 153   fHEventStat->GetXaxis()->SetBinLabel(7,"cls/2-goodtruth");
 154   fHEventStat->GetXaxis()->SetBinLabel(8,"cls/>3-truth");
 155   fOutputList->Add(fHEventStat);
 156
 157   fHTrkEffParGenPtEtaPhi = new TH3F("fHTrkEffParGenPtEtaPhi","Particle level truth Phi-Eta-p_{T} distribution of primary charged pions", 500, 0, 50, 60, -1.5, 1.5, 128, 0, 6.4);
 158   fHTrkEffParGenPtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 159   fHTrkEffParGenPtEtaPhi->GetYaxis()->SetTitle("#eta");
 160   fHTrkEffParGenPtEtaPhi->GetZaxis()->SetTitle("#phi");
 161   fOutputList->Add(fHTrkEffParGenPtEtaPhi);
 162
 163   fHTrkEffDetGenPtEtaPhi = new TH3F("fHTrkEffDetGenPtEtaPhi","Detector level truth Phi-Eta-p_{T} distribution of primary charged pions", 500, 0, 50, 60, -1.5, 1.5, 128, 0, 6.4);
 164   fHTrkEffDetGenPtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 165   fHTrkEffDetGenPtEtaPhi->GetYaxis()->SetTitle("#eta");
 166   fHTrkEffDetGenPtEtaPhi->GetZaxis()->SetTitle("#phi");
 167   fOutputList->Add(fHTrkEffDetGenPtEtaPhi);
 168
 169   fHTrkEffDetRecPtEtaPhi = new TH3F("fHTrkEffDetRecPtEtaPhi","Reconstructed track Phi-Eta-p_{T} distribution of primary charged pions", 500, 0, 50, 60, -1.5, 1.5, 128, 0, 6.4);
 170   fHTrkEffDetRecPtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 171   fHTrkEffDetRecPtEtaPhi->GetYaxis()->SetTitle("#eta");
 172   fHTrkEffDetRecPtEtaPhi->GetZaxis()->SetTitle("#phi");
 173   fOutputList->Add(fHTrkEffDetRecPtEtaPhi);
 174
 175   fHTrkEffDetRecFakePtEtaPhi = new TH3F("fHTrkEffDetRecFakePtEtaPhi","Reconstructed fake track Phi-Eta-p_{T} distribution of pions", 500, 0, 50, 60, -1.5, 1.5, 128, 0, 6.4);
 176   fHTrkEffDetRecFakePtEtaPhi->GetXaxis()->SetTitle("p_{T}^{gen} (GeV/c)");
 177   fHTrkEffDetRecFakePtEtaPhi->GetYaxis()->SetTitle("#eta");
 178   fHTrkEffDetRecFakePtEtaPhi->GetZaxis()->SetTitle("#phi");
 179   fOutputList->Add(fHTrkEffDetRecFakePtEtaPhi);
 180
 181   fHScaleFactor = new TH1F("fHScaleFactor", "Scale factor distribution without hadronic correction;Scale factor",100,0,10);
 182   fOutputList->Add(fHScaleFactor);
 183
 184   fHScaleFactor100HC = new TH1F("fHScaleFactor100HC", "Scale factor distribution with 100% hadronic correction;Scale factor",100,0,10);
 185   fOutputList->Add(fHScaleFactor100HC);
 186
 187   fHEOverPVsPt = new TH2F("fHEOverPVsPt", "E/P vs track p_{T}; p_{T} (GeV/c); E/P", 200 , 0, 4, 200, 0, 3.2);
 188   fOutputList->Add(fHEOverPVsPt);
 189
 190   fHEMCalResponsePion = new TH2F("fHEMCalResponsePion", "Pion E/P vs track p_{T}; p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);
 191   fOutputList->Add(fHEMCalResponsePion);
 192
 193   fHEMCalResponseElec = new TH2F("fHEMCalResponseElec", "Electron E/P vs track p_{T};  p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);
 194   fOutputList->Add(fHEMCalResponseElec);
 195
 196   fHEMCalResponseProton = new TH2F("fHEMCalResponseProton", "Proton E/P vs track p_{T};  p_{T} (GeV/c); E/P", 100 , 0, 4, 100, 0, 3.2);
 197   fOutputList->Add(fHEMCalResponseProton);
 198
 199   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);
 200   fOutputList->Add(fHEMCalRecdPhidEta);
 201
 202   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);
 203   fOutputList->Add(fHEMCalRecdPhidEtaP);
 204
 205   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);
 206   fOutputList->Add(fHEMCalRecdPhidEtaM);
 207
 208   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);
 209   fOutputList->Add(fHEMCalRecdPhidEta_Truth);
 210
 211   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);
 212   fOutputList->Add(fHEMCalRecdPhidEtaP_Truth);
 213
 214   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);
 215   fOutputList->Add(fHEMCalRecdPhidEtaM_Truth);
 216
 217   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);
 218   fOutputList->Add(fHEMCalRecdPhidEtaposEta);
 219
 220   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);
 221   fOutputList->Add(fHEMCalRecdPhidEtanegEta);
 222
 223   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);
 224   fOutputList->Add(fHPhotonEdiff100HC);
 225
 226   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);
 227   fOutputList->Add(fHPhotonEdiff70HC);
 228
 229   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);
 230   fOutputList->Add(fHPhotonEdiff30HC);
 231
 232   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);
 233   fOutputList->Add(fHPhotonEdiff0HC);
 234
 235   fHPhotonEVsClsE = new TH2F("fHPhotonEVsClsE","Cluster E vs. photon E_{Truth}; photon E_{Truth} (GeV); Cluster E (GeV)",500,0,5,500,0,5);
 236   fOutputList->Add(fHPhotonEVsClsE);
 237
 238   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.);
 239   fOutputList->Add(fHistEsub1Pch);
 240
 241   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.);
 242   fOutputList->Add(fHistEsub2Pch);
 243
 244   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);
 245   fOutputList->Add(fHistEsub1PchRat);
 246
 247   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);
 248   fOutputList->Add(fHistEsub2PchRat);
 249
 250   fTrackIndices = new TArrayI();
 251   fClusterIndices = new TArrayI();
 252
 253   fClusterArray = new TObjArray();
 254   fClusterArray->SetOwner(1);
 255
 256   PostData(1, fOutputList);
 257 }
 258
 259 //________________________________________________________________________
 260 void AliAnalysisTaskSOH::UserExec(Option_t *)
 261 {
 262   // Main loop, called for each event.
 263
 264   // Post output data.
 265   fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 266   if (!fESD) {
 267     printf("ERROR: fESD not available\n");
 268     return;
 269   }
 270
 271   fMC = MCEvent();
 272   if (!fMC) {
 273     printf("ERROR: fMC not available\n");
 274     return;
 275   }
 276
 277   fHEventStat->Fill(0.5);
 278
 279   if(fTrackIndices)
 280     fTrackIndices->Reset();
 281   if(fClusterIndices)
 282     fClusterIndices->Reset();
 283   if(fClusterArray)
 284     fClusterArray->Delete();
 285
 286   ProcessTrack();
 287   ProcessCluster();
 288   ProcessMc();
 289   ProcessScaleFactor();
 290
 291   PostData(1, fOutputList);
 292 }
 293
 294 //________________________________________________________________________
 295 void  AliAnalysisTaskSOH::ProcessTrack()
 296 {
 297   // Process track.
 298
 299   fTrackIndices->Set(fESD->GetNumberOfTracks());
 300   AliDebug(3,Form("%s:%d Selecting tracks",(char*)__FILE__,__LINE__));
 301
 302   Int_t isMth = 0;
 303   Int_t nTracks = 0;
 304
 305   Float_t ClsPos[3] = {-999,-999,-999};
 306   Double_t emcTrkpos[3] = {-999,-999,-999};
 307
 308   for(Int_t itr=0; itr<fESD->GetNumberOfTracks(); itr++)
 309   {
 310     AliESDtrack *esdtrack = fESD->GetTrack(itr);
 311     if(!esdtrack)continue;
 312     AliESDtrack *newTrack = GetAcceptTrack(esdtrack);
 313     if(!newTrack) continue;
 314
 315     Double_t clsE = -1;
 316     Int_t clsIndex = newTrack->GetEMCALcluster();
 317     if(newTrack->GetEMCALcluster()>-1)
 318     {
 319       AliESDCaloCluster *cluster = fESD->GetCaloCluster(clsIndex);
 320       if(IsGoodCluster(cluster))
 321       {
 322         isMth=1;
 323
 324         cluster->GetPosition(ClsPos);
 325         TVector3 VClsPos(ClsPos[0], ClsPos[1], ClsPos[2]);
 326
 327         AliEMCALTrack EMCTrk(*newTrack);
 328         if(!EMCTrk.PropagateToGlobal(ClsPos[0], ClsPos[1], ClsPos[2], 0.0, 0.0)) {
 329           continue;
 330         }
 331         EMCTrk.GetXYZ(emcTrkpos);
 332         TVector3 VemcTrkPos(emcTrkpos[0],emcTrkpos[1],emcTrkpos[2]);
 333
 334         Double_t dPhi = VClsPos.Phi() - VemcTrkPos.Phi();
 335         if (dPhi < -1*TMath::Pi()) dPhi += (2*TMath::Pi());
 336         else if (dPhi > TMath::Pi()) dPhi -= (2*TMath::Pi());
 337
 338         Double_t dEta = VClsPos.Eta() - VemcTrkPos.Eta();
 339
 340         fHEMCalRecdPhidEta->Fill(dEta, dPhi);
 341
 342         if((newTrack->GetLabel())>-1 && (newTrack->GetLabel()) < fMC->GetNumberOfTracks())
 343         {
 344           AliVParticle *vParticle = fMC->GetTrack(newTrack->GetLabel());
 345           if(IsGoodMcParticle(vParticle, newTrack->GetLabel()))
 346           {
 347             fHEMCalRecdPhidEta_Truth->Fill(dEta, dPhi);
 348             if(vParticle->Charge() > 0) fHEMCalRecdPhidEtaP_Truth->Fill(dEta, dPhi);
 349             if(vParticle->Charge() < 0) fHEMCalRecdPhidEtaM_Truth->Fill(dEta, dPhi);
 350           }
 351         }
 352
 353         if(esdtrack->Charge() > 0) {fHEMCalRecdPhidEtaP->Fill(dEta, dPhi);}
 354         if(esdtrack->Charge() < 0) {fHEMCalRecdPhidEtaM->Fill(dEta, dPhi);}
 355
 356         if(VemcTrkPos.Eta() > 0) fHEMCalRecdPhidEtaposEta->Fill(dEta, dPhi);
 357         if(VemcTrkPos.Eta() < 0) fHEMCalRecdPhidEtanegEta->Fill(dEta, dPhi);
 358
 359         clsE = cluster->E();
 360         if(newTrack->P()>0) fHEOverPVsPt->Fill(newTrack->Pt(),clsE/newTrack->P());
 361       }
 362
 363       Int_t ipart = newTrack->GetLabel();
 364       if(ipart>-1 && ipart<fMC->GetNumberOfTracks())
 365       {
 366         AliVParticle* vParticle = fMC->GetTrack(ipart);
 367         if(isMth && vParticle)
 368         {
 369           if(TMath::Abs(vParticle->PdgCode())==211)
 370           {
 371             fHEMCalResponsePion->Fill(newTrack->Pt(),clsE/newTrack->P());
 372           }
 373           if(TMath::Abs(vParticle->PdgCode())==11)
 374           {
 375             fHEMCalResponseElec->Fill(newTrack->Pt(),clsE/newTrack->P());
 376           }
 377           if(TMath::Abs(vParticle->PdgCode())==2212)
 378           {
 379             fHEMCalResponseProton->Fill(newTrack->Pt(),clsE/newTrack->P());
 380           }
 381         }
 382       }
 383     }
 384
 385  // fake and secondary tracks
 386     if(newTrack->GetLabel()<0 && newTrack->GetPID()==2) fHTrkEffDetRecFakePtEtaPhi->Fill(newTrack->Pt(), newTrack->Eta(), newTrack->Phi());
 387
 388  // Track Indices
 389     fTrackIndices->AddAt(itr,nTracks);
 390     nTracks++;
 391   }
 392
 393   fTrackIndices->Set(nTracks);
 394 }
 395
 396 //________________________________________________________________________
 397 void AliAnalysisTaskSOH::ProcessCluster()
 398 {
 399   // Process cluster.
 400
 401   Int_t nCluster = 0;
 402   TLorentzVector gamma;
 403   Double_t vertex[3] = {0, 0, 0};
 404   fESD->GetVertex()->GetXYZ(vertex);
 405   const Int_t nCaloClusters = fESD->GetNumberOfCaloClusters();
 406   fClusterIndices->Set(nCaloClusters);
 407
 408   for(Int_t itr=0; itr<nCaloClusters; itr++)
 409   {
 410     fHEventStat->Fill(1.5);
 411     AliESDCaloCluster *cluster = fESD->GetCaloCluster(itr);
 412     if(!IsGoodCluster(cluster)) continue;
 413     cluster->GetMomentum(gamma, vertex);
 414     if (gamma.Pt() < 0.15) continue;
 415     fHEventStat->Fill(2.5);
 416
 417     TArrayI *TrackLabels = cluster->GetTracksMatched();
 418
 419     if(TrackLabels->GetSize() == 1)
 420     {
 421       AliESDtrack *esdtrack = fESD->GetTrack(TrackLabels->GetAt(0));
 422       AliESDtrack *newTrack = GetAcceptTrack(esdtrack);
 423       if(newTrack && TMath::Abs(newTrack->Eta())<0.7)
 424       {
 425         Double_t Esub = newTrack->P();
 426         if (Esub > cluster->E()) Esub = cluster->E();
 427         fHistEsub1Pch->Fill(newTrack->P(), Esub);
 428         fHistEsub1PchRat->Fill(newTrack->P(), Esub/newTrack->P());
 429       }
 430     }
 431
 432     if(TrackLabels->GetSize() == 2)
 433     {
 434       AliESDtrack *esdtrack1 = fESD->GetTrack(TrackLabels->GetAt(0));
 435       AliESDtrack *esdtrack2 = fESD->GetTrack(TrackLabels->GetAt(1));
 436       AliESDtrack *newTrack1 = GetAcceptTrack(esdtrack1);
 437       AliESDtrack *newTrack2 = GetAcceptTrack(esdtrack2);
 438       if(newTrack1 && newTrack2 && TMath::Abs(newTrack1->Eta())<0.7  && TMath::Abs(newTrack2->Eta())<0.7)
 439       {
 440         Double_t Esub = newTrack1->P() + newTrack2->P();
 441         if (Esub > cluster->E()) Esub = cluster->E();
 442         fHistEsub2Pch->Fill(newTrack1->P() + newTrack2->P(), Esub);
 443         fHistEsub2PchRat->Fill(newTrack1->P() + newTrack2->P(), Esub/(newTrack1->P() + newTrack2->P()));
 444       }
 445       else if(newTrack1 && !(newTrack2) && TMath::Abs(newTrack1->Eta())<0.7)
 446       {
 447         Double_t Esub = newTrack1->P();
 448         if (Esub > cluster->E()) Esub = cluster->E();
 449         fHistEsub1Pch->Fill(newTrack1->P(), Esub);
 450         fHistEsub1PchRat->Fill(newTrack1->P(), Esub/newTrack1->P());
 451       }
 452       else if (!(newTrack1) && newTrack2 && TMath::Abs(newTrack2->Eta())<0.7)
 453       {
 454         Double_t Esub = newTrack2->P();
 455         if (Esub > cluster->E()) Esub = cluster->E();
 456         fHistEsub1Pch->Fill(newTrack2->P(), Esub);
 457         fHistEsub1PchRat->Fill(newTrack2->P(), Esub/newTrack2->P());
 458       }
 459       else {;}
 460     }
 461
 462     TArrayI *MCLabels = cluster->GetLabelsArray();
 463
 464     if(MCLabels->GetSize() == 0) fHEventStat->Fill(3.5);
 465     if(MCLabels->GetSize() == 1) fHEventStat->Fill(4.5);
 466     if(MCLabels->GetSize() == 2)
 467     {
 468       fHEventStat->Fill(5.5);
 469       AliVParticle* vParticle1 = fMC->GetTrack(MCLabels->GetAt(0));
 470       AliVParticle* vParticle2 = fMC->GetTrack(MCLabels->GetAt(1));
 471       if(IsGoodMcParticle(vParticle1, MCLabels->GetAt(0)) && IsGoodMcParticle(vParticle2, MCLabels->GetAt(1)))
 472       {
 473         fHEventStat->Fill(6.5);
 474         if((vParticle1->PdgCode()==22) && (vParticle2->PdgCode()==22)) {;}
 475         else if((vParticle1->PdgCode()!=22) && (vParticle2->PdgCode()!=22)) {;}
 476         else
 477         {
 478           fClusterIndices->AddAt(itr,nCluster);
 479           nCluster++;
 480         }
 481       }
 482     }
 483     if(MCLabels->GetSize() > 2) fHEventStat->Fill(7.5);
 484
 485     AliESDCaloCluster *newCluster = new AliESDCaloCluster(*cluster);
 486
 487     Double_t subE = 0;
 488     TArrayI arrayTrackMatched(fTrackIndices->GetSize());
 489     Int_t nGoodMatch = 0;
 490
 491     for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
 492     {
 493       AliESDtrack *trk = fESD->GetTrack(fTrackIndices->At(j));
 494       if(itr==trk->GetEMCALcluster())
 495       {
 496         arrayTrackMatched[nGoodMatch] = j;
 497         nGoodMatch ++;
 498         subE += trk->P();
 499       }
 500     }
 501
 502     arrayTrackMatched.Set(nGoodMatch);
 503     newCluster->AddTracksMatched(arrayTrackMatched);
 504
 505     Double_t clsE = newCluster->E();
 506     Double_t newE = clsE-subE;
 507     if(newE<0) newE = 0;
 508     newCluster->SetDispersion(newE);
 509     fClusterArray->Add(newCluster);
 510   }
 511
 512   fClusterIndices->Set(nCluster);
 513 }
 514 //________________________________________________________________________
 515 void AliAnalysisTaskSOH::ProcessMc()
 516 {
 517   // Process MC.
 518
 519   for(Int_t ipart=0; ipart<fMC->GetNumberOfTracks(); ipart++)
 520   {
 521     AliVParticle* vParticle = fMC->GetTrack(ipart);
 522     if(!IsGoodMcParticle(vParticle, ipart)) continue;
 523     Int_t pdgCode = vParticle->PdgCode();
 524
 525    //tracking effciency
 526     if(TMath::Abs(vParticle->Eta())<0.9)
 527     {
 528       if(TMath::Abs(pdgCode==211)) fHTrkEffParGenPtEtaPhi->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
 529       for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
 530       {
 531         AliESDtrack *esdtrack = fESD->GetTrack(fTrackIndices->At(j));
 532         if(esdtrack && esdtrack->GetLabel()==ipart)
 533         {
 534           if(TMath::Abs(pdgCode==211))
 535           {
 536             fHTrkEffDetGenPtEtaPhi->Fill(vParticle->Pt(), vParticle->Eta(), vParticle->Phi());
 537             fHTrkEffDetRecPtEtaPhi->Fill(esdtrack->Pt(), esdtrack->Eta(), esdtrack->Phi());
 538           }
 539
 540     //cluster E vs. truth photon energy
 541           for(Int_t k=0; k<fClusterIndices->GetSize(); k++)
 542           {
 543             AliESDCaloCluster *cluster = fESD->GetCaloCluster(fClusterIndices->At(k));
 544             Double_t clsE = cluster->E();
 545             TArrayI *MCLabels = cluster->GetLabelsArray();
 546             AliVParticle* vParticle1 = fMC->GetTrack(MCLabels->GetAt(0));
 547             AliVParticle* vParticle2 = fMC->GetTrack(MCLabels->GetAt(1));
 548
 549             if(vParticle1->PdgCode()==22 && vParticle2 == vParticle)
 550             {
 551               fHPhotonEdiff0HC->Fill(vParticle1->E(), (vParticle1->E() - clsE)/vParticle1->E());
 552               fHPhotonEVsClsE->Fill(vParticle1->E(), clsE);
 553
 554               if((clsE - 0.3*esdtrack->E())<0) fHPhotonEdiff30HC->Fill(vParticle1->E(), 1);
 555               else  fHPhotonEdiff30HC->Fill(vParticle1->E(), (vParticle1->E() + 0.3*esdtrack->E() - clsE)/vParticle1->E());
 556
 557               if((clsE - 0.7*esdtrack->E())<0) fHPhotonEdiff70HC->Fill(vParticle1->E(), 1);
 558               else  fHPhotonEdiff70HC->Fill(vParticle1->E(), (vParticle1->E() + 0.7*esdtrack->E() - clsE)/vParticle1->E());
 559
 560               if((clsE - esdtrack->E())<0) fHPhotonEdiff100HC->Fill(vParticle1->E(), 1);
 561               else  fHPhotonEdiff100HC->Fill(vParticle1->E(), (vParticle1->E() + esdtrack->E() - clsE)/vParticle1->E());
 562               continue;
 563             }
 564             if(vParticle2->PdgCode()==22 && vParticle1 == vParticle)
 565             {
 566               fHPhotonEdiff0HC->Fill(vParticle2->E(), (vParticle2->E() - clsE)/vParticle2->E());
 567               fHPhotonEVsClsE->Fill(vParticle2->E(), clsE);
 568
 569               if((clsE - 0.3*esdtrack->E())<0) fHPhotonEdiff30HC->Fill(vParticle2->E(), 1);
 570               else fHPhotonEdiff30HC->Fill(vParticle2->E(), (vParticle2->E() + 0.3*esdtrack->E() - clsE)/vParticle2->E());
 571
 572               if((clsE - 0.7*esdtrack->E())<0) fHPhotonEdiff70HC->Fill(vParticle2->E(), 1);
 573               else fHPhotonEdiff70HC->Fill(vParticle2->E(), (vParticle2->E() + 0.7*esdtrack->E() - clsE)/vParticle2->E());
 574
 575               if((clsE-esdtrack->E())<0) fHPhotonEdiff100HC->Fill(vParticle2->E(), 1);
 576               else fHPhotonEdiff100HC->Fill(vParticle2->E(), (vParticle2->E() + esdtrack->E() - clsE)/vParticle2->E());
 577             }
 578           }
 579           break;
 580         }
 581       }
 582     }
 583   }
 584 }
 585
 586 //________________________________________________________________________
 587 void AliAnalysisTaskSOH::ProcessScaleFactor()
 588 {
 589   // Scale factor.
 590
 591   const Double_t phiMax = 180 * TMath::DegToRad();
 592   const Double_t phiMin = 80 * TMath::DegToRad();
 593   const Double_t TPCArea= 2*TMath::Pi()*1.8;
 594   const Double_t EMCArea = (phiMax-phiMin)*1.4;
 595
 596   Double_t PtEMC = 0;
 597   Double_t PtTPC = 0;
 598
 599   for(Int_t j=0; j<fTrackIndices->GetSize(); j++)
 600   {
 601     AliESDtrack *trk = fESD->GetTrack(fTrackIndices->At(j));
 602     Double_t eta = trk->Eta();
 603     Double_t phi = trk->Phi();
 604     if(TMath::Abs(eta)<0.9) PtTPC += trk->Pt();
 605     if(TMath::Abs(eta)<0.7  && phi > phiMin && phi < phiMax ) PtEMC += trk->Pt();
 606   }
 607
 608   Double_t EtWithHadCorr = 0;
 609   Double_t EtWithoutHadCorr = 0;
 610   Double_t vertex[3] = {0, 0, 0};
 611   fESD->GetVertex()->GetXYZ(vertex);
 612   TLorentzVector gamma;
 613
 614   for(Int_t i=0; i<fClusterArray->GetEntriesFast(); i++)
 615   {
 616     AliESDCaloCluster *cluster = (AliESDCaloCluster*)fClusterArray->At(i);
 617     cluster->GetMomentum(gamma, vertex);
 618     Double_t sinTheta = TMath::Sqrt(1-TMath::Power(gamma.CosTheta(),2));
 619     EtWithoutHadCorr +=  cluster->E() * sinTheta;
 620     EtWithHadCorr += cluster->GetDispersion() * sinTheta;
 621   }
 622
 623   if(PtTPC>0)
 624   {
 625     fHScaleFactor->Fill((PtEMC+EtWithoutHadCorr)/EMCArea * TPCArea/PtTPC);
 626     fHScaleFactor100HC->Fill((PtEMC+EtWithHadCorr)/EMCArea * TPCArea/PtTPC);
 627   }
 628 }
 629
 630 //________________________________________________________________________
 631 AliESDtrack *AliAnalysisTaskSOH::GetAcceptTrack(AliESDtrack *esdtrack)
 632 {
 633   // Get accepted track.
 634
 635   static AliESDtrack newTrack;
 636   if(fEsdTrackCuts->AcceptTrack(esdtrack))
 637   {
 638     esdtrack->Copy(newTrack);
 639     newTrack.SetTRDQuality(0);
 640   }
 641   else if(fHybridTrackCuts1->AcceptTrack(esdtrack))
 642   {
 643     if(esdtrack->GetConstrainedParam())
 644     {
 645       esdtrack->Copy(newTrack);
 646       const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
 647       newTrack.Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
 648       newTrack.SetTRDQuality(1);
 649     }
 650     else
 651       return 0x0;
 652   }
 653   else if(fHybridTrackCuts2->AcceptTrack(esdtrack))
 654   {
 655     if(esdtrack->GetConstrainedParam())
 656     {
 657       esdtrack->Copy(newTrack);
 658       const AliExternalTrackParam* constrainParam = esdtrack->GetConstrainedParam();
 659       newTrack.Set(constrainParam->GetX(),constrainParam->GetAlpha(),constrainParam->GetParameter(),constrainParam->GetCovariance());
 660       newTrack.SetTRDQuality(2);
 661     }
 662     else
 663       return 0x0;
 664   }
 665   else
 666   {
 667     return 0x0;
 668   }
 669
 670   return &newTrack;
 671 }
 672
 673 //________________________________________________________________________
 674 Bool_t AliAnalysisTaskSOH::IsGoodMcParticle(AliVParticle* vParticle, Int_t ipart)
 675 {
 676   // Return true if good MC particle.
 677
 678   if(!vParticle) return kFALSE;
 679   if(!fMC->IsPhysicalPrimary(ipart)) return kFALSE;
 680   if (TMath::Abs(vParticle->Eta())>2) return kFALSE;
 681   if(vParticle->Pt()<0.15) return kFALSE;
 682   return kTRUE;
 683 }
 684
 685 //________________________________________________________________________
 686 Bool_t AliAnalysisTaskSOH::IsGoodCluster(AliESDCaloCluster *cluster)
 687 {
 688   // Return true if good cluster.
 689
 690   if(!cluster) return kFALSE;
 691   if (!cluster->IsEMCAL()) return kFALSE;
 692   return kTRUE;
 693 }
 694
 695 //________________________________________________________________________
 696 void AliAnalysisTaskSOH::Terminate(Option_t *)
 697 {
 698   // Terminate analysis.
 699 }