]>
Commit | Line | Data |
---|---|---|
1 | //_________________________________________________________________________ | |
2 | // Utility Class for transverse energy studies | |
3 | // Base class for ESD analysis | |
4 | // - reconstruction output | |
5 | // implementation file | |
6 | // | |
7 | //*-- Authors: Oystein Djuvsland (Bergen), David Silvermyr (ORNL) | |
8 | //_________________________________________________________________________ | |
9 | ||
10 | #include "AliAnalysisEtReconstructed.h" | |
11 | #include "AliAnalysisEtCuts.h" | |
12 | #include "AliESDtrack.h" | |
13 | #include "AliEMCALTrack.h" | |
14 | #include "AliESDCaloCluster.h" | |
15 | #include "TVector3.h" | |
16 | #include "TGeoGlobalMagField.h" | |
17 | #include "AliMagF.h" | |
18 | #include "AliVEvent.h" | |
19 | #include "AliESDEvent.h" | |
20 | #include "AliESDtrackCuts.h" | |
21 | #include "AliVParticle.h" | |
22 | #include "TDatabasePDG.h" | |
23 | #include "TList.h" | |
24 | #include "AliESDpid.h" | |
25 | #include <iostream> | |
26 | #include "TH2F.h" | |
27 | #include "AliAnalysisHadEtCorrections.h" | |
28 | #include "AliLog.h" | |
29 | #include "AliCentrality.h" | |
30 | ||
31 | ||
32 | ||
33 | ||
34 | using namespace std; | |
35 | ||
36 | ClassImp(AliAnalysisEtReconstructed); | |
37 | ||
38 | ||
39 | AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() : | |
40 | AliAnalysisEt() | |
41 | ,fCorrections(0) | |
42 | ,fPidCut(0) | |
43 | ,fClusterType(0) | |
44 | ,fHistChargedPionEnergyDeposit(0) | |
45 | ,fHistProtonEnergyDeposit(0) | |
46 | ,fHistAntiProtonEnergyDeposit(0) | |
47 | ,fHistChargedKaonEnergyDeposit(0) | |
48 | ,fHistMuonEnergyDeposit(0) | |
49 | ,fHistRemovedEnergy(0) | |
50 | ,fGeomCorrection(1.0) | |
51 | ,fEMinCorrection(1.0) | |
52 | { | |
53 | ||
54 | } | |
55 | ||
56 | AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed() | |
57 | { | |
58 | delete fCorrections; | |
59 | } | |
60 | ||
61 | Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev) | |
62 | { | |
63 | // analyse ESD event | |
64 | ResetEventValues(); | |
65 | ||
66 | if (!ev) { | |
67 | AliFatal("ERROR: Event does not exist"); | |
68 | return 0; | |
69 | } | |
70 | ||
71 | AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev); | |
72 | if (!event) { | |
73 | AliFatal("ERROR: ESD Event does not exist"); | |
74 | return 0; | |
75 | } | |
76 | ||
77 | Int_t cent = -1; | |
78 | if (fCentrality) | |
79 | { | |
80 | cent = fCentrality->GetCentralityClass10("V0M"); | |
81 | fCentClass = fCentrality->GetCentralityClass10("V0M"); | |
82 | } | |
83 | ||
84 | Double_t protonMass = fgProtonMass; | |
85 | ||
86 | //for PID | |
87 | AliESDpid pID; | |
88 | pID.MakePID(event); | |
89 | TObjArray* list = fEsdtrackCutsTPC->GetAcceptedTracks(event); | |
90 | ||
91 | Int_t nGoodTracks = list->GetEntries(); | |
92 | // printf("nGoodTracks %d nCaloClusters %d\n", nGoodTracks, event->GetNumberOfCaloClusters()); | |
93 | ||
94 | for (Int_t iTrack = 0; iTrack < nGoodTracks; iTrack++) | |
95 | { | |
96 | AliESDtrack *track = dynamic_cast<AliESDtrack*> (list->At(iTrack)); | |
97 | if (!track) | |
98 | { | |
99 | AliError(Form("ERROR: Could not get track %d", iTrack)); | |
100 | continue; | |
101 | } | |
102 | ||
103 | fMultiplicity++; | |
104 | ||
105 | ||
106 | Float_t nSigmaPion,nSigmaProton,nSigmaKaon,nSigmaElectron; | |
107 | nSigmaPion = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kPion)); | |
108 | nSigmaProton = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kProton)); | |
109 | nSigmaKaon = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kKaon)); | |
110 | nSigmaElectron = TMath::Abs(pID.NumberOfSigmasTPC(track,AliPID::kElectron)); | |
111 | /* | |
112 | bool isPion = (nSigmaPion<3.0 && nSigmaProton>2.0 && nSigmaKaon>2.0); | |
113 | bool isElectron = (nSigmaElectron<2.0 && nSigmaPion>4.0 && nSigmaProton>3.0 && nSigmaKaon>3.0); | |
114 | bool isKaon = (nSigmaPion>3.0 && nSigmaProton>2.0 && nSigmaKaon<2.0); | |
115 | bool isProton = (nSigmaPion>3.0 && nSigmaProton<2.0 && nSigmaKaon>2.0); | |
116 | */ | |
117 | ||
118 | Int_t nItsClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(0); | |
119 | Int_t nTPCClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(1); | |
120 | ||
121 | Float_t massPart = 0; | |
122 | ||
123 | const Double_t *pidWeights = track->PID(); | |
124 | Int_t maxpid = -1; | |
125 | Double_t maxpidweight = 0; | |
126 | ||
127 | if (pidWeights) | |
128 | { | |
129 | for (Int_t p =0; p < AliPID::kSPECIES; p++) | |
130 | { | |
131 | if (pidWeights[p] > maxpidweight) | |
132 | { | |
133 | maxpidweight = pidWeights[p]; | |
134 | maxpid = p; | |
135 | } | |
136 | } | |
137 | if (maxpid == AliPID::kProton) | |
138 | { | |
139 | //by definition of ET | |
140 | massPart = -protonMass*track->Charge(); | |
141 | } | |
142 | ||
143 | } | |
144 | ||
145 | Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart; | |
146 | ||
147 | fSparseTracks[0] = maxpid; | |
148 | fSparseTracks[1] = track->Charge(); | |
149 | fSparseTracks[2] = track->M(); | |
150 | fSparseTracks[3] = et; | |
151 | fSparseTracks[4] = track->Pt(); | |
152 | fSparseTracks[5] = track->Eta(); | |
153 | fSparseTracks[6] = cent; | |
154 | fSparseHistTracks->Fill(fSparseTracks); | |
155 | //printf("Rec track: iTrack %03d eta %4.3f phi %4.3f nITSCl %d nTPCCl %d\n", iTrack, track->Eta(), track->Phi(), nItsClusters, nTPCClusters); // tmp/debug printout | |
156 | ||
157 | if (TMath::Abs(track->Eta()) < fCuts->GetCommonEtaCut() && CheckGoodVertex(track) && nItsClusters > fCuts->GetReconstructedNItsClustersCut() && nTPCClusters > fCuts->GetReconstructedNTpcClustersCut() ) | |
158 | { | |
159 | fTotChargedEt += et; | |
160 | fChargedMultiplicity++; | |
161 | if (maxpid != -1) | |
162 | { | |
163 | if (maxpid == AliPID::kProton) | |
164 | { | |
165 | fProtonEt += et; | |
166 | } | |
167 | if (maxpid == AliPID::kPion) | |
168 | { | |
169 | fPionEt += et; | |
170 | } | |
171 | if (maxpid == AliPID::kKaon) | |
172 | { | |
173 | fChargedKaonEt += et; | |
174 | } | |
175 | if (maxpid == AliPID::kMuon) | |
176 | { | |
177 | fMuonEt += et; | |
178 | } | |
179 | if (maxpid == AliPID::kElectron) | |
180 | { | |
181 | fElectronEt += et; | |
182 | } | |
183 | } | |
184 | ||
185 | if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin) | |
186 | { | |
187 | fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart; | |
188 | if (maxpid != -1) | |
189 | { | |
190 | if (maxpid == AliPID::kProton) | |
191 | { | |
192 | fProtonEtAcc += et; | |
193 | } | |
194 | if (maxpid == AliPID::kPion) | |
195 | { | |
196 | fPionEtAcc += et; | |
197 | } | |
198 | if (maxpid == AliPID::kKaon) | |
199 | { | |
200 | fChargedKaonEtAcc += et; | |
201 | } | |
202 | if (maxpid == AliPID::kMuon) | |
203 | { | |
204 | fMuonEtAcc += et; | |
205 | } | |
206 | if (maxpid == AliPID::kElectron) | |
207 | { | |
208 | fElectronEtAcc += et; | |
209 | } | |
210 | } | |
211 | ||
212 | } | |
213 | } | |
214 | ||
215 | if (TrackHitsCalorimeter(track, event->GetMagneticField())) | |
216 | { | |
217 | Double_t phi = track->Phi(); | |
218 | Double_t pt = track->Pt(); | |
219 | // printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout | |
220 | if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt); | |
221 | else fHistPhivsPtNeg->Fill(phi, pt); | |
222 | } | |
223 | } | |
224 | ||
225 | for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++) | |
226 | { | |
227 | AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster); | |
228 | if (!cluster) | |
229 | { | |
230 | AliError(Form("ERROR: Could not get cluster %d", iCluster)); | |
231 | continue; | |
232 | } | |
233 | if (cluster->GetType() != fClusterType) continue; | |
234 | ||
235 | //if(cluster->GetTracksMatched() > 0) | |
236 | // printf("Rec Cluster: iCluster %03d E %4.3f type %.qd NCells %d, nmatched: %d, distance to closest: %f\n", iCluster, cluster->E(), (int)(cluster->GetType()), cluster->GetNCells(), cluster->GetNTracksMatched(), cluster->GetEmcCpvDistance()); // tmp/debug printout | |
237 | ||
238 | ||
239 | if (cluster->E() < fClusterEnergyCut) continue; | |
240 | ||
241 | Float_t pos[3]; | |
242 | ||
243 | cluster->GetPosition(pos); | |
244 | TVector3 cp(pos); | |
245 | ||
246 | Double_t distance = cluster->GetEmcCpvDistance(); | |
247 | Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex(); | |
248 | if ( cluster->IsEMCAL() ) { | |
249 | distance = CalcTrackClusterDistance(pos, &trackMatchedIndex, event); | |
250 | } | |
251 | ||
252 | fSparseClusters[0] = 0; | |
253 | fSparseClusters[1] = 0; | |
254 | fSparseClusters[2] = 0; | |
255 | fSparseClusters[6] = 0; | |
256 | fSparseClusters[7] = 0; | |
257 | fSparseClusters[8] = 0; | |
258 | fSparseClusters[9] = cent; | |
259 | fSparseClusters[10] = 0; | |
260 | ||
261 | ||
262 | if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex > -1) | |
263 | { | |
264 | AliVTrack *tmptrack = event->GetTrack(trackMatchedIndex); | |
265 | if (!tmptrack) | |
266 | { | |
267 | AliError("Error: track does not exist"); | |
268 | return -1; | |
269 | } | |
270 | const Double_t *pidWeights = tmptrack->PID(); | |
271 | ||
272 | Double_t maxpidweight = 0; | |
273 | Int_t maxpid = 0; | |
274 | Double_t massPart = 0; | |
275 | if (pidWeights) | |
276 | { | |
277 | for (Int_t p =0; p < AliPID::kSPECIES; p++) | |
278 | { | |
279 | if (pidWeights[p] > maxpidweight) | |
280 | { | |
281 | maxpidweight = pidWeights[p]; | |
282 | maxpid = p; | |
283 | } | |
284 | } | |
285 | if (maxpid == AliPID::kProton) | |
286 | { | |
287 | //by definition of ET | |
288 | massPart = -protonMass*tmptrack->Charge(); | |
289 | } | |
290 | } | |
291 | fSparseClusters[0] = maxpid; | |
292 | fSparseClusters[1] = tmptrack->Charge(); | |
293 | fSparseClusters[2] = tmptrack->M(); | |
294 | fSparseClusters[6] = tmptrack->E() * TMath::Sin(tmptrack->Theta()) + massPart;; | |
295 | fSparseClusters[7] = tmptrack->Pt(); | |
296 | fSparseClusters[8] = tmptrack->Eta(); | |
297 | } | |
298 | ||
299 | fSparseClusters[10] = distance; | |
300 | ||
301 | fHistTMDeltaR->Fill(distance); | |
302 | fHistTMDxDz->Fill(cluster->GetTrackDx(), cluster->GetTrackDz()); | |
303 | ||
304 | // Float_t clusteret = cluster->E() * TMath::Sin(cp.Theta()); | |
305 | ||
306 | Bool_t matched = false; | |
307 | ||
308 | if (cluster->IsEMCAL()) matched = distance < fTrackDistanceCut; | |
309 | else matched = (TMath::Abs(cluster->GetTrackDx()) < fTrackDxCut && TMath::Abs(cluster->GetTrackDz()) < fTrackDzCut); | |
310 | ||
311 | if (matched) | |
312 | { | |
313 | if (cluster->GetNTracksMatched() > 0 && trackMatchedIndex>=0) | |
314 | { | |
315 | AliVTrack *track = event->GetTrack(trackMatchedIndex); | |
316 | if (!track) { | |
317 | AliError("Error: track does not exist"); | |
318 | } | |
319 | else { | |
320 | const Double_t *pidWeights = track->PID(); | |
321 | ||
322 | Double_t maxpidweight = 0; | |
323 | Int_t maxpid = 0; | |
324 | ||
325 | if (pidWeights) | |
326 | { | |
327 | for (Int_t p =0; p < AliPID::kSPECIES; p++) | |
328 | { | |
329 | if (pidWeights[p] > maxpidweight) | |
330 | { | |
331 | maxpidweight = pidWeights[p]; | |
332 | maxpid = p; | |
333 | } | |
334 | } | |
335 | if (fCuts->GetHistMakeTreeDeposit() && fTreeDeposit) | |
336 | { | |
337 | fEnergyDeposited = cluster->E(); | |
338 | fEnergyTPC = track->E(); | |
339 | fCharge = track->Charge(); | |
340 | fParticlePid = maxpid; | |
341 | fPidProb = maxpidweight; | |
342 | AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track); | |
343 | if (!esdTrack) { | |
344 | AliError("Error: track does not exist"); | |
345 | } | |
346 | else { | |
347 | if (esdTrack) fTrackPassedCut = fEsdtrackCutsTPC->AcceptTrack(esdTrack); | |
348 | fTreeDeposit->Fill(); | |
349 | } | |
350 | } | |
351 | ||
352 | if (maxpidweight > fPidCut) | |
353 | { | |
354 | Float_t dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]); | |
355 | ||
356 | Float_t theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2; | |
357 | ||
358 | Float_t et = cluster->E() * TMath::Sin(theta); | |
359 | if (maxpid == AliPID::kProton) | |
360 | { | |
361 | ||
362 | if (track->Charge() == 1) | |
363 | { | |
364 | fBaryonEt += et; | |
365 | fHistProtonEnergyDeposit->Fill(cluster->E(), track->E()); | |
366 | } | |
367 | else if (track->Charge() == -1) | |
368 | { | |
369 | fAntiBaryonEt += et; | |
370 | fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E()); | |
371 | } | |
372 | } | |
373 | else if (maxpid == AliPID::kPion) | |
374 | { | |
375 | fMesonEt += et; | |
376 | fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E()); | |
377 | } | |
378 | else if (maxpid == AliPID::kKaon) | |
379 | { | |
380 | fMesonEt += et; | |
381 | fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E()); | |
382 | } | |
383 | else if (maxpid == AliPID::kMuon) | |
384 | { | |
385 | fHistMuonEnergyDeposit->Fill(cluster->E(), track->E()); | |
386 | } | |
387 | } | |
388 | } | |
389 | } | |
390 | } | |
391 | //continue; | |
392 | } // distance | |
393 | else | |
394 | { | |
395 | fSparseClusters[0] = AliPID::kPhoton; | |
396 | fSparseClusters[1] = 0; | |
397 | ||
398 | if (cluster->E() > fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue; | |
399 | if (cluster->E() < fClusterEnergyCut) continue; | |
400 | cluster->GetPosition(pos); | |
401 | ||
402 | // TODO: replace with TVector3, too lazy now... | |
403 | ||
404 | float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]); | |
405 | ||
406 | float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2; | |
407 | // float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) ); | |
408 | fTotNeutralEt += cluster->E() * TMath::Sin(theta); | |
409 | fNeutralMultiplicity++; | |
410 | } | |
411 | ||
412 | cluster->GetPosition(pos); | |
413 | ||
414 | // TODO: replace with TVector3 | |
415 | ||
416 | float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]); | |
417 | float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2; | |
418 | float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) ); | |
419 | fSparseClusters[3] = cluster->E() * TMath::Sin(theta); | |
420 | fSparseClusters[4] = cluster->E(); | |
421 | fSparseClusters[5] = eta; | |
422 | ||
423 | fSparseHistClusters->Fill(fSparseClusters); | |
424 | ||
425 | fMultiplicity++; | |
426 | } | |
427 | Double_t removedEnergy = GetChargedContribution(fNeutralMultiplicity) + GetNeutralContribution(fNeutralMultiplicity) - GetGammaContribution(fNeutralMultiplicity); | |
428 | fHistRemovedEnergy->Fill(removedEnergy); | |
429 | // std::cout << "fTotNeutralEt: " << fTotNeutralEt << ", Contribution from non-removed charged: " << GetChargedContribution(fNeutralMultiplicity) << ", neutral: " << GetNeutralContribution(fNeutralMultiplicity) << ", gammas: " << GetGammaContribution(fNeutralMultiplicity) << ", multiplicity: " << fNeutralMultiplicity<< std::endl; | |
430 | fTotNeutralEt = fGeomCorrection * fEMinCorrection * fTotNeutralEt - removedEnergy; | |
431 | fTotNeutralEtAcc = fTotNeutralEt/fGeomCorrection; | |
432 | fTotEt = fTotChargedEt + fTotNeutralEt; | |
433 | fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc; | |
434 | fSparseEt[0] = fTotEt; | |
435 | fSparseEt[1] = fTotNeutralEt; | |
436 | fSparseEt[2] = fTotChargedEtAcc; | |
437 | fSparseEt[3] = fMultiplicity; | |
438 | fSparseEt[4] = fNeutralMultiplicity; | |
439 | fSparseEt[5] = fChargedMultiplicity; | |
440 | fSparseEt[6] = cent; | |
441 | // Fill the histograms... | |
442 | FillHistograms(); | |
443 | ||
444 | return 0; | |
445 | } | |
446 | ||
447 | bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track) | |
448 | { // check vertex | |
449 | ||
450 | Float_t bxy = 999.; | |
451 | Float_t bz = 999.; | |
452 | if (!track) { | |
453 | AliError("ERROR: no track"); | |
454 | return kFALSE; | |
455 | } | |
456 | AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track); | |
457 | if (!esdTrack) { | |
458 | AliError("ERROR: no track"); | |
459 | return kFALSE; | |
460 | } | |
461 | esdTrack->GetImpactParametersTPC(bxy,bz); | |
462 | ||
463 | ||
464 | bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) && | |
465 | (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) && | |
466 | (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) && | |
467 | (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) && | |
468 | (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut()); | |
469 | ||
470 | return status; | |
471 | } | |
472 | ||
473 | void AliAnalysisEtReconstructed::Init() | |
474 | { // Init | |
475 | AliAnalysisEt::Init(); | |
476 | fPidCut = fCuts->GetReconstructedPidCut(); | |
477 | TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG)); | |
478 | if (!fCorrections) { | |
479 | cout<<"Warning! You have not set corrections. Your code will crash. You have to set the corrections."<<endl; | |
480 | } | |
481 | } | |
482 | ||
483 | bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField) | |
484 | { // propagate track to detector radius | |
485 | ||
486 | if (!track) { | |
487 | cout<<"Warning: track empty"<<endl; | |
488 | return kFALSE; | |
489 | } | |
490 | AliESDtrack *esdTrack= dynamic_cast<AliESDtrack*>(track); | |
491 | if (!esdTrack) { | |
492 | AliError("ERROR: no ESD track"); | |
493 | return kFALSE; | |
494 | } | |
495 | // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt()); | |
496 | ||
497 | Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField); | |
498 | ||
499 | // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt()); | |
500 | return prop && | |
501 | TMath::Abs(esdTrack->Eta()) < fEtaCutAcc && | |
502 | esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. && | |
503 | esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.; | |
504 | } | |
505 | ||
506 | void AliAnalysisEtReconstructed::FillOutputList(TList* list) | |
507 | { // add some extra histograms to the ones from base class | |
508 | AliAnalysisEt::FillOutputList(list); | |
509 | ||
510 | list->Add(fHistChargedPionEnergyDeposit); | |
511 | list->Add(fHistProtonEnergyDeposit); | |
512 | list->Add(fHistAntiProtonEnergyDeposit); | |
513 | list->Add(fHistChargedKaonEnergyDeposit); | |
514 | list->Add(fHistMuonEnergyDeposit); | |
515 | ||
516 | list->Add(fHistRemovedEnergy); | |
517 | } | |
518 | ||
519 | void AliAnalysisEtReconstructed::CreateHistograms() | |
520 | { // add some extra histograms to the ones from base class | |
521 | AliAnalysisEt::CreateHistograms(); | |
522 | ||
523 | Int_t nbinsEt = 1000; | |
524 | Double_t minEt = 0; | |
525 | Double_t maxEt = 10; | |
526 | ||
527 | // possibly change histogram limits | |
528 | if (fCuts) { | |
529 | nbinsEt = fCuts->GetHistNbinsParticleEt(); | |
530 | minEt = fCuts->GetHistMinParticleEt(); | |
531 | maxEt = fCuts->GetHistMaxParticleEt(); | |
532 | } | |
533 | ||
534 | TString histname; | |
535 | histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix; | |
536 | fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt); | |
537 | fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
538 | fHistChargedPionEnergyDeposit->SetYTitle("Energy of track"); | |
539 | ||
540 | histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix; | |
541 | fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt); | |
542 | fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
543 | fHistProtonEnergyDeposit->SetYTitle("Energy of track"); | |
544 | ||
545 | histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix; | |
546 | fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt); | |
547 | fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
548 | fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track"); | |
549 | ||
550 | histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix; | |
551 | fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt); | |
552 | fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
553 | fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track"); | |
554 | ||
555 | histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix; | |
556 | fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", nbinsEt, minEt, maxEt, nbinsEt, minEt, maxEt); | |
557 | fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
558 | fHistMuonEnergyDeposit->SetYTitle("Energy of track"); | |
559 | ||
560 | histname = "fHistRemovedEnergy" + fHistogramNameSuffix; | |
561 | fHistRemovedEnergy = new TH1F(histname.Data(), histname.Data(), 1000, 0, 20); | |
562 | //fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
563 | //fHistMuonEnergyDeposit->SetYTitle("Energy of track"); | |
564 | ||
565 | ||
566 | ||
567 | } | |
568 | ||
569 | Double_t | |
570 | AliAnalysisEtReconstructed::CalcTrackClusterDistance(const Float_t clsPos[3], | |
571 | Int_t *trkMatchId, | |
572 | const AliESDEvent *event) | |
573 | { // calculate distance between cluster and closest track | |
574 | ||
575 | Double_t trkPos[3] = {0,0,0}; | |
576 | ||
577 | Int_t bestTrkMatchId = -1; | |
578 | Double_t distance = 9999; // init to a big number | |
579 | ||
580 | Double_t dist = 0; | |
581 | Double_t distX = 0, distY = 0, distZ = 0; | |
582 | ||
583 | for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) { | |
584 | AliESDtrack *track = event->GetTrack(iTrack); | |
585 | if (!track) { | |
586 | AliError(Form("ERROR: Could not get track %d", iTrack)); | |
587 | continue; | |
588 | } | |
589 | ||
590 | // check for approx. eta and phi range before we propagate.. | |
591 | // TBD | |
592 | ||
593 | AliEMCALTrack *emctrack = new AliEMCALTrack(*track); | |
594 | if (!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) { | |
595 | continue; | |
596 | } | |
597 | emctrack->GetXYZ(trkPos); | |
598 | delete emctrack; | |
599 | ||
600 | distX = clsPos[0]-trkPos[0]; | |
601 | distY = clsPos[1]-trkPos[1]; | |
602 | distZ = clsPos[2]-trkPos[2]; | |
603 | dist = TMath::Sqrt(distX*distX + distY*distY + distZ*distZ); | |
604 | ||
605 | if (dist < distance) { | |
606 | distance = dist; | |
607 | bestTrkMatchId = iTrack; | |
608 | } | |
609 | } // iTrack | |
610 | ||
611 | // printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance); | |
612 | *trkMatchId = bestTrkMatchId; | |
613 | return distance; | |
614 | } | |
615 |