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