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   // Do base class initializations and if it fails -> bail out
 425   if (!fInitialized)
 426     AliAnalysisTaskEmcal::ExecOnce();
 427   if (!fInitialized)
 428     return;
 429
 430   // Init the analysis.
 431
 432   AliAnalysisTaskEmcal::ExecOnce();
 433
 434   if (!fInitialized)
 435     return;
 436
 437   if (dynamic_cast<AliAODEvent*>(InputEvent()))
 438     fEsdMode = kFALSE;
 439
 440   if (fEsdMode) { // optimization in case autobranch loading is off
 441     AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 442     if (fCaloName == "CaloClusters")
 443       am->LoadBranch("CaloClusters");
 444     if (fTracksName == "Tracks")
 445       am->LoadBranch("Tracks");
 446     am->LoadBranch("Centrality.");
 447   }
 448
 449   if (fEsdMode)
 450     fOutClusters = new TClonesArray("AliESDCaloCluster");
 451   else
 452     fOutClusters = new TClonesArray("AliAODCaloCluster");
 453
 454   fOutClusters->SetName(fOutCaloName);
 455
 456   // post output in event if not yet present
 457   if (!(InputEvent()->FindListObject(fOutCaloName))) {
 458     InputEvent()->AddObject(fOutClusters);
 459   }
 460   else {
 461     fInitialized = kFALSE;
 462     AliFatal(Form("%s: Container with same name %s already present. Aborting", GetName(), fOutCaloName.Data()));
 463     fInitialized = kFALSE;
 464     return;
 465   }
 466 }
 467
 468 //________________________________________________________________________
 469 void AliHadCorrTask::DoTrackLoop()
 470 {
 471   const Int_t Ntracks = fTracks->GetEntries();
 472   for (Int_t iTrk = 0; iTrk < Ntracks; ++iTrk) {
 473     Int_t NmatchClus = 0;
 474
 475     AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iTrk));
 476     if (!emctrack)
 477       continue;
 478     if (!emctrack->IsEMCAL())
 479       continue;
 480
 481     AliVTrack *track = emctrack->GetTrack();
 482     if (!track)
 483       continue;
 484     if (!AcceptTrack(track))
 485       continue;
 486
 487     if (track->GetTrackEtaOnEMCal() < fGeom->GetArm1EtaMin() + fEtaMatch ||
 488         track->GetTrackEtaOnEMCal() > fGeom->GetArm1EtaMax() - fEtaMatch ||
 489         track->GetTrackPhiOnEMCal() < fGeom->GetArm1PhiMin() * TMath::DegToRad() + fPhiMatch ||
 490         track->GetTrackPhiOnEMCal() > fGeom->GetArm1PhiMax() * TMath::DegToRad() - fPhiMatch)
 491       continue;
 492
 493     Int_t Nclus = emctrack->GetNumberOfMatchedObj();
 494
 495     for (Int_t iClus = 0; iClus < Nclus; ++iClus) {
 496       AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
 497       if (!emccluster)
 498         continue;
 499
 500       AliVCluster *cluster = emccluster->GetCluster();
 501       if (!cluster)
 502         continue;
 503       if (!AcceptCluster(cluster))
 504         continue;
 505
 506       Double_t etadiff = 999;
 507       Double_t phidiff = 999;
 508       AliPicoTrack::GetEtaPhiDiff(track, cluster, phidiff, etadiff);
 509
 510       if (TMath::Abs(phidiff) < fPhiMatch && TMath::Abs(etadiff) < fEtaMatch)
 511         NmatchClus++;
 512     }
 513     fHistNClusMatchCent->Fill(fCent, NmatchClus);
 514   }
 515 }
 516
 517 //________________________________________________________________________
 518 void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster, Double_t &totalTrkP, Int_t &Nmatches)
 519 {
 520   // Do the loop over matched tracks for cluster emccluster.
 521
 522   AliVCluster *cluster = emccluster->GetCluster();
 523   Int_t iClus = emccluster->IdInCollection();
 524   Double_t energyclus = cluster->E();
 525
 526   // loop over matched tracks
 527   const Int_t Ntrks = emccluster->GetNumberOfMatchedObj();
 528   for (Int_t i = 0; i < Ntrks; ++i) {
 529     Int_t    iTrack = emccluster->GetMatchedObjId(i);
 530     Double_t dR     = emccluster->GetMatchedObjDistance(i);
 531
 532     AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iTrack));
 533     if (!emctrack)
 534       continue;
 535     AliVTrack *track = emctrack->GetTrack();
 536     if (!track)
 537       continue;
 538     if (!AcceptTrack(track))
 539       continue;
 540
 541     // check if track also points to cluster
 542     if (fDoTrackClus && (track->GetEMCALcluster()) != iClus)
 543       continue;
 544
 545     Double_t etadiff = 999;
 546     Double_t phidiff = 999;
 547     AliPicoTrack::GetEtaPhiDiff(track, cluster, phidiff, etadiff);
 548
 549     Double_t mom       = track->P();
 550     Int_t    mombin    = GetMomBin(mom);
 551     Int_t    centbinch = fCentBin;
 552     if (track->Charge()<0)
 553       centbinch += 4;
 554
 555     if (fCreateHisto) {
 556       Int_t etabin = 0;
 557       if(track->Eta() > 0) etabin=1;
 558       fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
 559     }
 560
 561     Double_t etaCut   = 0.0;
 562     Double_t phiCutlo = 0.0;
 563     Double_t phiCuthi = 0.0;
 564
 565     if (fPhiMatch > 0) {
 566       phiCutlo = -fPhiMatch;
 567       phiCuthi = +fPhiMatch;
 568     } else {
 569       phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
 570       phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
 571     }
 572
 573     if (fEtaMatch > 0) {
 574       etaCut = fEtaMatch;
 575     } else {
 576       etaCut = GetEtaSigma(mombin);
 577     }
 578
 579     if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
 580       ++Nmatches;
 581       totalTrkP += track->P();
 582
 583       if (fCreateHisto) {
 584         if (fHadCorr > 1 && mombin > -1) {
 585           fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
 586         }
 587       }
 588     }
 589   }
 590 }
 591
 592 //________________________________________________________________________
 593 Bool_t AliHadCorrTask::Run()
 594 {
 595   // Run the hadronic correction
 596
 597   if (fCreateHisto)
 598     DoTrackLoop();
 599
 600   // delete output
 601   fOutClusters->Delete();
 602
 603    // loop over all clusters
 604   const Int_t Nclus = fCaloClusters->GetEntries();
 605   for (Int_t iClus = 0, clusCount=0; iClus < Nclus; ++iClus) {
 606     AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
 607     if (!emccluster)
 608       continue;
 609
 610     AliVCluster *cluster = emccluster->GetCluster();
 611     if (!cluster)
 612       continue;
 613     if (!AcceptCluster(cluster))
 614       continue;
 615
 616     Double_t energyclus = 0;
 617     if (fCreateHisto)
 618       fHistEbefore->Fill(fCent, cluster->E());
 619
 620     // apply correction / subtraction
 621     if (fHadCorr > 0) {
 622       // to subtract only the closest track set fHadCor to a %
 623       // to subtract all tracks within the cut set fHadCor to %+1
 624       if (fHadCorr > 1)
 625         energyclus = ApplyHadCorrAllTracks(emccluster, fHadCorr - 1);
 626       else
 627         energyclus = ApplyHadCorrOneTrack(emccluster, fHadCorr);
 628     }
 629
 630     if (energyclus < 0)
 631       energyclus = 0;
 632
 633     if (energyclus > 0) { // create corrected cluster
 634       AliVCluster *oc;
 635       if (fEsdMode) {
 636         AliESDCaloCluster *ec = dynamic_cast<AliESDCaloCluster*>(cluster);
 637         if (!ec)
 638           continue;
 639         oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*ec);
 640       } else {
 641         AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
 642         if (!ac)
 643           continue;
 644         oc = new ((*fOutClusters)[clusCount]) AliAODCaloCluster(*ac);
 645       }
 646       oc->SetE(energyclus);
 647
 648       ++clusCount;
 649
 650       if (fCreateHisto)
 651         fHistEafter->Fill(fCent, energyclus);
 652     }
 653   }
 654
 655   return kTRUE;
 656 }
 657
 658 //________________________________________________________________________
 659 Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Double_t hadCorr)
 660 {
 661   // Apply the hadronic correction with one track only.
 662
 663   AliVCluster *cluster = emccluster->GetCluster();
 664   Double_t energyclus  = cluster->E();
 665   Int_t    iMin        = emccluster->GetMatchedObjId();
 666   if (iMin < 0)
 667     return energyclus;
 668
 669   AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
 670   if (!emctrack)
 671     return energyclus;
 672
 673   // check if track also points to cluster
 674   Int_t cid = emctrack->GetMatchedObjId();
 675   if (fDoTrackClus && (cid!=emccluster->IdInCollection()))
 676     return energyclus;
 677
 678   AliVTrack *track = emctrack->GetTrack();
 679   if (!track)
 680     return energyclus;
 681
 682   Double_t mom = track->P();
 683   if (mom == 0)
 684     return energyclus;
 685
 686   Double_t dRmin      = emccluster->GetMatchedObjDistance();
 687   Double_t dEtaMin    = 1e9;
 688   Double_t dPhiMin    = 1e9;
 689   AliPicoTrack::GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
 690
 691   Int_t mombin = GetMomBin(mom);
 692   Int_t centbinch = fCentBin;
 693   if (track->Charge()<0)
 694     centbinch += 4;
 695
 696   // plot some histograms if switched on
 697   if (fCreateHisto) {
 698     Int_t etabin = 0;
 699     if(track->Eta() > 0)
 700       etabin = 1;
 701
 702     fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
 703
 704     if (mom > 0) {
 705       fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
 706       fHistEoPCent->Fill(fCent, energyclus / mom);
 707       fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
 708     }
 709   }
 710
 711   // define eta/phi cuts
 712   Double_t etaCut   = 0.0;
 713   Double_t phiCutlo = 0.0;
 714   Double_t phiCuthi = 0.0;
 715   if (fPhiMatch > 0) {
 716     phiCutlo = -fPhiMatch;
 717     phiCuthi = +fPhiMatch;
 718   }
 719   else {
 720     phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
 721     phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
 722   }
 723   if(fEtaMatch > 0) {
 724     etaCut = fEtaMatch;
 725   }
 726   else {
 727     etaCut = GetEtaSigma(mombin);
 728   }
 729
 730   // apply the correction if the track is in the eta/phi window
 731   if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
 732     energyclus -= hadCorr * mom;
 733   }
 734
 735   return energyclus;
 736 }
 737
 738 //________________________________________________________________________
 739 Double_t AliHadCorrTask::ApplyHadCorrAllTracks(AliEmcalParticle *emccluster, Double_t hadCorr)
 740 {
 741   // Apply the hadronic correction with all tracks.
 742
 743   AliVCluster *cluster = emccluster->GetCluster();
 744
 745   Double_t energyclus = cluster->E();
 746   Double_t cNcells = cluster->GetNCells();
 747
 748   Double_t totalTrkP  = 0.0; // count total track momentum
 749   Int_t    Nmatches   = 0;   // count total number of matches
 750
 751   // do the loop over the matched tracks and get the number of matches and the total momentum
 752   DoMatchedTracksLoop(emccluster, totalTrkP, Nmatches);
 753
 754   Double_t Esub = hadCorr * totalTrkP;
 755
 756   if (Esub > energyclus)
 757     Esub = energyclus;
 758
 759   // applying Peter's proposed algorithm
 760   // never subtract the full energy of the cluster
 761   Double_t clusEexcl = fEexclCell * cNcells;
 762   if (energyclus < clusEexcl)
 763     clusEexcl = energyclus;
 764   if ((energyclus - Esub) < clusEexcl)
 765     Esub = (energyclus - clusEexcl);
 766
 767   Double_t EoP = 0;
 768   if (totalTrkP>0)
 769     EoP = energyclus / totalTrkP;
 770
 771   // plot some histograms if switched on
 772   if (fCreateHisto) {
 773     fHistNclusvsCent->Fill(fCent);
 774     fHistNMatchCent->Fill(fCent, Nmatches);
 775     fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
 776
 777     if(Nmatches > 0)
 778       fHistNclusMatchvsCent->Fill(fCent);
 779
 780     if(Nmatches < 3)
 781       fHistNCellsEnergy[fCentBin][Nmatches]->Fill(energyclus, cNcells);
 782     else
 783       fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
 784
 785     if (EoP > 0) {
 786       fHistEoPCent->Fill(fCent, EoP);
 787       fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
 788     }
 789
 790     if (Nmatches == 1 && totalTrkP > 0) {
 791       Int_t iMin = emccluster->GetMatchedObjId();
 792       AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
 793       AliVTrack *track = emctrack->GetTrack();
 794       Int_t centbinchm = fCentBin;
 795       if (track->Charge()<0)
 796         centbinchm += 4;
 797
 798       fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
 799       fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
 800     }
 801   }
 802
 803   // apply the correction
 804   energyclus -= Esub;
 805
 806   return energyclus;
 807 }
 808