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