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