AliPhysics  f05a842 (f05a842)
 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 // Author: R.Reed, C.Loizides, S.Aiola
4 
6 
10 
11 // Actually registers the class with the base class
13 
14 //________________________________________________________________________
16  AliEmcalCorrectionComponent("AliEmcalCorrectionClusterHadronicCorrection"),
17  fPhiMatch(0.05),
18  fEtaMatch(0.025),
19  fDoTrackClus(kTRUE),
20  fHadCorr(0),
21  fEexclCell(0),
22  fDoExact(kFALSE),
23  fHistMatchEtaPhiAll(0),
24  fHistMatchEtaPhiAllTr(0),
25  fHistMatchEtaPhiAllCl(0),
26  fHistNclusvsCent(0),
27  fHistNclusMatchvsCent(0),
28  fHistEbefore(0),
29  fHistEafter(0),
30  fHistEoPCent(0),
31  fHistNMatchCent(0),
32  fHistNClusMatchCent(0)
33 {
34  // Default constructor
35  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
36 
37  for(Int_t i=0; i<8; i++) {
38  fHistEsubPch[i] = 0;
39  fHistEsubPchRat[i] = 0;
40  fHistEsubPchRatAll[i] = 0;
41 
42  if (i<4) {
43  fHistMatchEvsP[i] = 0;
44  fHistMatchdRvsEP[i] = 0;
45  fHistNMatchEnergy[i] = 0;
46 
51  fHistOversub[i] = 0;
52 
53  for(Int_t j=0; j<4; j++)
54  fHistNCellsEnergy[i][j] = 0;
55  }
56 
57  for(Int_t j=0; j<9; j++) {
58  for(Int_t k=0; k<2; k++)
59  fHistMatchEtaPhi[i][j][k] = 0;
60  }
61  }
62 }
63 
64 //________________________________________________________________________
66 {
67  // Destructor
68 }
69 
70 //________________________________________________________________________
72 {
73  // Initialization
74  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
76  // Do base class initializations and if it fails -> bail out
77  //AliAnalysisTaskEmcal::ExecOnce();
78  //if (!fInitialized) return;
79 
80  GetProperty("createHistos", fCreateHisto);
81  GetProperty("phiMatch", fPhiMatch);
82  GetProperty("etaMatch", fEtaMatch);
83  GetProperty("hadCorr", fHadCorr);
84  GetProperty("Eexcl", fEexclCell);
85  GetProperty("doTrackClus", fDoTrackClus);
86  Double_t trackPtMin = 0.15;
87  GetProperty("trackPtMin", trackPtMin);
88  Double_t clusterNonLinCorrEnergyMin = 0.15;
89  GetProperty("clusterNonLinCorrEnergyMin", clusterNonLinCorrEnergyMin);
90  Double_t clusterECut = 0.0;
91  GetProperty("clusterEMin", clusterECut);
92  Double_t clusterPtCut = 0.0;
93  GetProperty("clusterPtMin", clusterPtCut);
94 
96  fPartCont->SetParticlePtCut(trackPtMin);
98  fClusCont->SetClusNonLinCorrEnergyCut(clusterNonLinCorrEnergyMin);
99  fClusCont->SetClusECut(clusterECut);
100  fClusCont->SetClusPtCut(clusterPtCut);
101 
102  // Create my user objects.
103 
104  if (!fCreateHisto) return kTRUE;
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  return kTRUE;
231 }
232 
233 //________________________________________________________________________
235 {
236  // Run
237  AliDebug(3, Form("%s", __PRETTY_FUNCTION__));
239 
240  // Run the hadronic correction
241 
242  // loop over all clusters
243  fClusCont->ResetCurrentID();
244  AliVCluster *cluster = 0;
245 
246  while ((cluster = fClusCont->GetNextAcceptCluster())) {
247 
248  Double_t energyclus = 0;
249  if (fCreateHisto) {
250  fHistEbefore->Fill(fCent, cluster->GetNonLinCorrEnergy());
251  fHistNclusvsCent->Fill(fCent);
252  }
253 
254  // apply correction / subtraction
255  // to subtract only the closest track set fHadCor to a %
256  // to subtract all tracks within the cut set fHadCor to %+1
257  if (fHadCorr > 1) {
258  energyclus = ApplyHadCorrAllTracks(fClusCont->GetCurrentID(), fHadCorr - 1);
259  }
260  else if (fHadCorr > 0) {
261  energyclus = ApplyHadCorrOneTrack(fClusCont->GetCurrentID(), fHadCorr);
262  }
263  else {
264  energyclus = cluster->GetNonLinCorrEnergy();
265  }
266 
267  if (energyclus < 0) energyclus = 0;
268 
269  cluster->SetHadCorrEnergy(energyclus);
270 
271  if (fCreateHisto) fHistEafter->Fill(fCent, energyclus);
272 
273  }
274 
275  return kTRUE;
276 }
277 
278 //________________________________________________________________________
280 {
281  // Get momenum bin.
282 
283  UInt_t pbin=0;
284  if (p<0.5)
285  pbin=0;
286  else if (p>=0.5 && p<1.0)
287  pbin=1;
288  else if (p>=1.0 && p<1.5)
289  pbin=2;
290  else if (p>=1.5 && p<2.)
291  pbin=3;
292  else if (p>=2. && p<3.)
293  pbin=4;
294  else if (p>=3. && p<4.)
295  pbin=5;
296  else if (p>=4. && p<5.)
297  pbin=6;
298  else if (p>=5. && p<8.)
299  pbin=7;
300  else
301  pbin=8;
302 
303  return pbin;
304 }
305 
306 //________________________________________________________________________
308 {
309  // Get sigma in eta.
310 
311  Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
312  return 2.0*EtaSigma[pbin];
313 }
314 
315 //________________________________________________________________________
317 {
318  // Get phi mean.
319 
320  if (centbin==0) {
321  Double_t PhiMean[9]={0.036,
322  0.021,
323  0.0121,
324  0.0084,
325  0.0060,
326  0.0041,
327  0.0031,
328  0.0022,
329  0.001};
330  return PhiMean[pbin];
331  } else if (centbin==1) {
332  Double_t PhiMean[9]={0.036,
333  0.021,
334  0.0121,
335  0.0084,
336  0.0060,
337  0.0041,
338  0.0031,
339  0.0022,
340  0.001};
341  return PhiMean[pbin];
342  } else if (centbin==2) {
343  Double_t PhiMean[9]={0.036,
344  0.021,
345  0.0121,
346  0.0084,
347  0.0060,
348  0.0041,
349  0.0031,
350  0.0022,
351  0.001};
352  return PhiMean[pbin];
353  } else if (centbin==3) {
354  Double_t PhiMean[9]={0.036,
355  0.021,
356  0.0121,
357  0.0084,
358  0.0060,
359  0.0041,
360  0.0031,
361  0.0022,
362  0.001};
363  return PhiMean[pbin];
364  } else if (centbin==4) {
365  Double_t PhiMean[9]={0.036,
366  0.021,
367  0.0127,
368  0.0089,
369  0.0068,
370  0.0049,
371  0.0038,
372  0.0028,
373  0.0018};
374  return PhiMean[pbin]*(-1.);
375  } else if (centbin==5) {
376  Double_t PhiMean[9]={0.036,
377  0.021,
378  0.0127,
379  0.0089,
380  0.0068,
381  0.0048,
382  0.0038,
383  0.0028,
384  0.0018};
385  return PhiMean[pbin]*(-1.);
386  } else if (centbin==6) {
387  Double_t PhiMean[9]={0.036,
388  0.021,
389  0.0127,
390  0.0089,
391  0.0068,
392  0.0045,
393  0.0035,
394  0.0028,
395  0.0018};
396  return PhiMean[pbin]*(-1.);
397  } else if (centbin==7) {
398  Double_t PhiMean[9]={0.036,
399  0.021,
400  0.0127,
401  0.0089,
402  0.0068,
403  0.0043,
404  0.0035,
405  0.0028,
406  0.0018};
407  return PhiMean[pbin]*(-1.);
408  }
409 
410  return 0;
411 }
412 
413 //________________________________________________________________________
415 {
416  // Get phi sigma.
417 
418  if (centbin==0) {
419  Double_t PhiSigma[9]={0.0221,
420  0.0128,
421  0.0074,
422  0.0064,
423  0.0059,
424  0.0055,
425  0.0052,
426  0.0049,
427  0.0045};
428  return 2.*PhiSigma[pbin];
429  } else if (centbin==1) {
430  Double_t PhiSigma[9]={0.0217,
431  0.0120,
432  0.0076,
433  0.0066,
434  0.0062,
435  0.0058,
436  0.0054,
437  0.0054,
438  0.0045};
439  return 2.*PhiSigma[pbin];
440  } else if (centbin==2) {
441  Double_t PhiSigma[9]={0.0211,
442  0.0124,
443  0.0080,
444  0.0070,
445  0.0067,
446  0.0061,
447  0.0059,
448  0.0054,
449  0.0047};
450  return 2.*PhiSigma[pbin];
451  } else if (centbin==3) {
452  Double_t PhiSigma[9]={0.0215,
453  0.0124,
454  0.0082,
455  0.0073,
456  0.0069,
457  0.0064,
458  0.0060,
459  0.0055,
460  0.0047};
461  return 2.*PhiSigma[pbin];
462  } else if (centbin==4) {
463  Double_t PhiSigma[9]={0.0199,
464  0.0108,
465  0.0072,
466  0.0071,
467  0.0060,
468  0.0055,
469  0.0052,
470  0.0049,
471  0.0045};
472  return 2.*PhiSigma[pbin];
473  } else if (centbin==5) {
474  Double_t PhiSigma[9]={0.0200,
475  0.0110,
476  0.0074,
477  0.0071,
478  0.0064,
479  0.0059,
480  0.0055,
481  0.0052,
482  0.0045};
483  return 2.*PhiSigma[pbin];
484  } else if (centbin==6) {
485  Double_t PhiSigma[9]={0.0202,
486  0.0113,
487  0.0077,
488  0.0071,
489  0.0069,
490  0.0064,
491  0.0060,
492  0.0055,
493  0.0050};
494  return 2.*PhiSigma[pbin];
495  } else if (centbin==7) {
496  Double_t PhiSigma[9]={0.0205,
497  0.0113,
498  0.0080,
499  0.0074,
500  0.0078,
501  0.0067,
502  0.0062,
503  0.0055,
504  0.0050};
505  return 2.*PhiSigma[pbin];
506  }
507 
508  return 0;
509 }
510 
511 //________________________________________________________________________
513  Double_t &totalTrkP, Int_t &Nmatches, Double_t &trkPMCfrac, Int_t &NMCmatches)
514 {
515  // Do the loop over matched tracks for the cluster.
516 
517  AliVCluster* cluster = fClusCont->GetCluster(icluster);
518 
519  if (!cluster) return;
520 
521  // loop over matched tracks
522  Int_t Ntrks = Ntrks = cluster->GetNTracksMatched();
523  for (Int_t i = 0; i < Ntrks; ++i) {
524  AliVTrack* track = 0;
525 
526  if (fEsdMode) {
527  Int_t itrack = cluster->GetTrackMatchedIndex(i);
528  if (itrack >= 0) track = static_cast<AliVTrack*>(fPartCont->GetAcceptParticle(itrack));
529  }
530  else {
531  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(i));
532  UInt_t rejectionReason = 0;
533  if (!fPartCont->AcceptParticle(track, rejectionReason)) track = 0;
534  }
535 
536  if (!track) continue;
537 
538  Double_t etadiff = 999;
539  Double_t phidiff = 999;
540  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
541  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(etadiff, phidiff);
542 
543  // check if track also points to cluster
544  if (fDoTrackClus && (track->GetEMCALcluster() != icluster)) continue;
545 
546  Double_t mom = track->P();
547  UInt_t mombin = GetMomBin(mom);
548  Int_t centbinch = fCentBin;
549 
550  if (track->Charge() < 0) centbinch += fNcentBins;
551 
552  if (fCreateHisto) {
553  Int_t etabin = 0;
554  if (track->Eta() > 0) etabin=1;
555  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
556  fHistMatchEtaPhiAll->Fill(etadiff, phidiff);
557  }
558 
559  Double_t etaCut = 0.0;
560  Double_t phiCutlo = 0.0;
561  Double_t phiCuthi = 0.0;
562 
563  if (fPhiMatch > 0) {
564  phiCutlo = -fPhiMatch;
565  phiCuthi = +fPhiMatch;
566  }
567  else {
568  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
569  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
570  }
571 
572  if (fEtaMatch > 0) {
573  etaCut = fEtaMatch;
574  }
575  else {
576  etaCut = GetEtaSigma(mombin);
577  }
578 
579  if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
580  if (track->GetLabel() > fMinMCLabel) {
581  ++NMCmatches;
582  trkPMCfrac += mom;
583  }
584  ++Nmatches;
585  totalTrkP += mom;
586 
587  if (fCreateHisto) {
588  if (fHadCorr > 1) {
589  Double_t dphi = 0;
590  Double_t deta = 0;
591  GetEtaPhiDiff(track, cluster, dphi, deta);
592  Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
593  Double_t energyclus = cluster->GetNonLinCorrEnergy();
594  fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
595  }
596  }
597  }
598  }
599 
600  if (totalTrkP > 0) trkPMCfrac /= totalTrkP;
601 }
602 
603 //________________________________________________________________________
605 {
606  // Apply the hadronic correction with one track only.
607 
608  AliVCluster* cluster = fClusCont->GetCluster(icluster);
609 
610  Double_t energyclus = cluster->GetNonLinCorrEnergy();
611 
612  AliVTrack* track = 0;
613 
614  if (cluster->GetNTracksMatched() > 0) {
615  if (fEsdMode) {
616  Int_t itrack = cluster->GetTrackMatchedIndex(0);
617  if (itrack >= 0) track = static_cast<AliVTrack*>(fPartCont->GetAcceptParticle(itrack));
618  }
619  else {
620  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
621  UInt_t rejectionReason = 0;
622  if (!fPartCont->AcceptParticle(track, rejectionReason)) track = 0;
623  }
624  }
625 
626  if (!track || track->P() < 1e-6) return energyclus;
627 
628  Double_t mom = track->P();
629 
630  Double_t dEtaMin = 1e9;
631  Double_t dPhiMin = 1e9;
632  GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
633 
634  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(dEtaMin, dPhiMin);
635 
636  // check if track also points to cluster
637  Int_t cid = track->GetEMCALcluster();
638  if (fDoTrackClus && (cid != icluster)) return energyclus;
639 
640  UInt_t mombin = GetMomBin(mom);
641  Int_t centbinch = fCentBin;
642  if (track->Charge() < 0) centbinch += fNcentBins;
643 
644  // plot some histograms if switched on
645  if (fCreateHisto) {
646  Int_t etabin = 0;
647  if(track->Eta() > 0) etabin = 1;
648 
649  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
650  fHistMatchEtaPhiAll->Fill(dEtaMin, dPhiMin);
651 
652  if (mom > 0) {
653  Double_t etadiff = 0;
654  Double_t phidiff = 0;
655  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
656  Double_t dRmin = TMath::Sqrt(etadiff*etadiff + phidiff*phidiff);
657  fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
658  fHistEoPCent->Fill(fCent, energyclus / mom);
659  fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
660  }
661  }
662 
663  // define eta/phi cuts
664  Double_t etaCut = 0.0;
665  Double_t phiCutlo = 0.0;
666  Double_t phiCuthi = 0.0;
667  if (fPhiMatch > 0) {
668  phiCutlo = -fPhiMatch;
669  phiCuthi = +fPhiMatch;
670  }
671  else {
672  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
673  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
674  }
675  if (fEtaMatch > 0) {
676  etaCut = fEtaMatch;
677  }
678  else {
679  etaCut = GetEtaSigma(mombin);
680  }
681 
682  // apply the correction if the track is in the eta/phi window
683  if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
684  energyclus -= hadCorr * mom;
685  }
686 
687  return energyclus;
688 }
689 
690 //________________________________________________________________________
692 {
693  // Apply the hadronic correction with all tracks.
694 
695  AliVCluster* cluster = fClusCont->GetCluster(icluster);
696 
697  Double_t energyclus = cluster->GetNonLinCorrEnergy();
698  Double_t cNcells = cluster->GetNCells();
699 
700  Double_t totalTrkP = 0.0; // count total track momentum
701  Int_t Nmatches = 0; // count total number of matches
702 
703  Double_t trkPMCfrac = 0.0; // count total track momentum
704  Int_t NMCmatches = 0; // count total number of matches
705 
706  // do the loop over the matched tracks and get the number of matches and the total momentum
707  DoMatchedTracksLoop(icluster, totalTrkP, Nmatches, trkPMCfrac, NMCmatches);
708 
709  Double_t Esub = hadCorr * totalTrkP;
710 
711  if (Esub > energyclus) Esub = energyclus;
712 
713  // applying Peter's proposed algorithm
714  // never subtract the full energy of the cluster
715  Double_t clusEexcl = fEexclCell * cNcells;
716  if (energyclus < clusEexcl) clusEexcl = energyclus;
717  if ((energyclus - Esub) < clusEexcl) Esub = (energyclus - clusEexcl);
718 
719  // embedding
720  Double_t EsubMC = 0;
721  Double_t EsubBkg = 0;
722  Double_t EclusMC = 0;
723  Double_t EclusBkg = 0;
724  Double_t EclusCorr = 0;
725  Double_t EclusMCcorr = 0;
726  Double_t EclusBkgcorr = 0;
727  Double_t overSub = 0;
728  if (fIsEmbedded) {
729  EsubMC = hadCorr * totalTrkP * trkPMCfrac;
730  EsubBkg = hadCorr * totalTrkP - EsubMC;
731  EclusMC = energyclus * cluster->GetMCEnergyFraction();
732  EclusBkg = energyclus - EclusMC;
733 
734  if (energyclus > Esub)
735  EclusCorr = energyclus - Esub;
736 
737  if (EclusMC > EsubMC)
738  EclusMCcorr = EclusMC - EsubMC;
739 
740  if (EclusBkg > EsubBkg)
741  EclusBkgcorr = EclusBkg - EsubBkg;
742 
743  overSub = EclusMCcorr + EclusBkgcorr - EclusCorr;
744  }
745 
746  // plot some histograms if switched on
747  if (fCreateHisto) {
748  fHistNMatchCent->Fill(fCent, Nmatches);
749  fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
750 
751  if (Nmatches > 0) fHistNclusMatchvsCent->Fill(fCent);
752 
753  if (Nmatches < 3) {
754  fHistNCellsEnergy[fCentBin][Nmatches]->Fill(energyclus, cNcells);
755  }
756  else {
757  fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
758  }
759 
760  if (totalTrkP > 0) {
761  Double_t EoP = energyclus / totalTrkP;
762  fHistEoPCent->Fill(fCent, EoP);
763  fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
764 
765  fHistEsubPchRatAll[fCentBin]->Fill(totalTrkP, Esub / totalTrkP);
766 
767  if (Nmatches == 1) {
768  AliVTrack* track = 0;
769  if (fEsdMode) {
770  Int_t itrack = cluster->GetTrackMatchedIndex(0);
771  if (itrack >= 0) track = static_cast<AliVTrack*>(fPartCont->GetAcceptParticle(itrack));
772  }
773  else {
774  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
775  UInt_t rejectionReason = 0;
776  if (!fPartCont->AcceptParticle(track, rejectionReason)) track = 0;
777  }
778  if (track) {
779  Int_t centbinchm = fCentBin;
780  if (track->Charge() < 0) centbinchm += fNcentBins;
781  fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
782  fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
783  }
784  }
785 
786  if (fIsEmbedded) {
787  fHistOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
788 
789  if (cluster->GetMCEnergyFraction() > 0.95)
790  fHistOversubMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
791  else if (cluster->GetMCEnergyFraction() < 0.05)
792  fHistOversubNonMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
793 
794  if (trkPMCfrac < 0.05)
795  fHistNonEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
796  else if (trkPMCfrac > 0.95)
797  fHistEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
798  }
799  }
800  }
801 
802  if (fIsEmbedded && fDoExact) {
803  Esub -= overSub;
804  if (EclusBkgcorr + EclusMCcorr > 0) {
805  Double_t newfrac = EclusMCcorr / (EclusBkgcorr + EclusMCcorr);
806  cluster->SetMCEnergyFraction(newfrac);
807  }
808  }
809 
810  // apply the correction
811  energyclus -= Esub;
812 
813  return energyclus;
814 }
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 (only 1 match) ...
void AddContainer(inputObjectType type)
void SetParticlePtCut(Double_t cut)
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 > 95%)
TH2 * fHistOversub[4]
Over-subtracted energy / cluster energy (cluster MC energy fraction < 5%)
Int_t fCentBin
!event centrality bin
TH2 * fHistMatchEvsP[4]
deta vs. dphi of all cluster-track pairs (cl loop)
TH2 * fHistEsubPchRat[8]
Esub vs. total momentum of matched tracks (only 1 match)
TH1 * fHistEafter
average energy of clusters before 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
void SetClusNonLinCorrEnergyCut(Double_t cut)
Int_t fMinMCLabel
minimum MC label value for the tracks/clusters being considered MC particles
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
void GetProperty(std::string propertyName, T &property, bool requiredProperty=true, std::string correctionName="")
Retrieve property.
TH2 * fHistEoPCent
average energy of clusters after correction vs. centrality
virtual Bool_t AcceptParticle(const AliVParticle *vp, UInt_t &rejectionReason) const
TH2 * fHistNMatchEnergy[4]
cluster energy vs. track momentum of matched pairs
TH2 * fHistNonEmbTrackMatchesOversub[4]
Over-subtracted energy / cluster energy with embedded track matches (non-embedded matches < 5%) ...
AliVCluster * GetCluster(Int_t i) const
Double_t fMaxBinPt
max pt in histograms
TList * fOutput
! list of output histograms
TH1 * fHistEbefore
n clusters matched to some track vs. 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
TH2 * fHistEmbTrackMatchesOversub[4]
Esub/momentum of matched tracks vs. total momentum of matched tracks (all number of matches) ...
void SetClusPtCut(Double_t cut)
void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
TH2 * fHistMatchEtaPhiAllTr
deta vs. dphi of matched cluster-track pairs
bool Bool_t
Definition: External.C:53
TH1 * fHistEsubPch[8]
n clusters macthed to some track (tracks allowed to match more than one cluster)
void SetClusECut(Double_t cut)
TH2 * fHistOversubMCClusters[4]
Over-subtracted energy / cluster energy with non-embedded track matches (embedded matches < 5%) ...
AliParticleContainer * fPartCont
! pointer to the track/particle container
AliVCluster * GetNextAcceptCluster()
TH2 * fHistMatchEtaPhiAllCl
deta vs. dphi of all cluster-track pairs (tr loop)
Bool_t fIsEmbedded
trigger, embedded signal