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