AliPhysics  b0c77bb (b0c77bb)
AliEmcalCorrectionClusterHadronicCorrection.cxx
Go to the documentation of this file.
1 // AliEmcalCorrectionClusterHadronicCorrection
2 //
3 
5 
6 #include <map>
7 #include <utility>
8 
9 #include <TH2.h>
10 #include <TF1.h>
11 #include <TList.h>
12 
13 #include "AliClusterContainer.h"
14 #include "AliParticleContainer.h"
15 #include "AliVTrack.h"
16 
20 
21 // Actually registers the class with the base class
23 
28  AliEmcalCorrectionComponent("AliEmcalCorrectionClusterHadronicCorrection"),
29  fPhiMatch(0.05),
30  fEtaMatch(0.025),
31  fDoTrackClus(kTRUE),
32  fDoMomDepMatching(kFALSE),
33  fHadCorr(0),
34  fEexclCell(0),
35  fPlotOversubtractionHistograms(kFALSE),
36  fDoNotOversubtract(kFALSE),
37  fUseM02SubtractionScheme(kFALSE),
38  fUseConstantSubtraction(kFALSE),
39  fConstantSubtractionValue(0.),
40  fClusterContainerIndexMap(),
41  fParticleContainerIndexMap(),
42  fTrackToContainerMap(),
43  fHistMatchEtaPhiAll(0),
44  fHistMatchEtaPhiAllCl(0),
45  fHistNclusvsCent(0),
46  fHistNclusMatchvsCent(0),
47  fHistEbefore(0),
48  fHistEafter(0),
49  fHistEoPCent(0),
50  fHistNMatchCent(0),
51  fHistNClusMatchCent(0)
52 {
53  for(Int_t i=0; i<10; i++) {
54  fHistEsubPch[i] = 0;
55  fHistEsubPchRat[i] = 0;
56 
57  if (i<5) {
58  fHistEsubPchRatAll[i] = 0;
59 
60  fHistMatchEvsP[i] = 0;
61  fHistMatchdRvsEP[i] = 0;
62  fHistNMatchEnergy[i] = 0;
63 
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 
85 {
86 }
87 
92 {
93  // Initialization
95 
96  GetProperty("phiMatch", fPhiMatch);
97  GetProperty("etaMatch", fEtaMatch);
98  GetProperty("hadCorr", fHadCorr);
99  GetProperty("Eexcl", fEexclCell);
100  GetProperty("doTrackClus", fDoTrackClus);
101  GetProperty("doMomDepMatching", fDoMomDepMatching);
102  GetProperty("plotOversubtractionHistograms", fPlotOversubtractionHistograms);
103  GetProperty("doNotOversubtract", fDoNotOversubtract);
104  GetProperty("useM02SubtractionScheme", fUseM02SubtractionScheme);
105  GetProperty("useConstantSubtraction", fUseConstantSubtraction);
106  GetProperty("constantSubtractionValue", fConstantSubtractionValue);
107 
108  return kTRUE;
109 }
110 
115 {
117 
118  // Create my user objects.
119 
120  if (!fCreateHisto) return;
121 
122  fHistMatchEtaPhiAll = new TH2F("fHistMatchEtaPhiAll", "fHistMatchEtaPhiAll;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
124 
125  fHistMatchEtaPhiAllCl = new TH2F("fHistMatchEtaPhiAllCl", "fHistMatchEtaPhiAllCl;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
127 
128  TString name;
129  TString temp;
130  const Int_t nCentChBins = fNcentBins * 2;
131  for(Int_t icent=0; icent<nCentChBins; ++icent) {
132  for(Int_t ipt=0; ipt<9; ++ipt) {
133  for(Int_t ieta=0; ieta<2; ++ieta) {
134  name = Form("fHistMatchEtaPhi_%i_%i_%i",icent,ipt,ieta);
135  fHistMatchEtaPhi[icent][ipt][ieta] = new TH2F(name, name, fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
136  fHistMatchEtaPhi[icent][ipt][ieta]->SetXTitle("#Delta#eta");
137  fHistMatchEtaPhi[icent][ipt][ieta]->SetYTitle("#Delta#phi");
138  fOutput->Add(fHistMatchEtaPhi[icent][ipt][ieta]);
139  }
140  }
141 
142  name = Form("fHistEsubPch_%i",icent);
143  temp = Form("%s (Nmatches==1)",name.Data());
144  fHistEsubPch[icent]=new TH1F(name, temp, fNbins, fMinBinPt, fMaxBinPt);
145  fHistEsubPch[icent]->SetXTitle("#sum p (GeV) weighted with E_{sub}");
146  fOutput->Add(fHistEsubPch[icent]);
147 
148  name = Form("fHistEsubPchRat_%i",icent);
149  temp = Form("%s (Nmatches==1)",name.Data());
150  fHistEsubPchRat[icent]=new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
151  fHistEsubPchRat[icent]->SetXTitle("#Sigma p (GeV)");
152  fHistEsubPchRat[icent]->SetYTitle("E_{sub} / #sum p");
153  fOutput->Add(fHistEsubPchRat[icent]);
154 
155  if (icent<fNcentBins) {
156  for(Int_t itrk=0; itrk<4; ++itrk) {
157  name = Form("fHistNCellsEnergy_%i_%i",icent,itrk);
158  temp = Form("%s (Nmatches==%d);N_{cells};E_{clus} (GeV)",name.Data(),itrk);
159  fHistNCellsEnergy[icent][itrk] = new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, 101, -0.5, 100.5);
160  fOutput->Add(fHistNCellsEnergy[icent][itrk]);
161  }
162 
163  name = Form("fHistEsubPchRatAll_%i",icent);
164  temp = Form("%s (all Nmatches)",name.Data());
165  fHistEsubPchRatAll[icent]=new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
166  fHistEsubPchRatAll[icent]->SetXTitle("#Sigma p (GeV)");
167  fHistEsubPchRatAll[icent]->SetYTitle("E_{sub} / #sum p");
168  fOutput->Add(fHistEsubPchRatAll[icent]);
169 
170  name = Form("fHistMatchEvsP_%i",icent);
171  temp = Form("%s (all Nmatches)",name.Data());
172  fHistMatchEvsP[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
173  fHistMatchEvsP[icent]->SetXTitle("E_{clus} (GeV)");
174  fHistMatchEvsP[icent]->SetYTitle("E_{clus} / #sum p");
175  fOutput->Add(fHistMatchEvsP[icent]);
176 
177  name = Form("fHistMatchdRvsEP_%i",icent);
178  temp = Form("%s (all Nmatches)",name.Data());
179  fHistMatchdRvsEP[icent] = new TH2F(name, temp, fNbins, 0., 0.2, fNbins*2, 0., 10.);
180  fHistMatchdRvsEP[icent]->SetXTitle("#Delta R between track and cluster");
181  fHistMatchdRvsEP[icent]->SetYTitle("E_{clus} / p");
182  fOutput->Add(fHistMatchdRvsEP[icent]);
183 
184  name = Form("fHistNMatchEnergy_%i",icent);
185  fHistNMatchEnergy[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, 101, -0.5, 100.5);
186  fHistNMatchEnergy[icent]->SetXTitle("E_{clus} (GeV)");
187  fHistNMatchEnergy[icent]->SetYTitle("N_{matches}");
188  fOutput->Add(fHistNMatchEnergy[icent]);
189  }
190  }
191 
192  fHistNclusvsCent = new TH1F("Nclusvscent", "NclusVsCent; Cent (%)", 100, 0, 100);
193  fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent (all Nmatches); Cent (%)", 100, 0, 100);
194  fHistEbefore = new TH1F("Ebefore", "Ebefore; Cent (%); E_{clus} (GeV)", 100, 0, 100);
195  fHistEafter = new TH1F("Eafter", "Eafter; Cent (%); E_{clus} (GeV)", 100, 0, 100);
196  fHistEoPCent = new TH2F("EoPCent", "EoPCent; Cent (%); E_{clus} / #sum p", 100, 0, 100, fNbins*2, 0, 10);
197  fHistNMatchCent = new TH2F("NMatchesCent", "NMatchesCent; Cent (%); Nmatches", 100, 0, 100, 11, -0.5, 10.5);
198  fHistNClusMatchCent = new TH2F("NClusMatchesCent", "NClusMatchesCent; Cent (%); Nmatches", 100, 0, 100, 11, -0.5, 10.5);
199 
202  fOutput->Add(fHistEbefore);
203  fOutput->Add(fHistEafter);
204  fOutput->Add(fHistEoPCent);
205  fOutput->Add(fHistNMatchCent);
207 
209  for(Int_t icent=0; icent<fNcentBins; ++icent) {
210  name = Form("fHistEmbTrackMatchesOversub_%d",icent);
211  fHistEmbTrackMatchesOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
212  fHistEmbTrackMatchesOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
213  fHistEmbTrackMatchesOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
215 
216  name = Form("fHistNonEmbTrackMatchesOversub_%d",icent);
217  fHistNonEmbTrackMatchesOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
218  fHistNonEmbTrackMatchesOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
219  fHistNonEmbTrackMatchesOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
221 
222  name = Form("fHistOversubMCClusters_%d",icent);
223  fHistOversubMCClusters[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
224  fHistOversubMCClusters[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
225  fHistOversubMCClusters[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
226  fOutput->Add(fHistOversubMCClusters[icent]);
227 
228  name = Form("fHistOversubNonMCClusters_%d",icent);
229  fHistOversubNonMCClusters[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
230  fHistOversubNonMCClusters[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
231  fHistOversubNonMCClusters[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
232  fOutput->Add(fHistOversubNonMCClusters[icent]);
233 
234  name = Form("fHistOversub_%d",icent);
235  fHistOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
236  fHistOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
237  fHistOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
238  fOutput->Add(fHistOversub[icent]);
239  }
240  }
241  fOutput->SetOwner(kTRUE);
242 }
243 
248 {
251 }
252 
257 {
259 
260  // Build map from all available tracks to their corresponding particle containers.
261  // This is only required for AODs which store pointers to AliVTracks instead of
262  // indices.
263  if (!fEsdMode) {
265  }
266 
267  // Run the hadronic correction
268  // loop over all clusters
269  AliVCluster *cluster = 0;
270  AliClusterContainer * clusCont = 0;
271  TIter nextClusCont(&fClusterCollArray);
272  while ((clusCont = static_cast<AliClusterContainer*>(nextClusCont()))) {
273  auto clusItCont = clusCont->accepted_momentum();
274  for (AliClusterIterableMomentumContainer::iterator clusIterator = clusItCont.begin(); clusIterator != clusItCont.end(); ++clusIterator) {
275  cluster = static_cast<AliVCluster *>(clusIterator->second);
276 
277  Double_t energyclus = 0;
278  if (fCreateHisto) {
279  fHistEbefore->Fill(fCent, cluster->GetNonLinCorrEnergy());
280  fHistNclusvsCent->Fill(fCent);
281  }
282 
283  // apply correction / subtraction
284  // to subtract only the closest track set fHadCor to a %
285  // to subtract all tracks within the cut set fHadCor to %+1
286  if (fHadCorr > 1) {
287  energyclus = ApplyHadCorrAllTracks(fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, clusIterator.current_index()), fHadCorr - 1);
288  }
289  else if (fHadCorr > 0) {
290  energyclus = ApplyHadCorrOneTrack(fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, clusIterator.current_index()), fHadCorr);
291  }
292  else {
293  energyclus = cluster->GetNonLinCorrEnergy();
294  }
295 
296  if (energyclus < 0) energyclus = 0;
297 
298  cluster->SetHadCorrEnergy(energyclus);
299 
300  if (fCreateHisto) fHistEafter->Fill(fCent, energyclus);
301 
302  }
303  }
304 
305  return kTRUE;
306 }
307 
318 {
319  fTrackToContainerMap.clear();
320 
321  // Loop by index to avoid the possibility of iterators being invalidated.
322  AliParticleContainer * partCont = 0;
323  for (int i = 0; i < fParticleCollArray.GetEntriesFast(); i++)
324  {
325  partCont = GetParticleContainer(i);
326  if (!partCont) { AliErrorStream() << "Failed to retrieve particle container at index " << i << "\n"; }
327  for (int j = 0; j < partCont->GetNParticles(); j++)
328  {
329  AliVTrack * track = static_cast<AliVTrack *>(partCont->GetParticle(j));
330  fTrackToContainerMap.insert(std::make_pair(track, partCont));
331  }
332  }
333 }
334 
339 {
340  UInt_t pbin=0;
341  if (p<0.5)
342  pbin=0;
343  else if (p>=0.5 && p<1.0)
344  pbin=1;
345  else if (p>=1.0 && p<1.5)
346  pbin=2;
347  else if (p>=1.5 && p<2.)
348  pbin=3;
349  else if (p>=2. && p<3.)
350  pbin=4;
351  else if (p>=3. && p<4.)
352  pbin=5;
353  else if (p>=4. && p<5.)
354  pbin=6;
355  else if (p>=5. && p<8.)
356  pbin=7;
357  else
358  pbin=8;
359 
360  return pbin;
361 }
362 
367 {
368  Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
369  return 2.0*EtaSigma[pbin];
370 }
371 
376 {
377  if (centbin==0) {
378  Double_t PhiMean[9]={0.036,
379  0.021,
380  0.0121,
381  0.0084,
382  0.0060,
383  0.0041,
384  0.0031,
385  0.0022,
386  0.001};
387  return PhiMean[pbin];
388  } else if (centbin==1) {
389  Double_t PhiMean[9]={0.036,
390  0.021,
391  0.0121,
392  0.0084,
393  0.0060,
394  0.0041,
395  0.0031,
396  0.0022,
397  0.001};
398  return PhiMean[pbin];
399  } else if (centbin==2) {
400  Double_t PhiMean[9]={0.036,
401  0.021,
402  0.0121,
403  0.0084,
404  0.0060,
405  0.0041,
406  0.0031,
407  0.0022,
408  0.001};
409  return PhiMean[pbin];
410  } else if (centbin==3) {
411  Double_t PhiMean[9]={0.036,
412  0.021,
413  0.0121,
414  0.0084,
415  0.0060,
416  0.0041,
417  0.0031,
418  0.0022,
419  0.001};
420  return PhiMean[pbin];
421  } else if (centbin==4) {
422  Double_t PhiMean[9]={0.036,
423  0.021,
424  0.0127,
425  0.0089,
426  0.0068,
427  0.0049,
428  0.0038,
429  0.0028,
430  0.0018};
431  return PhiMean[pbin]*(-1.);
432  } else if (centbin==5) {
433  Double_t PhiMean[9]={0.036,
434  0.021,
435  0.0127,
436  0.0089,
437  0.0068,
438  0.0048,
439  0.0038,
440  0.0028,
441  0.0018};
442  return PhiMean[pbin]*(-1.);
443  } else if (centbin==6) {
444  Double_t PhiMean[9]={0.036,
445  0.021,
446  0.0127,
447  0.0089,
448  0.0068,
449  0.0045,
450  0.0035,
451  0.0028,
452  0.0018};
453  return PhiMean[pbin]*(-1.);
454  } else if (centbin==7) {
455  Double_t PhiMean[9]={0.036,
456  0.021,
457  0.0127,
458  0.0089,
459  0.0068,
460  0.0043,
461  0.0035,
462  0.0028,
463  0.0018};
464  return PhiMean[pbin]*(-1.);
465  }
466 
467  return 0;
468 }
469 
474 {
475  if (centbin==0) {
476  Double_t PhiSigma[9]={0.0221,
477  0.0128,
478  0.0074,
479  0.0064,
480  0.0059,
481  0.0055,
482  0.0052,
483  0.0049,
484  0.0045};
485  return 2.*PhiSigma[pbin];
486  } else if (centbin==1) {
487  Double_t PhiSigma[9]={0.0217,
488  0.0120,
489  0.0076,
490  0.0066,
491  0.0062,
492  0.0058,
493  0.0054,
494  0.0054,
495  0.0045};
496  return 2.*PhiSigma[pbin];
497  } else if (centbin==2) {
498  Double_t PhiSigma[9]={0.0211,
499  0.0124,
500  0.0080,
501  0.0070,
502  0.0067,
503  0.0061,
504  0.0059,
505  0.0054,
506  0.0047};
507  return 2.*PhiSigma[pbin];
508  } else if (centbin==3) {
509  Double_t PhiSigma[9]={0.0215,
510  0.0124,
511  0.0082,
512  0.0073,
513  0.0069,
514  0.0064,
515  0.0060,
516  0.0055,
517  0.0047};
518  return 2.*PhiSigma[pbin];
519  } else if (centbin==4) {
520  Double_t PhiSigma[9]={0.0199,
521  0.0108,
522  0.0072,
523  0.0071,
524  0.0060,
525  0.0055,
526  0.0052,
527  0.0049,
528  0.0045};
529  return 2.*PhiSigma[pbin];
530  } else if (centbin==5) {
531  Double_t PhiSigma[9]={0.0200,
532  0.0110,
533  0.0074,
534  0.0071,
535  0.0064,
536  0.0059,
537  0.0055,
538  0.0052,
539  0.0045};
540  return 2.*PhiSigma[pbin];
541  } else if (centbin==6) {
542  Double_t PhiSigma[9]={0.0202,
543  0.0113,
544  0.0077,
545  0.0071,
546  0.0069,
547  0.0064,
548  0.0060,
549  0.0055,
550  0.0050};
551  return 2.*PhiSigma[pbin];
552  } else if (centbin==7) {
553  Double_t PhiSigma[9]={0.0205,
554  0.0113,
555  0.0080,
556  0.0074,
557  0.0078,
558  0.0067,
559  0.0062,
560  0.0055,
561  0.0050};
562  return 2.*PhiSigma[pbin];
563  }
564 
565  return 0;
566 }
567 
572  Double_t &totalTrkP, Int_t &Nmatches, Double_t &trkPMCfrac, Int_t &NMCmatches)
573 {
574  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
575 
576  if (!cluster) return;
577  //These are the functional forms for the momentum dependent eta-phi track matching cut
578  //Values taken from the PCM analyses see:
579  //https://alice-notes.web.cern.ch/node/411 Fig. 46 for mom<3GeV these cuts are wider than the standard cuts
580  //these values were extracted in the 2012 pp8 TeV MonteCarlo and apparently showed robustness throughout different periods
581  TF1 funcPtDepEta("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
582  funcPtDepEta.SetParameters(0.04, 0.010, 2.5);
583  TF1 funcPtDepPhi("func", "[1] + 1 / pow(x + pow(1 / ([0] - [1]), 1 / [2]), [2])");
584  funcPtDepPhi.SetParameters(0.09, 0.015, 2.);
585 
586  // loop over matched tracks
587  Int_t Ntrks = cluster->GetNTracksMatched();
588  for (Int_t i = 0; i < Ntrks; ++i) {
589  AliVTrack* track = 0;
590 
591  if (fEsdMode) {
592  Int_t itrack = cluster->GetTrackMatchedIndex(i);
593  if (itrack >= 0) {
595  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
596  }
597  }
598  else {
599  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(i));
600  UInt_t rejectionReason = 0;
601  AliParticleContainer * partCont = fTrackToContainerMap.at(track);
602  if (!partCont) { AliErrorStream() << "Requested particle container not available!\n"; }
603  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
604  }
605 
606  if (!track) continue;
607 
608  Double_t etadiff = 999;
609  Double_t phidiff = 999;
610  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
611  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(etadiff, phidiff);
612 
613  // check if track also points to cluster
614  if (fDoTrackClus && (track->GetEMCALcluster() != icluster)) continue;
615 
616  Double_t mom = track->P();
617  UInt_t mombin = GetMomBin(mom);
618  Int_t centbinch = fCentBin;
619 
620  if (track->Charge() < 0) centbinch += fNcentBins;
621 
622  if (fCreateHisto) {
623  Int_t etabin = 0;
624  if (track->Eta() > 0) etabin=1;
625  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
626  fHistMatchEtaPhiAll->Fill(etadiff, phidiff);
627  }
628 
629  Double_t etaCut = 0.0;
630  Double_t phiCutlo = 0.0;
631  Double_t phiCuthi = 0.0;
632 
633  if (fPhiMatch > 0) {
634  phiCutlo = -fPhiMatch;
635  phiCuthi = +fPhiMatch;
636  }
637  else {
638  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
639  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
640  }
641 
642  if (fEtaMatch > 0) {
643  etaCut = fEtaMatch;
644  }
645  else {
646  etaCut = GetEtaSigma(mombin);
647  }
648  //Do momentum dependent track matching
649  if (fDoMomDepMatching) {
650  phiCutlo = -funcPtDepPhi.Eval(mom);
651  phiCuthi = +funcPtDepPhi.Eval(mom);
652  etaCut = funcPtDepEta.Eval(mom);
653  }
654 
655  if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
656  if (track->GetLabel() > fMinMCLabel) {
657  ++NMCmatches;
658  trkPMCfrac += mom;
659  }
660  ++Nmatches;
661  totalTrkP += mom;
662 
663  if (fCreateHisto) {
664  if (fHadCorr > 1) {
665  Double_t dphi = 0;
666  Double_t deta = 0;
667  GetEtaPhiDiff(track, cluster, dphi, deta);
668  Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
669  Double_t energyclus = cluster->GetNonLinCorrEnergy();
670  fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
671  }
672  }
673  }
674  }
675 
676  if (totalTrkP > 0) trkPMCfrac /= totalTrkP;
677 }
678 
683 {
684  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
685 
686  Double_t energyclus = cluster->GetNonLinCorrEnergy();
687 
688  AliVTrack* track = 0;
689 
690  if (cluster->GetNTracksMatched() > 0) {
691  if (fEsdMode) {
692  Int_t itrack = cluster->GetTrackMatchedIndex(0);
693  if (itrack >= 0) {
695  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
696  }
697  }
698  else {
699  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
700  UInt_t rejectionReason = 0;
701  AliParticleContainer * partCont = fTrackToContainerMap.at(track);
702  if (!partCont) { AliErrorStream() << "Requested particle container not available!\n"; }
703  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
704  }
705  }
706 
707  if (!track || track->P() < 1e-6) return energyclus;
708 
709  Double_t mom = track->P();
710 
711  Double_t dEtaMin = 1e9;
712  Double_t dPhiMin = 1e9;
713  GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
714 
715  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(dEtaMin, dPhiMin);
716 
717  // check if track also points to cluster
718  Int_t cid = track->GetEMCALcluster();
719  if (fDoTrackClus && (cid != icluster)) return energyclus;
720 
721  UInt_t mombin = GetMomBin(mom);
722  Int_t centbinch = fCentBin;
723  if (track->Charge() < 0) centbinch += fNcentBins;
724 
725  // plot some histograms if switched on
726  if (fCreateHisto) {
727  Int_t etabin = 0;
728  if(track->Eta() > 0) etabin = 1;
729 
730  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
731  fHistMatchEtaPhiAll->Fill(dEtaMin, dPhiMin);
732 
733  if (mom > 0) {
734  Double_t etadiff = 0;
735  Double_t phidiff = 0;
736  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
737  Double_t dRmin = TMath::Sqrt(etadiff*etadiff + phidiff*phidiff);
738  fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
739  fHistEoPCent->Fill(fCent, energyclus / mom);
740  fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
741  }
742  }
743 
744  // define eta/phi cuts
745  Double_t etaCut = 0.0;
746  Double_t phiCutlo = 0.0;
747  Double_t phiCuthi = 0.0;
748  if (fPhiMatch > 0) {
749  phiCutlo = -fPhiMatch;
750  phiCuthi = +fPhiMatch;
751  }
752  else {
753  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
754  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
755  }
756  if (fEtaMatch > 0) {
757  etaCut = fEtaMatch;
758  }
759  else {
760  etaCut = GetEtaSigma(mombin);
761  }
762 
763  // apply the correction if the track is in the eta/phi window
764  if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
765  energyclus -= hadCorr * mom;
766  }
767 
768  return energyclus;
769 }
770 
775 {
776  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
777 
778  Double_t energyclus = cluster->GetNonLinCorrEnergy();
779  Double_t cNcells = cluster->GetNCells();
780 
781  Double_t totalTrkP = 0.0; // count total track momentum
782  Int_t Nmatches = 0; // count total number of matches
783 
784  Double_t trkPMCfrac = 0.0; // count total track momentum
785  Int_t NMCmatches = 0; // count total number of matches
786 
787  // do the loop over the matched tracks and get the number of matches and the total momentum
788  DoMatchedTracksLoop(icluster, totalTrkP, Nmatches, trkPMCfrac, NMCmatches);
789 
790  Double_t Esub = hadCorr * totalTrkP;
791 
792  if (Esub > energyclus) Esub = energyclus;
793 
794  // applying Peter's proposed algorithm
795  // never subtract the full energy of the cluster
796  Double_t clusEexcl = fEexclCell * cNcells;
797  if (energyclus < clusEexcl) clusEexcl = energyclus;
798  if ((energyclus - Esub) < clusEexcl) Esub = (energyclus - clusEexcl);
799 
800  // If enabled, use cluster M02 value to determine how to correct cluster energy.
802  Esub = ComputeM02Subtraction(cluster, energyclus, Nmatches, totalTrkP, hadCorr);
803  }
804 
805  // embedding
806  Double_t EsubMC = 0;
807  Double_t EsubBkg = 0;
808  Double_t EclusMC = 0;
809  Double_t EclusBkg = 0;
810  Double_t EclusCorr = 0;
811  Double_t EclusMCcorr = 0;
812  Double_t EclusBkgcorr = 0;
813  Double_t overSub = 0;
815  EsubMC = hadCorr * totalTrkP * trkPMCfrac;
816  EsubBkg = hadCorr * totalTrkP - EsubMC;
817  EclusMC = energyclus * cluster->GetMCEnergyFraction();
818  EclusBkg = energyclus - EclusMC;
819 
820  if (energyclus > Esub)
821  EclusCorr = energyclus - Esub;
822 
823  if (EclusMC > EsubMC)
824  EclusMCcorr = EclusMC - EsubMC;
825 
826  if (EclusBkg > EsubBkg)
827  EclusBkgcorr = EclusBkg - EsubBkg;
828 
829  overSub = EclusMCcorr + EclusBkgcorr - EclusCorr;
830  }
831 
832  // plot some histograms if switched on
833  if (fCreateHisto) {
834  fHistNMatchCent->Fill(fCent, Nmatches);
835  fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
836 
837  if (Nmatches > 0) fHistNclusMatchvsCent->Fill(fCent);
838 
839  if (Nmatches < 3) {
840  fHistNCellsEnergy[fCentBin][Nmatches]->Fill(energyclus, cNcells);
841  }
842  else {
843  fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
844  }
845 
846  if (totalTrkP > 0) {
847  Double_t EoP = energyclus / totalTrkP;
848  fHistEoPCent->Fill(fCent, EoP);
849  fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
850 
851  fHistEsubPchRatAll[fCentBin]->Fill(totalTrkP, Esub / totalTrkP);
852 
853  if (Nmatches == 1) {
854  AliVTrack* track = 0;
855  if (fEsdMode) {
856  Int_t itrack = cluster->GetTrackMatchedIndex(0);
857  if (itrack >= 0) {
859  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
860  }
861  }
862  else {
863  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
864  UInt_t rejectionReason = 0;
865  AliParticleContainer * partCont = fTrackToContainerMap.at(track);
866  if (!partCont) { AliErrorStream() << "Requested particle container not available!\n"; }
867  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
868  }
869  if (track) {
870  Int_t centbinchm = fCentBin;
871  if (track->Charge() < 0) centbinchm += fNcentBins;
872  fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
873  fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
874  }
875  }
876 
878  fHistOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
879 
880  if (cluster->GetMCEnergyFraction() > 0.95)
881  fHistOversubMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
882  else if (cluster->GetMCEnergyFraction() < 0.05)
883  fHistOversubNonMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
884 
885  if (trkPMCfrac < 0.05)
886  fHistNonEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
887  else if (trkPMCfrac > 0.95)
888  fHistEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
889  }
890  }
891  }
892 
894  Esub -= overSub;
895  if (EclusBkgcorr + EclusMCcorr > 0) {
896  Double_t newfrac = EclusMCcorr / (EclusBkgcorr + EclusMCcorr);
897  cluster->SetMCEnergyFraction(newfrac);
898  }
899  }
900 
901  // apply the correction
902  energyclus -= Esub;
903 
904  return energyclus;
905 }
906 
910 Double_t AliEmcalCorrectionClusterHadronicCorrection::ComputeM02Subtraction(const AliVCluster* cluster, Double_t energyclus, Int_t Nmatches, Double_t totalTrkP, Double_t hadCorr)
911 {
912  Double_t Esub = 0.;
913  Double_t clusM02 = cluster->GetM02();
914 
915  // Note: The comments below give only a rough indication of what particle types appear in each selection.
916 
917  // For M02 in the single photon region, the signal is primarily: Single photons, single electrons, single MIPs
918  if (clusM02 > 0.1 && clusM02 < 0.4) {
919  if (Nmatches == 0) { // Single photon (modulo tracking efficiency)
920  Esub = 0;
921  }
922  else { // Single electron, single MIP
923  Esub = energyclus;
924  }
925  }
926 
927  // For large M02, the signal is primarily: Single hadronic shower, photon-photon overlap, photon-MIP overlap, MIP-MIP overlap,
928  // MIP-hadronic shower overlap, hadronic shower - hadronic shower overlap)
929  if (clusM02 > 0.4) {
930  if (Nmatches == 0) { // Single neutral hadronic shower, photon-photon overlap, photon-neutral had shower overlap, neutral had shower overlap (modulo tracking efficiency)
931  Esub = 0;
932  }
933  else if (Nmatches == 1) { // Single charged hadronic shower, photon-MIP overlap, MIP-MIP overlap, MIP-had shower overlap, had shower-shower overlap
936  }
937  else {
938  Esub = hadCorr * totalTrkP;
939  }
940  }
941  else if (Nmatches > 1) { // MIP-MIP overlap, had shower-shower overlap
942  Esub = energyclus;
943  }
944  }
945 
946  return Esub;
947 }
Int_t fNcentBins
How many centrality bins (this member copied from AliAnalysisTaskEmcal)
Bool_t fDoMomDepMatching
set the values fPhiMatch and fEtaMatch depending on the track momentum
AliParticleContainer * GetParticleContainer(Int_t i=0) const
double Double_t
Definition: External.C:58
Definition: External.C:236
TH2 * fHistMatchEtaPhiAll
!deta vs. dphi of matched cluster-track pairs
TH2 * fHistMatchEtaPhi[10][9][2]
!deta vs. dphi of matched cluster-track pairs
const AliClusterIterableMomentumContainer accepted_momentum() const
TH2 * fHistNonEmbTrackMatchesOversub[5]
!Over-subtracted energy / cluster energy with non-embedded track matches (embedded matches < 5%) ...
bidirectional stl iterator over the EMCAL iterable container
std::pair< int, U * > LocalIndexFromGlobalIndex(const int globalIndex) const
TH2 * fHistOversubMCClusters[5]
!Over-subtracted energy / cluster energy (cluster MC energy fraction > 95%)
Int_t GetNParticles() const
Bool_t fUseConstantSubtraction
Flag to perform constant rather than fractional subtract (only applicable if using M02 scheme) ...
static RegisterCorrectionComponent< AliEmcalCorrectionClusterHadronicCorrection > reg
TH2 * fHistOversubNonMCClusters[5]
!Over-subtracted energy / cluster energy (cluster MC energy fraction < 5%)
Bool_t fDoNotOversubtract
do not oversubtract energy from cluster (embedding only)
AliEmcalContainerIndexMap< AliClusterContainer, AliVCluster > fClusterContainerIndexMap
! Mapping between index and cluster containers
TH2 * fHistMatchEvsP[5]
!cluster energy vs. track momentum of matched pairs
Int_t fCentBin
! Event centrality bin
int GlobalIndexFromLocalIndex(const U *inputObject, const int localIndex) const
TH1 * fHistEsubPch[10]
!Esub vs. total momentum of matched tracks (only 1 match)
Container for particles within the EMCAL framework.
TObjArray fParticleCollArray
Particle/track collection array.
TObjArray fClusterCollArray
Cluster collection array.
TH1 * fHistEafter
!average energy of clusters after correction vs. centrality
TH2 * fHistEsubPchRatAll[5]
!Esub/momentum of matched tracks vs. total momentum of matched tracks (all number of matches) ...
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
Base class for correction components in the EMCal correction framework.
TH2 * fHistEsubPchRat[10]
!Esub/momentum of matched tracks vs. total momentum of matched tracks (only 1 match) ...
virtual AliVParticle * GetParticle(Int_t i=-1) const
Double_t fEexclCell
energy/cell that we cannot subtract from the clusters
static const AliEmcalContainerIndexMap< TClonesArray, AliVCluster > & GetEmcalContainerIndexMap()
Get the EMCal container utils associated with particle containers.
std::map< AliVTrack *, AliParticleContainer * > fTrackToContainerMap
! Mapping between AliVTracks and their respective particle containers. Needed for AODs only...
Int_t fMinMCLabel
Minimum MC label value for the tracks/clusters being considered MC particles.
V * GetObjectFromGlobalIndex(const int globalIndex) const
virtual Bool_t AcceptParticle(const AliVParticle *vp, UInt_t &rejectionReason) const
Double_t fMaxBinPt
Max pt in histograms.
Bool_t fPlotOversubtractionHistograms
compute and plot oversubtracted energy from embedded/signal matches (embedding only) ...
TList * fOutput
! List of output histograms
Double_t ComputeM02Subtraction(const AliVCluster *cluster, Double_t energyclus, Int_t Nmatches, Double_t totalTrkP, Double_t hadCorr)
Bool_t fCreateHisto
Flag to make some basic histograms.
Hadronic correction component in the EMCal correction framework.
TH1 * fHistEbefore
!average energy of clusters before correction vs. centrality
Double_t fCent
! Event centrality
void DoMatchedTracksLoop(Int_t icluster, Double_t &totalTrkP, Int_t &Nmatches, Double_t &trkPMCfrac, Int_t &NMCmatches)
Double_t fMinBinPt
Min pt in histograms.
void CopyMappingFrom(const AliEmcalContainerIndexMap< U2, V > &map, U *cont)
TH1 * fHistNclusMatchvsCent
!n clusters matched to some track vs. centrality
void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
bool Bool_t
Definition: External.C:53
AliEmcalContainerIndexMap< AliParticleContainer, AliVParticle > fParticleContainerIndexMap
! Mapping between index and particle containers
Bool_t fUseM02SubtractionScheme
Flag to enable hadronic correction scheme using cluster M02 value.
Container structure for EMCAL clusters.
TH2 * fHistOversub[5]
!Over-subtracted energy / cluster energy
TH2 * fHistEmbTrackMatchesOversub[5]
!Over-subtracted energy / cluster energy with embedded track matches (non-embedded matches < 5%) ...
static const AliEmcalContainerIndexMap< TClonesArray, AliVParticle > & GetEmcalContainerIndexMap()
Get the EMCal container utils associated with particle containers.
TH2 * fHistMatchEtaPhiAllCl
!deta vs. dphi of all cluster-track pairs (cl loop)
TH2 * fHistNClusMatchCent
!n clusters macthed to some track (tracks allowed to match more than one cluster) ...
bool GetProperty(std::string propertyName, T &property, bool requiredProperty=true, std::string correctionName="")
Retrieve property.
Double_t fConstantSubtractionValue
Value to be used for constant subtraction (only applicable if using constant subtraction in M02 schem...