AliPhysics  cdeda5a (cdeda5a)
 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  fClusterContainerIndexMap(),
30  fParticleContainerIndexMap(),
31  fHistMatchEtaPhiAll(0),
32  fHistMatchEtaPhiAllTr(0),
33  fHistMatchEtaPhiAllCl(0),
34  fHistNclusvsCent(0),
35  fHistNclusMatchvsCent(0),
36  fHistEbefore(0),
37  fHistEafter(0),
38  fHistEoPCent(0),
39  fHistNMatchCent(0),
40  fHistNClusMatchCent(0)
41 {
42  // Default constructor
43  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
44 
45  for(Int_t i=0; i<8; i++) {
46  fHistEsubPch[i] = 0;
47  fHistEsubPchRat[i] = 0;
48  fHistEsubPchRatAll[i] = 0;
49 
50  if (i<4) {
51  fHistMatchEvsP[i] = 0;
52  fHistMatchdRvsEP[i] = 0;
53  fHistNMatchEnergy[i] = 0;
54 
59  fHistOversub[i] = 0;
60 
61  for(Int_t j=0; j<4; j++)
62  fHistNCellsEnergy[i][j] = 0;
63  }
64 
65  for(Int_t j=0; j<9; j++) {
66  for(Int_t k=0; k<2; k++)
67  fHistMatchEtaPhi[i][j][k] = 0;
68  }
69  }
70 }
71 
72 //________________________________________________________________________
74 {
75  // Destructor
76 }
77 
78 //________________________________________________________________________
80 {
81  // Initialization
82  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
84  // Do base class initializations and if it fails -> bail out
85  //AliAnalysisTaskEmcal::ExecOnce();
86  //if (!fInitialized) return;
87 
88  GetProperty("createHistos", fCreateHisto);
89  GetProperty("phiMatch", fPhiMatch);
90  GetProperty("etaMatch", fEtaMatch);
91  GetProperty("hadCorr", fHadCorr);
92  GetProperty("Eexcl", fEexclCell);
93  GetProperty("doTrackClus", fDoTrackClus);
94 
95  if (!fEsdMode && fParticleCollArray.GetEntries() > 1) {
96  AliWarning("================================================================================");
97  AliWarning("== Added multiple particle containers when running with AOD!");
98  AliWarning("== Particle selection of the first particle container will be applied");
99  AliWarning("== to _ALL_ particles! If you need a different selection, then change");
100  AliWarning("== the order of adding the containers so that the desired container is first!");
101  AliWarning("================================================================================");
102  }
103 
104  return kTRUE;
105 }
106 
107 //________________________________________________________________________
109 {
110  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
112 
113  // Create my user objects.
114 
115  if (!fCreateHisto) return;
116 
117  fHistMatchEtaPhiAll = new TH2F("fHistMatchEtaPhiAll", "fHistMatchEtaPhiAll;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
119 
120  fHistMatchEtaPhiAllTr = new TH2F("fHistMatchEtaPhiAllTr", "fHistMatchEtaPhiAllTr;#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  name = Form("fHistEsubPchRatAll_%i",icent);
154  temp = Form("%s (all Nmatches)",name.Data());
155  fHistEsubPchRatAll[icent]=new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
156  fHistEsubPchRatAll[icent]->SetXTitle("#Sigma p (GeV)");
157  fHistEsubPchRatAll[icent]->SetYTitle("E_{sub} / #sum p");
158  fOutput->Add(fHistEsubPchRatAll[icent]);
159 
160  if (icent<fNcentBins) {
161  for(Int_t itrk=0; itrk<4; ++itrk) {
162  name = Form("fHistNCellsEnergy_%i_%i",icent,itrk);
163  temp = Form("%s (Nmatches==%d);N_{cells};E_{clus} (GeV)",name.Data(),itrk);
164  fHistNCellsEnergy[icent][itrk] = new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, 101, -0.5, 100.5);
165  fOutput->Add(fHistNCellsEnergy[icent][itrk]);
166  }
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 
206  if (fIsEmbedded) {
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 
242 //________________________________________________________________________
244 {
247 }
248 
249 //________________________________________________________________________
251 {
252  // Run
253  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
255 
256  // Run the hadronic correction
257 
258  // loop over all clusters
259  AliVCluster *cluster = 0;
260  AliClusterContainer * clusCont = 0;
261  TIter nextClusCont(&fClusterCollArray);
262  while ((clusCont = static_cast<AliClusterContainer*>(nextClusCont()))) {
263  auto clusItCont = clusCont->accepted_momentum();
264  for (AliClusterIterableMomentumContainer::iterator clusIterator = clusItCont.begin(); clusIterator != clusItCont.end(); ++clusIterator) {
265  cluster = static_cast<AliVCluster *>(clusIterator->second);
266 
267  Double_t energyclus = 0;
268  if (fCreateHisto) {
269  fHistEbefore->Fill(fCent, cluster->GetNonLinCorrEnergy());
270  fHistNclusvsCent->Fill(fCent);
271  }
272 
273  // apply correction / subtraction
274  // to subtract only the closest track set fHadCor to a %
275  // to subtract all tracks within the cut set fHadCor to %+1
276  if (fHadCorr > 1) {
277  energyclus = ApplyHadCorrAllTracks(fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, clusIterator.current_index()), fHadCorr - 1);
278  }
279  else if (fHadCorr > 0) {
280  energyclus = ApplyHadCorrOneTrack(fClusterContainerIndexMap.GlobalIndexFromLocalIndex(clusCont, clusIterator.current_index()), fHadCorr);
281  }
282  else {
283  energyclus = cluster->GetNonLinCorrEnergy();
284  }
285 
286  if (energyclus < 0) energyclus = 0;
287 
288  cluster->SetHadCorrEnergy(energyclus);
289 
290  if (fCreateHisto) fHistEafter->Fill(fCent, energyclus);
291 
292  }
293  }
294 
295  return kTRUE;
296 }
297 
298 //________________________________________________________________________
300 {
301  // Get momenum bin.
302 
303  UInt_t pbin=0;
304  if (p<0.5)
305  pbin=0;
306  else if (p>=0.5 && p<1.0)
307  pbin=1;
308  else if (p>=1.0 && p<1.5)
309  pbin=2;
310  else if (p>=1.5 && p<2.)
311  pbin=3;
312  else if (p>=2. && p<3.)
313  pbin=4;
314  else if (p>=3. && p<4.)
315  pbin=5;
316  else if (p>=4. && p<5.)
317  pbin=6;
318  else if (p>=5. && p<8.)
319  pbin=7;
320  else
321  pbin=8;
322 
323  return pbin;
324 }
325 
326 //________________________________________________________________________
328 {
329  // Get sigma in eta.
330 
331  Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
332  return 2.0*EtaSigma[pbin];
333 }
334 
335 //________________________________________________________________________
337 {
338  // Get phi mean.
339 
340  if (centbin==0) {
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==1) {
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==2) {
363  Double_t PhiMean[9]={0.036,
364  0.021,
365  0.0121,
366  0.0084,
367  0.0060,
368  0.0041,
369  0.0031,
370  0.0022,
371  0.001};
372  return PhiMean[pbin];
373  } else if (centbin==3) {
374  Double_t PhiMean[9]={0.036,
375  0.021,
376  0.0121,
377  0.0084,
378  0.0060,
379  0.0041,
380  0.0031,
381  0.0022,
382  0.001};
383  return PhiMean[pbin];
384  } else if (centbin==4) {
385  Double_t PhiMean[9]={0.036,
386  0.021,
387  0.0127,
388  0.0089,
389  0.0068,
390  0.0049,
391  0.0038,
392  0.0028,
393  0.0018};
394  return PhiMean[pbin]*(-1.);
395  } else if (centbin==5) {
396  Double_t PhiMean[9]={0.036,
397  0.021,
398  0.0127,
399  0.0089,
400  0.0068,
401  0.0048,
402  0.0038,
403  0.0028,
404  0.0018};
405  return PhiMean[pbin]*(-1.);
406  } else if (centbin==6) {
407  Double_t PhiMean[9]={0.036,
408  0.021,
409  0.0127,
410  0.0089,
411  0.0068,
412  0.0045,
413  0.0035,
414  0.0028,
415  0.0018};
416  return PhiMean[pbin]*(-1.);
417  } else if (centbin==7) {
418  Double_t PhiMean[9]={0.036,
419  0.021,
420  0.0127,
421  0.0089,
422  0.0068,
423  0.0043,
424  0.0035,
425  0.0028,
426  0.0018};
427  return PhiMean[pbin]*(-1.);
428  }
429 
430  return 0;
431 }
432 
433 //________________________________________________________________________
435 {
436  // Get phi sigma.
437 
438  if (centbin==0) {
439  Double_t PhiSigma[9]={0.0221,
440  0.0128,
441  0.0074,
442  0.0064,
443  0.0059,
444  0.0055,
445  0.0052,
446  0.0049,
447  0.0045};
448  return 2.*PhiSigma[pbin];
449  } else if (centbin==1) {
450  Double_t PhiSigma[9]={0.0217,
451  0.0120,
452  0.0076,
453  0.0066,
454  0.0062,
455  0.0058,
456  0.0054,
457  0.0054,
458  0.0045};
459  return 2.*PhiSigma[pbin];
460  } else if (centbin==2) {
461  Double_t PhiSigma[9]={0.0211,
462  0.0124,
463  0.0080,
464  0.0070,
465  0.0067,
466  0.0061,
467  0.0059,
468  0.0054,
469  0.0047};
470  return 2.*PhiSigma[pbin];
471  } else if (centbin==3) {
472  Double_t PhiSigma[9]={0.0215,
473  0.0124,
474  0.0082,
475  0.0073,
476  0.0069,
477  0.0064,
478  0.0060,
479  0.0055,
480  0.0047};
481  return 2.*PhiSigma[pbin];
482  } else if (centbin==4) {
483  Double_t PhiSigma[9]={0.0199,
484  0.0108,
485  0.0072,
486  0.0071,
487  0.0060,
488  0.0055,
489  0.0052,
490  0.0049,
491  0.0045};
492  return 2.*PhiSigma[pbin];
493  } else if (centbin==5) {
494  Double_t PhiSigma[9]={0.0200,
495  0.0110,
496  0.0074,
497  0.0071,
498  0.0064,
499  0.0059,
500  0.0055,
501  0.0052,
502  0.0045};
503  return 2.*PhiSigma[pbin];
504  } else if (centbin==6) {
505  Double_t PhiSigma[9]={0.0202,
506  0.0113,
507  0.0077,
508  0.0071,
509  0.0069,
510  0.0064,
511  0.0060,
512  0.0055,
513  0.0050};
514  return 2.*PhiSigma[pbin];
515  } else if (centbin==7) {
516  Double_t PhiSigma[9]={0.0205,
517  0.0113,
518  0.0080,
519  0.0074,
520  0.0078,
521  0.0067,
522  0.0062,
523  0.0055,
524  0.0050};
525  return 2.*PhiSigma[pbin];
526  }
527 
528  return 0;
529 }
530 
531 //________________________________________________________________________
533  Double_t &totalTrkP, Int_t &Nmatches, Double_t &trkPMCfrac, Int_t &NMCmatches)
534 {
535  // Do the loop over matched tracks for the cluster.
536 
537  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
538 
539  if (!cluster) return;
540 
541  // loop over matched tracks
542  Int_t Ntrks = Ntrks = cluster->GetNTracksMatched();
543  for (Int_t i = 0; i < Ntrks; ++i) {
544  AliVTrack* track = 0;
545 
546  if (fEsdMode) {
547  Int_t itrack = cluster->GetTrackMatchedIndex(i);
548  if (itrack >= 0) {
550  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
551  }
552  }
553  else {
554  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(i));
555  UInt_t rejectionReason = 0;
557  if (!partCont) { AliError("No particle container available!"); }
558  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
559  }
560 
561  if (!track) continue;
562 
563  Double_t etadiff = 999;
564  Double_t phidiff = 999;
565  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
566  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(etadiff, phidiff);
567 
568  // check if track also points to cluster
569  if (fDoTrackClus && (track->GetEMCALcluster() != icluster)) continue;
570 
571  Double_t mom = track->P();
572  UInt_t mombin = GetMomBin(mom);
573  Int_t centbinch = fCentBin;
574 
575  if (track->Charge() < 0) centbinch += fNcentBins;
576 
577  if (fCreateHisto) {
578  Int_t etabin = 0;
579  if (track->Eta() > 0) etabin=1;
580  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
581  fHistMatchEtaPhiAll->Fill(etadiff, phidiff);
582  }
583 
584  Double_t etaCut = 0.0;
585  Double_t phiCutlo = 0.0;
586  Double_t phiCuthi = 0.0;
587 
588  if (fPhiMatch > 0) {
589  phiCutlo = -fPhiMatch;
590  phiCuthi = +fPhiMatch;
591  }
592  else {
593  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
594  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
595  }
596 
597  if (fEtaMatch > 0) {
598  etaCut = fEtaMatch;
599  }
600  else {
601  etaCut = GetEtaSigma(mombin);
602  }
603 
604  if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
605  if (track->GetLabel() > fMinMCLabel) {
606  ++NMCmatches;
607  trkPMCfrac += mom;
608  }
609  ++Nmatches;
610  totalTrkP += mom;
611 
612  if (fCreateHisto) {
613  if (fHadCorr > 1) {
614  Double_t dphi = 0;
615  Double_t deta = 0;
616  GetEtaPhiDiff(track, cluster, dphi, deta);
617  Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
618  Double_t energyclus = cluster->GetNonLinCorrEnergy();
619  fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
620  }
621  }
622  }
623  }
624 
625  if (totalTrkP > 0) trkPMCfrac /= totalTrkP;
626 }
627 
628 //________________________________________________________________________
630 {
631  // Apply the hadronic correction with one track only.
632 
633  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
634 
635  Double_t energyclus = cluster->GetNonLinCorrEnergy();
636 
637  AliVTrack* track = 0;
638 
639  if (cluster->GetNTracksMatched() > 0) {
640  if (fEsdMode) {
641  Int_t itrack = cluster->GetTrackMatchedIndex(0);
642  if (itrack >= 0) {
644  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
645  }
646  }
647  else {
648  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
649  UInt_t rejectionReason = 0;
651  if (!partCont) { AliError("No particle container available!"); }
652  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
653  }
654  }
655 
656  if (!track || track->P() < 1e-6) return energyclus;
657 
658  Double_t mom = track->P();
659 
660  Double_t dEtaMin = 1e9;
661  Double_t dPhiMin = 1e9;
662  GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
663 
664  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(dEtaMin, dPhiMin);
665 
666  // check if track also points to cluster
667  Int_t cid = track->GetEMCALcluster();
668  if (fDoTrackClus && (cid != icluster)) return energyclus;
669 
670  UInt_t mombin = GetMomBin(mom);
671  Int_t centbinch = fCentBin;
672  if (track->Charge() < 0) centbinch += fNcentBins;
673 
674  // plot some histograms if switched on
675  if (fCreateHisto) {
676  Int_t etabin = 0;
677  if(track->Eta() > 0) etabin = 1;
678 
679  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
680  fHistMatchEtaPhiAll->Fill(dEtaMin, dPhiMin);
681 
682  if (mom > 0) {
683  Double_t etadiff = 0;
684  Double_t phidiff = 0;
685  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
686  Double_t dRmin = TMath::Sqrt(etadiff*etadiff + phidiff*phidiff);
687  fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
688  fHistEoPCent->Fill(fCent, energyclus / mom);
689  fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
690  }
691  }
692 
693  // define eta/phi cuts
694  Double_t etaCut = 0.0;
695  Double_t phiCutlo = 0.0;
696  Double_t phiCuthi = 0.0;
697  if (fPhiMatch > 0) {
698  phiCutlo = -fPhiMatch;
699  phiCuthi = +fPhiMatch;
700  }
701  else {
702  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
703  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
704  }
705  if (fEtaMatch > 0) {
706  etaCut = fEtaMatch;
707  }
708  else {
709  etaCut = GetEtaSigma(mombin);
710  }
711 
712  // apply the correction if the track is in the eta/phi window
713  if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
714  energyclus -= hadCorr * mom;
715  }
716 
717  return energyclus;
718 }
719 
720 //________________________________________________________________________
722 {
723  // Apply the hadronic correction with all tracks.
724 
725  AliVCluster* cluster = fClusterContainerIndexMap.GetObjectFromGlobalIndex(icluster);
726 
727  Double_t energyclus = cluster->GetNonLinCorrEnergy();
728  Double_t cNcells = cluster->GetNCells();
729 
730  Double_t totalTrkP = 0.0; // count total track momentum
731  Int_t Nmatches = 0; // count total number of matches
732 
733  Double_t trkPMCfrac = 0.0; // count total track momentum
734  Int_t NMCmatches = 0; // count total number of matches
735 
736  // do the loop over the matched tracks and get the number of matches and the total momentum
737  DoMatchedTracksLoop(icluster, totalTrkP, Nmatches, trkPMCfrac, NMCmatches);
738 
739  Double_t Esub = hadCorr * totalTrkP;
740 
741  if (Esub > energyclus) Esub = energyclus;
742 
743  // applying Peter's proposed algorithm
744  // never subtract the full energy of the cluster
745  Double_t clusEexcl = fEexclCell * cNcells;
746  if (energyclus < clusEexcl) clusEexcl = energyclus;
747  if ((energyclus - Esub) < clusEexcl) Esub = (energyclus - clusEexcl);
748 
749  // embedding
750  Double_t EsubMC = 0;
751  Double_t EsubBkg = 0;
752  Double_t EclusMC = 0;
753  Double_t EclusBkg = 0;
754  Double_t EclusCorr = 0;
755  Double_t EclusMCcorr = 0;
756  Double_t EclusBkgcorr = 0;
757  Double_t overSub = 0;
758  if (fIsEmbedded) {
759  EsubMC = hadCorr * totalTrkP * trkPMCfrac;
760  EsubBkg = hadCorr * totalTrkP - EsubMC;
761  EclusMC = energyclus * cluster->GetMCEnergyFraction();
762  EclusBkg = energyclus - EclusMC;
763 
764  if (energyclus > Esub)
765  EclusCorr = energyclus - Esub;
766 
767  if (EclusMC > EsubMC)
768  EclusMCcorr = EclusMC - EsubMC;
769 
770  if (EclusBkg > EsubBkg)
771  EclusBkgcorr = EclusBkg - EsubBkg;
772 
773  overSub = EclusMCcorr + EclusBkgcorr - EclusCorr;
774  }
775 
776  // plot some histograms if switched on
777  if (fCreateHisto) {
778  fHistNMatchCent->Fill(fCent, Nmatches);
779  fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
780 
781  if (Nmatches > 0) fHistNclusMatchvsCent->Fill(fCent);
782 
783  if (Nmatches < 3) {
784  fHistNCellsEnergy[fCentBin][Nmatches]->Fill(energyclus, cNcells);
785  }
786  else {
787  fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
788  }
789 
790  if (totalTrkP > 0) {
791  Double_t EoP = energyclus / totalTrkP;
792  fHistEoPCent->Fill(fCent, EoP);
793  fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
794 
795  fHistEsubPchRatAll[fCentBin]->Fill(totalTrkP, Esub / totalTrkP);
796 
797  if (Nmatches == 1) {
798  AliVTrack* track = 0;
799  if (fEsdMode) {
800  Int_t itrack = cluster->GetTrackMatchedIndex(0);
801  if (itrack >= 0) {
803  track = static_cast<AliVTrack*>(res.second->GetAcceptParticle(res.first));
804  }
805  }
806  else {
807  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
808  UInt_t rejectionReason = 0;
810  if (!partCont) { AliError("No particle container available!"); }
811  if (!partCont->AcceptParticle(track, rejectionReason)) track = 0;
812  }
813  if (track) {
814  Int_t centbinchm = fCentBin;
815  if (track->Charge() < 0) centbinchm += fNcentBins;
816  fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
817  fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
818  }
819  }
820 
821  if (fIsEmbedded) {
822  fHistOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
823 
824  if (cluster->GetMCEnergyFraction() > 0.95)
825  fHistOversubMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
826  else if (cluster->GetMCEnergyFraction() < 0.05)
827  fHistOversubNonMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
828 
829  if (trkPMCfrac < 0.05)
830  fHistNonEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
831  else if (trkPMCfrac > 0.95)
832  fHistEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
833  }
834  }
835  }
836 
837  if (fIsEmbedded && fDoExact) {
838  Esub -= overSub;
839  if (EclusBkgcorr + EclusMCcorr > 0) {
840  Double_t newfrac = EclusMCcorr / (EclusBkgcorr + EclusMCcorr);
841  cluster->SetMCEnergyFraction(newfrac);
842  }
843  }
844 
845  // apply the correction
846  energyclus -= Esub;
847 
848  return energyclus;
849 }
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) ...
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
const AliClusterIterableMomentumContainer accepted_momentum() const
bidirectional stl iterator over the EMCAL iterable container
std::pair< int, U * > LocalIndexFromGlobalIndex(const int globalIndex) const
static RegisterCorrectionComponent< AliEmcalCorrectionClusterHadronicCorrection > reg
AliEmcalContainerIndexMap< AliClusterContainer, AliVCluster > fClusterContainerIndexMap
! Mapping between index and cluster containers
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
int GlobalIndexFromLocalIndex(const U *inputObject, const int localIndex) const
TH2 * fHistMatchEvsP[4]
!cluster energy vs. track momentum of matched pairs
Container for particles within the EMCAL framework.
TObjArray fParticleCollArray
Particle/track collection array.
TH2 * fHistEsubPchRat[8]
!Esub/momentum of matched tracks vs. total momentum of matched tracks (only 1 match) ...
TObjArray fClusterCollArray
Cluster collection array.
TH1 * fHistEafter
!average energy of clusters after correction vs. centrality
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
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
TH2 * fHistNonEmbTrackMatchesOversub[4]
!Over-subtracted energy / cluster energy with non-embedded track matches (embedded matches < 5%) ...
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
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%) ...
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)
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)
AliEmcalContainerIndexMap< AliParticleContainer, AliVParticle > fParticleContainerIndexMap
! Mapping between index and particle containers
TH2 * fHistOversubMCClusters[4]
!Over-subtracted energy / cluster energy (cluster MC energy fraction > 95%)
Container structure for EMCAL clusters.
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.