]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGLF/totEt/AliAnalysisEtReconstructed.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fixing Oysteins old cross check histograms
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisEtReconstructed.cxx
CommitLineData
cf6522d1 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//_________________________________________________________________________
2fbf38ac 9
10#include "AliAnalysisEtReconstructed.h"
11#include "AliAnalysisEtCuts.h"
12#include "AliESDtrack.h"
ba136eb4 13#include "AliEMCALTrack.h"
2fbf38ac 14#include "AliESDCaloCluster.h"
15#include "TVector3.h"
ba136eb4 16#include "TGeoGlobalMagField.h"
17#include "AliMagF.h"
2fbf38ac 18#include "AliVEvent.h"
19#include "AliESDEvent.h"
b5821c13 20#include "AliESDtrackCuts.h"
2fbf38ac 21#include "AliVParticle.h"
87efb15c 22#include "TDatabasePDG.h"
23#include "TList.h"
b5821c13 24#include "AliESDpid.h"
2fbf38ac 25#include <iostream>
09fcb185 26#include "TH3F.h"
2fbf38ac 27#include "TH2F.h"
ef647350 28#include "TH2I.h"
29#include "TH1I.h"
30#include "TFile.h"
964c8159 31#include "AliAnalysisHadEtCorrections.h"
ef647350 32#include "AliAnalysisEtSelector.h"
0f6416f3 33#include "AliLog.h"
e9da35da 34#include "AliCentrality.h"
ef647350 35#include "AliPHOSGeoUtils.h"
36#include "AliPHOSGeometry.h"
d3ce32b8 37#include "AliAnalysisEtRecEffCorrection.h"
9a365626 38#include "AliESDpid.h"
0f6416f3 39
2fbf38ac 40
16abb579 41using namespace std;
42
43ClassImp(AliAnalysisEtReconstructed);
44
45
2fbf38ac 46AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() :
f61cec2f 47 AliAnalysisEt()
48 ,fCorrections(0)
49 ,fPidCut(0)
fb116385 50 ,nChargedHadronsMeasured(0)
51 ,nChargedHadronsTotal(0)
f61cec2f 52 ,fHistChargedPionEnergyDeposit(0)
53 ,fHistProtonEnergyDeposit(0)
54 ,fHistAntiProtonEnergyDeposit(0)
55 ,fHistChargedKaonEnergyDeposit(0)
56 ,fHistMuonEnergyDeposit(0)
57 ,fHistRemovedEnergy(0)
58 ,fGeomCorrection(1.0)
b2c10007 59 ,fEMinCorrection(1.0/0.687)
f61cec2f 60 ,fRecEffCorrection(1.0)
43dd5a38 61 ,fClusterPositionAccepted(0)
62 ,fClusterPositionAll(0)
63 ,fClusterPositionAcceptedEnergy(0)
64 ,fClusterPositionAllEnergy(0)
53302bfe 65 ,fClusterEnergy(0)
9a365626 66 ,fClusterEnergyCent(0)
d386205e 67 ,fClusterEnergyModifiedTrackMatchesCent(0)
9a365626 68 ,fClusterEnergyCentMatched(0)
69 ,fClusterEnergyCentNotMatched(0)
53302bfe 70 ,fClusterEt(0)
f61cec2f 71 ,fHistChargedEnergyRemoved(0)
72 ,fHistNeutralEnergyRemoved(0)
73 ,fHistGammaEnergyAdded(0)
ac610b08 74 ,fHistMatchedTracksEvspTvsCent(0)
75 ,fHistMatchedTracksEvspTvsCentEffCorr(0)
459e9c44 76 ,fHistMatchedTracksEvspTvsCentEffTMCorr(0)
3af79c0f 77 ,fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr(0)
5881f036 78 ,fHistMatchedTracksEvspTvsCentEffTMCorr500MeV(0)
6a152780 79 ,fHistFoundHadronsvsCent(0)
80 ,fHistNotFoundHadronsvsCent(0)
81 ,fHistFoundHadronsEtvsCent(0)
82 ,fHistNotFoundHadronsEtvsCent(0)
43dd5a38 83 ,fHistFoundHadronsvsCent500MeV(0)
84 ,fHistNotFoundHadronsvsCent500MeV(0)
85 ,fHistFoundHadronsEtvsCent500MeV(0)
86 ,fHistNotFoundHadronsEtvsCent500MeV(0)
d3ce32b8 87 ,fHistNominalRawEt(0)
88 ,fHistNominalNonLinHighEt(0)
89 ,fHistNominalNonLinLowEt(0)
90 ,fHistNominalEffHighEt(0)
91 ,fHistNominalEffLowEt(0)
43dd5a38 92 ,fHistTotRawEtEffCorr(0)
3e9c52ca 93 ,fHistTotRawEt(0)
43dd5a38 94 ,fHistTotRawEtEffCorr500MeV(0)
95 ,fHistTotAllRawEt(0)
96 ,fHistTotAllRawEtEffCorr(0)
9a365626 97 ,fHistNClustersPhosVsEmcal(0)
98 ,fHistClusterSizeVsCent(0)
99 ,fHistMatchedClusterSizeVsCent(0)
100 ,fHistTotAllRawEtVsTotalPt(0)
101 ,fHistTotAllRawEtVsTotalPtVsCent(0)
102 ,fHistTotMatchedRawEtVsTotalPtVsCent(0)
103 ,fHistPIDProtonsTrackMatchedDepositedVsNch(0)
104 ,fHistPIDAntiProtonsTrackMatchedDepositedVsNch(0)
0fcbed20 105 ,fHistPIDProtonsTrackMatchedDepositedVsNcl(0)
106 ,fHistPIDAntiProtonsTrackMatchedDepositedVsNcl(0)
9a365626 107 ,fHistPiKPTrackMatchedDepositedVsNch(0)
fb116385 108 //,
109 ,fHistPIDProtonsTrackMatchedDepositedVsNchNoEff(0)
110 ,fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff(0)
111 ,fHistPIDProtonsTrackMatchedDepositedVsNclNoEff(0)
112 ,fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff(0)
113 ,fHistPiKPTrackMatchedDepositedVsNchNoEff(0)
9a365626 114 ,fHistCentVsNchVsNclReco(0)
fb116385 115 ,fHistRawSignalReco(0)
116 ,fHistEffCorrSignalReco(0)
02d47689 117 ,fHistRecoRCorrVsPtVsCent(0)
2fbf38ac 118{
119
120}
121
e9da35da 122AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed()
f61cec2f 123{//destructor
e9da35da 124 delete fCorrections;
ef647350 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 */
d0c22dcc 129 delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */
130
131 delete fHistRemovedEnergy; // removed energy
43dd5a38 132 delete fClusterPositionAccepted;
133 delete fClusterPositionAll;
134 delete fClusterPositionAcceptedEnergy;
135 delete fClusterPositionAllEnergy;
53302bfe 136 delete fClusterEnergy;
9a365626 137 delete fClusterEnergyCent;
d386205e 138 delete fClusterEnergyModifiedTrackMatchesCent;
9a365626 139 delete fClusterEnergyCentMatched;
140 delete fClusterEnergyCentNotMatched;
53302bfe 141 delete fClusterEt;
311c6540 142 delete fHistChargedEnergyRemoved;
143 delete fHistNeutralEnergyRemoved;
144 delete fHistGammaEnergyAdded;
ac610b08 145 delete fHistMatchedTracksEvspTvsCent;
146 delete fHistMatchedTracksEvspTvsCentEffCorr;
459e9c44 147 delete fHistMatchedTracksEvspTvsCentEffTMCorr;
3af79c0f 148 delete fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr;
5881f036 149 delete fHistMatchedTracksEvspTvsCentEffTMCorr500MeV;
6a152780 150 delete fHistFoundHadronsvsCent;
151 delete fHistNotFoundHadronsvsCent;
152 delete fHistFoundHadronsEtvsCent;
153 delete fHistNotFoundHadronsEtvsCent;
43dd5a38 154 delete fHistFoundHadronsvsCent500MeV;
155 delete fHistNotFoundHadronsvsCent500MeV;
156 delete fHistFoundHadronsEtvsCent500MeV;
157 delete fHistNotFoundHadronsEtvsCent500MeV;
d3ce32b8 158 delete fHistNominalRawEt;
159 delete fHistNominalNonLinHighEt;
160 delete fHistNominalNonLinLowEt;
161 delete fHistNominalEffHighEt;
162 delete fHistNominalEffLowEt;
43dd5a38 163 delete fHistTotRawEtEffCorr;
3e9c52ca 164 delete fHistTotRawEt;
43dd5a38 165 delete fHistTotAllRawEt;
166 delete fHistTotAllRawEtEffCorr;
167 delete fHistTotRawEtEffCorr500MeV;
9a365626 168 delete fHistNClustersPhosVsEmcal;
169 delete fHistClusterSizeVsCent;
170 delete fHistMatchedClusterSizeVsCent;
171 delete fHistTotAllRawEtVsTotalPt;
172 delete fHistTotAllRawEtVsTotalPtVsCent;
173 delete fHistTotMatchedRawEtVsTotalPtVsCent;
174 delete fHistPIDProtonsTrackMatchedDepositedVsNch;
175 delete fHistPIDAntiProtonsTrackMatchedDepositedVsNch;
0fcbed20 176 delete fHistPIDProtonsTrackMatchedDepositedVsNcl;
177 delete fHistPIDAntiProtonsTrackMatchedDepositedVsNcl;
9a365626 178 delete fHistPiKPTrackMatchedDepositedVsNch;
fb116385 179 delete fHistPIDProtonsTrackMatchedDepositedVsNchNoEff;
180 delete fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff;
181 delete fHistPIDProtonsTrackMatchedDepositedVsNclNoEff;
182 delete fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff;
183 delete fHistPiKPTrackMatchedDepositedVsNchNoEff;
9a365626 184 delete fHistCentVsNchVsNclReco;
fb116385 185 delete fHistRawSignalReco;
186 delete fHistEffCorrSignalReco;
02d47689 187 delete fHistRecoRCorrVsPtVsCent;
cf6522d1 188}
189
190Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev)
e9da35da 191{
f61cec2f 192
193 //AliAnalysisEt::AnalyseEvent(ev);
e9da35da 194 // analyse ESD event
2fbf38ac 195 ResetEventValues();
e9da35da 196 if (!ev) {
197 AliFatal("ERROR: Event does not exist");
198 return 0;
199 }
200
2fbf38ac 201 AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev);
e9da35da 202 if (!event) {
203 AliFatal("ERROR: ESD Event does not exist");
204 return 0;
ec956c46 205 }
2aab9269 206 if(!fSelector){
207 AliFatal("ERROR: fSelector does not exist");
208 return 0;
209 }
ef647350 210 fSelector->SetEvent(event);
f61cec2f 211
743ce29f 212 Int_t cent = -1;
d3ce32b8 213 fCentrality = event->GetCentrality();
f6b36c54 214 if (fCentrality && cent)
743ce29f 215 {
d3ce32b8 216 cent = fCentrality->GetCentralityClass5("V0M");
217 fCentClass = fCentrality->GetCentralityClass5("V0M");
743ce29f 218 }
219
9a365626 220
221 //for PID
3af79c0f 222 //AliESDpid *pID = new AliESDpid();
223 //pID->MakePID(event);
9a365626 224 Float_t etPIDProtons = 0.0;
225 Float_t etPIDAntiProtons = 0.0;
226 Float_t etPiKPMatched = 0.0;
fb116385 227 Float_t etPIDProtonsNoEff = 0.0;
228 Float_t etPIDAntiProtonsNoEff = 0.0;
229 Float_t etPiKPMatchedNoEff = 0.0;
9a365626 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();
3af79c0f 246 //pID->MakeITSPID(track);
9a365626 247
248
249 }
250 }
251
5881f036 252 //TRefArray *caloClusters = fSelector->GetClusters();//just gets the correct set of clusters - does not apply any cuts
253 //Float_t fClusterMult = caloClusters->GetEntries();
09fcb185 254
d3ce32b8 255 Float_t nominalRawEt = 0;
5881f036 256 Float_t totEt500MeV = 0;
d3ce32b8 257 Float_t nonlinHighRawEt = 0;
258 Float_t nonlinLowRawEt = 0;
259 Float_t effHighRawEt = 0;
260 Float_t effLowRawEt = 0;
3e9c52ca 261 Float_t uncorrEt = 0;
f2a956c2 262 Float_t rawSignal = 0;
263 Float_t effCorrSignal = 0;
d3ce32b8 264
fb116385 265 nChargedHadronsMeasured = 0.0;
266 nChargedHadronsTotal = 0.0;
6a152780 267 Float_t nChargedHadronsEtMeasured = 0.0;
268 Float_t nChargedHadronsEtTotal = 0.0;
43dd5a38 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;
f809f93b 275 Float_t fTotRawEtEffCorr = 0.0;
43dd5a38 276 Float_t fTotAllRawEtEffCorr = 0.0;
9a365626 277 Int_t nPhosClusters = 0;
278 Int_t nEmcalClusters = 0;
d3ce32b8 279
fb116385 280
281 TRefArray *caloClusters = fSelector->GetClusters();
282 Int_t nCluster = caloClusters->GetEntries();
283
284 for (int iCluster = 0; iCluster < nCluster; iCluster++ )
2fbf38ac 285 {
fb116385 286 AliESDCaloCluster* cluster = ( AliESDCaloCluster* )caloClusters->At( iCluster );
2fbf38ac 287 if (!cluster)
288 {
e9da35da 289 AliError(Form("ERROR: Could not get cluster %d", iCluster));
2fbf38ac 290 continue;
291 }
ef647350 292 int x = 0;
f61cec2f 293 fCutFlow->Fill(x++);
9a365626 294 if(cluster->IsEMCAL()) nEmcalClusters++;
295 else nPhosClusters++;
86e7d5db 296 if(!fSelector->IsDetectorCluster(*cluster)) continue;
f61cec2f 297 fCutFlow->Fill(x++);
86e7d5db 298 if(!fSelector->PassMinEnergyCut(*cluster)) continue;
f61cec2f 299 fCutFlow->Fill(x++);
86e7d5db 300 if (!fSelector->PassDistanceToBadChannelCut(*cluster)) continue;
f61cec2f 301 fCutFlow->Fill(x++);
5881f036 302 if (!fSelector->CutGeometricalAcceptance(*cluster)) continue;
303 //fCutFlow->Fill(x++);
2fbf38ac 304 Float_t pos[3];
e9da35da 305
2fbf38ac 306 cluster->GetPosition(pos);
e9da35da 307 TVector3 cp(pos);
43dd5a38 308 fClusterPositionAll->Fill(cp.Phi(), cp.PseudoRapidity());
02d47689 309 Float_t fReconstructedE = cluster->E();
6280e578 310 Float_t lostEnergy = 0.0;
311 Float_t lostTrackPt = 0.0;
02d47689 312 fClusterPositionAllEnergy->Fill(cp.Phi(), cp.PseudoRapidity(),GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE);
87efb15c 313
5881f036 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
02d47689 315 fTotAllRawEt += TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
27d85daa 316 fTotAllRawEtEffCorr +=GetCorrectionModification(*cluster,0,0,cent)* CorrectForReconstructionEfficiency(*cluster,cent);
5881f036 317
02d47689 318 fClusterEnergyCent->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
09fcb185 319 Bool_t matched = kTRUE;//default to no track matched
6280e578 320 Bool_t countasmatched = kFALSE;
96e7f328 321 Bool_t correctedcluster = kFALSE;
6280e578 322
09fcb185 323 Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex();//find the index of the matched track
584f2478 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
02d47689 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);
6280e578 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
7ee3992c 335 // if(fReconstructedE - fsub* track->P() > 0.0){
96e7f328 336 //cout<<"match corrected"<<endl;
02d47689 337 fReconstructedE = fReconstructedE - fsub* track->P();
338 matched = kFALSE;
6280e578 339 countasmatched = kTRUE;
96e7f328 340 correctedcluster = kTRUE;
6280e578 341 lostEnergy = fsub* track->P();
342 lostTrackPt = track->Pt();
d386205e 343 fClusterEnergyModifiedTrackMatchesCent->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
02d47689 344 }
96e7f328 345// else{
346// cerr<<"match passed ";
347// cerr<<"E "<<fReconstructedE<<" fsubmeanhad "<<fsub<<" p "<< track->P();
348// if(correctedcluster) cout<<" corrected";
349// cerr<<endl;
350// }
02d47689 351 }
584f2478 352 }
09fcb185 353 }
e9da35da 354
e9da35da 355
356 if (matched)
357 {
f61cec2f 358
e9da35da 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 {
9a365626 366 totalMatchedPt +=track->Pt();
02d47689 367 fClusterEnergyCentMatched->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE,cent);
9a365626 368 fHistMatchedClusterSizeVsCent->Fill(cluster->GetNCells(),cent);
369
5881f036 370 float eff = fTmCorrections->TrackMatchingEfficiency(track->Pt(),cent);
584f2478 371 if(TMath::Abs(eff)<1e-5) eff = 1.0;
9a365626 372 //cout<<"pt "<<track->Pt()<<" eff "<<eff<<" total "<<nChargedHadronsTotal<<endl;
6a152780 373 nChargedHadronsMeasured++;
584f2478 374 nChargedHadronsTotal += 1/eff;
02d47689 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;
3e9c52ca 377 //One efficiency is the gamma efficiency and the other is the track matching efficiency.
02d47689 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;
9a365626 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;
02d47689 386 etPiKPMatchedNoEff +=TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
9a365626 387 if(isProton){
388 if(track->Charge()>0){
389 etPIDProtons += effCorrEt;
02d47689 390 etPIDProtonsNoEff +=TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
9a365626 391 }
392 else{
02d47689 393 etPIDAntiProtonsNoEff +=TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
9a365626 394 etPIDAntiProtons += effCorrEt;
395 }
396 }
02d47689 397 if(TMath::Sin(cp.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE>0.5){
43dd5a38 398 nChargedHadronsMeasured500MeV++;
399 nChargedHadronsTotal500MeV += 1/eff;
02d47689 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;
43dd5a38 402 }
5881f036 403 cluster->GetPosition(pos);
404 TVector3 p2(pos);
02d47689 405 uncorrEt += TMath::Sin(p2.Theta())*GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE;
96e7f328 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
6280e578 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 }
6280e578 421 const Double_t *pidWeights = track->PID();
422
423 Double_t maxpidweight = 0;
424 Int_t maxpid = 0;
425
426 if (pidWeights)
e9da35da 427 {
6280e578 428 for (Int_t p =0; p < AliPID::kSPECIES; p++)
e9da35da 429 {
6280e578 430 if (pidWeights[p] > maxpidweight)
e9da35da 431 {
6280e578 432 maxpidweight = pidWeights[p];
433 maxpid = p;
e9da35da 434 }
435 }
6280e578 436 if (fCuts->GetHistMakeTreeDeposit() && fDepositTree)
e9da35da 437 {
6280e578 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 }
e9da35da 451 }
6280e578 452
453 if (maxpidweight > fPidCut)
e9da35da 454 {
6280e578 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
02d47689 459 //Float_t et = fReconstructedE * TMath::Sin(theta);
e9da35da 460 if (maxpid == AliPID::kProton)
461 {
462
463 if (track->Charge() == 1)
464 {
02d47689 465 fHistProtonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
e9da35da 466 }
467 else if (track->Charge() == -1)
468 {
02d47689 469 fHistAntiProtonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
e9da35da 470 }
471 }
472 else if (maxpid == AliPID::kPion)
473 {
02d47689 474 fHistChargedPionEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
e9da35da 475 }
476 else if (maxpid == AliPID::kKaon)
477 {
02d47689 478 fHistChargedKaonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
e9da35da 479 }
480 else if (maxpid == AliPID::kMuon)
481 {
02d47689 482 fHistMuonEnergyDeposit->Fill(GetCorrectionModification(*cluster,0,0,cent)*fReconstructedE, track->E());
e9da35da 483 }
484 }
485 }
486 }
487 }
488 //continue;
ba136eb4 489 } // distance
6a152780 490 else{//these are clusters which were not track matched
f61cec2f 491 fCutFlow->Fill(x++);
492 //std::cout << x++ << std::endl;
6a152780 493
02d47689 494 //if (fReconstructedE > fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue;
495 //if (fReconstructedE < fClusterEnergyCut) continue;
6a152780 496 cluster->GetPosition(pos);
497
6280e578 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++;
e9da35da 539 }
ef647350 540 fMultiplicity++;
541 }
f61cec2f 542
fb116385 543
544 fHistRawSignalReco->Fill(rawSignal);
545 fHistEffCorrSignalReco->Fill(effCorrSignal);
546
9a365626 547 fHistNClustersPhosVsEmcal->Fill(nPhosClusters,nEmcalClusters,cent);
ef647350 548 fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity);
549 fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity);
550 fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity);
551 fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity);
f61cec2f 552
ef647350 553 fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity);
554 fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity);
2fbf38ac 555
8529c4c1 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);
e9da35da 562 fHistRemovedEnergy->Fill(removedEnergy);
f61cec2f 563
459e9c44 564 fTotNeutralEtAcc = fTotNeutralEt;
f809f93b 565 //fHistTotRawEtEffCorr->Fill(fTotNeutralEt,cent);
566 fHistTotRawEtEffCorr->Fill(fTotRawEtEffCorr,cent);
43dd5a38 567 fHistTotRawEt->Fill(fTotRawEt,cent);
568 fHistTotAllRawEt->Fill(fTotAllRawEt,cent);
9a365626 569 fHistTotAllRawEtVsTotalPt->Fill(fTotAllRawEt,totalPt);
570 fHistTotAllRawEtVsTotalPtVsCent->Fill(fTotAllRawEt,totalPt,cent);
571 fHistTotMatchedRawEtVsTotalPtVsCent->Fill(fTotAllRawEt,totalMatchedPt,cent);
43dd5a38 572 fHistTotAllRawEtEffCorr->Fill(fTotAllRawEtEffCorr,cent);
f809f93b 573 //cout<<"fTotAllRawEtEffCorr "<<fTotAllRawEtEffCorr<<" fTotAllRawEt "<<fTotAllRawEt<<" fTotRawEtEffCorr "<<fTotRawEtEffCorr<<"("<<fTotNeutralEt<<")"<<" fTotRawEt "<<fTotRawEt<<endl;
3e9c52ca 574 //cout<<"uncorr "<<uncorrEt<<" raw "<<nominalRawEt<<" tot raw "<<fTotNeutralEt;
8529c4c1 575 //cout<<" raw "<<fTotNeutralEt<<" removed "<<removedEnergy<<" etmin "<<GetMinEtCorrection(cent)<<" final ";
576 if(GetMinEtCorrection(cent)>0) fTotNeutralEt = (fTotNeutralEt - removedEnergy)/GetMinEtCorrection(cent);
577 //cout<<fTotNeutralEt<<endl;
3e9c52ca 578 //cout<<" tot corr "<<fTotNeutralEt<<endl;
2fbf38ac 579 fTotEt = fTotChargedEt + fTotNeutralEt;
f61cec2f 580// Fill the histograms...0
2fbf38ac 581 FillHistograms();
b2c10007 582 //std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl;
d3ce32b8 583 //cout<<"cent "<<cent<<" cluster mult "<<fClusterMult<<" fTotNeutralEt "<<fTotNeutralEt<<" nominalRawEt "<<nominalRawEt<<endl;
584 fHistNominalRawEt->Fill(nominalRawEt,cent);
43dd5a38 585 fHistTotRawEtEffCorr500MeV->Fill(totEt500MeV,cent);
d3ce32b8 586 fHistNominalNonLinHighEt->Fill(nonlinHighRawEt,cent);
587 fHistNominalNonLinLowEt->Fill(nonlinLowRawEt,cent);
588 fHistNominalEffHighEt->Fill(effHighRawEt,cent);
589 fHistNominalEffLowEt->Fill(effLowRawEt,cent);
6a152780 590 fHistFoundHadronsvsCent->Fill(nChargedHadronsMeasured,cent);
591 fHistNotFoundHadronsvsCent->Fill(nChargedHadronsTotal-nChargedHadronsMeasured,cent);
592 fHistFoundHadronsEtvsCent->Fill(nChargedHadronsEtMeasured,cent);
593 fHistNotFoundHadronsEtvsCent->Fill(nChargedHadronsEtTotal-nChargedHadronsEtMeasured,cent);
9a365626 594 //cout<<"found "<<nChargedHadronsMeasured<<" total "<<nChargedHadronsTotal<<" not found "<<nChargedHadronsTotal-nChargedHadronsMeasured<<" found "<< nChargedHadronsMeasured500MeV<<" not found "<< nChargedHadronsTotal500MeV-nChargedHadronsMeasured500MeV <<" total "<<nChargedHadronsTotal500MeV<<endl;
43dd5a38 595 fHistFoundHadronsvsCent500MeV->Fill(nChargedHadronsMeasured500MeV,cent);
596 fHistNotFoundHadronsvsCent500MeV->Fill(nChargedHadronsTotal500MeV-nChargedHadronsMeasured500MeV,cent);
597 fHistFoundHadronsEtvsCent500MeV->Fill(nChargedHadronsEtMeasured500MeV,cent);
598 fHistNotFoundHadronsEtvsCent500MeV->Fill(nChargedHadronsEtTotal500MeV-nChargedHadronsEtMeasured500MeV,cent);
6a152780 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;
9a365626 606 fHistPIDProtonsTrackMatchedDepositedVsNch->Fill(etPIDProtons,multiplicity);
607 fHistPIDAntiProtonsTrackMatchedDepositedVsNch->Fill(etPIDAntiProtons,multiplicity);
fb116385 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);
9a365626 615 fHistPiKPTrackMatchedDepositedVsNch->Fill(etPiKPMatched,multiplicity);
fb116385 616 fHistPiKPTrackMatchedDepositedVsNchNoEff->Fill(etPiKPMatchedNoEff,multiplicity);
3af79c0f 617 //delete pID;
2fbf38ac 618 return 0;
619}
620
621bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track)
f61cec2f 622{ // check vertex
2fbf38ac 623
624 Float_t bxy = 999.;
625 Float_t bz = 999.;
e9da35da 626 if (!track) {
627 AliError("ERROR: no track");
628 return kFALSE;
ec956c46 629 }
1b8c3d66 630 AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track);
e9da35da 631 if (!esdTrack) {
632 AliError("ERROR: no track");
633 return kFALSE;
1b8c3d66 634 }
635 esdTrack->GetImpactParametersTPC(bxy,bz);
ec956c46 636
2fbf38ac 637
e9da35da 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());
2fbf38ac 643
4998becf 644 return status;
2fbf38ac 645}
646
647void AliAnalysisEtReconstructed::Init()
f61cec2f 648{ // Init
2fbf38ac 649 AliAnalysisEt::Init();
83d0f02c 650 fPidCut = fCuts->GetReconstructedPidCut();
ba136eb4 651 TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG));
e9da35da 652 if (!fCorrections) {
653 cout<<"Warning! You have not set corrections. Your code will crash. You have to set the corrections."<<endl;
0f97be4c 654 }
2fbf38ac 655}
656
657bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField)
f61cec2f 658{ // propagate track to detector radius
e9da35da 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());
2fbf38ac 669
670 Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField);
671
cf6522d1 672 // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt());
f61cec2f 673 return prop && fSelector->CutGeometricalAcceptance(*esdTrack);
2fbf38ac 674}
675
87efb15c 676void AliAnalysisEtReconstructed::FillOutputList(TList* list)
f61cec2f 677{ // add some extra histograms to the ones from base class
87efb15c 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);
e9da35da 685
686 list->Add(fHistRemovedEnergy);
43dd5a38 687 list->Add(fClusterPositionAccepted);
688 list->Add(fClusterPositionAll);
689 list->Add(fClusterPositionAcceptedEnergy);
690 list->Add(fClusterPositionAllEnergy);
53302bfe 691 list->Add(fClusterEnergy);
9a365626 692 list->Add(fClusterEnergyCent);
d386205e 693 list->Add(fClusterEnergyModifiedTrackMatchesCent);
9a365626 694 list->Add(fClusterEnergyCentMatched);
695 list->Add(fClusterEnergyCentNotMatched);
53302bfe 696 list->Add(fClusterEt);
f61cec2f 697
ef647350 698 list->Add(fHistChargedEnergyRemoved);
699 list->Add(fHistNeutralEnergyRemoved);
700 list->Add(fHistGammaEnergyAdded);
ac610b08 701 list->Add(fHistMatchedTracksEvspTvsCent);
702 list->Add(fHistMatchedTracksEvspTvsCentEffCorr);
459e9c44 703 list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr);
3af79c0f 704 list->Add(fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr);
5881f036 705 list->Add(fHistMatchedTracksEvspTvsCentEffTMCorr500MeV);
6a152780 706 list->Add(fHistFoundHadronsvsCent);
707 list->Add(fHistNotFoundHadronsvsCent);
708 list->Add(fHistFoundHadronsEtvsCent);
709 list->Add(fHistNotFoundHadronsEtvsCent);
43dd5a38 710 list->Add(fHistFoundHadronsvsCent500MeV);
711 list->Add(fHistNotFoundHadronsvsCent500MeV);
712 list->Add(fHistFoundHadronsEtvsCent500MeV);
713 list->Add(fHistNotFoundHadronsEtvsCent500MeV);
d3ce32b8 714 list->Add(fHistNominalRawEt);
715 list->Add(fHistNominalNonLinHighEt);
716 list->Add(fHistNominalNonLinLowEt);
717 list->Add(fHistNominalEffHighEt);
718 list->Add(fHistNominalEffLowEt);
43dd5a38 719 list->Add(fHistTotRawEtEffCorr);
720 list->Add(fHistTotRawEtEffCorr500MeV);
721 list->Add(fHistTotAllRawEtEffCorr);
3e9c52ca 722 list->Add(fHistTotRawEt);
43dd5a38 723 list->Add(fHistTotAllRawEt);
9a365626 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);
0fcbed20 732 list->Add(fHistPIDProtonsTrackMatchedDepositedVsNcl);
733 list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNcl);
9a365626 734 list->Add(fHistPiKPTrackMatchedDepositedVsNch);
fb116385 735 list->Add(fHistPIDProtonsTrackMatchedDepositedVsNchNoEff);
736 list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNchNoEff);
737 list->Add(fHistPIDProtonsTrackMatchedDepositedVsNclNoEff);
738 list->Add(fHistPIDAntiProtonsTrackMatchedDepositedVsNclNoEff);
739 list->Add(fHistPiKPTrackMatchedDepositedVsNchNoEff);
9a365626 740 list->Add(fHistCentVsNchVsNclReco);
fb116385 741 list->Add(fHistRawSignalReco);
742 list->Add(fHistEffCorrSignalReco);
02d47689 743 list->Add(fHistRecoRCorrVsPtVsCent);
87efb15c 744}
745
746void AliAnalysisEtReconstructed::CreateHistograms()
f61cec2f 747{ // add some extra histograms to the ones from base class
87efb15c 748 AliAnalysisEt::CreateHistograms();
749
6fe6205f 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;
0fa8c632 755 Double_t minEt = 0;
6fe6205f 756 Double_t maxEt = 10*scale;
0fa8c632 757
758 // possibly change histogram limits
53302bfe 759// if (fCuts) {
760// nbinsEt = fCuts->GetHistNbinsParticleEt();
761// minEt = fCuts->GetHistMinParticleEt();
762// maxEt = fCuts->GetHistMaxParticleEt();
763// }
0fa8c632 764
87efb15c 765 TString histname;
766 histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix;
0fa8c632 767 fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c 768 fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
769 fHistChargedPionEnergyDeposit->SetYTitle("Energy of track");
e9da35da 770
87efb15c 771 histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632 772 fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c 773 fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
774 fHistProtonEnergyDeposit->SetYTitle("Energy of track");
e9da35da 775
87efb15c 776 histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632 777 fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c 778 fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
779 fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track");
e9da35da 780
87efb15c 781 histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632 782 fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c 783 fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
784 fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track");
e9da35da 785
87efb15c 786 histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix;
0fa8c632 787 fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt);
87efb15c 788 fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter");
789 fHistMuonEnergyDeposit->SetYTitle("Energy of track");
e9da35da 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
43dd5a38 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");
ef647350 815
53302bfe 816 histname = "fClusterEnergy" + fHistogramNameSuffix;
817 fClusterEnergy = new TH1F(histname.Data(), histname.Data(), 100, 0, 5);
9a365626 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
d386205e 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
9a365626 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");
53302bfe 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");
ef647350 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);
e9da35da 858
5881f036 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);
3af79c0f 862 fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr = new TH3F("fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr", "fHistPeripheralMatchedTracksEvspTvsCentEffTMCorr",100, 0, 3,100,0,3,20,-0.5,19.5);
5881f036 863 fHistMatchedTracksEvspTvsCentEffTMCorr500MeV = new TH3F("fHistMatchedTracksEvspTvsCentEffTMCorr500MeV", "fHistMatchedTracksEvspTvsCentEffTMCorr500MeV",100, 0, 3,100,0,3,20,-0.5,19.5);
9a365626 864
6fe6205f 865 float max = 200*scale;
9a365626 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);
5881f036 875
6fe6205f 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);
9a365626 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);
6fe6205f 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);
d3ce32b8 887
6fe6205f 888 maxEt = 500*scale;
d3ce32b8 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
6fe6205f 900 Float_t maxEtRange = 25*scale;
901 Float_t maxEtRangeHigh = 125*scale;
9a365626 902 Float_t minEtRange = 0;
903 Int_t nbinsMult = 100;
904 Float_t maxMult = 3000;
905 Float_t minMult = 0;
6fe6205f 906 Int_t nbinsCl = 250*scale;
907 Float_t maxCl = 500*scale;
9a365626 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);
0fcbed20 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);
9a365626 913 fHistPiKPTrackMatchedDepositedVsNch = new TH2F("fHistPiKPTrackMatchedDepositedVsNch","PiKP track matched",nbinsEt,minEtRange,maxEtRangeHigh,nbinsMult,minMult,maxMult);
fb116385 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
9a365626 922 fHistCentVsNchVsNclReco = new TH3F("fHistCentVsNchVsNclReco","Cent bin vs Nch Vs NCl",20,-0.5,19.5,nbinsMult,minMult,maxMult,nbinsCl,minCl,maxCl);
fb116385 923
924 fHistRawSignalReco = new TH1F("fHistRawSignalReco","fHistRawSignalReco",20,-0.5,19.5);
925 fHistEffCorrSignalReco = new TH1F("fHistEffCorrSignalReco","fHistEffCorrSignalReco",20,-0.5,19.5);
02d47689 926 fHistRecoRCorrVsPtVsCent = new TH3F("fHistRecoRCorrVsPtVsCent","fHistRecoRCorrVsPtVsCent",72,0,2,50,0,10,20,-0.5,19.5);
fb116385 927
d3ce32b8 928}
5881f036 929Double_t AliAnalysisEtReconstructed::ApplyModifiedCorrections(const AliESDCaloCluster& cluster,Int_t nonLinCorr, Int_t effCorr, Int_t cent)
d3ce32b8 930{
931 Float_t pos[3];
932 cluster.GetPosition(pos);
933 TVector3 cp(pos);
5881f036 934 Double_t corrEnergy = fReCorrections->CorrectedEnergy(cluster.E(),cent);
d3ce32b8 935
5881f036 936 Double_t factorNonLin = GetCorrectionModification(cluster, nonLinCorr,effCorr,cent);
d3ce32b8 937
43dd5a38 938 cout<<"Warning: This function should not get called!"<<endl;
d3ce32b8 939 //std::cout << "Original energy: " << cluster.E() << ", corrected energy: " << corrEnergy << std::endl;
940 return TMath::Sin(cp.Theta())*corrEnergy*factorNonLin;
941}
87efb15c 942
5881f036 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
d3ce32b8 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 }
5881f036 950 return cluster.E()*cent;
2aab9269 951}