- Welding section on absorber side (LHCVC2C_001)
[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 "AliMUONRecoDisplay.h"
63 #include "AliMUONRecoEvent.h"
64 #include "AliMUONRecoTrack.h"
65 #include "AliMC.h"
66
67 ClassImp(AliMUONRecoDisplay)
68
69 //-------------------------------------------------------------------
70 AliMUONRecoDisplay::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    }
86    fEvent = nevent;
87    fEmpty = kFALSE;
88    TBranch *branch = fTree->GetBranch("Event");
89    branch->SetAddress(&fEvReco);
90    
91    fEvGen  = new AliMUONRecoEvent();
92
93    TFile *galiceFile = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
94    galiceFile->cd();
95    if (nevent > gAlice->TreeE()->GetEntries() - 1) {
96       cout << "Event number out of range !\n";
97       gApplication->Terminate(0);
98    }
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;
109    // map this event
110    MapEvent(nevent);
111 }
112
113 //-------------------------------------------------------------------
114 AliMUONRecoDisplay::AliMUONRecoDisplay(const AliMUONRecoDisplay& rhs)
115   : AliDisplay(rhs)
116 {
117 // Protected copy constructor
118
119   Fatal("AliMUONRecoDisplay", "Not implemented.");
120 }
121
122 //-------------------------------------------------------------------
123 AliMUONRecoDisplay::~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 }
136
137 //-------------------------------------------------------------------
138 AliMUONRecoDisplay&  
139 AliMUONRecoDisplay::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           
150 //-------------------------------------------------------------------
151 Bool_t AliMUONRecoDisplay::Event(Int_t nevent)
152 {
153 // Go to event nevent
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 }
170
171 //-------------------------------------------------------------------
172 void 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;
177    fEvent = nevent;
178    fEmpty = kFALSE;
179    fFile->cd();
180    // check if the event is not empty and make fEvReco to point to it
181    fEmpty = !Event(nevent);
182    // just testing
183 //   if (!fEvReco) {
184 //      cout << "failed to load reco event ! Correct this.\n";
185 //      gApplication->Terminate(0);
186 //   }
187    if (fEmpty) return;
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
194    // connect MUON module
195    AliDetector *pMUON = gAlice->GetDetector("MUON");
196    if (!pMUON) {
197       cout << "MUON module not present.\n";
198       gApplication->Terminate(0);
199    }
200    // get the number of generated tracks
201    Int_t ntracks = (Int_t)pMUON->TreeH()->GetEntries();
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++) {
209       hit = (AliMUONHit *) pMUON->FirstHit(track);
210       if (!hit) continue;
211       particle = gAlice->GetMCApp()->Particle(hit->Track());
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
218          for (AliMUONHit *muonHit=(AliMUONHit*)pMUON->FirstHit(track);
219               muonHit;
220               muonHit=(AliMUONHit*)pMUON->NextHit()) {
221             ch = muonHit->Chamber() - 1;
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 //-------------------------------------------------------------------
234 void 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 
241    if (fEmpty) return;
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
249    TH2F *xygenFound[10]; 
250    TH2F *xygenLost[10]; 
251    TH2F *xyrecoGood[10];
252    TH2F *xyrecoFake[10];
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++) {
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]);
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) {
296                xygenFound[ch]->Fill(x,y);
297             } else {xygenLost[ch]->Fill(x,y);}
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) {
309                xyrecoGood[ch]->Fill(x,y);
310             } else {xyrecoFake[ch]->Fill(x,y);}
311          }
312       }
313    }
314    // finally plot them
315    for (ch=0; ch<10; ch++) {
316       canvas->cd(ch+1);   
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");
337    }
338    canvas->SetTitle("y vs. x for simulated and reconstructed tracks");
339    pad->cd();
340 }
341
342 //-------------------------------------------------------------------
343 void 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)
349    Int_t nevents = (Int_t)gAlice->TreeE()->GetEntries();
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
377       TFile *galiceFile = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
378       galiceFile->cd();
379       gAlice->GetEvent(event);
380       MapEvent(event);
381       if (fEmpty) {
382       // skip empty events
383          fEmpty = kFALSE;
384          continue;
385       }
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 //-------------------------------------------------------------------
426 void AliMUONRecoDisplay::ShowNextEvent(Int_t delta)
427 {
428 // overwritten from AliDisplay in order to get also mapping of next event
429    TFile *galiceFile = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
430    galiceFile->cd();
431    if (delta) {
432       gAlice->Clear();
433       Int_t currentEvent = gAlice->GetHeader()->GetEvent();
434       Int_t newEvent     = currentEvent + delta;
435       if (newEvent<0 || newEvent>(gAlice->TreeE()->GetEntries() - 1)) return;
436       Int_t nparticles = gAlice->GetEvent(newEvent);
437       cout << "Event : " << newEvent << " with " << nparticles << " particles\n";
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 /******************************************************************/     
461       MapEvent(newEvent);
462       fHighlited = -1;
463    }
464    LoadPoints();
465    fPad->cd();
466    Draw();
467    if (gROOT->GetListOfCanvases()->FindObject("xy")) XYPlot();
468 }
469 //-------------------------------------------------------------------
470 Bool_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
474    if (fEmpty) return kFALSE;
475    AliDetector *pMUON = gAlice->GetDetector("MUON");
476    Bool_t chHit[10];
477    Int_t ch;
478    for (ch=0; ch<10; ch++) chHit[ch] = kFALSE;
479    //loop hits
480    for (AliMUONHit *muonHit=(AliMUONHit*)pMUON->FirstHit(track);
481         muonHit;
482         muonHit=(AliMUONHit*)pMUON->NextHit()) {
483       ch = muonHit->Chamber() - 1;
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 //-------------------------------------------------------------------
500 void 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 //-------------------------------------------------------------------
525 void AliMUONRecoDisplay::SetDrawHits(Bool_t hits)
526 {
527 // Turns on/off Geant hits drawing
528    fDrawHits = hits;
529    DrawView(0,-90);
530 }
531
532 //-------------------------------------------------------------------
533 void 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 //-------------------------------------------------------------------
543 Int_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;
550    if (fEmpty) return -1;
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);
561       if (!gtrack) continue;
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);
574          if (dX<tolerance*0.144 && dY<tolerance*0.01) {// within tol*sigma resolution
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 //-------------------------------------------------------------------
604 void AliMUONRecoDisplay::Highlight(Int_t track)
605 {
606 // Highlight the specified track
607    if (fEmpty) return; 
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;
614 //   DrawView(15,-45,135);
615    fPad->cd();
616    Draw();
617 }
618
619 //-------------------------------------------------------------------
620 void AliMUONRecoDisplay::UnHighlight()
621 {
622 // Unhighlight a previous highlighted track
623    if (fHighlited < 0 || fEmpty) return;      // nothing to do
624    TPolyLine3D *line = (TPolyLine3D*)fPolyRecoList->UncheckedAt(fHighlited);
625    line->SetLineColor(kRed);
626    line->SetLineWidth(1);
627    fHighlited = -1;
628 //   DrawView(0,-90);
629    fPad->cd();
630    Draw();
631 }
632
633 //-------------------------------------------------------------------
634 void 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
694    if (fEmpty) return;
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 //-------------------------------------------------------------------
728 void 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
732    if (fEmpty) return;
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 //-------------------------------------------------------------------
753 TClonesArray* AliMUONRecoDisplay::MakePolyLines3D(TClonesArray *tracklist)
754 {
755 // Makes the list of polylines3D corresponding to the list of tracks
756    if (fEmpty) return 0;
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 //-------------------------------------------------------------------
803 void AliMUONRecoDisplay::PolyLineInfo(TClonesArray *line3Dlist)
804 {
805 // Prints information (x, y, z coordinates) for all constructed polylines
806    if (fEmpty) return;
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