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