]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG/EMCAL/AliHadCorrTask.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
AliAnalysisMuMu classes update
[u/mrichter/AliRoot.git] / PWG / EMCAL / AliHadCorrTask.cxx
CommitLineData
8417a2aa 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
27ClassImp(AliHadCorrTask)
28
29//________________________________________________________________________
30AliHadCorrTask::AliHadCorrTask() :
31 AliAnalysisTaskEmcal("AliHadCorrTask", kFALSE),
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//________________________________________________________________________
82AliHadCorrTask::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
116 fHistEmbTrackMatchesOversub[i] = 0;
117 fHistNonEmbTrackMatchesOversub[i] = 0;
118 fHistOversubMCClusters[i] = 0;
119 fHistOversubNonMCClusters[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
132 SetMakeGeneralHistograms(histo);
133
134 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
135}
136
137//________________________________________________________________________
138AliHadCorrTask::~AliHadCorrTask()
139{
140 // Destructor
141}
142
143//________________________________________________________________________
144UInt_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//________________________________________________________________________
172Double_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//________________________________________________________________________
181Double_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//________________________________________________________________________
279Double_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,
3030.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//________________________________________________________________________
377void AliHadCorrTask::UserCreateOutputObjects()
378{
379 // Create my user objects.
380
381 AliAnalysisTaskEmcal::UserCreateOutputObjects();
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);
391 fOutput->Add(fHistMatchEtaPhiAll);
392
393 fHistMatchEtaPhiAllTr = new TH2F("fHistMatchEtaPhiAllTr", "fHistMatchEtaPhiAllTr;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
394 fOutput->Add(fHistMatchEtaPhiAllTr);
395
396 fHistMatchEtaPhiAllCl = new TH2F("fHistMatchEtaPhiAllCl", "fHistMatchEtaPhiAllCl;#Delta#eta;#Delta#phi", fNbins, -0.1, 0.1, fNbins, -0.1, 0.1);
397 fOutput->Add(fHistMatchEtaPhiAllCl);
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
468 fOutput->Add(fHistNclusMatchvsCent);
469 fOutput->Add(fHistNclusvsCent);
470 fOutput->Add(fHistEbefore);
471 fOutput->Add(fHistEafter);
472 fOutput->Add(fHistEoPCent);
473 fOutput->Add(fHistNMatchCent);
474 fOutput->Add(fHistNClusMatchCent);
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}");
482 fOutput->Add(fHistEmbTrackMatchesOversub[icent]);
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}");
488 fOutput->Add(fHistNonEmbTrackMatchesOversub[icent]);
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//________________________________________________________________________
514void AliHadCorrTask::ExecOnce()
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
529 AliAnalysisTaskEmcal::ExecOnce();
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//________________________________________________________________________
553void AliHadCorrTask::DoTrackLoop()
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();
8417a2aa 585 if (!cluster) continue;
4dde7ef1 586 if (!cluster->IsEMCAL()) continue;
8417a2aa 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//________________________________________________________________________
600void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster,
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//________________________________________________________________________
683Bool_t AliHadCorrTask::Run()
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//________________________________________________________________________
749Double_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//________________________________________________________________________
830Double_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)
893 fHistNclusMatchvsCent->Fill(fCent);
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}