]>
Commit | Line | Data |
---|---|---|
2c0075de | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
15 | ||
88cb7938 | 16 | /* $Id$ */ |
2c0075de | 17 | |
c7ba256d | 18 | //Authors: Mihaela Gheata, Andrei Gheata 09/10/00 |
2c0075de | 19 | ////////////////////////////////////////////////////////////////////// |
20 | // // | |
21 | // AliMUONRecoDisplay // | |
22 | // // | |
23 | // This class subclasses AliDisplay and provides display of // | |
24 | // reconstructed tracks with following functionality : // | |
25 | // - front/top/side/3D display of MUON reconstructed tracks // | |
26 | // as polylines ; // | |
27 | // - context menu activated when the main pad is right-clicked // | |
28 | // The context menu contains following functions : // | |
29 | // * SetDrawHits() - switches on or off Geant hits ; // | |
30 | // * CutMomentum() - displays only tracks within Pmin - Pmax // | |
31 | // * ListTracks() - prints ID and momentum info. for all // | |
32 | // tracks within momentum range Pmin,Pmax ; // | |
33 | // * Highlight() - shows only one selected reco. track // | |
34 | // and its best matching Geant track; // | |
35 | // * UnHighlight() - self explaining; // | |
36 | // * RecoEfficiency() - compute reco. efficiency for all events// | |
37 | // from galice.root file; also fake track percentage; make // | |
38 | // plots for momentum precision // | |
39 | // * XYPlot() - make X-Y plots of reconstructed and // | |
40 | // generated tracks in all chambers // | |
41 | // // | |
42 | // Starting : generate and reconstruct events, then use the // | |
43 | // MUONrecodisplay.C macro // | |
44 | // // | |
45 | ////////////////////////////////////////////////////////////////////// | |
c7ba256d | 46 | |
c7ba256d | 47 | #include <TROOT.h> |
30178c30 | 48 | #include <TApplication.h> |
49 | #include <TFile.h> | |
50 | #include <TPolyLine3D.h> | |
51 | #include <TParticle.h> | |
52 | #include <TTree.h> | |
53 | #include <TH1.h> | |
54 | #include <TH2.h> | |
55 | #include <TCanvas.h> | |
56 | #include <TProfile.h> | |
88cb7938 | 57 | #include <TFolder.h> |
58 | #include <TClonesArray.h> | |
c7ba256d | 59 | #include <TSlider.h> |
c7ba256d | 60 | #include <TGeometry.h> |
88cb7938 | 61 | #include <TView.h> |
30178c30 | 62 | #include <Riostream.h> |
88cb7938 | 63 | |
30178c30 | 64 | #include "AliRun.h" |
65 | #include "AliDetector.h" | |
66 | #include "AliConfig.h" | |
67 | #include "AliHeader.h" | |
68 | #include "AliPoints.h" | |
69 | #include "AliMC.h" | |
88cb7938 | 70 | |
88cb7938 | 71 | #include "AliMUONRecoDisplay.h" |
345c4d04 | 72 | #include "AliMUONRecoEvent.h" |
73 | #include "AliMUONRecoTrack.h" | |
30178c30 | 74 | #include "AliMUONHit.h" |
8c343c7c | 75 | #include "AliLog.h" |
c7ba256d | 76 | |
77 | ClassImp(AliMUONRecoDisplay) | |
78 | ||
79 | //------------------------------------------------------------------- | |
80 | AliMUONRecoDisplay::AliMUONRecoDisplay(Int_t nevent) | |
30178c30 | 81 | : AliDisplay(750) |
c7ba256d | 82 | { |
83 | //************ Constructor of the reco. event display********** | |
84 | // get reconstructed event from file | |
85 | fFile = new TFile("tree_reco.root"); | |
86 | if (!fFile) { | |
87 | cout << "File tree_reco.root not found\n"; | |
88 | gApplication->Terminate(0); | |
89 | } | |
90 | fEvReco = 0; | |
91 | fTree = (TTree *) fFile->Get("TreeRecoEvent"); | |
92 | if (!fTree) { | |
93 | cout << "Tree of reconstructed events not found on file. Abort.\n"; | |
94 | gApplication->Terminate(0); | |
95 | } | |
2c0075de | 96 | fEvent = nevent; |
97 | fEmpty = kFALSE; | |
c7ba256d | 98 | TBranch *branch = fTree->GetBranch("Event"); |
99 | branch->SetAddress(&fEvReco); | |
100 | ||
101 | fEvGen = new AliMUONRecoEvent(); | |
102 | ||
7c85b268 | 103 | TFile *galiceFile = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); |
104 | galiceFile->cd(); | |
2c0075de | 105 | if (nevent > gAlice->TreeE()->GetEntries() - 1) { |
106 | cout << "Event number out of range !\n"; | |
107 | gApplication->Terminate(0); | |
108 | } | |
c7ba256d | 109 | gAlice->GetEvent(nevent); |
110 | ||
111 | ||
112 | fRecoTracks = 0; | |
113 | fGenTracks = 0; | |
114 | fPolyRecoList = 0; | |
115 | fPolyGenList = 0; | |
116 | fHighlited = -1; | |
117 | fMinMomentum = 0; | |
118 | fMaxMomentum = 999; | |
2c0075de | 119 | // map this event |
c7ba256d | 120 | MapEvent(nevent); |
121 | } | |
122 | ||
11ca64ac | 123 | //------------------------------------------------------------------- |
124 | AliMUONRecoDisplay::AliMUONRecoDisplay(const AliMUONRecoDisplay& rhs) | |
125 | : AliDisplay(rhs) | |
126 | { | |
127 | // Protected copy constructor | |
128 | ||
8c343c7c | 129 | AliFatal("Not implemented."); |
11ca64ac | 130 | } |
131 | ||
c7ba256d | 132 | //------------------------------------------------------------------- |
133 | AliMUONRecoDisplay::~AliMUONRecoDisplay() | |
134 | { | |
135 | // Destructor of display object | |
136 | if (fPolyRecoList) { | |
137 | fPolyRecoList->Delete(); | |
138 | delete fPolyRecoList; | |
139 | } | |
140 | if (fPolyGenList) { | |
141 | fPolyGenList->Delete(); | |
142 | delete fPolyGenList; | |
143 | } | |
144 | delete fEvGen; | |
145 | } | |
11ca64ac | 146 | |
147 | //------------------------------------------------------------------- | |
148 | AliMUONRecoDisplay& | |
149 | AliMUONRecoDisplay::operator=(const AliMUONRecoDisplay& rhs) | |
150 | { | |
151 | // Protected assignement operator | |
152 | ||
153 | if (this == &rhs) return *this; | |
154 | ||
8c343c7c | 155 | AliFatal("Not implemented."); |
11ca64ac | 156 | |
157 | return *this; | |
158 | } | |
159 | ||
2c0075de | 160 | //------------------------------------------------------------------- |
161 | Bool_t AliMUONRecoDisplay::Event(Int_t nevent) | |
162 | { | |
6a9bc541 | 163 | // Go to event nevent |
2c0075de | 164 | fEvent = nevent; |
165 | for (Int_t entry=0; entry<fTree->GetEntries(); entry++) { | |
166 | fTree->GetEntry(entry); | |
167 | if (fEvReco->GetNoEvent() == nevent) return kTRUE; | |
168 | } | |
169 | cout << "Event number " << nevent << " empty\n"; | |
170 | fRecoTracks = 0; | |
171 | fPolyRecoList = 0; | |
172 | fGenTracks = 0; | |
173 | fPolyGenList = 0; | |
174 | fHighlited = -1; | |
175 | fMinMomentum = 0; | |
176 | fMaxMomentum = 999; | |
177 | ||
178 | return kFALSE; | |
179 | } | |
c7ba256d | 180 | |
181 | //------------------------------------------------------------------- | |
182 | void AliMUONRecoDisplay::MapEvent(Int_t nevent) | |
183 | { | |
184 | // get generated event (nevent) from galice.root and corresponding | |
185 | // reconstructed event from tree_reco.root; | |
186 | cout << "mapping event " << nevent << endl; | |
2c0075de | 187 | fEvent = nevent; |
188 | fEmpty = kFALSE; | |
c7ba256d | 189 | fFile->cd(); |
2c0075de | 190 | // check if the event is not empty and make fEvReco to point to it |
191 | fEmpty = !Event(nevent); | |
c7ba256d | 192 | // just testing |
2c0075de | 193 | // if (!fEvReco) { |
194 | // cout << "failed to load reco event ! Correct this.\n"; | |
195 | // gApplication->Terminate(0); | |
196 | // } | |
197 | if (fEmpty) return; | |
c7ba256d | 198 | fRecoTracks = fEvReco->TracksPtr(); |
199 | fPolyRecoList = MakePolyLines3D(fRecoTracks); | |
200 | // clear previous event | |
201 | if (fEvGen) fEvGen->Clear(); | |
202 | fEvGen->SetNoEvent(nevent); | |
203 | // get list of particles | |
c7ba256d | 204 | // connect MUON module |
6a9bc541 | 205 | AliDetector *pMUON = gAlice->GetDetector("MUON"); |
206 | if (!pMUON) { | |
c7ba256d | 207 | cout << "MUON module not present.\n"; |
208 | gApplication->Terminate(0); | |
209 | } | |
210 | // get the number of generated tracks | |
88cb7938 | 211 | Int_t ntracks = (Int_t)pMUON->TreeH()->GetEntries(); |
c7ba256d | 212 | // Fill the fEvGen object |
213 | AliMUONRecoTrack *gtrack = 0; | |
214 | AliMUONHit *hit = 0; | |
215 | TParticle *particle; | |
216 | Int_t ch; | |
217 | // loop all tracks | |
218 | for (Int_t track=0; track<ntracks; track++) { | |
6a9bc541 | 219 | hit = (AliMUONHit *) pMUON->FirstHit(track); |
c7ba256d | 220 | if (!hit) continue; |
5d12ce38 | 221 | particle = gAlice->GetMCApp()->Particle(hit->Track()); |
c7ba256d | 222 | if (IsReconstructible(track) && TMath::Abs(particle->GetPdgCode())==13) { |
223 | gtrack = fEvGen->AddEmptyTrack(); | |
224 | gtrack->SetSign(TMath::Sign((Int_t)1, -particle->GetPdgCode())); | |
225 | // reset hits | |
226 | for (ch=0; ch<10; ch++) gtrack->SetHitPosition(ch,0,0,0); | |
227 | // loop all hits | |
6a9bc541 | 228 | for (AliMUONHit *muonHit=(AliMUONHit*)pMUON->FirstHit(track); |
c7ba256d | 229 | muonHit; |
6a9bc541 | 230 | muonHit=(AliMUONHit*)pMUON->NextHit()) { |
59f3ecff | 231 | ch = muonHit->Chamber() - 1; |
c7ba256d | 232 | if (ch<0 || ch>9) continue; |
233 | gtrack->SetHitPosition(ch, muonHit->X(), muonHit->Y(), muonHit->Z()); | |
234 | gtrack->SetMomReconstr(particle->Px(), particle->Py(), particle->Pz()); | |
235 | gtrack->SetVertexPos(particle->Vz()); | |
236 | gtrack->SetChi2r(0.0); | |
237 | } | |
238 | } | |
239 | } | |
240 | fGenTracks = fEvGen->TracksPtr(); | |
241 | fPolyGenList = MakePolyLines3D(fGenTracks); | |
242 | } | |
243 | //------------------------------------------------------------------- | |
244 | void AliMUONRecoDisplay::XYPlot() | |
245 | { | |
246 | // Plot reco. tracks hits in all chambers for current event: | |
247 | // - open blue squares : generated muons that were reconstructed accurate | |
248 | // - open cyan squares : generated muons that were not reconstructed | |
249 | // - filled green circles : reco. tracks (accurate) | |
250 | // - filled red squares : fake tracks | |
2c0075de | 251 | if (fEmpty) return; |
c7ba256d | 252 | Double_t kMaxRadius[10]; |
253 | kMaxRadius[0] = kMaxRadius[1] = 91.5; | |
254 | kMaxRadius[2] = kMaxRadius[3] = 122.5; | |
255 | kMaxRadius[4] = kMaxRadius[5] = 158.3; | |
256 | kMaxRadius[6] = kMaxRadius[7] = 260.0; | |
257 | kMaxRadius[8] = kMaxRadius[9] = 260.0; | |
258 | ||
7c85b268 | 259 | TH2F *xygenFound[10]; |
260 | TH2F *xygenLost[10]; | |
261 | TH2F *xyrecoGood[10]; | |
262 | TH2F *xyrecoFake[10]; | |
c7ba256d | 263 | Double_t x,y,r; |
264 | Int_t matches[500]; | |
265 | Int_t index, ch; | |
266 | ||
267 | TPad *pad = (TPad*)gROOT->GetSelectedPad(); | |
268 | TCanvas *canvas = new TCanvas("xy", "Reconstruction efficiency"); | |
269 | canvas->Clear(); | |
270 | canvas->cd(); | |
271 | canvas->Divide(3,4); | |
272 | canvas->cd(11); | |
273 | gPad->Delete(); | |
274 | canvas->cd(12); | |
275 | gPad->Delete(); | |
276 | canvas->cd(1); | |
277 | ||
278 | // Define histograms for x-y plots | |
279 | for (ch=0; ch<10; ch++) { | |
7c85b268 | 280 | xygenFound[ch] = new TH2F("xygen_found","",50,-kMaxRadius[ch],kMaxRadius[ch],50,-kMaxRadius[ch],kMaxRadius[ch]); |
281 | xygenLost[ch] = new TH2F("xygen_lost","",50,-kMaxRadius[ch],kMaxRadius[ch],50,-kMaxRadius[ch],kMaxRadius[ch]); | |
282 | xyrecoGood[ch] = new TH2F("xyreco_good","",50,-kMaxRadius[ch],kMaxRadius[ch],50,-kMaxRadius[ch],kMaxRadius[ch]); | |
283 | xyrecoFake[ch] = new TH2F("xyreco_fake","",50,-kMaxRadius[ch],kMaxRadius[ch],50,-kMaxRadius[ch],kMaxRadius[ch]); | |
c7ba256d | 284 | } |
285 | // find list of matching tracks | |
286 | fPrinted = kTRUE; // no need to print | |
287 | for (index=0; index<fRecoTracks->GetEntriesFast(); index++) { | |
288 | matches[index] = GetBestMatch(index); | |
289 | } | |
290 | // and fill them | |
291 | for (index=0; index<fGenTracks->GetEntriesFast(); index++) { // GEANT | |
292 | Bool_t wasreconst = kFALSE; | |
293 | for (Int_t i=0; i<fRecoTracks->GetEntriesFast(); i++) { | |
294 | if (matches[i] == index) { | |
295 | wasreconst = kTRUE; | |
296 | break; | |
297 | } | |
298 | } | |
299 | AliMUONRecoTrack *current = (AliMUONRecoTrack*)fGenTracks->UncheckedAt(index); | |
300 | for (ch=0; ch<10; ch++) { | |
301 | x = current->GetPosX(ch); | |
302 | y = current->GetPosY(ch); | |
303 | r = TMath::Sqrt(x*x +y*y); | |
304 | if (r >= 10) { | |
305 | if (wasreconst) { | |
7c85b268 | 306 | xygenFound[ch]->Fill(x,y); |
307 | } else {xygenLost[ch]->Fill(x,y);} | |
c7ba256d | 308 | } |
309 | } | |
310 | } | |
311 | for (index=0; index<fRecoTracks->GetEntriesFast(); index++) { // reco | |
312 | AliMUONRecoTrack *current = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(index); | |
313 | for (ch=0; ch<10; ch++) { | |
314 | x = current->GetPosX(ch); | |
315 | y = current->GetPosY(ch); | |
316 | r = TMath::Sqrt(x*x +y*y); | |
317 | if (r >= 10) { | |
318 | if (matches[index] >= 0) { | |
7c85b268 | 319 | xyrecoGood[ch]->Fill(x,y); |
320 | } else {xyrecoFake[ch]->Fill(x,y);} | |
c7ba256d | 321 | } |
322 | } | |
323 | } | |
324 | // finally plot them | |
325 | for (ch=0; ch<10; ch++) { | |
326 | canvas->cd(ch+1); | |
7c85b268 | 327 | xygenFound[ch]->SetMarkerColor(kBlue); |
328 | xygenFound[ch]->SetMarkerStyle(4); | |
329 | xygenFound[ch]->SetMarkerSize(0.5); | |
330 | xygenFound[ch]->SetStats(kFALSE); | |
331 | xygenFound[ch]->Draw(); | |
332 | xygenLost[ch]->SetMarkerColor(kCyan); | |
333 | xygenLost[ch]->SetMarkerStyle(4); | |
334 | xygenLost[ch]->SetMarkerSize(0.5); | |
335 | xygenLost[ch]->SetStats(kFALSE); | |
336 | xygenLost[ch]->Draw("SAME"); | |
337 | xyrecoGood[ch]->SetMarkerColor(kGreen); | |
338 | xyrecoGood[ch]->SetMarkerStyle(20); | |
339 | xyrecoGood[ch]->SetMarkerSize(0.4); | |
340 | xyrecoGood[ch]->SetStats(kFALSE); | |
341 | xyrecoGood[ch]->Draw("SAME"); | |
342 | xyrecoFake[ch]->SetMarkerColor(kRed); | |
343 | xyrecoFake[ch]->SetMarkerStyle(20); | |
344 | xyrecoFake[ch]->SetMarkerSize(0.5); | |
345 | xyrecoFake[ch]->SetStats(kFALSE); | |
346 | xyrecoFake[ch]->Draw("SAME"); | |
c7ba256d | 347 | } |
348 | canvas->SetTitle("y vs. x for simulated and reconstructed tracks"); | |
349 | pad->cd(); | |
350 | } | |
351 | ||
352 | //------------------------------------------------------------------- | |
353 | void AliMUONRecoDisplay::RecoEfficiency(Int_t first, Int_t last) | |
354 | { | |
355 | // Loop selected reconstructed events, compute efficiency as total number of | |
356 | // reconstructed tracks (accurate) over total number of generated | |
357 | // reconstructible muons. Make histogram for momentum precision and profile | |
358 | // of mean momentum precision versus total momentum (generated) | |
2c0075de | 359 | Int_t nevents = (Int_t)gAlice->TreeE()->GetEntries(); |
c7ba256d | 360 | if (last > nevents) last = nevents - 1; |
361 | if (first < 0) first = 0; | |
362 | nevents = last - first + 1; | |
363 | Int_t track; | |
364 | Float_t efficiency; | |
365 | Float_t generated=0, found=0, fake=0; // number of generated/found/fake tracks | |
366 | Double_t pgen, preco, dP; // dP=preco-pgen | |
367 | AliMUONRecoTrack *rtrack, *gtrack; // generated/reco. current tracks | |
368 | fPrinted = kTRUE; // no need to print | |
369 | Int_t currentEvent = gAlice->GetHeader()->GetEvent(); | |
370 | ||
371 | TH1F *pReso = new TH1F("Momentum precision", "dP = Prec - Pgen", 100, -5.0, 5.0); | |
372 | pReso->SetXTitle("dP [GeV/c]"); | |
373 | pReso->SetYTitle("dN/dP"); | |
374 | ||
375 | TProfile *dPp = new TProfile("", "dP vs. P", 50, 0, 100, 0, 5); | |
376 | dPp->SetXTitle("P [GeV/c]"); | |
377 | dPp->SetYTitle("<dP> [GeV/c]"); | |
378 | ||
379 | TPad *pad = (TPad*)gROOT->GetSelectedPad(); | |
380 | TCanvas *canvas = new TCanvas("c", "Reconstruction efficiency"); | |
381 | canvas->Divide(1,2); | |
382 | canvas->cd(1); | |
383 | ||
384 | // loop events | |
385 | for (Int_t event=first; event<=last; event++) { | |
386 | // get the reco. & gen. events | |
7c85b268 | 387 | TFile *galiceFile = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); |
388 | galiceFile->cd(); | |
c7ba256d | 389 | gAlice->GetEvent(event); |
390 | MapEvent(event); | |
2c0075de | 391 | if (fEmpty) { |
392 | // skip empty events | |
393 | fEmpty = kFALSE; | |
394 | continue; | |
395 | } | |
c7ba256d | 396 | generated += fGenTracks->GetEntriesFast(); |
397 | // loop reco. tracks | |
398 | for (track=0; track<fRecoTracks->GetEntriesFast(); track++) { | |
399 | rtrack = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(track); | |
400 | preco = rtrack->P(); | |
401 | // find correspondent generated track | |
402 | Int_t ind = GetBestMatch(track); | |
403 | if (ind >=0) { | |
404 | gtrack = (AliMUONRecoTrack*)fGenTracks->UncheckedAt(ind); | |
405 | pgen = gtrack->P(); | |
406 | found += 1; | |
407 | dP = preco - pgen; | |
408 | pReso->Fill(dP); | |
409 | dPp->Fill(pgen, TMath::Abs(dP)); | |
410 | } else { | |
411 | fake += 1; | |
412 | } | |
413 | } | |
414 | } | |
415 | efficiency = found/generated; | |
416 | cout << "=================================================================\n"; | |
417 | cout << "|| Reconstruction efficiency and momentum precision ||\n"; | |
418 | cout << "=================================================================\n"; | |
419 | cout << "Number of events processed " << nevents << endl; | |
420 | cout << "Total number of reconstructible muons : " << (Int_t)generated << endl; | |
421 | cout << "RECONSTRUCTION EFFICIENCY : " << efficiency*100 << " %" << endl; | |
422 | cout << "Fake track rate : " << 100*fake/generated << " %" << endl; | |
423 | cout << "Momentum precision fit : \n" << endl; | |
424 | pReso->Fit("gaus"); | |
425 | cout << "=================================================================\n"; | |
426 | ||
427 | canvas->cd(1); | |
428 | pReso->Draw(); | |
429 | canvas->cd(2); | |
430 | dPp->Draw("HIST"); | |
431 | pad->cd(); | |
432 | ShowNextEvent(currentEvent-last); | |
433 | } | |
434 | ||
435 | //------------------------------------------------------------------- | |
436 | void AliMUONRecoDisplay::ShowNextEvent(Int_t delta) | |
437 | { | |
438 | // overwritten from AliDisplay in order to get also mapping of next event | |
7c85b268 | 439 | TFile *galiceFile = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root"); |
440 | galiceFile->cd(); | |
c7ba256d | 441 | if (delta) { |
442 | gAlice->Clear(); | |
443 | Int_t currentEvent = gAlice->GetHeader()->GetEvent(); | |
444 | Int_t newEvent = currentEvent + delta; | |
2c0075de | 445 | if (newEvent<0 || newEvent>(gAlice->TreeE()->GetEntries() - 1)) return; |
c7ba256d | 446 | Int_t nparticles = gAlice->GetEvent(newEvent); |
447 | cout << "Event : " << newEvent << " with " << nparticles << " particles\n"; | |
88cb7938 | 448 | |
449 | /******************************************************************/ | |
450 | AliConfig* config = AliConfig::Instance(); | |
451 | TFolder* topfold = (TFolder*)config->GetTopFolder(); | |
452 | if (topfold == 0x0) | |
453 | { | |
8c343c7c | 454 | AliError("Can not get Alice top folder"); |
88cb7938 | 455 | return; |
456 | } | |
457 | TString fmdfoldname(config->GetDataFolderName()+"/"+"MUON"); | |
458 | TFolder* fmdfold = (TFolder*)topfold->FindObject(fmdfoldname); | |
459 | if (fmdfold == 0x0) | |
460 | { | |
8c343c7c | 461 | AliError("Can not get MUON folder"); |
88cb7938 | 462 | return; |
463 | } | |
464 | TTree* treeH = dynamic_cast<TTree*>(fmdfold->FindObject("TreeH")); | |
465 | if (treeH == 0x0) | |
466 | { | |
8c343c7c | 467 | AliError("Can not get TreeH"); |
88cb7938 | 468 | return; |
469 | } | |
470 | /******************************************************************/ | |
c7ba256d | 471 | MapEvent(newEvent); |
472 | fHighlited = -1; | |
473 | } | |
474 | LoadPoints(); | |
475 | fPad->cd(); | |
476 | Draw(); | |
477 | if (gROOT->GetListOfCanvases()->FindObject("xy")) XYPlot(); | |
478 | } | |
479 | //------------------------------------------------------------------- | |
30178c30 | 480 | Bool_t AliMUONRecoDisplay::IsReconstructible(Int_t track) const |
c7ba256d | 481 | { |
482 | // true if at least three hits in first 2 stations, 3 in last 2 stations | |
483 | // and one in station 3 | |
2c0075de | 484 | if (fEmpty) return kFALSE; |
6a9bc541 | 485 | AliDetector *pMUON = gAlice->GetDetector("MUON"); |
c7ba256d | 486 | Bool_t chHit[10]; |
487 | Int_t ch; | |
488 | for (ch=0; ch<10; ch++) chHit[ch] = kFALSE; | |
489 | //loop hits | |
6a9bc541 | 490 | for (AliMUONHit *muonHit=(AliMUONHit*)pMUON->FirstHit(track); |
c7ba256d | 491 | muonHit; |
6a9bc541 | 492 | muonHit=(AliMUONHit*)pMUON->NextHit()) { |
59f3ecff | 493 | ch = muonHit->Chamber() - 1; |
c7ba256d | 494 | if (ch<0 || ch>9) continue; |
495 | chHit[ch] = kTRUE; | |
496 | } | |
497 | Int_t nhits = 0; | |
498 | for (ch=0; ch<4; ch++) nhits += (chHit[ch])?1:0; | |
499 | if (nhits < 3) return kFALSE; | |
500 | nhits = 0; | |
501 | for (ch=4; ch<6; ch++) nhits+= (chHit[ch])?1:0; | |
502 | if (nhits < 1) return kFALSE; | |
503 | ||
504 | for (ch=7; ch<10; ch++) nhits+= (chHit[ch])?1:0; | |
505 | if (nhits < 3) return kFALSE; | |
506 | return kTRUE; | |
507 | } | |
508 | ||
509 | //------------------------------------------------------------------- | |
510 | void AliMUONRecoDisplay::DrawView(Float_t theta, Float_t phi, Float_t psi) | |
511 | { | |
512 | // ovewritten from base class to change the range for MUON | |
513 | gPad->SetCursor(kWatch); | |
514 | gPad->SetFillColor(1); | |
515 | gPad->Clear(); | |
516 | ||
517 | Int_t iret; | |
518 | TView *view = new TView(1); | |
519 | Float_t range = fRrange*fRangeSlider->GetMaximum(); | |
520 | view->SetRange(-range, -range, 0, range, range, 5*range); | |
521 | fZoomX0[0] = -1; | |
522 | fZoomY0[0] = -1; | |
523 | fZoomX1[0] = 1; | |
524 | fZoomY1[0] = 1; | |
525 | // Display Alice geometry | |
526 | gAlice->GetGeometry()->Draw("same"); | |
527 | //Loop on all detectors to add their products to the pad | |
528 | DrawHits(); | |
529 | // add itself to the list (last) | |
530 | AppendPad(); | |
531 | ||
532 | view->SetView(phi, theta, psi, iret); | |
533 | } | |
534 | //------------------------------------------------------------------- | |
535 | void AliMUONRecoDisplay::SetDrawHits(Bool_t hits) | |
536 | { | |
537 | // Turns on/off Geant hits drawing | |
538 | fDrawHits = hits; | |
539 | DrawView(0,-90); | |
540 | } | |
541 | ||
542 | //------------------------------------------------------------------- | |
543 | void AliMUONRecoDisplay::CutMomentum(Double_t min, Double_t max) | |
544 | { | |
545 | // Define momentum cut for displayed reconstructed tracks | |
546 | fMinMomentum = min; | |
547 | fMaxMomentum = max; | |
548 | if (fHighlited >= 0) UnHighlight(); | |
549 | DrawView(0,-90); | |
550 | } | |
551 | ||
552 | //------------------------------------------------------------------- | |
553 | Int_t AliMUONRecoDisplay::GetBestMatch(Int_t indr, Float_t tolerance) | |
554 | { | |
555 | // Find the index of best Geant track matching a reconstructed track : track | |
556 | // with maximum number of compatible hits (within tolerance*sigma bending and | |
557 | // non-bending resolution) and minimum number of fake hits; | |
558 | // If no match is found within a given tolerance, the method is called recursively | |
559 | // with increasing tolerance, until tolerance = 10; | |
2c0075de | 560 | if (fEmpty) return -1; |
c7ba256d | 561 | if (indr<0 || indr>=fRecoTracks->GetEntriesFast()) return -1; |
562 | AliMUONRecoTrack *rtrack = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(indr); | |
563 | AliMUONRecoTrack *gtrack = 0; | |
564 | Int_t bestMatch = -1; | |
565 | Int_t maxNcompat = 0; | |
566 | Int_t minNfakes = 10; | |
567 | Double_t xrhit,yrhit,radius,xghit,yghit,dX,dY; | |
568 | // loop over all Geant tracks | |
569 | for (Int_t indg=0; indg<fGenTracks->GetEntriesFast(); indg++) { | |
570 | gtrack = (AliMUONRecoTrack*)fGenTracks->UncheckedAt(indg); | |
2c0075de | 571 | if (!gtrack) continue; |
c7ba256d | 572 | Int_t ncompat = 0; // number of compat. hits for this track |
573 | Int_t nfake = 0; // number of fakes | |
574 | // loop chambers to find compatible hits | |
575 | for (Int_t ch=0; ch<10; ch++) { | |
576 | xrhit = rtrack->GetPosX(ch); | |
577 | yrhit = rtrack->GetPosY(ch); | |
578 | radius = TMath::Sqrt(xrhit*xrhit + yrhit*yrhit); | |
579 | if (radius<10) continue; // skip null hits | |
580 | xghit = gtrack->GetPosX(ch); | |
581 | yghit = gtrack->GetPosY(ch); | |
582 | dX = TMath::Abs(xghit-xrhit); | |
583 | dY = TMath::Abs(yghit-yrhit); | |
2c0075de | 584 | if (dX<tolerance*0.144 && dY<tolerance*0.01) {// within tol*sigma resolution |
c7ba256d | 585 | ncompat++; |
586 | continue; // compatible hit | |
587 | } else nfake++; // fake hit | |
588 | } | |
589 | if (ncompat && ncompat>=maxNcompat && nfake<minNfakes) { // this is best matching | |
590 | maxNcompat = ncompat; | |
591 | minNfakes = nfake; | |
592 | bestMatch = indg; | |
593 | } | |
594 | } | |
595 | if (bestMatch<0 && tolerance<=9.) bestMatch = GetBestMatch(indr, tolerance+=1); | |
596 | if (!fPrinted) { | |
597 | rtrack = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(indr); | |
598 | Int_t sign = rtrack->GetSign(); | |
599 | cout << "Reconstructed track : " << indr << "(" << sign << ")" << endl; | |
600 | rtrack->TrackInfo(); | |
601 | printf("Best matching Geant track within %i*sgm : %i\n", (Int_t)tolerance, bestMatch); | |
602 | if (bestMatch >=0) { | |
603 | gtrack = (AliMUONRecoTrack*)fGenTracks->UncheckedAt(bestMatch); | |
604 | gtrack->TrackInfo(); | |
605 | } | |
606 | cout << "-----------------------------------------------------------------\n"; | |
607 | } | |
608 | fPrinted = kTRUE; | |
609 | ||
610 | return bestMatch; | |
611 | } | |
612 | ||
613 | //------------------------------------------------------------------- | |
614 | void AliMUONRecoDisplay::Highlight(Int_t track) | |
615 | { | |
2c0075de | 616 | // Highlight the specified track |
617 | if (fEmpty) return; | |
c7ba256d | 618 | if (fHighlited >=0) UnHighlight(); |
619 | if (track<0 || track>fPolyRecoList->GetEntries()) return; | |
620 | TPolyLine3D *line = (TPolyLine3D*)fPolyRecoList->UncheckedAt(track); | |
621 | line->SetLineColor(kYellow); | |
622 | line->SetLineWidth(1); | |
623 | fHighlited = track; | |
2c0075de | 624 | // DrawView(15,-45,135); |
625 | fPad->cd(); | |
626 | Draw(); | |
c7ba256d | 627 | } |
628 | ||
629 | //------------------------------------------------------------------- | |
630 | void AliMUONRecoDisplay::UnHighlight() | |
631 | { | |
632 | // Unhighlight a previous highlighted track | |
2c0075de | 633 | if (fHighlited < 0 || fEmpty) return; // nothing to do |
c7ba256d | 634 | TPolyLine3D *line = (TPolyLine3D*)fPolyRecoList->UncheckedAt(fHighlited); |
635 | line->SetLineColor(kRed); | |
636 | line->SetLineWidth(1); | |
637 | fHighlited = -1; | |
2c0075de | 638 | // DrawView(0,-90); |
639 | fPad->cd(); | |
640 | Draw(); | |
c7ba256d | 641 | } |
642 | ||
643 | //------------------------------------------------------------------- | |
644 | void AliMUONRecoDisplay::DrawHits() | |
645 | { | |
646 | // Draw hits for all ALICE detectors. Overwrites the DrawHits() method of the | |
647 | // base class for reco. track drawing | |
648 | ||
649 | Float_t cutmin, cutmax, etamin, etamax, pmom, smin, smax, eta, theta, r; | |
650 | const Float_t kptcutmax = 2; | |
651 | const Float_t ketacutmax = 1.5; | |
652 | Float_t *pxyz; | |
653 | Int_t ntracks,track; | |
654 | TParticle *particle; | |
655 | TObjArray *points; | |
656 | AliPoints *pm; | |
657 | ||
658 | //Get cut slider | |
659 | smax = fCutSlider->GetMaximum(); | |
660 | smin = fCutSlider->GetMinimum(); | |
661 | cutmin = kptcutmax*smin; | |
662 | if (smax < 0.98) cutmax = kptcutmax*smax; | |
663 | else cutmax = 100000; | |
664 | ||
665 | //Get eta slider | |
666 | smax = fEtaSlider->GetMaximum(); | |
667 | smin = fEtaSlider->GetMinimum(); | |
668 | etamin = ketacutmax*(2*smin-1); | |
669 | etamax = ketacutmax*(2*smax-1); | |
670 | if (smin < 0.02) etamin = -1000; | |
671 | if (smax > 0.98) etamax = 1000; | |
672 | ||
673 | TIter next(gAlice->Modules()); | |
674 | AliModule *module; | |
675 | fHitsCuts = 0; | |
676 | if (fDrawHits) { | |
677 | // draw hits in all modules | |
678 | while((module = (AliModule*)next())) { | |
679 | if (!module->IsActive()) continue; | |
680 | points = module->Points(); | |
681 | if (!points) continue; | |
682 | ntracks = points->GetEntriesFast(); | |
683 | for (track=0;track<ntracks;track++) { | |
684 | pm = (AliPoints*)points->UncheckedAt(track); | |
685 | if (!pm) continue; | |
686 | particle = pm->GetParticle(); | |
687 | if (!particle) continue; | |
688 | pmom = particle->P(); | |
689 | if (pmom < cutmin) continue; | |
690 | if (pmom > cutmax) continue; | |
691 | // as a first approximation, take eta of first point | |
692 | pxyz = pm->GetP(); | |
693 | r = TMath::Sqrt(pxyz[0]*pxyz[0] + pxyz[1]*pxyz[1]); | |
694 | theta = TMath::ATan2(r,TMath::Abs(pxyz[2])); | |
695 | if(theta) eta = -TMath::Log(TMath::Tan(0.5*theta)); else eta = 1e10; | |
696 | if (pxyz[2] < 0) eta = -eta; | |
697 | if (eta < etamin || eta > etamax) continue; | |
698 | pm->Draw(); | |
699 | fHitsCuts += pm->GetN(); | |
700 | } | |
701 | } | |
702 | } | |
703 | // draw reconstructed tracks | |
2c0075de | 704 | if (fEmpty) return; |
c7ba256d | 705 | TPolyLine3D *line, *gline; |
706 | Int_t bestMatch; | |
707 | Double_t px,py,pz,p; | |
708 | AliMUONRecoTrack *rtrack; | |
709 | ||
710 | if (fHighlited >= 0) { | |
711 | line = (TPolyLine3D*)fPolyRecoList->UncheckedAt(fHighlited); | |
712 | line->Draw(); | |
713 | fPrinted = kFALSE; | |
714 | bestMatch = GetBestMatch(fHighlited); | |
715 | if (bestMatch>=0) { | |
716 | gline = (TPolyLine3D*)fPolyGenList->UncheckedAt(bestMatch); | |
717 | gline->SetLineColor(kRed); | |
718 | gline->SetLineWidth(2); | |
719 | gline->SetLineStyle(2); | |
720 | gline->Draw(); | |
721 | } | |
722 | } else { | |
723 | for (track=0; track<fPolyRecoList->GetEntries(); track++) { | |
724 | rtrack = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(track); | |
725 | px = rtrack->GetMomReconstr(0); | |
726 | py = rtrack->GetMomReconstr(1); | |
727 | pz = rtrack->GetMomReconstr(2); | |
728 | p = rtrack->P(); | |
729 | if (p>fMinMomentum && p<fMaxMomentum) { | |
730 | line = (TPolyLine3D*)fPolyRecoList->UncheckedAt(track); | |
731 | line->Draw(); | |
732 | } | |
733 | } | |
734 | } | |
735 | } | |
736 | ||
737 | //------------------------------------------------------------------- | |
738 | void AliMUONRecoDisplay::ListTracks() | |
739 | { | |
740 | // List momentum information of all reconstructed traccks within fPmin and fPmax | |
741 | // cuts, as well as their best matching Geant tracks | |
2c0075de | 742 | if (fEmpty) return; |
c7ba256d | 743 | cout << "================================================================\n"; |
744 | printf("Reconstructed tracks with momentum in range : %g , %g [GeV/c]\n", | |
745 | fMinMomentum, fMaxMomentum); | |
746 | cout << "----------------------------------------------------------------\n"; | |
747 | AliMUONRecoTrack *rtrack; | |
748 | Double_t p; | |
749 | Int_t sign; | |
750 | for (Int_t ind=0; ind<fRecoTracks->GetEntries(); ind++) { | |
751 | rtrack = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(ind); | |
752 | p = rtrack->P(); | |
753 | if (p>fMinMomentum && p<fMaxMomentum) { | |
754 | fPrinted = kFALSE; | |
755 | GetBestMatch(ind); | |
756 | sign = rtrack->GetSign(); | |
757 | } | |
758 | } | |
759 | cout << "================================================================\n"; | |
760 | } | |
761 | ||
762 | //------------------------------------------------------------------- | |
763 | TClonesArray* AliMUONRecoDisplay::MakePolyLines3D(TClonesArray *tracklist) | |
764 | { | |
765 | // Makes the list of polylines3D corresponding to the list of tracks | |
2c0075de | 766 | if (fEmpty) return 0; |
c7ba256d | 767 | if (tracklist!=fRecoTracks && tracklist!=fGenTracks) return 0; |
768 | Bool_t reco = (tracklist==fRecoTracks)?kTRUE:kFALSE; | |
769 | // make sure there is no other list in memory | |
770 | if (reco) { | |
771 | if (fPolyRecoList) { | |
772 | fPolyRecoList->Delete(); | |
773 | delete fPolyRecoList; | |
774 | fPolyRecoList = 0; | |
775 | } | |
776 | } else { | |
777 | if (fPolyGenList) { | |
778 | fPolyGenList->Delete(); | |
779 | delete fPolyGenList; | |
780 | fPolyGenList = 0; | |
781 | } | |
782 | } | |
783 | if (!tracklist->GetEntries()) return 0; | |
784 | ||
785 | AliMUONRecoTrack* track = 0; | |
786 | TClonesArray *polyLines3D = new TClonesArray("TPolyLine3D",1000); | |
787 | TClonesArray &polylist = *polyLines3D; | |
788 | TPolyLine3D *polyline = 0; | |
789 | ||
790 | // loop all tracks | |
791 | for (Int_t i=0; i<tracklist->GetEntries(); i++) { | |
792 | track = (AliMUONRecoTrack*)tracklist->UncheckedAt(i); | |
793 | Int_t ch = 0; | |
794 | Double_t x,y,z,r; | |
795 | polyline = new(polylist[i]) TPolyLine3D(2,""); | |
796 | polyline->SetLineColor(kRed); | |
797 | polyline->SetLineWidth(1); | |
798 | polyline->SetNextPoint(0,0,track->GetVertexPos()); // vertex point | |
799 | // loop chambers to fill TPolyLine3D objects | |
800 | for (ch=0; ch<10; ch++) { | |
801 | x = track->GetPosX(ch); | |
802 | y = track->GetPosY(ch); | |
803 | r = TMath::Sqrt(x*x + y*y); | |
804 | if (r < 10) continue; | |
805 | z = track->GetPosZ(ch); | |
806 | polyline->SetNextPoint(x,y,z); | |
807 | } | |
808 | } | |
809 | return polyLines3D; | |
810 | } | |
811 | ||
812 | //------------------------------------------------------------------- | |
813 | void AliMUONRecoDisplay::PolyLineInfo(TClonesArray *line3Dlist) | |
814 | { | |
815 | // Prints information (x, y, z coordinates) for all constructed polylines | |
2c0075de | 816 | if (fEmpty) return; |
c7ba256d | 817 | if (line3Dlist) { |
818 | TPolyLine3D *polyline = 0; | |
819 | for(Int_t trackIndex=0; trackIndex<line3Dlist->GetEntries(); trackIndex++) { | |
820 | polyline = (TPolyLine3D*)line3Dlist->UncheckedAt(trackIndex); | |
821 | polyline->ls(); | |
822 | Float_t *pl = polyline->GetP(); | |
823 | for (Int_t i=0; i<polyline->GetN() ;i++) { | |
824 | printf(" x[%d]=%g, y[%d]=%g, z[%d]=%g\n",i,pl[3*i],i,pl[3*i+1],i,pl[3*i+2]); | |
825 | } | |
826 | } | |
827 | } | |
828 | } |