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