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