Coverity fix: CID 22038
[u/mrichter/AliRoot.git] / PWG / EMCAL / 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 #include "AliEMCALGeometry.h"
22 #include "AliParticleContainer.h"
23 #include "AliClusterContainer.h"
24
25 #include "AliHadCorrTask.h"
26
27 ClassImp(AliHadCorrTask)
28
29 //________________________________________________________________________
30 AliHadCorrTask::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 //________________________________________________________________________
82 AliHadCorrTask::AliHadCorrTask(const char *name, Bool_t histo) : 
83   AliAnalysisTaskEmcal(name, histo),
84   fOutCaloName("CaloClustersCorr"),
85   fPhiMatch(0.05),
86   fEtaMatch(0.025),
87   fDoTrackClus(kTRUE),
88   fHadCorr(0),
89   fEexclCell(0),
90   fDoExact(kFALSE),
91   fEsdMode(kTRUE),
92   fOutClusters(0),
93   fHistMatchEtaPhiAll(0),
94   fHistMatchEtaPhiAllTr(0),
95   fHistMatchEtaPhiAllCl(0),
96   fHistNclusvsCent(0),
97   fHistNclusMatchvsCent(0),
98   fHistEbefore(0),
99   fHistEafter(0),
100   fHistEoPCent(0),
101   fHistNMatchCent(0),
102   fHistNClusMatchCent(0)
103 {
104   // Standard constructor.
105
106   for(Int_t i=0; i<8; i++) {
107     fHistEsubPch[i]    = 0;
108     fHistEsubPchRat[i] = 0;
109     fHistEsubPchRatAll[i] = 0;
110
111     if (i<4) {
112       fHistMatchEvsP[i]    = 0;
113       fHistMatchdRvsEP[i]  = 0;
114       fHistNMatchEnergy[i] = 0;
115
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 //________________________________________________________________________
138 AliHadCorrTask::~AliHadCorrTask()
139 {
140   // Destructor
141 }
142
143 //________________________________________________________________________
144 UInt_t AliHadCorrTask::GetMomBin(Double_t p) const
145 {
146   // Get momenum bin.
147
148   UInt_t pbin=0;
149   if (p<0.5) 
150     pbin=0;
151   else if (p>=0.5 && p<1.0) 
152     pbin=1;
153   else if (p>=1.0 && p<1.5) 
154     pbin=2;
155   else if (p>=1.5 && p<2.) 
156     pbin=3;
157   else if (p>=2. && p<3.) 
158     pbin=4;
159   else if (p>=3. && p<4.) 
160     pbin=5;
161   else if (p>=4. && p<5.) 
162     pbin=6;
163   else if (p>=5. && p<8.) 
164     pbin=7;
165   else 
166     pbin=8;
167
168   return pbin;
169 }
170
171 //________________________________________________________________________
172 Double_t AliHadCorrTask::GetEtaSigma(Int_t pbin) const
173 {
174   // Get sigma in eta.
175
176   Double_t EtaSigma[9]={0.0097,0.0075,0.0059,0.0055,0.0053,0.005,0.005,0.0045,0.0042};
177   return 2.0*EtaSigma[pbin];
178 }
179
180 //________________________________________________________________________
181 Double_t AliHadCorrTask::GetPhiMean(Int_t pbin, Int_t centbin) const
182 {
183   // Get phi mean.
184
185   if (centbin==0) { 
186     Double_t PhiMean[9]={0.036,
187                          0.021,
188                          0.0121,
189                          0.0084,
190                          0.0060,
191                          0.0041,
192                          0.0031,
193                          0.0022,
194                          0.001};
195     return PhiMean[pbin];
196   } else if (centbin==1) { 
197     Double_t PhiMean[9]={0.036,
198                          0.021,
199                          0.0121,
200                          0.0084,
201                          0.0060,
202                          0.0041,
203                          0.0031,
204                          0.0022,
205                          0.001};
206     return PhiMean[pbin];
207   } else if (centbin==2) { 
208     Double_t PhiMean[9]={0.036,
209                          0.021,
210                          0.0121,
211                          0.0084,
212                          0.0060,
213                          0.0041,
214                          0.0031,
215                          0.0022,
216                          0.001};
217     return PhiMean[pbin];
218   } else if (centbin==3) { 
219     Double_t PhiMean[9]={0.036,
220                          0.021,
221                          0.0121,
222                          0.0084,
223                          0.0060,
224                          0.0041,
225                          0.0031,
226                          0.0022,
227                          0.001};  
228     return PhiMean[pbin];
229   } else if (centbin==4) { 
230     Double_t PhiMean[9]={0.036,
231                          0.021,
232                          0.0127,
233                          0.0089,
234                          0.0068,
235                          0.0049,
236                          0.0038,
237                          0.0028,
238                          0.0018};
239     return PhiMean[pbin]*(-1.);
240   } else if (centbin==5) { 
241     Double_t PhiMean[9]={0.036,
242                          0.021,
243                          0.0127,
244                          0.0089,
245                          0.0068,
246                          0.0048,
247                          0.0038,
248                          0.0028,
249                          0.0018};
250     return PhiMean[pbin]*(-1.);
251   } else if (centbin==6) { 
252     Double_t PhiMean[9]={0.036,
253                          0.021,
254                          0.0127,
255                          0.0089,
256                          0.0068,
257                          0.0045,
258                          0.0035,
259                          0.0028,
260                          0.0018};
261     return PhiMean[pbin]*(-1.);
262   } else if (centbin==7) { 
263     Double_t PhiMean[9]={0.036,
264                          0.021,
265                          0.0127,
266                          0.0089,
267                          0.0068,
268                          0.0043,
269                          0.0035,
270                          0.0028,
271                          0.0018};
272     return PhiMean[pbin]*(-1.);
273   }
274
275   return 0;
276 }
277
278 //________________________________________________________________________
279 Double_t AliHadCorrTask::GetPhiSigma(Int_t pbin, Int_t centbin) const
280 {
281   // Get phi sigma.
282
283   if (centbin==0) { 
284     Double_t PhiSigma[9]={0.0221,
285                           0.0128,
286                           0.0074,
287                           0.0064,
288                           0.0059,
289                           0.0055,
290                           0.0052,
291                           0.0049,
292                           0.0045};
293     return 2.*PhiSigma[pbin];
294   } else if (centbin==1) { 
295     Double_t PhiSigma[9]={0.0217,
296                           0.0120,
297                           0.0076,
298                           0.0066,
299                           0.0062,
300                           0.0058,
301                           0.0054,
302                           0.0054,
303 0.0045};
304     return 2.*PhiSigma[pbin];
305   } else if (centbin==2) { 
306     Double_t PhiSigma[9]={0.0211,
307                           0.0124,
308                           0.0080,
309                           0.0070,
310                           0.0067,
311                           0.0061,
312                           0.0059,
313                           0.0054,
314                           0.0047};
315     return 2.*PhiSigma[pbin];
316   } else if (centbin==3) { 
317     Double_t PhiSigma[9]={0.0215,
318                           0.0124,
319                           0.0082,
320                           0.0073,
321                           0.0069,
322                           0.0064,
323                           0.0060,
324                           0.0055,
325                           0.0047};  
326     return 2.*PhiSigma[pbin];
327   } else if (centbin==4) { 
328     Double_t PhiSigma[9]={0.0199,
329                           0.0108,
330                           0.0072,
331                           0.0071,
332                           0.0060,
333                           0.0055,
334                           0.0052,
335                           0.0049,
336                           0.0045};
337     return 2.*PhiSigma[pbin];
338   } else if (centbin==5) { 
339     Double_t PhiSigma[9]={0.0200,
340                           0.0110,
341                           0.0074,
342                           0.0071,
343                           0.0064,
344                           0.0059,
345                           0.0055,
346                           0.0052,
347                           0.0045};
348     return 2.*PhiSigma[pbin];
349   } else if (centbin==6) { 
350     Double_t PhiSigma[9]={0.0202,
351                           0.0113,
352                           0.0077,
353                           0.0071,
354                           0.0069,
355                           0.0064,
356                           0.0060,
357                           0.0055,
358                           0.0050};
359     return 2.*PhiSigma[pbin];
360   } else if (centbin==7) { 
361     Double_t PhiSigma[9]={0.0205,
362                           0.0113,
363                           0.0080,
364                           0.0074,
365                           0.0078,
366                           0.0067,
367                           0.0062,
368                           0.0055,
369                           0.0050};
370     return 2.*PhiSigma[pbin];
371   }
372
373   return 0;
374 }
375
376 //________________________________________________________________________
377 void 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 //________________________________________________________________________
514 void 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 //________________________________________________________________________
553 void 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();
585       if (!cluster) continue;
586       if (!cluster->IsEMCAL()) continue;
587
588       Double_t etadiff = 999;
589       Double_t phidiff = 999;
590       AliPicoTrack::GetEtaPhiDiff(track, cluster, phidiff, etadiff);
591       fHistMatchEtaPhiAllTr->Fill(etadiff,phidiff);
592
593       if (TMath::Abs(phidiff) < fPhiMatch && TMath::Abs(etadiff) < fEtaMatch) NmatchClus++;
594     }
595     fHistNClusMatchCent->Fill(fCent, NmatchClus);
596   }
597 }
598
599 //________________________________________________________________________
600 void 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 //________________________________________________________________________
683 Bool_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 //________________________________________________________________________
749 Double_t AliHadCorrTask::ApplyHadCorrOneTrack(AliEmcalParticle *emccluster, Double_t hadCorr) 
750 {
751   // Apply the hadronic correction with one track only.
752
753   AliVCluster *cluster = emccluster->GetCluster();
754   if (!cluster) return 0;
755
756   Double_t energyclus  = cluster->E();
757   Int_t    iMin        = emccluster->GetMatchedObjId();
758   if (iMin < 0)
759     return energyclus;
760
761   AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
762   AliEmcalParticle *emctrack = static_cast<AliEmcalParticle*>(tracks->GetParticle(iMin));
763   if (!emctrack) return energyclus;
764
765   AliVTrack *track = emctrack->GetTrack();
766   if (!track) return energyclus;
767
768   Double_t mom = track->P();
769   if (mom < 1e-6) return energyclus;
770
771   Double_t dEtaMin    = 1e9;
772   Double_t dPhiMin    = 1e9;
773   AliPicoTrack::GetEtaPhiDiff(track, cluster, dPhiMin, dEtaMin);
774   if (fCreateHisto)
775     fHistMatchEtaPhiAllCl->Fill(dEtaMin, dPhiMin);
776
777   // check if track also points to cluster
778   Int_t cid = emctrack->GetMatchedObjId();
779   if (fDoTrackClus && (cid!=emccluster->IdInCollection())) return energyclus;
780
781   UInt_t mombin = GetMomBin(mom);
782   Int_t centbinch = fCentBin;
783   if (track->Charge()<0) centbinch += fNcentBins;
784
785   // plot some histograms if switched on
786   if (fCreateHisto) {
787     Int_t etabin = 0;
788     if(track->Eta() > 0) 
789       etabin = 1;
790             
791     fHistMatchEtaPhi[centbinch][mombin][etabin]->Fill(dEtaMin, dPhiMin);
792     fHistMatchEtaPhiAll->Fill(dEtaMin, dPhiMin);
793     
794     if (mom > 0) {
795       Double_t dRmin      = emccluster->GetMatchedObjDistance();
796       fHistMatchEvsP[fCentBin]->Fill(energyclus, energyclus / mom);
797       fHistEoPCent->Fill(fCent, energyclus / mom);
798       fHistMatchdRvsEP[fCentBin]->Fill(dRmin, energyclus / mom);
799     }
800   }
801            
802   // define eta/phi cuts
803   Double_t etaCut   = 0.0;
804   Double_t phiCutlo = 0.0;
805   Double_t phiCuthi = 0.0;
806   if (fPhiMatch > 0) {
807     phiCutlo = -fPhiMatch;
808     phiCuthi = +fPhiMatch;
809   }
810   else {
811     phiCutlo = GetPhiMean(mombin, centbinch) - GetPhiSigma(mombin, fCentBin);
812     phiCuthi = GetPhiMean(mombin, centbinch) + GetPhiSigma(mombin, fCentBin);
813   }
814   if(fEtaMatch > 0) {
815     etaCut = fEtaMatch;
816   }
817   else {
818     etaCut = GetEtaSigma(mombin);
819   }
820   
821   // apply the correction if the track is in the eta/phi window
822   if ((dPhiMin < phiCuthi && dPhiMin > phiCutlo) && TMath::Abs(dEtaMin) < etaCut) {
823     energyclus -= hadCorr * mom;
824   }
825
826   return energyclus;
827 }
828
829 //________________________________________________________________________
830 Double_t AliHadCorrTask::ApplyHadCorrAllTracks(AliEmcalParticle *emccluster, Double_t hadCorr) 
831 {
832   // Apply the hadronic correction with all tracks.
833
834   AliParticleContainer *tracks = static_cast<AliParticleContainer*>(fParticleCollArray.At(0));
835
836   AliVCluster *cluster = emccluster->GetCluster();
837   
838   Double_t energyclus = cluster->E();
839   Double_t cNcells = cluster->GetNCells();
840   
841   Double_t totalTrkP  = 0.0; // count total track momentum
842   Int_t    Nmatches   = 0;   // count total number of matches
843
844   Double_t trkPMCfrac   = 0.0; // count total track momentum
845   Int_t    NMCmatches   = 0;   // count total number of matches
846   
847   // do the loop over the matched tracks and get the number of matches and the total momentum
848   DoMatchedTracksLoop(emccluster, totalTrkP, Nmatches, trkPMCfrac, NMCmatches);
849
850   Double_t Esub = hadCorr * totalTrkP;
851         
852   if (Esub > energyclus) Esub = energyclus;
853         
854   // applying Peter's proposed algorithm
855   // never subtract the full energy of the cluster 
856   Double_t clusEexcl = fEexclCell * cNcells;
857   if (energyclus < clusEexcl) clusEexcl = energyclus;
858   if ((energyclus - Esub) < clusEexcl) Esub = (energyclus - clusEexcl);
859
860   // embedding
861   Double_t EsubMC       = 0;
862   Double_t EsubBkg      = 0;
863   Double_t EclusMC      = 0;
864   Double_t EclusBkg     = 0;
865   Double_t EclusCorr    = 0;
866   Double_t EclusMCcorr  = 0;
867   Double_t EclusBkgcorr = 0;
868   Double_t overSub = 0;
869   if (fIsEmbedded) {
870     EsubMC   = hadCorr * totalTrkP * trkPMCfrac;
871     EsubBkg  = hadCorr * totalTrkP - EsubMC;
872     EclusMC  = energyclus * cluster->GetMCEnergyFraction();
873     EclusBkg = energyclus - EclusMC;
874  
875     if (energyclus > Esub)
876       EclusCorr = energyclus - Esub;
877
878     if (EclusMC > EsubMC)
879       EclusMCcorr = EclusMC - EsubMC;
880   
881     if (EclusBkg > EsubBkg)
882       EclusBkgcorr = EclusBkg - EsubBkg;
883   
884     overSub = EclusMCcorr + EclusBkgcorr - EclusCorr;
885   }
886
887   // plot some histograms if switched on
888   if (fCreateHisto) {
889     fHistNMatchCent->Fill(fCent, Nmatches);
890     fHistNMatchEnergy[fCentBin]->Fill(energyclus, Nmatches);
891     
892     if(Nmatches > 0) 
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 }