]>
Commit | Line | Data |
---|---|---|
aa14645f | 1 | /************************************************************************** |
2 | * This file is property of and copyright by the ALICE HLT Project * | |
3 | * All rights reserved. * | |
4 | * * | |
5 | * Primary Authors: * | |
6 | * Artur Szostak <artursz@iafrica.com> * | |
7 | * * | |
8 | * Permission to use, copy, modify and distribute this software and its * | |
9 | * documentation strictly for non-commercial purposes is hereby granted * | |
10 | * without fee, provided that the above copyright notice appears in all * | |
11 | * copies and that both the copyright notice and this permission notice * | |
12 | * appear in the supporting documentation. The authors make no claims * | |
13 | * about the suitability of this software for any purpose. It is * | |
14 | * provided "as is" without express or implied warranty. * | |
15 | **************************************************************************/ | |
16 | ||
17 | // Simple macro to generate N-tuple for performance analysis. | |
18 | // This macro must be compiled and run like so from within the aliroot command | |
19 | // prompt: | |
20 | // root [0] gSystem->Load("libAliHLTMUON.so"); | |
21 | // root [0] .L MakeHitsTable.C+ | |
22 | ||
23 | #include "TVector3.h" | |
24 | #include "TSystem.h" | |
25 | #include "TROOT.h" | |
26 | #include "TParticle.h" | |
27 | #include "TArrayF.h" | |
28 | #include "TStopwatch.h" | |
29 | #include "TNtuple.h" | |
30 | #include "TFile.h" | |
31 | #include "TError.h" | |
32 | ||
33 | #include "AliLoader.h" | |
34 | #include "AliRunLoader.h" | |
35 | #include "AliStack.h" | |
36 | #include "AliRun.h" | |
37 | #include "AliMUON.h" | |
38 | #include "AliMUONMCDataInterface.h" | |
39 | #include "AliMUONHit.h" | |
40 | #include "AliHLTMUONRootifierComponent.h" | |
41 | #include "AliHLTMUONMansoTrack.h" | |
42 | ||
43 | #include <cstdlib> | |
44 | #include <vector> | |
45 | #include <iostream> | |
46 | using std::cout; | |
47 | using std::cerr; | |
48 | using std::endl; | |
49 | ||
50 | ||
51 | void MakeTriggerTable( | |
52 | Int_t firstEvent = 0, | |
53 | Int_t lastEvent = -1, | |
54 | const char* L0outputfile = "output_0x00000000.root" | |
55 | ) | |
56 | { | |
57 | gSystem->Load("libAliHLTMUON.so"); | |
58 | // Must pree load libAliHLTMUON.so before loading this macro and running it in compiled mode. | |
59 | ||
60 | TString fieldnames = "event:isprimary:pdgcode:sign:px:py:pz"; | |
61 | for (Int_t i = 11; i <= 14; i++) | |
62 | { | |
63 | fieldnames += ":x"; | |
64 | fieldnames += i; | |
65 | fieldnames += ":y"; | |
66 | fieldnames += i; | |
67 | fieldnames += ":z"; | |
68 | fieldnames += i; | |
69 | } | |
70 | fieldnames += ":L0sign:L0px:L0py:L0pz:L0fit"; | |
71 | for (Int_t i = 11; i <= 14; i++) | |
72 | { | |
73 | fieldnames += ":L0x"; | |
74 | fieldnames += i; | |
75 | fieldnames += ":L0y"; | |
76 | fieldnames += i; | |
77 | fieldnames += ":L0z"; | |
78 | fieldnames += i; | |
79 | } | |
80 | TNtuple* table = new TNtuple("triggertable", "Table of triggers.", fieldnames); | |
81 | ||
82 | TFile L0file(L0outputfile, "READ"); | |
83 | ||
84 | AliMUONMCDataInterface data; | |
85 | AliStack* stack; | |
86 | ||
87 | Int_t nEvents = data.NumberOfEvents(); | |
88 | if (lastEvent < 0) lastEvent = nEvents-1; | |
89 | for (Int_t event = firstEvent; event <= lastEvent; event++) | |
90 | { | |
91 | cout << "Processing event: " << event; | |
92 | ||
93 | char buf[1024]; | |
94 | char* str = &buf[0]; | |
95 | sprintf(str, "AliHLTMUONEvent;%d", event+1); | |
96 | AliHLTMUONEvent* L0event = dynamic_cast<AliHLTMUONEvent*>( L0file.Get(str) ); | |
97 | if (L0event == NULL) | |
98 | { | |
99 | cout << endl; | |
100 | cerr << "ERROR: Could not find " << str << " in " << L0outputfile << endl; | |
101 | continue; | |
102 | } | |
103 | ||
104 | data.GetEvent(event); | |
105 | ||
106 | stack = data.Stack(event); | |
107 | if (stack == NULL) | |
108 | { | |
109 | cout << endl; | |
110 | cerr << "ERROR: Stack is NULL" << endl; | |
111 | continue; | |
112 | } | |
113 | ||
114 | Int_t trackcount = data.NumberOfParticles(); | |
115 | cout << "\ttrackcount = " << trackcount << endl; | |
116 | ||
117 | TArrayF isPrimary(trackcount); | |
118 | TArrayF pdgcode(trackcount); | |
119 | TArrayF sign(trackcount); | |
120 | TArrayF px(trackcount); | |
121 | TArrayF py(trackcount); | |
122 | TArrayF pz(trackcount); | |
123 | TArrayF L0sign(trackcount); | |
124 | TArrayF L0px(trackcount); | |
125 | TArrayF L0py(trackcount); | |
126 | TArrayF L0pz(trackcount); | |
127 | TArrayF L0fit(trackcount); | |
128 | ||
129 | //cout << "Fill particle data" << endl; | |
130 | for (Int_t i = 0; i < trackcount; i++) | |
131 | { | |
132 | TParticle* p = data.Particle(i); | |
133 | isPrimary[i] = stack->IsPhysicalPrimary(i) ? 1 : 0; | |
134 | pdgcode[i] = p->GetPdgCode(); | |
135 | if (p->GetPDG() != NULL) | |
136 | { | |
137 | Double_t charge = p->GetPDG()->Charge(); | |
138 | if (charge > 0) | |
139 | sign[i] = 1; | |
140 | else if (charge < 0) | |
141 | sign[i] = -1; | |
142 | else | |
143 | sign[i] = 0; | |
144 | } | |
145 | else | |
146 | { | |
147 | sign[i] = 0; | |
148 | cerr << "ERROR: Could not get the PDG information." << endl; | |
149 | } | |
150 | px[i] = p->Px(); | |
151 | py[i] = p->Py(); | |
152 | pz[i] = p->Pz(); | |
153 | L0sign[i] = 0; | |
154 | L0px[i] = 0; | |
155 | L0py[i] = 0; | |
156 | L0pz[i] = 0; | |
157 | L0fit[i] = -1; | |
158 | } | |
159 | ||
160 | TArrayF hitX[4]; | |
161 | TArrayF hitY[4]; | |
162 | TArrayF hitZ[4]; | |
163 | TArrayF L0hitX[4]; | |
164 | TArrayF L0hitY[4]; | |
165 | TArrayF L0hitZ[4]; | |
166 | for (Int_t i = 0; i < 4; i++) | |
167 | { | |
168 | hitX[i].Set(trackcount); | |
169 | hitY[i].Set(trackcount); | |
170 | hitZ[i].Set(trackcount); | |
171 | L0hitX[i].Set(trackcount); | |
172 | L0hitY[i].Set(trackcount); | |
173 | L0hitZ[i].Set(trackcount); | |
174 | for (Int_t j = 0; j < trackcount; j++) | |
175 | { | |
176 | hitX[i][j] = 0; | |
177 | hitY[i][j] = 0; | |
178 | hitZ[i][j] = 0; | |
179 | L0hitX[i][j] = 0; | |
180 | L0hitY[i][j] = 0; | |
181 | L0hitZ[i][j] = 0; | |
182 | } | |
183 | } | |
184 | ||
185 | //cout << "Fill hits" << endl; | |
186 | for (Int_t i = 0; i < data.NumberOfTracks(); i++) | |
187 | for (Int_t j = 0; j < data.NumberOfHits(i); j++) | |
188 | { | |
189 | AliMUONHit* hit = data.Hit(i, j); | |
190 | Int_t chamber = hit->Chamber() - 1; | |
191 | if (chamber < 0 || chamber >= 14) | |
192 | { | |
193 | cerr << "ERROR: Chamber number " << chamber << " is out of range." << endl; | |
194 | continue; | |
195 | } | |
196 | ||
197 | //cout << "hit->Track() = " << hit->Track() << endl; | |
198 | if (hit->Track() < 0 || hit->Track() >= trackcount) | |
199 | { | |
200 | cerr << "ERROR: Track number " << hit->Track() << " is out of range. [0.." | |
201 | << trackcount << ")" << endl; | |
202 | continue; | |
203 | } | |
204 | Int_t particleindex = hit->Track(); | |
205 | ||
206 | if (chamber >= 10) | |
207 | { | |
208 | hitX[chamber-10][particleindex] = hit->X(); | |
209 | hitY[chamber-10][particleindex] = hit->Y(); | |
210 | hitZ[chamber-10][particleindex] = hit->Z(); | |
211 | } | |
212 | } // hits loop | |
213 | ||
214 | // Ok now go through the L0 information for this event and correlate | |
215 | // the L0 trigger records to the kine data. | |
216 | // This is done as follows: take a trigger record. Then for each kine track | |
217 | // calculate the distance between each hit of the kine track and corresponding | |
218 | // L0 trigger track fragment hit. | |
219 | // Then sum up the distances and divide by the number of hits matched. | |
220 | // This gives us a fit quality /Chi2 type parameter. The kine track that | |
221 | // has the smallest sum of differences is then chosen as the kine track | |
222 | // correlating to the current L0 trigger record being processed. | |
223 | for (Int_t i = 0; i < L0event->Array().GetEntriesFast(); i++) | |
224 | { | |
225 | if (L0event->Array().At(i)->IsA() != AliHLTMUONTriggerRecord::Class()) | |
226 | continue; | |
227 | AliHLTMUONTriggerRecord* trigrec = static_cast<AliHLTMUONTriggerRecord*>(L0event->Array().At(i)); | |
228 | ||
229 | // Find the best fitting kine track. | |
230 | Double_t bestFitQuality = 1e30; | |
231 | Int_t bestFitTrack = -1; | |
232 | for (Int_t kinetrack = 0; kinetrack < trackcount; kinetrack++) | |
233 | { | |
234 | Double_t sumOfDiffs = 0; | |
235 | Double_t hitsMatched = 0; | |
236 | for (Int_t j = 11; j <= 14; j++) | |
237 | { | |
238 | const TVector3& hit = trigrec->Hit(j); | |
239 | if (hit == TVector3(0,0,0)) continue; | |
240 | TVector3 hV(hitX[j-11][kinetrack], hitY[j-11][kinetrack], hitZ[j-11][kinetrack]); | |
241 | if (hV == TVector3(0,0,0)) continue; | |
242 | TVector3 diff = hV - hit; | |
243 | sumOfDiffs += diff.Mag2(); // Use the square of diff. More Chi^2 like properties. | |
244 | hitsMatched++; | |
245 | } | |
246 | ||
247 | // Now check the fit quality. | |
248 | if (hitsMatched <= 0) continue; | |
249 | Double_t fitQuality = sumOfDiffs / hitsMatched; | |
250 | if (fitQuality < bestFitQuality) | |
251 | { | |
252 | bestFitQuality = fitQuality; | |
253 | bestFitTrack = kinetrack; | |
254 | } | |
255 | } | |
256 | ||
257 | if (bestFitTrack < 0) | |
258 | { | |
259 | // No track was matched so we need to add a fake track record. | |
260 | Float_t args[7+4*3+5+4*3]; | |
261 | args[0] = event; | |
262 | args[1] = -1; | |
263 | args[2] = 0; | |
264 | args[3] = 0; | |
265 | args[4] = 0; | |
266 | args[5] = 0; | |
267 | args[6] = 0; | |
268 | for (Int_t j = 0; j < 4; j++) | |
269 | { | |
270 | args[7+j*3+0] = 0; | |
271 | args[7+j*3+1] = 0; | |
272 | args[7+j*3+2] = 0; | |
273 | } | |
274 | args[7+4*3+0] = trigrec->Sign(); | |
275 | args[7+4*3+1] = trigrec->Px(); | |
276 | args[7+4*3+2] = trigrec->Py(); | |
277 | args[7+4*3+3] = trigrec->Pz(); | |
278 | args[7+4*3+4] = -1; | |
279 | for (Int_t j = 0; j < 4; j++) | |
280 | { | |
281 | const TVector3& hit = trigrec->Hit(j+11); | |
282 | args[7+4*3+5+j*3+0] = hit.X(); | |
283 | args[7+4*3+5+j*3+1] = hit.Y(); | |
284 | args[7+4*3+5+j*3+2] = hit.Z(); | |
285 | } | |
286 | table->Fill(args); | |
287 | } | |
288 | else | |
289 | { | |
290 | // Fill the details about the L0 track to the best fitting track info. | |
291 | L0sign[bestFitTrack] = trigrec->Sign(); | |
292 | L0px[bestFitTrack] = trigrec->Px(); | |
293 | L0py[bestFitTrack] = trigrec->Py(); | |
294 | L0pz[bestFitTrack] = trigrec->Pz(); | |
295 | L0fit[bestFitTrack] = bestFitQuality; | |
296 | for (Int_t j = 11; j <= 14; j++) | |
297 | { | |
298 | const TVector3& hit = trigrec->Hit(j); | |
299 | if (hit == TVector3(0,0,0)) continue; | |
300 | L0hitX[j-11][bestFitTrack] = hit.X(); | |
301 | L0hitY[j-11][bestFitTrack] = hit.Y(); | |
302 | L0hitZ[j-11][bestFitTrack] = hit.Z(); | |
303 | } | |
304 | } | |
305 | } | |
306 | ||
307 | //cout << "Fill table" << endl; | |
308 | for (Int_t i = 0; i < trackcount; i++) | |
309 | { | |
310 | Float_t args[7+4*3+5+4*3]; | |
311 | args[0] = event; | |
312 | args[1] = isPrimary[i]; | |
313 | args[2] = pdgcode[i]; | |
314 | args[3] = sign[i]; | |
315 | args[4] = px[i]; | |
316 | args[5] = py[i]; | |
317 | args[6] = pz[i]; | |
318 | for (Int_t j = 0; j < 4; j++) | |
319 | { | |
320 | args[7+j*3+0] = hitX[j][i]; | |
321 | args[7+j*3+1] = hitY[j][i]; | |
322 | args[7+j*3+2] = hitZ[j][i]; | |
323 | } | |
324 | args[7+4*3+0] = L0sign[i]; | |
325 | args[7+4*3+1] = L0px[i]; | |
326 | args[7+4*3+2] = L0py[i]; | |
327 | args[7+4*3+3] = L0pz[i]; | |
328 | args[7+4*3+4] = L0fit[i]; | |
329 | for (Int_t j = 0; j < 4; j++) | |
330 | { | |
331 | args[7+4*3+5+j*3+0] = L0hitX[j][i]; | |
332 | args[7+4*3+5+j*3+1] = L0hitY[j][i]; | |
333 | args[7+4*3+5+j*3+2] = L0hitZ[j][i]; | |
334 | } | |
335 | table->Fill(args); | |
336 | } | |
337 | ||
338 | } // event loop | |
339 | ||
340 | TFile file("triggertable.root", "RECREATE"); | |
341 | file.cd(); | |
342 | table->Write(table->GetName(), TObject::kOverwrite); | |
343 | ||
344 | cout << "Done." << endl; | |
345 | } |