protections against failures in deleting event content
[u/mrichter/AliRoot.git] / PWGJE / EMCALJetTasks / AliHadCorrTask.cxx
CommitLineData
d6f334b5 1// $Id$
2//
3// Hadronic correction task.
4//
1f6c8e39 5// Author: R.Reed, C.Loizides, S.Aiola
d6f334b5 6
7#include <TChain.h>
0b777a09 8#include <TClonesArray.h>
35789a2d 9#include <TH1F.h>
10#include <TH2F.h>
d6f334b5 11#include <TList.h>
bcf52e4a 12
3a385c49 13#include "AliAnalysisManager.h"
d6f334b5 14#include "AliAODEvent.h"
1f6c8e39 15#include "AliAODCaloCluster.h"
d6f334b5 16#include "AliESDCaloCluster.h"
1f6c8e39 17#include "AliVTrack.h"
0b777a09 18#include "AliPicoTrack.h"
35789a2d 19#include "AliVEventHandler.h"
1f6c8e39 20#include "AliEmcalParticle.h"
bcf52e4a 21
35789a2d 22#include "AliHadCorrTask.h"
0b777a09 23
24ClassImp(AliHadCorrTask)
25
26//________________________________________________________________________
d6f334b5 27AliHadCorrTask::AliHadCorrTask() :
1f6c8e39 28 AliAnalysisTaskEmcal("AliHadCorrTask", kFALSE),
d6f334b5 29 fOutCaloName(),
a3f43e7e 30 fPhiMatch(0.05),
31 fEtaMatch(0.025),
32bf39af 32 fDoTrackClus(0),
f22cb876 33 fHadCorr(0),
1f6c8e39 34 fEexclCell(0),
f22cb876 35 fOutClusters(0),
0b777a09 36 fHistNclusvsCent(0),
37 fHistNclusMatchvsCent(0),
f22cb876 38 fHistEbefore(0),
4711c521 39 fHistEafter(0),
40 fHistEoPCent(0),
32bf39af 41 fHistNMatchCent(0),
a89e2005 42 fHistCentrality(0)
0b777a09 43{
d6f334b5 44 // Default constructor.
0b777a09 45
a3f43e7e 46 for(Int_t i=0; i<8; i++) {
1f6c8e39 47 fHistEsubPch[i] = 0;
4f833c06 48 fHistEsubPchRat[i] = 0;
1f6c8e39 49 for(Int_t j=0; j<4; j++) {
50 fHistNCellsEnergy[i][j] = 0;
51 }
a3f43e7e 52 if (i<4) {
1f6c8e39 53 fHistMatchEvsP[i] = 0;
54 fHistMatchdRvsEP[i] = 0;
55 fHistNMatchEnergy[i] = 0;
a3f43e7e 56 }
fbc9afc4 57 for(Int_t j=0; j<9; j++) {
1f6c8e39 58 for(Int_t k=0; k<2; k++) {
59 fHistMatchEtaPhi[i][j][k] = 0;
60 }
a3f43e7e 61 }
62 }
0b777a09 63}
64
d6f334b5 65//________________________________________________________________________
5d845887 66AliHadCorrTask::AliHadCorrTask(const char *name, Bool_t histo) :
1f6c8e39 67 AliAnalysisTaskEmcal(name, histo),
5d845887 68 fOutCaloName("CaloClustersCorr"),
69 fPhiMatch(0.05),
70 fEtaMatch(0.025),
a89e2005 71 fDoTrackClus(1),
5d845887 72 fHadCorr(0),
1f6c8e39 73 fEexclCell(0),
5d845887 74 fOutClusters(0),
5d845887 75 fHistNclusvsCent(0),
76 fHistNclusMatchvsCent(0),
77 fHistEbefore(0),
78 fHistEafter(0),
79 fHistEoPCent(0),
80 fHistNMatchCent(0),
a89e2005 81 fHistCentrality(0)
5d845887 82{
83 // Standard constructor.
84
85 for(Int_t i=0; i<8; i++) {
86 if (i<4) {
87 fHistMatchEvsP[i] = 0;
88 fHistMatchdRvsEP[i] = 0;
89 for(Int_t j=0; j<3; j++) {
1f6c8e39 90 fHistNCellsEnergy[i][j] = 0;
5d845887 91 }
92 }
4f833c06 93 fHistEsubPch[i] = 0;
1f6c8e39 94 fHistEsubPchRat[i] = 0;
5d845887 95 for(Int_t j=0; j<9; j++) {
1f6c8e39 96 for(Int_t k=0; k<2; k++) {
97 fHistMatchEtaPhi[i][j][k] = 0;
98 }
5d845887 99 }
100 }
101
102 fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
d6f334b5 103}
0b777a09 104
105//________________________________________________________________________
106AliHadCorrTask::~AliHadCorrTask()
107{
108 // Destructor
109}
110
111//________________________________________________________________________
d6f334b5 112Int_t AliHadCorrTask::GetMomBin(Double_t p) const
113{
114 // Get momenum bin.
115
116 Int_t pbin=-1;
4f833c06 117 if (p<0.5)
d6f334b5 118 pbin=0;
fbc9afc4 119 else if (p>=0.5 && p<1.0)
d6f334b5 120 pbin=1;
fbc9afc4 121 else if (p>=1.0 && p<1.5)
d6f334b5 122 pbin=2;
fbc9afc4 123 else if (p>=1.5 && p<2.)
d6f334b5 124 pbin=3;
fbc9afc4 125 else if (p>=2. && p<3.)
d6f334b5 126 pbin=4;
fbc9afc4 127 else if (p>=3. && p<4.)
128 pbin=5;
129 else if (p>=4. && p<5.)
130 pbin=6;
131 else if (p>=5. && p<8.)
132 pbin=7;
133 else if (p>=8.)
134 pbin=8;
d6f334b5 135
136 return pbin;
137}
138
139//________________________________________________________________________
f6e8a928 140Double_t AliHadCorrTask::GetEtaSigma(Int_t pbin) const
141{
4f833c06 142 // Get sigma in eta.
143
1f6c8e39 144 Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
f6e8a928 145 return 2.0*EtaSigma[pbin];
146}
1f6c8e39 147
f6e8a928 148//________________________________________________________________________
149Double_t AliHadCorrTask::GetPhiMean(Int_t pbin, Int_t centbin) const
150{
4f833c06 151 // Get phi mean.
152
153 if (centbin==0) {
f6e8a928 154 Double_t PhiMean[9]={0.036,
155 0.021,
156 0.0121,
157 0.0084,
158 0.0060,
159 0.0041,
160 0.0031,
161 0.0022,
162 0.001};
163 return PhiMean[pbin];
4f833c06 164 } else if (centbin==1) {
f6e8a928 165 Double_t PhiMean[9]={0.036,
166 0.021,
167 0.0121,
168 0.0084,
169 0.0060,
170 0.0041,
171 0.0031,
172 0.0022,
173 0.001};
174 return PhiMean[pbin];
4f833c06 175 } else if (centbin==2) {
f6e8a928 176 Double_t PhiMean[9]={0.036,
177 0.021,
178 0.0121,
179 0.0084,
180 0.0060,
181 0.0041,
182 0.0031,
183 0.0022,
184 0.001};
185 return PhiMean[pbin];
4f833c06 186 } else if (centbin==3) {
f6e8a928 187 Double_t PhiMean[9]={0.036,
188 0.021,
189 0.0121,
190 0.0084,
191 0.0060,
192 0.0041,
193 0.0031,
194 0.0022,
195 0.001};
196 return PhiMean[pbin];
4f833c06 197 } else if (centbin==4) {
f6e8a928 198 Double_t PhiMean[9]={0.036,
199 0.021,
200 0.0127,
201 0.0089,
202 0.0068,
203 0.0049,
204 0.0038,
205 0.0028,
206 0.0018};
207 return PhiMean[pbin]*(-1.);
4f833c06 208 } else if (centbin==5) {
f6e8a928 209 Double_t PhiMean[9]={0.036,
210 0.021,
211 0.0127,
212 0.0089,
213 0.0068,
214 0.0048,
215 0.0038,
216 0.0028,
217 0.0018};
218 return PhiMean[pbin]*(-1.);
4f833c06 219 } else if (centbin==6) {
f6e8a928 220 Double_t PhiMean[9]={0.036,
221 0.021,
222 0.0127,
223 0.0089,
224 0.0068,
225 0.0045,
226 0.0035,
227 0.0028,
228 0.0018};
229 return PhiMean[pbin]*(-1.);
4f833c06 230 } else if (centbin==7) {
f6e8a928 231 Double_t PhiMean[9]={0.036,
232 0.021,
233 0.0127,
234 0.0089,
235 0.0068,
236 0.0043,
237 0.0035,
238 0.0028,
239 0.0018};
240 return PhiMean[pbin]*(-1.);
241 }
242
243 return 0;
f6e8a928 244}
1f6c8e39 245
f6e8a928 246//________________________________________________________________________
247Double_t AliHadCorrTask::GetPhiSigma(Int_t pbin, Int_t centbin) const
248{
4f833c06 249 // Get phi sigma.
250
251 if (centbin==0) {
f6e8a928 252 Double_t PhiSigma[9]={0.0221,
253 0.0128,
254 0.0074,
255 0.0064,
256 0.0059,
257 0.0055,
258 0.0052,
259 0.0049,
260 0.0045};
261 return 2.*PhiSigma[pbin];
4f833c06 262 } else if (centbin==1) {
f6e8a928 263 Double_t PhiSigma[9]={0.0217,
264 0.0120,
265 0.0076,
266 0.0066,
267 0.0062,
268 0.0058,
269 0.0054,
270 0.0054,
2710.0045};
272 return 2.*PhiSigma[pbin];
4f833c06 273 } else if (centbin==2) {
f6e8a928 274 Double_t PhiSigma[9]={0.0211,
275 0.0124,
276 0.0080,
277 0.0070,
278 0.0067,
279 0.0061,
280 0.0059,
281 0.0054,
282 0.0047};
283 return 2.*PhiSigma[pbin];
4f833c06 284 } else if (centbin==3) {
f6e8a928 285 Double_t PhiSigma[9]={0.0215,
286 0.0124,
287 0.0082,
288 0.0073,
289 0.0069,
290 0.0064,
291 0.0060,
292 0.0055,
293 0.0047};
294 return 2.*PhiSigma[pbin];
4f833c06 295 } else if (centbin==4) {
f6e8a928 296 Double_t PhiSigma[9]={0.0199,
297 0.0108,
298 0.0072,
299 0.0071,
300 0.0060,
301 0.0055,
302 0.0052,
303 0.0049,
304 0.0045};
305 return 2.*PhiSigma[pbin];
4f833c06 306 } else if (centbin==5) {
f6e8a928 307 Double_t PhiSigma[9]={0.0200,
308 0.0110,
309 0.0074,
310 0.0071,
311 0.0064,
312 0.0059,
313 0.0055,
314 0.0052,
315 0.0045};
316 return 2.*PhiSigma[pbin];
4f833c06 317 } else if (centbin==6) {
f6e8a928 318 Double_t PhiSigma[9]={0.0202,
319 0.0113,
320 0.0077,
321 0.0071,
322 0.0069,
323 0.0064,
324 0.0060,
325 0.0055,
326 0.0050};
327 return 2.*PhiSigma[pbin];
4f833c06 328 } else if (centbin==7) {
f6e8a928 329 Double_t PhiSigma[9]={0.0205,
330 0.0113,
331 0.0080,
332 0.0074,
333 0.0078,
334 0.0067,
335 0.0062,
336 0.0055,
337 0.0050};
338 return 2.*PhiSigma[pbin];
339 }
340
341 return 0;
f6e8a928 342}
343
344//________________________________________________________________________
0b777a09 345void AliHadCorrTask::UserCreateOutputObjects()
346{
d6f334b5 347 // Create my user objects.
0b777a09 348
35789a2d 349 AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
35789a2d 350 if (!handler) {
181a0a09 351 AliError(Form("%s: Input handler not available!", GetName()));
35789a2d 352 return;
353 }
354
355 if (handler->InheritsFrom("AliESDInputHandler")) {
356 fOutClusters = new TClonesArray("AliESDCaloCluster");
d6f334b5 357 } else if (handler->InheritsFrom("AliAODInputHandler")) {
35789a2d 358 fOutClusters = new TClonesArray("AliAODCaloCluster");
d6f334b5 359 } else {
181a0a09 360 AliError(Form("%s: Input handler not recognized!", GetName()));
35789a2d 361 return;
362 }
0b777a09 363 fOutClusters->SetName(fOutCaloName);
bcf52e4a 364
5d845887 365 if (!fCreateHisto)
366 return;
367
0b777a09 368 OpenFile(1);
1f6c8e39 369 fOutput = new TList();
370 fOutput->SetOwner();
0b777a09 371
1f6c8e39 372 TString name;
bcf52e4a 373
a3f43e7e 374 for(Int_t icent=0; icent<8; ++icent) {
fbc9afc4 375 for(Int_t ipt=0; ipt<9; ++ipt) {
1f6c8e39 376 for(Int_t ieta=0; ieta<2; ++ieta) {
1f6c8e39 377 name = Form("fHistMatchEtaPhi_%i_%i_%i",icent,ipt,ieta);
378 fHistMatchEtaPhi[icent][ipt][ieta] = new TH2F(name, name, 200, -0.1, 0.1, 400, -0.2, 0.2);
379 fOutput->Add(fHistMatchEtaPhi[icent][ipt][ieta]);
380 }
4711c521 381 }
382
1f6c8e39 383 for(Int_t itrk=0; itrk<4; ++itrk) {
1f6c8e39 384 name = Form("fHistNCellsEnergy_%i_%i",icent,itrk);
385 fHistNCellsEnergy[icent][itrk] = new TH2F(name, name, 1000, 0, 100, 101, -0.5, 100.5);
386 fOutput->Add(fHistNCellsEnergy[icent][itrk]);
387 }
388
389 name = Form("fHistEsubPch_%i",icent);
390 fHistEsubPch[icent]=new TH1F(name, name, 400, 0., 100.);
391 fOutput->Add(fHistEsubPch[icent]);
392
393 name = Form("fHistEsubPchRat_%i",icent);
394 fHistEsubPchRat[icent]=new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
395 fOutput->Add(fHistEsubPchRat[icent]);
1f6c8e39 396
4f833c06 397 if (icent<4) {
1f6c8e39 398 name = Form("fHistMatchEvsP_%i",icent);
a3f43e7e 399 fHistMatchEvsP[icent] = new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
1f6c8e39 400 fOutput->Add(fHistMatchEvsP[icent]);
0b777a09 401
a3f43e7e 402 name = Form("fHistMatchdRvsEP_%i",icent);
403 fHistMatchdRvsEP[icent] = new TH2F(name, name, 1000, 0., 1., 1000, 0., 10.);
1f6c8e39 404 fOutput->Add(fHistMatchdRvsEP[icent]);
405
406 name = Form("fHistNMatchEnergy_%i",icent);
407 fHistNMatchEnergy[icent] = new TH2F(name, name, 1000, 0, 100, 101, -0.5, 100.5);
408 fOutput->Add(fHistNMatchEnergy[icent]);
a3f43e7e 409 }
d6f334b5 410 }
0b777a09 411
32bf39af 412 fHistCentrality = new TH1F("fHistCentrality", "Centrality", 100, 0, 100);
d6f334b5 413 fHistNclusvsCent = new TH1F("Nclusvscent", "NclusVsCent", 100, 0, 100);
414 fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent", 100, 0, 100);
a89e2005 415 fHistEbefore = new TH1F("Ebefore", "Ebefore", 100, 0, 100);
416 fHistEafter = new TH1F("Eafter", "Eafter", 100, 0, 100);
d6f334b5 417 fHistEoPCent = new TH2F("EoPCent", "EoPCent", 100, 0, 100, 1000, 0, 10);
1f6c8e39 418 fHistNMatchCent = new TH2F("NMatchesCent", "NMatchesCent", 100, 0, 100, 11, -0.5, 10.5);
1f6c8e39 419
420 fOutput->Add(fHistNclusMatchvsCent);
421 fOutput->Add(fHistNclusvsCent);
422 fOutput->Add(fHistEbefore);
423 fOutput->Add(fHistEafter);
424 fOutput->Add(fHistEoPCent);
425 fOutput->Add(fHistNMatchCent);
1f6c8e39 426 fOutput->Add(fHistCentrality);
1f6c8e39 427
428 PostData(1, fOutput);
0b777a09 429}
430
1f6c8e39 431//________________________________________________________________________
1f6c8e39 432void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster, Double_t &totalTrkP, Int_t &Nmatches)
433{
c7d52659 434 // Do the loop over matched tracks for cluster emccluster.
1f6c8e39 435
436 AliVCluster *cluster = emccluster->GetCluster();
437 Int_t iClus = emccluster->IdInCollection();
438 Double_t energyclus = cluster->E();
439
440 // loop over matched tracks
c7d52659 441 const Int_t Ntrks = emccluster->GetNumberOfMatchedObj();
1f6c8e39 442 for (Int_t i = 0; i < Ntrks; ++i) {
443 Int_t iTrack = emccluster->GetMatchedObjId(i);
444 Double_t dR = emccluster->GetMatchedObjDistance(i);
445
4f833c06 446 AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iTrack));
1f6c8e39 447 if (!emctrack)
448 continue;
1f6c8e39 449 AliVTrack *track = emctrack->GetTrack();
450 if (!track)
451 continue;
452 if (!AcceptTrack(track))
453 continue;
454
875b35fc 455 // check if track also points to cluster
456 if (fDoTrackClus && (track->GetEMCALcluster()) != iClus)
457 continue;
458
1f6c8e39 459 Double_t etadiff = 999;
460 Double_t phidiff = 999;
461 AliPicoTrack::GetEtaPhiDiff(track, cluster, phidiff, etadiff);
462
463 Double_t mom = track->P();
464 Int_t mombin = GetMomBin(mom);
465 Int_t centbinch = fCentBin;
c7d52659 466 if (track->Charge()<0)
1f6c8e39 467 centbinch += 4;
468
469 if (fCreateHisto) {
470 Int_t etabin = 0;
471 if(track->Eta() > 0) etabin=1;
1f6c8e39 472 fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
473 }
474
c7d52659 475 Double_t etaCut = 0.0;
476 Double_t phiCutlo = 0.0;
477 Double_t phiCuthi = 0.0;
1f6c8e39 478
4f833c06 479 if (fPhiMatch > 0) {
c7d52659 480 phiCutlo = -fPhiMatch;
875b35fc 481 phiCuthi = +fPhiMatch;
482 } else {
c7d52659 483 phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
484 phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
1f6c8e39 485 }
486
487 if (fEtaMatch > 0) {
c7d52659 488 etaCut = fEtaMatch;
875b35fc 489 } else {
c7d52659 490 etaCut = GetEtaSigma(mombin);
3a385c49 491 }
3a385c49 492
c7d52659 493 if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
875b35fc 494 ++Nmatches;
495 totalTrkP += track->P();
496
497 if (fCreateHisto) {
498 if (fHadCorr > 1 && mombin > -1) {
499 fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
500 }
1f6c8e39 501 }
502 }
503 }
3a385c49 504}
505
0b777a09 506//________________________________________________________________________
1f6c8e39 507Bool_t AliHadCorrTask::Run()
0b777a09 508{
1f6c8e39 509 // Run the hadronic correction
0b777a09 510
d6f334b5 511 // post output in event if not yet present
71c1dd76 512 fOutClusters->Delete();
0b777a09 513 if (!(InputEvent()->FindListObject(fOutCaloName)))
514 InputEvent()->AddObject(fOutClusters);
d6f334b5 515
516 // delete output
517 fOutClusters->Delete();
0b777a09 518
d6f334b5 519 // esd or aod mode
520 Bool_t esdMode = kTRUE;
521 if (dynamic_cast<AliAODEvent*>(InputEvent()))
522 esdMode = kFALSE;
0b777a09 523
c7d52659 524 if (esdMode) { // optimization in case autobranch loading is off
525 AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
526 if (fCaloName == "CaloClusters")
527 am->LoadBranch("CaloClusters");
528 if (fTracksName == "Tracks")
529 am->LoadBranch("Tracks");
9cf91fa5 530 am->LoadBranch("Centrality.");
c7d52659 531 }
0b777a09 532
85835c5b 533 if (fCreateHisto) {
1f6c8e39 534 fHistCentrality->Fill(fCent);
85835c5b 535 }
32bf39af 536
85835c5b 537 // loop over all clusters
538 const Int_t Nclus = fCaloClusters->GetEntries();
a89e2005 539 for (Int_t iClus = 0, clusCount=0; iClus < Nclus; ++iClus) {
4f833c06 540 AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
1f6c8e39 541 if (!emccluster)
d6f334b5 542 continue;
543
1f6c8e39 544 AliVCluster *cluster = emccluster->GetCluster();
545 if (!cluster)
546 continue;
547 if (!AcceptCluster(cluster))
d6f334b5 548 continue;
549
1f6c8e39 550 Double_t energyclus = 0;
62206c47 551 if (fCreateHisto)
552 fHistEbefore->Fill(fCent, cluster->E());
32bf39af 553
d6f334b5 554 // apply correction / subtraction
bcf52e4a 555 if (fHadCorr > 0) {
d6f334b5 556 // to subtract only the closest track set fHadCor to a %
557 // to subtract all tracks within the cut set fHadCor to %+1
1f6c8e39 558 if (fHadCorr > 1)
559 energyclus = ApplyHadCorrAllTracks(emccluster, fHadCorr - 1);
560 else
561 energyclus = ApplyHadCorrOneTrack(emccluster, fHadCorr);
0b777a09 562 }
bcf52e4a 563
564 if (energyclus < 0)
d6f334b5 565 energyclus = 0;
d6f334b5 566
bcf52e4a 567 if (energyclus > 0) { // create corrected cluster
d6f334b5 568 AliVCluster *oc;
4f833c06 569 if (esdMode) {
570 AliESDCaloCluster *ec = dynamic_cast<AliESDCaloCluster*>(cluster);
571 if (!ec)
572 continue;
573 oc = new ((*fOutClusters)[clusCount]) AliESDCaloCluster(*ec);
574 } else {
575 AliAODCaloCluster *ac = dynamic_cast<AliAODCaloCluster*>(cluster);
576 if (!ac)
577 continue;
578 oc = new ((*fOutClusters)[clusCount]) AliAODCaloCluster(*ac);
579 }
d6f334b5 580 oc->SetE(energyclus);
1f6c8e39 581
d6f334b5 582 ++clusCount;
1f6c8e39 583
584 if (fCreateHisto)
585 fHistEafter->Fill(fCent, energyclus);
d6f334b5 586 }
0b777a09 587 }
1f6c8e39 588
589 return kTRUE;
590}
bcf52e4a 591
1f6c8e39 592//________________________________________________________________________
593Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Double_t hadCorr)
594{
4f833c06 595 // Apply the hadronic correction with one track only.
596
1f6c8e39 597 AliVCluster *cluster = emccluster->GetCluster();
181a0a09 598 Double_t energyclus = cluster->E();
599 Int_t iMin = emccluster->GetMatchedObjId();
1f6c8e39 600 if (iMin < 0)
601 return energyclus;
602
4f833c06 603 AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
1f6c8e39 604 if (!emctrack)
605 return energyclus;
606
181a0a09 607 // check if track also points to cluster
608 Int_t cid = emctrack->GetMatchedObjId();
875b35fc 609 if (fDoTrackClus && (cid!=emccluster->IdInCollection()))
c7d52659 610 return energyclus;
611
1f6c8e39 612 AliVTrack *track = emctrack->GetTrack();
613 if (!track)
614 return energyclus;
615
616 Double_t mom = track->P();
617 if (mom == 0)
618 return energyclus;
619
620 Double_t dRmin = emccluster->GetMatchedObjDistance();
621 Double_t dEtaMin = 1e9;
622 Double_t dPhiMin = 1e9;
623 AliPicoTrack::GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
624
625 Int_t mombin = GetMomBin(mom);
626 Int_t centbinch = fCentBin;
c7d52659 627 if (track->Charge()<0)
1f6c8e39 628 centbinch += 4;
629
181a0a09 630 // plot some histograms if switched on
1f6c8e39 631 if (fCreateHisto) {
632 Int_t etabin = 0;
633 if(track->Eta() > 0)
634 etabin = 1;
635
636 fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
637
638 if (mom > 0) {
639 fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
640 fHistEoPCent->Fill(fCent, energyclus / mom);
641 fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
642 }
643 }
644
181a0a09 645 // define eta/phi cuts
c7d52659 646 Double_t etaCut = 0.0;
647 Double_t phiCutlo = 0.0;
648 Double_t phiCuthi = 0.0;
1f6c8e39 649 if (fPhiMatch > 0) {
c7d52659 650 phiCutlo = -fPhiMatch;
875b35fc 651 phiCuthi = +fPhiMatch;
1f6c8e39 652 }
653 else {
c7d52659 654 phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
655 phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
1f6c8e39 656 }
657 if(fEtaMatch > 0) {
c7d52659 658 etaCut = fEtaMatch;
1f6c8e39 659 }
660 else {
c7d52659 661 etaCut = GetEtaSigma(mombin);
1f6c8e39 662 }
663
181a0a09 664 // apply the correction if the track is in the eta/phi window
c7d52659 665 if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
875b35fc 666 energyclus -= hadCorr * mom;
1f6c8e39 667 }
668
669 return energyclus;
0b777a09 670}
4711c521 671
0b777a09 672//________________________________________________________________________
1f6c8e39 673Double_t AliHadCorrTask::ApplyHadCorrAllTracks(AliEmcalParticle *emccluster, Double_t hadCorr)
0b777a09 674{
4f833c06 675 // Apply the hadronic correction with all tracks.
676
1f6c8e39 677 AliVCluster *cluster = emccluster->GetCluster();
678
679 Double_t energyclus = cluster->E();
680 Double_t cNcells = cluster->GetNCells();
681
682 Double_t totalTrkP = 0.0; // count total track momentum
683 Int_t Nmatches = 0; // count total number of matches
684
181a0a09 685 // do the loop over the matched tracks and get the number of matches and the total momentum
1f6c8e39 686 DoMatchedTracksLoop(emccluster, totalTrkP, Nmatches);
687
a9880dcd 688 if(fCreateHisto && Nmatches == 0) {
689 fHistNclusvsCent->Fill(fCent);
690 fHistNCellsEnergy[fCentBin][0]->Fill(energyclus, cNcells);
691 }
692
f54a0dd2 693 if (totalTrkP <= 0)
694 return energyclus;
695
1f6c8e39 696 Double_t Esub = hadCorr * totalTrkP;
697
62206c47 698 Double_t clusEexcl = fEexclCell * cNcells;
699
700 if (energyclus < clusEexcl)
701 clusEexcl = energyclus;
702
703 if (Esub > energyclus)
704 Esub = energyclus;
705
181a0a09 706 // applying Peter's proposed algorithm
707 // never subtract the full energy of the cluster
62206c47 708 if ((energyclus - Esub) < clusEexcl)
709 Esub = (energyclus - clusEexcl);
710
f54a0dd2 711 Double_t EoP = energyclus / totalTrkP;
1f6c8e39 712
181a0a09 713 // plot some histograms if switched on
714 if (fCreateHisto) {
1f6c8e39 715 fHistNMatchCent->Fill(fCent, Nmatches);
716 fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
717
718 if(Nmatches > 0)
719 fHistNclusMatchvsCent->Fill(fCent);
720
a9880dcd 721 if(Nmatches == 1)
1f6c8e39 722 fHistNCellsEnergy[fCentBin][1]->Fill(energyclus, cNcells);
723 else if(Nmatches == 2)
724 fHistNCellsEnergy[fCentBin][2]->Fill(energyclus, cNcells);
a9880dcd 725 else if(Nmatches > 2)
1f6c8e39 726 fHistNCellsEnergy[fCentBin][3]->Fill(energyclus, cNcells);
727
728 if (EoP > 0) {
729 fHistEoPCent->Fill(fCent, EoP);
730 fHistMatchEvsP[fCentBin]->Fill(energyclus, EoP);
731 }
732
733 if (Nmatches == 1) {
734 Int_t iMin = emccluster->GetMatchedObjId();
4f833c06 735 AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
1f6c8e39 736 AliVTrack *track = emctrack->GetTrack();
737 Int_t centbinchm = fCentBin;
c7d52659 738 if (track->Charge()<0)
1f6c8e39 739 centbinchm += 4;
740
f54a0dd2 741 fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
1f6c8e39 742 fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
743 }
744 }
1f6c8e39 745
181a0a09 746 // apply the correction
1f6c8e39 747 energyclus -= Esub;
748
749 return energyclus;
0b777a09 750}
751