]>
Commit | Line | Data |
---|---|---|
0e676e3f | 1 | // Author: Artur Szostak |
2 | // Email: artur@alice.phy.uct.ac.za | artursz@iafrica.com | |
3 | ||
4 | #include "AliHLTMUONTracker.h" | |
5 | #include "Tracking/MansoTracker.hpp" | |
6 | #include "AliMUONDataInterface.h" | |
7 | #include "AliRunLoader.h" | |
8 | #include "AliLog.h" | |
9 | #include "AliESD.h" | |
10 | #include "AliESDMuonTrack.h" | |
11 | ||
12 | ClassImp(AliHLTMUONTracker) | |
13 | ||
14 | ||
15 | AliHLTMUONTracker::AliHLTMUONTracker(AliRunLoader* runloader) : AliTracker() | |
16 | { | |
17 | // Creates the the MicrodHLT object and its associated data source and sink | |
18 | // objects. The MicrodHLT object is then initialised by hooking to these objects. | |
19 | ||
20 | AliDebug(2, Form("Called for object 0x%X and with runloader = 0x%X", (ULong_t)this, (ULong_t)runloader)); | |
21 | ||
22 | // Create the dHLT objects. | |
23 | fdHLT = new AliMUONHLT::MicrodHLT(); | |
24 | fTriggers = new AliMUONHLT::TriggerSource(); | |
25 | fClusters = new AliMUONHLT::ClusterSource(); | |
26 | fTracks = new AliMUONHLT::TrackSink(); | |
27 | ||
28 | // Hook up all the objects. | |
29 | fdHLT->SetTriggerSource(fTriggers); | |
30 | fdHLT->SetClusterSource(fClusters); | |
31 | fdHLT->SetTrackSink(fTracks); | |
32 | } | |
33 | ||
34 | ||
35 | AliHLTMUONTracker::~AliHLTMUONTracker() | |
36 | { | |
37 | // Deletes all dHLT objects created in the constructor. | |
38 | ||
39 | AliDebug(2, Form("Called for object 0x%X", (ULong_t)this)); | |
40 | delete fTracks; | |
41 | delete fClusters; | |
42 | delete fTriggers; | |
43 | delete fdHLT; | |
44 | } | |
45 | ||
46 | ||
47 | Int_t AliHLTMUONTracker::LoadClusters(TTree* data) | |
48 | { | |
49 | // Fills the trigger and cluster data source from the current runloader. | |
50 | // The runloader must be open and initialised properly. | |
51 | ||
52 | AliDebug(1, Form("Called for object 0x%X and with data tree = 0x%X", (ULong_t)this, (ULong_t)data)); | |
53 | ||
54 | // Initialise the MUON data interface to use current runloader. | |
55 | AliMUONDataInterface di; | |
56 | di.UseCurrentRunLoader(); | |
57 | AliDebug(2, Form("Loading for event %d", di.CurrentEvent())); | |
58 | ||
59 | // Load the trigger records. | |
d6576e27 | 60 | fTriggers->DataToUse(AliMUONHLT::TriggerSource::FromLocalTriggers); |
0e676e3f | 61 | fTriggers->FillFrom(&di, di.CurrentEvent()); |
62 | ||
63 | #ifndef LOG_NO_DEBUG | |
64 | TString str = "====================== Loaded Triggers ========================\nTrig#\tSign\tPt\t\tX1\t\tY1\t\tX2\t\tY2\n"; | |
65 | for (fTriggers->GetFirstEvent(); fTriggers->MoreEvents(); fTriggers->GetNextEvent()) | |
66 | for (fTriggers->GetFirstBlock(); fTriggers->MoreBlocks(); fTriggers->GetNextBlock()) | |
67 | for (fTriggers->GetFirstTrigger(); fTriggers->MoreTriggers(); fTriggers->GetNextTrigger()) | |
68 | { | |
69 | const AliMUONHLT::TriggerRecord* trig = fTriggers->GetTrigger(); | |
70 | str += trig->TriggerNumber(); | |
71 | str += "\t"; | |
72 | str += trig->ParticleSign(); | |
73 | str += "\t"; | |
74 | str += trig->Pt(); | |
75 | str += "\t"; | |
76 | str += trig->Station1Point().fX; | |
77 | str += "\t"; | |
78 | str += trig->Station1Point().fY; | |
79 | str += "\t"; | |
80 | str += trig->Station2Point().fX; | |
81 | str += "\t"; | |
82 | str += trig->Station2Point().fY; | |
83 | str += "\n"; | |
84 | } | |
85 | AliDebug(5, str); | |
86 | #endif // LOG_NO_DEBUG | |
87 | ||
88 | // Load cluster points (reconstructed hits) | |
d6576e27 | 89 | fClusters->DataToUse(AliMUONHLT::ClusterSource::FromRawClusters); |
0e676e3f | 90 | fClusters->FillFrom(&di, di.CurrentEvent()); |
91 | ||
92 | #ifndef LOG_NO_DEBUG | |
93 | str = "====================== Loaded Clusters ========================\nChamber\tX\t\tY\n"; | |
94 | for (fClusters->GetFirstEvent(); fClusters->MoreEvents(); fClusters->GetNextEvent()) | |
95 | for (fClusters->GetFirstBlock(); fClusters->MoreBlocks(); fClusters->GetNextBlock()) | |
96 | for (fClusters->GetFirstCluster(); fClusters->MoreClusters(); fClusters->GetNextCluster()) | |
97 | { | |
98 | Float_t x, y; | |
99 | fClusters->FetchCluster(x, y); | |
100 | str += fClusters->Chamber(); | |
101 | str += "\t"; | |
102 | str += x; | |
103 | str += "\t"; | |
104 | str += y; | |
105 | str += "\n"; | |
106 | } | |
107 | AliDebug(5, str); | |
108 | #endif // LOG_NO_DEBUG | |
109 | ||
110 | return 0; | |
111 | } | |
112 | ||
113 | ||
114 | void AliHLTMUONTracker::UnloadClusters() | |
115 | { | |
116 | // Frees the triggers and clusters loaded in LoadClusters(). | |
117 | ||
118 | AliDebug(1, Form("Called for object 0x%X", (ULong_t)this)); | |
119 | ||
120 | // Release internal arrays. | |
121 | fTriggers->Clear(); | |
122 | fClusters->Clear(); | |
123 | ||
124 | return; | |
125 | } | |
126 | ||
127 | ||
128 | const AliMUONHLT::TriggerRecord* | |
129 | AliHLTMUONTracker::FindTriggerRecord(const AliMUONHLT::Track* track) const | |
130 | { | |
131 | // Finds the corresponding trigger record object for the given track. | |
132 | // This is doen by matching up the trigger record number with the trigger | |
133 | // record ID number in the track. | |
134 | ||
135 | fTriggers->GetEvent( fTracks->CurrentEvent() ); | |
136 | for (fTriggers->GetFirstBlock(); fTriggers->MoreBlocks(); fTriggers->GetNextBlock()) | |
137 | for (fTriggers->GetFirstTrigger(); fTriggers->MoreTriggers(); fTriggers->GetNextTrigger()) | |
138 | { | |
139 | const AliMUONHLT::TriggerRecord* trigrec = fTriggers->GetTrigger(); | |
140 | if ( trigrec->TriggerNumber() == track->TriggerID() ) | |
141 | return trigrec; | |
142 | } | |
143 | return NULL; | |
144 | } | |
145 | ||
146 | ||
147 | void AliHLTMUONTracker::LeastSquaresFit(const Double_t x[4], const Double_t y[4], Double_t& m, Double_t& c) const | |
148 | { | |
149 | // Least squares fit for a 4 point line: y = m*x + c | |
150 | ||
151 | Double_t n = 4.; // Number of data points. | |
152 | ||
153 | // Compute the sums. | |
154 | Double_t sumx, sumy, sumxx, sumxy; | |
155 | sumx = sumy = sumxx = sumxy = 0.; | |
156 | for (Int_t i = 0; i < 4; i++) | |
157 | { | |
158 | sumx += x[i]; | |
159 | sumy += y[i]; | |
160 | sumxx += x[i] * x[i]; | |
161 | sumxy += x[i] * y[i]; | |
162 | } | |
163 | // Now compute the denominator then m and c parameters. | |
164 | Double_t denom = n*sumxx - sumx*sumx; | |
165 | m = (n*sumxy - sumx*sumy) / denom; | |
166 | c = (sumy*sumxx - sumx*sumxy) / denom; | |
167 | } | |
168 | ||
169 | ||
170 | Double_t AliHLTMUONTracker::ComputeChi2(const AliMUONHLT::Track* track) const | |
171 | { | |
172 | // Computes the Chi^2 for the line fit of the found track fragment. | |
173 | ||
174 | const AliMUONHLT::TriggerRecord* trigrec = FindTriggerRecord(track); | |
175 | if (trigrec == NULL) return -1.; | |
176 | AliDebug(10, Form("Found trigger #%d, particle sign = %d, pt = %f", | |
177 | trigrec->TriggerNumber(), | |
178 | trigrec->ParticleSign(), | |
179 | trigrec->Pt() | |
180 | ) | |
181 | ); | |
182 | ||
183 | // Initialise X, Y and Z coordinate arrays. | |
184 | Double_t st4x, st4y, st4z, st5x, st5y, st5z; | |
185 | if (track->Hit(6).fX == 0. && track->Hit(6).fY == 0.) | |
186 | { | |
187 | st4x = track->Hit(7).fX; | |
188 | st4y = track->Hit(7).fY; | |
189 | st4z = dHLT::Tracking::MansoTracker::GetZ8(); | |
190 | } | |
191 | else | |
192 | { | |
193 | st4x = track->Hit(6).fX; | |
194 | st4y = track->Hit(6).fY; | |
195 | st4z = dHLT::Tracking::MansoTracker::GetZ7(); | |
196 | } | |
197 | if (track->Hit(8).fX == 0. && track->Hit(8).fY == 0.) | |
198 | { | |
199 | st5x = track->Hit(9).fX; | |
200 | st5y = track->Hit(9).fY; | |
201 | st5z = dHLT::Tracking::MansoTracker::GetZ10(); | |
202 | } | |
203 | else | |
204 | { | |
205 | st5x = track->Hit(8).fX; | |
206 | st5y = track->Hit(8).fY; | |
207 | st5z = dHLT::Tracking::MansoTracker::GetZ9(); | |
208 | } | |
209 | ||
210 | Double_t x[4] = {st4x, st5x, | |
211 | trigrec->Station1Point().fX, trigrec->Station2Point().fX | |
212 | }; | |
213 | Double_t y[4] = {st4y, st5y, | |
214 | trigrec->Station1Point().fY, trigrec->Station2Point().fY | |
215 | }; | |
216 | Double_t z[4] = {st4z, st5z, | |
217 | dHLT::Tracking::MansoTracker::GetZ11(), | |
218 | dHLT::Tracking::MansoTracker::GetZ13() | |
219 | }; | |
220 | ||
221 | // Fit a line to the x, y, z data points. | |
222 | Double_t mx, cx; // for x = mx*z + cx | |
223 | Double_t my, cy; // for y = my*z + cy | |
224 | LeastSquaresFit(z, x, mx, cx); | |
225 | AliDebug(11, Form("Fitted line to xz plane: m = %f c = %f", mx, cx)); | |
226 | LeastSquaresFit(z, y, my, cy); | |
227 | AliDebug(11, Form("Fitted line to yz plane: m = %f c = %f", my, cy)); | |
228 | ||
229 | // Now compute the residual, square it and add it to the total chi2. | |
230 | Double_t chi2 = 0.; | |
231 | for (Int_t i = 0; i < 4; i++) | |
232 | { | |
233 | Double_t rx = x[i] - (mx*z[i] + cx); | |
234 | Double_t ry = y[i] - (my*z[i] + cy); | |
235 | AliDebug(15, Form("Real x = %f, fitted x = %f for z = %f", x[i], (mx*z[i] + cx), z[i])); | |
236 | AliDebug(15, Form("Real y = %f, fitted y = %f for z = %f", y[i], (my*z[i] + cy), z[i])); | |
237 | chi2 += rx*rx + ry*ry; | |
238 | AliDebug(15, Form("Running total for Chi2 = %f", chi2)); | |
239 | } | |
240 | ||
241 | return chi2; | |
242 | } | |
243 | ||
244 | ||
245 | Int_t AliHLTMUONTracker::Clusters2Tracks(AliESD* event) | |
246 | { | |
247 | // Runs the dHLT tracking algorithm via the MicrodHLT fdHLT object. The found tracks are | |
248 | // filled by MicrodHLT in fTracks, which then need to be unpacked into the AliRoot ESD. | |
249 | ||
250 | AliDebug(1, Form("Called for object 0x%X and ESD object = 0x%X", (ULong_t)this, (ULong_t)event)); | |
251 | ||
252 | // Run the dHLT tracking algorithm. | |
253 | fdHLT->Run(); | |
254 | AliDebug(1, "Finished running dHLT."); | |
255 | ||
256 | // Unpack the output and fill the ESD. | |
257 | #ifndef LOG_NO_DEBUG | |
258 | TString str = "====================== Found Tracks ========================\nID\tSign\tP\t\tPt\t\tChamber\tX\t\tY\n"; | |
259 | #endif // LOG_NO_DEBUG | |
260 | for (fTracks->GetFirstEvent(); fTracks->MoreEvents(); fTracks->GetNextEvent()) | |
261 | for (fTracks->GetFirstBlock(); fTracks->MoreBlocks(); fTracks->GetNextBlock()) | |
262 | for (fTracks->GetFirstTrack(); fTracks->MoreTracks(); fTracks->GetNextTrack()) | |
263 | { | |
264 | const AliMUONHLT::Track* track = fTracks->GetTrack(); | |
265 | ||
266 | #ifndef LOG_NO_DEBUG | |
267 | // Build some debug logging information. | |
268 | str += track->TriggerID(); | |
269 | str += "\t"; | |
270 | str += track->ParticleSign(); | |
271 | str += "\t"; | |
272 | str += track->P(); | |
273 | str += "\t"; | |
274 | str += track->Pt(); | |
275 | str += "\t7\t"; | |
276 | str += track->Hit(6).fX; | |
277 | str += "\t"; | |
278 | str += track->Hit(6).fY; | |
279 | str += "\n\t\t\t\t\t\t8\t"; | |
280 | str += track->Hit(7).fX; | |
281 | str += "\t"; | |
282 | str += track->Hit(7).fY; | |
283 | str += "\n\t\t\t\t\t\t9\t"; | |
284 | str += track->Hit(8).fX; | |
285 | str += "\t"; | |
286 | str += track->Hit(8).fY; | |
287 | str += "\n\t\t\t\t\t\t10\t"; | |
288 | str += track->Hit(9).fX; | |
289 | str += "\t"; | |
290 | str += track->Hit(9).fY; | |
291 | str += "\n"; | |
292 | #endif // LOG_NO_DEBUG | |
293 | ||
294 | Double_t z = dHLT::Tracking::MansoTracker::GetZ7(); | |
295 | Double_t x = track->Hit(6).fX; | |
296 | if (track->Hit(6).fX == 0. && track->Hit(6).fY == 0.) | |
297 | { | |
298 | z = dHLT::Tracking::MansoTracker::GetZ8(); | |
299 | x = track->Hit(7).fX; | |
300 | }; | |
301 | ||
302 | Double_t p = track->P(); | |
303 | Double_t pt = track->Pt(); | |
304 | ||
305 | // Using the approximation: px / pz = x / z , and pz^2 = p^2 - pt^2 | |
306 | Double_t pz2 = TMath::Abs(p*p - pt*pt); | |
307 | Double_t px2 = -1.; | |
308 | if (z != 0.) | |
309 | px2 = pz2 * x*x / (z*z); | |
310 | Double_t pyz = -1.; | |
311 | if (p*p - px2 >= 0.) | |
312 | pyz = TMath::Sqrt(p*p - px2); | |
313 | Double_t px = -1.; | |
314 | if (px2 >= 0.) | |
315 | px = TMath::Sqrt(px2); | |
316 | Double_t py = -1.; | |
317 | if (pt*pt - px2 >= 0.) | |
318 | py = TMath::Sqrt(pt*pt - px2); // pt^2 = px^2 + py^2 | |
319 | Double_t pz = -1.; | |
320 | if (pz2 >= 0.) | |
321 | pz = TMath::Sqrt(pz2); | |
322 | ||
323 | Double_t invmom = -1.0; | |
324 | if (pyz != 0.) | |
325 | invmom = 1. / pyz * track->ParticleSign(); | |
326 | Double_t thetax = TMath::ATan2(px, pz); | |
327 | Double_t thetay = TMath::ATan2(py, pz); | |
328 | ||
329 | Double_t chi2 = ComputeChi2(track); | |
330 | ||
331 | // Create MUON track, fill it and add it to the ESD. | |
332 | AliESDMuonTrack mt; | |
333 | mt.SetInverseBendingMomentum(invmom); | |
334 | mt.SetThetaX(thetax); | |
335 | mt.SetThetaY(thetay); | |
336 | ||
337 | // Our algorithm assumes the particle originates from the origin. | |
338 | mt.SetZ(0.0); | |
339 | mt.SetBendingCoor(0.0); | |
340 | mt.SetNonBendingCoor(0.0); | |
341 | ||
342 | mt.SetChi2(chi2); | |
343 | mt.SetNHit(2); // Always 2, thats the way Manso's algorithm works. | |
344 | mt.SetMatchTrigger( chi2 != -1. ? 1 : 0 ); | |
345 | mt.SetChi2MatchTrigger(chi2); | |
346 | ||
347 | AliDebug(2, Form( | |
348 | "Adding AliESDMuonTrack with inverse bending momentum" | |
349 | " = %f, theta X = %f, theta Y = %f, chi2 = %f", | |
350 | invmom, thetax, thetay, chi2) | |
351 | ); | |
352 | ||
353 | event->AddMuonTrack(&mt); | |
354 | } | |
355 | ||
356 | #ifndef LOG_NO_DEBUG | |
357 | AliDebug(4, str); | |
358 | #endif // LOG_NO_DEBUG | |
359 | ||
360 | return 0; | |
361 | } | |
362 |