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