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