.so cleanup: removed from gSystem->Load()
[u/mrichter/AliRoot.git] / HLT / MUON / macros / MakeTrackTable.C
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] .x 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 #include "TMath.h"
33
34 #include "AliCDBManager.h"
35 #include "AliLoader.h"
36 #include "AliRunLoader.h"
37 #include "AliStack.h"
38 #include "AliRun.h"
39 #include "AliMUON.h"
40 #include "AliMUONMCDataInterface.h"
41 #include "AliMUONHit.h"
42 #include "AliHLTMUONEvent.h"
43 #include "AliHLTMUONRootifierComponent.h"
44 #include "AliHLTMUONMansoTrack.h"
45
46 #include <cstdlib>
47 #include <vector>
48 #include <iostream>
49 using std::cout;
50 using std::cerr;
51 using std::endl;
52
53
54 void MakeTrackTable(
55                 Int_t firstEvent = 0,
56                 Int_t lastEvent = -1,
57                 const char* dHLToutputfile = "output.root",
58                 Float_t maxSigma = 4., // 4 standard deviations
59                 Float_t sigmaX = 0.1,  // 1 mm resolution
60                 Float_t sigmaY = 0.01, // 100 micron resolution
61                 Float_t sigmaZ = 0.02,  // 200 microns resolution
62                 Float_t sigmaXtrg = 0.5,  // 5 mm resolution
63                 Float_t sigmaYtrg = 0.5,  // 5 mm resolution
64                 Float_t sigmaZtrg = 0.02  // 2 microns resolution
65         )
66 {
67         // Setup the CDB default storage and run number if nothing was set.
68         AliCDBManager* cdbManager = AliCDBManager::Instance();
69         if (cdbManager == NULL)
70         {
71                 cerr << "ERROR: Global CDB manager object does not exist." << endl;
72                 return;
73         }
74         if (cdbManager->GetDefaultStorage() == NULL)
75         {
76                 cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
77         }
78         if (cdbManager->GetRun() == -1)
79         {
80                 cdbManager->SetRun(0);
81         }
82         
83         gSystem->Load("libAliHLTMUON");
84         // Must pree load libAliHLTMUON.so before loading this macro and running it in compiled mode.
85
86         TString fieldnames = "event:isprimary:hastrack:cantrigger:pdgcode:sign:px:py:pz";
87         for (Int_t i = 0; i <= 14; i++)
88         {
89                 fieldnames += ":x";
90                 fieldnames += i;
91                 fieldnames += ":y";
92                 fieldnames += i;
93                 fieldnames += ":z";
94                 fieldnames += i;
95         }
96         fieldnames += ":dHLTsign:dHLTpx:dHLTpy:dHLTpz:dHLTfit";
97         for (Int_t i = 1; i <= 14; i++)
98         {
99                 fieldnames += ":dHLTx";
100                 fieldnames += i;
101                 fieldnames += ":dHLTy";
102                 fieldnames += i;
103                 fieldnames += ":dHLTz";
104                 fieldnames += i;
105         }
106         TNtuple* table = new TNtuple("tracktable", "Table of tracks.", fieldnames);
107         
108         TFile dHLTfile(dHLToutputfile, "READ");
109
110         AliMUONMCDataInterface data;
111         AliStack* stack;
112         
113         Int_t nEvents = data.NumberOfEvents();
114         if (lastEvent < 0) lastEvent = nEvents-1;
115         for (Int_t event = firstEvent; event <= lastEvent; event++)
116         {
117                 cout << "Processing event: " << event;
118                 
119                 char buf[1024];
120                 char* str = &buf[0];
121                 sprintf(str, "AliHLTMUONEvent;%d", event+1);
122                 AliHLTMUONEvent* dHLTevent = dynamic_cast<AliHLTMUONEvent*>( dHLTfile.Get(str) );
123                 if (dHLTevent == NULL)
124                 {
125                         cout << endl;
126                         cerr << "ERROR: Could not find " << str << " in " << dHLToutputfile << endl;
127                         continue;
128                 }
129                 
130                 data.GetEvent(event);
131                 
132                 stack = data.Stack(event);
133                 if (stack == NULL)
134                 {
135                         cout << endl;
136                         cerr << "ERROR: Stack is NULL" << endl;
137                         continue;
138                 }
139                 
140                 Int_t trackcount = data.NumberOfParticles();
141                 cout << "\ttrackcount = " << trackcount << endl;
142                 
143                 TArrayF isPrimary(trackcount);
144                 TArrayF hasTrack(trackcount);
145                 TArrayI hitsInTrigger(trackcount);
146                 TArrayF pdgcode(trackcount);
147                 TArrayF sign(trackcount);
148                 TArrayF px(trackcount);
149                 TArrayF py(trackcount);
150                 TArrayF pz(trackcount);
151                 TArrayF vx(trackcount);
152                 TArrayF vy(trackcount);
153                 TArrayF vz(trackcount);
154                 TArrayF dHLTsign(trackcount);
155                 TArrayF dHLTpx(trackcount);
156                 TArrayF dHLTpy(trackcount);
157                 TArrayF dHLTpz(trackcount);
158                 TArrayF dHLTfit(trackcount);
159                 
160                 //cout << "Fill particle data" << endl;
161                 for (Int_t i = 0; i < trackcount; i++)
162                 {
163                         TParticle* p = data.Particle(i);
164                         //isPrimary[i] = (i < stack->GetNprimary()) ? 1 : 0;
165                         isPrimary[i] = stack->IsPhysicalPrimary(i) ? 1 : 0;
166                         hasTrack[i] = 0;
167                         hitsInTrigger[i] = 0;
168                         pdgcode[i] = p->GetPdgCode();
169                         if (p->GetPDG() != NULL)
170                         {
171                                 Double_t charge = p->GetPDG()->Charge();
172                                 if (charge > 0)
173                                         sign[i] = 1;
174                                 else if (charge < 0)
175                                         sign[i] = -1;
176                                 else
177                                         sign[i] = 0;
178                         }
179                         else
180                         {
181                                 sign[i] = 0;
182                                 cerr << "ERROR: Could not get the PDG information." << endl;
183                         }
184                         px[i] = p->Px();
185                         py[i] = p->Py();
186                         pz[i] = p->Pz();
187                         vx[i] = p->Vx();
188                         vy[i] = p->Vy();
189                         vz[i] = p->Vz();
190                         dHLTsign[i] = 0;
191                         dHLTpx[i] = 0;
192                         dHLTpy[i] = 0;
193                         dHLTpz[i] = 0;
194                         dHLTfit[i] = -1;
195                 }
196                 
197                 TArrayF hitX[14];
198                 TArrayF hitY[14];
199                 TArrayF hitZ[14];
200                 TArrayF dHLThitX[14];
201                 TArrayF dHLThitY[14];
202                 TArrayF dHLThitZ[14];
203                 for (Int_t i = 0; i < 14; i++)
204                 {
205                         hitX[i].Set(trackcount);
206                         hitY[i].Set(trackcount);
207                         hitZ[i].Set(trackcount);
208                         dHLThitX[i].Set(trackcount);
209                         dHLThitY[i].Set(trackcount);
210                         dHLThitZ[i].Set(trackcount);
211                         for (Int_t j = 0; j < trackcount; j++)
212                         {
213                                 hitX[i][j] = 0;
214                                 hitY[i][j] = 0;
215                                 hitZ[i][j] = 0;
216                                 dHLThitX[i][j] = 0;
217                                 dHLThitY[i][j] = 0;
218                                 dHLThitZ[i][j] = 0;
219                         }
220                 }
221
222                 //cout << "Fill hits" << endl;
223                 for (Int_t i = 0; i < data.NumberOfTracks(); i++)
224                 for (Int_t j = 0; j < data.NumberOfHits(i); j++)
225                 {
226                         AliMUONHit* hit = data.Hit(i, j);
227                         Int_t chamber = hit->Chamber() - 1;
228                         if (chamber < 0 || chamber >= 14)
229                         {
230                                 cerr << "ERROR: Chamber number " << chamber << " is out of range." << endl;
231                                 continue;
232                         }
233                         
234                         //cout << "hit->Track() = " << hit->Track() << endl;
235                         if (hit->Track() < 0 || hit->Track() >= trackcount)
236                         {
237                                 cerr << "ERROR: Track number " << hit->Track() << " is out of range. [0.."
238                                         << trackcount << ")" << endl;
239                                 continue;
240                         }
241                         Int_t particleindex = hit->Track();
242
243                         hitX[chamber][particleindex] = hit->X();
244                         hitY[chamber][particleindex] = hit->Y();
245                         hitZ[chamber][particleindex] = hit->Z();
246                                         
247                         if (chamber == 0)
248                         {
249                                 hasTrack[particleindex] = 1;
250                         }
251                         if (chamber >= 10)
252                         {
253                                 hitsInTrigger[particleindex]++;
254                         }
255                 } // hits loop
256                 
257                 // Ok now go through the dHLT information for this event and correlate
258                 // the dHLT tracks to the kine data.
259                 // This is done as follows: take a dHLT track. Then for each kine track calculate
260                 // the distance between each hit of the kine track and corresponding dHLT track.
261                 // Then sum up the distances and divide by the number of hits matched.
262                 // This gives us a fit quality /Chi2 type parameter. The kine track that
263                 // has the smallest sum of differences is then chosen as the kine track
264                 // correlating to the current dHLT track being processed.
265                 for (Int_t i = 0; i < dHLTevent->Array().GetEntriesFast(); i++)
266                 {
267                         if (dHLTevent->Array().At(i)->IsA() != AliHLTMUONMansoTrack::Class())
268                                 continue;
269                         AliHLTMUONMansoTrack* mtrack = static_cast<AliHLTMUONMansoTrack*>(dHLTevent->Array().At(i));
270                         
271                         // Find the best fitting kine track.
272                         Double_t bestFitQuality = 1e30;
273                         Int_t bestFitTrack = -1;
274                         for (Int_t kinetrack = 0; kinetrack < trackcount; kinetrack++)
275                         {
276                                 Double_t sumOfDiffs = 0;
277                                 Double_t hitsMatched = 0;
278                                 for (Int_t j = 7; j <= 10; j++)
279                                 {
280                                         const AliHLTMUONRecHit* hit = mtrack->Hit(j);
281                                         if (hit == NULL) continue;
282                                         if (hitX[j-1][kinetrack] == 0 && hitY[j-1][kinetrack] == 0 && hitZ[j-1][kinetrack] == 0)
283                                                 continue;
284                                         
285                                         TVector3 hV(hitX[j-1][kinetrack], hitY[j-1][kinetrack], hitZ[j-1][kinetrack]);
286                                         TVector3 diff = hV - hit->Coordinate();
287                                         Double_t diffX = diff.X() / sigmaX;
288                                         Double_t diffY = diff.Y() / sigmaY;
289                                         Double_t diffZ = diff.Z() / sigmaZ;
290                                         if (diffX > maxSigma) continue;
291                                         if (diffY > maxSigma) continue;
292                                         if (diffZ > maxSigma) continue;
293                                         sumOfDiffs += diffX*diffX + diffY*diffY + diffZ*diffZ;
294                                         hitsMatched++;
295                                 }
296                                 if (mtrack->TriggerRecord() != NULL)
297                                 {
298                                         for (Int_t j = 11; j <= 14; j++)
299                                         {
300                                                 const TVector3& hit = mtrack->TriggerRecord()->Hit(j);
301                                                 if (hit == TVector3(0,0,0)) continue;
302                                                 if (hitX[j-1][kinetrack] == 0 && hitY[j-1][kinetrack] == 0 && hitZ[j-1][kinetrack] == 0)
303                                                         continue;
304                                                 
305                                                 TVector3 hV(hitX[j-1][kinetrack], hitY[j-1][kinetrack], hitZ[j-1][kinetrack]);
306                                                 TVector3 diff = hV - hit;
307                                                 Double_t diffX = diff.X() / sigmaXtrg;
308                                                 Double_t diffY = diff.Y() / sigmaYtrg;
309                                                 Double_t diffZ = diff.Z() / sigmaZtrg;
310                                                 if (diffX > maxSigma) continue;
311                                                 if (diffY > maxSigma) continue;
312                                                 if (diffZ > maxSigma) continue;
313                                                 sumOfDiffs += diffX*diffX + diffY*diffY + diffZ*diffZ;
314                                                 hitsMatched++;
315                                         }
316                                 }
317                                 
318                                 // Now check the fit quality.
319                                 if (hitsMatched <= 0) continue;
320                                 Double_t fitQuality = TMath::Sqrt(sumOfDiffs) / hitsMatched;
321                                 if (fitQuality < bestFitQuality)
322                                 {
323                                         bestFitQuality = fitQuality;
324                                         bestFitTrack = kinetrack;
325                                 }
326                         }
327                         
328                         if (bestFitTrack < 0)
329                         {
330                                 // No track was matched so we need to add a fake track record.
331                                 Float_t args[12+14*3+5+14*3];
332                                 args[0] = event;
333                                 args[1] = -1;
334                                 args[2] = 0;
335                                 args[3] = 0;
336                                 args[4] = 0;
337                                 args[5] = 0;
338                                 args[6] = 0;
339                                 args[7] = 0;
340                                 args[8] = 0;
341                                 args[9] = 0;
342                                 args[10] = 0;
343                                 args[11] = 0;
344                                 for (Int_t j = 0; j < 14; j++)
345                                 {                       
346                                         args[12+j*3+0] = 0;
347                                         args[12+j*3+1] = 0;
348                                         args[12+j*3+2] = 0;
349                                 }
350                                 args[12+14*3+0] = mtrack->Sign();
351                                 args[12+14*3+1] = mtrack->Px();
352                                 args[12+14*3+2] = mtrack->Py();
353                                 args[12+14*3+3] = mtrack->Pz();
354                                 args[12+14*3+4] = -1;
355                                 for (Int_t j = 0; j < 14; j++)
356                                 {                       
357                                         args[12+14*3+5+j*3+0] = 0;
358                                         args[12+14*3+5+j*3+1] = 0;
359                                         args[12+14*3+5+j*3+2] = 0;
360                                 }
361                                 for (Int_t j = 6; j < 10; j++)
362                                 {
363                                         const AliHLTMUONRecHit* hit = mtrack->Hit(j+1);
364                                         if (hit == NULL) continue;
365                                         args[12+14*3+5+j*3+0] = hit->X();
366                                         args[12+14*3+5+j*3+1] = hit->Y();
367                                         args[12+14*3+5+j*3+2] = hit->Z();
368                                 }
369                                 if (mtrack->TriggerRecord() != NULL)
370                                 {
371                                         for (Int_t j = 10; j < 14; j++)
372                                         {
373                                                 const TVector3& hit = mtrack->TriggerRecord()->Hit(j+1);
374                                                 args[12+14*3+5+j*3+0] = hit.X();
375                                                 args[12+14*3+5+j*3+1] = hit.Y();
376                                                 args[12+14*3+5+j*3+2] = hit.Z();
377                                         }
378                                 }
379                                 table->Fill(args);
380                         }
381                         else
382                         {
383                                 // Fill the details about the dHLT track to the best fitting track info.
384                                 dHLTsign[bestFitTrack] = mtrack->Sign();
385                                 dHLTpx[bestFitTrack] = mtrack->Px();
386                                 dHLTpy[bestFitTrack] = mtrack->Py();
387                                 dHLTpz[bestFitTrack] = mtrack->Pz();
388                                 dHLTfit[bestFitTrack] = bestFitQuality;
389                                 for (Int_t j = 7; j <= 10; j++)
390                                 {
391                                         const AliHLTMUONRecHit* hit = mtrack->Hit(j);
392                                         if (hit == NULL) continue;
393                                         dHLThitX[j-1][bestFitTrack] = hit->X();
394                                         dHLThitY[j-1][bestFitTrack] = hit->Y();
395                                         dHLThitZ[j-1][bestFitTrack] = hit->Z();
396                                 }
397                                 if (mtrack->TriggerRecord() != NULL)
398                                 {
399                                         for (Int_t j = 11; j <= 14; j++)
400                                         {
401                                                 const TVector3& hit = mtrack->TriggerRecord()->Hit(j);
402                                                 if (hit == TVector3(0,0,0)) continue;
403                                                 dHLThitX[j-1][bestFitTrack] = hit.X();
404                                                 dHLThitY[j-1][bestFitTrack] = hit.Y();
405                                                 dHLThitZ[j-1][bestFitTrack] = hit.Z();
406                                         }
407                                 }
408                         }
409                 }
410                 
411                 //cout << "Fill table" << endl;
412                 for (Int_t i = 0; i < trackcount; i++)
413                 {
414                         Float_t args[12+14*3+5+14*3];
415                         args[0] = event;
416                         args[1] = isPrimary[i];
417                         args[2] = hasTrack[i];
418                         args[3] = (hitsInTrigger[i] > 2) ? 1 : 0;
419                         args[4] = pdgcode[i];
420                         args[5] = sign[i];
421                         args[6] = px[i];
422                         args[7] = py[i];
423                         args[8] = pz[i];
424                         args[9] = vx[i];
425                         args[10] = vy[i];
426                         args[11] = vz[i];
427                         for (Int_t j = 0; j < 14; j++)
428                         {                       
429                                 args[12+j*3+0] = hitX[j][i];
430                                 args[12+j*3+1] = hitY[j][i];
431                                 args[12+j*3+2] = hitZ[j][i];
432                         }
433                         args[12+14*3+0] = dHLTsign[i];
434                         args[12+14*3+1] = dHLTpx[i];
435                         args[12+14*3+2] = dHLTpy[i];
436                         args[12+14*3+3] = dHLTpz[i];
437                         args[12+14*3+4] = dHLTfit[i];
438                         for (Int_t j = 0; j < 14; j++)
439                         {                       
440                                 args[12+14*3+5+j*3+0] = dHLThitX[j][i];
441                                 args[12+14*3+5+j*3+1] = dHLThitY[j][i];
442                                 args[12+14*3+5+j*3+2] = dHLThitZ[j][i];
443                         }
444                         table->Fill(args);
445                 }
446                 
447         } // event loop
448
449         TFile file("tracktable.root", "RECREATE");
450         file.cd();
451         table->Write(table->GetName(), TObject::kOverwrite);
452
453         cout << "Done." << endl;
454 }