AliPhysics  41af4b0 (41af4b0)
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 <TList.h>
11 
12 #include "AliClusterContainer.h"
13 #include "AliParticleContainer.h"
14 #include "AliVTrack.h"
15 
19 
20 // Actually registers the class with the base class
22 
27  AliEmcalCorrectionComponent("AliEmcalCorrectionClusterHadronicCorrection"),
28  fPhiMatch(0.05),
29  fEtaMatch(0.025),
30  fDoTrackClus(kTRUE),
31  fHadCorr(0),
32  fEexclCell(0),
33  fPlotOversubtractionHistograms(kFALSE),
34  fDoNotOversubtract(kFALSE),
35  fUseM02SubtractionScheme(kFALSE),
36  fUseConstantSubtraction(kFALSE),
37  fConstantSubtractionValue(0.),
38  fClusterContainerIndexMap(),
39  fParticleContainerIndexMap(),
40  fTrackToContainerMap(),
41  fHistMatchEtaPhiAll(0),
42  fHistMatchEtaPhiAllCl(0),
43  fHistNclusvsCent(0),
44  fHistNclusMatchvsCent(0),
45  fHistEbefore(0),
46  fHistEafter(0),
47  fHistEoPCent(0),
48  fHistNMatchCent(0),
49  fHistNClusMatchCent(0)
50 {
51  for(Int_t i=0; i<10; i++) {
52  fHistEsubPch[i] = 0;
53  fHistEsubPchRat[i] = 0;
54 
55  if (i<5) {
56  fHistEsubPchRatAll[i] = 0;
57 
58  fHistMatchEvsP[i] = 0;
59  fHistMatchdRvsEP[i] = 0;
60  fHistNMatchEnergy[i] = 0;
61 
66  fHistOversub[i] = 0;
67 
68  for(Int_t j=0; j<4; j++)
69  fHistNCellsEnergy[i][j] = 0;
70  }
71 
72  for(Int_t j=0; j<9; j++) {
73  for(Int_t k=0; k<2; k++)
74  fHistMatchEtaPhi[i][j][k] = 0;
75  }
76  }
77 }
78 
83 {
84 }
85 
90 {
91  // Initialization
93 
94  GetProperty("phiMatch", fPhiMatch);
95  GetProperty("etaMatch", fEtaMatch);
96  GetProperty("hadCorr", fHadCorr);
97  GetProperty("Eexcl", fEexclCell);
98  GetProperty("doTrackClus", fDoTrackClus);
99  GetProperty("plotOversubtractionHistograms", fPlotOversubtractionHistograms);
100  GetProperty("doNotOversubtract", fDoNotOversubtract);
101  GetProperty("useM02SubtractionScheme", fUseM02SubtractionScheme);
102  GetProperty("useConstantSubtraction", fUseConstantSubtraction);
103  GetProperty("constantSubtractionValue", fConstantSubtractionValue);
104 
105  return kTRUE;
106 }
107 
112 {
114 
115  // Create my user objects.
116 
117  if (!fCreateHisto) return;
118 
119  fHistMatchEtaPhiAll = new TH2F("fHistMatchEtaPhiAll", "fHistMatchEtaPhiAll;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
121 
122  fHistMatchEtaPhiAllCl = new TH2F("fHistMatchEtaPhiAllCl", "fHistMatchEtaPhiAllCl;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
124 
125  TString name;
126  TString temp;
127  const Int_t nCentChBins = fNcentBins * 2;
128  for(Int_t icent=0; icent<nCentChBins; ++icent) {
129  for(Int_t ipt=0; ipt<9; ++ipt) {
130  for(Int_t ieta=0; ieta<2; ++ieta) {
131  name = Form("fHistMatchEtaPhi_%i_%i_%i",icent,ipt,ieta);
132  fHistMatchEtaPhi[icent][ipt][ieta] = new TH2F(name, name, fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
133  fHistMatchEtaPhi[icent][ipt][ieta]->SetXTitle("#Delta#eta");
134  fHistMatchEtaPhi[icent][ipt][ieta]->SetYTitle("#Delta#phi");
135  fOutput->Add(fHistMatchEtaPhi[icent][ipt][ieta]);
136  }
137  }
138 
139  name = Form("fHistEsubPch_%i",icent);
140  temp = Form("%s (Nmatches==1)",name.Data());
141  fHistEsubPch[icent]=new TH1F(name, temp, fNbins, fMinBinPt, fMaxBinPt);
142  fHistEsubPch[icent]->SetXTitle("#sum p (GeV) weighted with E_{sub}");
143  fOutput->Add(fHistEsubPch[icent]);
144 
145  name = Form("fHistEsubPchRat_%i",icent);
146  temp = Form("%s (Nmatches==1)",name.Data());
147  fHistEsubPchRat[icent]=new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
148  fHistEsubPchRat[icent]->SetXTitle("#Sigma p (GeV)");
149  fHistEsubPchRat[icent]->SetYTitle("E_{sub} / #sum p");
150  fOutput->Add(fHistEsubPchRat[icent]);
151 
152  if (icent<fNcentBins) {
153  for(Int_t itrk=0; itrk<4; ++itrk) {
154  name = Form("fHistNCellsEnergy_%i_%i",icent,itrk);
155  temp = Form("%s (Nmatches==%d);N_{cells};E_{clus} (GeV)",name.Data(),itrk);
156  fHistNCellsEnergy[icent][itrk] = new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, 101, -0.5, 100.5);
157  fOutput->Add(fHistNCellsEnergy[icent][itrk]);
158  }
159 
160  name = Form("fHistEsubPchRatAll_%i",icent);
161  temp = Form("%s (all Nmatches)",name.Data());
162  fHistEsubPchRatAll[icent]=new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
163  fHistEsubPchRatAll[icent]->SetXTitle("#Sigma p (GeV)");
164  fHistEsubPchRatAll[icent]->SetYTitle("E_{sub} / #sum p");
165  fOutput->Add(fHistEsubPchRatAll[icent]);
166 
167  name = Form("fHistMatchEvsP_%i",icent);
168  temp = Form("%s (all Nmatches)",name.Data());
169  fHistMatchEvsP[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
170  fHistMatchEvsP[icent]->SetXTitle("E_{clus} (GeV)");
171  fHistMatchEvsP[icent]->SetYTitle("E_{clus} / #sum p");
172  fOutput->Add(fHistMatchEvsP[icent]);
173 
174  name = Form("fHistMatchdRvsEP_%i",icent);
175  temp = Form("%s (all Nmatches)",name.Data());
176  fHistMatchdRvsEP[icent] = new TH2F(name, temp, fNbins, 0., 0.2, fNbins*2, 0., 10.);
177  fHistMatchdRvsEP[icent]->SetXTitle("#Delta R between track and cluster");
178  fHistMatchdRvsEP[icent]->SetYTitle("E_{clus} / p");
179  fOutput->Add(fHistMatchdRvsEP[icent]);
180 
181  name = Form("fHistNMatchEnergy_%i",icent);
182  fHistNMatchEnergy[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, 101, -0.5, 100.5);
183  fHistNMatchEnergy[icent]->SetXTitle("E_{clus} (GeV)");
184  fHistNMatchEnergy[icent]->SetYTitle("N_{matches}");
185  fOutput->Add(fHistNMatchEnergy[icent]);
186  }
187  }
188 
189  fHistNclusvsCent = new TH1F("Nclusvscent", "NclusVsCent; Cent (%)", 100, 0, 100);
190  fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent (all Nmatches); Cent (%)", 100, 0, 100);
191  fHistEbefore = new TH1F("Ebefore", "Ebefore; Cent (%); E_{clus} (GeV)", 100, 0, 100);
192  fHistEafter = new TH1F("Eafter", "Eafter; Cent (%); E_{clus} (GeV)", 100, 0, 100);
193  fHistEoPCent = new TH2F("EoPCent", "EoPCent; Cent (%); E_{clus} / #sum p", 100, 0, 100, fNbins*2, 0, 10);
194  fHistNMatchCent = new TH2F("NMatchesCent", "NMatchesCent; Cent (%); Nmatches", 100, 0, 100, 11, -0.5, 10.5);
195  fHistNClusMatchCent = new TH2F("NClusMatchesCent", "NClusMatchesCent; Cent (%); Nmatches", 100, 0, 100, 11, -0.5, 10.5);
196 
199  fOutput->Add(fHistEbefore);
200  fOutput->Add(fHistEafter);
201  fOutput->Add(fHistEoPCent);
202  fOutput->Add(fHistNMatchCent);
204 
206  for(Int_t icent=0; icent<fNcentBins; ++icent) {
207  name = Form("fHistEmbTrackMatchesOversub_%d",icent);
208  fHistEmbTrackMatchesOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
209  fHistEmbTrackMatchesOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
210  fHistEmbTrackMatchesOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
212 
213  name = Form("fHistNonEmbTrackMatchesOversub_%d",icent);
214  fHistNonEmbTrackMatchesOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
215  fHistNonEmbTrackMatchesOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
216  fHistNonEmbTrackMatchesOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
218 
219  name = Form("fHistOversubMCClusters_%d",icent);
220  fHistOversubMCClusters[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
221  fHistOversubMCClusters[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
222  fHistOversubMCClusters[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
223  fOutput->Add(fHistOversubMCClusters[icent]);
224 
225  name = Form("fHistOversubNonMCClusters_%d",icent);
226  fHistOversubNonMCClusters[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
227  fHistOversubNonMCClusters[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
228  fHistOversubNonMCClusters[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
229  fOutput->Add(fHistOversubNonMCClusters[icent]);
230 
231  name = Form("fHistOversub_%d",icent);
232  fHistOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
233  fHistOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
234  fHistOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
235  fOutput->Add(fHistOversub[icent]);
236  }
237  }
238  fOutput->SetOwner(kTRUE);
239 }
240 
245 {
248 }
249 
254 {
256 
257  // Build map from all available tracks to their corresponding particle containers.
258  // This is only required for AODs which store pointers to AliVTracks instead of
259  // indices.
260  if (!fEsdMode) {
262  }
263 
264  // Run the hadronic correction
265  // loop over all clusters
266  AliVCluster *cluster = 0;
267  AliClusterContainer * clusCont = 0;
268  TIter nextClusCont(&fClusterCollArray);
269  while ((clusCont = static_cast<AliClusterContainer*>(nextClusCont()))) {
270  auto clusItCont = clusCont->accepted_momentum();
271  for (AliClusterIterableMomentumContainer::iterator clusIterator = clusItCont.begin(); clusIterator != clusItCont.end(); ++clusIterator) {
272  cluster = static_cast<AliVCluster *>(clusIterator->second);
273 
274  Double_t energyclus = 0;
275  if (fCreateHisto) {
276  fHistEbefore->Fill(fCent, cluster->GetNonLinCorrEnergy());
277  fHistNclusvsCent->Fill(fCent);
278  }
279 
280  // apply correction / subtraction
281  // to subtract only the closest track set fHadCor to a %
282  // to subtract all tracks within the cut set fHadCor to %+1
283  if (fHadCorr > 1) {
284  energyclus = ApplyHadCorrAllTracks(fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, clusIterator.current_index()), fHadCorr - 1);
285  }
286  else if (fHadCorr > 0) {
287  energyclus = ApplyHadCorrOneTrack(fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, clusIterator.current_index()), fHadCorr);
288  }
289  else {
290  energyclus = cluster->GetNonLinCorrEnergy();
291  }
292 
293  if (energyclus < 0) energyclus = 0;
294 
295  cluster->SetHadCorrEnergy(energyclus);
296 
297  if (fCreateHisto) fHistEafter->Fill(fCent, energyclus);
298 
299  }
300  }
301 
302  return kTRUE;
303 }
304 
315 {
316  fTrackToContainerMap.clear();
317 
318  // Loop by index to avoid the possibility of iterators being invalidated.
319  AliParticleContainer * partCont = 0;
320  for (int i = 0; i < fParticleCollArray.GetEntriesFast(); i++)
321  {
322  partCont = GetParticleContainer(i);
323  if (!partCont) { AliErrorStream() << "Failed to retrieve particle container at index " << i << "\n"; }
324  for (int j = 0; j < partCont->GetNParticles(); j++)
325  {
326  AliVTrack * track = static_cast<AliVTrack *>(partCont->GetParticle(j));
327  fTrackToContainerMap.insert(std::make_pair(track, partCont));
328  }
329  }
330 }
331 
336 {
337  UInt_t pbin=0;
338  if (p<0.5)
339  pbin=0;
340  else if (p>=0.5 && p<1.0)
341  pbin=1;
342  else if (p>=1.0 && p<1.5)
343  pbin=2;
344  else if (p>=1.5 && p<2.)
345  pbin=3;
346  else if (p>=2. && p<3.)
347  pbin=4;
348  else if (p>=3. && p<4.)
349  pbin=5;
350  else if (p>=4. && p<5.)
351  pbin=6;
352  else if (p>=5. && p<8.)
353  pbin=7;
354  else
355  pbin=8;
356 
357  return pbin;
358 }
359 
364 {
365  Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
366  return 2.0*EtaSigma[pbin];
367 }
368 
373 {
374  if (centbin==0) {
375  Double_t PhiMean[9]={0.036,
376  0.021,
377  0.0121,
378  0.0084,
379  0.0060,
380  0.0041,
381  0.0031,
382  0.0022,
383  0.001};
384  return PhiMean[pbin];
385  } else if (centbin==1) {
386  Double_t PhiMean[9]={0.036,
387  0.021,
388  0.0121,
389  0.0084,
390  0.0060,
391  0.0041,
392  0.0031,
393  0.0022,
394  0.001};
395  return PhiMean[pbin];
396  } else if (centbin==2) {
397  Double_t PhiMean[9]={0.036,
398  0.021,
399  0.0121,
400  0.0084,
401  0.0060,
402  0.0041,
403  0.0031,
404  0.0022,
405  0.001};
406  return PhiMean[pbin];
407  } else if (centbin==3) {
408  Double_t PhiMean[9]={0.036,
409  0.021,
410  0.0121,
411  0.0084,
412  0.0060,
413  0.0041,
414  0.0031,
415  0.0022,
416  0.001};
417  return PhiMean[pbin];
418  } else if (centbin==4) {
419  Double_t PhiMean[9]={0.036,
420  0.021,
421  0.0127,
422  0.0089,
423  0.0068,
424  0.0049,
425  0.0038,
426  0.0028,
427  0.0018};
428  return PhiMean[pbin]*(-1.);
429  } else if (centbin==5) {
430  Double_t PhiMean[9]={0.036,
431  0.021,
432  0.0127,
433  0.0089,
434  0.0068,
435  0.0048,
436  0.0038,
437  0.0028,
438  0.0018};
439  return PhiMean[pbin]*(-1.);
440  } else if (centbin==6) {
441  Double_t PhiMean[9]={0.036,
442  0.021,
443  0.0127,
444  0.0089,
445  0.0068,
446  0.0045,
447  0.0035,
448  0.0028,
449  0.0018};
450  return PhiMean[pbin]*(-1.);
451  } else if (centbin==7) {
452  Double_t PhiMean[9]={0.036,
453  0.021,
454  0.0127,
455  0.0089,
456  0.0068,
457  0.0043,
458  0.0035,
459  0.0028,
460  0.0018};
461  return PhiMean[pbin]*(-1.);
462  }
463 
464  return 0;
465 }
466 
471 {
472  if (centbin==0) {
473  Double_t PhiSigma[9]={0.0221,
474  0.0128,
475  0.0074,
476  0.0064,
477  0.0059,
478  0.0055,
479  0.0052,
480  0.0049,
481  0.0045};
482  return 2.*PhiSigma[pbin];
483  } else if (centbin==1) {
484  Double_t PhiSigma[9]={0.0217,
485  0.0120,
486  0.0076,
487  0.0066,
488  0.0062,
489  0.0058,
490  0.0054,
491  0.0054,
492  0.0045};
493  return 2.*PhiSigma[pbin];
494  } else if (centbin==2) {
495  Double_t PhiSigma[9]={0.0211,
496  0.0124,
497  0.0080,
498  0.0070,
499  0.0067,
500  0.0061,
501  0.0059,
502  0.0054,
503  0.0047};
504  return 2.*PhiSigma[pbin];
505  } else if (centbin==3) {
506  Double_t PhiSigma[9]={0.0215,
507  0.0124,
508  0.0082,
509  0.0073,
510  0.0069,
511  0.0064,
512  0.0060,
513  0.0055,
514  0.0047};
515  return 2.*PhiSigma[pbin];
516  } else if (centbin==4) {
517  Double_t PhiSigma[9]={0.0199,
518  0.0108,
519  0.0072,
520  0.0071,
521  0.0060,
522  0.0055,
523  0.0052,
524  0.0049,
525  0.0045};
526  return 2.*PhiSigma[pbin];
527  } else if (centbin==5) {
528  Double_t PhiSigma[9]={0.0200,
529  0.0110,
530  0.0074,
531  0.0071,
532  0.0064,
533  0.0059,
534  0.0055,
535  0.0052,
536  0.0045};
537  return 2.*PhiSigma[pbin];
538  } else if (centbin==6) {
539  Double_t PhiSigma[9]={0.0202,
540  0.0113,
541  0.0077,
542  0.0071,
543  0.0069,
544  0.0064,
545  0.0060,
546  0.0055,
547  0.0050};
548  return 2.*PhiSigma[pbin];
549  } else if (centbin==7) {
550  Double_t PhiSigma[9]={0.0205,
551  0.0113,
552  0.0080,
553  0.0074,
554  0.0078,
555  0.0067,
556  0.0062,
557  0.0055,
558  0.0050};
559  return 2.*PhiSigma[pbin];
560  }
561 
562  return 0;
563 }
564 
569  Double_t &totalTrkP, Int_t &Nmatches, Double_t &trkPMCfrac, Int_t &NMCmatches)
570 {
571  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
572 
573  if (!cluster) return;
574 
575  // loop over matched tracks
576  Int_t Ntrks = cluster->GetNTracksMatched();
577  for (Int_t i = 0; i < Ntrks; ++i) {
578  AliVTrack* track = 0;
579 
580  if (fEsdMode) {
581  Int_t itrack = cluster->GetTrackMatchedIndex(i);
582  if (itrack >= 0) {
584  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
585  }
586  }
587  else {
588  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(i));
589  UInt_t rejectionReason = 0;
590  AliParticleContainer * partCont = fTrackToContainerMap.at(track);
591  if (!partCont) { AliErrorStream() << "Requested particle container not available!\n"; }
592  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
593  }
594 
595  if (!track) continue;
596 
597  Double_t etadiff = 999;
598  Double_t phidiff = 999;
599  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
600  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(etadiff, phidiff);
601 
602  // check if track also points to cluster
603  if (fDoTrackClus && (track->GetEMCALcluster() != icluster)) continue;
604 
605  Double_t mom = track->P();
606  UInt_t mombin = GetMomBin(mom);
607  Int_t centbinch = fCentBin;
608 
609  if (track->Charge() < 0) centbinch += fNcentBins;
610 
611  if (fCreateHisto) {
612  Int_t etabin = 0;
613  if (track->Eta() > 0) etabin=1;
614  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
615  fHistMatchEtaPhiAll->Fill(etadiff, phidiff);
616  }
617 
618  Double_t etaCut = 0.0;
619  Double_t phiCutlo = 0.0;
620  Double_t phiCuthi = 0.0;
621 
622  if (fPhiMatch > 0) {
623  phiCutlo = -fPhiMatch;
624  phiCuthi = +fPhiMatch;
625  }
626  else {
627  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
628  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
629  }
630 
631  if (fEtaMatch > 0) {
632  etaCut = fEtaMatch;
633  }
634  else {
635  etaCut = GetEtaSigma(mombin);
636  }
637 
638  if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
639  if (track->GetLabel() > fMinMCLabel) {
640  ++NMCmatches;
641  trkPMCfrac += mom;
642  }
643  ++Nmatches;
644  totalTrkP += mom;
645 
646  if (fCreateHisto) {
647  if (fHadCorr > 1) {
648  Double_t dphi = 0;
649  Double_t deta = 0;
650  GetEtaPhiDiff(track, cluster, dphi, deta);
651  Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
652  Double_t energyclus = cluster->GetNonLinCorrEnergy();
653  fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
654  }
655  }
656  }
657  }
658 
659  if (totalTrkP > 0) trkPMCfrac /= totalTrkP;
660 }
661 
666 {
667  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
668 
669  Double_t energyclus = cluster->GetNonLinCorrEnergy();
670 
671  AliVTrack* track = 0;
672 
673  if (cluster->GetNTracksMatched() > 0) {
674  if (fEsdMode) {
675  Int_t itrack = cluster->GetTrackMatchedIndex(0);
676  if (itrack >= 0) {
678  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
679  }
680  }
681  else {
682  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
683  UInt_t rejectionReason = 0;
684  AliParticleContainer * partCont = fTrackToContainerMap.at(track);
685  if (!partCont) { AliErrorStream() << "Requested particle container not available!\n"; }
686  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
687  }
688  }
689 
690  if (!track || track->P() < 1e-6) return energyclus;
691 
692  Double_t mom = track->P();
693 
694  Double_t dEtaMin = 1e9;
695  Double_t dPhiMin = 1e9;
696  GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
697 
698  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(dEtaMin, dPhiMin);
699 
700  // check if track also points to cluster
701  Int_t cid = track->GetEMCALcluster();
702  if (fDoTrackClus && (cid != icluster)) return energyclus;
703 
704  UInt_t mombin = GetMomBin(mom);
705  Int_t centbinch = fCentBin;
706  if (track->Charge() < 0) centbinch += fNcentBins;
707 
708  // plot some histograms if switched on
709  if (fCreateHisto) {
710  Int_t etabin = 0;
711  if(track->Eta() > 0) etabin = 1;
712 
713  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
714  fHistMatchEtaPhiAll->Fill(dEtaMin, dPhiMin);
715 
716  if (mom > 0) {
717  Double_t etadiff = 0;
718  Double_t phidiff = 0;
719  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
720  Double_t dRmin = TMath::Sqrt(etadiff*etadiff + phidiff*phidiff);
721  fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
722  fHistEoPCent->Fill(fCent, energyclus / mom);
723  fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
724  }
725  }
726 
727  // define eta/phi cuts
728  Double_t etaCut = 0.0;
729  Double_t phiCutlo = 0.0;
730  Double_t phiCuthi = 0.0;
731  if (fPhiMatch > 0) {
732  phiCutlo = -fPhiMatch;
733  phiCuthi = +fPhiMatch;
734  }
735  else {
736  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
737  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
738  }
739  if (fEtaMatch > 0) {
740  etaCut = fEtaMatch;
741  }
742  else {
743  etaCut = GetEtaSigma(mombin);
744  }
745 
746  // apply the correction if the track is in the eta/phi window
747  if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
748  energyclus -= hadCorr * mom;
749  }
750 
751  return energyclus;
752 }
753 
758 {
759  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
760 
761  Double_t energyclus = cluster->GetNonLinCorrEnergy();
762  Double_t cNcells = cluster->GetNCells();
763 
764  Double_t totalTrkP = 0.0; // count total track momentum
765  Int_t Nmatches = 0; // count total number of matches
766 
767  Double_t trkPMCfrac = 0.0; // count total track momentum
768  Int_t NMCmatches = 0; // count total number of matches
769 
770  // do the loop over the matched tracks and get the number of matches and the total momentum
771  DoMatchedTracksLoop(icluster, totalTrkP, Nmatches, trkPMCfrac, NMCmatches);
772 
773  Double_t Esub = hadCorr * totalTrkP;
774 
775  if (Esub > energyclus) Esub = energyclus;
776 
777  // applying Peter's proposed algorithm
778  // never subtract the full energy of the cluster
779  Double_t clusEexcl = fEexclCell * cNcells;
780  if (energyclus < clusEexcl) clusEexcl = energyclus;
781  if ((energyclus - Esub) < clusEexcl) Esub = (energyclus - clusEexcl);
782 
783  // If enabled, use cluster M02 value to determine how to correct cluster energy.
785  Esub = ComputeM02Subtraction(cluster, energyclus, Nmatches, totalTrkP, hadCorr);
786  }
787 
788  // embedding
789  Double_t EsubMC = 0;
790  Double_t EsubBkg = 0;
791  Double_t EclusMC = 0;
792  Double_t EclusBkg = 0;
793  Double_t EclusCorr = 0;
794  Double_t EclusMCcorr = 0;
795  Double_t EclusBkgcorr = 0;
796  Double_t overSub = 0;
798  EsubMC = hadCorr * totalTrkP * trkPMCfrac;
799  EsubBkg = hadCorr * totalTrkP - EsubMC;
800  EclusMC = energyclus * cluster->GetMCEnergyFraction();
801  EclusBkg = energyclus - EclusMC;
802 
803  if (energyclus > Esub)
804  EclusCorr = energyclus - Esub;
805 
806  if (EclusMC > EsubMC)
807  EclusMCcorr = EclusMC - EsubMC;
808 
809  if (EclusBkg > EsubBkg)
810  EclusBkgcorr = EclusBkg - EsubBkg;
811 
812  overSub = EclusMCcorr + EclusBkgcorr - EclusCorr;
813  }
814 
815  // plot some histograms if switched on
816  if (fCreateHisto) {
817  fHistNMatchCent->Fill(fCent, Nmatches);
818  fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
819 
820  if (Nmatches > 0) fHistNclusMatchvsCent->Fill(fCent);
821 
822  if (Nmatches < 3) {
823  fHistNCellsEnergy[fCentBin][Nmatches]->Fill(energyclus, cNcells);
824  }
825  else {
826  fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
827  }
828 
829  if (totalTrkP > 0) {
830  Double_t EoP = energyclus / totalTrkP;
831  fHistEoPCent->Fill(fCent, EoP);
832  fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
833 
834  fHistEsubPchRatAll[fCentBin]->Fill(totalTrkP, Esub / totalTrkP);
835 
836  if (Nmatches == 1) {
837  AliVTrack* track = 0;
838  if (fEsdMode) {
839  Int_t itrack = cluster->GetTrackMatchedIndex(0);
840  if (itrack >= 0) {
842  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
843  }
844  }
845  else {
846  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
847  UInt_t rejectionReason = 0;
848  AliParticleContainer * partCont = fTrackToContainerMap.at(track);
849  if (!partCont) { AliErrorStream() << "Requested particle container not available!\n"; }
850  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
851  }
852  if (track) {
853  Int_t centbinchm = fCentBin;
854  if (track->Charge() < 0) centbinchm += fNcentBins;
855  fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
856  fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
857  }
858  }
859 
861  fHistOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
862 
863  if (cluster->GetMCEnergyFraction() > 0.95)
864  fHistOversubMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
865  else if (cluster->GetMCEnergyFraction() < 0.05)
866  fHistOversubNonMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
867 
868  if (trkPMCfrac < 0.05)
869  fHistNonEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
870  else if (trkPMCfrac > 0.95)
871  fHistEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
872  }
873  }
874  }
875 
877  Esub -= overSub;
878  if (EclusBkgcorr + EclusMCcorr > 0) {
879  Double_t newfrac = EclusMCcorr / (EclusBkgcorr + EclusMCcorr);
880  cluster->SetMCEnergyFraction(newfrac);
881  }
882  }
883 
884  // apply the correction
885  energyclus -= Esub;
886 
887  return energyclus;
888 }
889 
893 Double_t AliEmcalCorrectionClusterHadronicCorrection::ComputeM02Subtraction(const AliVCluster* cluster, Double_t energyclus, Int_t Nmatches, Double_t totalTrkP, Double_t hadCorr)
894 {
895  Double_t Esub = 0.;
896  Double_t clusM02 = cluster->GetM02();
897 
898  // Note: The comments below give only a rough indication of what particle types appear in each selection.
899 
900  // For M02 in the single photon region, the signal is primarily: Single photons, single electrons, single MIPs
901  if (clusM02 > 0.1 && clusM02 < 0.4) {
902  if (Nmatches == 0) { // Single photon (modulo tracking efficiency)
903  Esub = 0;
904  }
905  else { // Single electron, single MIP
906  Esub = energyclus;
907  }
908  }
909 
910  // For large M02, the signal is primarily: Single hadronic shower, photon-photon overlap, photon-MIP overlap, MIP-MIP overlap,
911  // MIP-hadronic shower overlap, hadronic shower - hadronic shower overlap)
912  if (clusM02 > 0.4) {
913  if (Nmatches == 0) { // Single neutral hadronic shower, photon-photon overlap, photon-neutral had shower overlap, neutral had shower overlap (modulo tracking efficiency)
914  Esub = 0;
915  }
916  else if (Nmatches == 1) { // Single charged hadronic shower, photon-MIP overlap, MIP-MIP overlap, MIP-had shower overlap, had shower-shower overlap
919  }
920  else {
921  Esub = hadCorr * totalTrkP;
922  }
923  }
924  else if (Nmatches > 1) { // MIP-MIP overlap, had shower-shower overlap
925  Esub = energyclus;
926  }
927  }
928 
929  return Esub;
930 }
Int_t fNcentBins
How many centrality bins (this member copied from AliAnalysisTaskEmcal)
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...