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