AliPhysics  8bb951a (8bb951a)
 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 //________________________________________________________________________
141 UInt_t AliHadCorrTask::GetMomBin(Double_t p) const
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 //________________________________________________________________________
169 Double_t AliHadCorrTask::GetEtaSigma(Int_t pbin) const
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 //________________________________________________________________________
178 Double_t AliHadCorrTask::GetPhiMean(Int_t pbin, Int_t centbin) const
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 //________________________________________________________________________
276 Double_t AliHadCorrTask::GetPhiSigma(Int_t pbin, Int_t centbin) const
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 (!fInitialized) 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  }
708  }
709 
710  return kTRUE;
711 }
712 
713 //________________________________________________________________________
714 Double_t AliHadCorrTask::ApplyHadCorrOneTrack(Int_t icluster, Double_t hadCorr)
715 {
716  // Apply the hadronic correction with one track only.
717 
720 
721  AliVCluster* cluster = clusters->GetCluster(icluster);
722 
723  Double_t energyclus = cluster->GetNonLinCorrEnergy();
724 
725  AliVTrack* track = 0;
726 
727  if (cluster->GetNTracksMatched() > 0) {
728  if (fEsdMode) {
729  Int_t itrack = cluster->GetTrackMatchedIndex(0);
730  if (itrack >= 0) track = static_cast<AliVTrack*>(tracks->GetAcceptParticle(itrack));
731  }
732  else {
733  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
734  UInt_t rejectionReason = 0;
735  if (!tracks->AcceptParticle(track, rejectionReason)) track = 0;
736  }
737  }
738 
739  if (!track || track->P() < 1e-6) return energyclus;
740 
741  Double_t mom = track->P();
742 
743  Double_t dEtaMin = 1e9;
744  Double_t dPhiMin = 1e9;
745  GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
746 
747  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(dEtaMin, dPhiMin);
748 
749  // check if track also points to cluster
750  Int_t cid = track->GetEMCALcluster();
751  if (fDoTrackClus && (cid != icluster)) return energyclus;
752 
753  UInt_t mombin = GetMomBin(mom);
754  Int_t centbinch = fCentBin;
755  if (track->Charge() < 0) centbinch += fNcentBins;
756 
757  // plot some histograms if switched on
758  if (fCreateHisto) {
759  Int_t etabin = 0;
760  if(track->Eta() > 0) etabin = 1;
761 
762  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
763  fHistMatchEtaPhiAll->Fill(dEtaMin, dPhiMin);
764 
765  if (mom > 0) {
766  Double_t etadiff = 0;
767  Double_t phidiff = 0;
768  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
769  Double_t dRmin = TMath::Sqrt(etadiff*etadiff + phidiff*phidiff);
770  fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
771  fHistEoPCent->Fill(fCent, energyclus / mom);
772  fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
773  }
774  }
775 
776  // define eta/phi cuts
777  Double_t etaCut = 0.0;
778  Double_t phiCutlo = 0.0;
779  Double_t phiCuthi = 0.0;
780  if (fPhiMatch > 0) {
781  phiCutlo = -fPhiMatch;
782  phiCuthi = +fPhiMatch;
783  }
784  else {
785  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
786  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
787  }
788  if (fEtaMatch > 0) {
789  etaCut = fEtaMatch;
790  }
791  else {
792  etaCut = GetEtaSigma(mombin);
793  }
794 
795  // apply the correction if the track is in the eta/phi window
796  if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
797  energyclus -= hadCorr * mom;
798  }
799 
800  return energyclus;
801 }
802 
803 //________________________________________________________________________
804 Double_t AliHadCorrTask::ApplyHadCorrAllTracks(Int_t icluster, Double_t hadCorr)
805 {
806  // Apply the hadronic correction with all tracks.
807 
810 
811  AliVCluster* cluster = clusters->GetCluster(icluster);
812 
813  Double_t energyclus = cluster->GetNonLinCorrEnergy();
814  Double_t cNcells = cluster->GetNCells();
815 
816  Double_t totalTrkP = 0.0; // count total track momentum
817  Int_t Nmatches = 0; // count total number of matches
818 
819  Double_t trkPMCfrac = 0.0; // count total track momentum
820  Int_t NMCmatches = 0; // count total number of matches
821 
822  // do the loop over the matched tracks and get the number of matches and the total momentum
823  DoMatchedTracksLoop(icluster, totalTrkP, Nmatches, trkPMCfrac, NMCmatches);
824 
825  Double_t Esub = hadCorr * totalTrkP;
826 
827  if (Esub > energyclus) Esub = energyclus;
828 
829  // applying Peter's proposed algorithm
830  // never subtract the full energy of the cluster
831  Double_t clusEexcl = fEexclCell * cNcells;
832  if (energyclus < clusEexcl) clusEexcl = energyclus;
833  if ((energyclus - Esub) < clusEexcl) Esub = (energyclus - clusEexcl);
834 
835  // embedding
836  Double_t EsubMC = 0;
837  Double_t EsubBkg = 0;
838  Double_t EclusMC = 0;
839  Double_t EclusBkg = 0;
840  Double_t EclusCorr = 0;
841  Double_t EclusMCcorr = 0;
842  Double_t EclusBkgcorr = 0;
843  Double_t overSub = 0;
844  if (fIsEmbedded) {
845  EsubMC = hadCorr * totalTrkP * trkPMCfrac;
846  EsubBkg = hadCorr * totalTrkP - EsubMC;
847  EclusMC = energyclus * cluster->GetMCEnergyFraction();
848  EclusBkg = energyclus - EclusMC;
849 
850  if (energyclus > Esub)
851  EclusCorr = energyclus - Esub;
852 
853  if (EclusMC > EsubMC)
854  EclusMCcorr = EclusMC - EsubMC;
855 
856  if (EclusBkg > EsubBkg)
857  EclusBkgcorr = EclusBkg - EsubBkg;
858 
859  overSub = EclusMCcorr + EclusBkgcorr - EclusCorr;
860  }
861 
862  // plot some histograms if switched on
863  if (fCreateHisto) {
864  fHistNMatchCent->Fill(fCent, Nmatches);
865  fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
866 
867  if (Nmatches > 0) fHistNclusMatchvsCent->Fill(fCent);
868 
869  if (Nmatches < 3) {
870  fHistNCellsEnergy[fCentBin][Nmatches]->Fill(energyclus, cNcells);
871  }
872  else {
873  fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
874  }
875 
876  if (totalTrkP > 0) {
877  Double_t EoP = energyclus / totalTrkP;
878  fHistEoPCent->Fill(fCent, EoP);
879  fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
880 
881  fHistEsubPchRatAll[fCentBin]->Fill(totalTrkP, Esub / totalTrkP);
882 
883  if (Nmatches == 1) {
884  AliVTrack* track = 0;
885  if (fEsdMode) {
886  Int_t itrack = cluster->GetTrackMatchedIndex(0);
887  if (itrack >= 0) track = static_cast<AliVTrack*>(tracks->GetAcceptParticle(itrack));
888  }
889  else {
890  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
891  UInt_t rejectionReason = 0;
892  if (!tracks->AcceptParticle(track, rejectionReason)) track = 0;
893  }
894  if (track) {
895  Int_t centbinchm = fCentBin;
896  if (track->Charge() < 0) centbinchm += fNcentBins;
897  fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
898  fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
899  }
900  }
901 
902  if (fIsEmbedded) {
903  fHistOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
904 
905  if (cluster->GetMCEnergyFraction() > 0.95)
906  fHistOversubMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
907  else if (cluster->GetMCEnergyFraction() < 0.05)
908  fHistOversubNonMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
909 
910  if (trkPMCfrac < 0.05)
911  fHistNonEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
912  else if (trkPMCfrac > 0.95)
913  fHistEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
914  }
915  }
916  }
917 
918  if (fIsEmbedded && fDoExact) {
919  Esub -= overSub;
920  if (EclusBkgcorr + EclusMCcorr > 0) {
921  Double_t newfrac = EclusMCcorr / (EclusBkgcorr + EclusMCcorr);
922  cluster->SetMCEnergyFraction(newfrac);
923  }
924  }
925 
926  // apply the correction
927  energyclus -= Esub;
928 
929  return energyclus;
930 }
TH2 * fHistNonEmbTrackMatchesOversub[4]
Over-subtracted energy / cluster energy with embedded track matches (non-embedded matches < 5%) ...
TH2 * fHistNClusMatchCent
n matches vs. centraity
TString fOutCaloName
TH2 * fHistOversubNonMCClusters[4]
Over-subtracted energy / cluster energy (cluster MC energy fraction > 95%)
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.
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
TList * fOutput
!output list
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 * 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_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
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
void SetClassName(const char *clname)
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
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)
TH2 * fHistEsubPchRat[8]
Esub vs. total momentum of matched tracks (only 1 match)
TH1 * fHistEbefore
n clusters matched to some track vs. centrality
Bool_t fInitialized
whether or not the task has been already initialized
TH1 * fHistNclusvsCent
matching distance vs. E/P
Container structure for EMCAL clusters.
AliVCluster * GetNextAcceptCluster()
void ResetCurrentID(Int_t i=-1)
Int_t fNbins
no. of pt bins