MUON/AliMUONRecoDisplay.cxx

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