New package for reconstructed tracks (A. Gheata):
authorgosset <gosset@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Nov 2000 11:24:10 +0000 (11:24 +0000)
committergosset <gosset@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Nov 2000 11:24:10 +0000 (11:24 +0000)
* tree output in the file "tree_reco.root"
* display events and make histograms from this file

MUON/AliMUONEventReconstructor.cxx
MUON/AliMUONEventReconstructor.h
MUON/AliMUONRecoDisplay.cxx [new file with mode: 0644]
MUON/AliMUONRecoDisplay.h [new file with mode: 0644]
MUON/AliMUONRecoEvent.cxx [new file with mode: 0644]
MUON/AliMUONRecoEvent.h [new file with mode: 0644]
MUON/MUONLinkDef.h
MUON/MUONreco.C
MUON/MUONrecodisplay.C [new file with mode: 0644]
MUON/Makefile
MUON/ReadTree.C [new file with mode: 0644]

index 391f28f384c28863316db5e0073e7cf307940b66..fb97ac1924757d4a8edef8e93aaa56673ef24390 100644 (file)
 
 /*
 $Log$
+Revision 1.18  2000/10/26 12:47:03  gosset
+Real distance between chambers of each station taken into account
+for the reconstruction parameters "fSegmentMaxDistBending[5]"
+
 Revision 1.17  2000/10/24 09:26:20  gosset
 Comments updated
 
@@ -131,6 +135,7 @@ Addition of files for track reconstruction in C++
 #include "AliMagF.h"
 #include "AliRun.h"
 #include "TParticle.h"
+#include "AliMUONRecoEvent.h"
 
 //************* Defaults parameters for reconstruction
 static const Double_t kDefaultMinBendingMomentum = 3.0;
@@ -192,6 +197,12 @@ AliMUONEventReconstructor::AliMUONEventReconstructor(void)
     gAlice->Field()->Dump();
     cout << endl;
   }
+  
+  // Initialize to 0 pointers to RecoEvent, tree and tree file
+  fRecoEvent = 0;
+  fEventTree = 0;
+  fTreeFile  = 0;
+  
   return;
 }
 
@@ -210,6 +221,12 @@ AliMUONEventReconstructor & AliMUONEventReconstructor::operator=(const AliMUONEv
 AliMUONEventReconstructor::~AliMUONEventReconstructor(void)
 {
   // Destructor for class AliMUONEventReconstructor
+  if (fTreeFile) {
+     fTreeFile->Close();
+     delete fTreeFile;
+  }
+//  if (fEventTree) delete fEventTree;
+  if (fRecoEvent) delete fRecoEvent;
   delete fHitsForRecPtr; // Correct destruction of everything ???? or delete [] ????
   for (Int_t st = 0; st < kMaxMuonTrackingStations; st++)
     delete fSegmentsPtr[st]; // Correct destruction of everything ????
@@ -1355,3 +1372,30 @@ void AliMUONEventReconstructor::EventDump(void)
   return;
 }
 
+void AliMUONEventReconstructor::FillEvent()
+{
+// Create a new AliMUONRecoEvent, fill its track list, then add it as a
+// leaf in the Event branch of TreeRecoEvent tree
+   cout << "Enter FillEvent() ...\n";
+
+   if (!fRecoEvent) {
+      fRecoEvent = new AliMUONRecoEvent();
+   } else {
+      fRecoEvent->Clear();
+   }
+   //save current directory
+   TDirectory *current =  gDirectory;
+   if (!fTreeFile)  fTreeFile  = new TFile("tree_reco.root", "RECREATE");
+   if (!fEventTree) fEventTree = new TTree("TreeRecoEvent", "MUON reconstructed events");
+   if (fRecoEvent->MakeDumpTracks(fRecTracksPtr)) {
+      if (fPrintLevel > 1) fRecoEvent->EventInfo();
+      TBranch *branch = fEventTree->GetBranch("Event");
+      if (!branch) branch = fEventTree->Branch("Event", "AliMUONRecoEvent", &fRecoEvent, 64000,1);
+      branch->SetAutoDelete();
+      fTreeFile->cd();
+      fEventTree->Fill();
+      fTreeFile->Write();
+   }
+   // restore directory
+   current->cd();
+}
index 6af4408fa9fe5ed65b4f0f1ff58ab6afa4d19895..a1d8f96ce015ab81355cc0ddd56c927a15917bb0 100644 (file)
@@ -13,6 +13,7 @@ class AliMUONSegment;
 class TClonesArray;
 class TFile;
 class TTree;
+class AliMUONRecoEvent;
 
 // Constants which should be elsewhere ????
 const Int_t kMaxMuonTrackingChambers = 10;
@@ -75,6 +76,7 @@ class AliMUONEventReconstructor : public TObject {
   Double_t GetBendingMomentumFromImpactParam(Double_t ImpactParam);
   void EventReconstruct(void);
   void EventDump(void);  // dump reconstructed event
+  void FillEvent();      // fill and write tree of reconstructed events
 
  protected:
 
@@ -126,6 +128,11 @@ class AliMUONEventReconstructor : public TObject {
   TClonesArray *fRecTrackHitsPtr; // pointer to array of hits on reconstructed tracks
   Int_t fNRecTrackHits; // number of hits on reconstructed tracks
 
+  // Objects needed for tree writing
+  AliMUONRecoEvent *fRecoEvent; // the reconstructed event
+  TTree            *fEventTree; // tree of reconstructed events
+  TFile            *fTreeFile;  // file where the tree is outputed
+
   // Functions
   void ResetHitsForRec(void);
   void MakeEventToBeReconstructed(void);
diff --git a/MUON/AliMUONRecoDisplay.cxx b/MUON/AliMUONRecoDisplay.cxx
new file mode 100644 (file)
index 0000000..e390183
--- /dev/null
@@ -0,0 +1,678 @@
+//Authors: Mihaela Gheata, Andrei Gheata 09/10/00
+
+#include <iostream.h>
+#include <AliRun.h>
+#include <TClonesArray.h>
+#include "AliMUONRecoEvent.h"
+#include "AliMUONRecoDisplay.h"
+#include <TROOT.h>
+#include <AliPoints.h>
+#include <TSlider.h>
+#include <TView.h>
+#include <TGeometry.h>
+
+ClassImp(AliMUONRecoDisplay)
+
+//-------------------------------------------------------------------
+AliMUONRecoDisplay::AliMUONRecoDisplay(Int_t nevent)
+                  :AliDisplay(750)
+{
+//************ Constructor of the reco. event display**********
+   // get reconstructed event from file
+   fFile = new TFile("tree_reco.root");
+   if (!fFile) {
+      cout << "File tree_reco.root not found\n";
+      gApplication->Terminate(0);
+   }
+   fEvReco = 0;
+   fTree = (TTree *) fFile->Get("TreeRecoEvent");
+   if (!fTree) {
+      cout << "Tree of reconstructed events not found on file. Abort.\n";
+      gApplication->Terminate(0);
+   }
+   if (nevent > fTree->GetEntries()) {
+      cout << "Event number out of range. Aborting.\n";
+      gApplication->Terminate(0);
+   }
+   TBranch *branch = fTree->GetBranch("Event");
+   branch->SetAddress(&fEvReco);
+   
+   fEvGen  = new AliMUONRecoEvent();
+
+   TFile *galice_file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+   galice_file->cd();
+   gAlice->GetEvent(nevent);
+
+
+   fRecoTracks = 0;
+   fGenTracks  = 0;
+   fPolyRecoList = 0;
+   fPolyGenList  = 0;
+   fHighlited   = -1;
+   fMinMomentum = 0;
+   fMaxMomentum = 999;
+   MapEvent(nevent);
+}
+
+//-------------------------------------------------------------------
+AliMUONRecoDisplay::~AliMUONRecoDisplay()
+{
+// Destructor of display object
+   if (fPolyRecoList) {
+      fPolyRecoList->Delete();
+      delete fPolyRecoList;
+   }
+   if (fPolyGenList) {
+      fPolyGenList->Delete();
+      delete fPolyGenList;
+   }
+   delete fEvGen;
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::MapEvent(Int_t nevent)
+{
+// get generated event (nevent) from galice.root and corresponding
+// reconstructed event from tree_reco.root; 
+   cout << "mapping event " << nevent << endl;
+   fFile->cd();
+   fTree->GetEntry(nevent);
+   // just testing
+   if (!fEvReco) {
+      cout << "failed to load reco event ! Correct this.\n";
+      gApplication->Terminate(0);
+   }
+   fRecoTracks   = fEvReco->TracksPtr();
+   fPolyRecoList = MakePolyLines3D(fRecoTracks);   
+   // clear previous event
+   if (fEvGen) fEvGen->Clear();
+   fEvGen->SetNoEvent(nevent);
+   // get list of particles
+   TClonesArray *Particles = gAlice->Particles();
+   // connect MUON module
+   AliDetector *MUON = gAlice->GetDetector("MUON");
+   if (!MUON) {
+      cout << "MUON module not present.\n";
+      gApplication->Terminate(0);
+   }
+   // get the number of generated tracks
+   Int_t ntracks = (Int_t)gAlice->TreeH()->GetEntries();
+   // Fill the fEvGen object
+   AliMUONRecoTrack *gtrack = 0;
+   AliMUONHit *hit = 0;
+   TParticle *particle;
+   Int_t ch;
+   // loop all tracks
+   for (Int_t track=0; track<ntracks; track++) {
+      hit = (AliMUONHit *) MUON->FirstHit(track);
+      if (!hit) continue;
+      particle = (TParticle *) Particles->UncheckedAt(hit->Track());
+      if (IsReconstructible(track) && TMath::Abs(particle->GetPdgCode())==13) {
+         gtrack = fEvGen->AddEmptyTrack();
+        gtrack->SetSign(TMath::Sign((Int_t)1, -particle->GetPdgCode()));
+        // reset hits
+        for (ch=0; ch<10; ch++) gtrack->SetHitPosition(ch,0,0,0);
+        // loop all hits
+        for (AliMUONHit *muonHit=(AliMUONHit*)MUON->FirstHit(track);
+             muonHit;
+             muonHit=(AliMUONHit*)MUON->NextHit()) {
+           ch = muonHit->fChamber - 1;
+           if (ch<0 || ch>9) continue;
+           gtrack->SetHitPosition(ch, muonHit->X(),  muonHit->Y(),  muonHit->Z());
+           gtrack->SetMomReconstr(particle->Px(), particle->Py(), particle->Pz());
+           gtrack->SetVertexPos(particle->Vz());
+           gtrack->SetChi2r(0.0);
+        } 
+      }
+   }
+   fGenTracks   = fEvGen->TracksPtr();
+   fPolyGenList = MakePolyLines3D(fGenTracks);
+}
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::XYPlot()
+{
+// Plot reco. tracks hits in all chambers for current event:
+//  - open blue squares : generated muons that were reconstructed accurate
+//  - open cyan squares : generated muons that were not reconstructed
+//  - filled green circles : reco. tracks (accurate)
+//  - filled red squares   : fake tracks 
+   Double_t kMaxRadius[10];
+   kMaxRadius[0] =  kMaxRadius[1] = 91.5;
+   kMaxRadius[2] =  kMaxRadius[3] = 122.5;
+   kMaxRadius[4] =  kMaxRadius[5] = 158.3;
+   kMaxRadius[6] =  kMaxRadius[7] = 260.0;
+   kMaxRadius[8] =  kMaxRadius[9] = 260.0;
+
+   TH2F *xygen_found[10]; 
+   TH2F *xygen_lost[10]; 
+   TH2F *xyreco_good[10];
+   TH2F *xyreco_fake[10];
+   Double_t x,y,r;
+   Int_t matches[500];
+   Int_t index, ch;
+
+   TPad *pad = (TPad*)gROOT->GetSelectedPad();
+   TCanvas *canvas = new TCanvas("xy", "Reconstruction efficiency");
+   canvas->Clear();
+   canvas->cd();
+   canvas->Divide(3,4);
+   canvas->cd(11);
+   gPad->Delete();
+   canvas->cd(12);
+   gPad->Delete();
+   canvas->cd(1);
+   
+   // Define histograms for x-y plots
+   for (ch=0; ch<10; ch++) {
+      xygen_found[ch] = new TH2F("xygen_found","",50,-kMaxRadius[ch],kMaxRadius[ch],50,-kMaxRadius[ch],kMaxRadius[ch]);
+      xygen_lost[ch] = new TH2F("xygen_lost","",50,-kMaxRadius[ch],kMaxRadius[ch],50,-kMaxRadius[ch],kMaxRadius[ch]);
+      xyreco_good[ch] = new TH2F("xyreco_good","",50,-kMaxRadius[ch],kMaxRadius[ch],50,-kMaxRadius[ch],kMaxRadius[ch]);
+      xyreco_fake[ch] = new TH2F("xyreco_fake","",50,-kMaxRadius[ch],kMaxRadius[ch],50,-kMaxRadius[ch],kMaxRadius[ch]);
+   }
+   // find list of matching tracks
+   fPrinted = kTRUE; // no need to print
+   for (index=0; index<fRecoTracks->GetEntriesFast(); index++) {
+      matches[index] = GetBestMatch(index);  
+   }
+   // and fill them
+   for (index=0; index<fGenTracks->GetEntriesFast(); index++) {        // GEANT
+      Bool_t wasreconst = kFALSE;
+      for (Int_t i=0; i<fRecoTracks->GetEntriesFast(); i++) {
+         if (matches[i] == index) {
+           wasreconst = kTRUE;
+           break;
+        }
+      }
+      AliMUONRecoTrack *current = (AliMUONRecoTrack*)fGenTracks->UncheckedAt(index);
+      for (ch=0; ch<10; ch++) {
+         x = current->GetPosX(ch);
+         y = current->GetPosY(ch);
+         r = TMath::Sqrt(x*x +y*y);
+         if (r >= 10) {
+            if (wasreconst) {
+               xygen_found[ch]->Fill(x,y);
+            } else {xygen_lost[ch]->Fill(x,y);}
+         }
+      }
+   }
+   for (index=0; index<fRecoTracks->GetEntriesFast(); index++) {       // reco
+      AliMUONRecoTrack *current = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(index);
+      for (ch=0; ch<10; ch++) {
+         x = current->GetPosX(ch);
+         y = current->GetPosY(ch);
+         r = TMath::Sqrt(x*x +y*y);
+         if (r >= 10) {
+            if (matches[index] >= 0) {
+               xyreco_good[ch]->Fill(x,y);
+            } else {xyreco_fake[ch]->Fill(x,y);}
+         }
+      }
+   }
+   // finally plot them
+   for (ch=0; ch<10; ch++) {
+      canvas->cd(ch+1);   
+      xygen_found[ch]->SetMarkerColor(kBlue);
+      xygen_found[ch]->SetMarkerStyle(4);
+      xygen_found[ch]->SetMarkerSize(0.5);
+      xygen_found[ch]->SetStats(kFALSE);
+      xygen_found[ch]->Draw();
+      xygen_lost[ch]->SetMarkerColor(kCyan);
+      xygen_lost[ch]->SetMarkerStyle(4);
+      xygen_lost[ch]->SetMarkerSize(0.5);
+      xygen_lost[ch]->SetStats(kFALSE);
+      xygen_lost[ch]->Draw("SAME");
+      xyreco_good[ch]->SetMarkerColor(kGreen);
+      xyreco_good[ch]->SetMarkerStyle(20);
+      xyreco_good[ch]->SetMarkerSize(0.4);
+      xyreco_good[ch]->SetStats(kFALSE);
+      xyreco_good[ch]->Draw("SAME");
+      xyreco_fake[ch]->SetMarkerColor(kRed);
+      xyreco_fake[ch]->SetMarkerStyle(20);
+      xyreco_fake[ch]->SetMarkerSize(0.5);
+      xyreco_fake[ch]->SetStats(kFALSE);
+      xyreco_fake[ch]->Draw("SAME");
+   }
+   canvas->SetTitle("y vs. x for simulated and reconstructed tracks");
+   pad->cd();
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::RecoEfficiency(Int_t first, Int_t last)
+{
+// Loop selected reconstructed events, compute efficiency as total number of 
+// reconstructed tracks (accurate) over total number of generated
+// reconstructible muons. Make histogram for momentum precision and profile
+// of mean momentum precision versus total momentum (generated)
+   Int_t nevents = (Int_t)fTree->GetEntries();
+   if (last > nevents) last = nevents - 1;
+   if (first < 0) first = 0;
+   nevents = last - first + 1;
+   Int_t track;
+   Float_t efficiency;
+   Float_t generated=0, found=0, fake=0; // number of generated/found/fake tracks
+   Double_t pgen, preco, dP;             // dP=preco-pgen
+   AliMUONRecoTrack *rtrack, *gtrack;    // generated/reco. current tracks
+   fPrinted = kTRUE;                     // no need to print
+   Int_t currentEvent = gAlice->GetHeader()->GetEvent();
+   
+   TH1F *pReso = new TH1F("Momentum precision", "dP = Prec - Pgen", 100, -5.0, 5.0);
+   pReso->SetXTitle("dP [GeV/c]");
+   pReso->SetYTitle("dN/dP");
+   
+   TProfile *dPp = new TProfile("", "dP vs. P", 50, 0, 100, 0, 5);
+   dPp->SetXTitle("P [GeV/c]");
+   dPp->SetYTitle("<dP> [GeV/c]"); 
+
+   TPad *pad = (TPad*)gROOT->GetSelectedPad();
+   TCanvas *canvas = new TCanvas("c", "Reconstruction efficiency");
+   canvas->Divide(1,2);
+   canvas->cd(1);
+
+   // loop events
+   for (Int_t event=first; event<=last; event++) {
+      // get the reco. & gen. events
+      TFile *galice_file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+      galice_file->cd();
+      gAlice->GetEvent(event);
+      MapEvent(event);
+      generated += fGenTracks->GetEntriesFast();
+      // loop reco. tracks
+      for (track=0; track<fRecoTracks->GetEntriesFast(); track++) {
+         rtrack = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(track);
+         preco  = rtrack->P();
+        // find correspondent generated track
+        Int_t ind = GetBestMatch(track);
+        if (ind >=0) {
+           gtrack = (AliMUONRecoTrack*)fGenTracks->UncheckedAt(ind);
+           pgen  = gtrack->P();
+           found += 1;
+           dP = preco - pgen;
+           pReso->Fill(dP);
+           dPp->Fill(pgen, TMath::Abs(dP));
+        } else {
+           fake += 1;
+        }
+      }
+   }
+   efficiency = found/generated;
+   cout << "=================================================================\n";
+   cout << "||       Reconstruction efficiency and momentum precision      ||\n";
+   cout << "=================================================================\n";
+   cout << "Number of events processed              " << nevents << endl;
+   cout << "Total number of reconstructible muons : " << (Int_t)generated << endl;
+   cout << "RECONSTRUCTION EFFICIENCY :             " << efficiency*100 << " %" << endl;
+   cout << "Fake track rate :                       " << 100*fake/generated << " %" << endl;
+   cout << "Momentum precision fit : \n" << endl;
+   pReso->Fit("gaus");
+   cout << "=================================================================\n";    
+   
+   canvas->cd(1);
+   pReso->Draw();
+   canvas->cd(2);
+   dPp->Draw("HIST");
+   pad->cd();
+   ShowNextEvent(currentEvent-last);
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::ShowNextEvent(Int_t delta)
+{
+// overwritten from AliDisplay in order to get also mapping of next event
+   TFile *galice_file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+   galice_file->cd();
+   if (delta) {
+      gAlice->Clear();
+      Int_t currentEvent = gAlice->GetHeader()->GetEvent();
+      Int_t newEvent     = currentEvent + delta;
+      Int_t nparticles = gAlice->GetEvent(newEvent);
+      cout << "Event : " << newEvent << " with " << nparticles << " particles\n";
+      if (!gAlice->TreeH()) return;
+      MapEvent(newEvent);
+      fHighlited = -1;
+   }
+   LoadPoints();
+   fPad->cd();
+   Draw();
+   if (gROOT->GetListOfCanvases()->FindObject("xy")) XYPlot();
+}
+//-------------------------------------------------------------------
+Bool_t AliMUONRecoDisplay::IsReconstructible(Int_t track)
+{
+// true if at least three hits in first 2 stations, 3 in last 2 stations
+// and one in station 3
+   AliDetector *MUON = gAlice->GetDetector("MUON");
+   Bool_t chHit[10];
+   Int_t ch;
+   for (ch=0; ch<10; ch++) chHit[ch] = kFALSE;
+   //loop hits
+   for (AliMUONHit *muonHit=(AliMUONHit*)MUON->FirstHit(track);
+        muonHit;
+       muonHit=(AliMUONHit*)MUON->NextHit()) {
+      ch = muonHit->fChamber - 1;
+      if (ch<0 || ch>9) continue;
+      chHit[ch] = kTRUE;
+   }
+   Int_t nhits = 0;
+   for (ch=0; ch<4; ch++) nhits += (chHit[ch])?1:0;
+   if (nhits < 3) return kFALSE;
+   nhits = 0;
+   for (ch=4; ch<6; ch++) nhits+= (chHit[ch])?1:0;
+   if (nhits < 1) return kFALSE;
+   
+   for (ch=7; ch<10; ch++) nhits+= (chHit[ch])?1:0;
+   if (nhits < 3) return kFALSE;
+   return kTRUE;
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::DrawView(Float_t theta, Float_t phi, Float_t psi)
+{
+// ovewritten from base class to change the range for MUON
+   gPad->SetCursor(kWatch);
+   gPad->SetFillColor(1);
+   gPad->Clear();
+   
+   Int_t iret;
+   TView *view = new TView(1);
+   Float_t range = fRrange*fRangeSlider->GetMaximum();
+   view->SetRange(-range, -range, 0, range, range, 5*range);
+   fZoomX0[0] = -1;
+   fZoomY0[0] = -1;
+   fZoomX1[0] =  1;
+   fZoomY1[0] =  1;
+   // Display Alice geometry
+   gAlice->GetGeometry()->Draw("same");
+   //Loop on all detectors to add their products to the pad
+   DrawHits();
+   // add itself to the list (last)
+   AppendPad();
+   
+   view->SetView(phi, theta, psi, iret);
+}
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::SetDrawHits(Bool_t hits)
+{
+// Turns on/off Geant hits drawing
+   fDrawHits = hits;
+   DrawView(0,-90);
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::CutMomentum(Double_t min, Double_t max)
+{
+// Define momentum cut for displayed reconstructed tracks
+   fMinMomentum = min;
+   fMaxMomentum = max;
+   if (fHighlited >= 0) UnHighlight();
+   DrawView(0,-90);
+}
+
+//-------------------------------------------------------------------
+Int_t AliMUONRecoDisplay::GetBestMatch(Int_t indr, Float_t tolerance)
+{
+// Find the index of best Geant track matching a reconstructed track : track
+//     with maximum number of compatible hits (within tolerance*sigma bending and
+//     non-bending resolution) and minimum number of fake hits;
+// If no match is found within a given tolerance, the method is called recursively
+//      with increasing tolerance, until tolerance = 10;
+   if (indr<0 || indr>=fRecoTracks->GetEntriesFast()) return -1;
+   AliMUONRecoTrack *rtrack = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(indr);
+   AliMUONRecoTrack *gtrack = 0;
+   Int_t bestMatch = -1;
+   Int_t maxNcompat = 0;
+   Int_t minNfakes = 10;
+   Double_t xrhit,yrhit,radius,xghit,yghit,dX,dY;
+// loop over all Geant tracks
+   for (Int_t indg=0; indg<fGenTracks->GetEntriesFast(); indg++) {
+      gtrack = (AliMUONRecoTrack*)fGenTracks->UncheckedAt(indg);
+       if (!gtrack) continue;
+      Int_t ncompat = 0;      // number of compat. hits for this track
+      Int_t nfake = 0;        // number of fakes
+      // loop chambers to find compatible hits
+      for (Int_t ch=0; ch<10; ch++) {
+         xrhit = rtrack->GetPosX(ch);
+         yrhit = rtrack->GetPosY(ch);
+         radius = TMath::Sqrt(xrhit*xrhit + yrhit*yrhit);
+         if (radius<10) continue; // skip null hits
+         xghit = gtrack->GetPosX(ch);
+         yghit = gtrack->GetPosY(ch);
+         dX = TMath::Abs(xghit-xrhit);
+         dY = TMath::Abs(yghit-yrhit);
+         if (dX<tolerance*0.144 && dY<tolerance*0.01) {      // within 3 sigma resolution
+            ncompat++;
+            continue;      // compatible hit
+         } else nfake++;      // fake hit
+      }
+      if (ncompat && ncompat>=maxNcompat && nfake<minNfakes) { // this is best matching
+         maxNcompat = ncompat;
+         minNfakes = nfake;
+         bestMatch = indg;
+      }
+   }
+   if (bestMatch<0 && tolerance<=9.) bestMatch = GetBestMatch(indr, tolerance+=1);
+   if (!fPrinted) {
+      rtrack = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(indr);
+      Int_t sign = rtrack->GetSign();
+      cout << "Reconstructed track : " << indr << "(" << sign << ")" << endl;
+      rtrack->TrackInfo();
+      printf("Best matching Geant track within %i*sgm : %i\n", (Int_t)tolerance, bestMatch);
+      if (bestMatch >=0) {
+         gtrack = (AliMUONRecoTrack*)fGenTracks->UncheckedAt(bestMatch);
+         gtrack->TrackInfo();
+      }
+      cout << "-----------------------------------------------------------------\n";
+   }
+   fPrinted = kTRUE;
+
+   return bestMatch;
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::Highlight(Int_t track)
+{
+// Highlight the specified track 
+   if (fHighlited >=0) UnHighlight();
+   if (track<0 || track>fPolyRecoList->GetEntries()) return;
+   TPolyLine3D *line = (TPolyLine3D*)fPolyRecoList->UncheckedAt(track);
+   line->SetLineColor(kYellow);
+   line->SetLineWidth(1);
+   fHighlited = track;
+   DrawView(15,-45,135);
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::UnHighlight()
+{
+// Unhighlight a previous highlighted track
+   if (fHighlited < 0) return;      // nothing to do
+   TPolyLine3D *line = (TPolyLine3D*)fPolyRecoList->UncheckedAt(fHighlited);
+   line->SetLineColor(kRed);
+   line->SetLineWidth(1);
+   fHighlited = -1;
+   DrawView(0,-90);
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::DrawHits()
+{
+//    Draw hits for all ALICE detectors. Overwrites the DrawHits() method of the
+//     base class for reco. track drawing
+
+   Float_t cutmin, cutmax, etamin, etamax, pmom, smin, smax, eta, theta, r;
+   const Float_t kptcutmax = 2;
+   const Float_t ketacutmax = 1.5;
+   Float_t *pxyz;
+   Int_t ntracks,track;
+   TParticle *particle;
+   TObjArray *points;
+   AliPoints *pm;
+      
+   //Get cut slider
+   smax   = fCutSlider->GetMaximum();
+   smin   = fCutSlider->GetMinimum();
+   cutmin = kptcutmax*smin;
+   if (smax < 0.98) cutmax = kptcutmax*smax;
+   else             cutmax = 100000;
+   
+   //Get eta slider
+   smax   = fEtaSlider->GetMaximum();
+   smin   = fEtaSlider->GetMinimum();
+   etamin = ketacutmax*(2*smin-1);
+   etamax = ketacutmax*(2*smax-1);
+   if (smin < 0.02) etamin = -1000;
+   if (smax > 0.98) etamax =  1000;
+      
+   TIter next(gAlice->Modules());
+   AliModule *module;
+   fHitsCuts = 0;
+   if (fDrawHits) {
+      // draw hits in all modules
+      while((module = (AliModule*)next())) {
+         if (!module->IsActive()) continue;
+         points = module->Points();
+         if (!points) continue;
+         ntracks = points->GetEntriesFast();
+         for (track=0;track<ntracks;track++) {
+            pm = (AliPoints*)points->UncheckedAt(track);
+            if (!pm) continue;
+            particle = pm->GetParticle();
+            if (!particle) continue;
+            pmom = particle->P();
+            if (pmom < cutmin) continue;
+            if (pmom > cutmax) continue;
+            // as a first approximation, take eta of first point
+            pxyz  = pm->GetP();
+            r     = TMath::Sqrt(pxyz[0]*pxyz[0] + pxyz[1]*pxyz[1]);
+            theta = TMath::ATan2(r,TMath::Abs(pxyz[2]));
+            if(theta) eta = -TMath::Log(TMath::Tan(0.5*theta)); else eta = 1e10;
+            if (pxyz[2] < 0) eta = -eta;
+            if (eta < etamin || eta > etamax) continue;
+            pm->Draw();
+            fHitsCuts += pm->GetN();
+         }
+      }
+   }
+   // draw reconstructed tracks
+   TPolyLine3D *line, *gline;
+   Int_t bestMatch;
+   Double_t px,py,pz,p;
+   AliMUONRecoTrack *rtrack;
+
+   if (fHighlited >= 0) {
+      line = (TPolyLine3D*)fPolyRecoList->UncheckedAt(fHighlited);
+      line->Draw();
+      fPrinted = kFALSE;
+      bestMatch = GetBestMatch(fHighlited);
+      if (bestMatch>=0) {
+         gline = (TPolyLine3D*)fPolyGenList->UncheckedAt(bestMatch);
+         gline->SetLineColor(kRed);
+         gline->SetLineWidth(2);
+         gline->SetLineStyle(2);
+         gline->Draw();
+      }
+   } else {
+      for (track=0; track<fPolyRecoList->GetEntries(); track++) {
+         rtrack = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(track);
+         px = rtrack->GetMomReconstr(0);
+         py = rtrack->GetMomReconstr(1);
+         pz = rtrack->GetMomReconstr(2);
+         p  = rtrack->P();
+         if (p>fMinMomentum && p<fMaxMomentum) {
+            line = (TPolyLine3D*)fPolyRecoList->UncheckedAt(track);
+            line->Draw();
+         }
+      }
+   }
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::ListTracks()
+{
+// List momentum information of all reconstructed traccks within fPmin and fPmax
+//     cuts, as well as their best matching Geant tracks
+   cout << "================================================================\n";
+   printf("Reconstructed tracks with momentum in range : %g , %g [GeV/c]\n",
+         fMinMomentum, fMaxMomentum);
+   cout << "----------------------------------------------------------------\n";
+   AliMUONRecoTrack *rtrack;
+   Double_t p;
+   Int_t sign;
+   for (Int_t ind=0; ind<fRecoTracks->GetEntries(); ind++) {
+      rtrack = (AliMUONRecoTrack*)fRecoTracks->UncheckedAt(ind);
+      p  = rtrack->P();
+      if (p>fMinMomentum && p<fMaxMomentum) {
+         fPrinted = kFALSE;
+         GetBestMatch(ind);
+         sign = rtrack->GetSign();
+      }
+   }      
+   cout << "================================================================\n";
+}
+
+//-------------------------------------------------------------------
+TClonesArray* AliMUONRecoDisplay::MakePolyLines3D(TClonesArray *tracklist)
+{
+// Makes the list of polylines3D corresponding to the list of tracks
+   if (tracklist!=fRecoTracks && tracklist!=fGenTracks) return 0;
+   Bool_t reco = (tracklist==fRecoTracks)?kTRUE:kFALSE;
+   // make sure there is no other list in memory
+   if (reco) {
+      if (fPolyRecoList) {
+         fPolyRecoList->Delete();
+         delete fPolyRecoList;
+         fPolyRecoList = 0;
+      }
+   } else {
+      if (fPolyGenList) {
+         fPolyGenList->Delete();
+         delete fPolyGenList;
+         fPolyGenList = 0;
+      }      
+   }
+   if (!tracklist->GetEntries()) return 0;
+   
+   AliMUONRecoTrack* track = 0;
+   TClonesArray *polyLines3D = new TClonesArray("TPolyLine3D",1000);
+   TClonesArray &polylist = *polyLines3D;
+   TPolyLine3D *polyline = 0;
+   
+   // loop all tracks
+   for (Int_t i=0; i<tracklist->GetEntries(); i++) {
+      track = (AliMUONRecoTrack*)tracklist->UncheckedAt(i);
+      Int_t ch = 0;
+      Double_t x,y,z,r;
+      polyline = new(polylist[i]) TPolyLine3D(2,"");
+      polyline->SetLineColor(kRed);
+      polyline->SetLineWidth(1);
+      polyline->SetNextPoint(0,0,track->GetVertexPos()); // vertex point
+      // loop chambers to fill TPolyLine3D objects
+      for (ch=0; ch<10; ch++) {
+         x = track->GetPosX(ch);
+         y = track->GetPosY(ch);
+         r = TMath::Sqrt(x*x + y*y);
+         if (r < 10) continue;
+         z = track->GetPosZ(ch);
+         polyline->SetNextPoint(x,y,z);
+      }      
+   }
+   return polyLines3D;
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoDisplay::PolyLineInfo(TClonesArray *line3Dlist)
+{
+// Prints information (x, y, z coordinates) for all constructed polylines
+   if (line3Dlist) {
+      TPolyLine3D *polyline = 0;
+      for(Int_t trackIndex=0; trackIndex<line3Dlist->GetEntries(); trackIndex++) {
+         polyline = (TPolyLine3D*)line3Dlist->UncheckedAt(trackIndex);
+         polyline->ls();
+         Float_t *pl = polyline->GetP();
+         for (Int_t i=0; i<polyline->GetN() ;i++) {
+            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]);
+         }
+      }
+   }
+}
+
+
diff --git a/MUON/AliMUONRecoDisplay.h b/MUON/AliMUONRecoDisplay.h
new file mode 100644 (file)
index 0000000..aa36e4a
--- /dev/null
@@ -0,0 +1,88 @@
+// Authors : M.Gheata, A.Gheata 09/10/00
+#ifndef MUON_RECDISPLAY
+#define MUON_RECDISPLAY
+
+//////////////////////////////////////////////////////////////////////
+//                                                                  //
+// AliMUONRecoDisplay                                              //
+//                                                                 //
+// This class subclasses AliDisplay and provides display of         //
+// reconstructed tracks with following functionality :                     //
+//     - front/top/side/3D display of MUON reconstructed tracks    //
+//        as polylines ;                                            //
+//     - context menu activated when the main pad is right-clicked //
+//     The context menu contains following functions :             //
+//     * SetDrawHits() - switches on or off Geant hits ;           //
+//     * CutMomentum() - displays only tracks within Pmin - Pmax   //
+//     * ListTracks()  - prints ID and momentum info. for all      //
+//     tracks within momentum range Pmin,Pmax ;                    //
+//     * Highlight()   - shows only one selected reco. track       //
+//     and its best matching Geant track;                          //
+//     * UnHighlight() - self explaining;                          //
+//     * RecoEfficiency() - compute reco. efficiency for all events//
+//        from galice.root file; also fake track percentage; make   //
+//        plots for momentum precision                              //
+//      * XYPlot()      - make X-Y plots of reconstructed and       //
+//        generated tracks in all chambers                          //
+//                                                                 //
+//      Starting : generate and reconstruct events, then use the    //
+//                 MUONrecodisplay.C macro                          //
+//                                                                  //
+//////////////////////////////////////////////////////////////////////
+
+#include <TApplication.h>
+#include <TROOT.h>
+#include <TFile.h>
+#include <TPolyLine3D.h>
+#include <TParticle.h>
+#include <AliDisplay.h>
+#include <TTree.h>
+#include <TH1.h>
+#include <TH2.h>
+#include <TCanvas.h>
+#include <TProfile.h>
+#include <AliDetector.h>
+#include "AliMUONHit.h"
+
+
+class AliMUONRecoDisplay:public AliDisplay {
+
+private:
+//methods
+   Int_t              GetBestMatch(Int_t indr, Float_t tolerance=3.0);
+   TClonesArray*      MakePolyLines3D(TClonesArray *tracklist);
+   void               MapEvent(Int_t nevent);
+   Bool_t             IsReconstructible(Int_t track);
+//data members
+   AliMUONRecoEvent  *fEvGen;                   // Geant event
+   AliMUONRecoEvent  *fEvReco;                  // reconstructed event
+   TFile             *fFile;                    // file with reco. event tree
+   TTree             *fTree;                    // tree with reco. events
+   TClonesArray      *fPolyRecoList;            // list of TPolyLine3D's for reco. tracks
+   TClonesArray      *fPolyGenList;             // list of TPolyLine3D's for generated tracks
+   TClonesArray      *fRecoTracks;              // list of reco tracks
+   TClonesArray      *fGenTracks;               // list of GEANT tracks
+   Int_t              fHighlited;               // index of current highlited track
+   Double_t           fMinMomentum;             // min. cut of momentum
+   Double_t           fMaxMomentum;             // max. cut of momentum
+   Bool_t             fPrinted;                        // tracks info switch
+
+public:
+   AliMUONRecoDisplay(Int_t nevent=0);
+   virtual             ~AliMUONRecoDisplay();
+   virtual void       DrawHits();
+   virtual void       DrawView(Float_t theta, Float_t phi, Float_t psi = 0);
+   virtual void       SetDrawHits(Bool_t hits = kTRUE);        // *MENU*
+   virtual void       ShowNextEvent(Int_t delta = 1);
+   void               ListTracks();                            // *MENU*
+   void               Highlight(Int_t track=0);                // *MENU*
+   void               UnHighlight();                           // *MENU*
+   void               CutMomentum(Double_t min=0, Double_t max=999);   // *MENU*
+   void               PolyLineInfo(TClonesArray *line3Dlist);
+   void               RecoEfficiency(Int_t first=0, Int_t last=10000);  // *MENU*
+   void               XYPlot();                                 // *MENU*
+   
+   ClassDef(AliMUONRecoDisplay,0)      // MUON reco. event display
+};
+
+#endif
diff --git a/MUON/AliMUONRecoEvent.cxx b/MUON/AliMUONRecoEvent.cxx
new file mode 100644 (file)
index 0000000..92325b2
--- /dev/null
@@ -0,0 +1,236 @@
+//Authors: Mihaela Gheata, Andrei Gheata 09/10/00
+
+////////////////////////////////////////////////////////////////////
+//                                                                //
+// AliMUONRecoEvent and  AliMUONRecoTrack and AliMUONRecoDisplay  //
+//                                                                //
+// This header contains all classes needed for storing MUON       //
+// reconstructed events . The tree of reconstructed               //
+// events is constructed and filled during the FillEvent() method //
+// of AliMUONEventReconstructor class, when all reconstruction    //
+// information is available.                                      //
+//                                                                //
+////////////////////////////////////////////////////////////////////
+
+#include <iostream.h>
+#include <AliRun.h>
+#include <TClonesArray.h>
+#include "AliMUONRecoEvent.h"
+#include "AliMUONTrack.h"
+#include "AliMUONTrackParam.h"
+#include "AliMUONHitForRec.h"
+#include "AliMUONTrackHit.h"
+
+ClassImp(AliMUONRecoTrack)
+ClassImp(AliMUONRecoEvent)
+
+//-------------------------------------------------------------------
+AliMUONRecoEvent::AliMUONRecoEvent(Int_t eventNo) 
+{
+// Reconstructed event constructor
+   fTracks     = new TClonesArray("AliMUONRecoTrack",200);
+   fNevr       = eventNo;
+   fNtracks = 0;
+}
+
+//-------------------------------------------------------------------
+AliMUONRecoEvent::~AliMUONRecoEvent() 
+{
+// Destructor of AliMUONRecoEvent
+   fTracks->Delete();
+   delete fTracks;
+   fTracks = 0;
+}
+
+//-------------------------------------------------------------------
+AliMUONRecoTrack* AliMUONRecoEvent::AddEmptyTrack()
+{
+// Add a empty AliMUONRecoTrackObject to the track list
+   TClonesArray &dumptracks = *fTracks;
+   return (new(dumptracks[fNtracks++])AliMUONRecoTrack(kTRUE));
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoEvent::Clear(Option_t *option)
+{
+// Clears all track pointers from the list
+//   fTracks->Clear(option);
+   fTracks->Delete();
+   fNtracks=0;
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoEvent::EventInfo()
+{
+// Prints reconstructed event information
+   cout << "*********************Reco Dumper**************************" << endl;
+   cout << "Event number : " << fNevr << endl;
+   cout << "   Number of tracks : " << fNtracks << endl;
+   AliMUONRecoTrack *currentTrack =0;
+   Int_t trackIndex = 0;
+   for(trackIndex=0; trackIndex<fNtracks; trackIndex++) {
+      currentTrack = (AliMUONRecoTrack*)fTracks->UncheckedAt(trackIndex);
+      cout << "Track : " << trackIndex << endl;
+      cout << "   Sign : " << currentTrack->GetSign() << endl;
+      cout << "   Vertex position    : " << currentTrack->GetVertexPos() << endl;
+      Double_t momreco[3];
+      for (Int_t axis=0; axis<3; axis++) {
+         momreco[axis] = currentTrack->GetMomReconstr(axis);
+      }
+      cout << "   Reconstructed mom. : " << "Px=" << momreco[0] << "Py=" << momreco[1] << "Pz=" << momreco[2] << endl;
+      cout << "   Chi squared        : " << currentTrack->GetChi2r() << endl;
+      cout << "   Hits positions     : " << endl;
+      Double_t xhit, yhit, zhit;
+      for (Int_t chamber=0; chamber<10; chamber++) {
+         xhit = currentTrack->GetPosX(chamber);
+         yhit = currentTrack->GetPosY(chamber);
+         zhit = currentTrack->GetPosZ(chamber);
+//         cout <<"      chamber" << chamber << " X=" << xhit << " Y=" << yhit << " Z=" << zhit << endl;
+      }         
+   }
+   cout << "**********************************************************" << endl;
+}
+
+//-------------------------------------------------------------------
+Bool_t AliMUONRecoEvent::MakeDumpTracks(TClonesArray *tracksPtr)
+{
+// This method takes the pointer of the list of reconstructed tracks from
+// AliMUONEventReconstructor and fill the reconstructed AliMUONRecoEvent
+// fields.
+       cout << "Enter MakeDumpTracks..." << endl;
+   Int_t nTracks = tracksPtr->GetEntriesFast();
+   if (nTracks == 0) {
+      cout << "AliMUONRecoEvent::MakeDumpTracks: Number of tracks is zero !" << endl;
+      return kFALSE;
+   }
+   cout << tracksPtr << endl;
+   if (!tracksPtr) {
+      cout << "AliMUONRecoEvent::MakeDumpTracks() : You should call SetRecoTracksPtr() first..." << endl;
+      return kFALSE;
+   }
+       // Get event number
+   Int_t no_event = gAlice->GetHeader()->GetEvent();
+   tracksPtr->Compress();  // simple loop
+   AliMUONRecoTrack *currentTrack;
+   Int_t trackIndex = 0, nTrackHits = 0;
+   Double_t z,bendingSlope, nonBendingSlope, pYZ;
+   Double_t pX, pY, pZ;                        // reconstructed momentum components
+   Int_t isign;                                        // charge sign
+   AliMUONTrack *track = 0;
+   AliMUONTrackParam *trackParam = 0;
+   // Fill event number and number of tracks
+   fNevr = no_event;
+   // Loop over reconstructed tracks
+   for (trackIndex=0; trackIndex<nTracks; trackIndex++) {
+      currentTrack = AddEmptyTrack();
+      track = (AliMUONTrack*) ((*tracksPtr)[trackIndex]);
+      nTrackHits = track->GetNTrackHits();
+      trackParam = track->GetTrackParamAtVertex();
+      bendingSlope = trackParam->GetBendingSlope();
+      nonBendingSlope = trackParam->GetNonBendingSlope();
+      z = trackParam->GetZ();
+      pYZ = 1/TMath::Abs(trackParam->GetInverseBendingMomentum());
+      pZ = pYZ/TMath::Sqrt(1+bendingSlope*bendingSlope);
+      pX = pZ * nonBendingSlope;
+      pY = pZ * bendingSlope;
+      if (trackParam->GetInverseBendingMomentum()<0) isign=-1; else isign=1;
+      currentTrack->SetVertexPos(z);
+      currentTrack->SetMomReconstr(pX,pY,pZ);
+      currentTrack->SetSign(isign);
+//        currentTrack->SetChi2r(trackParam->GetChi2());
+      currentTrack->SetChi2r(0);
+      AliMUONTrackHit *trackHit;
+      Double_t xhit,yhit,zhit;
+         // Loop over track hits
+      for (Int_t trackHitIndex = 0; trackHitIndex < nTrackHits; trackHitIndex++) {
+         trackHit = (AliMUONTrackHit*) (*(track->GetTrackHitsPtr()))[trackHitIndex];
+         xhit = trackHit->GetHitForRecPtr()->GetNonBendingCoor();
+         yhit = trackHit->GetHitForRecPtr()->GetBendingCoor();
+         zhit = trackHit->GetHitForRecPtr()->GetZ();
+         if (trackHitIndex >= 0 && trackHitIndex < 10) {
+            currentTrack->SetHitPosition(trackHitIndex,xhit,yhit,zhit);
+         } else { cout << "track " << trackIndex << " hit out of range" << endl;} 
+      }
+   }
+   cout << "Leave MakeDumpTracks..." << endl;
+   return kTRUE;
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoEvent::Streamer(TBuffer &R__b)
+{
+// Streams an object of class AliMUONRecoEvent
+   if (R__b.IsReading()) {
+      Version_t R__v = R__b.ReadVersion(); if (R__v) { }
+      TObject::Streamer(R__b);
+      R__b >> fNevr;
+      R__b >> fNtracks;
+      fTracks->Clear();
+      fTracks->Streamer(R__b);
+   } else {
+      cout << "...writing event to file...\n";
+      R__b.WriteVersion(AliMUONRecoEvent::IsA());
+      TObject::Streamer(R__b);
+      R__b << fNevr;
+      R__b << fNtracks;
+      fTracks->Streamer(R__b);
+   }
+}
+
+//-------------------------------------------------------------------
+AliMUONRecoTrack::AliMUONRecoTrack(Bool_t active)
+{
+//Constructor of AliMUONRecoTrack
+   fSign  = 0;
+   fZvr   = 0.0;
+   fChi2r = 0.0;
+   if (active) {
+       for (Int_t axis=0; axis<3; axis++) {
+       fPr[axis] = 0.0;
+       }
+       for (Int_t chamber=0; chamber<10; chamber++) {
+       fPosX[chamber] = 0.0;
+       fPosY[chamber] = 0.0;
+       fPosZ[chamber] = 0.0;
+       }
+   }
+}
+
+//-------------------------------------------------------------------
+const Double_t AliMUONRecoTrack::Phi()
+{
+// Return trach phi angle
+       return TMath::ATan2(fPr[2], fPr[1]);
+}
+
+//-------------------------------------------------------------------
+const Double_t AliMUONRecoTrack::Theta()
+{
+// Return trach theta angle
+   return TMath::ACos(fPr[3] / P());
+}
+
+//-------------------------------------------------------------------
+void AliMUONRecoTrack::SetMomReconstr(Double_t px, Double_t py, Double_t pz)
+{
+// Set the track reconstructed momentum 
+   fPr[0] = px;
+   fPr[1] = py;
+   fPr[2] = pz;            
+} 
+   
+//-------------------------------------------------------------------          
+void AliMUONRecoTrack::SetHitPosition(Int_t chamber, Double_t x, Double_t y, Double_t z)
+{
+// Set hit coordinates in chamber[0..9]
+   fPosX[chamber] = x;
+   fPosY[chamber] = y;
+   fPosZ[chamber] = z;
+}
+//-------------------------------------------------------------------          
+void AliMUONRecoTrack::TrackInfo()
+{
+// Prints momentum info for this track
+   cout << "Px=" << GetMomReconstr(0) << " Py=" << GetMomReconstr(1) <<
+          " Pz=" << GetMomReconstr(2) << " P=" << P() << endl;
+}
diff --git a/MUON/AliMUONRecoEvent.h b/MUON/AliMUONRecoEvent.h
new file mode 100644 (file)
index 0000000..96cf2d3
--- /dev/null
@@ -0,0 +1,102 @@
+// Authors : M.Gheata, A.Gheata 09/10/00
+#ifndef MUON_RECEVENT
+#define MUON_RECEVENT
+
+////////////////////////////////////////////////////////////////////
+//                                                                //
+// AliMUONRecoEvent and AliMUONRecoTrack                          //
+//                                                                //
+// This header contains all classes needed for storing MUON       //
+// reconstructed events in a tree                                 //
+//                                                                //
+////////////////////////////////////////////////////////////////////
+
+#include <TObject.h>
+#include <TFile.h>
+#include <TParticle.h>
+#include <AliDetector.h>
+#include "AliMUONHit.h"
+
+class AliMUONRecoTrack;
+
+
+/////////////////////////////////////////////////////////////////////
+//                                                                 //
+// AliMUONRecoEvent                                                //
+//                                                                 //
+// This class handles an array of reconstructed tracks.            //
+// It provides :                                                   //
+//     - filling the tracks array according the information stored//
+//     in AliMUONEventReconstructor class ;                       //
+//     - printing event and track informations : event number,    //
+//     number of tracks, hits positions, reconstr. momentum.      //
+//                                                                 //
+/////////////////////////////////////////////////////////////////////
+
+
+class AliMUONRecoEvent:public TObject {
+
+private:
+   Int_t             fNevr;               // event number
+   Int_t             fNtracks;            // number of tracks
+   TClonesArray      *fTracks;            // list of AliMUONRecoTracks
+   
+public:
+                     AliMUONRecoEvent(Int_t eventNo = 0);
+   virtual           ~AliMUONRecoEvent();
+   AliMUONRecoTrack* AddEmptyTrack();
+   void              Clear(Option_t *option = "");
+   void              EventInfo();
+   Int_t             GetNoEvent()  const {return fNevr;}
+   Int_t             GetNoTracks() const {return fNtracks;}
+   Bool_t            MakeDumpTracks(TClonesArray *tracksPtr);
+   void              SetNoEvent(Int_t event)   {fNevr = event;}
+   void              SetNoTracks(Int_t ntracks) {fNtracks = ntracks;} 
+   TClonesArray*     TracksPtr() {return fTracks;}
+
+   ClassDef(AliMUONRecoEvent,1)        // Reconstructed event for MUON module
+};
+
+////////////////////////////////////////////////////////////////////
+//                                                                //
+// AliMUONRecoTrack                                               //
+//                                                                //
+// This class represents a reconstructed muon track.              //
+//                                                                //
+////////////////////////////////////////////////////////////////////
+
+class AliMUONRecoTrack:public TObject {
+
+private:
+    Int_t       fSign;                  // charge sign
+    Double_t   fZvr;                   // z of track vertex point
+    Double_t   fChi2r;                 // chi squared for reco. track
+    Double_t   fPr[3];                 // reconstr. momentum (same as in vertex)
+    Double_t   fPosX[10];              // hit X position in all chambers
+    Double_t   fPosY[10];              // hit Y position in all chambers    
+    Double_t   fPosZ[10];              // hit Z position in all chambers
+
+public:
+   AliMUONRecoTrack() { }
+    AliMUONRecoTrack(Bool_t active);
+    virtual        ~AliMUONRecoTrack() { }     //desctructor
+    const Double_t GetChi2r() const {return fChi2r;};
+    const Double_t GetMomReconstr(Int_t axis) const {return fPr[axis];};
+    const Int_t    GetSign() const {return fSign;};
+    const Double_t GetPosX(Int_t chamber) const {return fPosX[chamber];};
+    const Double_t GetPosY(Int_t chamber) const {return fPosY[chamber];};
+    const Double_t GetPosZ(Int_t chamber) const {return fPosZ[chamber];};
+    const Double_t GetVertexPos() { return fZvr;};
+    const Double_t P() {return TMath::Sqrt(fPr[0]*fPr[0] + fPr[1]*fPr[1] + fPr[2]*fPr[2]);};
+    const Double_t Phi();
+    void           SetChi2r(Double_t chi) { fChi2r = chi;};
+    void           SetHitPosition(Int_t chamber, Double_t x, Double_t y, Double_t z);
+    void           SetMomReconstr(Double_t px, Double_t py, Double_t pz);
+    void           SetSign(Int_t sign) {fSign = sign;};
+    void           SetVertexPos(Double_t zvr) {fZvr = zvr;};
+    const Double_t Theta();
+    void           TrackInfo();
+    ClassDef(AliMUONRecoTrack,1)       // A reconstructed muon track
+};
+
+#endif
\ No newline at end of file
index 90ca4942f4b44051a1d5dc3e59d90721db3a1287..3ab275b78719245ccce9af119e359b5d81d70f81 100644 (file)
@@ -48,6 +48,9 @@
 #pragma link C++ class  AliMUONSegmentationSlatN;
 #pragma link C++ class  AliMUONConstants;
 #pragma link C++ class  AliMUONClusterInput;
+#pragma link C++ class  AliMUONRecoEvent-;
+#pragma link C++ class  AliMUONRecoTrack;
+#pragma link C++ class  AliMUONRecoDisplay;
 #endif
 
 
index 802d15d304efb04d4407fb0940e56b3d38bc2d76..cbe7db48fc3e49770f034ee40daaf48dc09ca020 100644 (file)
@@ -65,6 +65,8 @@ void MUONreco (Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t RecGeantHits = 0
       if (Reco->GetBkgGeantFile())Reco->NextBkgGeantEvent();
     }
     Reco->EventReconstruct();
+    // Write this event in a tree
+    Reco->FillEvent();
     // Dump current event
     Reco->EventDump();
   } // Event loop
diff --git a/MUON/MUONrecodisplay.C b/MUON/MUONrecodisplay.C
new file mode 100644 (file)
index 0000000..374690d
--- /dev/null
@@ -0,0 +1,62 @@
+#include <iostream.h>
+
+void MUONrecodisplay(Int_t evNumber=0)
+{
+////////////////////////////////////////////////////////////////////
+// The display pops up a context menu when the right-button is    //
+//     clicked on the main pad. The following functions are available //
+//     * SetDrawHits() - switches on or off Geant hits ;                               //
+//     * CutMomentum() - displays only tracks within Pmin - Pmax       //
+//     * ListTracks()  - prints ID and momentum info. for all                          //
+//     tracks within momentum range Pmin,Pmax;                                                         //
+//     * Highlight()   - shows only one selected reco. track                           //
+//     and its best matching Geant track;                                                                              //
+//     *UnHighlight()  - self explaining;                                                                              //
+//     *SetDrawHits()  - switch on or off Geant hits                                                   //
+////////////////////////////////////////////////////////////////////
+   if (evNumber<0) return;
+
+
+//Dynamically link some shared libs
+   if (gClassTable->GetID("AliRun") < 0) {
+      gROOT->LoadMacro("loadlibs.C");
+       loadlibs();
+   }
+    
+   
+// Connect the Root Galice file containing ...
+    
+   TFile *galice_file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+   if (!galice_file) galice_file = new TFile("galice.root");
+   if (!galice_file) {
+      cout << "File galice.root not found\n";
+      return;
+   }
+     
+// Get AliRun object from file or create it if not on file
+   if (!gAlice) {
+       gAlice = (AliRun*)galice_file->Get("gAlice");
+       if (gAlice) {
+         cout << "AliRun object found on file\n";
+      } else {
+         cout << "AliRun object not found on file !!!\n";
+         return;
+      }
+   }
+    
+   AliMUONRecoDisplay *display = new AliMUONRecoDisplay(evNumber);   
+
+   display->DisableDetector("ITS");
+   display->DisableDetector("TPC");
+   display->DisableDetector("TOF");
+   display->DisableDetector("RICH");
+   display->DisableDetector("ZDC");
+   display->DisableDetector("CASTOR");
+   display->DisableDetector("TRD");
+   display->DisableDetector("FMD");
+   display->DisableDetector("PHOS");
+   display->DisableDetector("PMD");
+
+   display->ShowNextEvent(0);
+
+}
index c68cd23a534cae1f6c00f16eef439f084b8ba529..9a3cccd9e2e2a2e63f44e265649a1da369fc8d29 100644 (file)
@@ -34,7 +34,8 @@ SRCS         = AliMUONChamber.cxx AliMUONChamberTrigger.cxx \
                AliMUONTriggerConstants.cxx  AliMUONConstants.cxx \
               AliMUONClusterInput.cxx  \
               AliMUONSegmentationSlatModule.cxx AliMUONSegmentationSlatModuleN.cxx \
-              AliMUONSegmentationSlat.cxx AliMUONSegmentationSlatN.cxx 
+              AliMUONSegmentationSlat.cxx AliMUONSegmentationSlatN.cxx \
+              AliMUONRecoEvent.cxx AliMUONRecoDisplay.cxx
 
 # C++ Headers
 
diff --git a/MUON/ReadTree.C b/MUON/ReadTree.C
new file mode 100644 (file)
index 0000000..c1eac1d
--- /dev/null
@@ -0,0 +1,55 @@
+#include <iostream.h>
+
+
+//***************************************************************************
+
+void ReadTree()
+{
+// This is a basic example on how the tree of reconstructed events
+// should be accessed
+//Dynamically link some shared libs
+   if (gClassTable->GetID("AliRun") < 0) {
+      gROOT->LoadMacro("loadlibs.C");
+       loadlibs();
+   }
+    
+   
+// Connect the Root Galice file containing ...
+   TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("tree_reco.root");
+   if (!file) file = new TFile("tree_reco.root");
+   if (!file) {
+      cout << "File tree_reco.root not found\n";
+      return;
+   }
+    
+   TFile *galice_file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
+   if (!galice_file) galice_file = new TFile("galice.root");
+   if (!galice_file) {
+      cout << "File galice.root not found\n";
+      return;
+   }
+     
+// Get AliRun object from file or create it if not on file
+   if (!gAlice) {
+       gAlice = (AliRun*)galice_file->Get("gAlice");
+       if (gAlice) {
+         cout << "AliRun object found on file\n";
+      } else {
+         cout << "AliRun object not found on file !!!\n";
+         return;
+      }
+   }
+    
+   TTree *tree = (TTree*)file->Get("TreeRecoEvent");
+   TBranch *branch = tree->GetBranch("Event");
+   AliMUONRecoEvent *event = 0;
+   branch->SetAddress(&event);
+   Int_t nentries = (Int_t)tree->GetEntries(); // no. of reco events on file
+
+   for (Int_t evNumber=0; evNumber<nentries; evNumber++) {
+      tree->GetEntry(evNumber);
+      cout << "Event : " << event->GetNoEvent() << " with" << event->GetNoTracks() << " tracks\n";
+// print reconstr. event information
+//      event->EventInfo();
+   }
+}