]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/MUON/macros/MakeTriggerTable.C
Adding a protection in the case of DAQ FXS failure
[u/mrichter/AliRoot.git] / HLT / MUON / macros / MakeTriggerTable.C
CommitLineData
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>
46using std::cout;
47using std::cerr;
48using std::endl;
49
50
51void 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}