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