]>
Commit | Line | Data |
---|---|---|
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> |
26 | #include "TH2F.h" | |
ef647350 | 27 | #include "TH2I.h" |
28 | #include "TH1I.h" | |
29 | #include "TFile.h" | |
964c8159 | 30 | #include "AliAnalysisHadEtCorrections.h" |
ef647350 | 31 | #include "AliAnalysisEtSelector.h" |
0f6416f3 | 32 | #include "AliLog.h" |
e9da35da | 33 | #include "AliCentrality.h" |
ef647350 | 34 | #include "AliPHOSGeoUtils.h" |
35 | #include "AliPHOSGeometry.h" | |
0f6416f3 | 36 | |
2fbf38ac | 37 | |
16abb579 | 38 | using namespace std; |
39 | ||
40 | ClassImp(AliAnalysisEtReconstructed); | |
41 | ||
42 | ||
2fbf38ac | 43 | AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() : |
f61cec2f | 44 | AliAnalysisEt() |
45 | ,fCorrections(0) | |
46 | ,fPidCut(0) | |
47 | ,fHistChargedPionEnergyDeposit(0) | |
48 | ,fHistProtonEnergyDeposit(0) | |
49 | ,fHistAntiProtonEnergyDeposit(0) | |
50 | ,fHistChargedKaonEnergyDeposit(0) | |
51 | ,fHistMuonEnergyDeposit(0) | |
52 | ,fHistRemovedEnergy(0) | |
53 | ,fGeomCorrection(1.0) | |
b2c10007 | 54 | ,fEMinCorrection(1.0/0.687) |
f61cec2f | 55 | ,fRecEffCorrection(1.0) |
56 | ,fClusterPosition(0) | |
57 | ,fHistChargedEnergyRemoved(0) | |
58 | ,fHistNeutralEnergyRemoved(0) | |
59 | ,fHistGammaEnergyAdded(0) | |
2fbf38ac | 60 | { |
61 | ||
62 | } | |
63 | ||
e9da35da | 64 | AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed() |
f61cec2f | 65 | {//destructor |
e9da35da | 66 | delete fCorrections; |
ef647350 | 67 | delete fHistChargedPionEnergyDeposit; /** Energy deposited in calorimeter by charged pions */ |
68 | delete fHistProtonEnergyDeposit; /** Energy deposited in calorimeter by protons */ | |
69 | delete fHistAntiProtonEnergyDeposit; /** Energy deposited in calorimeter by anti-protons */ | |
70 | delete fHistChargedKaonEnergyDeposit; /** Energy deposited in calorimeter by charged kaons */ | |
d0c22dcc | 71 | delete fHistMuonEnergyDeposit; /** Energy deposited in calorimeter by muons */ |
72 | ||
73 | delete fHistRemovedEnergy; // removed energy | |
311c6540 | 74 | delete fClusterPosition; |
75 | delete fHistChargedEnergyRemoved; | |
76 | delete fHistNeutralEnergyRemoved; | |
77 | delete fHistGammaEnergyAdded; | |
d0c22dcc | 78 | |
cf6522d1 | 79 | } |
80 | ||
81 | Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev) | |
e9da35da | 82 | { |
f61cec2f | 83 | |
84 | //AliAnalysisEt::AnalyseEvent(ev); | |
e9da35da | 85 | // analyse ESD event |
2fbf38ac | 86 | ResetEventValues(); |
e9da35da | 87 | if (!ev) { |
88 | AliFatal("ERROR: Event does not exist"); | |
89 | return 0; | |
90 | } | |
91 | ||
2fbf38ac | 92 | AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev); |
e9da35da | 93 | if (!event) { |
94 | AliFatal("ERROR: ESD Event does not exist"); | |
95 | return 0; | |
ec956c46 | 96 | } |
ef647350 | 97 | fSelector->SetEvent(event); |
f61cec2f | 98 | |
743ce29f | 99 | Int_t cent = -1; |
f6b36c54 | 100 | if (fCentrality && cent) |
743ce29f | 101 | { |
102 | cent = fCentrality->GetCentralityClass10("V0M"); | |
ef647350 | 103 | fCentClass = fCentrality->GetCentralityClass10("V0M"); |
743ce29f | 104 | } |
105 | ||
2fbf38ac | 106 | for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++) |
107 | { | |
108 | AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster); | |
109 | if (!cluster) | |
110 | { | |
e9da35da | 111 | AliError(Form("ERROR: Could not get cluster %d", iCluster)); |
2fbf38ac | 112 | continue; |
113 | } | |
ef647350 | 114 | int x = 0; |
f61cec2f | 115 | fCutFlow->Fill(x++); |
86e7d5db | 116 | if(!fSelector->IsDetectorCluster(*cluster)) continue; |
f61cec2f | 117 | fCutFlow->Fill(x++); |
86e7d5db | 118 | if(!fSelector->PassMinEnergyCut(*cluster)) continue; |
f61cec2f | 119 | fCutFlow->Fill(x++); |
86e7d5db | 120 | if (!fSelector->PassDistanceToBadChannelCut(*cluster)) continue; |
f61cec2f | 121 | fCutFlow->Fill(x++); |
e9da35da | 122 | |
2fbf38ac | 123 | Float_t pos[3]; |
e9da35da | 124 | |
2fbf38ac | 125 | cluster->GetPosition(pos); |
e9da35da | 126 | TVector3 cp(pos); |
f61cec2f | 127 | |
e9da35da | 128 | Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex(); |
87efb15c | 129 | |
e9da35da | 130 | Bool_t matched = false; |
131 | ||
f61cec2f | 132 | |
86e7d5db | 133 | matched = !fSelector->PassTrackMatchingCut(*cluster); |
e9da35da | 134 | |
135 | if (matched) | |
136 | { | |
f61cec2f | 137 | |
e9da35da | 138 | if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0) |
139 | { | |
140 | AliVTrack *track = event->GetTrack(trackMatchedIndex); | |
141 | if (!track) { | |
142 | AliError("Error: track does not exist"); | |
143 | } | |
144 | else { | |
145 | const Double_t *pidWeights = track->PID(); | |
146 | ||
147 | Double_t maxpidweight = 0; | |
148 | Int_t maxpid = 0; | |
149 | ||
150 | if (pidWeights) | |
151 | { | |
152 | for (Int_t p =0; p < AliPID::kSPECIES; p++) | |
153 | { | |
154 | if (pidWeights[p] > maxpidweight) | |
155 | { | |
156 | maxpidweight = pidWeights[p]; | |
157 | maxpid = p; | |
158 | } | |
159 | } | |
f61cec2f | 160 | if (fCuts->GetHistMakeTreeDeposit() && fDepositTree) |
e9da35da | 161 | { |
162 | fEnergyDeposited = cluster->E(); | |
f61cec2f | 163 | fMomentumTPC = track->P(); |
e9da35da | 164 | fCharge = track->Charge(); |
165 | fParticlePid = maxpid; | |
166 | fPidProb = maxpidweight; | |
167 | AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track); | |
168 | if (!esdTrack) { | |
169 | AliError("Error: track does not exist"); | |
170 | } | |
171 | else { | |
172 | if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack); | |
f61cec2f | 173 | fDepositTree->Fill(); |
e9da35da | 174 | } |
175 | } | |
176 | ||
177 | if (maxpidweight > fPidCut) | |
178 | { | |
f61cec2f | 179 | //Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]); |
e9da35da | 180 | |
f61cec2f | 181 | //Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2; |
e9da35da | 182 | |
f61cec2f | 183 | //Float_t et = cluster->E() * TMath::Sin(theta); |
e9da35da | 184 | if (maxpid == AliPID::kProton) |
185 | { | |
186 | ||
187 | if (track->Charge() == 1) | |
188 | { | |
e9da35da | 189 | fHistProtonEnergyDeposit->Fill(cluster->E(), track->E()); |
190 | } | |
191 | else if (track->Charge() == -1) | |
192 | { | |
e9da35da | 193 | fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E()); |
194 | } | |
195 | } | |
196 | else if (maxpid == AliPID::kPion) | |
197 | { | |
e9da35da | 198 | fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E()); |
199 | } | |
200 | else if (maxpid == AliPID::kKaon) | |
201 | { | |
e9da35da | 202 | fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E()); |
203 | } | |
204 | else if (maxpid == AliPID::kMuon) | |
205 | { | |
206 | fHistMuonEnergyDeposit->Fill(cluster->E(), track->E()); | |
207 | } | |
208 | } | |
209 | } | |
210 | } | |
211 | } | |
212 | //continue; | |
ba136eb4 | 213 | } // distance |
e9da35da | 214 | else |
215 | { | |
f61cec2f | 216 | fCutFlow->Fill(x++); |
217 | //std::cout << x++ << std::endl; | |
e9da35da | 218 | |
ef647350 | 219 | //if (cluster->E() > fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue; |
220 | //if (cluster->E() < fClusterEnergyCut) continue; | |
e9da35da | 221 | cluster->GetPosition(pos); |
f61cec2f | 222 | |
223 | TVector3 p2(pos); | |
224 | ||
225 | fClusterPosition->Fill(p2.Phi(), p2.PseudoRapidity()); | |
2fbf38ac | 226 | |
b2c10007 | 227 | fTotNeutralEt += CalculateTransverseEnergy(*cluster); |
e9da35da | 228 | fNeutralMultiplicity++; |
229 | } | |
ef647350 | 230 | fMultiplicity++; |
231 | } | |
f61cec2f | 232 | |
ef647350 | 233 | fChargedEnergyRemoved = GetChargedContribution(fNeutralMultiplicity); |
234 | fNeutralEnergyRemoved = GetNeutralContribution(fNeutralMultiplicity); | |
235 | fHistChargedEnergyRemoved->Fill(fChargedEnergyRemoved, fNeutralMultiplicity); | |
236 | fHistNeutralEnergyRemoved->Fill(fNeutralEnergyRemoved, fNeutralMultiplicity); | |
f61cec2f | 237 | |
ef647350 | 238 | fGammaEnergyAdded = GetGammaContribution(fNeutralMultiplicity); |
239 | fHistGammaEnergyAdded->Fill(fGammaEnergyAdded, fNeutralMultiplicity); | |
2fbf38ac | 240 | |
b2c10007 | 241 | Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) + GetGammaContribution(fNeutralMultiplicity) + GetSecondaryContribution(fNeutralMultiplicity); |
e9da35da | 242 | fHistRemovedEnergy->Fill(removedEnergy); |
f61cec2f | 243 | |
ef647350 | 244 | fTotNeutralEt = fGeomCorrection * fEMinCorrection * (fTotNeutralEt - removedEnergy); |
2fbf38ac | 245 | fTotEt = fTotChargedEt + fTotNeutralEt; |
f61cec2f | 246 | // Fill the histograms...0 |
2fbf38ac | 247 | FillHistograms(); |
b2c10007 | 248 | //std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl; |
2fbf38ac | 249 | return 0; |
250 | } | |
251 | ||
252 | bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track) | |
f61cec2f | 253 | { // check vertex |
2fbf38ac | 254 | |
255 | Float_t bxy = 999.; | |
256 | Float_t bz = 999.; | |
e9da35da | 257 | if (!track) { |
258 | AliError("ERROR: no track"); | |
259 | return kFALSE; | |
ec956c46 | 260 | } |
1b8c3d66 | 261 | AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track); |
e9da35da | 262 | if (!esdTrack) { |
263 | AliError("ERROR: no track"); | |
264 | return kFALSE; | |
1b8c3d66 | 265 | } |
266 | esdTrack->GetImpactParametersTPC(bxy,bz); | |
ec956c46 | 267 | |
2fbf38ac | 268 | |
e9da35da | 269 | bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) && |
270 | (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) && | |
271 | (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) && | |
272 | (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) && | |
273 | (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut()); | |
2fbf38ac | 274 | |
4998becf | 275 | return status; |
2fbf38ac | 276 | } |
277 | ||
278 | void AliAnalysisEtReconstructed::Init() | |
f61cec2f | 279 | { // Init |
2fbf38ac | 280 | AliAnalysisEt::Init(); |
83d0f02c | 281 | fPidCut = fCuts->GetReconstructedPidCut(); |
ba136eb4 | 282 | TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG)); |
e9da35da | 283 | if (!fCorrections) { |
284 | cout<<"Warning! You have not set corrections. Your code will crash. You have to set the corrections."<<endl; | |
0f97be4c | 285 | } |
2fbf38ac | 286 | } |
287 | ||
288 | bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField) | |
f61cec2f | 289 | { // propagate track to detector radius |
2fbf38ac | 290 | |
e9da35da | 291 | if (!track) { |
292 | cout<<"Warning: track empty"<<endl; | |
293 | return kFALSE; | |
294 | } | |
295 | AliESDtrack *esdTrack= dynamic_cast<AliESDtrack*>(track); | |
296 | if (!esdTrack) { | |
297 | AliError("ERROR: no ESD track"); | |
298 | return kFALSE; | |
299 | } | |
300 | // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt()); | |
2fbf38ac | 301 | |
302 | Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField); | |
303 | ||
cf6522d1 | 304 | // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt()); |
f61cec2f | 305 | return prop && fSelector->CutGeometricalAcceptance(*esdTrack); |
2fbf38ac | 306 | } |
307 | ||
87efb15c | 308 | void AliAnalysisEtReconstructed::FillOutputList(TList* list) |
f61cec2f | 309 | { // add some extra histograms to the ones from base class |
87efb15c | 310 | AliAnalysisEt::FillOutputList(list); |
311 | ||
312 | list->Add(fHistChargedPionEnergyDeposit); | |
313 | list->Add(fHistProtonEnergyDeposit); | |
314 | list->Add(fHistAntiProtonEnergyDeposit); | |
315 | list->Add(fHistChargedKaonEnergyDeposit); | |
316 | list->Add(fHistMuonEnergyDeposit); | |
e9da35da | 317 | |
318 | list->Add(fHistRemovedEnergy); | |
ef647350 | 319 | list->Add(fClusterPosition); |
f61cec2f | 320 | |
ef647350 | 321 | list->Add(fHistChargedEnergyRemoved); |
322 | list->Add(fHistNeutralEnergyRemoved); | |
323 | list->Add(fHistGammaEnergyAdded); | |
87efb15c | 324 | } |
325 | ||
326 | void AliAnalysisEtReconstructed::CreateHistograms() | |
f61cec2f | 327 | { // add some extra histograms to the ones from base class |
87efb15c | 328 | AliAnalysisEt::CreateHistograms(); |
329 | ||
0fa8c632 | 330 | Int_t nbinsEt = 1000; |
331 | Double_t minEt = 0; | |
332 | Double_t maxEt = 10; | |
333 | ||
334 | // possibly change histogram limits | |
335 | if (fCuts) { | |
e9da35da | 336 | nbinsEt = fCuts->GetHistNbinsParticleEt(); |
337 | minEt = fCuts->GetHistMinParticleEt(); | |
338 | maxEt = fCuts->GetHistMaxParticleEt(); | |
0fa8c632 | 339 | } |
340 | ||
87efb15c | 341 | TString histname; |
342 | histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix; | |
0fa8c632 | 343 | fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt); |
87efb15c | 344 | fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); |
345 | fHistChargedPionEnergyDeposit->SetYTitle("Energy of track"); | |
e9da35da | 346 | |
87efb15c | 347 | histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix; |
0fa8c632 | 348 | fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt); |
87efb15c | 349 | fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); |
350 | fHistProtonEnergyDeposit->SetYTitle("Energy of track"); | |
e9da35da | 351 | |
87efb15c | 352 | histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix; |
0fa8c632 | 353 | fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt); |
87efb15c | 354 | fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); |
355 | fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track"); | |
e9da35da | 356 | |
87efb15c | 357 | histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix; |
0fa8c632 | 358 | fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt); |
87efb15c | 359 | fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); |
360 | fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track"); | |
e9da35da | 361 | |
87efb15c | 362 | histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix; |
0fa8c632 | 363 | fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt); |
87efb15c | 364 | fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); |
365 | fHistMuonEnergyDeposit->SetYTitle("Energy of track"); | |
e9da35da | 366 | |
367 | histname = "fHistRemovedEnergy" + fHistogramNameSuffix; | |
368 | fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20); | |
369 | //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
370 | //fHistMuonEnergyDeposit->SetYTitle("Energy of track"); | |
371 | ||
ef647350 | 372 | histname = "fClusterPosition" + fHistogramNameSuffix; |
373 | fClusterPosition = new TH2D(histname.Data(), "Position of accepted neutral clusters",1000, -2.0, -.5, 1000, -.13 , 0.13); | |
374 | fClusterPosition->SetXTitle("Energy deposited in calorimeter"); | |
375 | fClusterPosition->SetYTitle("Energy of track"); | |
376 | ||
377 | ||
378 | histname = "fHistChargedEnergyRemoved" + fHistogramNameSuffix; | |
379 | fHistChargedEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5); | |
380 | ||
381 | histname = "fHistNeutralEnergyRemoved" + fHistogramNameSuffix; | |
382 | fHistNeutralEnergyRemoved = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5); | |
383 | ||
384 | histname = "fHistGammaEnergyAdded" + fHistogramNameSuffix; | |
385 | fHistGammaEnergyAdded = new TH2D(histname.Data(), histname.Data(), 1000, .0, 30, 100, -0.5 , 99.5); | |
e9da35da | 386 | |
87efb15c | 387 | |
f20da103 | 388 | } |