PWGJE/EMCALJetTasks/AliHadCorrTask.cxx

   1 // $Id$
   2 //
   3 // Hadronic correction task.
   4 //
   5 // Author: R.Reed, C.Loizides, S.Aiola
   6
   7 #include <TChain.h>
   8 #include <TClonesArray.h>
   9 #include <TH1F.h>
  10 #include <TH2F.h>
  11 #include <TList.h>
  12
  13 #include "AliAnalysisManager.h"
  14 #include "AliAODEvent.h"
  15 #include "AliAODCaloCluster.h"
  16 #include "AliESDCaloCluster.h"
  17 #include "AliVTrack.h"
  18 #include "AliPicoTrack.h"
  19 #include "AliVEventHandler.h"
  20 #include "AliEmcalParticle.h"
  21 #include "AliEMCALGeometry.h"
  22
  23 #include "AliHadCorrTask.h"
  24
  25 ClassImp(AliHadCorrTask)
  26
  27 //________________________________________________________________________
  28 AliHadCorrTask::AliHadCorrTask() :
  29   AliAnalysisTaskEmcal("AliHadCorrTask", kFALSE),
  30   fOutCaloName(),
  31   fPhiMatch(0.05),
  32   fEtaMatch(0.025),
  33   fDoTrackClus(0),
  34   fHadCorr(0),
  35   fEexclCell(0),
  36   fEsdMode(kTRUE),
  37   fOutClusters(0),
  38   fHistNclusvsCent(0),
  39   fHistNclusMatchvsCent(0),
  40   fHistEbefore(0),
  41   fHistEafter(0),
  42   fHistEoPCent(0),
  43   fHistNMatchCent(0),
  44   fHistNClusMatchCent(0)
  45 {
  46   // Default constructor.
  47
  48   for(Int_t i=0; i<8; i++) {
  49       fHistEsubPch[i]    = 0;
  50       fHistEsubPchRat[i] = 0;
  51     for(Int_t j=0; j<4; j++) {
  52       fHistNCellsEnergy[i][j] = 0;
  53     }
  54     if (i<4) {
  55       fHistMatchEvsP[i]    = 0;
  56       fHistMatchdRvsEP[i]  = 0;
  57       fHistNMatchEnergy[i] = 0;
  58     }
  59     for(Int_t j=0; j<9; j++) {
  60       for(Int_t k=0; k<2; k++) {
  61         fHistMatchEtaPhi[i][j][k] = 0;
  62       }
  63     }
  64   }
  65 }
  66
  67 //________________________________________________________________________
  68 AliHadCorrTask::AliHadCorrTask(const char *name, Bool_t histo) :
  69   AliAnalysisTaskEmcal(name, histo),
  70   fOutCaloName("CaloClustersCorr"),
  71   fPhiMatch(0.05),
  72   fEtaMatch(0.025),
  73   fDoTrackClus(1),
  74   fHadCorr(0),
  75   fEexclCell(0),
  76   fEsdMode(kTRUE),
  77   fOutClusters(0),
  78   fHistNclusvsCent(0),
  79   fHistNclusMatchvsCent(0),
  80   fHistEbefore(0),
  81   fHistEafter(0),
  82   fHistEoPCent(0),
  83   fHistNMatchCent(0),
  84   fHistNClusMatchCent(0)
  85 {
  86   // Standard constructor.
  87
  88   for(Int_t i=0; i<8; i++) {
  89     if (i<4) {
  90       fHistMatchEvsP[i]   = 0;
  91       fHistMatchdRvsEP[i] = 0;
  92       for(Int_t j=0; j<3; j++) {
  93         fHistNCellsEnergy[i][j] = 0;
  94       }
  95     }
  96     fHistEsubPch[i]    = 0;
  97     fHistEsubPchRat[i] = 0;
  98     for(Int_t j=0; j<9; j++) {
  99       for(Int_t k=0; k<2; k++) {
 100         fHistMatchEtaPhi[i][j][k] = 0;
 101       }
 102     }
 103   }
 104
 105   SetMakeGeneralHistograms(histo);
 106
 107   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
 108 }
 109
 110 //________________________________________________________________________
 111 AliHadCorrTask::~AliHadCorrTask()
 112 {
 113   // Destructor
 114 }
 115
 116 //________________________________________________________________________
 117 Int_t AliHadCorrTask::GetMomBin(Double_t p) const
 118 {
 119   // Get momenum bin.
 120
 121   Int_t pbin=-1;
 122   if (p<0.5)
 123     pbin=0;
 124   else if (p>=0.5 && p<1.0)
 125     pbin=1;
 126   else if (p>=1.0 && p<1.5)
 127     pbin=2;
 128   else if (p>=1.5 && p<2.)
 129     pbin=3;
 130   else if (p>=2. && p<3.)
 131     pbin=4;
 132   else if (p>=3. && p<4.)
 133     pbin=5;
 134   else if (p>=4. && p<5.)
 135     pbin=6;
 136   else if (p>=5. && p<8.)
 137     pbin=7;
 138   else if (p>=8.)
 139     pbin=8;
 140
 141   return pbin;
 142 }
 143
 144 //________________________________________________________________________
 145 Double_t AliHadCorrTask::GetEtaSigma(Int_t pbin) const
 146 {
 147   // Get sigma in eta.
 148
 149   Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
 150   return 2.0*EtaSigma[pbin];
 151 }
 152
 153 //________________________________________________________________________
 154 Double_t AliHadCorrTask::GetPhiMean(Int_t pbin, Int_t centbin) const
 155 {
 156   // Get phi mean.
 157
 158   if (centbin==0) {
 159     Double_t PhiMean[9]={0.036,
 160                          0.021,
 161                          0.0121,
 162                          0.0084,
 163                          0.0060,
 164                          0.0041,
 165                          0.0031,
 166                          0.0022,
 167                          0.001};
 168     return PhiMean[pbin];
 169   } else if (centbin==1) {
 170     Double_t PhiMean[9]={0.036,
 171                          0.021,
 172                          0.0121,
 173                          0.0084,
 174                          0.0060,
 175                          0.0041,
 176                          0.0031,
 177                          0.0022,
 178                          0.001};
 179     return PhiMean[pbin];
 180   } else if (centbin==2) {
 181     Double_t PhiMean[9]={0.036,
 182                          0.021,
 183                          0.0121,
 184                          0.0084,
 185                          0.0060,
 186                          0.0041,
 187                          0.0031,
 188                          0.0022,
 189                          0.001};
 190     return PhiMean[pbin];
 191   } else if (centbin==3) {
 192     Double_t PhiMean[9]={0.036,
 193                          0.021,
 194                          0.0121,
 195                          0.0084,
 196                          0.0060,
 197                          0.0041,
 198                          0.0031,
 199                          0.0022,
 200                          0.001};
 201     return PhiMean[pbin];
 202   } else if (centbin==4) {
 203     Double_t PhiMean[9]={0.036,
 204                          0.021,
 205                          0.0127,
 206                          0.0089,
 207                          0.0068,
 208                          0.0049,
 209                          0.0038,
 210                          0.0028,
 211                          0.0018};
 212     return PhiMean[pbin]*(-1.);
 213   } else if (centbin==5) {
 214     Double_t PhiMean[9]={0.036,
 215                          0.021,
 216                          0.0127,
 217                          0.0089,
 218                          0.0068,
 219                          0.0048,
 220                          0.0038,
 221                          0.0028,
 222                          0.0018};
 223     return PhiMean[pbin]*(-1.);
 224   } else if (centbin==6) {
 225     Double_t PhiMean[9]={0.036,
 226                          0.021,
 227                          0.0127,
 228                          0.0089,
 229                          0.0068,
 230                          0.0045,
 231                          0.0035,
 232                          0.0028,
 233                          0.0018};
 234     return PhiMean[pbin]*(-1.);
 235   } else if (centbin==7) {
 236     Double_t PhiMean[9]={0.036,
 237                          0.021,
 238                          0.0127,
 239                          0.0089,
 240                          0.0068,
 241                          0.0043,
 242                          0.0035,
 243                          0.0028,
 244                          0.0018};
 245     return PhiMean[pbin]*(-1.);
 246   }
 247
 248   return 0;
 249 }
 250
 251 //________________________________________________________________________
 252 Double_t AliHadCorrTask::GetPhiSigma(Int_t pbin, Int_t centbin) const
 253 {
 254   // Get phi sigma.
 255
 256   if (centbin==0) {
 257     Double_t PhiSigma[9]={0.0221,
 258                           0.0128,
 259                           0.0074,
 260                           0.0064,
 261                           0.0059,
 262                           0.0055,
 263                           0.0052,
 264                           0.0049,
 265                           0.0045};
 266     return 2.*PhiSigma[pbin];
 267   } else if (centbin==1) {
 268     Double_t PhiSigma[9]={0.0217,
 269                           0.0120,
 270                           0.0076,
 271                           0.0066,
 272                           0.0062,
 273                           0.0058,
 274                           0.0054,
 275                           0.0054,
 276 0.0045};
 277     return 2.*PhiSigma[pbin];
 278   } else if (centbin==2) {
 279     Double_t PhiSigma[9]={0.0211,
 280                           0.0124,
 281                           0.0080,
 282                           0.0070,
 283                           0.0067,
 284                           0.0061,
 285                           0.0059,
 286                           0.0054,
 287                           0.0047};
 288     return 2.*PhiSigma[pbin];
 289   } else if (centbin==3) {
 290     Double_t PhiSigma[9]={0.0215,
 291                           0.0124,
 292                           0.0082,
 293                           0.0073,
 294                           0.0069,
 295                           0.0064,
 296                           0.0060,
 297                           0.0055,
 298                           0.0047};
 299     return 2.*PhiSigma[pbin];
 300   } else if (centbin==4) {
 301     Double_t PhiSigma[9]={0.0199,
 302                           0.0108,
 303                           0.0072,
 304                           0.0071,
 305                           0.0060,
 306                           0.0055,
 307                           0.0052,
 308                           0.0049,
 309                           0.0045};
 310     return 2.*PhiSigma[pbin];
 311   } else if (centbin==5) {
 312     Double_t PhiSigma[9]={0.0200,
 313                           0.0110,
 314                           0.0074,
 315                           0.0071,
 316                           0.0064,
 317                           0.0059,
 318                           0.0055,
 319                           0.0052,
 320                           0.0045};
 321     return 2.*PhiSigma[pbin];
 322   } else if (centbin==6) {
 323     Double_t PhiSigma[9]={0.0202,
 324                           0.0113,
 325                           0.0077,
 326                           0.0071,
 327                           0.0069,
 328                           0.0064,
 329                           0.0060,
 330                           0.0055,
 331                           0.0050};
 332     return 2.*PhiSigma[pbin];
 333   } else if (centbin==7) {
 334     Double_t PhiSigma[9]={0.0205,
 335                           0.0113,
 336                           0.0080,
 337                           0.0074,
 338                           0.0078,
 339                           0.0067,
 340                           0.0062,
 341                           0.0055,
 342                           0.0050};
 343     return 2.*PhiSigma[pbin];
 344   }
 345
 346   return 0;
 347 }
 348
 349 //________________________________________________________________________
 350 void AliHadCorrTask::UserCreateOutputObjects()
 351 {
 352   // Create my user objects.
 353
 354   if (!fCreateHisto)
 355     return;
 356
 357   AliAnalysisTaskEmcal::UserCreateOutputObjects();
 358
 359   TString name;
 360
 361   for(Int_t icent=0; icent<8; ++icent) {
 362     for(Int_t ipt=0; ipt<9; ++ipt) {
 363       for(Int_t ieta=0; ieta<2; ++ieta) {
 364         name = Form("fHistMatchEtaPhi_%i_%i_%i",icent,ipt,ieta);
 365         fHistMatchEtaPhi[icent][ipt][ieta] = new TH2F(name, name, 200, -0.1, 0.1, 400, -0.2, 0.2);
 366         fOutput->Add(fHistMatchEtaPhi[icent][ipt][ieta]);
 367       }
 368     }
 369
 370     name = Form("fHistEsubPch_%i",icent);
 371     fHistEsubPch[icent]=new TH1F(name, name, 400, 0., 100.);
 372     fOutput->Add(fHistEsubPch[icent]);
 373
 374     name = Form("fHistEsubPchRat_%i",icent);
 375     fHistEsubPchRat[icent]=new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
 376     fOutput->Add(fHistEsubPchRat[icent]);
 377
 378     if (icent<4) {
 379       for(Int_t itrk=0; itrk<4; ++itrk) {
 380         name = Form("fHistNCellsEnergy_%i_%i",icent,itrk);
 381         fHistNCellsEnergy[icent][itrk]  = new TH2F(name, name, 1000, 0, 100, 101, -0.5, 100.5);
 382         fOutput->Add(fHistNCellsEnergy[icent][itrk]);
 383       }
 384
 385       name = Form("fHistMatchEvsP_%i",icent);
 386       fHistMatchEvsP[icent] = new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
 387       fOutput->Add(fHistMatchEvsP[icent]);
 388
 389       name = Form("fHistMatchdRvsEP_%i",icent);
 390       fHistMatchdRvsEP[icent] = new TH2F(name, name, 1000, 0., 1., 1000, 0., 10.);
 391       fOutput->Add(fHistMatchdRvsEP[icent]);
 392
 393       name = Form("fHistNMatchEnergy_%i",icent);
 394       fHistNMatchEnergy[icent]  = new TH2F(name, name, 1000, 0, 100, 101, -0.5, 100.5);
 395       fOutput->Add(fHistNMatchEnergy[icent]);
 396     }
 397   }
 398
 399   fHistNclusvsCent      = new TH1F("Nclusvscent",      "NclusVsCent",      100, 0, 100);
 400   fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent", 100, 0, 100);
 401   fHistEbefore          = new TH1F("Ebefore",          "Ebefore",          100, 0, 100);
 402   fHistEafter           = new TH1F("Eafter",           "Eafter",           100, 0, 100);
 403   fHistEoPCent          = new TH2F("EoPCent",          "EoPCent",          100, 0, 100, 1000, 0,   10);
 404   fHistNMatchCent       = new TH2F("NMatchesCent",     "NMatchesCent",     100, 0, 100, 11, -0.5, 10.5);
 405   fHistNClusMatchCent   = new TH2F("NClusMatchesCent", "NClusMatchesCent", 100, 0, 100, 11, -0.5, 10.5);
 406
 407   fOutput->Add(fHistNclusMatchvsCent);
 408   fOutput->Add(fHistNclusvsCent);
 409   fOutput->Add(fHistEbefore);
 410   fOutput->Add(fHistEafter);
 411   fOutput->Add(fHistEoPCent);
 412   fOutput->Add(fHistNMatchCent);
 413   fOutput->Add(fHistNClusMatchCent);
 414
 415   PostData(1, fOutput);
 416 }
 417
 418 //________________________________________________________________________
 419 void AliHadCorrTask::ExecOnce()
 420 {
 421   // Init the analysis.
 422
 423   if (dynamic_cast<AliAODEvent*>(InputEvent()))
 424     fEsdMode = kFALSE;
 425
 426   if (fEsdMode) { // optimization in case autobranch loading is off
 427     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 428     if (fCaloName == "CaloClusters")
 429       am->LoadBranch("CaloClusters");
 430     if (fTracksName == "Tracks")
 431       am->LoadBranch("Tracks");
 432     am->LoadBranch("Centrality.");
 433   }
 434
 435   if (fEsdMode)
 436     fOutClusters = new TClonesArray("AliESDCaloCluster");
 437   else
 438     fOutClusters = new TClonesArray("AliAODCaloCluster");
 439
 440   fOutClusters->SetName(fOutCaloName);
 441
 442   // post output in event if not yet present
 443   if (!(InputEvent()->FindListObject(fOutCaloName))) {
 444     InputEvent()->AddObject(fOutClusters);
 445   }
 446   else {
 447     AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fOutCaloName.Data()));
 448     return;
 449   }
 450
 451   AliAnalysisTaskEmcal::ExecOnce();
 452 }
 453
 454 //________________________________________________________________________
 455 void AliHadCorrTask::DoTrackLoop()
 456 {
 457   const Int_t Ntracks = fTracks->GetEntries();
 458   for (Int_t iTrk = 0; iTrk < Ntracks; ++iTrk) {
 459     Int_t NmatchClus = 0;
 460
 461     AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iTrk));
 462     if (!emctrack)
 463       continue;
 464     if (!emctrack->IsEMCAL())
 465       continue;
 466
 467     AliVTrack *track = emctrack->GetTrack();
 468     if (!track)
 469       continue;
 470     if (!AcceptTrack(track))
 471       continue;
 472
 473     if (track->GetTrackEtaOnEMCal() < fGeom->GetArm1EtaMin() + fEtaMatch ||
 474         track->GetTrackEtaOnEMCal() > fGeom->GetArm1EtaMax() - fEtaMatch ||
 475         track->GetTrackPhiOnEMCal() < fGeom->GetArm1PhiMin() * TMath::DegToRad() + fPhiMatch ||
 476         track->GetTrackPhiOnEMCal() > fGeom->GetArm1PhiMax() * TMath::DegToRad() - fPhiMatch)
 477       continue;
 478
 479     Int_t Nclus = emctrack->GetNumberOfMatchedObj();
 480
 481     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
 482       AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
 483       if (!emccluster)
 484         continue;
 485
 486       AliVCluster *cluster = emccluster->GetCluster();
 487       if (!cluster)
 488         continue;
 489       if (!AcceptCluster(cluster))
 490         continue;
 491
 492       Double_t etadiff = 999;
 493       Double_t phidiff = 999;
 494       AliPicoTrack::GetEtaPhiDiff(track, cluster, phidiff, etadiff);
 495
 496       if (TMath::Abs(phidiff) < fPhiMatch && TMath::Abs(etadiff) < fEtaMatch)
 497         NmatchClus++;
 498     }
 499     fHistNClusMatchCent->Fill(fCent, NmatchClus);
 500   }
 501 }
 502
 503 //________________________________________________________________________
 504 void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster, Double_t &totalTrkP, Int_t &Nmatches)
 505 {
 506   // Do the loop over matched tracks for cluster emccluster.
 507
 508   AliVCluster *cluster = emccluster->GetCluster();
 509   Int_t iClus = emccluster->IdInCollection();
 510   Double_t energyclus = cluster->E();
 511
 512   // loop over matched tracks
 513   const Int_t Ntrks = emccluster->GetNumberOfMatchedObj();
 514   for (Int_t i = 0; i < Ntrks; ++i) {
 515     Int_t    iTrack = emccluster->GetMatchedObjId(i);
 516     Double_t dR     = emccluster->GetMatchedObjDistance(i);
 517
 518     AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iTrack));
 519     if (!emctrack)
 520       continue;
 521     AliVTrack *track = emctrack->GetTrack();
 522     if (!track)
 523       continue;
 524     if (!AcceptTrack(track))
 525       continue;
 526
 527     // check if track also points to cluster
 528     if (fDoTrackClus && (track->GetEMCALcluster()) != iClus)
 529       continue;
 530
 531     Double_t etadiff = 999;
 532     Double_t phidiff = 999;
 533     AliPicoTrack::GetEtaPhiDiff(track, cluster, phidiff, etadiff);
 534
 535     Double_t mom       = track->P();
 536     Int_t    mombin    = GetMomBin(mom);
 537     Int_t    centbinch = fCentBin;
 538     if (track->Charge()<0)
 539       centbinch += 4;
 540
 541     if (fCreateHisto) {
 542       Int_t etabin = 0;
 543       if(track->Eta() > 0) etabin=1;
 544       fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
 545     }
 546
 547     Double_t etaCut   = 0.0;
 548     Double_t phiCutlo = 0.0;
 549     Double_t phiCuthi = 0.0;
 550
 551     if (fPhiMatch > 0) {
 552       phiCutlo = -fPhiMatch;
 553       phiCuthi = +fPhiMatch;
 554     } else {
 555       phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
 556       phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
 557     }
 558
 559     if (fEtaMatch > 0) {
 560       etaCut = fEtaMatch;
 561     } else {
 562       etaCut = GetEtaSigma(mombin);
 563     }
 564
 565     if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
 566       ++Nmatches;
 567       totalTrkP += track->P();
 568
 569       if (fCreateHisto) {
 570         if (fHadCorr > 1 && mombin > -1) {
 571           fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
 572         }
 573       }
 574     }
 575   }
 576 }
 577
 578 //________________________________________________________________________
 579 Bool_t AliHadCorrTask::Run()
 580 {
 581   // Run the hadronic correction
 582
 583   if (fCreateHisto)
 584     DoTrackLoop();
 585
 586   // delete output
 587   fOutClusters->Delete();
 588
 589    // loop over all clusters
 590   const Int_t Nclus = fCaloClusters->GetEntries();
 591   for (Int_t iClus = 0, clusCount=0; iClus < Nclus; ++iClus) {
 592     AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
 593     if (!emccluster)
 594       continue;
 595
 596     AliVCluster *cluster = emccluster->GetCluster();
 597     if (!cluster)
 598       continue;
 599     if (!AcceptCluster(cluster))
 600       continue;
 601
 602     Double_t energyclus = 0;
 603     if (fCreateHisto)
 604       fHistEbefore->Fill(fCent, cluster->E());
 605
 606     // apply correction / subtraction
 607     if (fHadCorr > 0) {
 608       // to subtract only the closest track set fHadCor to a %
 609       // to subtract all tracks within the cut set fHadCor to %+1
 610       if (fHadCorr > 1)
 611         energyclus = ApplyHadCorrAllTracks(emccluster, fHadCorr - 1);
 612       else
 613         energyclus = ApplyHadCorrOneTrack(emccluster, fHadCorr);
 614     }
 615
 616     if (energyclus < 0)
 617       energyclus = 0;
 618
 619     if (energyclus > 0) { // create corrected cluster
 620       AliVCluster *oc;
 621       if (fEsdMode) {
 622         AliESDCaloCluster *ec = dynamic_cast<AliESDCaloCluster*>(cluster);
 623         if (!ec)
 624           continue;
 625         oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*ec);
 626       } else {
 627         AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
 628         if (!ac)
 629           continue;
 630         oc = new ((*fOutClusters)[clusCount]) AliAODCaloCluster(*ac);
 631       }
 632       oc->SetE(energyclus);
 633
 634       ++clusCount;
 635
 636       if (fCreateHisto)
 637         fHistEafter->Fill(fCent, energyclus);
 638     }
 639   }
 640
 641   return kTRUE;
 642 }
 643
 644 //________________________________________________________________________
 645 Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Double_t hadCorr)
 646 {
 647   // Apply the hadronic correction with one track only.
 648
 649   AliVCluster *cluster = emccluster->GetCluster();
 650   Double_t energyclus  = cluster->E();
 651   Int_t    iMin        = emccluster->GetMatchedObjId();
 652   if (iMin < 0)
 653     return energyclus;
 654
 655   AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
 656   if (!emctrack)
 657     return energyclus;
 658
 659   // check if track also points to cluster
 660   Int_t cid = emctrack->GetMatchedObjId();
 661   if (fDoTrackClus && (cid!=emccluster->IdInCollection()))
 662     return energyclus;
 663
 664   AliVTrack *track = emctrack->GetTrack();
 665   if (!track)
 666     return energyclus;
 667
 668   Double_t mom = track->P();
 669   if (mom == 0)
 670     return energyclus;
 671
 672   Double_t dRmin      = emccluster->GetMatchedObjDistance();
 673   Double_t dEtaMin    = 1e9;
 674   Double_t dPhiMin    = 1e9;
 675   AliPicoTrack::GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
 676
 677   Int_t mombin = GetMomBin(mom);
 678   Int_t centbinch = fCentBin;
 679   if (track->Charge()<0)
 680     centbinch += 4;
 681
 682   // plot some histograms if switched on
 683   if (fCreateHisto) {
 684     Int_t etabin = 0;
 685     if(track->Eta() > 0)
 686       etabin = 1;
 687
 688     fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
 689
 690     if (mom > 0) {
 691       fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
 692       fHistEoPCent->Fill(fCent, energyclus / mom);
 693       fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
 694     }
 695   }
 696
 697   // define eta/phi cuts
 698   Double_t etaCut   = 0.0;
 699   Double_t phiCutlo = 0.0;
 700   Double_t phiCuthi = 0.0;
 701   if (fPhiMatch > 0) {
 702     phiCutlo = -fPhiMatch;
 703     phiCuthi = +fPhiMatch;
 704   }
 705   else {
 706     phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
 707     phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
 708   }
 709   if(fEtaMatch > 0) {
 710     etaCut = fEtaMatch;
 711   }
 712   else {
 713     etaCut = GetEtaSigma(mombin);
 714   }
 715
 716   // apply the correction if the track is in the eta/phi window
 717   if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
 718     energyclus -= hadCorr * mom;
 719   }
 720
 721   return energyclus;
 722 }
 723
 724 //________________________________________________________________________
 725 Double_t AliHadCorrTask::ApplyHadCorrAllTracks(AliEmcalParticle *emccluster, Double_t hadCorr)
 726 {
 727   // Apply the hadronic correction with all tracks.
 728
 729   AliVCluster *cluster = emccluster->GetCluster();
 730
 731   Double_t energyclus = cluster->E();
 732   Double_t cNcells = cluster->GetNCells();
 733
 734   Double_t totalTrkP  = 0.0; // count total track momentum
 735   Int_t    Nmatches   = 0;   // count total number of matches
 736
 737   // do the loop over the matched tracks and get the number of matches and the total momentum
 738   DoMatchedTracksLoop(emccluster, totalTrkP, Nmatches);
 739
 740   Double_t Esub = hadCorr * totalTrkP;
 741
 742   if (Esub > energyclus)
 743     Esub = energyclus;
 744
 745   // applying Peter's proposed algorithm
 746   // never subtract the full energy of the cluster
 747   Double_t clusEexcl = fEexclCell * cNcells;
 748   if (energyclus < clusEexcl)
 749     clusEexcl = energyclus;
 750   if ((energyclus - Esub) < clusEexcl)
 751     Esub = (energyclus - clusEexcl);
 752
 753   Double_t EoP = 0;
 754   if (totalTrkP>0)
 755     EoP = energyclus / totalTrkP;
 756
 757   // plot some histograms if switched on
 758   if (fCreateHisto) {
 759     fHistNclusvsCent->Fill(fCent);
 760     fHistNMatchCent->Fill(fCent, Nmatches);
 761     fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
 762
 763     if(Nmatches > 0)
 764       fHistNclusMatchvsCent->Fill(fCent);
 765
 766     if(Nmatches < 3)
 767       fHistNCellsEnergy[fCentBin][Nmatches]->Fill(energyclus, cNcells);
 768     else
 769       fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
 770
 771     if (EoP > 0) {
 772       fHistEoPCent->Fill(fCent, EoP);
 773       fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
 774     }
 775
 776     if (Nmatches == 1 && totalTrkP > 0) {
 777       Int_t iMin = emccluster->GetMatchedObjId();
 778       AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
 779       AliVTrack *track = emctrack->GetTrack();
 780       Int_t centbinchm = fCentBin;
 781       if (track->Charge()<0)
 782         centbinchm += 4;
 783
 784       fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
 785       fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
 786     }
 787   }
 788
 789   // apply the correction
 790   energyclus -= Esub;
 791
 792   return energyclus;
 793 }
 794