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