AliPhysics  ec7afe5 (ec7afe5)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliEmcalCorrectionClusterHadronicCorrection.cxx
Go to the documentation of this file.
1 // AliEmcalCorrectionClusterHadronicCorrection
2 //
3 
5 
6 #include <TH2.h>
7 #include <TList.h>
8 
9 #include "AliClusterContainer.h"
10 #include "AliParticleContainer.h"
11 #include "AliVTrack.h"
12 
16 
17 // Actually registers the class with the base class
19 
24  AliEmcalCorrectionComponent("AliEmcalCorrectionClusterHadronicCorrection"),
25  fPhiMatch(0.05),
26  fEtaMatch(0.025),
27  fDoTrackClus(kTRUE),
28  fHadCorr(0),
29  fEexclCell(0),
30  fDoExact(kFALSE),
31  fClusterContainerIndexMap(),
32  fParticleContainerIndexMap(),
33  fHistMatchEtaPhiAll(0),
34  fHistMatchEtaPhiAllCl(0),
35  fHistNclusvsCent(0),
36  fHistNclusMatchvsCent(0),
37  fHistEbefore(0),
38  fHistEafter(0),
39  fHistEoPCent(0),
40  fHistNMatchCent(0),
41  fHistNClusMatchCent(0)
42 {
43  for(Int_t i=0; i<10; i++) {
44  fHistEsubPch[i] = 0;
45  fHistEsubPchRat[i] = 0;
46 
47  if (i<5) {
48  fHistEsubPchRatAll[i] = 0;
49 
50  fHistMatchEvsP[i] = 0;
51  fHistMatchdRvsEP[i] = 0;
52  fHistNMatchEnergy[i] = 0;
53 
58  fHistOversub[i] = 0;
59 
60  for(Int_t j=0; j<4; j++)
61  fHistNCellsEnergy[i][j] = 0;
62  }
63 
64  for(Int_t j=0; j<9; j++) {
65  for(Int_t k=0; k<2; k++)
66  fHistMatchEtaPhi[i][j][k] = 0;
67  }
68  }
69 }
70 
75 {
76 }
77 
82 {
83  // Initialization
85 
86  GetProperty("createHistos", fCreateHisto);
87  GetProperty("phiMatch", fPhiMatch);
88  GetProperty("etaMatch", fEtaMatch);
89  GetProperty("hadCorr", fHadCorr);
90  GetProperty("Eexcl", fEexclCell);
91  GetProperty("doTrackClus", fDoTrackClus);
92 
93  if (!fEsdMode && fParticleCollArray.GetEntries() > 1) {
94  AliWarning("================================================================================");
95  AliWarning("== Added multiple particle containers when running with AOD!");
96  AliWarning("== Particle selection of the first particle container will be applied");
97  AliWarning("== to _ALL_ particles! If you need a different selection, then change");
98  AliWarning("== the order of adding the containers so that the desired container is first!");
99  AliWarning("================================================================================");
100  }
101 
102  return kTRUE;
103 }
104 
109 {
111 
112  // Create my user objects.
113 
114  if (!fCreateHisto) return;
115 
116  fHistMatchEtaPhiAll = new TH2F("fHistMatchEtaPhiAll", "fHistMatchEtaPhiAll;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
118 
119  fHistMatchEtaPhiAllCl = new TH2F("fHistMatchEtaPhiAllCl", "fHistMatchEtaPhiAllCl;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
121 
122  TString name;
123  TString temp;
124  const Int_t nCentChBins = fNcentBins * 2;
125  for(Int_t icent=0; icent<nCentChBins; ++icent) {
126  for(Int_t ipt=0; ipt<9; ++ipt) {
127  for(Int_t ieta=0; ieta<2; ++ieta) {
128  name = Form("fHistMatchEtaPhi_%i_%i_%i",icent,ipt,ieta);
129  fHistMatchEtaPhi[icent][ipt][ieta] = new TH2F(name, name, fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
130  fHistMatchEtaPhi[icent][ipt][ieta]->SetXTitle("#Delta#eta");
131  fHistMatchEtaPhi[icent][ipt][ieta]->SetYTitle("#Delta#phi");
132  fOutput->Add(fHistMatchEtaPhi[icent][ipt][ieta]);
133  }
134  }
135 
136  name = Form("fHistEsubPch_%i",icent);
137  temp = Form("%s (Nmatches==1)",name.Data());
138  fHistEsubPch[icent]=new TH1F(name, temp, fNbins, fMinBinPt, fMaxBinPt);
139  fHistEsubPch[icent]->SetXTitle("#sum p (GeV) weighted with E_{sub}");
140  fOutput->Add(fHistEsubPch[icent]);
141 
142  name = Form("fHistEsubPchRat_%i",icent);
143  temp = Form("%s (Nmatches==1)",name.Data());
144  fHistEsubPchRat[icent]=new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
145  fHistEsubPchRat[icent]->SetXTitle("#Sigma p (GeV)");
146  fHistEsubPchRat[icent]->SetYTitle("E_{sub} / #sum p");
147  fOutput->Add(fHistEsubPchRat[icent]);
148 
149  if (icent<fNcentBins) {
150  for(Int_t itrk=0; itrk<4; ++itrk) {
151  name = Form("fHistNCellsEnergy_%i_%i",icent,itrk);
152  temp = Form("%s (Nmatches==%d);N_{cells};E_{clus} (GeV)",name.Data(),itrk);
153  fHistNCellsEnergy[icent][itrk] = new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, 101, -0.5, 100.5);
154  fOutput->Add(fHistNCellsEnergy[icent][itrk]);
155  }
156 
157  name = Form("fHistEsubPchRatAll_%i",icent);
158  temp = Form("%s (all Nmatches)",name.Data());
159  fHistEsubPchRatAll[icent]=new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
160  fHistEsubPchRatAll[icent]->SetXTitle("#Sigma p (GeV)");
161  fHistEsubPchRatAll[icent]->SetYTitle("E_{sub} / #sum p");
162  fOutput->Add(fHistEsubPchRatAll[icent]);
163 
164  name = Form("fHistMatchEvsP_%i",icent);
165  temp = Form("%s (all Nmatches)",name.Data());
166  fHistMatchEvsP[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
167  fHistMatchEvsP[icent]->SetXTitle("E_{clus} (GeV)");
168  fHistMatchEvsP[icent]->SetYTitle("E_{clus} / #sum p");
169  fOutput->Add(fHistMatchEvsP[icent]);
170 
171  name = Form("fHistMatchdRvsEP_%i",icent);
172  temp = Form("%s (all Nmatches)",name.Data());
173  fHistMatchdRvsEP[icent] = new TH2F(name, temp, fNbins, 0., 0.2, fNbins*2, 0., 10.);
174  fHistMatchdRvsEP[icent]->SetXTitle("#Delta R between track and cluster");
175  fHistMatchdRvsEP[icent]->SetYTitle("E_{clus} / p");
176  fOutput->Add(fHistMatchdRvsEP[icent]);
177 
178  name = Form("fHistNMatchEnergy_%i",icent);
179  fHistNMatchEnergy[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, 101, -0.5, 100.5);
180  fHistNMatchEnergy[icent]->SetXTitle("E_{clus} (GeV)");
181  fHistNMatchEnergy[icent]->SetYTitle("N_{matches}");
182  fOutput->Add(fHistNMatchEnergy[icent]);
183  }
184  }
185 
186  fHistNclusvsCent = new TH1F("Nclusvscent", "NclusVsCent; Cent (%)", 100, 0, 100);
187  fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent (all Nmatches); Cent (%)", 100, 0, 100);
188  fHistEbefore = new TH1F("Ebefore", "Ebefore; Cent (%); E_{clus} (GeV)", 100, 0, 100);
189  fHistEafter = new TH1F("Eafter", "Eafter; Cent (%); E_{clus} (GeV)", 100, 0, 100);
190  fHistEoPCent = new TH2F("EoPCent", "EoPCent; Cent (%); E_{clus} / #sum p", 100, 0, 100, fNbins*2, 0, 10);
191  fHistNMatchCent = new TH2F("NMatchesCent", "NMatchesCent; Cent (%); Nmatches", 100, 0, 100, 11, -0.5, 10.5);
192  fHistNClusMatchCent = new TH2F("NClusMatchesCent", "NClusMatchesCent; Cent (%); Nmatches", 100, 0, 100, 11, -0.5, 10.5);
193 
196  fOutput->Add(fHistEbefore);
197  fOutput->Add(fHistEafter);
198  fOutput->Add(fHistEoPCent);
199  fOutput->Add(fHistNMatchCent);
201 
202  if (fIsEmbedded) {
203  for(Int_t icent=0; icent<fNcentBins; ++icent) {
204  name = Form("fHistEmbTrackMatchesOversub_%d",icent);
205  fHistEmbTrackMatchesOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
206  fHistEmbTrackMatchesOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
207  fHistEmbTrackMatchesOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
209 
210  name = Form("fHistNonEmbTrackMatchesOversub_%d",icent);
211  fHistNonEmbTrackMatchesOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
212  fHistNonEmbTrackMatchesOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
213  fHistNonEmbTrackMatchesOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
215 
216  name = Form("fHistOversubMCClusters_%d",icent);
217  fHistOversubMCClusters[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
218  fHistOversubMCClusters[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
219  fHistOversubMCClusters[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
220  fOutput->Add(fHistOversubMCClusters[icent]);
221 
222  name = Form("fHistOversubNonMCClusters_%d",icent);
223  fHistOversubNonMCClusters[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
224  fHistOversubNonMCClusters[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
225  fHistOversubNonMCClusters[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
226  fOutput->Add(fHistOversubNonMCClusters[icent]);
227 
228  name = Form("fHistOversub_%d",icent);
229  fHistOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
230  fHistOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
231  fHistOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
232  fOutput->Add(fHistOversub[icent]);
233  }
234  }
235  fOutput->SetOwner(kTRUE);
236 }
237 
242 {
245 }
246 
251 {
253 
254  // Run the hadronic correction
255  // loop over all clusters
256  AliVCluster *cluster = 0;
257  AliClusterContainer * clusCont = 0;
258  TIter nextClusCont(&fClusterCollArray);
259  while ((clusCont = static_cast<AliClusterContainer*>(nextClusCont()))) {
260  auto clusItCont = clusCont->accepted_momentum();
261  for (AliClusterIterableMomentumContainer::iterator clusIterator = clusItCont.begin(); clusIterator != clusItCont.end(); ++clusIterator) {
262  cluster = static_cast<AliVCluster *>(clusIterator->second);
263 
264  Double_t energyclus = 0;
265  if (fCreateHisto) {
266  fHistEbefore->Fill(fCent, cluster->GetNonLinCorrEnergy());
267  fHistNclusvsCent->Fill(fCent);
268  }
269 
270  // apply correction / subtraction
271  // to subtract only the closest track set fHadCor to a %
272  // to subtract all tracks within the cut set fHadCor to %+1
273  if (fHadCorr > 1) {
274  energyclus = ApplyHadCorrAllTracks(fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, clusIterator.current_index()), fHadCorr - 1);
275  }
276  else if (fHadCorr > 0) {
277  energyclus = ApplyHadCorrOneTrack(fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, clusIterator.current_index()), fHadCorr);
278  }
279  else {
280  energyclus = cluster->GetNonLinCorrEnergy();
281  }
282 
283  if (energyclus < 0) energyclus = 0;
284 
285  cluster->SetHadCorrEnergy(energyclus);
286 
287  if (fCreateHisto) fHistEafter->Fill(fCent, energyclus);
288 
289  }
290  }
291 
292  return kTRUE;
293 }
294 
299 {
300  UInt_t pbin=0;
301  if (p<0.5)
302  pbin=0;
303  else if (p>=0.5 && p<1.0)
304  pbin=1;
305  else if (p>=1.0 && p<1.5)
306  pbin=2;
307  else if (p>=1.5 && p<2.)
308  pbin=3;
309  else if (p>=2. && p<3.)
310  pbin=4;
311  else if (p>=3. && p<4.)
312  pbin=5;
313  else if (p>=4. && p<5.)
314  pbin=6;
315  else if (p>=5. && p<8.)
316  pbin=7;
317  else
318  pbin=8;
319 
320  return pbin;
321 }
322 
327 {
328  Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
329  return 2.0*EtaSigma[pbin];
330 }
331 
336 {
337  if (centbin==0) {
338  Double_t PhiMean[9]={0.036,
339  0.021,
340  0.0121,
341  0.0084,
342  0.0060,
343  0.0041,
344  0.0031,
345  0.0022,
346  0.001};
347  return PhiMean[pbin];
348  } else if (centbin==1) {
349  Double_t PhiMean[9]={0.036,
350  0.021,
351  0.0121,
352  0.0084,
353  0.0060,
354  0.0041,
355  0.0031,
356  0.0022,
357  0.001};
358  return PhiMean[pbin];
359  } else if (centbin==2) {
360  Double_t PhiMean[9]={0.036,
361  0.021,
362  0.0121,
363  0.0084,
364  0.0060,
365  0.0041,
366  0.0031,
367  0.0022,
368  0.001};
369  return PhiMean[pbin];
370  } else if (centbin==3) {
371  Double_t PhiMean[9]={0.036,
372  0.021,
373  0.0121,
374  0.0084,
375  0.0060,
376  0.0041,
377  0.0031,
378  0.0022,
379  0.001};
380  return PhiMean[pbin];
381  } else if (centbin==4) {
382  Double_t PhiMean[9]={0.036,
383  0.021,
384  0.0127,
385  0.0089,
386  0.0068,
387  0.0049,
388  0.0038,
389  0.0028,
390  0.0018};
391  return PhiMean[pbin]*(-1.);
392  } else if (centbin==5) {
393  Double_t PhiMean[9]={0.036,
394  0.021,
395  0.0127,
396  0.0089,
397  0.0068,
398  0.0048,
399  0.0038,
400  0.0028,
401  0.0018};
402  return PhiMean[pbin]*(-1.);
403  } else if (centbin==6) {
404  Double_t PhiMean[9]={0.036,
405  0.021,
406  0.0127,
407  0.0089,
408  0.0068,
409  0.0045,
410  0.0035,
411  0.0028,
412  0.0018};
413  return PhiMean[pbin]*(-1.);
414  } else if (centbin==7) {
415  Double_t PhiMean[9]={0.036,
416  0.021,
417  0.0127,
418  0.0089,
419  0.0068,
420  0.0043,
421  0.0035,
422  0.0028,
423  0.0018};
424  return PhiMean[pbin]*(-1.);
425  }
426 
427  return 0;
428 }
429 
434 {
435  if (centbin==0) {
436  Double_t PhiSigma[9]={0.0221,
437  0.0128,
438  0.0074,
439  0.0064,
440  0.0059,
441  0.0055,
442  0.0052,
443  0.0049,
444  0.0045};
445  return 2.*PhiSigma[pbin];
446  } else if (centbin==1) {
447  Double_t PhiSigma[9]={0.0217,
448  0.0120,
449  0.0076,
450  0.0066,
451  0.0062,
452  0.0058,
453  0.0054,
454  0.0054,
455  0.0045};
456  return 2.*PhiSigma[pbin];
457  } else if (centbin==2) {
458  Double_t PhiSigma[9]={0.0211,
459  0.0124,
460  0.0080,
461  0.0070,
462  0.0067,
463  0.0061,
464  0.0059,
465  0.0054,
466  0.0047};
467  return 2.*PhiSigma[pbin];
468  } else if (centbin==3) {
469  Double_t PhiSigma[9]={0.0215,
470  0.0124,
471  0.0082,
472  0.0073,
473  0.0069,
474  0.0064,
475  0.0060,
476  0.0055,
477  0.0047};
478  return 2.*PhiSigma[pbin];
479  } else if (centbin==4) {
480  Double_t PhiSigma[9]={0.0199,
481  0.0108,
482  0.0072,
483  0.0071,
484  0.0060,
485  0.0055,
486  0.0052,
487  0.0049,
488  0.0045};
489  return 2.*PhiSigma[pbin];
490  } else if (centbin==5) {
491  Double_t PhiSigma[9]={0.0200,
492  0.0110,
493  0.0074,
494  0.0071,
495  0.0064,
496  0.0059,
497  0.0055,
498  0.0052,
499  0.0045};
500  return 2.*PhiSigma[pbin];
501  } else if (centbin==6) {
502  Double_t PhiSigma[9]={0.0202,
503  0.0113,
504  0.0077,
505  0.0071,
506  0.0069,
507  0.0064,
508  0.0060,
509  0.0055,
510  0.0050};
511  return 2.*PhiSigma[pbin];
512  } else if (centbin==7) {
513  Double_t PhiSigma[9]={0.0205,
514  0.0113,
515  0.0080,
516  0.0074,
517  0.0078,
518  0.0067,
519  0.0062,
520  0.0055,
521  0.0050};
522  return 2.*PhiSigma[pbin];
523  }
524 
525  return 0;
526 }
527 
532  Double_t &totalTrkP, Int_t &Nmatches, Double_t &trkPMCfrac, Int_t &NMCmatches)
533 {
534  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
535 
536  if (!cluster) return;
537 
538  // loop over matched tracks
539  Int_t Ntrks = Ntrks = cluster->GetNTracksMatched();
540  for (Int_t i = 0; i < Ntrks; ++i) {
541  AliVTrack* track = 0;
542 
543  if (fEsdMode) {
544  Int_t itrack = cluster->GetTrackMatchedIndex(i);
545  if (itrack >= 0) {
547  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
548  }
549  }
550  else {
551  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(i));
552  UInt_t rejectionReason = 0;
554  if (!partCont) { AliError("No particle container available!"); }
555  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
556  }
557 
558  if (!track) continue;
559 
560  Double_t etadiff = 999;
561  Double_t phidiff = 999;
562  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
563  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(etadiff, phidiff);
564 
565  // check if track also points to cluster
566  if (fDoTrackClus && (track->GetEMCALcluster() != icluster)) continue;
567 
568  Double_t mom = track->P();
569  UInt_t mombin = GetMomBin(mom);
570  Int_t centbinch = fCentBin;
571 
572  if (track->Charge() < 0) centbinch += fNcentBins;
573 
574  if (fCreateHisto) {
575  Int_t etabin = 0;
576  if (track->Eta() > 0) etabin=1;
577  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
578  fHistMatchEtaPhiAll->Fill(etadiff, phidiff);
579  }
580 
581  Double_t etaCut = 0.0;
582  Double_t phiCutlo = 0.0;
583  Double_t phiCuthi = 0.0;
584 
585  if (fPhiMatch > 0) {
586  phiCutlo = -fPhiMatch;
587  phiCuthi = +fPhiMatch;
588  }
589  else {
590  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
591  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
592  }
593 
594  if (fEtaMatch > 0) {
595  etaCut = fEtaMatch;
596  }
597  else {
598  etaCut = GetEtaSigma(mombin);
599  }
600 
601  if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
602  if (track->GetLabel() > fMinMCLabel) {
603  ++NMCmatches;
604  trkPMCfrac += mom;
605  }
606  ++Nmatches;
607  totalTrkP += mom;
608 
609  if (fCreateHisto) {
610  if (fHadCorr > 1) {
611  Double_t dphi = 0;
612  Double_t deta = 0;
613  GetEtaPhiDiff(track, cluster, dphi, deta);
614  Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
615  Double_t energyclus = cluster->GetNonLinCorrEnergy();
616  fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
617  }
618  }
619  }
620  }
621 
622  if (totalTrkP > 0) trkPMCfrac /= totalTrkP;
623 }
624 
629 {
630  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
631 
632  Double_t energyclus = cluster->GetNonLinCorrEnergy();
633 
634  AliVTrack* track = 0;
635 
636  if (cluster->GetNTracksMatched() > 0) {
637  if (fEsdMode) {
638  Int_t itrack = cluster->GetTrackMatchedIndex(0);
639  if (itrack >= 0) {
641  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
642  }
643  }
644  else {
645  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
646  UInt_t rejectionReason = 0;
648  if (!partCont) { AliError("No particle container available!"); }
649  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
650  }
651  }
652 
653  if (!track || track->P() < 1e-6) return energyclus;
654 
655  Double_t mom = track->P();
656 
657  Double_t dEtaMin = 1e9;
658  Double_t dPhiMin = 1e9;
659  GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
660 
661  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(dEtaMin, dPhiMin);
662 
663  // check if track also points to cluster
664  Int_t cid = track->GetEMCALcluster();
665  if (fDoTrackClus && (cid != icluster)) return energyclus;
666 
667  UInt_t mombin = GetMomBin(mom);
668  Int_t centbinch = fCentBin;
669  if (track->Charge() < 0) centbinch += fNcentBins;
670 
671  // plot some histograms if switched on
672  if (fCreateHisto) {
673  Int_t etabin = 0;
674  if(track->Eta() > 0) etabin = 1;
675 
676  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
677  fHistMatchEtaPhiAll->Fill(dEtaMin, dPhiMin);
678 
679  if (mom > 0) {
680  Double_t etadiff = 0;
681  Double_t phidiff = 0;
682  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
683  Double_t dRmin = TMath::Sqrt(etadiff*etadiff + phidiff*phidiff);
684  fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
685  fHistEoPCent->Fill(fCent, energyclus / mom);
686  fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
687  }
688  }
689 
690  // define eta/phi cuts
691  Double_t etaCut = 0.0;
692  Double_t phiCutlo = 0.0;
693  Double_t phiCuthi = 0.0;
694  if (fPhiMatch > 0) {
695  phiCutlo = -fPhiMatch;
696  phiCuthi = +fPhiMatch;
697  }
698  else {
699  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
700  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
701  }
702  if (fEtaMatch > 0) {
703  etaCut = fEtaMatch;
704  }
705  else {
706  etaCut = GetEtaSigma(mombin);
707  }
708 
709  // apply the correction if the track is in the eta/phi window
710  if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
711  energyclus -= hadCorr * mom;
712  }
713 
714  return energyclus;
715 }
716 
721 {
722  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
723 
724  Double_t energyclus = cluster->GetNonLinCorrEnergy();
725  Double_t cNcells = cluster->GetNCells();
726 
727  Double_t totalTrkP = 0.0; // count total track momentum
728  Int_t Nmatches = 0; // count total number of matches
729 
730  Double_t trkPMCfrac = 0.0; // count total track momentum
731  Int_t NMCmatches = 0; // count total number of matches
732 
733  // do the loop over the matched tracks and get the number of matches and the total momentum
734  DoMatchedTracksLoop(icluster, totalTrkP, Nmatches, trkPMCfrac, NMCmatches);
735 
736  Double_t Esub = hadCorr * totalTrkP;
737 
738  if (Esub > energyclus) Esub = energyclus;
739 
740  // applying Peter's proposed algorithm
741  // never subtract the full energy of the cluster
742  Double_t clusEexcl = fEexclCell * cNcells;
743  if (energyclus < clusEexcl) clusEexcl = energyclus;
744  if ((energyclus - Esub) < clusEexcl) Esub = (energyclus - clusEexcl);
745 
746  // embedding
747  Double_t EsubMC = 0;
748  Double_t EsubBkg = 0;
749  Double_t EclusMC = 0;
750  Double_t EclusBkg = 0;
751  Double_t EclusCorr = 0;
752  Double_t EclusMCcorr = 0;
753  Double_t EclusBkgcorr = 0;
754  Double_t overSub = 0;
755  if (fIsEmbedded) {
756  EsubMC = hadCorr * totalTrkP * trkPMCfrac;
757  EsubBkg = hadCorr * totalTrkP - EsubMC;
758  EclusMC = energyclus * cluster->GetMCEnergyFraction();
759  EclusBkg = energyclus - EclusMC;
760 
761  if (energyclus > Esub)
762  EclusCorr = energyclus - Esub;
763 
764  if (EclusMC > EsubMC)
765  EclusMCcorr = EclusMC - EsubMC;
766 
767  if (EclusBkg > EsubBkg)
768  EclusBkgcorr = EclusBkg - EsubBkg;
769 
770  overSub = EclusMCcorr + EclusBkgcorr - EclusCorr;
771  }
772 
773  // plot some histograms if switched on
774  if (fCreateHisto) {
775  fHistNMatchCent->Fill(fCent, Nmatches);
776  fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
777 
778  if (Nmatches > 0) fHistNclusMatchvsCent->Fill(fCent);
779 
780  if (Nmatches < 3) {
781  fHistNCellsEnergy[fCentBin][Nmatches]->Fill(energyclus, cNcells);
782  }
783  else {
784  fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
785  }
786 
787  if (totalTrkP > 0) {
788  Double_t EoP = energyclus / totalTrkP;
789  fHistEoPCent->Fill(fCent, EoP);
790  fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
791 
792  fHistEsubPchRatAll[fCentBin]->Fill(totalTrkP, Esub / totalTrkP);
793 
794  if (Nmatches == 1) {
795  AliVTrack* track = 0;
796  if (fEsdMode) {
797  Int_t itrack = cluster->GetTrackMatchedIndex(0);
798  if (itrack >= 0) {
800  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
801  }
802  }
803  else {
804  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
805  UInt_t rejectionReason = 0;
807  if (!partCont) { AliError("No particle container available!"); }
808  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
809  }
810  if (track) {
811  Int_t centbinchm = fCentBin;
812  if (track->Charge() < 0) centbinchm += fNcentBins;
813  fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
814  fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
815  }
816  }
817 
818  if (fIsEmbedded) {
819  fHistOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
820 
821  if (cluster->GetMCEnergyFraction() > 0.95)
822  fHistOversubMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
823  else if (cluster->GetMCEnergyFraction() < 0.05)
824  fHistOversubNonMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
825 
826  if (trkPMCfrac < 0.05)
827  fHistNonEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
828  else if (trkPMCfrac > 0.95)
829  fHistEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
830  }
831  }
832  }
833 
834  if (fIsEmbedded && fDoExact) {
835  Esub -= overSub;
836  if (EclusBkgcorr + EclusMCcorr > 0) {
837  Double_t newfrac = EclusMCcorr / (EclusBkgcorr + EclusMCcorr);
838  cluster->SetMCEnergyFraction(newfrac);
839  }
840  }
841 
842  // apply the correction
843  energyclus -= Esub;
844 
845  return energyclus;
846 }
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%)
static RegisterCorrectionComponent< AliEmcalCorrectionClusterHadronicCorrection > reg
TH2 * fHistOversubNonMCClusters[5]
!Over-subtracted energy / cluster energy (cluster MC energy fraction < 5%)
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) ...
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.
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.
TList * fOutput
! List of output histograms
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)
ClassImp(AliAnalysisTaskCRC) AliAnalysisTaskCRC
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
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.
Bool_t fIsEmbedded
Trigger, embedded signal.