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