]> git.uio.no Git - u/mrichter/AliRoot.git/blob - PWGLF/totEt/AliAnalysisEtReconstructed.cxx
Merge branch 'TPCdev'
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEtReconstructed.cxx
1 //_________________________________________________________________________
2 //  Utility Class for transverse energy studies
3 //  Base class for ESD analysis
4 //  - reconstruction output
5 //  implementation file
6 //
7 //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL)
8 //_________________________________________________________________________
9
10 #include "AliAnalysisEtReconstructed.h"
11 #include "AliAnalysisEtCuts.h"
12 #include "AliESDtrack.h"
13 #include "AliEMCALTrack.h"
14 #include "AliESDCaloCluster.h"
15 #include "TVector3.h"
16 #include "TGeoGlobalMagField.h"
17 #include "AliMagF.h"
18 #include "AliVEvent.h"
19 #include "AliESDEvent.h"
20 #include "AliESDtrackCuts.h"
21 #include "AliVParticle.h"
22 #include "TDatabasePDG.h"
23 #include "TList.h"
24 #include "AliESDpid.h"
25 #include <iostream>
26 #include "TH3F.h"
27 #include "TH2F.h"
28 #include "TH2I.h"
29 #include "TH1I.h"
30 #include "TFile.h"
31 #include "AliAnalysisHadEtCorrections.h"
32 #include "AliAnalysisEtSelector.h"
33 #include "AliLog.h"
34 #include "AliCentrality.h"
35 #include "AliPHOSGeoUtils.h"
36 #include "AliPHOSGeometry.h"
37 #include "AliAnalysisEtRecEffCorrection.h"
38 #include "AliESDpid.h"
39
40
41 using namespace std;
42
43 ClassImp(AliAnalysisEtReconstructed);
44
45
46 AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
47         AliAnalysisEt()
48         ,fCorrections(0)
49         ,fPidCut(0)
50         ,nChargedHadronsMeasured(0)
51         ,nChargedHadronsTotal(0)
52         ,fHistChargedPionEnergyDeposit(0)
53         ,fHistProtonEnergyDeposit(0)
54         ,fHistAntiProtonEnergyDeposit(0)
55         ,fHistChargedKaonEnergyDeposit(0)
56         ,fHistMuonEnergyDeposit(0)
57         ,fHistRemovedEnergy(0)
58         ,fGeomCorrection(1.0)
59         ,fEMinCorrection(1.0/0.687)
60         ,fRecEffCorrection(1.0)
61         ,fClusterPositionAccepted(0)
62         ,fClusterPositionAll(0)
63         ,fClusterPositionAcceptedEnergy(0)
64         ,fClusterPositionAllEnergy(0)
65         ,fClusterEnergy(0)
66         ,fClusterEnergyCent(0)
67         ,fClusterEnergyCentMatched(0)
68         ,fClusterEnergyCentNotMatched(0)
69         ,fClusterEt(0)
70         ,fHistChargedEnergyRemoved(0)
71         ,fHistNeutralEnergyRemoved(0)
72         ,fHistGammaEnergyAdded(0)
73         ,fHistMatchedTracksEvspTvsCent(0)
74         ,fHistMatchedTracksEvspTvsCentEffCorr(0)
75         ,fHistMatchedTracksEvspTvsCentEffTMCorr(0)
76         ,fHistMatchedTracksEvspTvsCentEffTMCorr500MeV(0)
77         ,fHistFoundHadronsvsCent(0)
78         ,fHistNotFoundHadronsvsCent(0)
79         ,fHistFoundHadronsEtvsCent(0)
80         ,fHistNotFoundHadronsEtvsCent(0)
81         ,fHistFoundHadronsvsCent500MeV(0)
82         ,fHistNotFoundHadronsvsCent500MeV(0)
83         ,fHistFoundHadronsEtvsCent500MeV(0)
84         ,fHistNotFoundHadronsEtvsCent500MeV(0)
85         ,fHistNominalRawEt(0)
86         ,fHistNominalNonLinHighEt(0)
87         ,fHistNominalNonLinLowEt(0)
88         ,fHistNominalEffHighEt(0)
89         ,fHistNominalEffLowEt(0)
90         ,fHistTotRawEtEffCorr(0)
91         ,fHistTotRawEt(0)
92         ,fHistTotRawEtEffCorr500MeV(0)
93         ,fHistTotAllRawEt(0)
94         ,fHistTotAllRawEtEffCorr(0)
95         ,fHistNClustersPhosVsEmcal(0)
96         ,fHistClusterSizeVsCent(0)
97         ,fHistMatchedClusterSizeVsCent(0)
98         ,fHistTotAllRawEtVsTotalPt(0)
99         ,fHistTotAllRawEtVsTotalPtVsCent(0)
100         ,fHistTotMatchedRawEtVsTotalPtVsCent(0)
101         ,fHistPIDProtonsTrackMatchedDepositedVsNch(0)
102         ,fHistPIDAntiProtonsTrackMatchedDepositedVsNch(0)
103         ,fHistPIDProtonsTrackMatchedDepositedVsNcl(0)
104         ,fHistPIDAntiProtonsTrackMatchedDepositedVsNcl(0)
105         ,fHistPiKPTrackMatchedDepositedVsNch(0)
106                                                       //,
107         ,fHistPIDProtonsTrackMatchedDepositedVsNchNoEff(0)
108         ,fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff(0)
109         ,fHistPIDProtonsTrackMatchedDepositedVsNclNoEff(0)
110         ,fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff(0)
111         ,fHistPiKPTrackMatchedDepositedVsNchNoEff(0)
112         ,fHistCentVsNchVsNclReco(0)
113         ,fHistRawSignalReco(0)
114         ,fHistEffCorrSignalReco(0)
115 {
116
117 }
118
119 AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
120 {//destructor
121     delete fCorrections;
122     delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */
123     delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */
124     delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */
125     delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */
126     delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
127
128     delete fHistRemovedEnergy; // removed energy
129     delete fClusterPositionAccepted;
130     delete fClusterPositionAll;
131     delete fClusterPositionAcceptedEnergy;
132     delete fClusterPositionAllEnergy;
133     delete fClusterEnergy;
134     delete fClusterEnergyCent;
135     delete fClusterEnergyCentMatched;
136     delete fClusterEnergyCentNotMatched;
137     delete fClusterEt;
138     delete fHistChargedEnergyRemoved;
139     delete fHistNeutralEnergyRemoved;
140     delete fHistGammaEnergyAdded;
141     delete fHistMatchedTracksEvspTvsCent;
142     delete fHistMatchedTracksEvspTvsCentEffCorr;
143     delete fHistMatchedTracksEvspTvsCentEffTMCorr;
144     delete fHistMatchedTracksEvspTvsCentEffTMCorr500MeV;
145     delete fHistFoundHadronsvsCent;
146     delete fHistNotFoundHadronsvsCent;
147     delete fHistFoundHadronsEtvsCent;
148     delete fHistNotFoundHadronsEtvsCent;
149     delete fHistFoundHadronsvsCent500MeV;
150     delete fHistNotFoundHadronsvsCent500MeV;
151     delete fHistFoundHadronsEtvsCent500MeV;
152     delete fHistNotFoundHadronsEtvsCent500MeV;
153     delete fHistNominalRawEt;
154     delete fHistNominalNonLinHighEt;
155     delete fHistNominalNonLinLowEt;
156     delete fHistNominalEffHighEt;
157     delete fHistNominalEffLowEt;
158     delete fHistTotRawEtEffCorr;
159     delete fHistTotRawEt;
160     delete fHistTotAllRawEt;
161     delete fHistTotAllRawEtEffCorr;
162     delete fHistTotRawEtEffCorr500MeV;
163     delete fHistNClustersPhosVsEmcal;
164     delete fHistClusterSizeVsCent;
165     delete fHistMatchedClusterSizeVsCent;
166     delete fHistTotAllRawEtVsTotalPt;
167     delete fHistTotAllRawEtVsTotalPtVsCent;
168     delete fHistTotMatchedRawEtVsTotalPtVsCent;
169     delete fHistPIDProtonsTrackMatchedDepositedVsNch;
170     delete fHistPIDAntiProtonsTrackMatchedDepositedVsNch;
171     delete fHistPIDProtonsTrackMatchedDepositedVsNcl;
172     delete fHistPIDAntiProtonsTrackMatchedDepositedVsNcl;
173     delete fHistPiKPTrackMatchedDepositedVsNch;
174     delete fHistPIDProtonsTrackMatchedDepositedVsNchNoEff;
175     delete fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff;
176     delete fHistPIDProtonsTrackMatchedDepositedVsNclNoEff;
177     delete fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff;
178     delete fHistPiKPTrackMatchedDepositedVsNchNoEff;
179     delete fHistCentVsNchVsNclReco;
180     delete fHistRawSignalReco;
181     delete fHistEffCorrSignalReco;
182 }
183
184 Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
185 {
186
187     //AliAnalysisEt::AnalyseEvent(ev);
188     // analyse ESD event
189     ResetEventValues();
190     if (!ev) {
191         AliFatal("ERROR: Event does not exist");
192         return 0;
193     }
194
195     AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
196     if (!event) {
197         AliFatal("ERROR: ESD Event does not exist");
198         return 0;
199     }
200     if(!fSelector){
201         AliFatal("ERROR: fSelector does not exist");
202         return 0;
203     }
204     fSelector->SetEvent(event);
205     
206     Int_t cent = -1;
207     fCentrality = event->GetCentrality();
208     if (fCentrality && cent)
209     {
210         cent = fCentrality->GetCentralityClass5("V0M");
211         fCentClass = fCentrality->GetCentralityClass5("V0M");
212     }
213
214
215   //for PID
216   AliESDpid *pID = new AliESDpid();
217   pID->MakePID(event);
218   Float_t etPIDProtons = 0.0;
219   Float_t etPIDAntiProtons = 0.0;
220   Float_t etPiKPMatched = 0.0;
221   Float_t etPIDProtonsNoEff = 0.0;
222   Float_t etPIDAntiProtonsNoEff = 0.0;
223   Float_t etPiKPMatchedNoEff = 0.0;
224   Float_t multiplicity = fEsdtrackCutsTPC->GetReferenceMultiplicity(event,kTRUE);
225
226
227     Float_t totalMatchedPt = 0.0;
228     Float_t totalPt = 0.0;
229     TObjArray* list  = fEsdtrackCutsTPC->GetAcceptedTracks(event);
230     Int_t nGoodTracks = list->GetEntries();
231     for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++){
232         AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack));
233         if (!track)
234           {
235             Printf("ERROR: Could not get track %d", iTrack);
236             continue;
237           }
238         else{
239           totalPt +=track->Pt();
240           pID->MakeITSPID(track);
241
242
243         }
244     }
245
246     //TRefArray *caloClusters = fSelector->GetClusters();//just gets the correct set of clusters - does not apply any cuts
247     //Float_t fClusterMult = caloClusters->GetEntries();
248
249     Float_t nominalRawEt = 0;
250     Float_t totEt500MeV = 0;
251     Float_t nonlinHighRawEt = 0;
252     Float_t nonlinLowRawEt = 0;
253     Float_t effHighRawEt = 0;
254     Float_t effLowRawEt = 0;
255     Float_t uncorrEt = 0;
256     Float_t rawSignal;
257     Float_t effCorrSignal;
258
259     nChargedHadronsMeasured = 0.0;
260     nChargedHadronsTotal = 0.0;
261     Float_t nChargedHadronsEtMeasured = 0.0;
262     Float_t nChargedHadronsEtTotal = 0.0;
263     Float_t nChargedHadronsMeasured500MeV = 0.0;
264     Float_t nChargedHadronsTotal500MeV = 0.0;
265     Float_t nChargedHadronsEtMeasured500MeV = 0.0;
266     Float_t nChargedHadronsEtTotal500MeV = 0.0;
267     Float_t fTotAllRawEt = 0.0;
268     Float_t fTotRawEt = 0.0;
269     Float_t fTotRawEtEffCorr = 0.0;
270     Float_t fTotAllRawEtEffCorr = 0.0;
271     Int_t nPhosClusters = 0;
272     Int_t nEmcalClusters = 0;
273
274
275     TRefArray *caloClusters = fSelector->GetClusters();
276     Int_t nCluster = caloClusters->GetEntries();
277
278     for (int iCluster = 0; iCluster < nCluster; iCluster++ )
279     {
280         AliESDCaloCluster* cluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
281         if (!cluster)
282         {
283             AliError(Form("ERROR: Could not get cluster %d", iCluster));
284             continue;
285         }
286         int x = 0;
287         fCutFlow->Fill(x++);
288         if(cluster->IsEMCAL()) nEmcalClusters++;
289         else nPhosClusters++;
290         if(!fSelector->IsDetectorCluster(*cluster)) continue;
291         fCutFlow->Fill(x++);
292         if(!fSelector->PassMinEnergyCut(*cluster)) continue;
293         fCutFlow->Fill(x++);
294         if (!fSelector->PassDistanceToBadChannelCut(*cluster)) continue;
295         fCutFlow->Fill(x++);
296         if (!fSelector->CutGeometricalAcceptance(*cluster)) continue;
297         //fCutFlow->Fill(x++);
298         Float_t pos[3];
299
300         cluster->GetPosition(pos);
301         TVector3 cp(pos);
302         fClusterPositionAll->Fill(cp.Phi(), cp.PseudoRapidity());
303         fClusterPositionAllEnergy->Fill(cp.Phi(), cp.PseudoRapidity(),cluster->E());
304
305         //if(TMath::Abs(cp.Eta())> fCuts->fCuts->GetGeometryEmcalEtaAccCut() || cp.Phi() >  fCuts->GetGeometryEmcalPhiAccMaxCut()*TMath::Pi()/180. ||  cp.Phi() >  fCuts->GetGeometryEmcalPhiAccMinCut()*TMath::Pi()/180.) continue;//Do not accept if cluster is not in the acceptance
306         fTotAllRawEt += TMath::Sin(cp.Theta())*cluster->E();
307         fTotAllRawEtEffCorr += CorrectForReconstructionEfficiency(*cluster,cent);
308
309         fClusterEnergyCent->Fill(cluster->E(),cent);
310
311         Bool_t matched = kTRUE;//default to no track matched
312         Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();//find the index of the matched track
313         matched = !(fSelector->PassTrackMatchingCut(*cluster));//PassTrackMatchingCut is false if there is a matched track
314         if(matched){//if the track match is good (, is the track good?
315           if(trackMatchedIndex < 0) matched=kFALSE;//If the index is bad, don't count it
316           if(matched){
317             AliESDtrack *track = event->GetTrack(trackMatchedIndex);
318             //if this is a good track, accept track will return true.  The track matched is good, so not track matched is false
319             matched = fEsdtrackCutsTPC->AcceptTrack(track);//If the track is bad, don't count it
320           }
321         }
322
323
324         if (matched)
325         {
326           
327             if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0)
328             {
329                 AliVTrack *track = event->GetTrack(trackMatchedIndex);
330                 if (!track) {
331                     AliError("Error: track does not exist");
332                 }
333                 else {
334                   totalMatchedPt +=track->Pt();
335                   fClusterEnergyCentMatched->Fill(cluster->E(),cent);
336                   fHistMatchedClusterSizeVsCent->Fill(cluster->GetNCells(),cent);
337
338                   float eff = fTmCorrections->TrackMatchingEfficiency(track->Pt(),cent);
339                   if(TMath::Abs(eff)<1e-5) eff = 1.0;
340                   //cout<<"pt "<<track->Pt()<<" eff "<<eff<<" total "<<nChargedHadronsTotal<<endl;
341                   nChargedHadronsMeasured++;
342                   nChargedHadronsTotal += 1/eff;
343                   Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,cent);
344                   nChargedHadronsEtMeasured+= TMath::Sin(cp.Theta())*cluster->E();
345                   //One efficiency is the gamma efficiency and the other is the track matching efficiency.
346                   nChargedHadronsEtTotal+= 1/eff *TMath::Sin(cp.Theta())*cluster->E();
347                   //cout<<"nFound "<<1<<" nFoundTotal "<<1/eff<<" etMeas "<<TMath::Sin(cp.Theta())*cluster->E()<<" ET total "<< 1/eff *TMath::Sin(cp.Theta())*cluster->E()<<endl;
348
349                   Float_t nSigmaPion = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion); 
350                   Float_t nSigmaProton = fPIDResponse->NumberOfSigmasTPC(track, AliPID::kProton); 
351                   bool isProton = (nSigmaPion>3.0 && nSigmaProton<3.0 && track->Pt()<0.9);
352                   //cout<<"NSigmaProton "<<nSigmaProton<<endl;
353                   etPiKPMatched += effCorrEt;
354                   etPiKPMatchedNoEff  +=TMath::Sin(cp.Theta())*cluster->E();
355                   if(isProton){
356                     if(track->Charge()>0){
357                       etPIDProtons += effCorrEt;
358                       etPIDProtonsNoEff +=TMath::Sin(cp.Theta())*cluster->E();
359                     }
360                     else{
361                       etPIDAntiProtonsNoEff +=TMath::Sin(cp.Theta())*cluster->E();
362                       etPIDAntiProtons += effCorrEt;
363                     }
364                   }
365                   if(TMath::Sin(cp.Theta())*cluster->E()>0.5){
366                     nChargedHadronsMeasured500MeV++;
367                     nChargedHadronsTotal500MeV += 1/eff;
368                     nChargedHadronsEtMeasured500MeV+= TMath::Sin(cp.Theta())*cluster->E();
369                     nChargedHadronsEtTotal500MeV+= 1/eff *TMath::Sin(cp.Theta())*cluster->E();
370                   }
371                   fHistMatchedTracksEvspTvsCent->Fill(track->P(),TMath::Sin(cp.Theta())*cluster->E(),cent);
372                   fHistMatchedTracksEvspTvsCentEffCorr->Fill(track->P(),effCorrEt,cent);
373                   //Weighed by the number of tracks we didn't find
374                   fHistMatchedTracksEvspTvsCentEffTMCorr->Fill(track->P(), effCorrEt,cent, (1/eff-1) );
375                   cluster->GetPosition(pos);      
376                   TVector3 p2(pos);
377                   uncorrEt += TMath::Sin(p2.Theta())*cluster->E();
378                   if(uncorrEt>=0.5) fHistMatchedTracksEvspTvsCentEffTMCorr500MeV->Fill(track->P(), effCorrEt,cent, (1/eff-1) );
379                     const Double_t *pidWeights = track->PID();
380
381                     Double_t maxpidweight = 0;
382                     Int_t maxpid = 0;
383
384                     if (pidWeights)
385                     {
386                         for (Int_t p =0; p < AliPID::kSPECIES; p++)
387                         {
388                             if (pidWeights[p] > maxpidweight)
389                             {
390                                 maxpidweight = pidWeights[p];
391                                 maxpid = p;
392                             }
393                         }
394                         if (fCuts->GetHistMakeTreeDeposit() && fDepositTree)
395                         {
396                             fEnergyDeposited = cluster->E();
397                             fMomentumTPC = track->P();
398                             fCharge = track->Charge();
399                             fParticlePid = maxpid;
400                             fPidProb = maxpidweight;
401                             AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
402                             if (!esdTrack) {
403                                 AliError("Error: track does not exist");
404                             }
405                             else {
406                                 if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack);
407                                 fDepositTree->Fill();
408                             }
409                         }
410
411                         if (maxpidweight > fPidCut)
412                         {
413                             //Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]);
414
415                             //Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2;
416
417                             //Float_t et = cluster->E() * TMath::Sin(theta);
418                             if (maxpid == AliPID::kProton)
419                             {
420
421                                 if (track->Charge() == 1)
422                                 {
423                                     fHistProtonEnergyDeposit->Fill(cluster->E(), track->E());
424                                 }
425                                 else if (track->Charge() == -1)
426                                 {
427                                     fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E());
428                                 }
429                             }
430                             else if (maxpid == AliPID::kPion)
431                             {
432                                 fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E());
433                             }
434                             else if (maxpid == AliPID::kKaon)
435                             {
436                                 fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E());
437                             }
438                             else if (maxpid == AliPID::kMuon)
439                             {
440                                 fHistMuonEnergyDeposit->Fill(cluster->E(), track->E());
441                             }
442                         }
443                     }
444                 }
445             }
446             //continue;
447         } // distance
448         else{//these are clusters which were not track matched
449           fCutFlow->Fill(x++);
450           //std::cout << x++ << std::endl;
451           
452           //if (cluster->E() >  fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
453           //if (cluster->E() < fClusterEnergyCut) continue;
454           cluster->GetPosition(pos);
455           
456             TVector3 p2(pos);
457             
458             fClusterPositionAccepted->Fill(p2.Phi(), p2.PseudoRapidity());
459             fClusterPositionAcceptedEnergy->Fill(p2.Phi(), p2.PseudoRapidity(),cluster->E());
460             fClusterEnergy->Fill(cluster->E());
461             fClusterEnergyCentNotMatched->Fill(cluster->E(),cent);
462             fHistClusterSizeVsCent->Fill(cluster->GetNCells(),cent);
463             fClusterEt->Fill(TMath::Sin(p2.Theta())*cluster->E());
464             uncorrEt += TMath::Sin(p2.Theta())*cluster->E();
465             float myuncorrEt = TMath::Sin(p2.Theta())*cluster->E();
466             fTotRawEt += myuncorrEt;
467
468             Double_t effCorrEt = CorrectForReconstructionEfficiency(*cluster,cent);
469             rawSignal += myuncorrEt;
470             effCorrSignal +=effCorrEt;
471             //cout<<"cluster energy "<<cluster->E()<<" eff corr Et "<<effCorrEt<<endl;
472             fTotRawEtEffCorr += effCorrEt;
473             fTotNeutralEt += effCorrEt;
474             nominalRawEt += effCorrEt;
475             if(myuncorrEt>=0.5){
476               totEt500MeV += effCorrEt;
477               //cout<<"test "<<myuncorrEt<<"> 0.5"<<endl;
478             }
479             else{
480               //cout<<"test "<<myuncorrEt<<"< 0.5"<<endl;
481             }
482             nonlinHighRawEt += effCorrEt*GetCorrectionModification(*cluster,1,0,cent);
483             nonlinLowRawEt += effCorrEt*GetCorrectionModification(*cluster,-1,0,cent);
484             effHighRawEt += effCorrEt*GetCorrectionModification(*cluster,0,1,cent);
485             effLowRawEt += effCorrEt*GetCorrectionModification(*cluster,0,-1,cent);
486             fNeutralMultiplicity++;
487         }
488         fMultiplicity++;
489     }
490     
491
492     fHistRawSignalReco->Fill(rawSignal);
493     fHistEffCorrSignalReco->Fill(effCorrSignal);
494
495     fHistNClustersPhosVsEmcal->Fill(nPhosClusters,nEmcalClusters,cent);
496     fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
497     fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
498     fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
499     fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
500     
501     fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
502     fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
503
504     Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) + GetGammaContribution(fNeutralMultiplicity) + GetSecondaryContribution(fNeutralMultiplicity);
505     fHistRemovedEnergy->Fill(removedEnergy);
506     
507     fTotNeutralEtAcc = fTotNeutralEt;
508     //fHistTotRawEtEffCorr->Fill(fTotNeutralEt,cent);
509     fHistTotRawEtEffCorr->Fill(fTotRawEtEffCorr,cent);
510     fHistTotRawEt->Fill(fTotRawEt,cent);
511     fHistTotAllRawEt->Fill(fTotAllRawEt,cent);
512     fHistTotAllRawEtVsTotalPt->Fill(fTotAllRawEt,totalPt);
513     fHistTotAllRawEtVsTotalPtVsCent->Fill(fTotAllRawEt,totalPt,cent);
514     fHistTotMatchedRawEtVsTotalPtVsCent->Fill(fTotAllRawEt,totalMatchedPt,cent);
515     fHistTotAllRawEtEffCorr->Fill(fTotAllRawEtEffCorr,cent);
516     //cout<<"fTotAllRawEtEffCorr "<<fTotAllRawEtEffCorr<<" fTotAllRawEt "<<fTotAllRawEt<<" fTotRawEtEffCorr "<<fTotRawEtEffCorr<<"("<<fTotNeutralEt<<")"<<" fTotRawEt "<<fTotRawEt<<endl;
517     //cout<<"uncorr "<<uncorrEt<<" raw "<<nominalRawEt<<" tot raw "<<fTotNeutralEt;
518     fTotNeutralEt = fGeomCorrection * fEMinCorrection * (fTotNeutralEt - removedEnergy);
519     //cout<<" tot corr "<<fTotNeutralEt<<endl;
520     fTotEt = fTotChargedEt + fTotNeutralEt;
521 // Fill the histograms...0
522     FillHistograms();
523     //std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
524     //cout<<"cent "<<cent<<" cluster mult "<<fClusterMult<<" fTotNeutralEt "<<fTotNeutralEt<<" nominalRawEt "<<nominalRawEt<<endl;
525     fHistNominalRawEt->Fill(nominalRawEt,cent);
526     fHistTotRawEtEffCorr500MeV->Fill(totEt500MeV,cent);
527     fHistNominalNonLinHighEt->Fill(nonlinHighRawEt,cent);
528     fHistNominalNonLinLowEt->Fill(nonlinLowRawEt,cent);
529     fHistNominalEffHighEt->Fill(effHighRawEt,cent);
530     fHistNominalEffLowEt->Fill(effLowRawEt,cent);
531     fHistFoundHadronsvsCent->Fill(nChargedHadronsMeasured,cent);
532     fHistNotFoundHadronsvsCent->Fill(nChargedHadronsTotal-nChargedHadronsMeasured,cent);
533     fHistFoundHadronsEtvsCent->Fill(nChargedHadronsEtMeasured,cent);
534     fHistNotFoundHadronsEtvsCent->Fill(nChargedHadronsEtTotal-nChargedHadronsEtMeasured,cent);
535     //cout<<"found "<<nChargedHadronsMeasured<<" total "<<nChargedHadronsTotal<<" not found "<<nChargedHadronsTotal-nChargedHadronsMeasured<<" found "<< nChargedHadronsMeasured500MeV<<" not found "<< nChargedHadronsTotal500MeV-nChargedHadronsMeasured500MeV <<" total "<<nChargedHadronsTotal500MeV<<endl;
536     fHistFoundHadronsvsCent500MeV->Fill(nChargedHadronsMeasured500MeV,cent);
537     fHistNotFoundHadronsvsCent500MeV->Fill(nChargedHadronsTotal500MeV-nChargedHadronsMeasured500MeV,cent);
538     fHistFoundHadronsEtvsCent500MeV->Fill(nChargedHadronsEtMeasured500MeV,cent);
539     fHistNotFoundHadronsEtvsCent500MeV->Fill(nChargedHadronsEtTotal500MeV-nChargedHadronsEtMeasured500MeV,cent);
540 //     cout<<"Number of hadrons measured:  "<<nChargedHadronsMeasured<<" Estimated total number of hadrons "<<nChargedHadronsTotal<<" ET in track matched hadrons "<<
541 //       nChargedHadronsEtMeasured;
542 //     if(nChargedHadronsMeasured>0)cout<<" ("<<nChargedHadronsEtMeasured/nChargedHadronsMeasured<<") ";
543 //     cout<<" ET in all hadrons ";
544 //     cout<<nChargedHadronsEtTotal;
545 //     if(nChargedHadronsTotal>0) cout<<" ("<<nChargedHadronsEtTotal/nChargedHadronsTotal<<") ";
546 //     cout<<endl;
547     fHistPIDProtonsTrackMatchedDepositedVsNch->Fill(etPIDProtons,multiplicity);
548     fHistPIDAntiProtonsTrackMatchedDepositedVsNch->Fill(etPIDAntiProtons,multiplicity);
549     fHistPIDProtonsTrackMatchedDepositedVsNcl->Fill(etPIDProtons,nCluster);
550     fHistPIDAntiProtonsTrackMatchedDepositedVsNcl->Fill(etPIDAntiProtons,nCluster);
551     fHistPIDProtonsTrackMatchedDepositedVsNchNoEff->Fill(etPIDProtonsNoEff,multiplicity);
552     fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff->Fill(etPIDAntiProtonsNoEff,multiplicity);
553     fHistPIDProtonsTrackMatchedDepositedVsNclNoEff->Fill(etPIDProtonsNoEff,nCluster);
554     fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff->Fill(etPIDAntiProtonsNoEff,nCluster);
555     fHistCentVsNchVsNclReco->Fill(cent,multiplicity,nCluster);
556     fHistPiKPTrackMatchedDepositedVsNch->Fill(etPiKPMatched,multiplicity);
557     fHistPiKPTrackMatchedDepositedVsNchNoEff->Fill(etPiKPMatchedNoEff,multiplicity);
558     delete pID;
559     return 0;
560 }
561
562 bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
563 { // check vertex
564
565     Float_t bxy = 999.;
566     Float_t bz = 999.;
567     if (!track) {
568         AliError("ERROR: no track");
569         return kFALSE;
570     }
571     AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
572     if (!esdTrack) {
573         AliError("ERROR: no track");
574         return kFALSE;
575     }
576     esdTrack->GetImpactParametersTPC(bxy,bz);
577
578
579     bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
580                   (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
581                   (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
582                   (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
583                   (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
584
585     return status;
586 }
587
588 void AliAnalysisEtReconstructed::Init()
589 { // Init
590     AliAnalysisEt::Init();
591     fPidCut = fCuts->GetReconstructedPidCut();
592     TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
593     if (!fCorrections) {
594         cout<<"Warning!  You have not set corrections.  Your code will crash.  You have to set the corrections."<<endl;
595     }
596 }
597
598 bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
599 { // propagate track to detector radius
600     if (!track) {
601         cout<<"Warning: track empty"<<endl;
602         return kFALSE;
603     }
604     AliESDtrack *esdTrack= dynamic_cast<AliESDtrack*>(track);
605     if (!esdTrack) {
606         AliError("ERROR: no ESD track");
607         return kFALSE;
608     }
609     // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
610
611     Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
612
613     // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
614     return prop && fSelector->CutGeometricalAcceptance(*esdTrack);
615 }
616
617 void AliAnalysisEtReconstructed::FillOutputList(TList* list)
618 { // add some extra histograms to the ones from base class
619     AliAnalysisEt::FillOutputList(list);
620
621     list->Add(fHistChargedPionEnergyDeposit);
622     list->Add(fHistProtonEnergyDeposit);
623     list->Add(fHistAntiProtonEnergyDeposit);
624     list->Add(fHistChargedKaonEnergyDeposit);
625     list->Add(fHistMuonEnergyDeposit);
626
627     list->Add(fHistRemovedEnergy);
628     list->Add(fClusterPositionAccepted);
629     list->Add(fClusterPositionAll);
630     list->Add(fClusterPositionAcceptedEnergy);
631     list->Add(fClusterPositionAllEnergy);
632     list->Add(fClusterEnergy);
633     list->Add(fClusterEnergyCent);
634     list->Add(fClusterEnergyCentMatched);
635     list->Add(fClusterEnergyCentNotMatched);
636     list->Add(fClusterEt);
637     
638     list->Add(fHistChargedEnergyRemoved);
639     list->Add(fHistNeutralEnergyRemoved);
640     list->Add(fHistGammaEnergyAdded);
641     list->Add(fHistMatchedTracksEvspTvsCent);
642     list->Add(fHistMatchedTracksEvspTvsCentEffCorr);
643     list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr);
644     list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr500MeV);
645     list->Add(fHistFoundHadronsvsCent);
646     list->Add(fHistNotFoundHadronsvsCent);
647     list->Add(fHistFoundHadronsEtvsCent);
648     list->Add(fHistNotFoundHadronsEtvsCent);
649     list->Add(fHistFoundHadronsvsCent500MeV);
650     list->Add(fHistNotFoundHadronsvsCent500MeV);
651     list->Add(fHistFoundHadronsEtvsCent500MeV);
652     list->Add(fHistNotFoundHadronsEtvsCent500MeV);
653     list->Add(fHistNominalRawEt);
654     list->Add(fHistNominalNonLinHighEt);
655     list->Add(fHistNominalNonLinLowEt);
656     list->Add(fHistNominalEffHighEt);
657     list->Add(fHistNominalEffLowEt);
658     list->Add(fHistTotRawEtEffCorr);
659     list->Add(fHistTotRawEtEffCorr500MeV);
660     list->Add(fHistTotAllRawEtEffCorr);
661     list->Add(fHistTotRawEt);
662     list->Add(fHistTotAllRawEt);
663     list->Add(fHistNClustersPhosVsEmcal);
664     list->Add(fHistClusterSizeVsCent);
665     list->Add(fHistMatchedClusterSizeVsCent);
666     list->Add(fHistTotAllRawEtVsTotalPt);
667     list->Add(fHistTotAllRawEtVsTotalPtVsCent);
668     list->Add(fHistTotMatchedRawEtVsTotalPtVsCent);
669     list->Add(fHistPIDProtonsTrackMatchedDepositedVsNch);
670     list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNch);
671     list->Add(fHistPIDProtonsTrackMatchedDepositedVsNcl);
672     list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNcl);
673     list->Add(fHistPiKPTrackMatchedDepositedVsNch);
674     list->Add(fHistPIDProtonsTrackMatchedDepositedVsNchNoEff);
675     list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff);
676     list->Add(fHistPIDProtonsTrackMatchedDepositedVsNclNoEff);
677     list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff);
678     list->Add(fHistPiKPTrackMatchedDepositedVsNchNoEff);
679     list->Add(fHistCentVsNchVsNclReco);
680     list->Add(fHistRawSignalReco);
681     list->Add(fHistEffCorrSignalReco);
682 }
683
684 void AliAnalysisEtReconstructed::CreateHistograms()
685 { // add some extra histograms to the ones from base class
686     AliAnalysisEt::CreateHistograms();
687
688     Int_t nbinsEt = 1000;
689     Double_t minEt = 0;
690     Double_t maxEt = 10;
691
692     // possibly change histogram limits
693 //     if (fCuts) {
694 //         nbinsEt = fCuts->GetHistNbinsParticleEt();
695 //         minEt = fCuts->GetHistMinParticleEt();
696 //         maxEt = fCuts->GetHistMaxParticleEt();
697 //     }
698
699     TString histname;
700     histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
701     fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
702     fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
703     fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
704
705     histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
706     fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
707     fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
708     fHistProtonEnergyDeposit->SetYTitle("Energy of track");
709
710     histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
711     fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
712     fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
713     fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
714
715     histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
716     fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
717     fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
718     fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
719
720     histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
721     fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
722     fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
723     fHistMuonEnergyDeposit->SetYTitle("Energy of track");
724
725     histname = "fHistRemovedEnergy" + fHistogramNameSuffix;
726     fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20);
727     //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
728     //fHistMuonEnergyDeposit->SetYTitle("Energy of track");
729
730     histname = "fClusterPositionAccepted" + fHistogramNameSuffix;
731     fClusterPositionAccepted = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
732     fClusterPositionAccepted->SetXTitle("#phi");
733     fClusterPositionAccepted->SetYTitle("#eta");
734
735     histname = "fClusterPositionAll" + fHistogramNameSuffix;
736     fClusterPositionAll = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
737     fClusterPositionAll->SetXTitle("#phi");
738     fClusterPositionAll->SetYTitle("#eta");
739
740     histname = "fClusterPositionAcceptedEnergy" + fHistogramNameSuffix;
741     fClusterPositionAcceptedEnergy = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
742     fClusterPositionAcceptedEnergy->SetXTitle("#phi");
743     fClusterPositionAcceptedEnergy->SetYTitle("#eta");
744
745     histname = "fClusterPositionAllEnergy" + fHistogramNameSuffix;
746     fClusterPositionAllEnergy = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
747     fClusterPositionAllEnergy->SetXTitle("#phi");
748     fClusterPositionAllEnergy->SetYTitle("#eta");
749
750     histname = "fClusterEnergy" + fHistogramNameSuffix;
751     fClusterEnergy = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
752     fClusterEnergy->SetYTitle("Number of clusters");
753     fClusterEnergy->SetXTitle("Energy of cluster");
754
755     histname = "fClusterEnergyCent" + fHistogramNameSuffix;
756     fClusterEnergyCent = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
757     fClusterEnergyCent->SetXTitle("Energy of cluster");
758     fClusterEnergyCent->SetYTitle("Centrality Bin");
759     fClusterEnergyCent->SetZTitle("Number of clusters");
760
761     histname = "fClusterEnergyCentMatched" + fHistogramNameSuffix;
762     fClusterEnergyCentMatched = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
763     fClusterEnergyCentMatched->SetXTitle("Energy of cluster");
764     fClusterEnergyCentMatched->SetYTitle("Centrality Bin");
765     fClusterEnergyCentMatched->SetZTitle("Number of Clusters");
766
767     histname = "fClusterEnergyCentNotMatched" + fHistogramNameSuffix;
768     fClusterEnergyCentNotMatched = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
769     fClusterEnergyCentNotMatched->SetXTitle("Energy of cluster");
770     fClusterEnergyCentNotMatched->SetYTitle("Centrality Bin");
771     fClusterEnergyCentNotMatched->SetZTitle("Number of clusters");
772
773     histname = "fClusterEt" + fHistogramNameSuffix;
774     fClusterEt = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
775     fClusterEt->SetXTitle("Number of clusters");
776     fClusterEt->SetYTitle("E_{T} of cluster");
777
778     histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix;
779     fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
780
781     histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix;
782     fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
783
784     histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix;
785     fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
786
787     fHistMatchedTracksEvspTvsCent = new TH3F("fHistMatchedTracksEvspTvsCent", "fHistMatchedTracksEvspTvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
788     fHistMatchedTracksEvspTvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffCorr", "fHistMatchedTracksEvspTvsCentEffCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
789     fHistMatchedTracksEvspTvsCentEffTMCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr", "fHistMatchedTracksEvspTvsCentEffTMCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
790     fHistMatchedTracksEvspTvsCentEffTMCorr500MeV = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr500MeV", "fHistMatchedTracksEvspTvsCentEffTMCorr500MeV",100, 0, 3,100,0,3,20,-0.5,19.5);
791
792     float max = 200;
793     if(fHistogramNameSuffix.Contains("P")){max = 100;}
794     fHistFoundHadronsvsCent = new TH2F("fHistFoundHadronsvsCent","fHistFoundHadronsvsCent",100,0,max,20,-0.5,19.5);
795        fHistNotFoundHadronsvsCent = new TH2F("fHistNotFoundHadronsvsCent","fHistNotFoundHadronsvsCent",100,0,max,20,-0.5,19.5);
796     fHistFoundHadronsEtvsCent = new TH2F("fHistFoundHadronsEtvsCent","fHistFoundHadronsEtvsCent",100,0,max,20,-0.5,19.5);
797        fHistNotFoundHadronsEtvsCent = new TH2F("fHistNotFoundHadronsEtvsCent","fHistNotFoundHadronsEtvsCent",100,0,max,20,-0.5,19.5);
798     fHistFoundHadronsvsCent500MeV = new TH2F("fHistFoundHadronsvsCent500MeV","fHistFoundHadronsvsCent500MeV",100,0,max,20,-0.5,19.5);
799     fHistNotFoundHadronsvsCent500MeV = new TH2F("fHistNotFoundHadronsvsCent500MeV","fHistNotFoundHadronsvsCent500MeV",100,0,max,20,-0.5,19.5);
800     fHistFoundHadronsEtvsCent500MeV = new TH2F("fHistFoundHadronsEtvsCent500MeV","fHistFoundHadronsEtvsCent500MeV",100,0,max,20,-0.5,19.5);
801     fHistNotFoundHadronsEtvsCent500MeV = new TH2F("fHistNotFoundHadronsEtvsCent500MeV","fHistNotFoundHadronsEtvsCent500MeV",100,0,max,20,-0.5,19.5);
802
803     fHistTotRawEtEffCorr = new TH2F("fHistTotRawEtEffCorr","fHistTotRawEtEffCorr",250,0,250,20,-0.5,19.5);
804     fHistTotRawEt = new TH2F("fHistTotRawEt","fHistTotRawEt",250,0,250,20,-0.5,19.5);
805     fHistTotRawEtEffCorr500MeV = new TH2F("fHistTotRawEtEffCorr500MeV","fHistTotRawEtEffCorr500MeV",250,0,250,20,-0.5,19.5);
806     fHistTotAllRawEt = new TH2F("fHistTotAllRawEt","fHistTotAllRawEt",250,0,250,20,-0.5,19.5);
807     fHistTotAllRawEtEffCorr = new TH2F("fHistTotAllRawEtEffCorr","fHistTotAllRawEtEffCorr",250,0,250,20,-0.5,19.5);
808     fHistNClustersPhosVsEmcal = new TH3F("fHistNClustersPhosVsEmcal","fHistNClustersPhosVsEmcal",50,0,50,250,0,250,20,-0.5,19);
809     fHistClusterSizeVsCent = new TH2F("fHistClusterSizeVsCent","fHistClusterSizeVsCent",10,0.5,10.5,20,-0.5,19.5);
810     fHistMatchedClusterSizeVsCent = new TH2F("fHistMatchedClusterSizeVsCent","fHistMatchedClusterSizeVsCent",10,0.5,10.5,20,-0.5,19.5);
811     fHistTotAllRawEtVsTotalPt = new TH2F("fHistTotAllRawEtVsTotalPt","fHistTotAllRawEtVsTotalPt",125,0,250,200,0,2000);
812     fHistTotAllRawEtVsTotalPtVsCent = new TH3F("fHistTotAllRawEtVsTotalPtVsCent","fHistTotAllRawEtVsTotalPtVsCent",125,0,250,200,0,2000,20,-0.5,19.5);
813     fHistTotMatchedRawEtVsTotalPtVsCent = new TH3F("fHistTotMatchedRawEtVsTotalPtVsCent","fHistTotMatchedRawEtVsTotalPtVsCent",250,0,250,100,0,200,20,-0.5,19.5);
814     
815     maxEt = 500;
816     histname = "fHistNominalRawEt" + fHistogramNameSuffix;
817     fHistNominalRawEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
818     histname = "fHistNominalNonLinHighEt" + fHistogramNameSuffix;
819     fHistNominalNonLinHighEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
820     histname = "fHistNominalNonLinLowEt" + fHistogramNameSuffix;
821     fHistNominalNonLinLowEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
822     histname = "fHistNominalEffHighEt" + fHistogramNameSuffix;
823     fHistNominalEffHighEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
824     histname = "fHistNominalEffLowEt" + fHistogramNameSuffix;
825     fHistNominalEffLowEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
826
827       Float_t maxEtRange = 25;
828       Float_t maxEtRangeHigh = 125;
829       Float_t minEtRange = 0;
830       Int_t nbinsMult = 100;
831       Float_t maxMult = 3000;
832       Float_t minMult = 0;
833       Int_t nbinsCl = 250;
834       Float_t maxCl = 500;
835       Float_t minCl = 0;
836     fHistPIDProtonsTrackMatchedDepositedVsNch = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNch","PID'd protons deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
837     fHistPIDAntiProtonsTrackMatchedDepositedVsNch = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNch","PID'd #bar{p} E_{T} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
838     fHistPIDProtonsTrackMatchedDepositedVsNcl = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNcl","PID'd protons deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
839     fHistPIDAntiProtonsTrackMatchedDepositedVsNcl = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNcl","PID'd #bar{p} E_{T} deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
840     fHistPiKPTrackMatchedDepositedVsNch = new TH2F("fHistPiKPTrackMatchedDepositedVsNch","PiKP track matched",nbinsEt,minEtRange,maxEtRangeHigh,nbinsMult,minMult,maxMult);
841
842     fHistPIDProtonsTrackMatchedDepositedVsNchNoEff = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNchNoEff","PID'd protons deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
843     fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff","PID'd #bar{p} E_{T} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
844     fHistPIDProtonsTrackMatchedDepositedVsNclNoEff = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNclNoEff","PID'd protons deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
845     fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff","PID'd #bar{p} E_{T} deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
846     fHistPiKPTrackMatchedDepositedVsNchNoEff = new TH2F("fHistPiKPTrackMatchedDepositedVsNchNoEff","PiKP track matched",nbinsEt,minEtRange,maxEtRangeHigh,nbinsMult,minMult,maxMult);
847
848
849     fHistCentVsNchVsNclReco = new TH3F("fHistCentVsNchVsNclReco","Cent bin vs Nch Vs NCl",20,-0.5,19.5,nbinsMult,minMult,maxMult,nbinsCl,minCl,maxCl);
850
851    fHistRawSignalReco = new TH1F("fHistRawSignalReco","fHistRawSignalReco",20,-0.5,19.5);
852    fHistEffCorrSignalReco = new TH1F("fHistEffCorrSignalReco","fHistEffCorrSignalReco",20,-0.5,19.5);
853
854 }
855 Double_t AliAnalysisEtReconstructed::ApplyModifiedCorrections(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t cent)
856 {
857   Float_t pos[3];
858   cluster.GetPosition(pos);
859   TVector3 cp(pos);
860   Double_t corrEnergy = fReCorrections->CorrectedEnergy(cluster.E(),cent);
861   
862   Double_t factorNonLin = GetCorrectionModification(cluster, nonLinCorr,effCorr,cent);
863
864     cout<<"Warning:  This function should not get called!"<<endl;
865   //std::cout << "Original energy: " << cluster.E() << ", corrected energy: " << corrEnergy << std::endl;
866   return TMath::Sin(cp.Theta())*corrEnergy*factorNonLin;
867 }
868
869 Double_t AliAnalysisEtReconstructed::GetCorrectionModification(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t cent){//nonLinCorr 0 = nominal 1 = high -1 = low, effCorr  0 = nominal 1 = high -1 = low
870   if(nonLinCorr==0){
871     cout<<"Warning:  This function should not get called!"<<endl;//this statement is basically here to avoid a compilation warning
872   }
873   if(effCorr==0){
874     cout<<"Warning:  This function should not get called!"<<endl;//this statement is basically here to avoid a compilation warning
875   }
876   return cluster.E()*cent;
877 }