AliPhysics  608b256 (608b256)
AliHadCorrTask.cxx
Go to the documentation of this file.
1 //
2 // Hadronic correction task.
3 //
4 // Author: R.Reed, C.Loizides, S.Aiola
5 
6 #include <TChain.h>
7 #include <TClonesArray.h>
8 #include <TH1F.h>
9 #include <TH2F.h>
10 #include <TList.h>
11 
12 #include "AliAnalysisManager.h"
13 #include "AliAODEvent.h"
14 #include "AliAODCaloCluster.h"
15 #include "AliESDCaloCluster.h"
16 #include "AliVTrack.h"
17 #include "AliVEventHandler.h"
18 #include "AliEMCALGeometry.h"
19 #include "AliParticleContainer.h"
20 #include "AliClusterContainer.h"
21 
22 #include "AliHadCorrTask.h"
23 
24 ClassImp(AliHadCorrTask)
25 
26 //________________________________________________________________________
29  fOutCaloName(),
30  fPhiMatch(0.05),
31  fEtaMatch(0.025),
32  fDoTrackClus(kTRUE),
33  fHadCorr(0),
34  fEexclCell(0),
35  fDoExact(kFALSE),
36  fEsdMode(kTRUE),
37  fOutClusters(0),
38  fHistMatchEtaPhiAll(0),
39  fHistMatchEtaPhiAllTr(0),
40  fHistMatchEtaPhiAllCl(0),
41  fHistNclusvsCent(0),
42  fHistNclusMatchvsCent(0),
43  fHistEbefore(0),
44  fHistEafter(0),
45  fHistEoPCent(0),
46  fHistNMatchCent(0),
47  fHistNClusMatchCent(0)
48 {
49  // Default constructor.
50 
51  for(Int_t i=0; i<10; i++) {
52  fHistEsubPch[i] = 0;
53  fHistEsubPchRat[i] = 0;
54 
55  if (i<5) {
56  fHistEsubPchRatAll[i] = 0;
57 
58  fHistMatchEvsP[i] = 0;
59  fHistMatchdRvsEP[i] = 0;
60  fHistNMatchEnergy[i] = 0;
61 
62  fHistEmbTrackMatchesOversub[i] = 0;
63  fHistNonEmbTrackMatchesOversub[i] = 0;
64  fHistOversubMCClusters[i] = 0;
65  fHistOversubNonMCClusters[i] = 0;
66  fHistOversub[i] = 0;
67 
68  for(Int_t j=0; j<4; j++)
69  fHistNCellsEnergy[i][j] = 0;
70  }
71 
72  for(Int_t j=0; j<9; j++) {
73  for(Int_t k=0; k<2; k++)
74  fHistMatchEtaPhi[i][j][k] = 0;
75  }
76  }
77 }
78 
79 //________________________________________________________________________
80 AliHadCorrTask::AliHadCorrTask(const char *name, Bool_t histo) :
81  AliAnalysisTaskEmcal(name, histo),
82  fOutCaloName("CaloClustersCorr"),
83  fPhiMatch(0.05),
84  fEtaMatch(0.025),
85  fDoTrackClus(kTRUE),
86  fHadCorr(0),
87  fEexclCell(0),
88  fDoExact(kFALSE),
89  fEsdMode(kTRUE),
90  fOutClusters(0),
91  fHistMatchEtaPhiAll(0),
92  fHistMatchEtaPhiAllTr(0),
93  fHistMatchEtaPhiAllCl(0),
94  fHistNclusvsCent(0),
95  fHistNclusMatchvsCent(0),
96  fHistEbefore(0),
97  fHistEafter(0),
98  fHistEoPCent(0),
99  fHistNMatchCent(0),
100  fHistNClusMatchCent(0)
101 {
102  // Standard constructor.
103 
104  for(Int_t i=0; i<10; i++) {
105  fHistEsubPch[i] = 0;
106  fHistEsubPchRat[i] = 0;
107 
108  if (i<5) {
109  fHistEsubPchRatAll[i] = 0;
110 
111  fHistMatchEvsP[i] = 0;
112  fHistMatchdRvsEP[i] = 0;
113  fHistNMatchEnergy[i] = 0;
114 
117  fHistOversubMCClusters[i] = 0;
119  fHistOversub[i] = 0;
120 
121  for(Int_t j=0; j<4; j++)
122  fHistNCellsEnergy[i][j] = 0;
123  }
124 
125  for(Int_t j=0; j<9; j++) {
126  for(Int_t k=0; k<2; k++)
127  fHistMatchEtaPhi[i][j][k] = 0;
128  }
129  }
130 
132 
133  fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
134 }
135 
136 //________________________________________________________________________
138 {
139  // Destructor
140 }
141 
142 //________________________________________________________________________
144 {
145  // Get momenum bin.
146 
147  UInt_t pbin=0;
148  if (p<0.5)
149  pbin=0;
150  else if (p>=0.5 && p<1.0)
151  pbin=1;
152  else if (p>=1.0 && p<1.5)
153  pbin=2;
154  else if (p>=1.5 && p<2.)
155  pbin=3;
156  else if (p>=2. && p<3.)
157  pbin=4;
158  else if (p>=3. && p<4.)
159  pbin=5;
160  else if (p>=4. && p<5.)
161  pbin=6;
162  else if (p>=5. && p<8.)
163  pbin=7;
164  else
165  pbin=8;
166 
167  return pbin;
168 }
169 
170 //________________________________________________________________________
172 {
173  // Get sigma in eta.
174 
175  Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
176  return 2.0*EtaSigma[pbin];
177 }
178 
179 //________________________________________________________________________
181 {
182  // Get phi mean.
183 
184  if (centbin==0) {
185  Double_t PhiMean[9]={0.036,
186  0.021,
187  0.0121,
188  0.0084,
189  0.0060,
190  0.0041,
191  0.0031,
192  0.0022,
193  0.001};
194  return PhiMean[pbin];
195  } else if (centbin==1) {
196  Double_t PhiMean[9]={0.036,
197  0.021,
198  0.0121,
199  0.0084,
200  0.0060,
201  0.0041,
202  0.0031,
203  0.0022,
204  0.001};
205  return PhiMean[pbin];
206  } else if (centbin==2) {
207  Double_t PhiMean[9]={0.036,
208  0.021,
209  0.0121,
210  0.0084,
211  0.0060,
212  0.0041,
213  0.0031,
214  0.0022,
215  0.001};
216  return PhiMean[pbin];
217  } else if (centbin==3) {
218  Double_t PhiMean[9]={0.036,
219  0.021,
220  0.0121,
221  0.0084,
222  0.0060,
223  0.0041,
224  0.0031,
225  0.0022,
226  0.001};
227  return PhiMean[pbin];
228  } else if (centbin==4) {
229  Double_t PhiMean[9]={0.036,
230  0.021,
231  0.0127,
232  0.0089,
233  0.0068,
234  0.0049,
235  0.0038,
236  0.0028,
237  0.0018};
238  return PhiMean[pbin]*(-1.);
239  } else if (centbin==5) {
240  Double_t PhiMean[9]={0.036,
241  0.021,
242  0.0127,
243  0.0089,
244  0.0068,
245  0.0048,
246  0.0038,
247  0.0028,
248  0.0018};
249  return PhiMean[pbin]*(-1.);
250  } else if (centbin==6) {
251  Double_t PhiMean[9]={0.036,
252  0.021,
253  0.0127,
254  0.0089,
255  0.0068,
256  0.0045,
257  0.0035,
258  0.0028,
259  0.0018};
260  return PhiMean[pbin]*(-1.);
261  } else if (centbin==7) {
262  Double_t PhiMean[9]={0.036,
263  0.021,
264  0.0127,
265  0.0089,
266  0.0068,
267  0.0043,
268  0.0035,
269  0.0028,
270  0.0018};
271  return PhiMean[pbin]*(-1.);
272  }
273 
274  return 0;
275 }
276 
277 //________________________________________________________________________
279 {
280  // Get phi sigma.
281 
282  if (centbin==0) {
283  Double_t PhiSigma[9]={0.0221,
284  0.0128,
285  0.0074,
286  0.0064,
287  0.0059,
288  0.0055,
289  0.0052,
290  0.0049,
291  0.0045};
292  return 2.*PhiSigma[pbin];
293  } else if (centbin==1) {
294  Double_t PhiSigma[9]={0.0217,
295  0.0120,
296  0.0076,
297  0.0066,
298  0.0062,
299  0.0058,
300  0.0054,
301  0.0054,
302 0.0045};
303  return 2.*PhiSigma[pbin];
304  } else if (centbin==2) {
305  Double_t PhiSigma[9]={0.0211,
306  0.0124,
307  0.0080,
308  0.0070,
309  0.0067,
310  0.0061,
311  0.0059,
312  0.0054,
313  0.0047};
314  return 2.*PhiSigma[pbin];
315  } else if (centbin==3) {
316  Double_t PhiSigma[9]={0.0215,
317  0.0124,
318  0.0082,
319  0.0073,
320  0.0069,
321  0.0064,
322  0.0060,
323  0.0055,
324  0.0047};
325  return 2.*PhiSigma[pbin];
326  } else if (centbin==4) {
327  Double_t PhiSigma[9]={0.0199,
328  0.0108,
329  0.0072,
330  0.0071,
331  0.0060,
332  0.0055,
333  0.0052,
334  0.0049,
335  0.0045};
336  return 2.*PhiSigma[pbin];
337  } else if (centbin==5) {
338  Double_t PhiSigma[9]={0.0200,
339  0.0110,
340  0.0074,
341  0.0071,
342  0.0064,
343  0.0059,
344  0.0055,
345  0.0052,
346  0.0045};
347  return 2.*PhiSigma[pbin];
348  } else if (centbin==6) {
349  Double_t PhiSigma[9]={0.0202,
350  0.0113,
351  0.0077,
352  0.0071,
353  0.0069,
354  0.0064,
355  0.0060,
356  0.0055,
357  0.0050};
358  return 2.*PhiSigma[pbin];
359  } else if (centbin==7) {
360  Double_t PhiSigma[9]={0.0205,
361  0.0113,
362  0.0080,
363  0.0074,
364  0.0078,
365  0.0067,
366  0.0062,
367  0.0055,
368  0.0050};
369  return 2.*PhiSigma[pbin];
370  }
371 
372  return 0;
373 }
374 
375 //________________________________________________________________________
377 {
378  // Create my user objects.
379 
381 
382  if (!fCreateHisto) return;
383 
384  TString name;
385  TString temp;
386 
387  const Int_t nCentChBins = fNcentBins * 2;
388 
389  fHistMatchEtaPhiAll = new TH2F("fHistMatchEtaPhiAll", "fHistMatchEtaPhiAll;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
391 
392  fHistMatchEtaPhiAllTr = new TH2F("fHistMatchEtaPhiAllTr", "fHistMatchEtaPhiAllTr;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
394 
395  fHistMatchEtaPhiAllCl = new TH2F("fHistMatchEtaPhiAllCl", "fHistMatchEtaPhiAllCl;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
397 
398  for(Int_t icent=0; icent<nCentChBins; ++icent) {
399  for(Int_t ipt=0; ipt<9; ++ipt) {
400  for(Int_t ieta=0; ieta<2; ++ieta) {
401  name = Form("fHistMatchEtaPhi_%i_%i_%i",icent,ipt,ieta);
402  fHistMatchEtaPhi[icent][ipt][ieta] = new TH2F(name, name, fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
403  fHistMatchEtaPhi[icent][ipt][ieta]->SetXTitle("#Delta#eta");
404  fHistMatchEtaPhi[icent][ipt][ieta]->SetYTitle("#Delta#phi");
405  fOutput->Add(fHistMatchEtaPhi[icent][ipt][ieta]);
406  }
407  }
408 
409  name = Form("fHistEsubPch_%i",icent);
410  temp = Form("%s (Nmatches==1)",name.Data());
411  fHistEsubPch[icent]=new TH1F(name, temp, fNbins, fMinBinPt, fMaxBinPt);
412  fHistEsubPch[icent]->SetXTitle("#sum p (GeV) weighted with E_{sub}");
413  fOutput->Add(fHistEsubPch[icent]);
414 
415  name = Form("fHistEsubPchRat_%i",icent);
416  temp = Form("%s (Nmatches==1)",name.Data());
417  fHistEsubPchRat[icent]=new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
418  fHistEsubPchRat[icent]->SetXTitle("#Sigma p (GeV)");
419  fHistEsubPchRat[icent]->SetYTitle("E_{sub} / #sum p");
420  fOutput->Add(fHistEsubPchRat[icent]);
421 
422  if (icent<fNcentBins) {
423  for(Int_t itrk=0; itrk<4; ++itrk) {
424  name = Form("fHistNCellsEnergy_%i_%i",icent,itrk);
425  temp = Form("%s (Nmatches==%d);N_{cells};E_{clus} (GeV)",name.Data(),itrk);
426  fHistNCellsEnergy[icent][itrk] = new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, 101, -0.5, 100.5);
427  fOutput->Add(fHistNCellsEnergy[icent][itrk]);
428  }
429 
430  name = Form("fHistEsubPchRatAll_%i",icent);
431  temp = Form("%s (all Nmatches)",name.Data());
432  fHistEsubPchRatAll[icent]=new TH2F(name, temp, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
433  fHistEsubPchRatAll[icent]->SetXTitle("#Sigma p (GeV)");
434  fHistEsubPchRatAll[icent]->SetYTitle("E_{sub} / #sum p");
435  fOutput->Add(fHistEsubPchRatAll[icent]);
436 
437  name = Form("fHistMatchEvsP_%i",icent);
438  temp = Form("%s (all Nmatches)",name.Data());
439  fHistMatchEvsP[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins*2, 0., 10.);
440  fHistMatchEvsP[icent]->SetXTitle("E_{clus} (GeV)");
441  fHistMatchEvsP[icent]->SetYTitle("E_{clus} / #sum p");
442  fOutput->Add(fHistMatchEvsP[icent]);
443 
444  name = Form("fHistMatchdRvsEP_%i",icent);
445  temp = Form("%s (all Nmatches)",name.Data());
446  fHistMatchdRvsEP[icent] = new TH2F(name, temp, fNbins, 0., 0.2, fNbins*2, 0., 10.);
447  fHistMatchdRvsEP[icent]->SetXTitle("#Delta R between track and cluster");
448  fHistMatchdRvsEP[icent]->SetYTitle("E_{clus} / p");
449  fOutput->Add(fHistMatchdRvsEP[icent]);
450 
451  name = Form("fHistNMatchEnergy_%i",icent);
452  fHistNMatchEnergy[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, 101, -0.5, 100.5);
453  fHistNMatchEnergy[icent]->SetXTitle("E_{clus} (GeV)");
454  fHistNMatchEnergy[icent]->SetYTitle("N_{matches}");
455  fOutput->Add(fHistNMatchEnergy[icent]);
456  }
457  }
458 
459  fHistNclusvsCent = new TH1F("Nclusvscent", "NclusVsCent; Cent (%)", 100, 0, 100);
460  fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent (all Nmatches); Cent (%)", 100, 0, 100);
461  fHistEbefore = new TH1F("Ebefore", "Ebefore; Cent (%); E_{clus} (GeV)", 100, 0, 100);
462  fHistEafter = new TH1F("Eafter", "Eafter; Cent (%); E_{clus} (GeV)", 100, 0, 100);
463  fHistEoPCent = new TH2F("EoPCent", "EoPCent; Cent (%); E_{clus} / #sum p", 100, 0, 100, fNbins*2, 0, 10);
464  fHistNMatchCent = new TH2F("NMatchesCent", "NMatchesCent; Cent (%); Nmatches", 100, 0, 100, 11, -0.5, 10.5);
465  fHistNClusMatchCent = new TH2F("NClusMatchesCent", "NClusMatchesCent; Cent (%); Nmatches", 100, 0, 100, 11, -0.5, 10.5);
466 
469  fOutput->Add(fHistEbefore);
470  fOutput->Add(fHistEafter);
471  fOutput->Add(fHistEoPCent);
472  fOutput->Add(fHistNMatchCent);
474 
475  if (fIsEmbedded) {
476  for(Int_t icent=0; icent<fNcentBins; ++icent) {
477  name = Form("fHistEmbTrackMatchesOversub_%d",icent);
478  fHistEmbTrackMatchesOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
479  fHistEmbTrackMatchesOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
480  fHistEmbTrackMatchesOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
482 
483  name = Form("fHistNonEmbTrackMatchesOversub_%d",icent);
484  fHistNonEmbTrackMatchesOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
485  fHistNonEmbTrackMatchesOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
486  fHistNonEmbTrackMatchesOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
488 
489  name = Form("fHistOversubMCClusters_%d",icent);
490  fHistOversubMCClusters[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
491  fHistOversubMCClusters[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
492  fHistOversubMCClusters[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
493  fOutput->Add(fHistOversubMCClusters[icent]);
494 
495  name = Form("fHistOversubNonMCClusters_%d",icent);
496  fHistOversubNonMCClusters[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
497  fHistOversubNonMCClusters[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
498  fHistOversubNonMCClusters[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
499  fOutput->Add(fHistOversubNonMCClusters[icent]);
500 
501  name = Form("fHistOversub_%d",icent);
502  fHistOversub[icent] = new TH2F(name, name, fNbins, fMinBinPt, fMaxBinPt, fNbins, 0, 1.2);
503  fHistOversub[icent]->GetXaxis()->SetTitle("E_{clus}^{raw} (GeV)");
504  fHistOversub[icent]->GetYaxis()->SetTitle("E_{oversub} / E_{clus}^{raw}");
505  fOutput->Add(fHistOversub[icent]);
506  }
507  }
508 
509  PostData(1, fOutput);
510 }
511 
512 //________________________________________________________________________
514 {
515  // Initialize the analysis.
516 
518  if (!tracks) {
519  AliError(Form("%s: This task needs a particle container!", GetName()));
520  return;
521  }
522 
523  TClass trackClass(tracks->GetClassName());
524  if (!trackClass.InheritsFrom("AliVTrack")) {
525  tracks->SetClassName("AliVTrack"); // enforce only AliVTrack and derived classes
526  }
527 
529  if (!clusters) {
530  AliError(Form("%s: This task needs a cluster container!", GetName()));
531  return;
532  }
533 
534  TClass clusterClass(clusters->GetClassName());
535  if (!clusterClass.InheritsFrom("AliVCluster")) {
536  clusters->SetClassName("AliVCluster"); // enforce only AliVCluster and derived classes
537  }
538 
539  // Do base class initializations and if it fails -> bail out
541  if (!fLocalInitialized) return;
542 
543  if (dynamic_cast<AliAODEvent*>(InputEvent())) fEsdMode = kFALSE;
544 
545  if (!fOutCaloName.IsNull()) {
546  if (fEsdMode) {
547  fOutClusters = new TClonesArray("AliESDCaloCluster");
548  }
549  else {
550  fOutClusters = new TClonesArray("AliAODCaloCluster");
551  }
552 
553  fOutClusters->SetName(fOutCaloName);
555  }
556 }
557 
558 //________________________________________________________________________
560  Double_t &totalTrkP, Int_t &Nmatches, Double_t &trkPMCfrac, Int_t &NMCmatches)
561 {
562  // Do the loop over matched tracks for the cluster.
563 
566 
567  AliVCluster* cluster = clusters->GetCluster(icluster);
568 
569  if (!cluster) return;
570 
571  // loop over matched tracks
572  Int_t Ntrks = cluster->GetNTracksMatched();
573  for (Int_t i = 0; i < Ntrks; ++i) {
574  AliVTrack* track = 0;
575 
576  if (fEsdMode) {
577  Int_t itrack = cluster->GetTrackMatchedIndex(i);
578  if (itrack >= 0) track = static_cast<AliVTrack*>(tracks->GetAcceptParticle(itrack));
579  }
580  else {
581  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(i));
582  UInt_t rejectionReason = 0;
583  if (!tracks->AcceptParticle(track, rejectionReason)) track = 0;
584  }
585 
586  if (!track) continue;
587 
588  Double_t etadiff = 999;
589  Double_t phidiff = 999;
590  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
591  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(etadiff, phidiff);
592 
593  // check if track also points to cluster
594  if (fDoTrackClus && (track->GetEMCALcluster() != icluster)) continue;
595 
596  Double_t mom = track->P();
597  UInt_t mombin = GetMomBin(mom);
598  Int_t centbinch = fCentBin;
599 
600  if (track->Charge() < 0) centbinch += fNcentBins;
601 
602  if (fCreateHisto) {
603  Int_t etabin = 0;
604  if (track->Eta() > 0) etabin=1;
605  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
606  fHistMatchEtaPhiAll->Fill(etadiff, phidiff);
607  }
608 
609  Double_t etaCut = 0.0;
610  Double_t phiCutlo = 0.0;
611  Double_t phiCuthi = 0.0;
612 
613  if (fPhiMatch > 0) {
614  phiCutlo = -fPhiMatch;
615  phiCuthi = +fPhiMatch;
616  }
617  else {
618  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
619  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
620  }
621 
622  if (fEtaMatch > 0) {
623  etaCut = fEtaMatch;
624  }
625  else {
626  etaCut = GetEtaSigma(mombin);
627  }
628 
629  if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
630  if (track->GetLabel() > fMinMCLabel) {
631  ++NMCmatches;
632  trkPMCfrac += mom;
633  }
634  ++Nmatches;
635  totalTrkP += mom;
636 
637  if (fCreateHisto) {
638  if (fHadCorr > 1) {
639  Double_t dphi = 0;
640  Double_t deta = 0;
641  GetEtaPhiDiff(track, cluster, dphi, deta);
642  Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
643  Double_t energyclus = cluster->GetNonLinCorrEnergy();
644  fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
645  }
646  }
647  }
648  }
649 
650  if (totalTrkP > 0) trkPMCfrac /= totalTrkP;
651 }
652 
653 //________________________________________________________________________
655 {
656  // Run the hadronic correction
657 
659 
660  // delete output
661  if (fOutClusters) fOutClusters->Delete();
662 
663  Int_t clusCount = 0;
664 
665  // loop over all clusters
666  clusters->ResetCurrentID();
667  AliVCluster *cluster = 0;
668  while ((cluster = clusters->GetNextAcceptCluster())) {
669 
670  Double_t energyclus = 0;
671  if (fCreateHisto) {
672  fHistEbefore->Fill(fCent, cluster->GetNonLinCorrEnergy());
673  fHistNclusvsCent->Fill(fCent);
674  }
675 
676  // apply correction / subtraction
677  // to subtract only the closest track set fHadCor to a %
678  // to subtract all tracks within the cut set fHadCor to %+1
679  if (fHadCorr > 1) {
680  energyclus = ApplyHadCorrAllTracks(clusters->GetCurrentID(), fHadCorr - 1);
681  }
682  else if (fHadCorr > 0) {
683  energyclus = ApplyHadCorrOneTrack(clusters->GetCurrentID(), fHadCorr);
684  }
685  else {
686  energyclus = cluster->GetNonLinCorrEnergy();
687  }
688 
689  if (energyclus < 0) energyclus = 0;
690 
691  cluster->SetHadCorrEnergy(energyclus);
692 
693  if (fCreateHisto) fHistEafter->Fill(fCent, energyclus);
694 
695  if (fOutClusters && energyclus > 0) { // create corrected cluster
696  AliVCluster *oc;
697  if (fEsdMode) {
698  AliESDCaloCluster *ec = dynamic_cast<AliESDCaloCluster*>(cluster);
699  if (!ec) continue;
700  oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*ec);
701  }
702  else {
703  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
704  if (!ac) continue;
705  oc = new ((*fOutClusters)[clusCount]) AliAODCaloCluster(*ac);
706  }
707  ++clusCount;
708  oc->SetE(energyclus);
709  oc->SetNonLinCorrEnergy(cluster->GetNonLinCorrEnergy()); //get the non linearity corrected energy of the clusters. Works only if the clustermaker was ran before
710  oc->SetHadCorrEnergy(energyclus); //same as the default energy field of this specific copy of the cluster container
711  }
712  }
713 
714  return kTRUE;
715 }
716 
717 //________________________________________________________________________
719 {
720  // Apply the hadronic correction with one track only.
721 
724 
725  AliVCluster* cluster = clusters->GetCluster(icluster);
726 
727  Double_t energyclus = cluster->GetNonLinCorrEnergy();
728 
729  AliVTrack* track = 0;
730 
731  if (cluster->GetNTracksMatched() > 0) {
732  if (fEsdMode) {
733  Int_t itrack = cluster->GetTrackMatchedIndex(0);
734  if (itrack >= 0) track = static_cast<AliVTrack*>(tracks->GetAcceptParticle(itrack));
735  }
736  else {
737  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
738  UInt_t rejectionReason = 0;
739  if (!tracks->AcceptParticle(track, rejectionReason)) track = 0;
740  }
741  }
742 
743  if (!track || track->P() < 1e-6) return energyclus;
744 
745  Double_t mom = track->P();
746 
747  Double_t dEtaMin = 1e9;
748  Double_t dPhiMin = 1e9;
749  GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
750 
751  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(dEtaMin, dPhiMin);
752 
753  // check if track also points to cluster
754  Int_t cid = track->GetEMCALcluster();
755  if (fDoTrackClus && (cid != icluster)) return energyclus;
756 
757  UInt_t mombin = GetMomBin(mom);
758  Int_t centbinch = fCentBin;
759  if (track->Charge() < 0) centbinch += fNcentBins;
760 
761  // plot some histograms if switched on
762  if (fCreateHisto) {
763  Int_t etabin = 0;
764  if(track->Eta() > 0) etabin = 1;
765 
766  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
767  fHistMatchEtaPhiAll->Fill(dEtaMin, dPhiMin);
768 
769  if (mom > 0) {
770  Double_t etadiff = 0;
771  Double_t phidiff = 0;
772  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
773  Double_t dRmin = TMath::Sqrt(etadiff*etadiff + phidiff*phidiff);
774  fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
775  fHistEoPCent->Fill(fCent, energyclus / mom);
776  fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
777  }
778  }
779 
780  // define eta/phi cuts
781  Double_t etaCut = 0.0;
782  Double_t phiCutlo = 0.0;
783  Double_t phiCuthi = 0.0;
784  if (fPhiMatch > 0) {
785  phiCutlo = -fPhiMatch;
786  phiCuthi = +fPhiMatch;
787  }
788  else {
789  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
790  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
791  }
792  if (fEtaMatch > 0) {
793  etaCut = fEtaMatch;
794  }
795  else {
796  etaCut = GetEtaSigma(mombin);
797  }
798 
799  // apply the correction if the track is in the eta/phi window
800  if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
801  energyclus -= hadCorr * mom;
802  }
803 
804  return energyclus;
805 }
806 
807 //________________________________________________________________________
809 {
810  // Apply the hadronic correction with all tracks.
811 
814 
815  AliVCluster* cluster = clusters->GetCluster(icluster);
816 
817  Double_t energyclus = cluster->GetNonLinCorrEnergy();
818  Double_t cNcells = cluster->GetNCells();
819 
820  Double_t totalTrkP = 0.0; // count total track momentum
821  Int_t Nmatches = 0; // count total number of matches
822 
823  Double_t trkPMCfrac = 0.0; // count total track momentum
824  Int_t NMCmatches = 0; // count total number of matches
825 
826  // do the loop over the matched tracks and get the number of matches and the total momentum
827  DoMatchedTracksLoop(icluster, totalTrkP, Nmatches, trkPMCfrac, NMCmatches);
828 
829  Double_t Esub = hadCorr * totalTrkP;
830 
831  if (Esub > energyclus) Esub = energyclus;
832 
833  // applying Peter's proposed algorithm
834  // never subtract the full energy of the cluster
835  Double_t clusEexcl = fEexclCell * cNcells;
836  if (energyclus < clusEexcl) clusEexcl = energyclus;
837  if ((energyclus - Esub) < clusEexcl) Esub = (energyclus - clusEexcl);
838 
839  // embedding
840  Double_t EsubMC = 0;
841  Double_t EsubBkg = 0;
842  Double_t EclusMC = 0;
843  Double_t EclusBkg = 0;
844  Double_t EclusCorr = 0;
845  Double_t EclusMCcorr = 0;
846  Double_t EclusBkgcorr = 0;
847  Double_t overSub = 0;
848  if (fIsEmbedded) {
849  EsubMC = hadCorr * totalTrkP * trkPMCfrac;
850  EsubBkg = hadCorr * totalTrkP - EsubMC;
851  EclusMC = energyclus * cluster->GetMCEnergyFraction();
852  EclusBkg = energyclus - EclusMC;
853 
854  if (energyclus > Esub)
855  EclusCorr = energyclus - Esub;
856 
857  if (EclusMC > EsubMC)
858  EclusMCcorr = EclusMC - EsubMC;
859 
860  if (EclusBkg > EsubBkg)
861  EclusBkgcorr = EclusBkg - EsubBkg;
862 
863  overSub = EclusMCcorr + EclusBkgcorr - EclusCorr;
864  }
865 
866  // plot some histograms if switched on
867  if (fCreateHisto) {
868  fHistNMatchCent->Fill(fCent, Nmatches);
869  fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
870 
871  if (Nmatches > 0) fHistNclusMatchvsCent->Fill(fCent);
872 
873  if (Nmatches < 3) {
874  fHistNCellsEnergy[fCentBin][Nmatches]->Fill(energyclus, cNcells);
875  }
876  else {
877  fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
878  }
879 
880  if (totalTrkP > 0) {
881  Double_t EoP = energyclus / totalTrkP;
882  fHistEoPCent->Fill(fCent, EoP);
883  fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
884 
885  fHistEsubPchRatAll[fCentBin]->Fill(totalTrkP, Esub / totalTrkP);
886 
887  if (Nmatches == 1) {
888  AliVTrack* track = 0;
889  if (fEsdMode) {
890  Int_t itrack = cluster->GetTrackMatchedIndex(0);
891  if (itrack >= 0) track = static_cast<AliVTrack*>(tracks->GetAcceptParticle(itrack));
892  }
893  else {
894  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
895  UInt_t rejectionReason = 0;
896  if (!tracks->AcceptParticle(track, rejectionReason)) track = 0;
897  }
898  if (track) {
899  Int_t centbinchm = fCentBin;
900  if (track->Charge() < 0) centbinchm += fNcentBins;
901  fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
902  fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
903  }
904  }
905 
906  if (fIsEmbedded) {
907  fHistOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
908 
909  if (cluster->GetMCEnergyFraction() > 0.95)
910  fHistOversubMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
911  else if (cluster->GetMCEnergyFraction() < 0.05)
912  fHistOversubNonMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
913 
914  if (trkPMCfrac < 0.05)
915  fHistNonEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
916  else if (trkPMCfrac > 0.95)
917  fHistEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
918  }
919  }
920  }
921 
922  if (fIsEmbedded && fDoExact) {
923  Esub -= overSub;
924  if (EclusBkgcorr + EclusMCcorr > 0) {
925  Double_t newfrac = EclusMCcorr / (EclusBkgcorr + EclusMCcorr);
926  cluster->SetMCEnergyFraction(newfrac);
927  }
928  }
929 
930  // apply the correction
931  energyclus -= Esub;
932 
933  return energyclus;
934 }
TH2 * fHistOversubMCClusters[5]
Over-subtracted energy / cluster energy with non-embedded track matches (embedded matches < 5%) ...
TH2 * fHistMatchdRvsEP[5]
n cells vs. cluster energy
TH2 * fHistNClusMatchCent
n matches vs. centraity
TString fOutCaloName
double Double_t
Definition: External.C:58
Definition: External.C:236
Double_t fHadCorr
Bool_t Run()
Run function. This is the core function of the analysis and contains the user code. Therefore users have to implement this function.
Double_t GetEtaSigma(Int_t pbin) const
Base task in the EMCAL framework.
Bool_t fLocalInitialized
whether or not the task has been already initialized
TH1 * fHistEafter
average energy of clusters before correction vs. centrality
TH2 * fHistEsubPchRat[10]
Esub vs. total momentum of matched tracks (only 1 match)
static void GetEtaPhiDiff(const AliVTrack *t, const AliVCluster *v, Double_t &phidiff, Double_t &etadiff)
Calculate and difference between a track (t) and a cluster (c).
Double_t fMinBinPt
min pt in histograms
Double_t GetPhiMean(Int_t pbin, Int_t centbin) const
Int_t fCentBin
!event centrality bin
Double_t ApplyHadCorrAllTracks(Int_t icluster, Double_t hadCorr)
void ExecOnce()
Perform steps needed to initialize the analysis.
UInt_t GetMomBin(Double_t pt) const
TH2 * fHistNMatchEnergy[5]
cluster energy vs. track momentum of matched pairs
Double_t fPhiMatch
TClonesArray * fOutClusters
ESD/AOD mode.
Container for particles within the EMCAL framework.
const TString & GetClassName() const
Bool_t fIsEmbedded
trigger, embedded signal
TH2 * fHistEsubPchRatAll[5]
Esub/momentum of matched tracks vs. total momentum of matched tracks (only 1 match) ...
TH2 * fHistMatchEtaPhiAll
deta vs. dphi of matched cluster-track pairs
AliParticleContainer * GetParticleContainer(Int_t i=0) const
Get particle container attached to this task.
TH2 * fHistMatchEtaPhi[10][9][2]
output cluster collection
TH2 * fHistEmbTrackMatchesOversub[5]
Esub/momentum of matched tracks vs. total momentum of matched tracks (all number of matches) ...
Double_t fEtaMatch
TH2 * fHistNMatchCent
E/P vs. centrality.
TH2 * fHistOversubNonMCClusters[5]
Over-subtracted energy / cluster energy (cluster MC energy fraction > 95%)
Double_t ApplyHadCorrOneTrack(Int_t icluster, Double_t hadCorr)
int Int_t
Definition: External.C:63
unsigned int UInt_t
Definition: External.C:33
Int_t GetCurrentID() const
Double_t fEexclCell
Double_t GetPhiSigma(Int_t pbin, Int_t centbin) const
void DoMatchedTracksLoop(Int_t icluster, Double_t &totalTrkP, Int_t &Nmatches, Double_t &trkPMCfrac, Int_t &NMCmatches)
Int_t fNcentBins
how many centrality bins
AliClusterContainer * GetClusterContainer(Int_t i=0) const
Get cluster container attached to this task.
virtual AliVParticle * GetAcceptParticle(Int_t i=-1) const
TH1 * fHistNclusMatchvsCent
n clusters vs. centrality
Double_t fCent
!event centrality
virtual Bool_t AcceptParticle(const AliVParticle *vp, UInt_t &rejectionReason) const
AliVCluster * GetCluster(Int_t i) const
void SetClassName(const char *clname)
TH2 * fHistMatchEtaPhiAllCl
deta vs. dphi of all cluster-track pairs (tr loop)
TH2 * fHistOversub[5]
Over-subtracted energy / cluster energy (cluster MC energy fraction < 5%)
Int_t fMinMCLabel
minimum MC label value for the tracks/clusters being considered MC particles
TH2 * fHistEoPCent
average energy of clusters after correction vs. centrality
TH2 * fHistNonEmbTrackMatchesOversub[5]
Over-subtracted energy / cluster energy with embedded track matches (non-embedded matches < 5%) ...
AliEmcalList * fOutput
!output list
Double_t fMaxBinPt
max pt in histograms
TH1 * fHistEsubPch[10]
n clusters macthed to some track (tracks allowed to match more than one cluster)
TH2 * fHistMatchEtaPhiAllTr
deta vs. dphi of matched cluster-track pairs
void UserCreateOutputObjects()
virtual ~AliHadCorrTask()
Bool_t fCreateHisto
whether or not create histograms
void SetMakeGeneralHistograms(Bool_t g)
Enable general histograms.
TH2 * fHistNCellsEnergy[5][4]
n matches vs. cluster energy
virtual void ExecOnce()
Perform steps needed to initialize the analysis.
void AddObjectToEvent(TObject *obj, Bool_t attempt=kFALSE)
Add object to event.
void UserCreateOutputObjects()
Main initialization function on the worker.
bool Bool_t
Definition: External.C:53
TH2 * fHistMatchEvsP[5]
deta vs. dphi of all cluster-track pairs (cl loop)
TH1 * fHistEbefore
n clusters matched to some track vs. centrality
TH1 * fHistNclusvsCent
matching distance vs. E/P
Container structure for EMCAL clusters.
AliVCluster * GetNextAcceptCluster()
void ResetCurrentID(Int_t i=-1)
Reset the iterator to a given index.
Int_t fNbins
no. of pt bins