.so cleanup: more gSystem->Load()
[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:
b0635849 20// root [0] gSystem->Load("libAliHLTMUON");
d8cf7391 21// root [0] .x MakeHitsTable.C+
aa14645f 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
846e6f41 33#include "AliCDBManager.h"
aa14645f 34#include "AliLoader.h"
35#include "AliRunLoader.h"
36#include "AliStack.h"
37#include "AliRun.h"
38#include "AliMUON.h"
39#include "AliMUONMCDataInterface.h"
40#include "AliMUONHit.h"
3b376c57 41#include "AliHLTMUONEvent.h"
aa14645f 42#include "AliHLTMUONRootifierComponent.h"
43#include "AliHLTMUONMansoTrack.h"
44
45#include <cstdlib>
46#include <vector>
47#include <iostream>
48using std::cout;
49using std::cerr;
50using std::endl;
51
52
53void MakeTriggerTable(
54 Int_t firstEvent = 0,
55 Int_t lastEvent = -1,
d8cf7391 56 const char* L0outputfile = "output.root",
af5746f2 57 Float_t maxSigma = 4., // 4 standard deviations
58 Float_t sigmaXtrg = 0.5, // 5 mm resolution
59 Float_t sigmaYtrg = 0.5, // 5 mm resolution
60 Float_t sigmaZtrg = 0.02 // 2 microns resolution
aa14645f 61 )
62{
846e6f41 63 // Setup the CDB default storage and run number if nothing was set.
64 AliCDBManager* cdbManager = AliCDBManager::Instance();
65 if (cdbManager == NULL)
66 {
67 cerr << "ERROR: Global CDB manager object does not exist." << endl;
68 return;
69 }
70 if (cdbManager->GetDefaultStorage() == NULL)
71 {
162637e4 72 cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
846e6f41 73 }
74 if (cdbManager->GetRun() == -1)
75 {
76 cdbManager->SetRun(0);
77 }
78
4070f709 79 gSystem->Load("libAliHLTMUON");
aa14645f 80 // Must pree load libAliHLTMUON.so before loading this macro and running it in compiled mode.
81
82 TString fieldnames = "event:isprimary:pdgcode:sign:px:py:pz";
83 for (Int_t i = 11; i <= 14; i++)
84 {
85 fieldnames += ":x";
86 fieldnames += i;
87 fieldnames += ":y";
88 fieldnames += i;
89 fieldnames += ":z";
90 fieldnames += i;
91 }
92 fieldnames += ":L0sign:L0px:L0py:L0pz:L0fit";
93 for (Int_t i = 11; i <= 14; i++)
94 {
95 fieldnames += ":L0x";
96 fieldnames += i;
97 fieldnames += ":L0y";
98 fieldnames += i;
99 fieldnames += ":L0z";
100 fieldnames += i;
101 }
102 TNtuple* table = new TNtuple("triggertable", "Table of triggers.", fieldnames);
103
104 TFile L0file(L0outputfile, "READ");
105
106 AliMUONMCDataInterface data;
107 AliStack* stack;
108
109 Int_t nEvents = data.NumberOfEvents();
110 if (lastEvent < 0) lastEvent = nEvents-1;
111 for (Int_t event = firstEvent; event <= lastEvent; event++)
112 {
113 cout << "Processing event: " << event;
114
115 char buf[1024];
116 char* str = &buf[0];
117 sprintf(str, "AliHLTMUONEvent;%d", event+1);
118 AliHLTMUONEvent* L0event = dynamic_cast<AliHLTMUONEvent*>( L0file.Get(str) );
119 if (L0event == NULL)
120 {
121 cout << endl;
122 cerr << "ERROR: Could not find " << str << " in " << L0outputfile << endl;
123 continue;
124 }
125
126 data.GetEvent(event);
127
128 stack = data.Stack(event);
129 if (stack == NULL)
130 {
131 cout << endl;
132 cerr << "ERROR: Stack is NULL" << endl;
133 continue;
134 }
135
136 Int_t trackcount = data.NumberOfParticles();
137 cout << "\ttrackcount = " << trackcount << endl;
138
139 TArrayF isPrimary(trackcount);
140 TArrayF pdgcode(trackcount);
141 TArrayF sign(trackcount);
142 TArrayF px(trackcount);
143 TArrayF py(trackcount);
144 TArrayF pz(trackcount);
145 TArrayF L0sign(trackcount);
146 TArrayF L0px(trackcount);
147 TArrayF L0py(trackcount);
148 TArrayF L0pz(trackcount);
149 TArrayF L0fit(trackcount);
150
151 //cout << "Fill particle data" << endl;
152 for (Int_t i = 0; i < trackcount; i++)
153 {
154 TParticle* p = data.Particle(i);
155 isPrimary[i] = stack->IsPhysicalPrimary(i) ? 1 : 0;
156 pdgcode[i] = p->GetPdgCode();
157 if (p->GetPDG() != NULL)
158 {
159 Double_t charge = p->GetPDG()->Charge();
160 if (charge > 0)
161 sign[i] = 1;
162 else if (charge < 0)
163 sign[i] = -1;
164 else
165 sign[i] = 0;
166 }
167 else
168 {
169 sign[i] = 0;
170 cerr << "ERROR: Could not get the PDG information." << endl;
171 }
172 px[i] = p->Px();
173 py[i] = p->Py();
174 pz[i] = p->Pz();
175 L0sign[i] = 0;
176 L0px[i] = 0;
177 L0py[i] = 0;
178 L0pz[i] = 0;
179 L0fit[i] = -1;
180 }
181
182 TArrayF hitX[4];
183 TArrayF hitY[4];
184 TArrayF hitZ[4];
185 TArrayF L0hitX[4];
186 TArrayF L0hitY[4];
187 TArrayF L0hitZ[4];
188 for (Int_t i = 0; i < 4; i++)
189 {
190 hitX[i].Set(trackcount);
191 hitY[i].Set(trackcount);
192 hitZ[i].Set(trackcount);
193 L0hitX[i].Set(trackcount);
194 L0hitY[i].Set(trackcount);
195 L0hitZ[i].Set(trackcount);
196 for (Int_t j = 0; j < trackcount; j++)
197 {
198 hitX[i][j] = 0;
199 hitY[i][j] = 0;
200 hitZ[i][j] = 0;
201 L0hitX[i][j] = 0;
202 L0hitY[i][j] = 0;
203 L0hitZ[i][j] = 0;
204 }
205 }
206
207 //cout << "Fill hits" << endl;
208 for (Int_t i = 0; i < data.NumberOfTracks(); i++)
209 for (Int_t j = 0; j < data.NumberOfHits(i); j++)
210 {
211 AliMUONHit* hit = data.Hit(i, j);
212 Int_t chamber = hit->Chamber() - 1;
213 if (chamber < 0 || chamber >= 14)
214 {
215 cerr << "ERROR: Chamber number " << chamber << " is out of range." << endl;
216 continue;
217 }
218
219 //cout << "hit->Track() = " << hit->Track() << endl;
220 if (hit->Track() < 0 || hit->Track() >= trackcount)
221 {
222 cerr << "ERROR: Track number " << hit->Track() << " is out of range. [0.."
223 << trackcount << ")" << endl;
224 continue;
225 }
226 Int_t particleindex = hit->Track();
227
228 if (chamber >= 10)
229 {
230 hitX[chamber-10][particleindex] = hit->X();
231 hitY[chamber-10][particleindex] = hit->Y();
232 hitZ[chamber-10][particleindex] = hit->Z();
233 }
234 } // hits loop
235
236 // Ok now go through the L0 information for this event and correlate
237 // the L0 trigger records to the kine data.
238 // This is done as follows: take a trigger record. Then for each kine track
239 // calculate the distance between each hit of the kine track and corresponding
240 // L0 trigger track fragment hit.
241 // Then sum up the distances and divide by the number of hits matched.
242 // This gives us a fit quality /Chi2 type parameter. The kine track that
243 // has the smallest sum of differences is then chosen as the kine track
244 // correlating to the current L0 trigger record being processed.
245 for (Int_t i = 0; i < L0event->Array().GetEntriesFast(); i++)
246 {
247 if (L0event->Array().At(i)->IsA() != AliHLTMUONTriggerRecord::Class())
248 continue;
249 AliHLTMUONTriggerRecord* trigrec = static_cast<AliHLTMUONTriggerRecord*>(L0event->Array().At(i));
250
251 // Find the best fitting kine track.
252 Double_t bestFitQuality = 1e30;
253 Int_t bestFitTrack = -1;
254 for (Int_t kinetrack = 0; kinetrack < trackcount; kinetrack++)
255 {
256 Double_t sumOfDiffs = 0;
257 Double_t hitsMatched = 0;
258 for (Int_t j = 11; j <= 14; j++)
259 {
260 const TVector3& hit = trigrec->Hit(j);
261 if (hit == TVector3(0,0,0)) continue;
262 TVector3 hV(hitX[j-11][kinetrack], hitY[j-11][kinetrack], hitZ[j-11][kinetrack]);
263 if (hV == TVector3(0,0,0)) continue;
264 TVector3 diff = hV - hit;
af5746f2 265 Double_t diffX = diff.X() / sigmaXtrg;
266 Double_t diffY = diff.Y() / sigmaYtrg;
267 Double_t diffZ = diff.Z() / sigmaZtrg;
268 if (diffX > maxSigma) continue;
269 if (diffY > maxSigma) continue;
270 if (diffZ > maxSigma) continue;
271 sumOfDiffs += diffX*diffX + diffY*diffY + diffZ*diffZ;
aa14645f 272 hitsMatched++;
273 }
274
275 // Now check the fit quality.
276 if (hitsMatched <= 0) continue;
277 Double_t fitQuality = sumOfDiffs / hitsMatched;
278 if (fitQuality < bestFitQuality)
279 {
280 bestFitQuality = fitQuality;
281 bestFitTrack = kinetrack;
282 }
283 }
284
285 if (bestFitTrack < 0)
286 {
287 // No track was matched so we need to add a fake track record.
288 Float_t args[7+4*3+5+4*3];
289 args[0] = event;
290 args[1] = -1;
291 args[2] = 0;
292 args[3] = 0;
293 args[4] = 0;
294 args[5] = 0;
295 args[6] = 0;
296 for (Int_t j = 0; j < 4; j++)
297 {
298 args[7+j*3+0] = 0;
299 args[7+j*3+1] = 0;
300 args[7+j*3+2] = 0;
301 }
302 args[7+4*3+0] = trigrec->Sign();
303 args[7+4*3+1] = trigrec->Px();
304 args[7+4*3+2] = trigrec->Py();
305 args[7+4*3+3] = trigrec->Pz();
306 args[7+4*3+4] = -1;
307 for (Int_t j = 0; j < 4; j++)
308 {
309 const TVector3& hit = trigrec->Hit(j+11);
310 args[7+4*3+5+j*3+0] = hit.X();
311 args[7+4*3+5+j*3+1] = hit.Y();
312 args[7+4*3+5+j*3+2] = hit.Z();
313 }
314 table->Fill(args);
315 }
316 else
317 {
318 // Fill the details about the L0 track to the best fitting track info.
319 L0sign[bestFitTrack] = trigrec->Sign();
320 L0px[bestFitTrack] = trigrec->Px();
321 L0py[bestFitTrack] = trigrec->Py();
322 L0pz[bestFitTrack] = trigrec->Pz();
323 L0fit[bestFitTrack] = bestFitQuality;
324 for (Int_t j = 11; j <= 14; j++)
325 {
326 const TVector3& hit = trigrec->Hit(j);
327 if (hit == TVector3(0,0,0)) continue;
328 L0hitX[j-11][bestFitTrack] = hit.X();
329 L0hitY[j-11][bestFitTrack] = hit.Y();
330 L0hitZ[j-11][bestFitTrack] = hit.Z();
331 }
332 }
333 }
334
335 //cout << "Fill table" << endl;
336 for (Int_t i = 0; i < trackcount; i++)
337 {
338 Float_t args[7+4*3+5+4*3];
339 args[0] = event;
340 args[1] = isPrimary[i];
341 args[2] = pdgcode[i];
342 args[3] = sign[i];
343 args[4] = px[i];
344 args[5] = py[i];
345 args[6] = pz[i];
346 for (Int_t j = 0; j < 4; j++)
347 {
348 args[7+j*3+0] = hitX[j][i];
349 args[7+j*3+1] = hitY[j][i];
350 args[7+j*3+2] = hitZ[j][i];
351 }
352 args[7+4*3+0] = L0sign[i];
353 args[7+4*3+1] = L0px[i];
354 args[7+4*3+2] = L0py[i];
355 args[7+4*3+3] = L0pz[i];
356 args[7+4*3+4] = L0fit[i];
357 for (Int_t j = 0; j < 4; j++)
358 {
359 args[7+4*3+5+j*3+0] = L0hitX[j][i];
360 args[7+4*3+5+j*3+1] = L0hitY[j][i];
361 args[7+4*3+5+j*3+2] = L0hitZ[j][i];
362 }
363 table->Fill(args);
364 }
365
366 } // event loop
367
368 TFile file("triggertable.root", "RECREATE");
369 file.cd();
370 table->Write(table->GetName(), TObject::kOverwrite);
371
372 cout << "Done." << endl;
373}