AliPhysics  a0db429 (a0db429)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
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  if (!tracks->AcceptParticle(track)) track = 0;
581  }
582 
583  if (!track) continue;
584 
585  Double_t etadiff = 999;
586  Double_t phidiff = 999;
587  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
588  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(etadiff, phidiff);
589 
590  // check if track also points to cluster
591  if (fDoTrackClus && (track->GetEMCALcluster() != icluster)) continue;
592 
593  Double_t mom = track->P();
594  UInt_t mombin = GetMomBin(mom);
595  Int_t centbinch = fCentBin;
596 
597  if (track->Charge() < 0) centbinch += fNcentBins;
598 
599  if (fCreateHisto) {
600  Int_t etabin = 0;
601  if (track->Eta() > 0) etabin=1;
602  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
603  fHistMatchEtaPhiAll->Fill(etadiff, phidiff);
604  }
605 
606  Double_t etaCut = 0.0;
607  Double_t phiCutlo = 0.0;
608  Double_t phiCuthi = 0.0;
609 
610  if (fPhiMatch > 0) {
611  phiCutlo = -fPhiMatch;
612  phiCuthi = +fPhiMatch;
613  }
614  else {
615  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
616  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
617  }
618 
619  if (fEtaMatch > 0) {
620  etaCut = fEtaMatch;
621  }
622  else {
623  etaCut = GetEtaSigma(mombin);
624  }
625 
626  if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
627  if (track->GetLabel() > fMinMCLabel) {
628  ++NMCmatches;
629  trkPMCfrac += mom;
630  }
631  ++Nmatches;
632  totalTrkP += mom;
633 
634  if (fCreateHisto) {
635  if (fHadCorr > 1) {
636  Double_t dphi = 0;
637  Double_t deta = 0;
638  GetEtaPhiDiff(track, cluster, dphi, deta);
639  Double_t dR = TMath::Sqrt(dphi*dphi + deta*deta);
640  Double_t energyclus = cluster->GetNonLinCorrEnergy();
641  fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
642  }
643  }
644  }
645  }
646 
647  if (totalTrkP > 0) trkPMCfrac /= totalTrkP;
648 }
649 
650 //________________________________________________________________________
652 {
653  // Run the hadronic correction
654 
656 
657  // delete output
658  if (fOutClusters) fOutClusters->Delete();
659 
660  Int_t clusCount = 0;
661 
662  // loop over all clusters
663  clusters->ResetCurrentID();
664  AliVCluster *cluster = 0;
665  while ((cluster = clusters->GetNextAcceptCluster())) {
666 
667  Double_t energyclus = 0;
668  if (fCreateHisto) {
669  fHistEbefore->Fill(fCent, cluster->GetNonLinCorrEnergy());
670  fHistNclusvsCent->Fill(fCent);
671  }
672 
673  // apply correction / subtraction
674  // to subtract only the closest track set fHadCor to a %
675  // to subtract all tracks within the cut set fHadCor to %+1
676  if (fHadCorr > 1) {
677  energyclus = ApplyHadCorrAllTracks(clusters->GetCurrentID(), fHadCorr - 1);
678  }
679  else if (fHadCorr > 0) {
680  energyclus = ApplyHadCorrOneTrack(clusters->GetCurrentID(), fHadCorr);
681  }
682  else {
683  energyclus = cluster->GetNonLinCorrEnergy();
684  }
685 
686  if (energyclus < 0) energyclus = 0;
687 
688  cluster->SetHadCorrEnergy(energyclus);
689 
690  if (fCreateHisto) fHistEafter->Fill(fCent, energyclus);
691 
692  if (fOutClusters && energyclus > 0) { // create corrected cluster
693  AliVCluster *oc;
694  if (fEsdMode) {
695  AliESDCaloCluster *ec = dynamic_cast<AliESDCaloCluster*>(cluster);
696  if (!ec) continue;
697  oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*ec);
698  }
699  else {
700  AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
701  if (!ac) continue;
702  oc = new ((*fOutClusters)[clusCount]) AliAODCaloCluster(*ac);
703  }
704  ++clusCount;
705  oc->SetE(energyclus);
706  }
707  }
708 
709  return kTRUE;
710 }
711 
712 //________________________________________________________________________
713 Double_t AliHadCorrTask::ApplyHadCorrOneTrack(Int_t icluster, Double_t hadCorr)
714 {
715  // Apply the hadronic correction with one track only.
716 
719 
720  AliVCluster* cluster = clusters->GetCluster(icluster);
721 
722  Double_t energyclus = cluster->GetNonLinCorrEnergy();
723 
724  AliVTrack* track = 0;
725 
726  if (cluster->GetNTracksMatched() > 0) {
727  if (fEsdMode) {
728  Int_t itrack = cluster->GetTrackMatchedIndex(0);
729  if (itrack >= 0) track = static_cast<AliVTrack*>(tracks->GetAcceptParticle(itrack));
730  }
731  else {
732  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
733  if (!tracks->AcceptParticle(track)) track = 0;
734  }
735  }
736 
737  if (!track || track->P() < 1e-6) return energyclus;
738 
739  Double_t mom = track->P();
740 
741  Double_t dEtaMin = 1e9;
742  Double_t dPhiMin = 1e9;
743  GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
744 
745  if (fCreateHisto) fHistMatchEtaPhiAllCl->Fill(dEtaMin, dPhiMin);
746 
747  // check if track also points to cluster
748  Int_t cid = track->GetEMCALcluster();
749  if (fDoTrackClus && (cid != icluster)) return energyclus;
750 
751  UInt_t mombin = GetMomBin(mom);
752  Int_t centbinch = fCentBin;
753  if (track->Charge() < 0) centbinch += fNcentBins;
754 
755  // plot some histograms if switched on
756  if (fCreateHisto) {
757  Int_t etabin = 0;
758  if(track->Eta() > 0) etabin = 1;
759 
760  fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
761  fHistMatchEtaPhiAll->Fill(dEtaMin, dPhiMin);
762 
763  if (mom > 0) {
764  Double_t etadiff = 0;
765  Double_t phidiff = 0;
766  GetEtaPhiDiff(track, cluster, phidiff, etadiff);
767  Double_t dRmin = TMath::Sqrt(etadiff*etadiff + phidiff*phidiff);
768  fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
769  fHistEoPCent->Fill(fCent, energyclus / mom);
770  fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
771  }
772  }
773 
774  // define eta/phi cuts
775  Double_t etaCut = 0.0;
776  Double_t phiCutlo = 0.0;
777  Double_t phiCuthi = 0.0;
778  if (fPhiMatch > 0) {
779  phiCutlo = -fPhiMatch;
780  phiCuthi = +fPhiMatch;
781  }
782  else {
783  phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
784  phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
785  }
786  if (fEtaMatch > 0) {
787  etaCut = fEtaMatch;
788  }
789  else {
790  etaCut = GetEtaSigma(mombin);
791  }
792 
793  // apply the correction if the track is in the eta/phi window
794  if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
795  energyclus -= hadCorr * mom;
796  }
797 
798  return energyclus;
799 }
800 
801 //________________________________________________________________________
802 Double_t AliHadCorrTask::ApplyHadCorrAllTracks(Int_t icluster, Double_t hadCorr)
803 {
804  // Apply the hadronic correction with all tracks.
805 
808 
809  AliVCluster* cluster = clusters->GetCluster(icluster);
810 
811  Double_t energyclus = cluster->GetNonLinCorrEnergy();
812  Double_t cNcells = cluster->GetNCells();
813 
814  Double_t totalTrkP = 0.0; // count total track momentum
815  Int_t Nmatches = 0; // count total number of matches
816 
817  Double_t trkPMCfrac = 0.0; // count total track momentum
818  Int_t NMCmatches = 0; // count total number of matches
819 
820  // do the loop over the matched tracks and get the number of matches and the total momentum
821  DoMatchedTracksLoop(icluster, totalTrkP, Nmatches, trkPMCfrac, NMCmatches);
822 
823  Double_t Esub = hadCorr * totalTrkP;
824 
825  if (Esub > energyclus) Esub = energyclus;
826 
827  // applying Peter's proposed algorithm
828  // never subtract the full energy of the cluster
829  Double_t clusEexcl = fEexclCell * cNcells;
830  if (energyclus < clusEexcl) clusEexcl = energyclus;
831  if ((energyclus - Esub) < clusEexcl) Esub = (energyclus - clusEexcl);
832 
833  // embedding
834  Double_t EsubMC = 0;
835  Double_t EsubBkg = 0;
836  Double_t EclusMC = 0;
837  Double_t EclusBkg = 0;
838  Double_t EclusCorr = 0;
839  Double_t EclusMCcorr = 0;
840  Double_t EclusBkgcorr = 0;
841  Double_t overSub = 0;
842  if (fIsEmbedded) {
843  EsubMC = hadCorr * totalTrkP * trkPMCfrac;
844  EsubBkg = hadCorr * totalTrkP - EsubMC;
845  EclusMC = energyclus * cluster->GetMCEnergyFraction();
846  EclusBkg = energyclus - EclusMC;
847 
848  if (energyclus > Esub)
849  EclusCorr = energyclus - Esub;
850 
851  if (EclusMC > EsubMC)
852  EclusMCcorr = EclusMC - EsubMC;
853 
854  if (EclusBkg > EsubBkg)
855  EclusBkgcorr = EclusBkg - EsubBkg;
856 
857  overSub = EclusMCcorr + EclusBkgcorr - EclusCorr;
858  }
859 
860  // plot some histograms if switched on
861  if (fCreateHisto) {
862  fHistNMatchCent->Fill(fCent, Nmatches);
863  fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
864 
865  if (Nmatches > 0) fHistNclusMatchvsCent->Fill(fCent);
866 
867  if (Nmatches < 3) {
868  fHistNCellsEnergy[fCentBin][Nmatches]->Fill(energyclus, cNcells);
869  }
870  else {
871  fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
872  }
873 
874  if (totalTrkP > 0) {
875  Double_t EoP = energyclus / totalTrkP;
876  fHistEoPCent->Fill(fCent, EoP);
877  fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
878 
879  fHistEsubPchRatAll[fCentBin]->Fill(totalTrkP, Esub / totalTrkP);
880 
881  if (Nmatches == 1) {
882  AliVTrack* track = 0;
883  if (fEsdMode) {
884  Int_t itrack = cluster->GetTrackMatchedIndex(0);
885  if (itrack >= 0) track = static_cast<AliVTrack*>(tracks->GetAcceptParticle(itrack));
886  }
887  else {
888  track = static_cast<AliVTrack*>(cluster->GetTrackMatched(0));
889  if (!tracks->AcceptParticle(track)) track = 0;
890  }
891  if (track) {
892  Int_t centbinchm = fCentBin;
893  if (track->Charge() < 0) centbinchm += fNcentBins;
894  fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
895  fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
896  }
897  }
898 
899  if (fIsEmbedded) {
900  fHistOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
901 
902  if (cluster->GetMCEnergyFraction() > 0.95)
903  fHistOversubMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
904  else if (cluster->GetMCEnergyFraction() < 0.05)
905  fHistOversubNonMCClusters[fCentBin]->Fill(energyclus, overSub / energyclus);
906 
907  if (trkPMCfrac < 0.05)
908  fHistNonEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
909  else if (trkPMCfrac > 0.95)
910  fHistEmbTrackMatchesOversub[fCentBin]->Fill(energyclus, overSub / energyclus);
911  }
912  }
913  }
914 
915  if (fIsEmbedded && fDoExact) {
916  Esub -= overSub;
917  if (EclusBkgcorr + EclusMCcorr > 0) {
918  Double_t newfrac = EclusMCcorr / (EclusBkgcorr + EclusMCcorr);
919  cluster->SetMCEnergyFraction(newfrac);
920  }
921  }
922 
923  // apply the correction
924  energyclus -= Esub;
925 
926  return energyclus;
927 }
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
void SetClassName(const char *clname)
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 GetPhiMean(Int_t pbin, Int_t centbin) const
Int_t fCentBin
event centrality
AliVCluster * GetNextAcceptCluster(Int_t i=-1)
TH2 * fHistNCellsEnergy[4][4]
n matches vs. cluster energy
TH2 * fHistOversub[4]
Over-subtracted energy / cluster energy (cluster MC energy fraction < 5%)
Bool_t AcceptParticle(AliVParticle *vp)
Double_t ApplyHadCorrAllTracks(Int_t icluster, Double_t hadCorr)
UInt_t GetMomBin(Double_t pt) const
TList * fOutput
x-section from pythia header
Double_t fPhiMatch
TClonesArray * fOutClusters
ESD/AOD mode.
const TString & GetClassName() const
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
AliVParticle * GetAcceptParticle(Int_t i)
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)
AliClusterContainer * GetClusterContainer(Int_t i=0) const
TH1 * fHistNclusMatchvsCent
n clusters vs. centrality
Double_t fCent
trigger patch info array
TH2 * fHistMatchdRvsEP[4]
n cells vs. cluster energy
AliVCluster * GetCluster(Int_t i) const
TH2 * fHistMatchEtaPhiAllCl
deta vs. dphi of all cluster-track pairs (tr loop)
void SetClassName(const char *clname)
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
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()
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
TH1 * fHistNclusvsCent
matching distance vs. E/P
void ResetCurrentID(Int_t i=0)