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