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