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