25c6daab1b18edc314a7c421f36e08b6a059bc79
[u/mrichter/AliRoot.git] / HLT / MUON / macros / MakeTriggerTable.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
33 #include "AliCDBManager.h"
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"
41 #include "AliHLTMUONEvent.h"
42 #include "AliHLTMUONRootifierComponent.h"
43 #include "AliHLTMUONMansoTrack.h"
44
45 #include <cstdlib>
46 #include <vector>
47 #include <iostream>
48 using std::cout;
49 using std::cerr;
50 using std::endl;
51
52
53 void MakeTriggerTable(
54                 Int_t firstEvent = 0,
55                 Int_t lastEvent = -1,
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
61         )
62 {
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         {
72                 cdbManager->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
73         }
74         if (cdbManager->GetRun() == -1)
75         {
76                 cdbManager->SetRun(0);
77         }
78         
79         gSystem->Load("libAliHLTMUON");
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;
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;
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 }