]>
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" | |
20 | #include "AliVParticle.h" | |
87efb15c | 21 | #include "TDatabasePDG.h" |
22 | #include "TList.h" | |
2fbf38ac | 23 | #include <iostream> |
24 | #include "TH2F.h" | |
25 | ||
26 | AliAnalysisEtReconstructed::AliAnalysisEtReconstructed() : | |
27 | AliAnalysisEt() | |
87efb15c | 28 | ,fNTpcClustersCut(0) |
29 | ,fNItsClustersCut(0) | |
2fbf38ac | 30 | ,fTrackDistanceCut(0) |
87efb15c | 31 | ,fPidCut(0) |
2fbf38ac | 32 | ,fClusterType(0) |
87efb15c | 33 | ,fHistChargedPionEnergyDeposit(0) |
34 | ,fHistProtonEnergyDeposit(0) | |
35 | ,fHistAntiProtonEnergyDeposit(0) | |
36 | ,fHistChargedKaonEnergyDeposit(0) | |
37 | ,fHistMuonEnergyDeposit(0) | |
2fbf38ac | 38 | { |
39 | ||
40 | } | |
41 | ||
cf6522d1 | 42 | AliAnalysisEtReconstructed::~AliAnalysisEtReconstructed() |
2fbf38ac | 43 | { |
cf6522d1 | 44 | } |
45 | ||
46 | Int_t AliAnalysisEtReconstructed::AnalyseEvent(AliVEvent* ev) | |
47 | { // analyse ESD event | |
2fbf38ac | 48 | ResetEventValues(); |
49 | AliESDEvent *event = dynamic_cast<AliESDEvent*>(ev); | |
50 | ||
87efb15c | 51 | Double_t protonMass = fPdgDB->GetParticle("proton")->Mass(); |
52 | ||
2fbf38ac | 53 | for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) |
54 | { | |
55 | AliVParticle *track = event->GetTrack(iTrack); | |
56 | if (!track) | |
57 | { | |
58 | Printf("ERROR: Could not get track %d", iTrack); | |
59 | continue; | |
60 | } | |
61 | ||
62 | fMultiplicity++; | |
63 | ||
64 | Int_t nItsClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(0); | |
65 | Int_t nTPCClusters = dynamic_cast<AliESDtrack*>(track)->GetNcls(1); | |
66 | ||
67 | Float_t massPart = 0; | |
68 | ||
69 | const Double_t *pidWeights = track->PID(); | |
87efb15c | 70 | Int_t maxpid = -1; |
71 | Double_t maxpidweight = 0; | |
72 | ||
2fbf38ac | 73 | if (pidWeights) |
74 | { | |
2fbf38ac | 75 | for (Int_t p =0; p < AliPID::kSPECIES; p++) |
76 | { | |
77 | if (pidWeights[p] > maxpidweight) | |
78 | { | |
79 | maxpidweight = pidWeights[p]; | |
80 | maxpid = p; | |
81 | } | |
82 | } | |
83 | if (maxpid == AliPID::kProton) | |
84 | { | |
87efb15c | 85 | //by definition of ET |
86 | massPart = -protonMass*track->Charge(); | |
2fbf38ac | 87 | } |
88 | ||
89 | } | |
90 | ||
91 | Double_t et = track->E() * TMath::Sin(track->Theta()) + massPart; | |
cf6522d1 | 92 | // 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 |
2fbf38ac | 93 | |
4998becf | 94 | if (TMath::Abs(track->Eta()) < fCuts->GetCommonEtaCut() && CheckGoodVertex(track) && nItsClusters > fCuts->GetReconstructedNItsClustersCut() && nTPCClusters > fCuts->GetReconstructedNTpcClustersCut() ) |
2fbf38ac | 95 | { |
96 | fTotChargedEt += et; | |
97 | fChargedMultiplicity++; | |
87efb15c | 98 | if (maxpid != -1) |
99 | { | |
100 | if (maxpid == AliPID::kPion) | |
101 | { | |
102 | fProtonEt += et; | |
103 | } | |
104 | if (maxpid == AliPID::kKaon) | |
105 | { | |
106 | fChargedKaonEt += et; | |
107 | } | |
108 | if (maxpid == AliPID::kMuon) | |
109 | { | |
110 | fMuonEt += et; | |
111 | } | |
112 | if (maxpid == AliPID::kElectron) | |
113 | { | |
114 | fElectronEt += et; | |
115 | } | |
116 | } | |
2fbf38ac | 117 | |
118 | if (TMath::Abs(track->Eta()) < fEtaCutAcc && track->Phi() < fPhiCutAccMax && track->Phi() > fPhiCutAccMin) | |
119 | { | |
120 | fTotChargedEtAcc += track->E()*TMath::Sin(track->Theta()) + massPart; | |
87efb15c | 121 | if (maxpid != -1) |
122 | { | |
123 | if (maxpid == AliPID::kPion) | |
124 | { | |
125 | fProtonEtAcc += et; | |
126 | } | |
127 | if (maxpid == AliPID::kKaon) | |
128 | { | |
129 | fChargedKaonEtAcc += et; | |
130 | } | |
131 | if (maxpid == AliPID::kMuon) | |
132 | { | |
133 | fMuonEtAcc += et; | |
134 | } | |
135 | if (maxpid == AliPID::kElectron) | |
136 | { | |
137 | fElectronEtAcc += et; | |
138 | } | |
139 | } | |
140 | ||
2fbf38ac | 141 | } |
142 | } | |
143 | ||
2fbf38ac | 144 | if (TrackHitsCalorimeter(track, event->GetMagneticField())) |
145 | { | |
13b0d3c1 | 146 | Double_t phi = track->Phi(); |
147 | Double_t pt = track->Pt(); | |
cf6522d1 | 148 | // printf("Rec track hit: iTrack %03d phi %4.3f pt %4.3f\n", iTrack, phi, pt); // tmp/debug printout |
13b0d3c1 | 149 | if (track->Charge() > 0) fHistPhivsPtPos->Fill(phi, pt); |
150 | else fHistPhivsPtNeg->Fill(phi, pt); | |
2fbf38ac | 151 | } |
152 | } | |
153 | ||
154 | for (Int_t iCluster = 0; iCluster < event->GetNumberOfCaloClusters(); iCluster++) | |
155 | { | |
156 | AliESDCaloCluster* cluster = event->GetCaloCluster(iCluster); | |
157 | if (!cluster) | |
158 | { | |
159 | Printf("ERROR: Could not get cluster %d", iCluster); | |
160 | continue; | |
161 | } | |
162 | ||
87efb15c | 163 | if (cluster->GetType() != fClusterType) continue; |
ba136eb4 | 164 | //if(cluster->GetTracksMatched() > 0) |
87efb15c | 165 | //printf("Rec Cluster: iCluster %03d E %4.3f type %d NCells %d, nmatched: %d, distance to closest: %f\n", iCluster, cluster->E(), (int)(cluster->GetType()), cluster->GetNCells(), cluster->GetNTracksMatched(), cluster->GetEmcCpvDistance()); // tmp/debug printout |
166 | ||
167 | ||
2fbf38ac | 168 | if (cluster->E() < fClusterEnergyCut) continue; |
169 | Float_t pos[3]; | |
170 | TVector3 cp(pos); | |
171 | cluster->GetPosition(pos); | |
87efb15c | 172 | |
ba136eb4 | 173 | Double_t distance = cluster->GetEmcCpvDistance(); |
174 | Int_t trackMatchedIndex = cluster->GetTrackMatchedIndex(); | |
175 | if ( cluster->IsEMCAL() ) { | |
176 | distance = CalcTrackClusterDistance(pos, &trackMatchedIndex, event); | |
177 | } | |
178 | ||
179 | fHistTMDeltaR->Fill(distance); | |
180 | if (distance < fTrackDistanceCut) | |
2fbf38ac | 181 | { |
ba136eb4 | 182 | if (cluster->GetNTracksMatched() == 1 && trackMatchedIndex>0) |
87efb15c | 183 | { |
ba136eb4 | 184 | AliVTrack *track = event->GetTrack(trackMatchedIndex); |
87efb15c | 185 | const Double_t *pidWeights = track->PID(); |
186 | ||
187 | Double_t maxpidweight = 0; | |
188 | Int_t maxpid = 0; | |
189 | ||
190 | if (pidWeights) | |
191 | { | |
192 | for (Int_t p =0; p < AliPID::kSPECIES; p++) | |
193 | { | |
194 | if (pidWeights[p] > maxpidweight) | |
195 | { | |
196 | maxpidweight = pidWeights[p]; | |
197 | maxpid = p; | |
198 | } | |
199 | } | |
200 | if(maxpidweight > fPidCut) | |
201 | { | |
202 | if(maxpid == AliPID::kProton) | |
203 | { | |
204 | if(track->Charge() == 1) | |
205 | { | |
206 | fHistProtonEnergyDeposit->Fill(cluster->E(), track->E()); | |
207 | } | |
208 | else if(track->Charge() == -1) | |
209 | { | |
210 | fHistAntiProtonEnergyDeposit->Fill(cluster->E(), track->E()); | |
211 | } | |
212 | } | |
213 | else if(maxpid == AliPID::kPion) | |
214 | { | |
215 | fHistChargedPionEnergyDeposit->Fill(cluster->E(), track->E()); | |
216 | } | |
217 | else if(maxpid == AliPID::kKaon) | |
218 | { | |
219 | fHistChargedKaonEnergyDeposit->Fill(cluster->E(), track->E()); | |
220 | } | |
221 | else if(maxpid == AliPID::kMuon) | |
222 | { | |
223 | fHistMuonEnergyDeposit->Fill(cluster->E(), track->E()); | |
224 | } | |
225 | } | |
226 | } | |
227 | } | |
228 | ||
2fbf38ac | 229 | continue; |
ba136eb4 | 230 | } // distance |
2fbf38ac | 231 | |
4998becf | 232 | if (cluster->E() > fSingleCellEnergyCut && cluster->GetNCells() == fCuts->GetCommonSingleCell()) continue; |
2fbf38ac | 233 | |
234 | cluster->GetPosition(pos); | |
235 | ||
236 | // TODO: replace with TVector3, too lazy now... | |
237 | ||
238 | float dist = TMath::Sqrt(pos[0]*pos[0] + pos[1]*pos[1]); | |
239 | ||
240 | float theta = TMath::ATan(pos[2]/dist)+TMath::Pi()/2; | |
241 | // float eta = TMath::Log(TMath::Abs( TMath::Tan( 0.5 * theta ) ) ); | |
242 | fTotNeutralEt += cluster->E() * TMath::Sin(theta); | |
13b0d3c1 | 243 | fNeutralMultiplicity++; |
2fbf38ac | 244 | |
245 | fMultiplicity++; | |
2fbf38ac | 246 | } |
247 | ||
248 | fTotNeutralEtAcc = fTotNeutralEt; | |
249 | fTotEt = fTotChargedEt + fTotNeutralEt; | |
250 | fTotEtAcc = fTotChargedEtAcc + fTotNeutralEtAcc; | |
251 | ||
2fbf38ac | 252 | // Fill the histograms... |
253 | FillHistograms(); | |
254 | ||
255 | return 0; | |
256 | } | |
257 | ||
258 | bool AliAnalysisEtReconstructed::CheckGoodVertex(AliVParticle* track) | |
cf6522d1 | 259 | { // check vertex |
2fbf38ac | 260 | |
261 | Float_t bxy = 999.; | |
262 | Float_t bz = 999.; | |
263 | dynamic_cast<AliESDtrack*>(track)->GetImpactParametersTPC(bxy,bz); | |
264 | ||
4998becf | 265 | bool status = (TMath::Abs(track->Xv()) < fCuts->GetReconstructedVertexXCut()) && |
266 | (TMath::Abs(track->Yv()) < fCuts->GetReconstructedVertexYCut()) && | |
267 | (TMath::Abs(track->Zv()) < fCuts->GetReconstructedVertexZCut()) && | |
268 | (TMath::Abs(bxy) < fCuts->GetReconstructedIPxyCut()) && | |
269 | (TMath::Abs(bz) < fCuts->GetReconstructedIPzCut()); | |
2fbf38ac | 270 | |
4998becf | 271 | return status; |
2fbf38ac | 272 | } |
273 | ||
274 | void AliAnalysisEtReconstructed::Init() | |
cf6522d1 | 275 | { // Init |
2fbf38ac | 276 | AliAnalysisEt::Init(); |
87efb15c | 277 | fNItsClustersCut = fCuts->GetReconstructedNItsClustersCut(); |
278 | fNTpcClustersCut = fCuts->GetReconstructedNTpcClustersCut(); | |
83d0f02c | 279 | fPidCut = fCuts->GetReconstructedPidCut(); |
ba136eb4 | 280 | TGeoGlobalMagField::Instance()->SetField(new AliMagF("Maps","Maps", 1., 1., AliMagF::k5kG)); |
281 | ||
2fbf38ac | 282 | } |
283 | ||
284 | bool AliAnalysisEtReconstructed::TrackHitsCalorimeter(AliVParticle* track, Double_t magField) | |
cf6522d1 | 285 | { // propagate track to detector radius |
2fbf38ac | 286 | |
287 | AliESDtrack *esdTrack = dynamic_cast<AliESDtrack*>(track); | |
cf6522d1 | 288 | // Printf("Propagating track: eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt()); |
2fbf38ac | 289 | |
290 | Bool_t prop = esdTrack->PropagateTo(fDetectorRadius, magField); | |
291 | ||
cf6522d1 | 292 | // if (prop) Printf("Track propagated, eta: %f, phi: %f, pt: %f", esdTrack->Eta(), esdTrack->Phi(), esdTrack->Pt()); |
2fbf38ac | 293 | return prop && |
294 | TMath::Abs(esdTrack->Eta()) < fEtaCutAcc && | |
295 | esdTrack->Phi() > fPhiCutAccMin*TMath::Pi()/180. && | |
296 | esdTrack->Phi() < fPhiCutAccMax*TMath::Pi()/180.; | |
297 | } | |
298 | ||
87efb15c | 299 | void AliAnalysisEtReconstructed::FillOutputList(TList* list) |
300 | { | |
301 | AliAnalysisEt::FillOutputList(list); | |
302 | ||
303 | list->Add(fHistChargedPionEnergyDeposit); | |
304 | list->Add(fHistProtonEnergyDeposit); | |
305 | list->Add(fHistAntiProtonEnergyDeposit); | |
306 | list->Add(fHistChargedKaonEnergyDeposit); | |
307 | list->Add(fHistMuonEnergyDeposit); | |
308 | } | |
309 | ||
310 | void AliAnalysisEtReconstructed::CreateHistograms() | |
311 | { | |
312 | AliAnalysisEt::CreateHistograms(); | |
313 | ||
314 | TString histname; | |
315 | histname = "fHistChargedPionEnergyDeposit" + fHistogramNameSuffix; | |
316 | fHistChargedPionEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #pi^{+/-}", 1000, 0, 10, 1000, 0, 10); | |
317 | fHistChargedPionEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
318 | fHistChargedPionEnergyDeposit->SetYTitle("Energy of track"); | |
319 | ||
320 | histname = "fHistProtonEnergyDeposit" + fHistogramNameSuffix; | |
321 | fHistProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by protons", 1000, 0, 10, 1000, 0, 10); | |
322 | fHistProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
323 | fHistProtonEnergyDeposit->SetYTitle("Energy of track"); | |
324 | ||
325 | histname = "fHistAntiProtonEnergyDeposit" + fHistogramNameSuffix; | |
326 | fHistAntiProtonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by anti-protons", 1000, 0, 10, 1000, 0, 10); | |
327 | fHistAntiProtonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
328 | fHistAntiProtonEnergyDeposit->SetYTitle("Energy of track"); | |
329 | ||
330 | histname = "fHistChargedKaonEnergyDeposit" + fHistogramNameSuffix; | |
331 | fHistChargedKaonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by K^{+/-}", 1000, 0, 10, 1000, 0, 10); | |
332 | fHistChargedKaonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
333 | fHistChargedKaonEnergyDeposit->SetYTitle("Energy of track"); | |
334 | ||
335 | histname = "fHistMuonEnergyDeposit" + fHistogramNameSuffix; | |
336 | fHistMuonEnergyDeposit = new TH2F(histname.Data(), "Energy deposited by #mu^{+/-}", 1000, 0, 10, 1000, 0, 10); | |
337 | fHistMuonEnergyDeposit->SetXTitle("Energy deposited in calorimeter"); | |
338 | fHistMuonEnergyDeposit->SetYTitle("Energy of track"); | |
339 | ||
340 | ||
341 | } | |
ba136eb4 | 342 | |
343 | Double_t | |
344 | AliAnalysisEtReconstructed::CalcTrackClusterDistance(Float_t clsPos[3], | |
345 | Int_t *trkMatchId, | |
346 | AliESDEvent *event) | |
347 | { // calculate distance between cluster and closest track | |
348 | ||
349 | Double_t trkPos[3] = {0,0,0}; | |
350 | ||
351 | Int_t bestTrkMatchId = -1; | |
352 | Double_t distance = 9999; // init to a big number | |
353 | ||
354 | Double_t dist = 0; | |
355 | Double_t distX = 0, distY = 0, distZ = 0; | |
356 | ||
357 | for (Int_t iTrack = 0; iTrack < event->GetNumberOfTracks(); iTrack++) { | |
358 | AliESDtrack *track = event->GetTrack(iTrack); | |
359 | if (!track) { | |
360 | Printf("ERROR: Could not get track %d", iTrack); | |
361 | continue; | |
362 | } | |
363 | ||
364 | // check for approx. eta and phi range before we propagate.. | |
365 | // TBD | |
366 | ||
367 | AliEMCALTrack *emctrack = new AliEMCALTrack(*track); | |
368 | if (!emctrack->PropagateToGlobal(clsPos[0],clsPos[1],clsPos[2],0.,0.) ) { | |
369 | continue; | |
370 | } | |
371 | emctrack->GetXYZ(trkPos); | |
372 | if (emctrack) delete emctrack; | |
373 | ||
374 | distX = clsPos[0]-trkPos[0]; | |
375 | distY = clsPos[1]-trkPos[1]; | |
376 | distZ = clsPos[2]-trkPos[2]; | |
377 | dist = TMath::Sqrt(distX*distX + distY*distY + distZ*distZ); | |
378 | ||
379 | if (dist < distance) { | |
380 | distance = dist; | |
381 | bestTrkMatchId = iTrack; | |
382 | } | |
383 | } // iTrack | |
384 | ||
385 | // printf("CalcTrackClusterDistance: bestTrkMatch %d origTrkMatch %d distance %f\n", bestTrkMatchId, *trkMatchId, distance); | |
386 | *trkMatchId = bestTrkMatchId; | |
387 | return distance; | |
388 | } |