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