MUON/AliMUONRecoDisplay.cxx

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