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