protections against failures in deleting event content
[u/mrichter/AliRoot.git] / PWGJE / 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   fHistCentrality(0)
43 {
44   // Default constructor.
45
46   for(Int_t i=0; i<8; i++) {
47       fHistEsubPch[i]    = 0;
48       fHistEsubPchRat[i] = 0;
49     for(Int_t j=0; j<4; j++) {
50       fHistNCellsEnergy[i][j] = 0;
51     }
52     if (i<4) {
53       fHistMatchEvsP[i]    = 0;
54       fHistMatchdRvsEP[i]  = 0;
55       fHistNMatchEnergy[i] = 0;
56     }
57     for(Int_t j=0; j<9; j++) {
58       for(Int_t k=0; k<2; k++) {
59         fHistMatchEtaPhi[i][j][k] = 0;
60       }
61     }
62   } 
63 }
64
65 //________________________________________________________________________
66 AliHadCorrTask::AliHadCorrTask(const char *name, Bool_t histo) : 
67   AliAnalysisTaskEmcal(name, histo),
68   fOutCaloName("CaloClustersCorr"),
69   fPhiMatch(0.05),
70   fEtaMatch(0.025),
71   fDoTrackClus(1),
72   fHadCorr(0),
73   fEexclCell(0),
74   fOutClusters(0),
75   fHistNclusvsCent(0),
76   fHistNclusMatchvsCent(0),
77   fHistEbefore(0),
78   fHistEafter(0),
79   fHistEoPCent(0),
80   fHistNMatchCent(0),
81   fHistCentrality(0)
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++) {
90         fHistNCellsEnergy[i][j] = 0;
91       }
92     }
93     fHistEsubPch[i]    = 0;
94     fHistEsubPchRat[i] = 0;
95     for(Int_t j=0; j<9; j++) {
96       for(Int_t k=0; k<2; k++) {
97         fHistMatchEtaPhi[i][j][k] = 0;
98       }
99     }
100   } 
101
102   fBranchNames="ESD:AliESDRun.,AliESDHeader.,PrimaryVertex.";
103 }
104
105 //________________________________________________________________________
106 AliHadCorrTask::~AliHadCorrTask()
107 {
108   // Destructor
109 }
110
111 //________________________________________________________________________
112 Int_t AliHadCorrTask::GetMomBin(Double_t p) const
113 {
114   // Get momenum bin.
115
116   Int_t pbin=-1;
117   if (p<0.5) 
118     pbin=0;
119   else if (p>=0.5 && p<1.0) 
120     pbin=1;
121   else if (p>=1.0 && p<1.5) 
122     pbin=2;
123   else if (p>=1.5 && p<2.) 
124     pbin=3;
125   else if (p>=2. && p<3.) 
126     pbin=4;
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;
135
136   return pbin;
137 }
138
139 //________________________________________________________________________
140 Double_t AliHadCorrTask::GetEtaSigma(Int_t pbin) const
141 {
142   // Get sigma in eta.
143
144   Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
145   return 2.0*EtaSigma[pbin];
146 }
147
148 //________________________________________________________________________
149 Double_t AliHadCorrTask::GetPhiMean(Int_t pbin, Int_t centbin) const
150 {
151   // Get phi mean.
152
153   if (centbin==0) { 
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];
164   } else if (centbin==1) { 
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];
175   } else if (centbin==2) { 
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];
186   } else if (centbin==3) { 
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];
197   } else if (centbin==4) { 
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.);
208   } else if (centbin==5) { 
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.);
219   } else if (centbin==6) { 
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.);
230   } else if (centbin==7) { 
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;
244 }
245
246 //________________________________________________________________________
247 Double_t AliHadCorrTask::GetPhiSigma(Int_t pbin, Int_t centbin) const
248 {
249   // Get phi sigma.
250
251   if (centbin==0) { 
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];
262   } else if (centbin==1) { 
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,
271 0.0045};
272     return 2.*PhiSigma[pbin];
273   } else if (centbin==2) { 
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];
284   } else if (centbin==3) { 
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];
295   } else if (centbin==4) { 
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];
306   } else if (centbin==5) { 
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];
317   } else if (centbin==6) { 
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];
328   } else if (centbin==7) { 
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;
342 }
343
344 //________________________________________________________________________
345 void AliHadCorrTask::UserCreateOutputObjects()
346 {
347   // Create my user objects.
348
349   AliVEventHandler* handler = AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler();
350   if (!handler) {
351     AliError(Form("%s: Input handler not available!", GetName()));
352     return;
353   }
354  
355   if (handler->InheritsFrom("AliESDInputHandler")) {
356     fOutClusters = new TClonesArray("AliESDCaloCluster");
357   } else if (handler->InheritsFrom("AliAODInputHandler")) {
358     fOutClusters = new TClonesArray("AliAODCaloCluster");
359   } else {
360     AliError(Form("%s: Input handler not recognized!", GetName()));
361     return;
362   }
363   fOutClusters->SetName(fOutCaloName);
364
365   if (!fCreateHisto)
366     return;
367
368   OpenFile(1);
369   fOutput = new TList();
370   fOutput->SetOwner();
371
372   TString name;
373
374   for(Int_t icent=0; icent<8; ++icent) {
375     for(Int_t ipt=0; ipt<9; ++ipt) {
376       for(Int_t ieta=0; ieta<2; ++ieta) {
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       }
381     }
382
383     for(Int_t itrk=0; itrk<4; ++itrk) {
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]);
396     
397     if (icent<4) {
398       name = Form("fHistMatchEvsP_%i",icent);
399       fHistMatchEvsP[icent] = new TH2F(name, name, 400, 0., 200., 1000, 0., 10.);
400       fOutput->Add(fHistMatchEvsP[icent]);
401
402       name = Form("fHistMatchdRvsEP_%i",icent);
403       fHistMatchdRvsEP[icent] = new TH2F(name, name, 1000, 0., 1., 1000, 0., 10.);
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]);
409     }
410   }
411
412   fHistCentrality       = new TH1F("fHistCentrality",  "Centrality",       100, 0, 100);
413   fHistNclusvsCent      = new TH1F("Nclusvscent",      "NclusVsCent",      100, 0, 100);
414   fHistNclusMatchvsCent = new TH1F("NclusMatchvscent", "NclusMatchVsCent", 100, 0, 100);
415   fHistEbefore          = new TH1F("Ebefore",          "Ebefore",          100, 0, 100);
416   fHistEafter           = new TH1F("Eafter",           "Eafter",           100, 0, 100);
417   fHistEoPCent          = new TH2F("EoPCent",          "EoPCent",          100, 0, 100, 1000, 0,   10);
418   fHistNMatchCent       = new TH2F("NMatchesCent",     "NMatchesCent",     100, 0, 100, 11, -0.5, 10.5);
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);
426   fOutput->Add(fHistCentrality);
427
428   PostData(1, fOutput);
429 }
430
431 //________________________________________________________________________
432 void AliHadCorrTask::DoMatchedTracksLoop(AliEmcalParticle *emccluster, Double_t &totalTrkP, Int_t &Nmatches) 
433 {
434   // Do the loop over matched tracks for cluster emccluster.
435
436   AliVCluster *cluster = emccluster->GetCluster();
437   Int_t iClus = emccluster->IdInCollection();
438   Double_t energyclus = cluster->E();
439
440   // loop over matched tracks
441   const Int_t Ntrks = emccluster->GetNumberOfMatchedObj();
442   for (Int_t i = 0; i < Ntrks; ++i) {
443     Int_t    iTrack = emccluster->GetMatchedObjId(i);
444     Double_t dR     = emccluster->GetMatchedObjDistance(i);
445     
446     AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iTrack));
447     if (!emctrack)
448       continue;
449     AliVTrack *track = emctrack->GetTrack();
450     if (!track)
451       continue;
452     if (!AcceptTrack(track))
453       continue;
454
455     // check if track also points to cluster
456     if (fDoTrackClus && (track->GetEMCALcluster()) != iClus)
457       continue;
458
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;
466     if (track->Charge()<0) 
467       centbinch += 4;
468
469     if (fCreateHisto) {
470       Int_t etabin = 0;
471       if(track->Eta() > 0) etabin=1;
472       fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(etadiff, phidiff);
473     }
474     
475     Double_t etaCut   = 0.0;
476     Double_t phiCutlo = 0.0;
477     Double_t phiCuthi = 0.0;
478
479     if (fPhiMatch > 0) {
480       phiCutlo = -fPhiMatch;
481       phiCuthi = +fPhiMatch;
482     } else {
483       phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
484       phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
485     }
486
487     if (fEtaMatch > 0) {
488       etaCut = fEtaMatch;
489     } else {
490       etaCut = GetEtaSigma(mombin);
491     }
492
493     if ((phidiff < phiCuthi && phidiff > phiCutlo) && TMath::Abs(etadiff) < etaCut) {
494       ++Nmatches;
495       totalTrkP += track->P();
496
497       if (fCreateHisto) {
498         if (fHadCorr > 1 && mombin > -1) {
499           fHistMatchdRvsEP[fCentBin]->Fill(dR, energyclus / mom);
500         }
501       }
502     }
503   }
504 }
505
506 //________________________________________________________________________
507 Bool_t AliHadCorrTask::Run() 
508 {
509   // Run the hadronic correction
510
511   // post output in event if not yet present
512   fOutClusters->Delete();
513   if (!(InputEvent()->FindListObject(fOutCaloName)))
514     InputEvent()->AddObject(fOutClusters);
515   
516   // delete output
517   fOutClusters->Delete();
518
519   // esd or aod mode
520   Bool_t esdMode = kTRUE;
521   if (dynamic_cast<AliAODEvent*>(InputEvent()))
522       esdMode = kFALSE;
523
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");
530     am->LoadBranch("Centrality.");      
531   }
532
533   if (fCreateHisto) {
534     fHistCentrality->Fill(fCent);
535   }
536
537    // loop over all clusters
538   const Int_t Nclus = fCaloClusters->GetEntries();
539   for (Int_t iClus = 0, clusCount=0; iClus < Nclus; ++iClus) {
540     AliEmcalParticle *emccluster = static_cast<AliEmcalParticle*>(fCaloClusters->At(iClus));
541     if (!emccluster)
542       continue;
543
544     AliVCluster *cluster = emccluster->GetCluster();
545     if (!cluster)
546       continue;
547     if (!AcceptCluster(cluster))
548       continue;
549
550     Double_t energyclus = 0;
551     if (fCreateHisto)
552       fHistEbefore->Fill(fCent, cluster->E());
553   
554     // apply correction / subtraction
555     if (fHadCorr > 0) {
556       // to subtract only the closest track set fHadCor to a %
557       // to subtract all tracks within the cut set fHadCor to %+1
558       if (fHadCorr > 1)
559         energyclus = ApplyHadCorrAllTracks(emccluster, fHadCorr - 1);   
560       else 
561         energyclus = ApplyHadCorrOneTrack(emccluster, fHadCorr);        
562     }
563
564     if (energyclus < 0) 
565       energyclus = 0;
566
567     if (energyclus > 0) { // create corrected cluster
568       AliVCluster *oc;
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       }
580       oc->SetE(energyclus);
581
582       ++clusCount;
583
584       if (fCreateHisto)
585         fHistEafter->Fill(fCent, energyclus);
586     }
587   }
588   
589   return kTRUE;
590 }
591
592 //________________________________________________________________________
593 Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Double_t hadCorr) 
594 {
595   // Apply the hadronic correction with one track only.
596
597   AliVCluster *cluster = emccluster->GetCluster();
598   Double_t energyclus  = cluster->E();
599   Int_t    iMin        = emccluster->GetMatchedObjId();
600   if (iMin < 0)
601     return energyclus;
602
603   AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
604   if (!emctrack)
605     return energyclus;
606
607   // check if track also points to cluster
608   Int_t cid = emctrack->GetMatchedObjId();
609   if (fDoTrackClus && (cid!=emccluster->IdInCollection())) 
610     return energyclus;
611
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;
627   if (track->Charge()<0) 
628     centbinch += 4;
629
630   // plot some histograms if switched on
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            
645   // define eta/phi cuts
646   Double_t etaCut   = 0.0;
647   Double_t phiCutlo = 0.0;
648   Double_t phiCuthi = 0.0;
649   if (fPhiMatch > 0) {
650     phiCutlo = -fPhiMatch;
651     phiCuthi = +fPhiMatch;
652   }
653   else {
654     phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
655     phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
656   }
657   if(fEtaMatch > 0) {
658     etaCut = fEtaMatch;
659   }
660   else {
661     etaCut = GetEtaSigma(mombin);
662   }
663   
664   // apply the correction if the track is in the eta/phi window
665   if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
666     energyclus -= hadCorr * mom;
667   }
668
669   return energyclus;
670 }
671
672 //________________________________________________________________________
673 Double_t AliHadCorrTask::ApplyHadCorrAllTracks(AliEmcalParticle *emccluster, Double_t hadCorr) 
674 {
675   // Apply the hadronic correction with all tracks.
676
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   
685   // do the loop over the matched tracks and get the number of matches and the total momentum
686   DoMatchedTracksLoop(emccluster, totalTrkP, Nmatches);
687
688   if(fCreateHisto && Nmatches == 0) {
689     fHistNclusvsCent->Fill(fCent);
690     fHistNCellsEnergy[fCentBin][0]->Fill(energyclus, cNcells);
691   }
692
693   if (totalTrkP <= 0)
694     return energyclus;
695
696   Double_t Esub = hadCorr * totalTrkP;
697
698   Double_t clusEexcl = fEexclCell * cNcells;
699
700   if (energyclus < clusEexcl) 
701     clusEexcl = energyclus;
702         
703   if (Esub > energyclus) 
704     Esub = energyclus;
705         
706   // applying Peter's proposed algorithm
707   // never subtract the full energy of the cluster 
708   if ((energyclus - Esub) < clusEexcl) 
709     Esub = (energyclus - clusEexcl);
710
711   Double_t EoP = energyclus / totalTrkP;
712
713   // plot some histograms if switched on
714   if (fCreateHisto) {
715     fHistNMatchCent->Fill(fCent, Nmatches);
716     fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
717     
718     if(Nmatches > 0) 
719       fHistNclusMatchvsCent->Fill(fCent);
720       
721     if(Nmatches == 1)
722       fHistNCellsEnergy[fCentBin][1]->Fill(energyclus, cNcells);        
723     else if(Nmatches == 2)
724       fHistNCellsEnergy[fCentBin][2]->Fill(energyclus, cNcells);
725     else if(Nmatches > 2)
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();
735       AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(fTracks->At(iMin));
736       AliVTrack *track = emctrack->GetTrack();
737       Int_t centbinchm = fCentBin;
738       if (track->Charge()<0) 
739         centbinchm += 4;
740       
741       fHistEsubPchRat[centbinchm]->Fill(totalTrkP, Esub / totalTrkP);
742       fHistEsubPch[centbinchm]->Fill(totalTrkP, Esub);
743     } 
744   }
745
746   // apply the correction
747   energyclus -= Esub;
748
749   return energyclus;
750 }
751