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