Modifications needed by the HBT analysis (P.Skowronski)
[u/mrichter/AliRoot.git] / MUON / AliMUONRecoDisplay.cxx
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 /* $Id$ */
17
18 //Authors: Mihaela Gheata, Andrei Gheata 09/10/00
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 //////////////////////////////////////////////////////////////////////
46
47 #include <Riostream.h>
48
49 #include <TROOT.h>
50 #include <TFolder.h>
51 #include <TClonesArray.h>
52 #include <TGeometry.h>
53 #include <TSlider.h>
54 #include <TGeometry.h>
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"
64 #include "AliMC.h"
65
66 ClassImp(AliMUONRecoDisplay)
67
68 //-------------------------------------------------------------------
69 AliMUONRecoDisplay::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    }
85    fEvent = nevent;
86    fEmpty = kFALSE;
87    TBranch *branch = fTree->GetBranch("Event");
88    branch->SetAddress(&fEvReco);
89    
90    fEvGen  = new AliMUONRecoEvent();
91
92    TFile *galiceFile = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
93    galiceFile->cd();
94    if (nevent > gAlice->TreeE()->GetEntries() - 1) {
95       cout << "Event number out of range !\n";
96       gApplication->Terminate(0);
97    }
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;
108    // map this event
109    MapEvent(nevent);
110 }
111
112 //-------------------------------------------------------------------
113 AliMUONRecoDisplay::~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 }
126 //-------------------------------------------------------------------
127 Bool_t AliMUONRecoDisplay::Event(Int_t nevent)
128 {
129 // Go to event nevent
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 }
146
147 //-------------------------------------------------------------------
148 void 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;
153    fEvent = nevent;
154    fEmpty = kFALSE;
155    fFile->cd();
156    // check if the event is not empty and make fEvReco to point to it
157    fEmpty = !Event(nevent);
158    // just testing
159 //   if (!fEvReco) {
160 //      cout << "failed to load reco event ! Correct this.\n";
161 //      gApplication->Terminate(0);
162 //   }
163    if (fEmpty) return;
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
170    // connect MUON module
171    AliDetector *pMUON = gAlice->GetDetector("MUON");
172    if (!pMUON) {
173       cout << "MUON module not present.\n";
174       gApplication->Terminate(0);
175    }
176    // get the number of generated tracks
177    Int_t ntracks = (Int_t)pMUON->TreeH()->GetEntries();
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++) {
185       hit = (AliMUONHit *) pMUON->FirstHit(track);
186       if (!hit) continue;
187       particle = gAlice->GetMCApp()->Particle(hit->Track());
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
194          for (AliMUONHit *muonHit=(AliMUONHit*)pMUON->FirstHit(track);
195               muonHit;
196               muonHit=(AliMUONHit*)pMUON->NextHit()) {
197             ch = muonHit->Chamber() - 1;
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 //-------------------------------------------------------------------
210 void 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 
217    if (fEmpty) return;
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
225    TH2F *xygenFound[10]; 
226    TH2F *xygenLost[10]; 
227    TH2F *xyrecoGood[10];
228    TH2F *xyrecoFake[10];
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++) {
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]);
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) {
272                xygenFound[ch]->Fill(x,y);
273             } else {xygenLost[ch]->Fill(x,y);}
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) {
285                xyrecoGood[ch]->Fill(x,y);
286             } else {xyrecoFake[ch]->Fill(x,y);}
287          }
288       }
289    }
290    // finally plot them
291    for (ch=0; ch<10; ch++) {
292       canvas->cd(ch+1);   
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");
313    }
314    canvas->SetTitle("y vs. x for simulated and reconstructed tracks");
315    pad->cd();
316 }
317
318 //-------------------------------------------------------------------
319 void 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)
325    Int_t nevents = (Int_t)gAlice->TreeE()->GetEntries();
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
353       TFile *galiceFile = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
354       galiceFile->cd();
355       gAlice->GetEvent(event);
356       MapEvent(event);
357       if (fEmpty) {
358       // skip empty events
359          fEmpty = kFALSE;
360          continue;
361       }
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 //-------------------------------------------------------------------
402 void AliMUONRecoDisplay::ShowNextEvent(Int_t delta)
403 {
404 // overwritten from AliDisplay in order to get also mapping of next event
405    TFile *galiceFile = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
406    galiceFile->cd();
407    if (delta) {
408       gAlice->Clear();
409       Int_t currentEvent = gAlice->GetHeader()->GetEvent();
410       Int_t newEvent     = currentEvent + delta;
411       if (newEvent<0 || newEvent>(gAlice->TreeE()->GetEntries() - 1)) return;
412       Int_t nparticles = gAlice->GetEvent(newEvent);
413       cout << "Event : " << newEvent << " with " << nparticles << " particles\n";
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 /******************************************************************/     
437       MapEvent(newEvent);
438       fHighlited = -1;
439    }
440    LoadPoints();
441    fPad->cd();
442    Draw();
443    if (gROOT->GetListOfCanvases()->FindObject("xy")) XYPlot();
444 }
445 //-------------------------------------------------------------------
446 Bool_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
450    if (fEmpty) return kFALSE;
451    AliDetector *pMUON = gAlice->GetDetector("MUON");
452    Bool_t chHit[10];
453    Int_t ch;
454    for (ch=0; ch<10; ch++) chHit[ch] = kFALSE;
455    //loop hits
456    for (AliMUONHit *muonHit=(AliMUONHit*)pMUON->FirstHit(track);
457         muonHit;
458         muonHit=(AliMUONHit*)pMUON->NextHit()) {
459       ch = muonHit->Chamber() - 1;
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 //-------------------------------------------------------------------
476 void 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 //-------------------------------------------------------------------
501 void AliMUONRecoDisplay::SetDrawHits(Bool_t hits)
502 {
503 // Turns on/off Geant hits drawing
504    fDrawHits = hits;
505    DrawView(0,-90);
506 }
507
508 //-------------------------------------------------------------------
509 void 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 //-------------------------------------------------------------------
519 Int_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;
526    if (fEmpty) return -1;
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);
537       if (!gtrack) continue;
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);
550          if (dX<tolerance*0.144 && dY<tolerance*0.01) {// within tol*sigma resolution
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 //-------------------------------------------------------------------
580 void AliMUONRecoDisplay::Highlight(Int_t track)
581 {
582 // Highlight the specified track
583    if (fEmpty) return; 
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;
590 //   DrawView(15,-45,135);
591    fPad->cd();
592    Draw();
593 }
594
595 //-------------------------------------------------------------------
596 void AliMUONRecoDisplay::UnHighlight()
597 {
598 // Unhighlight a previous highlighted track
599    if (fHighlited < 0 || fEmpty) return;      // nothing to do
600    TPolyLine3D *line = (TPolyLine3D*)fPolyRecoList->UncheckedAt(fHighlited);
601    line->SetLineColor(kRed);
602    line->SetLineWidth(1);
603    fHighlited = -1;
604 //   DrawView(0,-90);
605    fPad->cd();
606    Draw();
607 }
608
609 //-------------------------------------------------------------------
610 void 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
670    if (fEmpty) return;
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 //-------------------------------------------------------------------
704 void 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
708    if (fEmpty) return;
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 //-------------------------------------------------------------------
729 TClonesArray* AliMUONRecoDisplay::MakePolyLines3D(TClonesArray *tracklist)
730 {
731 // Makes the list of polylines3D corresponding to the list of tracks
732    if (fEmpty) return 0;
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 //-------------------------------------------------------------------
779 void AliMUONRecoDisplay::PolyLineInfo(TClonesArray *line3Dlist)
780 {
781 // Prints information (x, y, z coordinates) for all constructed polylines
782    if (fEmpty) return;
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