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