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