AliPhysics  cda3415 (cda3415)
 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 
20 //________________________________________________________________________
22  AliEmcalCorrectionComponent("AliEmcalCorrectionClusterHadronicCorrection"),
23  fPhiMatch(0.05),
24  fEtaMatch(0.025),
25  fDoTrackClus(kTRUE),
26  fHadCorr(0),
27  fEexclCell(0),
28  fDoExact(kFALSE),
29  fHistMatchEtaPhiAll(0),
30  fHistMatchEtaPhiAllTr(0),
31  fHistMatchEtaPhiAllCl(0),
32  fHistNclusvsCent(0),
33  fHistNclusMatchvsCent(0),
34  fHistEbefore(0),
35  fHistEafter(0),
36  fHistEoPCent(0),
37  fHistNMatchCent(0),
38  fHistNClusMatchCent(0)
39 {
40  // Default constructor
41  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
42 
43  for(Int_t i=0; i<8; i++) {
44  fHistEsubPch[i] = 0;
45  fHistEsubPchRat[i] = 0;
46  fHistEsubPchRatAll[i] = 0;
47 
48  if (i<4) {
49  fHistMatchEvsP[i] = 0;
50  fHistMatchdRvsEP[i] = 0;
51  fHistNMatchEnergy[i] = 0;
52 
57  fHistOversub[i] = 0;
58 
59  for(Int_t j=0; j<4; j++)
60  fHistNCellsEnergy[i][j] = 0;
61  }
62 
63  for(Int_t j=0; j<9; j++) {
64  for(Int_t k=0; k<2; k++)
65  fHistMatchEtaPhi[i][j][k] = 0;
66  }
67  }
68 }
69 
70 //________________________________________________________________________
72 {
73  // Destructor
74 }
75 
76 //________________________________________________________________________
78 {
79  // Initialization
80  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
82  // Do base class initializations and if it fails -> bail out
83  //AliAnalysisTaskEmcal::ExecOnce();
84  //if (!fInitialized) return;
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  return kTRUE;
94 }
95 
96 //________________________________________________________________________
98 {
99  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
101 
102  // Create my user objects.
103 
104  if (!fCreateHisto) return;
105 
106  fHistMatchEtaPhiAll = new TH2F("fHistMatchEtaPhiAll", "fHistMatchEtaPhiAll;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
108 
109  fHistMatchEtaPhiAllTr = new TH2F("fHistMatchEtaPhiAllTr", "fHistMatchEtaPhiAllTr;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
111 
112  fHistMatchEtaPhiAllCl = new TH2F("fHistMatchEtaPhiAllCl", "fHistMatchEtaPhiAllCl;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
114 
115  TString name;
116  TString temp;
117  const Int_t nCentChBins = fNcentBins * 2;
118  for(Int_t icent=0; icent<nCentChBins; ++icent) {
119  for(Int_t ipt=0; ipt<9; ++ipt) {
120  for(Int_t ieta=0; ieta<2; ++ieta) {
121  name = Form("fHistMatchEtaPhi_%i_%i_%i",icent,ipt,ieta);
122  fHistMatchEtaPhi[icent][ipt][ieta] = new TH2F(name, name, fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
123  fHistMatchEtaPhi[icent][ipt][ieta]->SetXTitle("#Delta#eta");
124  fHistMatchEtaPhi[icent][ipt][ieta]->SetYTitle("#Delta#phi");
125  fOutput->Add(fHistMatchEtaPhi[icent][ipt][ieta]);
126  }
127  }
128 
129  name = Form("fHistEsubPch_%i",icent);
130  temp = Form("%s (Nmatches==1)",name.Data());
131  fHistEsubPch[icent]=new TH1F(name, temp, fNbins, fMinBinPt, fMaxBinPt);
132  fHistEsubPch[icent]->SetXTitle("#sum p (GeV) weighted with E_{sub}");
133  fOutput->Add(fHistEsubPch[icent]);
134 
135  name = Form("fHistEsubPchRat_%i",icent);
136  temp = Form("%s (Nmatches==1)",name.Data());
137  fHistEsubPchRat[icent]=new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
138  fHistEsubPchRat[icent]->SetXTitle("#Sigma p (GeV)");
139  fHistEsubPchRat[icent]->SetYTitle("E_{sub} / #sum p");
140  fOutput->Add(fHistEsubPchRat[icent]);
141 
142  name = Form("fHistEsubPchRatAll_%i",icent);
143  temp = Form("%s (all Nmatches)",name.Data());
144  fHistEsubPchRatAll[icent]=new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
145  fHistEsubPchRatAll[icent]->SetXTitle("#Sigma p (GeV)");
146  fHistEsubPchRatAll[icent]->SetYTitle("E_{sub} / #sum p");
147  fOutput->Add(fHistEsubPchRatAll[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("fHistMatchEvsP_%i",icent);
158  temp = Form("%s (all Nmatches)",name.Data());
159  fHistMatchEvsP[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
160  fHistMatchEvsP[icent]->SetXTitle("E_{clus} (GeV)");
161  fHistMatchEvsP[icent]->SetYTitle("E_{clus} / #sum p");
162  fOutput->Add(fHistMatchEvsP[icent]);
163 
164  name = Form("fHistMatchdRvsEP_%i",icent);
165  temp = Form("%s (all Nmatches)",name.Data());
166  fHistMatchdRvsEP[icent] = new TH2F(name, temp, fNbins, 0., 0.2, fNbins*2, 0., 10.);
167  fHistMatchdRvsEP[icent]->SetXTitle("#Delta R between track and cluster");
168  fHistMatchdRvsEP[icent]->SetYTitle("E_{clus} / p");
169  fOutput->Add(fHistMatchdRvsEP[icent]);
170 
171  name = Form("fHistNMatchEnergy_%i",icent);
172  fHistNMatchEnergy[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, 101, -0.5, 100.5);
173  fHistNMatchEnergy[icent]->SetXTitle("E_{clus} (GeV)");
174  fHistNMatchEnergy[icent]->SetYTitle("N_{matches}");
175  fOutput->Add(fHistNMatchEnergy[icent]);
176  }
177  }
178 
179  fHistNclusvsCent = new TH1F("Nclusvscent", "NclusVsCent; Cent (%)", 100, 0, 100);
180  fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent (all Nmatches); Cent (%)", 100, 0, 100);
181  fHistEbefore = new TH1F("Ebefore", "Ebefore; Cent (%); E_{clus} (GeV)", 100, 0, 100);
182  fHistEafter = new TH1F("Eafter", "Eafter; Cent (%); E_{clus} (GeV)", 100, 0, 100);
183  fHistEoPCent = new TH2F("EoPCent", "EoPCent; Cent (%); E_{clus} / #sum p", 100, 0, 100, fNbins*2, 0, 10);
184  fHistNMatchCent = new TH2F("NMatchesCent", "NMatchesCent; Cent (%); Nmatches", 100, 0, 100, 11, -0.5, 10.5);
185  fHistNClusMatchCent = new TH2F("NClusMatchesCent", "NClusMatchesCent; Cent (%); Nmatches", 100, 0, 100, 11, -0.5, 10.5);
186 
189  fOutput->Add(fHistEbefore);
190  fOutput->Add(fHistEafter);
191  fOutput->Add(fHistEoPCent);
192  fOutput->Add(fHistNMatchCent);
194 
195  if (fIsEmbedded) {
196  for(Int_t icent=0; icent<fNcentBins; ++icent) {
197  name = Form("fHistEmbTrackMatchesOversub_%d",icent);
198  fHistEmbTrackMatchesOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
199  fHistEmbTrackMatchesOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
200  fHistEmbTrackMatchesOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
202 
203  name = Form("fHistNonEmbTrackMatchesOversub_%d",icent);
204  fHistNonEmbTrackMatchesOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
205  fHistNonEmbTrackMatchesOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
206  fHistNonEmbTrackMatchesOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
208 
209  name = Form("fHistOversubMCClusters_%d",icent);
210  fHistOversubMCClusters[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
211  fHistOversubMCClusters[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
212  fHistOversubMCClusters[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
213  fOutput->Add(fHistOversubMCClusters[icent]);
214 
215  name = Form("fHistOversubNonMCClusters_%d",icent);
216  fHistOversubNonMCClusters[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
217  fHistOversubNonMCClusters[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
218  fHistOversubNonMCClusters[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
219  fOutput->Add(fHistOversubNonMCClusters[icent]);
220 
221  name = Form("fHistOversub_%d",icent);
222  fHistOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
223  fHistOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
224  fHistOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
225  fOutput->Add(fHistOversub[icent]);
226  }
227  }
228  fOutput->SetOwner(kTRUE);
229 }
230 
231 //________________________________________________________________________
233 {
234  // Run
235  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
237 
238  // Run the hadronic correction
239 
240  // loop over all clusters
241  fClusCont->ResetCurrentID();
242  AliVCluster *cluster = 0;
243 
244  while ((cluster = fClusCont->GetNextAcceptCluster())) {
245 
246  Double_t energyclus = 0;
247  if (fCreateHisto) {
248  fHistEbefore->Fill(fCent, cluster->GetNonLinCorrEnergy());
249  fHistNclusvsCent->Fill(fCent);
250  }
251 
252  // apply correction / subtraction
253  // to subtract only the closest track set fHadCor to a %
254  // to subtract all tracks within the cut set fHadCor to %+1
255  if (fHadCorr > 1) {
256  energyclus = ApplyHadCorrAllTracks(fClusCont->GetCurrentID(), fHadCorr - 1);
257  }
258  else if (fHadCorr > 0) {
259  energyclus = ApplyHadCorrOneTrack(fClusCont->GetCurrentID(), fHadCorr);
260  }
261  else {
262  energyclus = cluster->GetNonLinCorrEnergy();
263  }
264 
265  if (energyclus < 0) energyclus = 0;
266 
267  cluster->SetHadCorrEnergy(energyclus);
268 
269  if (fCreateHisto) fHistEafter->Fill(fCent, energyclus);
270 
271  }
272 
273  return kTRUE;
274 }
275 
276 //________________________________________________________________________
278 {
279  // Get momenum bin.
280 
281  UInt_t pbin=0;
282  if (p<0.5)
283  pbin=0;
284  else if (p>=0.5 && p<1.0)
285  pbin=1;
286  else if (p>=1.0 && p<1.5)
287  pbin=2;
288  else if (p>=1.5 && p<2.)
289  pbin=3;
290  else if (p>=2. && p<3.)
291  pbin=4;
292  else if (p>=3. && p<4.)
293  pbin=5;
294  else if (p>=4. && p<5.)
295  pbin=6;
296  else if (p>=5. && p<8.)
297  pbin=7;
298  else
299  pbin=8;
300 
301  return pbin;
302 }
303 
304 //________________________________________________________________________
306 {
307  // Get sigma in eta.
308 
309  Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
310  return 2.0*EtaSigma[pbin];
311 }
312 
313 //________________________________________________________________________
315 {
316  // Get phi mean.
317 
318  if (centbin==0) {
319  Double_t PhiMean[9]={0.036,
320  0.021,
321  0.0121,
322  0.0084,
323  0.0060,
324  0.0041,
325  0.0031,
326  0.0022,
327  0.001};
328  return PhiMean[pbin];
329  } else if (centbin==1) {
330  Double_t PhiMean[9]={0.036,
331  0.021,
332  0.0121,
333  0.0084,
334  0.0060,
335  0.0041,
336  0.0031,
337  0.0022,
338  0.001};
339  return PhiMean[pbin];
340  } else if (centbin==2) {
341  Double_t PhiMean[9]={0.036,
342  0.021,
343  0.0121,
344  0.0084,
345  0.0060,
346  0.0041,
347  0.0031,
348  0.0022,
349  0.001};
350  return PhiMean[pbin];
351  } else if (centbin==3) {
352  Double_t PhiMean[9]={0.036,
353  0.021,
354  0.0121,
355  0.0084,
356  0.0060,
357  0.0041,
358  0.0031,
359  0.0022,
360  0.001};
361  return PhiMean[pbin];
362  } else if (centbin==4) {
363  Double_t PhiMean[9]={0.036,
364  0.021,
365  0.0127,
366  0.0089,
367  0.0068,
368  0.0049,
369  0.0038,
370  0.0028,
371  0.0018};
372  return PhiMean[pbin]*(-1.);
373  } else if (centbin==5) {
374  Double_t PhiMean[9]={0.036,
375  0.021,
376  0.0127,
377  0.0089,
378  0.0068,
379  0.0048,
380  0.0038,
381  0.0028,
382  0.0018};
383  return PhiMean[pbin]*(-1.);
384  } else if (centbin==6) {
385  Double_t PhiMean[9]={0.036,
386  0.021,
387  0.0127,
388  0.0089,
389  0.0068,
390  0.0045,
391  0.0035,
392  0.0028,
393  0.0018};
394  return PhiMean[pbin]*(-1.);
395  } else if (centbin==7) {
396  Double_t PhiMean[9]={0.036,
397  0.021,
398  0.0127,
399  0.0089,
400  0.0068,
401  0.0043,
402  0.0035,
403  0.0028,
404  0.0018};
405  return PhiMean[pbin]*(-1.);
406  }
407 
408  return 0;
409 }
410 
411 //________________________________________________________________________
413 {
414  // Get phi sigma.
415 
416  if (centbin==0) {
417  Double_t PhiSigma[9]={0.0221,
418  0.0128,
419  0.0074,
420  0.0064,
421  0.0059,
422  0.0055,
423  0.0052,
424  0.0049,
425  0.0045};
426  return 2.*PhiSigma[pbin];
427  } else if (centbin==1) {
428  Double_t PhiSigma[9]={0.0217,
429  0.0120,
430  0.0076,
431  0.0066,
432  0.0062,
433  0.0058,
434  0.0054,
435  0.0054,
436  0.0045};
437  return 2.*PhiSigma[pbin];
438  } else if (centbin==2) {
439  Double_t PhiSigma[9]={0.0211,
440  0.0124,
441  0.0080,
442  0.0070,
443  0.0067,
444  0.0061,
445  0.0059,
446  0.0054,
447  0.0047};
448  return 2.*PhiSigma[pbin];
449  } else if (centbin==3) {
450  Double_t PhiSigma[9]={0.0215,
451  0.0124,
452  0.0082,
453  0.0073,
454  0.0069,
455  0.0064,
456  0.0060,
457  0.0055,
458  0.0047};
459  return 2.*PhiSigma[pbin];
460  } else if (centbin==4) {
461  Double_t PhiSigma[9]={0.0199,
462  0.0108,
463  0.0072,
464  0.0071,
465  0.0060,
466  0.0055,
467  0.0052,
468  0.0049,
469  0.0045};
470  return 2.*PhiSigma[pbin];
471  } else if (centbin==5) {
472  Double_t PhiSigma[9]={0.0200,
473  0.0110,
474  0.0074,
475  0.0071,
476  0.0064,
477  0.0059,
478  0.0055,
479  0.0052,
480  0.0045};
481  return 2.*PhiSigma[pbin];
482  } else if (centbin==6) {
483  Double_t PhiSigma[9]={0.0202,
484  0.0113,
485  0.0077,
486  0.0071,
487  0.0069,
488  0.0064,
489  0.0060,
490  0.0055,
491  0.0050};
492  return 2.*PhiSigma[pbin];
493  } else if (centbin==7) {
494  Double_t PhiSigma[9]={0.0205,
495  0.0113,
496  0.0080,
497  0.0074,
498  0.0078,
499  0.0067,
500  0.0062,
501  0.0055,
502  0.0050};
503  return 2.*PhiSigma[pbin];
504  }
505 
506  return 0;
507 }
508 
509 //________________________________________________________________________
511  Double_t &totalTrkP, Int_t &Nmatches, Double_t &trkPMCfrac, Int_t &NMCmatches)
512 {
513  // Do the loop over matched tracks for the cluster.
514 
515  AliVCluster* cluster = fClusCont->GetCluster(icluster);
516 
517  if (!cluster) return;
518 
519  // loop over matched tracks
520  Int_t Ntrks = Ntrks = cluster->GetNTracksMatched();
521  for (Int_t i = 0; i < Ntrks; ++i) {
522  AliVTrack* track = 0;
523 
524  if (fEsdMode) {
525  Int_t itrack = cluster->GetTrackMatchedIndex(i);
526  if (itrack >= 0) track = static_cast<AliVTrack*>(fPartCont->GetAcceptParticle(itrack));
527  }
528  else {
529  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(i));
530  UInt_t rejectionReason = 0;
531  if (!fPartCont->AcceptParticle(track, rejectionReason)) track = 0;
532  }
533 
534  if (!track) continue;
535 
536  Double_t etadiff = 999;
537  Double_t phidiff = 999;
538  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
539  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(etadiff, phidiff);
540 
541  // check if track also points to cluster
542  if (fDoTrackClus && (track->GetEMCALcluster() != icluster)) continue;
543 
544  Double_t mom = track->P();
545  UInt_t mombin = GetMomBin(mom);
546  Int_t centbinch = fCentBin;
547 
548  if (track->Charge() < 0) centbinch += fNcentBins;
549 
550  if (fCreateHisto) {
551  Int_t etabin = 0;
552  if (track->Eta() > 0) etabin=1;
553  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
554  fHistMatchEtaPhiAll->Fill(etadiff, phidiff);
555  }
556 
557  Double_t etaCut = 0.0;
558  Double_t phiCutlo = 0.0;
559  Double_t phiCuthi = 0.0;
560 
561  if (fPhiMatch > 0) {
562  phiCutlo = -fPhiMatch;
563  phiCuthi = +fPhiMatch;
564  }
565  else {
566  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
567  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
568  }
569 
570  if (fEtaMatch > 0) {
571  etaCut = fEtaMatch;
572  }
573  else {
574  etaCut = GetEtaSigma(mombin);
575  }
576 
577  if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
578  if (track->GetLabel() > fMinMCLabel) {
579  ++NMCmatches;
580  trkPMCfrac += mom;
581  }
582  ++Nmatches;
583  totalTrkP += mom;
584 
585  if (fCreateHisto) {
586  if (fHadCorr > 1) {
587  Double_t dphi = 0;
588  Double_t deta = 0;
589  GetEtaPhiDiff(track, cluster, dphi, deta);
590  Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
591  Double_t energyclus = cluster->GetNonLinCorrEnergy();
592  fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
593  }
594  }
595  }
596  }
597 
598  if (totalTrkP > 0) trkPMCfrac /= totalTrkP;
599 }
600 
601 //________________________________________________________________________
603 {
604  // Apply the hadronic correction with one track only.
605 
606  AliVCluster* cluster = fClusCont->GetCluster(icluster);
607 
608  Double_t energyclus = cluster->GetNonLinCorrEnergy();
609 
610  AliVTrack* track = 0;
611 
612  if (cluster->GetNTracksMatched() > 0) {
613  if (fEsdMode) {
614  Int_t itrack = cluster->GetTrackMatchedIndex(0);
615  if (itrack >= 0) track = static_cast<AliVTrack*>(fPartCont->GetAcceptParticle(itrack));
616  }
617  else {
618  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
619  UInt_t rejectionReason = 0;
620  if (!fPartCont->AcceptParticle(track, rejectionReason)) track = 0;
621  }
622  }
623 
624  if (!track || track->P() < 1e-6) return energyclus;
625 
626  Double_t mom = track->P();
627 
628  Double_t dEtaMin = 1e9;
629  Double_t dPhiMin = 1e9;
630  GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
631 
632  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(dEtaMin, dPhiMin);
633 
634  // check if track also points to cluster
635  Int_t cid = track->GetEMCALcluster();
636  if (fDoTrackClus && (cid != icluster)) return energyclus;
637 
638  UInt_t mombin = GetMomBin(mom);
639  Int_t centbinch = fCentBin;
640  if (track->Charge() < 0) centbinch += fNcentBins;
641 
642  // plot some histograms if switched on
643  if (fCreateHisto) {
644  Int_t etabin = 0;
645  if(track->Eta() > 0) etabin = 1;
646 
647  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
648  fHistMatchEtaPhiAll->Fill(dEtaMin, dPhiMin);
649 
650  if (mom > 0) {
651  Double_t etadiff = 0;
652  Double_t phidiff = 0;
653  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
654  Double_t dRmin = TMath::Sqrt(etadiff*etadiff + phidiff*phidiff);
655  fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
656  fHistEoPCent->Fill(fCent, energyclus / mom);
657  fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
658  }
659  }
660 
661  // define eta/phi cuts
662  Double_t etaCut = 0.0;
663  Double_t phiCutlo = 0.0;
664  Double_t phiCuthi = 0.0;
665  if (fPhiMatch > 0) {
666  phiCutlo = -fPhiMatch;
667  phiCuthi = +fPhiMatch;
668  }
669  else {
670  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
671  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
672  }
673  if (fEtaMatch > 0) {
674  etaCut = fEtaMatch;
675  }
676  else {
677  etaCut = GetEtaSigma(mombin);
678  }
679 
680  // apply the correction if the track is in the eta/phi window
681  if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
682  energyclus -= hadCorr * mom;
683  }
684 
685  return energyclus;
686 }
687 
688 //________________________________________________________________________
690 {
691  // Apply the hadronic correction with all tracks.
692 
693  AliVCluster* cluster = fClusCont->GetCluster(icluster);
694 
695  Double_t energyclus = cluster->GetNonLinCorrEnergy();
696  Double_t cNcells = cluster->GetNCells();
697 
698  Double_t totalTrkP = 0.0; // count total track momentum
699  Int_t Nmatches = 0; // count total number of matches
700 
701  Double_t trkPMCfrac = 0.0; // count total track momentum
702  Int_t NMCmatches = 0; // count total number of matches
703 
704  // do the loop over the matched tracks and get the number of matches and the total momentum
705  DoMatchedTracksLoop(icluster, totalTrkP, Nmatches, trkPMCfrac, NMCmatches);
706 
707  Double_t Esub = hadCorr * totalTrkP;
708 
709  if (Esub > energyclus) Esub = energyclus;
710 
711  // applying Peter's proposed algorithm
712  // never subtract the full energy of the cluster
713  Double_t clusEexcl = fEexclCell * cNcells;
714  if (energyclus < clusEexcl) clusEexcl = energyclus;
715  if ((energyclus - Esub) < clusEexcl) Esub = (energyclus - clusEexcl);
716 
717  // embedding
718  Double_t EsubMC = 0;
719  Double_t EsubBkg = 0;
720  Double_t EclusMC = 0;
721  Double_t EclusBkg = 0;
722  Double_t EclusCorr = 0;
723  Double_t EclusMCcorr = 0;
724  Double_t EclusBkgcorr = 0;
725  Double_t overSub = 0;
726  if (fIsEmbedded) {
727  EsubMC = hadCorr * totalTrkP * trkPMCfrac;
728  EsubBkg = hadCorr * totalTrkP - EsubMC;
729  EclusMC = energyclus * cluster->GetMCEnergyFraction();
730  EclusBkg = energyclus - EclusMC;
731 
732  if (energyclus > Esub)
733  EclusCorr = energyclus - Esub;
734 
735  if (EclusMC > EsubMC)
736  EclusMCcorr = EclusMC - EsubMC;
737 
738  if (EclusBkg > EsubBkg)
739  EclusBkgcorr = EclusBkg - EsubBkg;
740 
741  overSub = EclusMCcorr + EclusBkgcorr - EclusCorr;
742  }
743 
744  // plot some histograms if switched on
745  if (fCreateHisto) {
746  fHistNMatchCent->Fill(fCent, Nmatches);
747  fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
748 
749  if (Nmatches > 0) fHistNclusMatchvsCent->Fill(fCent);
750 
751  if (Nmatches < 3) {
752  fHistNCellsEnergy[fCentBin][Nmatches]->Fill(energyclus, cNcells);
753  }
754  else {
755  fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
756  }
757 
758  if (totalTrkP > 0) {
759  Double_t EoP = energyclus / totalTrkP;
760  fHistEoPCent->Fill(fCent, EoP);
761  fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
762 
763  fHistEsubPchRatAll[fCentBin]->Fill(totalTrkP, Esub / totalTrkP);
764 
765  if (Nmatches == 1) {
766  AliVTrack* track = 0;
767  if (fEsdMode) {
768  Int_t itrack = cluster->GetTrackMatchedIndex(0);
769  if (itrack >= 0) track = static_cast<AliVTrack*>(fPartCont->GetAcceptParticle(itrack));
770  }
771  else {
772  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
773  UInt_t rejectionReason = 0;
774  if (!fPartCont->AcceptParticle(track, rejectionReason)) track = 0;
775  }
776  if (track) {
777  Int_t centbinchm = fCentBin;
778  if (track->Charge() < 0) centbinchm += fNcentBins;
779  fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
780  fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
781  }
782  }
783 
784  if (fIsEmbedded) {
785  fHistOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
786 
787  if (cluster->GetMCEnergyFraction() > 0.95)
788  fHistOversubMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
789  else if (cluster->GetMCEnergyFraction() < 0.05)
790  fHistOversubNonMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
791 
792  if (trkPMCfrac < 0.05)
793  fHistNonEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
794  else if (trkPMCfrac > 0.95)
795  fHistEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
796  }
797  }
798  }
799 
800  if (fIsEmbedded && fDoExact) {
801  Esub -= overSub;
802  if (EclusBkgcorr + EclusMCcorr > 0) {
803  Double_t newfrac = EclusMCcorr / (EclusBkgcorr + EclusMCcorr);
804  cluster->SetMCEnergyFraction(newfrac);
805  }
806  }
807 
808  // apply the correction
809  energyclus -= Esub;
810 
811  return energyclus;
812 }
Int_t fNcentBins
How many centrality bins (this member copied from AliAnalysisTaskEmcal)
TH2 * fHistEsubPchRatAll[8]
!Esub/momentum of matched tracks vs. total momentum of matched tracks (all number of matches) ...
double Double_t
Definition: External.C:58
Definition: External.C:236
TH2 * fHistMatchEtaPhiAll
!deta vs. dphi of matched cluster-track pairs
static RegisterCorrectionComponent< AliEmcalCorrectionClusterHadronicCorrection > reg
TH2 * fHistOversubNonMCClusters[4]
!Over-subtracted energy / cluster energy (cluster MC energy fraction < 5%)
TH2 * fHistOversub[4]
!Over-subtracted energy / cluster energy
Int_t fCentBin
! Event centrality bin
TH2 * fHistMatchEvsP[4]
!cluster energy vs. track momentum of matched pairs
TH2 * fHistEsubPchRat[8]
!Esub/momentum of matched tracks vs. total momentum of matched tracks (only 1 match) ...
TH1 * fHistEafter
!average energy of clusters after correction vs. centrality
AliClusterContainer * fClusCont
Pointer to the cluster container.
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.
Double_t fEexclCell
energy/cell that we cannot subtract from the clusters
Int_t fMinMCLabel
Minimum MC label value for the tracks/clusters being considered MC particles.
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
virtual Bool_t AcceptParticle(const AliVParticle *vp, UInt_t &rejectionReason) const
TH2 * fHistNonEmbTrackMatchesOversub[4]
!Over-subtracted energy / cluster energy with non-embedded track matches (embedded matches < 5%) ...
AliVCluster * GetCluster(Int_t i) const
Double_t fMaxBinPt
Max pt in histograms.
TList * fOutput
! List of output histograms
Bool_t fCreateHisto
Default YAML configuration.
Hadronic correction component in the EMCal correction framework.
TH1 * fHistEbefore
!average energy of clusters before correction vs. centrality
Double_t fCent
! Event centrality
TH2 * fHistMatchEtaPhi[8][9][2]
!deta vs. dphi of matched cluster-track pairs
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.
TH2 * fHistEmbTrackMatchesOversub[4]
!Over-subtracted energy / cluster energy with embedded track matches (non-embedded matches < 5%) ...
TH1 * fHistNclusMatchvsCent
!n clusters matched to some track vs. centrality
void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
TH2 * fHistMatchEtaPhiAllTr
!deta vs. dphi of all cluster-track pairs (tr loop)
bool Bool_t
Definition: External.C:53
TH1 * fHistEsubPch[8]
!Esub vs. total momentum of matched tracks (only 1 match)
TH2 * fHistOversubMCClusters[4]
!Over-subtracted energy / cluster energy (cluster MC energy fraction > 95%)
AliParticleContainer * fPartCont
Pointer to the track/particle container.
AliVCluster * GetNextAcceptCluster()
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.