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