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