1 /**************************************************************************
2 * This file is property of and copyright by the ALICE HLT Project *
3 * All rights reserved. *
6 * Artur Szostak <artursz@iafrica.com> *
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 **************************************************************************/
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
20 // root [0] gSystem->Load("libAliHLTMUON");
21 // root [0] .x MakeHitsTable.C+
26 #include "TParticle.h"
28 #include "TStopwatch.h"
33 #include "AliCDBManager.h"
34 #include "AliLoader.h"
35 #include "AliRunLoader.h"
39 #include "AliMUONMCDataInterface.h"
40 #include "AliMUONHit.h"
41 #include "AliHLTMUONEvent.h"
42 #include "AliHLTMUONRootifierComponent.h"
43 #include "AliHLTMUONMansoTrack.h"
53 void MakeTriggerTable(
56 const char* L0outputfile = "output.root",
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
63 // Setup the CDB default storage and run number if nothing was set.
64 AliCDBManager* cdbManager = AliCDBManager::Instance();
65 if (cdbManager == NULL)
67 cerr << "ERROR: Global CDB manager object does not exist." << endl;
70 if (cdbManager->GetDefaultStorage() == NULL)
72 cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
74 if (cdbManager->GetRun() == -1)
76 cdbManager->SetRun(0);
79 gSystem->Load("libAliHLTMUON");
80 // Must pree load libAliHLTMUON.so before loading this macro and running it in compiled mode.
82 TString fieldnames = "event:isprimary:pdgcode:sign:px:py:pz";
83 for (Int_t i = 11; i <= 14; i++)
92 fieldnames += ":L0sign:L0px:L0py:L0pz:L0fit";
93 for (Int_t i = 11; i <= 14; i++)
102 TNtuple* table = new TNtuple("triggertable", "Table of triggers.", fieldnames);
104 TFile L0file(L0outputfile, "READ");
106 AliMUONMCDataInterface data;
109 Int_t nEvents = data.NumberOfEvents();
110 if (lastEvent < 0) lastEvent = nEvents-1;
111 for (Int_t event = firstEvent; event <= lastEvent; event++)
113 cout << "Processing event: " << event;
117 sprintf(str, "AliHLTMUONEvent;%d", event+1);
118 AliHLTMUONEvent* L0event = dynamic_cast<AliHLTMUONEvent*>( L0file.Get(str) );
122 cerr << "ERROR: Could not find " << str << " in " << L0outputfile << endl;
126 data.GetEvent(event);
128 stack = data.Stack(event);
132 cerr << "ERROR: Stack is NULL" << endl;
136 Int_t trackcount = data.NumberOfParticles();
137 cout << "\ttrackcount = " << trackcount << endl;
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);
151 //cout << "Fill particle data" << endl;
152 for (Int_t i = 0; i < trackcount; i++)
154 TParticle* p = data.Particle(i);
155 isPrimary[i] = stack->IsPhysicalPrimary(i) ? 1 : 0;
156 pdgcode[i] = p->GetPdgCode();
157 if (p->GetPDG() != NULL)
159 Double_t charge = p->GetPDG()->Charge();
170 cerr << "ERROR: Could not get the PDG information." << endl;
188 for (Int_t i = 0; i < 4; i++)
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++)
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++)
211 AliMUONHit* hit = data.Hit(i, j);
212 Int_t chamber = hit->Chamber() - 1;
213 if (chamber < 0 || chamber >= 14)
215 cerr << "ERROR: Chamber number " << chamber << " is out of range." << endl;
219 //cout << "hit->Track() = " << hit->Track() << endl;
220 if (hit->Track() < 0 || hit->Track() >= trackcount)
222 cerr << "ERROR: Track number " << hit->Track() << " is out of range. [0.."
223 << trackcount << ")" << endl;
226 Int_t particleindex = hit->Track();
230 hitX[chamber-10][particleindex] = hit->X();
231 hitY[chamber-10][particleindex] = hit->Y();
232 hitZ[chamber-10][particleindex] = hit->Z();
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++)
247 if (L0event->Array().At(i)->IsA() != AliHLTMUONTriggerRecord::Class())
249 AliHLTMUONTriggerRecord* trigrec = static_cast<AliHLTMUONTriggerRecord*>(L0event->Array().At(i));
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++)
256 Double_t sumOfDiffs = 0;
257 Double_t hitsMatched = 0;
258 for (Int_t j = 11; j <= 14; j++)
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;
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;
275 // Now check the fit quality.
276 if (hitsMatched <= 0) continue;
277 Double_t fitQuality = sumOfDiffs / hitsMatched;
278 if (fitQuality < bestFitQuality)
280 bestFitQuality = fitQuality;
281 bestFitTrack = kinetrack;
285 if (bestFitTrack < 0)
287 // No track was matched so we need to add a fake track record.
288 Float_t args[7+4*3+5+4*3];
296 for (Int_t j = 0; j < 4; j++)
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();
307 for (Int_t j = 0; j < 4; j++)
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();
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++)
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();
335 //cout << "Fill table" << endl;
336 for (Int_t i = 0; i < trackcount; i++)
338 Float_t args[7+4*3+5+4*3];
340 args[1] = isPrimary[i];
341 args[2] = pdgcode[i];
346 for (Int_t j = 0; j < 4; j++)
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];
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++)
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];
368 TFile file("triggertable.root", "RECREATE");
370 table->Write(table->GetName(), TObject::kOverwrite);
372 cout << "Done." << endl;