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