MUON/AliMUONRecoDisplay.cxx

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