]> git.uio.no Git - u/mrichter/AliRoot.git/blame_incremental - PWGLF/totEt/AliAnalysisEtReconstructed.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame_incremental
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (non-incremental) | history | HEAD
Fixing Oysteins old cross check histograms
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEtReconstructed.cxx
... / ...
CommitLineData
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
41using namespace std;
42
43ClassImp(AliAnalysisEtReconstructed);
44
45
46AliAnalysisEtReconstructed::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
122AliAnalysisEtReconstructed::~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
190Int_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(cent) + GetNeutralContribution(cent) + GetGammaContribution(cent) + GetSecondaryContribution(cent);//fNeutralMultiplicity
557 Double_t removedEnergy = GetHadronContribution(cent) + GetNeutronContribution(cent) + GetKaonContribution(cent) + GetSecondaryContribution(cent);//fNeutralMultiplicity
558 //cout<<" centbin "<<cent<<" removed energy ch "<< GetHadronContribution(cent)<<" n " << GetNeutronContribution(cent) <<" kaon "<< GetKaonContribution(cent)<<" secondary "<< GetSecondaryContribution(cent);//fNeutralMultiplicity
559 //cout<<" test min et "<<fTmCorrections->GetMinEtCorrection(cent);
560 //cout<<" test neutral "<<fTmCorrections->NeutralContr(cent);
561 //cout<<" test neutron "<<fTmCorrections->NeutralContr(cent);
562 fHistRemovedEnergy->Fill(removedEnergy);
563
564 fTotNeutralEtAcc = fTotNeutralEt;
565 //fHistTotRawEtEffCorr->Fill(fTotNeutralEt,cent);
566 fHistTotRawEtEffCorr->Fill(fTotRawEtEffCorr,cent);
567 fHistTotRawEt->Fill(fTotRawEt,cent);
568 fHistTotAllRawEt->Fill(fTotAllRawEt,cent);
569 fHistTotAllRawEtVsTotalPt->Fill(fTotAllRawEt,totalPt);
570 fHistTotAllRawEtVsTotalPtVsCent->Fill(fTotAllRawEt,totalPt,cent);
571 fHistTotMatchedRawEtVsTotalPtVsCent->Fill(fTotAllRawEt,totalMatchedPt,cent);
572 fHistTotAllRawEtEffCorr->Fill(fTotAllRawEtEffCorr,cent);
573 //cout<<"fTotAllRawEtEffCorr "<<fTotAllRawEtEffCorr<<" fTotAllRawEt "<<fTotAllRawEt<<" fTotRawEtEffCorr "<<fTotRawEtEffCorr<<"("<<fTotNeutralEt<<")"<<" fTotRawEt "<<fTotRawEt<<endl;
574 //cout<<"uncorr "<<uncorrEt<<" raw "<<nominalRawEt<<" tot raw "<<fTotNeutralEt;
575 //cout<<" raw "<<fTotNeutralEt<<" removed "<<removedEnergy<<" etmin "<<GetMinEtCorrection(cent)<<" final ";
576 if(GetMinEtCorrection(cent)>0) fTotNeutralEt = (fTotNeutralEt - removedEnergy)/GetMinEtCorrection(cent);
577 //cout<<fTotNeutralEt<<endl;
578 //cout<<" tot corr "<<fTotNeutralEt<<endl;
579 fTotEt = fTotChargedEt + fTotNeutralEt;
580// Fill the histograms...0
581 FillHistograms();
582 //std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
583 //cout<<"cent "<<cent<<" cluster mult "<<fClusterMult<<" fTotNeutralEt "<<fTotNeutralEt<<" nominalRawEt "<<nominalRawEt<<endl;
584 fHistNominalRawEt->Fill(nominalRawEt,cent);
585 fHistTotRawEtEffCorr500MeV->Fill(totEt500MeV,cent);
586 fHistNominalNonLinHighEt->Fill(nonlinHighRawEt,cent);
587 fHistNominalNonLinLowEt->Fill(nonlinLowRawEt,cent);
588 fHistNominalEffHighEt->Fill(effHighRawEt,cent);
589 fHistNominalEffLowEt->Fill(effLowRawEt,cent);
590 fHistFoundHadronsvsCent->Fill(nChargedHadronsMeasured,cent);
591 fHistNotFoundHadronsvsCent->Fill(nChargedHadronsTotal-nChargedHadronsMeasured,cent);
592 fHistFoundHadronsEtvsCent->Fill(nChargedHadronsEtMeasured,cent);
593 fHistNotFoundHadronsEtvsCent->Fill(nChargedHadronsEtTotal-nChargedHadronsEtMeasured,cent);
594 //cout<<"found "<<nChargedHadronsMeasured<<" total "<<nChargedHadronsTotal<<" not found "<<nChargedHadronsTotal-nChargedHadronsMeasured<<" found "<< nChargedHadronsMeasured500MeV<<" not found "<< nChargedHadronsTotal500MeV-nChargedHadronsMeasured500MeV <<" total "<<nChargedHadronsTotal500MeV<<endl;
595 fHistFoundHadronsvsCent500MeV->Fill(nChargedHadronsMeasured500MeV,cent);
596 fHistNotFoundHadronsvsCent500MeV->Fill(nChargedHadronsTotal500MeV-nChargedHadronsMeasured500MeV,cent);
597 fHistFoundHadronsEtvsCent500MeV->Fill(nChargedHadronsEtMeasured500MeV,cent);
598 fHistNotFoundHadronsEtvsCent500MeV->Fill(nChargedHadronsEtTotal500MeV-nChargedHadronsEtMeasured500MeV,cent);
599// cout<<"Number of hadrons measured: "<<nChargedHadronsMeasured<<" Estimated total number of hadrons "<<nChargedHadronsTotal<<" ET in track matched hadrons "<<
600// nChargedHadronsEtMeasured;
601// if(nChargedHadronsMeasured>0)cout<<" ("<<nChargedHadronsEtMeasured/nChargedHadronsMeasured<<") ";
602// cout<<" ET in all hadrons ";
603// cout<<nChargedHadronsEtTotal;
604// if(nChargedHadronsTotal>0) cout<<" ("<<nChargedHadronsEtTotal/nChargedHadronsTotal<<") ";
605// cout<<endl;
606 fHistPIDProtonsTrackMatchedDepositedVsNch->Fill(etPIDProtons,multiplicity);
607 fHistPIDAntiProtonsTrackMatchedDepositedVsNch->Fill(etPIDAntiProtons,multiplicity);
608 fHistPIDProtonsTrackMatchedDepositedVsNcl->Fill(etPIDProtons,nCluster);
609 fHistPIDAntiProtonsTrackMatchedDepositedVsNcl->Fill(etPIDAntiProtons,nCluster);
610 fHistPIDProtonsTrackMatchedDepositedVsNchNoEff->Fill(etPIDProtonsNoEff,multiplicity);
611 fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff->Fill(etPIDAntiProtonsNoEff,multiplicity);
612 fHistPIDProtonsTrackMatchedDepositedVsNclNoEff->Fill(etPIDProtonsNoEff,nCluster);
613 fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff->Fill(etPIDAntiProtonsNoEff,nCluster);
614 fHistCentVsNchVsNclReco->Fill(cent,multiplicity,nCluster);
615 fHistPiKPTrackMatchedDepositedVsNch->Fill(etPiKPMatched,multiplicity);
616 fHistPiKPTrackMatchedDepositedVsNchNoEff->Fill(etPiKPMatchedNoEff,multiplicity);
617 //delete pID;
618 return 0;
619}
620
621bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
622{ // check vertex
623
624 Float_t bxy = 999.;
625 Float_t bz = 999.;
626 if (!track) {
627 AliError("ERROR: no track");
628 return kFALSE;
629 }
630 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
631 if (!esdTrack) {
632 AliError("ERROR: no track");
633 return kFALSE;
634 }
635 esdTrack->GetImpactParametersTPC(bxy,bz);
636
637
638 bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) &&
639 (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) &&
640 (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) &&
641 (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) &&
642 (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut());
643
644 return status;
645}
646
647void AliAnalysisEtReconstructed::Init()
648{ // Init
649 AliAnalysisEt::Init();
650 fPidCut = fCuts->GetReconstructedPidCut();
651 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
652 if (!fCorrections) {
653 cout<<"Warning! You have not set corrections. Your code will crash. You have to set the corrections."<<endl;
654 }
655}
656
657bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
658{ // propagate track to detector radius
659 if (!track) {
660 cout<<"Warning: track empty"<<endl;
661 return kFALSE;
662 }
663 AliESDtrack *esdTrack= dynamic_cast<AliESDtrack*>(track);
664 if (!esdTrack) {
665 AliError("ERROR: no ESD track");
666 return kFALSE;
667 }
668 // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
669
670 Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
671
672 // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
673 return prop && fSelector->CutGeometricalAcceptance(*esdTrack);
674}
675
676void AliAnalysisEtReconstructed::FillOutputList(TList* list)
677{ // add some extra histograms to the ones from base class
678 AliAnalysisEt::FillOutputList(list);
679
680 list->Add(fHistChargedPionEnergyDeposit);
681 list->Add(fHistProtonEnergyDeposit);
682 list->Add(fHistAntiProtonEnergyDeposit);
683 list->Add(fHistChargedKaonEnergyDeposit);
684 list->Add(fHistMuonEnergyDeposit);
685
686 list->Add(fHistRemovedEnergy);
687 list->Add(fClusterPositionAccepted);
688 list->Add(fClusterPositionAll);
689 list->Add(fClusterPositionAcceptedEnergy);
690 list->Add(fClusterPositionAllEnergy);
691 list->Add(fClusterEnergy);
692 list->Add(fClusterEnergyCent);
693 list->Add(fClusterEnergyModifiedTrackMatchesCent);
694 list->Add(fClusterEnergyCentMatched);
695 list->Add(fClusterEnergyCentNotMatched);
696 list->Add(fClusterEt);
697
698 list->Add(fHistChargedEnergyRemoved);
699 list->Add(fHistNeutralEnergyRemoved);
700 list->Add(fHistGammaEnergyAdded);
701 list->Add(fHistMatchedTracksEvspTvsCent);
702 list->Add(fHistMatchedTracksEvspTvsCentEffCorr);
703 list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr);
704 list->Add(fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr);
705 list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr500MeV);
706 list->Add(fHistFoundHadronsvsCent);
707 list->Add(fHistNotFoundHadronsvsCent);
708 list->Add(fHistFoundHadronsEtvsCent);
709 list->Add(fHistNotFoundHadronsEtvsCent);
710 list->Add(fHistFoundHadronsvsCent500MeV);
711 list->Add(fHistNotFoundHadronsvsCent500MeV);
712 list->Add(fHistFoundHadronsEtvsCent500MeV);
713 list->Add(fHistNotFoundHadronsEtvsCent500MeV);
714 list->Add(fHistNominalRawEt);
715 list->Add(fHistNominalNonLinHighEt);
716 list->Add(fHistNominalNonLinLowEt);
717 list->Add(fHistNominalEffHighEt);
718 list->Add(fHistNominalEffLowEt);
719 list->Add(fHistTotRawEtEffCorr);
720 list->Add(fHistTotRawEtEffCorr500MeV);
721 list->Add(fHistTotAllRawEtEffCorr);
722 list->Add(fHistTotRawEt);
723 list->Add(fHistTotAllRawEt);
724 list->Add(fHistNClustersPhosVsEmcal);
725 list->Add(fHistClusterSizeVsCent);
726 list->Add(fHistMatchedClusterSizeVsCent);
727 list->Add(fHistTotAllRawEtVsTotalPt);
728 list->Add(fHistTotAllRawEtVsTotalPtVsCent);
729 list->Add(fHistTotMatchedRawEtVsTotalPtVsCent);
730 list->Add(fHistPIDProtonsTrackMatchedDepositedVsNch);
731 list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNch);
732 list->Add(fHistPIDProtonsTrackMatchedDepositedVsNcl);
733 list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNcl);
734 list->Add(fHistPiKPTrackMatchedDepositedVsNch);
735 list->Add(fHistPIDProtonsTrackMatchedDepositedVsNchNoEff);
736 list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff);
737 list->Add(fHistPIDProtonsTrackMatchedDepositedVsNclNoEff);
738 list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff);
739 list->Add(fHistPiKPTrackMatchedDepositedVsNchNoEff);
740 list->Add(fHistCentVsNchVsNclReco);
741 list->Add(fHistRawSignalReco);
742 list->Add(fHistEffCorrSignalReco);
743 list->Add(fHistRecoRCorrVsPtVsCent);
744}
745
746void AliAnalysisEtReconstructed::CreateHistograms()
747{ // add some extra histograms to the ones from base class
748 AliAnalysisEt::CreateHistograms();
749
750 Float_t scale = 1;//scale up histograms if EMCal 2011 so we have the right range
751 if(fDataSet==2011 && !fHistogramNameSuffix.Contains("P")){
752 scale = 2.5;
753 }
754 Int_t nbinsEt = 1000*scale;
755 Double_t minEt = 0;
756 Double_t maxEt = 10*scale;
757
758 // possibly change histogram limits
759// if (fCuts) {
760// nbinsEt = fCuts->GetHistNbinsParticleEt();
761// minEt = fCuts->GetHistMinParticleEt();
762// maxEt = fCuts->GetHistMaxParticleEt();
763// }
764
765 TString histname;
766 histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
767 fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
768 fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
769 fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
770
771 histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
772 fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
773 fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
774 fHistProtonEnergyDeposit->SetYTitle("Energy of track");
775
776 histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
777 fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
778 fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
779 fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
780
781 histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
782 fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
783 fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
784 fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
785
786 histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
787 fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
788 fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
789 fHistMuonEnergyDeposit->SetYTitle("Energy of track");
790
791 histname = "fHistRemovedEnergy" + fHistogramNameSuffix;
792 fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20);
793 //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
794 //fHistMuonEnergyDeposit->SetYTitle("Energy of track");
795
796 histname = "fClusterPositionAccepted" + fHistogramNameSuffix;
797 fClusterPositionAccepted = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
798 fClusterPositionAccepted->SetXTitle("#phi");
799 fClusterPositionAccepted->SetYTitle("#eta");
800
801 histname = "fClusterPositionAll" + fHistogramNameSuffix;
802 fClusterPositionAll = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
803 fClusterPositionAll->SetXTitle("#phi");
804 fClusterPositionAll->SetYTitle("#eta");
805
806 histname = "fClusterPositionAcceptedEnergy" + fHistogramNameSuffix;
807 fClusterPositionAcceptedEnergy = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
808 fClusterPositionAcceptedEnergy->SetXTitle("#phi");
809 fClusterPositionAcceptedEnergy->SetYTitle("#eta");
810
811 histname = "fClusterPositionAllEnergy" + fHistogramNameSuffix;
812 fClusterPositionAllEnergy = new TH2D(histname.Data(), "Position of accepted neutral clusters",300, -TMath::Pi(),TMath::Pi(), 100, -0.7 , 0.7);
813 fClusterPositionAllEnergy->SetXTitle("#phi");
814 fClusterPositionAllEnergy->SetYTitle("#eta");
815
816 histname = "fClusterEnergy" + fHistogramNameSuffix;
817 fClusterEnergy = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
818 fClusterEnergy->SetYTitle("Number of clusters");
819 fClusterEnergy->SetXTitle("Energy of cluster");
820
821 histname = "fClusterEnergyCent" + fHistogramNameSuffix;
822 fClusterEnergyCent = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
823 fClusterEnergyCent->SetXTitle("Energy of cluster");
824 fClusterEnergyCent->SetYTitle("Centrality Bin");
825 fClusterEnergyCent->SetZTitle("Number of clusters");
826
827 histname = "fClusterEnergyModifiedTrackMatchesCent" + fHistogramNameSuffix;
828 fClusterEnergyModifiedTrackMatchesCent = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
829 fClusterEnergyModifiedTrackMatchesCent->SetXTitle("Energy of cluster");
830 fClusterEnergyModifiedTrackMatchesCent->SetYTitle("Centrality Bin");
831 fClusterEnergyModifiedTrackMatchesCent->SetZTitle("Number of clusters");
832
833 histname = "fClusterEnergyCentMatched" + fHistogramNameSuffix;
834 fClusterEnergyCentMatched = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
835 fClusterEnergyCentMatched->SetXTitle("Energy of cluster");
836 fClusterEnergyCentMatched->SetYTitle("Centrality Bin");
837 fClusterEnergyCentMatched->SetZTitle("Number of Clusters");
838
839 histname = "fClusterEnergyCentNotMatched" + fHistogramNameSuffix;
840 fClusterEnergyCentNotMatched = new TH2F(histname.Data(), histname.Data(), 100, 0, 5,20,-0.5,19.5);
841 fClusterEnergyCentNotMatched->SetXTitle("Energy of cluster");
842 fClusterEnergyCentNotMatched->SetYTitle("Centrality Bin");
843 fClusterEnergyCentNotMatched->SetZTitle("Number of clusters");
844
845 histname = "fClusterEt" + fHistogramNameSuffix;
846 fClusterEt = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
847 fClusterEt->SetXTitle("Number of clusters");
848 fClusterEt->SetYTitle("E_{T} of cluster");
849
850 histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix;
851 fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
852
853 histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix;
854 fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
855
856 histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix;
857 fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5);
858
859 fHistMatchedTracksEvspTvsCent = new TH3F("fHistMatchedTracksEvspTvsCent", "fHistMatchedTracksEvspTvsCent",100, 0, 3,100,0,3,20,-0.5,19.5);
860 fHistMatchedTracksEvspTvsCentEffCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffCorr", "fHistMatchedTracksEvspTvsCentEffCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
861 fHistMatchedTracksEvspTvsCentEffTMCorr = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr", "fHistMatchedTracksEvspTvsCentEffTMCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
862 fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr = new TH3F("fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr", "fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
863 fHistMatchedTracksEvspTvsCentEffTMCorr500MeV = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr500MeV", "fHistMatchedTracksEvspTvsCentEffTMCorr500MeV",100, 0, 3,100,0,3,20,-0.5,19.5);
864
865 float max = 200*scale;
866 if(fHistogramNameSuffix.Contains("P")){max = 100;}
867 fHistFoundHadronsvsCent = new TH2F("fHistFoundHadronsvsCent","fHistFoundHadronsvsCent",100,0,max,20,-0.5,19.5);
868 fHistNotFoundHadronsvsCent = new TH2F("fHistNotFoundHadronsvsCent","fHistNotFoundHadronsvsCent",100,0,max,20,-0.5,19.5);
869 fHistFoundHadronsEtvsCent = new TH2F("fHistFoundHadronsEtvsCent","fHistFoundHadronsEtvsCent",100,0,max,20,-0.5,19.5);
870 fHistNotFoundHadronsEtvsCent = new TH2F("fHistNotFoundHadronsEtvsCent","fHistNotFoundHadronsEtvsCent",100,0,max,20,-0.5,19.5);
871 fHistFoundHadronsvsCent500MeV = new TH2F("fHistFoundHadronsvsCent500MeV","fHistFoundHadronsvsCent500MeV",100,0,max,20,-0.5,19.5);
872 fHistNotFoundHadronsvsCent500MeV = new TH2F("fHistNotFoundHadronsvsCent500MeV","fHistNotFoundHadronsvsCent500MeV",100,0,max,20,-0.5,19.5);
873 fHistFoundHadronsEtvsCent500MeV = new TH2F("fHistFoundHadronsEtvsCent500MeV","fHistFoundHadronsEtvsCent500MeV",100,0,max,20,-0.5,19.5);
874 fHistNotFoundHadronsEtvsCent500MeV = new TH2F("fHistNotFoundHadronsEtvsCent500MeV","fHistNotFoundHadronsEtvsCent500MeV",100,0,max,20,-0.5,19.5);
875
876 fHistTotRawEtEffCorr = new TH2F("fHistTotRawEtEffCorr","fHistTotRawEtEffCorr",250,0,250*scale,20,-0.5,19.5);
877 fHistTotRawEt = new TH2F("fHistTotRawEt","fHistTotRawEt",250,0,250*scale,20,-0.5,19.5);
878 fHistTotRawEtEffCorr500MeV = new TH2F("fHistTotRawEtEffCorr500MeV","fHistTotRawEtEffCorr500MeV",250,0,250*scale,20,-0.5,19.5);
879 fHistTotAllRawEt = new TH2F("fHistTotAllRawEt","fHistTotAllRawEt",250,0,250*scale,20,-0.5,19.5);
880 fHistTotAllRawEtEffCorr = new TH2F("fHistTotAllRawEtEffCorr","fHistTotAllRawEtEffCorr",250,0,250*scale,20,-0.5,19.5);
881 fHistNClustersPhosVsEmcal = new TH3F("fHistNClustersPhosVsEmcal","fHistNClustersPhosVsEmcal",50,0,50,250*scale,0,250*scale,20,-0.5,19);
882 fHistClusterSizeVsCent = new TH2F("fHistClusterSizeVsCent","fHistClusterSizeVsCent",10,0.5,10.5,20,-0.5,19.5);
883 fHistMatchedClusterSizeVsCent = new TH2F("fHistMatchedClusterSizeVsCent","fHistMatchedClusterSizeVsCent",10,0.5,10.5,20,-0.5,19.5);
884 fHistTotAllRawEtVsTotalPt = new TH2F("fHistTotAllRawEtVsTotalPt","fHistTotAllRawEtVsTotalPt",125,0,250*scale,200,0,2000);
885 fHistTotAllRawEtVsTotalPtVsCent = new TH3F("fHistTotAllRawEtVsTotalPtVsCent","fHistTotAllRawEtVsTotalPtVsCent",125,0,250*scale,200,0,2000,20,-0.5,19.5);
886 fHistTotMatchedRawEtVsTotalPtVsCent = new TH3F("fHistTotMatchedRawEtVsTotalPtVsCent","fHistTotMatchedRawEtVsTotalPtVsCent",250,0,250*scale,100,0,200,20,-0.5,19.5);
887
888 maxEt = 500*scale;
889 histname = "fHistNominalRawEt" + fHistogramNameSuffix;
890 fHistNominalRawEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
891 histname = "fHistNominalNonLinHighEt" + fHistogramNameSuffix;
892 fHistNominalNonLinHighEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
893 histname = "fHistNominalNonLinLowEt" + fHistogramNameSuffix;
894 fHistNominalNonLinLowEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
895 histname = "fHistNominalEffHighEt" + fHistogramNameSuffix;
896 fHistNominalEffHighEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
897 histname = "fHistNominalEffLowEt" + fHistogramNameSuffix;
898 fHistNominalEffLowEt = new TH2D(histname.Data(), histname.Data(),nbinsEt,minEt,maxEt,20,-0.5,19.5);
899
900 Float_t maxEtRange = 25*scale;
901 Float_t maxEtRangeHigh = 125*scale;
902 Float_t minEtRange = 0;
903 Int_t nbinsMult = 100;
904 Float_t maxMult = 3000;
905 Float_t minMult = 0;
906 Int_t nbinsCl = 250*scale;
907 Float_t maxCl = 500*scale;
908 Float_t minCl = 0;
909 fHistPIDProtonsTrackMatchedDepositedVsNch = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNch","PID'd protons deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
910 fHistPIDAntiProtonsTrackMatchedDepositedVsNch = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNch","PID'd #bar{p} E_{T} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
911 fHistPIDProtonsTrackMatchedDepositedVsNcl = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNcl","PID'd protons deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
912 fHistPIDAntiProtonsTrackMatchedDepositedVsNcl = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNcl","PID'd #bar{p} E_{T} deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
913 fHistPiKPTrackMatchedDepositedVsNch = new TH2F("fHistPiKPTrackMatchedDepositedVsNch","PiKP track matched",nbinsEt,minEtRange,maxEtRangeHigh,nbinsMult,minMult,maxMult);
914
915 fHistPIDProtonsTrackMatchedDepositedVsNchNoEff = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNchNoEff","PID'd protons deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
916 fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff","PID'd #bar{p} E_{T} deposited in calorimeter vs multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsMult,minMult,maxMult);
917 fHistPIDProtonsTrackMatchedDepositedVsNclNoEff = new TH2F("fHistPIDProtonsTrackMatchedDepositedVsNclNoEff","PID'd protons deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
918 fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff = new TH2F("fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff","PID'd #bar{p} E_{T} deposited in calorimeter vs cluster multiplicity",nbinsEt,minEtRange,maxEtRange,nbinsCl,minCl,maxCl);
919 fHistPiKPTrackMatchedDepositedVsNchNoEff = new TH2F("fHistPiKPTrackMatchedDepositedVsNchNoEff","PiKP track matched",nbinsEt,minEtRange,maxEtRangeHigh,nbinsMult,minMult,maxMult);
920
921
922 fHistCentVsNchVsNclReco = new TH3F("fHistCentVsNchVsNclReco","Cent bin vs Nch Vs NCl",20,-0.5,19.5,nbinsMult,minMult,maxMult,nbinsCl,minCl,maxCl);
923
924 fHistRawSignalReco = new TH1F("fHistRawSignalReco","fHistRawSignalReco",20,-0.5,19.5);
925 fHistEffCorrSignalReco = new TH1F("fHistEffCorrSignalReco","fHistEffCorrSignalReco",20,-0.5,19.5);
926 fHistRecoRCorrVsPtVsCent = new TH3F("fHistRecoRCorrVsPtVsCent","fHistRecoRCorrVsPtVsCent",72,0,2,50,0,10,20,-0.5,19.5);
927
928}
929Double_t AliAnalysisEtReconstructed::ApplyModifiedCorrections(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t cent)
930{
931 Float_t pos[3];
932 cluster.GetPosition(pos);
933 TVector3 cp(pos);
934 Double_t corrEnergy = fReCorrections->CorrectedEnergy(cluster.E(),cent);
935
936 Double_t factorNonLin = GetCorrectionModification(cluster, nonLinCorr,effCorr,cent);
937
938 cout<<"Warning: This function should not get called!"<<endl;
939 //std::cout << "Original energy: " << cluster.E() << ", corrected energy: " << corrEnergy << std::endl;
940 return TMath::Sin(cp.Theta())*corrEnergy*factorNonLin;
941}
942
943Double_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
944 if(nonLinCorr==0){
945 cout<<"Warning: This function should not get called!"<<endl;//this statement is basically here to avoid a compilation warning
946 }
947 if(effCorr==0){
948 cout<<"Warning: This function should not get called!"<<endl;//this statement is basically here to avoid a compilation warning
949 }
950 return cluster.E()*cent;
951}