]> git.uio.no Git - u/mrichter/AliRoot.git/commitdiff
Splitting Draw option in cluaster finder AZ (Sacha)
authormartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 4 Nov 2005 13:19:18 +0000 (13:19 +0000)
committermartinez <martinez@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 4 Nov 2005 13:19:18 +0000 (13:19 +0000)
MUON/AliMUONClusterDrawAZ.cxx [new file with mode: 0644]
MUON/AliMUONClusterDrawAZ.h [new file with mode: 0644]
MUON/AliMUONClusterFinderAZ.cxx
MUON/AliMUONClusterFinderAZ.h
MUON/AliMUONReconstructor.cxx
MUON/MUONrecLinkDef.h
MUON/libMUONrec.pkg

diff --git a/MUON/AliMUONClusterDrawAZ.cxx b/MUON/AliMUONClusterDrawAZ.cxx
new file mode 100644 (file)
index 0000000..cebbc27
--- /dev/null
@@ -0,0 +1,612 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+// Cluster drawing object for AZ cluster finder 
+
+#include <stdlib.h>
+#include <Riostream.h>
+//#include <TROOT.h>
+#include <TCanvas.h>
+#include <TLine.h>
+//#include <TTree.h>
+#include <TH2.h>
+#include <TView.h>
+#include <TStyle.h>
+
+#include "AliMUONClusterDrawAZ.h"
+#include "AliMUONClusterFinderAZ.h"
+#include "AliHeader.h"
+#include "AliRun.h"
+#include "AliMUON.h"
+//#include "AliMUONChamber.h"
+#include "AliMUONDigit.h"
+#include "AliMUONHit.h"
+#include "AliMUONRawCluster.h"
+//#include "AliMUONClusterInput.h"
+#include "AliMUONPixel.h"
+//#include "AliMC.h"
+#include "AliMUONLoader.h"
+#include "AliLog.h"
+
+ClassImp(AliMUONClusterDrawAZ)
+//_____________________________________________________________________________
+AliMUONClusterDrawAZ::AliMUONClusterDrawAZ()
+  : TObject()
+{
+// Default constructor
+  fFind = NULL; fData = NULL;
+  for (Int_t i=0; i<4; i++) fHist[i] = NULL;
+}
+
+//_____________________________________________________________________________
+AliMUONClusterDrawAZ::AliMUONClusterDrawAZ(AliMUONClusterFinderAZ *clusFinder)
+  : TObject()
+{
+// Constructor
+  fFind = clusFinder;
+  for (Int_t i=0; i<4; i++) fHist[i] = NULL;
+  fDebug = 1; 
+  fEvent = fChamber = 0;
+  Init();
+}
+
+//_____________________________________________________________________________
+AliMUONClusterDrawAZ::~AliMUONClusterDrawAZ()
+{
+  // Destructor
+}
+
+//_____________________________________________________________________________
+AliMUONClusterDrawAZ::AliMUONClusterDrawAZ(const AliMUONClusterDrawAZ& rhs)
+  : TObject(rhs)
+{
+// Protected copy constructor
+
+  AliFatal("Not implemented.");
+}
+
+
+//_____________________________________________________________________________
+AliMUONClusterDrawAZ&  
+AliMUONClusterDrawAZ::operator=(const AliMUONClusterDrawAZ& rhs)
+{
+// Protected assignement operator
+
+  if (this == &rhs) return *this;
+
+  AliFatal("Not implemented.");
+  return *this;  
+}    
+
+//_____________________________________________________________________________
+void AliMUONClusterDrawAZ::Init()
+{
+  // Initialization
+
+  TCanvas *c1 = new TCanvas("c1","Clusters",0,0,600,700);
+  //c1->SetFillColor(10);
+  c1->Divide(1,2);
+  new TCanvas("c2","Mlem",700,0,600,350);
+
+  // Get pointer to Alice detectors
+  //AliMUON *muon  = (AliMUON*) gAlice->GetModule("MUON");
+  //if (!muon) return;
+  //Loaders
+  AliRunLoader *rl = AliRunLoader::GetRunLoader();
+  AliLoader *gime = rl->GetLoader("MUONLoader");
+  fData = ((AliMUONLoader*)gime)->GetMUONData();
+
+  gime->LoadHits("READ"); 
+  gime->LoadRecPoints("READ"); 
+}
+
+//_____________________________________________________________________________
+Bool_t AliMUONClusterDrawAZ::FindEvCh(Int_t nev, Int_t ch)
+{
+  // Find requested event and chamber (skip the ones before the selected)
+
+  if (nev < fEvent) return kFALSE;
+  else if (nev == fEvent && ch < fChamber) return kFALSE;
+  fEvent = nev;
+  fChamber = ch;
+  return kTRUE;
+}
+
+//_____________________________________________________________________________
+void AliMUONClusterDrawAZ::DrawCluster()
+{
+  // Draw preclusters
+
+  TCanvas *c1 = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c1");
+
+  cout << " nev         " << fEvent << endl;
+    
+  char hName[4];
+  for (Int_t cath = 0; cath < 2; cath++) {
+    // Build histograms
+    if (fHist[cath*2]) {fHist[cath*2]->Delete(); fHist[cath*2] = 0;}
+    if (fHist[cath*2+1]) {fHist[cath*2+1]->Delete(); fHist[cath*2+1] = 0;}
+    if (fFind->GetNPads(cath) == 0) continue; // cluster on one cathode only
+    Float_t wxMin = 999, wxMax = 0, wyMin = 999, wyMax = 0; 
+    Int_t minDx = 0, maxDx = 0, minDy = 0, maxDy = 0;
+    for (Int_t i = 0; i < fFind->GetNPads(0)+fFind->GetNPads(1); i++) {
+      if (fFind->GetIJ(0,i) != cath) continue;
+      if (fFind->GetXyq(3,i) < wxMin) { wxMin = fFind->GetXyq(3,i); minDx = i; }
+      if (fFind->GetXyq(3,i) > wxMax) { wxMax = fFind->GetXyq(3,i); maxDx = i; }
+      if (fFind->GetXyq(4,i) < wyMin) { wyMin = fFind->GetXyq(4,i); minDy = i; }
+      if (fFind->GetXyq(4,i) > wyMax) { wyMax = fFind->GetXyq(4,i); maxDy = i; }
+    }
+    cout << minDx << " " << maxDx << " " << minDy << " " << maxDy << endl;
+    Int_t nx, ny, padSize;
+    Float_t xmin = 9999, xmax = -9999, ymin = 9999, ymax = -9999;
+    if (TMath::Nint(fFind->GetXyq(3,minDx)*1000) == TMath::Nint(fFind->GetXyq(3,maxDx)*1000) &&
+       TMath::Nint(fFind->GetXyq(4,minDy)*1000) == TMath::Nint(fFind->GetXyq(4,maxDy)*1000)) {
+      // the same segmentation
+      cout << " Same" << endl;
+      cout << fFind->GetXyq(3,minDx) << " " << fFind->GetXyq(3,maxDx) << " " 
+          << fFind->GetXyq(4,minDy) << " " << fFind->GetXyq(4,maxDy) << endl;
+      for (Int_t i = 0; i < fFind->GetNPads(0)+fFind->GetNPads(1); i++) {
+       if (fFind->GetIJ(0,i) != cath) continue;
+       if (fFind->GetXyq(0,i) < xmin) xmin = fFind->GetXyq(0,i);
+       if (fFind->GetXyq(0,i) > xmax) xmax = fFind->GetXyq(0,i);
+       if (fFind->GetXyq(1,i) < ymin) ymin = fFind->GetXyq(1,i);
+       if (fFind->GetXyq(1,i) > ymax) ymax = fFind->GetXyq(1,i);
+      }
+      xmin -= fFind->GetXyq(3,minDx); xmax += fFind->GetXyq(3,minDx);
+      ymin -= fFind->GetXyq(4,minDy); ymax += fFind->GetXyq(4,minDy);
+      nx = TMath::Nint ((xmax-xmin)/wxMin/2);
+      ny = TMath::Nint ((ymax-ymin)/wyMin/2);
+      cout << xmin << " " << xmax << " " << nx << " " << ymin << " " << ymax << " " << ny << endl;
+      sprintf(hName,"h%d",cath*2);
+      fHist[cath*2] = new TH2D(hName,"cluster",nx,xmin,xmax,ny,ymin,ymax);
+      for (Int_t i = 0; i < fFind->GetNPads(0)+fFind->GetNPads(1); i++) {
+       if (fFind->GetIJ(0,i) != cath) continue;
+       fHist[cath*2]->Fill(fFind->GetXyq(0,i),fFind->GetXyq(1,i),fFind->GetXyq(2,i));
+       //cout << fXyq[0][i] << fXyq[1][i] << fXyq[2][i] << endl;
+      }
+    } else {
+      // different segmentation in the cluster
+      cout << " Different" << endl;
+      cout << fFind->GetXyq(3,minDx) << " " << fFind->GetXyq(3,maxDx) << " " 
+          << fFind->GetXyq(4,minDy) << " " << fFind->GetXyq(4,maxDy) << endl;
+      Int_t nOK = 0;
+      Int_t indx, locMin, locMax;
+      if (TMath::Nint(fFind->GetXyq(3,minDx)*1000) != TMath::Nint(fFind->GetXyq(3,maxDx)*1000)) {
+       // different segmentation along x
+       indx = 0;
+       locMin = minDx;
+       locMax = maxDx;
+      } else {
+       // different segmentation along y
+       indx = 1;
+       locMin = minDy;
+       locMax = maxDy;
+      }
+      Int_t loc = locMin;
+      for (Int_t i = 0; i < 2; i++) {
+       // loop over different pad sizes
+       if (i > 0) loc = locMax;
+       padSize = TMath::Nint(fFind->GetXyq(indx+3,loc)*1000);
+       xmin = 9999; xmax = -9999; ymin = 9999; ymax = -9999;
+       for (Int_t j = 0; j < fFind->GetNPads(0)+fFind->GetNPads(1); j++) {
+         if (fFind->GetIJ(0,j) != cath) continue;
+         if (TMath::Nint(fFind->GetXyq(indx+3,j)*1000) != padSize) continue;
+         nOK++;
+         xmin = TMath::Min (xmin,fFind->GetXyq(0,j));
+         xmax = TMath::Max (xmax,fFind->GetXyq(0,j));
+         ymin = TMath::Min (ymin,fFind->GetXyq(1,j));
+         ymax = TMath::Max (ymax,fFind->GetXyq(1,j));
+       }
+       xmin -= fFind->GetXyq(3,loc); xmax += fFind->GetXyq(3,loc);
+       ymin -= fFind->GetXyq(4,loc); ymax += fFind->GetXyq(4,loc);
+       nx = TMath::Nint ((xmax-xmin)/fFind->GetXyq(3,loc)/2);
+       ny = TMath::Nint ((ymax-ymin)/fFind->GetXyq(4,loc)/2);
+       sprintf(hName,"h%d",cath*2+i);
+       fHist[cath*2+i] = new TH2D(hName,"cluster",nx,xmin,xmax,ny,ymin,ymax);
+       for (Int_t j = 0; j < fFind->GetNPads(0)+fFind->GetNPads(1); j++) {
+         if (fFind->GetIJ(0,j) != cath) continue;
+         if (TMath::Nint(fFind->GetXyq(indx+3,j)*1000) != padSize) continue;
+         fHist[cath*2+i]->Fill(fFind->GetXyq(0,j),fFind->GetXyq(1,j),fFind->GetXyq(2,j));
+       }
+      } // for (Int_t i=0;
+      if (nOK != fFind->GetNPads(cath)) cout << " *** Too many segmentations: nPads, nOK " 
+                                            << fFind->GetNPads(cath) << " " << nOK << endl;
+    } // if (TMath::Nint(fFind->GetXyq(3,minDx)*1000)
+  } // for (Int_t cath = 0;
+       
+  // Draw histograms and coordinates
+  for (Int_t cath = 0; cath < 2; cath++) {
+    if (cath == 0) ModifyHistos();
+    if (fFind->GetNPads(cath) == 0) continue; // cluster on one cathode only
+    c1->cd(cath+1);
+    gPad->SetTheta(55);
+    gPad->SetPhi(30);
+    Double_t x, y, x0, y0, r1 = 999, r2 = 0;
+    if (fHist[cath*2+1]) {
+      // 
+      x0 = fHist[cath*2]->GetXaxis()->GetXmin() - 1000*TMath::Cos(30*TMath::Pi()/180);
+      y0 = fHist[cath*2]->GetYaxis()->GetXmin() - 1000*TMath::Sin(30*TMath::Pi()/180);
+      r1 = 0;
+      Int_t ihist=cath*2;
+      for (Int_t iy = 1; iy <= fHist[ihist]->GetNbinsY(); iy++) {
+       y = fHist[ihist]->GetYaxis()->GetBinCenter(iy) 
+         + fHist[ihist]->GetYaxis()->GetBinWidth(iy);
+       for (Int_t ix = 1; ix <= fHist[ihist]->GetNbinsX(); ix++) {
+         if (fHist[ihist]->GetCellContent(ix,iy) > 0.1) {
+           x = fHist[ihist]->GetXaxis()->GetBinCenter(ix)
+             + fHist[ihist]->GetXaxis()->GetBinWidth(ix);
+           r1 = TMath::Max (r1,TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)));
+         }
+       }
+      }
+      ihist = cath*2 + 1 ;
+      for (Int_t iy = 1; iy <= fHist[ihist]->GetNbinsY(); iy++) {
+       y = fHist[ihist]->GetYaxis()->GetBinCenter(iy)
+         + fHist[ihist]->GetYaxis()->GetBinWidth(iy);
+       for (Int_t ix = 1; ix <= fHist[ihist]->GetNbinsX(); ix++) {
+         if (fHist[ihist]->GetCellContent(ix,iy) > 0.1) {
+           x = fHist[ihist]->GetXaxis()->GetBinCenter(ix)
+             + fHist[ihist]->GetXaxis()->GetBinWidth(ix);
+           r2 = TMath::Max (r2,TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)));
+         }
+       }
+      }
+      cout << r1 << " " << r2 << endl;
+    } // if (fHist[cath*2+1])
+    if (r1 > r2) {
+      //fHist[cath*2]->Draw("lego1");
+      fHist[cath*2]->Draw("lego1Fb");
+      //if (fHist[cath*2+1]) fHist[cath*2+1]->Draw("lego1SameAxisBb");
+      if (fHist[cath*2+1]) fHist[cath*2+1]->Draw("lego1SameAxisBbFb");
+    } else {
+      //fHist[cath*2+1]->Draw("lego1");
+      fHist[cath*2+1]->Draw("lego1Fb");
+      //fHist[cath*2]->Draw("lego1SameAxisBb");
+      fHist[cath*2]->Draw("lego1SameAxisFbBb");
+    }
+    c1->Update();
+  } // for (Int_t cath = 0;
+
+  // Draw simulated and reconstructed hits 
+  DrawHits();
+}
+
+//_____________________________________________________________________________
+void AliMUONClusterDrawAZ::DrawHits()
+{
+  // Draw simulated and reconstructed hits 
+
+  TView *view[2] = { 0x0, 0x0 };
+  Double_t p1[3]={0}, p2[3], xNDC[6];
+  TLine *line[99] = {0};
+  TCanvas *c1 = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c1");
+  if (c1) {
+    c1->cd(1);
+    view[0] = c1->Pad()->GetView();
+    c1->cd(2);
+    view[1] = c1->Pad()->GetView();
+  }
+  fData->SetTreeAddress("H RC");
+
+  TH2D *hist = fHist[0] ? fHist[0] : fHist[2];
+  p2[2] = hist->GetMaximum();
+
+  // Draw simulated hits
+  cout << " *** Simulated hits *** " << endl;
+  Int_t ntracks = (Int_t) fData->GetNtracks();
+  fnMu = 0;
+  Int_t ix, iy, iok, nLine = 0;
+  TClonesArray *hits = NULL;
+  for (Int_t i = 0; i < ntracks; i++) {
+    fData->GetTrack(i);
+    hits = fData->Hits();
+    Int_t nhits = (Int_t) hits->GetEntriesFast();
+    AliMUONHit* mHit;
+    for (Int_t ihit = 0; ihit < nhits; ihit++) {
+      mHit = (AliMUONHit*) hits->UncheckedAt(ihit);
+      if (mHit->Chamber() != fChamber+1) continue;  // chamber number
+      if (TMath::Abs(mHit->Z()-fFind->GetZpad()) > 1) continue; // different slat
+      p2[0] = p1[0] = mHit->X();        // x-pos of hit
+      p2[1] = p1[1] = mHit->Y();        // y-pos
+      if (p1[0] < hist->GetXaxis()->GetXmin() || 
+         p1[0] > hist->GetXaxis()->GetXmax()) continue;
+      if (p1[1] < hist->GetYaxis()->GetXmin() || 
+         p1[1] > hist->GetYaxis()->GetXmax()) continue;
+      // Check if track comes thru pads with signal
+      iok = 0;
+      for (Int_t ihist = 0; ihist < 4; ihist++) {
+       if (!fHist[ihist]) continue;
+       ix = fHist[ihist]->GetXaxis()->FindBin(p1[0]);
+       iy = fHist[ihist]->GetYaxis()->FindBin(p1[1]);
+       if (fHist[ihist]->GetCellContent(ix,iy) > 0.5) {iok = 1; break;}
+      }
+      if (!iok) continue;
+      gStyle->SetLineColor(1);
+      if (TMath::Abs((Int_t)mHit->Particle()) == 13) {
+       gStyle->SetLineColor(4);
+       if (fnMu < 2) {
+         fxyMu[fnMu][0] = p1[0];
+         fxyMu[fnMu++][1] = p1[1];
+       }
+      }            
+      if (fDebug) printf(" X=%10.4f, Y=%10.4f, Z=%10.4f\n",p1[0],p1[1],mHit->Z());
+      if (view[0] || view[1]) {
+       // Take into account track angles
+       p2[0] += mHit->Tlength() * TMath::Sin(mHit->Theta()/180*TMath::Pi()) 
+                                * TMath::Cos(mHit->Phi()/180*TMath::Pi()) / 2;
+       p2[1] += mHit->Tlength() * TMath::Sin(mHit->Theta()/180*TMath::Pi()) 
+                                * TMath::Sin(mHit->Phi()/180*TMath::Pi()) / 2;
+       for (Int_t ipad = 1; ipad < 3; ipad++) {
+         c1->cd(ipad);
+         view[ipad-1]->WCtoNDC(p1, &xNDC[0]);
+         view[ipad-1]->WCtoNDC(p2, &xNDC[3]);
+         //c1->DrawLine(xpad[0],xpad[1],xpad[3],xpad[4]);
+         line[nLine] = new TLine(xNDC[0],xNDC[1],xNDC[3],xNDC[4]);
+         line[nLine++]->Draw();
+       }
+      }
+    } // for (Int_t ihit = 0; ihit < nhits;
+  } // for (Int_t i = 0; i < ntracks;
+
+  // Draw reconstructed coordinates
+  fData->GetRawClusters();
+  TClonesArray *rawclust = fData->RawClusters(fChamber);
+  AliMUONRawCluster *mRaw;
+  gStyle->SetLineColor(3);
+  cout << " *** Reconstructed hits *** " << endl;
+  if (rawclust) {
+    for (Int_t i = 0; i < rawclust ->GetEntriesFast(); i++) {
+      mRaw = (AliMUONRawCluster*)rawclust ->UncheckedAt(i);
+      if (TMath::Abs(mRaw->GetZ(0)-fFind->GetZpad()) > 1) continue; // different slat
+      p2[0] = p1[0] = mRaw->GetX(0);        // x-pos of hit
+      p2[1] = p1[1] = mRaw->GetY(0);        // y-pos
+      if (p1[0] < hist->GetXaxis()->GetXmin() || 
+         p1[0] > hist->GetXaxis()->GetXmax()) continue;
+      if (p1[1] < hist->GetYaxis()->GetXmin() || 
+         p1[1] > hist->GetYaxis()->GetXmax()) continue;
+      /*
+       treeD->GetEvent(cath);
+       cout << mRaw->fMultiplicity[0] << mRaw->fMultiplicity[1] << endl;
+       for (Int_t j=0; j<mRaw->fMultiplicity[cath]; j++) {
+       Int_t digit = mRaw->fIndexMap[j][cath];
+       cout << ((AliMUONDigit*)fMuonDigits->UncheckedAt(digit))->Signal() << endl;
+       }
+      */
+      // Check if track comes thru pads with signal
+      iok = 0;
+      for (Int_t ihist = 0; ihist < 4; ihist++) {
+       if (!fHist[ihist]) continue;
+       ix = fHist[ihist]->GetXaxis()->FindBin(p1[0]);
+       iy = fHist[ihist]->GetYaxis()->FindBin(p1[1]);
+       if (fHist[ihist]->GetCellContent(ix,iy) > 0.5) {iok = 1; break;}
+      }
+      if (!iok) continue;
+      if (fDebug) printf(" X=%10.4f, Y=%10.4f, Z=%10.4f\n",p1[0],p1[1],mRaw->GetZ(0));
+      if (view[0] || view[1]) {
+       for (Int_t ipad = 1; ipad < 3; ipad++) {
+         c1->cd(ipad);
+         view[ipad-1]->WCtoNDC(p1, &xNDC[0]);
+         view[ipad-1]->WCtoNDC(p2, &xNDC[3]);
+         line[nLine] = new TLine(xNDC[0],xNDC[1],xNDC[3],xNDC[4]);
+         line[nLine++]->Draw();
+       }
+      }
+    } // for (Int_t i = 0; i < rawclust ->GetEntries();
+  } // if (rawclust)
+  c1->Update();
+}
+
+//_____________________________________________________________________________
+Int_t AliMUONClusterDrawAZ::Next()
+{
+  // What to do next?
+  // File
+  FILE *lun = 0;
+  //lun = fopen("pull.dat","w");
+
+  for (Int_t i = 0; i < fnMu; i++) {
+    // Check again if muon comes thru the used pads (due to extra splitting)
+    for (Int_t j = 0; j < fFind->GetNPads(0)+fFind->GetNPads(1); j++) {
+      if (TMath::Abs(fxyMu[i][0]-fFind->GetXyq(0,j))<fFind->GetXyq(3,j) && 
+         TMath::Abs(fxyMu[i][1]-fFind->GetXyq(1,j))<fFind->GetXyq(4,j)) {
+       if (fDebug) printf("%12.3e %12.3e %12.3e %12.3e\n",fxyMu[i][2],fxyMu[i][3],fxyMu[i][4],fxyMu[i][5]);
+       if (lun) fprintf(lun,"%4d %2d %12.3e %12.3e %12.3e %12.3e\n",fEvent,fChamber,fxyMu[i][2],fxyMu[i][3],fxyMu[i][4],fxyMu[i][5]);
+       break;
+      }
+    }
+  } // for (Int_t i=0; i<fnMu;
+
+  // What's next?
+  char command[8];
+  cout << " What is next? " << endl;
+  command[0] = ' '; 
+  gets(command);
+  if (command[0] == 'n' || command[0] == 'N') { fEvent++; fChamber = 0; } // next event
+  else if (command[0] == 'q' || command[0] == 'Q') { if (lun) fclose(lun); } // exit display 
+  else if (command[0] == 'c' || command[0] == 'C') sscanf(command+1,"%d",&fChamber); // new chamber
+  else if (command[0] == 'e' || command[0] == 'E') { sscanf(command+1,"%d",&fEvent); fChamber = 0; } // new event
+  else return 1; // Next precluster
+  return 0;
+}
+
+
+//_____________________________________________________________________________
+void AliMUONClusterDrawAZ::ModifyHistos(void)
+{
+  // Modify histograms to bring them to (approximately) the same size
+  Int_t nhist = 0;
+  Float_t hlim[4][4], hbin[4][4]; // first index - xmin, xmax, ymin, ymax
+
+  Float_t binMin[4] = {999,999,999,999};
+
+  for (Int_t i = 0; i < 4; i++) {
+    if (!fHist[i]) {
+      hlim[0][i] = hlim[2][i] = 999;
+      hlim[1][i] = hlim[3][i] = -999;
+      continue;
+    }
+    hlim[0][i] = fHist[i]->GetXaxis()->GetXmin(); // xmin
+    hlim[1][i] = fHist[i]->GetXaxis()->GetXmax(); // xmax
+    hlim[2][i] = fHist[i]->GetYaxis()->GetXmin(); // ymin
+    hlim[3][i] = fHist[i]->GetYaxis()->GetXmax(); // ymax
+    hbin[0][i] = hbin[1][i] = fHist[i]->GetXaxis()->GetBinWidth(1);
+    hbin[2][i] = hbin[3][i] = fHist[i]->GetYaxis()->GetBinWidth(1);
+    binMin[0] = TMath::Min(binMin[0],hbin[0][i]);
+    binMin[2] = TMath::Min(binMin[2],hbin[2][i]);
+    nhist++;
+  }
+  binMin[1] = binMin[0];
+  binMin[3] = binMin[2];
+  cout << " Nhist: " << nhist << endl;
+
+  // Adjust histo limits for cathode with different segmentation
+  for (Int_t i = 0; i < 4; i+=2) {
+    if (!fHist[i+1]) continue;
+    Int_t imin, imax, i1 = i + 1;
+    for (Int_t lim = 0; lim < 4; lim++) {
+      while (1) {
+       if (hlim[lim][i] < hlim[lim][i1]) {
+         imin = i;
+         imax = i1;
+       } else {
+         imin = i1;
+         imax = i;
+       }
+       if (TMath::Abs(hlim[lim][imin]-hlim[lim][imax])<0.01*binMin[lim]) break;
+       if (lim == 0 || lim == 2) {
+         // find lower limit
+         hlim[lim][imax] -= hbin[lim][imax];
+       } else {
+         // find upper limit
+         hlim[lim][imin] += hbin[lim][imin];
+       }
+      } // while (1)
+    }
+  }
+    
+
+  Int_t imnmx = 0, nExtra = 0;
+  for (Int_t lim = 0; lim < 4; lim++) {
+    if (lim == 0 || lim == 2) imnmx = TMath::LocMin(4,hlim[lim]); // find lower limit
+    else imnmx = TMath::LocMax(4,hlim[lim]); // find upper limit
+
+    // Adjust histogram limit
+    for (Int_t i = 0; i < 4; i++) {
+      if (!fHist[i]) continue;
+      nExtra = TMath::Nint ((hlim[lim][imnmx]-hlim[lim][i]) / hbin[lim][i]);
+      hlim[lim][i] += nExtra * hbin[lim][i];
+    }
+  }
+    
+  // Rebuild histograms 
+  TH2D *hist = 0;
+  Int_t nx, ny;
+  Double_t x, y, cont, cmax=0;
+  char hName[4];
+  for (Int_t ihist = 0; ihist < 4; ihist++) {
+    if (!fHist[ihist]) continue;
+    nx = TMath::Nint((hlim[1][ihist]-hlim[0][ihist])/hbin[0][ihist]);
+    ny = TMath::Nint((hlim[3][ihist]-hlim[2][ihist])/hbin[2][ihist]);
+    cout << ihist << " " << hlim[0][ihist] << " " << hlim[1][ihist] << " " << nx;
+    cout << " " << hlim[2][ihist] << " " << hlim[3][ihist] << " " << ny << endl;
+    sprintf(hName,"hh%d",ihist);
+    hist =  new TH2D(hName,"hist",nx,hlim[0][ihist],hlim[1][ihist],ny,hlim[2][ihist],hlim[3][ihist]);
+    for (Int_t i=1; i<=fHist[ihist]->GetNbinsX(); i++) {
+      x = fHist[ihist]->GetXaxis()->GetBinCenter(i);
+      for (Int_t j=1; j<=fHist[ihist]->GetNbinsY(); j++) {
+       y = fHist[ihist]->GetYaxis()->GetBinCenter(j);
+       cont = fHist[ihist]->GetCellContent(i,j);
+       hist->Fill(x,y,cont);
+      }
+    }
+    cmax = TMath::Max (cmax,hist->GetMaximum());
+    sprintf(hName,"%s%d",fHist[ihist]->GetName(),ihist);
+    fHist[ihist]->Delete();
+    fHist[ihist] = new TH2D(*hist);
+    fHist[ihist]->SetName(hName);
+    fHist[ihist]->SetNdivisions(505,"Z");
+    hist->Delete(); 
+  }
+  if (fDebug) printf("%f \n",cmax);
+
+  for (Int_t ihist = 0; ihist < 4; ihist++) {
+    if (!fHist[ihist]) continue;
+    fHist[ihist]->SetMaximum(cmax);
+    fHist[ihist]->SetMinimum(0);
+  }
+}
+
+//_____________________________________________________________________________
+void AliMUONClusterDrawAZ::AdjustHist(Double_t *xylim, const AliMUONPixel *pixPtr)
+{
+  // Adjust histogram limits for pixel drawing
+
+  Float_t xypads[4];
+  if (fHist[0]) {
+    xypads[0] = fHist[0]->GetXaxis()->GetXmin();
+    xypads[1] = -fHist[0]->GetXaxis()->GetXmax();
+    xypads[2] = fHist[0]->GetYaxis()->GetXmin();
+    xypads[3] = -fHist[0]->GetYaxis()->GetXmax();
+    for (Int_t i = 0; i < 4; i++) {
+      while(1) {
+       if (xylim[i] < xypads[i]) break;
+       xylim[i] -= 2*pixPtr->Size(i/2);
+      }
+    }
+  } 
+}
+
+//_____________________________________________________________________________
+void AliMUONClusterDrawAZ::DrawHist(const char* canvas, TH2D *hist)
+{
+  // Draw histogram in given canvas 
+
+  Int_t ix = 0;
+  //((TCanvas*)gROOT->FindObject("c2"))->cd();
+  ((TCanvas*)gROOT->FindObject(canvas))->cd();
+  gPad->SetTheta(55);
+  gPad->SetPhi(30);
+  hist->Draw("lego1Fb");
+  gPad->Update();
+  gets((char*)&ix);
+}
+
+//_____________________________________________________________________________
+void AliMUONClusterDrawAZ::FillMuon(Int_t nfit, const Double_t *parOk, const Double_t *errOk)
+{
+  // Fill muon information
+
+  Int_t indx, imax;
+  Double_t cmax, rad;
+  for (Int_t i = 0; i < fnMu; i++) {
+    cmax = fxyMu[i][6];
+    for (Int_t j = 0; j < nfit; j++) {
+      indx = j<2 ? j*2 : j*2+1;  
+      rad = (fxyMu[i][0]-parOk[indx])*(fxyMu[i][0]-parOk[indx]) +
+            (fxyMu[i][1]-parOk[indx+1])*(fxyMu[i][1]-parOk[indx+1]);
+      if (rad < cmax) {
+       cmax = rad; 
+       imax = indx;
+       fxyMu[i][6] = cmax;
+       fxyMu[i][2] = parOk[imax] - fxyMu[i][0];
+       fxyMu[i][4] = parOk[imax+1] - fxyMu[i][1];
+       fxyMu[i][3] = errOk[imax];
+       fxyMu[i][5] = errOk[imax+1];
+      }
+    }      
+  }
+}
+
diff --git a/MUON/AliMUONClusterDrawAZ.h b/MUON/AliMUONClusterDrawAZ.h
new file mode 100644 (file)
index 0000000..7a00734
--- /dev/null
@@ -0,0 +1,55 @@
+#ifndef ALIMUONCLUSTERDRAWAZ_H
+#define ALIMUONCLUSTERDRAWAZ_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/// \ingroup rec
+/// \class AliMUONClusterDrawAZ
+/// \brief Cluster drawing object for AZ cluster finder in MUON arm of ALICE
+
+#include "AliMUONClusterDrawAZ.h"
+
+class TH2D;
+class AliMUONData;
+class AliMUONPixel;
+class AliMUONClusterFinderAZ;
+
+class AliMUONClusterDrawAZ : public TObject 
+{
+public:
+  AliMUONClusterDrawAZ(); // default constructor
+  AliMUONClusterDrawAZ(AliMUONClusterFinderAZ *clusFinder); // Constructor
+  virtual ~AliMUONClusterDrawAZ(); // Destructor
+
+  void     DrawCluster(); // draw precluster
+  void     AdjustHist(Double_t *xylim, const AliMUONPixel *pixPtr);
+  void     DrawHist(const char* canvas, TH2D *hist); // draw histogram in canvas
+  Int_t    Next(); // commands for drawing
+  Bool_t   FindEvCh(Int_t nev, Int_t ch); // find requested event and chamber
+  void     FillMuon(Int_t nfit, const Double_t *parOk, const Double_t *errOk); // fill muon info
+  void     ResetMuon() { fxyMu[0][6] = fxyMu[1][6] = 9999; } // reset muons
+
+protected:
+  AliMUONClusterDrawAZ(const AliMUONClusterDrawAZ& rhs);
+  AliMUONClusterDrawAZ& operator=(const AliMUONClusterDrawAZ& rhs);
+
+private:
+  AliMUONData *fData; //! pointer to muon data container
+  AliMUONClusterFinderAZ* fFind; //! pointer to ClusterFinder
+  TH2D*      fHist[4]; // ! histograms
+  Int_t      fnMu; // ! number of muons passing thru the selected area
+  Double_t   fxyMu[2][7]; // ! muon information
+  Int_t      fEvent; // ! current event
+  Int_t      fChamber; //! current chamber
+  Int_t      fDebug; // ! debug level
+
+  // Functions
+
+  void   Init(); // initialization
+  void   ModifyHistos(); // modify histograms
+  void   DrawHits(); // draw simulated and reconstructed hits
+
+ClassDef(AliMUONClusterDrawAZ,0) // cluster drawing for MUON arm of ALICE
+};
+
+#endif
index 58695abdc76452660580676fc3245f055eb07880..3ed43b5ffe2cff2a2b02e8a78fdbc3675abf2ac2 100644 (file)
 
 /* $Id$ */
 
-// Clusterizer class developped by A. Zinchenko (Dubna)
+// Clusterizer class developed by A. Zinchenko (Dubna)
 
 #include <stdlib.h>
 #include <Riostream.h>
-#include <TROOT.h>
-#include <TCanvas.h>
-#include <TLine.h>
-#include <TTree.h>
+//#include <TROOT.h>
 #include <TH2.h>
-#include <TView.h>
-#include <TStyle.h>
 #include <TMinuit.h>
 #include <TMatrixD.h>
 
 #include "AliMUONClusterFinderAZ.h"
+#include "AliMUONClusterDrawAZ.h"
 #include "AliHeader.h"
 #include "AliRun.h"
 #include "AliMUON.h"
-#include "AliMUONChamber.h"
+//#include "AliMUONChamber.h"
 #include "AliMUONDigit.h"
-#include "AliMUONHit.h"
-#include "AliMUONChamber.h"
+//#include "AliMUONHit.h"
 #include "AliMUONRawCluster.h"
 #include "AliMUONClusterInput.h"
 #include "AliMUONPixel.h"
-#include "AliMC.h"
-#include "AliMUONLoader.h"
+//#include "AliMC.h"
+//#include "AliMUONLoader.h"
 #include "AliLog.h"
 
 ClassImp(AliMUONClusterFinderAZ)
@@ -52,22 +47,19 @@ ClassImp(AliMUONClusterFinderAZ)
 //FILE *lun1 = fopen("nxny.dat","w");
 
 //_____________________________________________________________________________
-AliMUONClusterFinderAZ::AliMUONClusterFinderAZ(Bool_t draw, Int_t iReco)
+AliMUONClusterFinderAZ::AliMUONClusterFinderAZ(Bool_t draw)
   : AliMUONClusterFinderVS()
 {
 // Constructor
-  for (Int_t i=0; i<4; i++) {fHist[i] = 0;}
-  //fMuonDigits = 0;
   fSegmentation[1] = fSegmentation[0] = 0; 
-  //AZ fgClusterFinder = 0x0;
-  //AZ fgMinuit = 0x0; 
   if (!fgClusterFinder) fgClusterFinder = this;
   if (!fgMinuit) fgMinuit = new TMinuit(8);
-  fDraw = draw;
-  fReco = iReco;
   fPixArray = new TObjArray(20); 
   fDebug = 0; //0;
-  if (draw) fDebug = 1;
+  if (draw) {
+    fDebug = 1;
+    fDraw = new AliMUONClusterDrawAZ(this);
+  }
   AliWarning("*** Running AZ cluster finder ***");
 }
 
@@ -85,6 +77,7 @@ AliMUONClusterFinderAZ::~AliMUONClusterFinderAZ()
 {
   // Destructor
   delete fgMinuit; fgMinuit = 0; delete fPixArray; fPixArray = 0;
+  delete fDraw;
 }
 
 //_____________________________________________________________________________
@@ -92,54 +85,41 @@ void AliMUONClusterFinderAZ::FindRawClusters()
 {
 // To provide the same interface as in AliMUONClusterFinderVS
 
-  ResetRawClusters(); //AZ
+  ResetRawClusters(); 
   EventLoop (gAlice->GetHeader()->GetEvent(), AliMUONClusterInput::Instance()->Chamber());
 }
 
 //_____________________________________________________________________________
-void AliMUONClusterFinderAZ::EventLoop(Int_t nev=0, Int_t ch=0)
+void AliMUONClusterFinderAZ::EventLoop(Int_t nev, Int_t ch)
 {
 // Loop over digits
   
-  static Int_t nev0 = -1, ch0 = -1;
-  if (fDraw) {
-    // Find requested event and chamber 
-    if (nev < nev0) return;
-    else if (nev == nev0 && ch < ch0) return;
-  }
-  nev0 = nev;
-  ch0 = ch;
+  if (fDraw && !fDraw->FindEvCh(nev, ch)) return;
 
   AliMUON *pMuon = (AliMUON*) gAlice->GetModule("MUON");
   AliMUONChamber *iChamber = &(pMuon->Chamber(ch));
-  //fSegmentation[0] = iChamber->SegmentationModel(1);
-  //fSegmentation[1] = iChamber->SegmentationModel(2);
   fResponse = iChamber->ResponseModel();
   fSegmentation[0] = AliMUONClusterInput::Instance()->Segmentation2(0);
   fSegmentation[1] = AliMUONClusterInput::Instance()->Segmentation2(1);
   //AZ fResponse = AliMUONClusterInput::Instance()->Response();
     
-  Int_t ndigits[2]={9,9}, nShown[2]={0};
-  for (Int_t i=0; i<2; i++) {
-    for (Int_t j=0; j<fgkDim; j++) {fUsed[i][j]=kFALSE;}
+  Int_t ndigits[2] = {9,9}, nShown[2] = {0};
+  for (Int_t i = 0; i < 2; i++) {
+    for (Int_t j = 0; j < fgkDim; j++) { fUsed[i][j] = kFALSE; }
   }
 
 next:
-  if (ndigits[0] == nShown[0] && ndigits[1] == nShown[1]) {
-    // No more clusters
-    if (fReco) return;
-    //AZ ch0++;
-    return; // next chamber
-  }
+  if (ndigits[0] == nShown[0] && ndigits[1] == nShown[1]) return;
+
   Float_t xpad, ypad, zpad, zpad0;
   Bool_t first = kTRUE;
   if (fDebug) cout << " *** Event # " << nev << " chamber: " << ch << endl;
   fnPads[0] = fnPads[1] = 0;
-  for (Int_t i=0; i<fgkDim; i++) {fPadIJ[1][i] = 0;}
+  for (Int_t i = 0; i < fgkDim; i++) { fPadIJ[1][i] = 0; }
 
-  for (Int_t iii = 0; iii<2; iii++) { 
+  for (Int_t iii = 0; iii < 2; iii++) { 
     Int_t cath = TMath::Odd(iii);
-    ndigits[cath] = AliMUONClusterInput::Instance()->NDigits(cath); //AZ
+    ndigits[cath] = AliMUONClusterInput::Instance()->NDigits(cath); 
     if (!ndigits[0] && !ndigits[1]) return;
     if (ndigits[cath] == 0) continue;
     if (fDebug) cout << " ndigits: " << ndigits[cath] << " " << cath << endl;
@@ -194,7 +174,7 @@ next:
   } // for (Int_t iii = 0;
 
   fZpad = zpad0;
-  if (fDraw) DrawCluster(nev0, ch0); 
+  if (fDraw) fDraw->DrawCluster(); 
 
   // Use MLEM for cluster finder
   Int_t nMax = 1, localMax[100], maxPos[100];
@@ -220,446 +200,13 @@ next:
       }
     }
   }
-  if (fReco || Next(nev0, ch0)) goto next;
-}
-
-//_____________________________________________________________________________
-void AliMUONClusterFinderAZ::DrawCluster(Int_t nev0, Int_t ch0)
-{
-  // Draw preclusters
-  TCanvas *c1 = 0;
-  TView *view = 0;
-  TH2F *hist = 0;
-  Double_t p1[3]={0}, p2[3];
-  if (!gPad) {
-    c1 = new TCanvas("c1","Clusters",0,0,600,700);
-    //c1->SetFillColor(10);
-    c1->Divide(1,2);
-    new TCanvas("c2","Mlem",700,0,600,350);
-  } else {
-    c1 = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c1");
-  }    
-
-  Int_t ntracks = 0;
-
-  // Get pointer to Alice detectors
-  AliMUON *muon  = (AliMUON*) gAlice->GetModule("MUON");
-  if (!muon) return;
-
-  //Loaders
-  AliRunLoader *rl = AliRunLoader::GetRunLoader();
-  AliLoader *gime = rl->GetLoader("MUONLoader");
-  AliMUONData *data = ((AliMUONLoader*)gime)->GetMUONData();
-
-  gime->LoadHits("READ"); 
-  TTree *treeH = gime->TreeH();
-  ntracks = (Int_t) treeH->GetEntries();
-  cout << " nev         " << nev0 << " ntracks " << ntracks << endl;
-  gime->LoadRecPoints("READ"); 
-  TTree *treeR = data->TreeR();
-  if (treeR) {
-    data->SetTreeAddress("RC");
-    data->GetRawClusters();
-  }
-    
-  TLine *line[99]={0};
-  Int_t nLine = 0;
-  char hName[4];
-  for (Int_t cath = 0; cath<2; cath++) {
-    // Build histograms
-    if (fHist[cath*2]) {fHist[cath*2]->Delete(); fHist[cath*2] = 0;}
-    if (fHist[cath*2+1]) {fHist[cath*2+1]->Delete(); fHist[cath*2+1] = 0;}
-    if (fnPads[cath] == 0) continue; // cluster on one cathode only
-    Float_t wxMin=999, wxMax=0, wyMin=999, wyMax=0; 
-    Int_t minDx=0, maxDx=0, minDy=0, maxDy=0;
-    for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
-      if (fPadIJ[0][i] != cath) continue;
-      if (fXyq[3][i] < wxMin) {wxMin = fXyq[3][i]; minDx = i;}
-      if (fXyq[3][i] > wxMax) {wxMax = fXyq[3][i]; maxDx = i;}
-      if (fXyq[4][i] < wyMin) {wyMin = fXyq[4][i]; minDy = i;}
-      if (fXyq[4][i] > wyMax) {wyMax = fXyq[4][i]; maxDy = i;}
-    }
-    cout << minDx << maxDx << minDy << maxDy << endl;
-    Int_t nx, ny, padSize;
-    Float_t xmin=9999, xmax=-9999, ymin=9999, ymax=-9999;
-    if (TMath::Nint(fXyq[3][minDx]*1000) == TMath::Nint(fXyq[3][maxDx]*1000) &&
-       TMath::Nint(fXyq[4][minDy]*1000) == TMath::Nint(fXyq[4][maxDy]*1000)) {
-      // the same segmentation
-      cout << " Same" << endl;
-      cout << fXyq[3][minDx] << " " << fXyq[3][maxDx] << " " << fXyq[4][minDy] << " " << fXyq[4][maxDy] << endl;
-      for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
-       if (fPadIJ[0][i] != cath) continue;
-       if (fXyq[0][i] < xmin) xmin = fXyq[0][i];
-       if (fXyq[0][i] > xmax) xmax = fXyq[0][i];
-       if (fXyq[1][i] < ymin) ymin = fXyq[1][i];
-       if (fXyq[1][i] > ymax) ymax = fXyq[1][i];
-      }
-      xmin -= fXyq[3][minDx]; xmax += fXyq[3][minDx];
-      ymin -= fXyq[4][minDy]; ymax += fXyq[4][minDy];
-      nx = TMath::Nint ((xmax-xmin)/wxMin/2);
-      ny = TMath::Nint ((ymax-ymin)/wyMin/2);
-      cout << xmin << " " << xmax << " " << nx << " " << ymin << " " << ymax << " " << ny << endl;
-      sprintf(hName,"h%d",cath*2);
-      fHist[cath*2] = new TH2F(hName,"cluster",nx,xmin,xmax,ny,ymin,ymax);
-      //cout << fHist[cath*2] << " " << fnPads[cath] << endl;
-      for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
-       if (fPadIJ[0][i] != cath) continue;
-       fHist[cath*2]->Fill(fXyq[0][i],fXyq[1][i],fXyq[2][i]);
-       //cout << fXyq[0][i] << fXyq[1][i] << fXyq[2][i] << endl;
-      }
-    } else {
-      // different segmentation in the cluster
-      cout << " Different" << endl;
-      cout << fXyq[3][minDx] << " " << fXyq[3][maxDx] << " " << fXyq[4][minDy] << " " << fXyq[4][maxDy] << endl;
-      Int_t nOK = 0;
-      Int_t indx, locMin, locMax;
-      if (TMath::Nint(fXyq[3][minDx]*1000) != TMath::Nint(fXyq[3][maxDx]*1000)) {
-       // different segmentation along x
-       indx = 0;
-       locMin = minDx;
-       locMax = maxDx;
-      } else {
-       // different segmentation along y
-       indx = 1;
-       locMin = minDy;
-       locMax = maxDy;
-      }
-      Int_t loc = locMin;
-      for (Int_t i=0; i<2; i++) {
-       // loop over different pad sizes
-       if (i>0) loc = locMax;
-       padSize = TMath::Nint(fXyq[indx+3][loc]*1000);
-       xmin = 9999; xmax = -9999; ymin = 9999; ymax = -9999;
-       for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
-         if (fPadIJ[0][j] != cath) continue;
-         if (TMath::Nint(fXyq[indx+3][j]*1000) != padSize) continue;
-         nOK++;
-         xmin = TMath::Min (xmin,fXyq[0][j]);
-         xmax = TMath::Max (xmax,fXyq[0][j]);
-         ymin = TMath::Min (ymin,fXyq[1][j]);
-         ymax = TMath::Max (ymax,fXyq[1][j]);
-       }
-       xmin -= fXyq[3][loc]; xmax += fXyq[3][loc];
-       ymin -= fXyq[4][loc]; ymax += fXyq[4][loc];
-       nx = TMath::Nint ((xmax-xmin)/fXyq[3][loc]/2);
-       ny = TMath::Nint ((ymax-ymin)/fXyq[4][loc]/2);
-       sprintf(hName,"h%d",cath*2+i);
-       fHist[cath*2+i] = new TH2F(hName,"cluster",nx,xmin,xmax,ny,ymin,ymax);
-       for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
-         if (fPadIJ[0][j] != cath) continue;
-         if (TMath::Nint(fXyq[indx+3][j]*1000) != padSize) continue;
-         fHist[cath*2+i]->Fill(fXyq[0][j],fXyq[1][j],fXyq[2][j]);
-       }
-      } // for (Int_t i=0;
-      if (nOK != fnPads[cath]) cout << " *** Too many segmentations: nPads, nOK " << fnPads[cath] << " " << nOK << endl;
-    } // if (TMath::Nint(fXyq[3][minDx]*1000)
-  } // for (Int_t cath = 0;
-       
-  // Draw histograms and coordinates
-  for (Int_t cath=0; cath<2; cath++) {
-    if (cath == 0) ModifyHistos();
-    if (fnPads[cath] == 0) continue; // cluster on one cathode only
-    if (fDraw) {
-      c1->cd(cath+1);
-      gPad->SetTheta(55);
-      gPad->SetPhi(30);
-      Double_t x, y, x0, y0, r1=999, r2=0;
-      if (fHist[cath*2+1]) {
-       // 
-       x0 = fHist[cath*2]->GetXaxis()->GetXmin() - 1000*TMath::Cos(30*TMath::Pi()/180);
-       y0 = fHist[cath*2]->GetYaxis()->GetXmin() - 1000*TMath::Sin(30*TMath::Pi()/180);
-       r1 = 0;
-       Int_t ihist=cath*2;
-       for (Int_t iy=1; iy<=fHist[ihist]->GetNbinsY(); iy++) {
-         y = fHist[ihist]->GetYaxis()->GetBinCenter(iy) 
-           + fHist[ihist]->GetYaxis()->GetBinWidth(iy);
-         for (Int_t ix=1; ix<=fHist[ihist]->GetNbinsX(); ix++) {
-           if (fHist[ihist]->GetCellContent(ix,iy) > 0.1) {
-             x = fHist[ihist]->GetXaxis()->GetBinCenter(ix)
-               + fHist[ihist]->GetXaxis()->GetBinWidth(ix);
-             r1 = TMath::Max (r1,TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)));
-           }
-         }
-       }
-       ihist = cath*2 + 1 ;
-       for (Int_t iy=1; iy<=fHist[ihist]->GetNbinsY(); iy++) {
-         y = fHist[ihist]->GetYaxis()->GetBinCenter(iy)
-           + fHist[ihist]->GetYaxis()->GetBinWidth(iy);
-         for (Int_t ix=1; ix<=fHist[ihist]->GetNbinsX(); ix++) {
-           if (fHist[ihist]->GetCellContent(ix,iy) > 0.1) {
-             x = fHist[ihist]->GetXaxis()->GetBinCenter(ix)
-               + fHist[ihist]->GetXaxis()->GetBinWidth(ix);
-             r2 = TMath::Max (r2,TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)));
-           }
-         }
-       }
-       cout << r1 << " " << r2 << endl;
-      } // if (fHist[cath*2+1])
-      if (r1 > r2) {
-       //fHist[cath*2]->Draw("lego1");
-       fHist[cath*2]->Draw("lego1Fb");
-       //if (fHist[cath*2+1]) fHist[cath*2+1]->Draw("lego1SameAxisBb");
-       if (fHist[cath*2+1]) fHist[cath*2+1]->Draw("lego1SameAxisBbFb");
-      } else {
-       //fHist[cath*2+1]->Draw("lego1");
-       fHist[cath*2+1]->Draw("lego1Fb");
-       //fHist[cath*2]->Draw("lego1SameAxisBb");
-       fHist[cath*2]->Draw("lego1SameAxisFbBb");
-      }
-      c1->Update();
-    } // if (fDraw)
-  } // for (Int_t cath = 0;
-
-  // Draw generated hits
-  Double_t xNDC[6];
-  hist = fHist[0] ? fHist[0] : fHist[2];
-  p2[2] = hist->GetMaximum();
-  view = 0;
-  if (c1) view = c1->Pad()->GetView();
-  cout << " *** GEANT hits *** " << endl;
-  fnMu = 0;
-  Int_t ix, iy, iok;
-  for (Int_t i=0; i<ntracks; i++) {
-    treeH->GetEvent(i);
-    for (AliMUONHit* mHit=(AliMUONHit*)muon->FirstHit(-1); 
-        mHit;
-        mHit=(AliMUONHit*)muon->NextHit()) {
-      if (mHit->Chamber() != ch0+1) continue;  // chamber number
-      if (TMath::Abs(mHit->Z()-fZpad) > 1) continue; // different slat
-      p2[0] = p1[0] = mHit->X();        // x-pos of hit
-      p2[1] = p1[1] = mHit->Y();        // y-pos
-      if (p1[0] < hist->GetXaxis()->GetXmin() || 
-         p1[0] > hist->GetXaxis()->GetXmax()) continue;
-      if (p1[1] < hist->GetYaxis()->GetXmin() || 
-         p1[1] > hist->GetYaxis()->GetXmax()) continue;
-      // Check if track comes thru pads with signal
-      iok = 0;
-      for (Int_t ihist=0; ihist<4; ihist++) {
-       if (!fHist[ihist]) continue;
-       ix = fHist[ihist]->GetXaxis()->FindBin(p1[0]);
-       iy = fHist[ihist]->GetYaxis()->FindBin(p1[1]);
-       if (fHist[ihist]->GetCellContent(ix,iy) > 0.5) {iok = 1; break;}
-      }
-      if (!iok) continue;
-      gStyle->SetLineColor(1);
-      if (TMath::Abs((Int_t)mHit->Particle()) == 13) {
-       gStyle->SetLineColor(4);
-       if (fnMu < 2) {
-         fxyMu[fnMu][0] = p1[0];
-         fxyMu[fnMu++][1] = p1[1];
-       }
-      }            
-      if (fDebug) printf(" X=%10.4f, Y=%10.4f, Z=%10.4f\n",p1[0],p1[1],mHit->Z());
-      if (view) {
-       // Take into account track angles
-       p2[0] += mHit->Tlength() * TMath::Sin(mHit->Theta()/180*TMath::Pi()) 
-                                * TMath::Cos(mHit->Phi()/180*TMath::Pi()) / 2;
-       p2[1] += mHit->Tlength() * TMath::Sin(mHit->Theta()/180*TMath::Pi()) 
-                                * TMath::Sin(mHit->Phi()/180*TMath::Pi()) / 2;
-       view->WCtoNDC(p1, &xNDC[0]);
-       view->WCtoNDC(p2, &xNDC[3]);
-       for (Int_t ipad=1; ipad<3; ipad++) {
-         c1->cd(ipad);
-         //c1->DrawLine(xpad[0],xpad[1],xpad[3],xpad[4]);
-         line[nLine] = new TLine(xNDC[0],xNDC[1],xNDC[3],xNDC[4]);
-         line[nLine++]->Draw();
-       }
-      }
-    } // for (AliMUONHit* mHit=
-  } // for (Int_t i=0; i<ntracks;
-
-  // Draw reconstructed coordinates
-  TClonesArray *listMUONrawclust = data->RawClusters(ch0);
-  AliMUONRawCluster *mRaw;
-  gStyle->SetLineColor(3);
-  cout << " *** Reconstructed hits *** " << endl;
-  if (listMUONrawclust) {
-    for (Int_t i=0; i<listMUONrawclust ->GetEntries(); i++) {
-      mRaw = (AliMUONRawCluster*)listMUONrawclust ->UncheckedAt(i);
-      if (TMath::Abs(mRaw->GetZ(0)-fZpad) > 1) continue; // different slat
-      p2[0] = p1[0] = mRaw->GetX(0);        // x-pos of hit
-      p2[1] = p1[1] = mRaw->GetY(0);        // y-pos
-      if (p1[0] < hist->GetXaxis()->GetXmin() || 
-         p1[0] > hist->GetXaxis()->GetXmax()) continue;
-      if (p1[1] < hist->GetYaxis()->GetXmin() || 
-         p1[1] > hist->GetYaxis()->GetXmax()) continue;
-      /*
-       treeD->GetEvent(cath);
-       cout << mRaw->fMultiplicity[0] << mRaw->fMultiplicity[1] << endl;
-       for (Int_t j=0; j<mRaw->fMultiplicity[cath]; j++) {
-       Int_t digit = mRaw->fIndexMap[j][cath];
-       cout << ((AliMUONDigit*)fMuonDigits->UncheckedAt(digit))->Signal() << endl;
-       }
-      */
-      // Check if track comes thru pads with signal
-      iok = 0;
-      for (Int_t ihist=0; ihist<4; ihist++) {
-       if (!fHist[ihist]) continue;
-       ix = fHist[ihist]->GetXaxis()->FindBin(p1[0]);
-       iy = fHist[ihist]->GetYaxis()->FindBin(p1[1]);
-       if (fHist[ihist]->GetCellContent(ix,iy) > 0.5) {iok = 1; break;}
-      }
-      if (!iok) continue;
-      if (fDebug) printf(" X=%10.4f, Y=%10.4f, Z=%10.4f\n",p1[0],p1[1],mRaw->GetZ(0));
-      if (view) {
-       view->WCtoNDC(p1, &xNDC[0]);
-       view->WCtoNDC(p2, &xNDC[3]);
-       for (Int_t ipad=1; ipad<3; ipad++) {
-         c1->cd(ipad);
-         line[nLine] = new TLine(xNDC[0],xNDC[1],xNDC[3],xNDC[4]);
-         line[nLine++]->Draw();
-       }
-      }
-    } // for (Int_t i=0; i<listMUONrawclust ->GetEntries();
-  } // if (listMUONrawclust)
-  c1->Update();
-}
-
-//_____________________________________________________________________________
-Int_t AliMUONClusterFinderAZ::Next(Int_t &nev0, Int_t &ch0)
-{
-  // What to do next?
-  // File
-  FILE *lun = 0;
-  //lun = fopen("pull.dat","w");
-
-  for (Int_t i=0; i<fnMu; i++) {
-    // Check again if muon come thru the used pads (due to extra splitting)
-    for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
-      if (TMath::Abs(fxyMu[i][0]-fXyq[0][j])<fXyq[3][j] && 
-         TMath::Abs(fxyMu[i][1]-fXyq[1][j])<fXyq[4][j]) {
-       if (fDebug) printf("%12.3e %12.3e %12.3e %12.3e\n",fxyMu[i][2],fxyMu[i][3],fxyMu[i][4],fxyMu[i][5]);
-       if (lun) fprintf(lun,"%4d %2d %12.3e %12.3e %12.3e %12.3e\n",nev0,ch0,fxyMu[i][2],fxyMu[i][3],fxyMu[i][4],fxyMu[i][5]);
-       break;
-      }
-    }
-  } // for (Int_t i=0; i<fnMu;
-
-  // What's next?
-  char command[8];
-  cout << " What is next? " << endl;
-  command[0] = ' '; 
-  gets(command);
-  if (command[0] == 'n' || command[0] == 'N') {nev0++; ch0 = 0; } // next event
-  else if (command[0] == 'q' || command[0] == 'Q') { if (lun) fclose(lun); } // exit display 
-  else if (command[0] == 'c' || command[0] == 'C') sscanf(command+1,"%d",&ch0); // new chamber
-  else if (command[0] == 'e' || command[0] == 'E') { sscanf(command+1,"%d",&nev0); ch0 = 0; } // new event
-  else return 1; // Next precluster
-  return 0;
-}
-
-
-//_____________________________________________________________________________
-void AliMUONClusterFinderAZ::ModifyHistos(void)
-{
-  // Modify histograms to bring them to (approximately) the same size
-  Int_t nhist = 0;
-  Float_t hlim[4][4], hbin[4][4]; // first index - xmin, xmax, ymin, ymax
-
-  Float_t binMin[4] = {999,999,999,999};
-
-  for (Int_t i = 0; i < 4; i++) {
-    if (!fHist[i]) {
-      hlim[0][i] = hlim[2][i] = 999;
-      hlim[1][i] = hlim[3][i] = -999;
-      continue;
-    }
-    hlim[0][i] = fHist[i]->GetXaxis()->GetXmin(); // xmin
-    hlim[1][i] = fHist[i]->GetXaxis()->GetXmax(); // xmax
-    hlim[2][i] = fHist[i]->GetYaxis()->GetXmin(); // ymin
-    hlim[3][i] = fHist[i]->GetYaxis()->GetXmax(); // ymax
-    hbin[0][i] = hbin[1][i] = fHist[i]->GetXaxis()->GetBinWidth(1);
-    hbin[2][i] = hbin[3][i] = fHist[i]->GetYaxis()->GetBinWidth(1);
-    binMin[0] = TMath::Min(binMin[0],hbin[0][i]);
-    binMin[2] = TMath::Min(binMin[2],hbin[2][i]);
-    nhist++;
-  }
-  binMin[1] = binMin[0];
-  binMin[3] = binMin[2];
-  cout << " Nhist: " << nhist << endl;
-
-  // Adjust histo limits for cathode with different segmentation
-  for (Int_t i = 0; i < 4; i+=2) {
-    if (!fHist[i+1]) continue;
-    Int_t imin, imax, i1 = i + 1;
-    for (Int_t lim = 0; lim < 4; lim++) {
-      while (1) {
-       if (hlim[lim][i] < hlim[lim][i1]) {
-         imin = i;
-         imax = i1;
-       } else {
-         imin = i1;
-         imax = i;
-       }
-       if (TMath::Abs(hlim[lim][imin]-hlim[lim][imax])<0.01*binMin[lim]) break;
-       if (lim == 0 || lim == 2) {
-         // find lower limit
-         hlim[lim][imax] -= hbin[lim][imax];
-       } else {
-         // find upper limit
-         hlim[lim][imin] += hbin[lim][imin];
-       }
-      } // while (1)
-    }
-  }
-    
-
-  Int_t imnmx = 0, nExtra = 0;
-  for (Int_t lim = 0; lim < 4; lim++) {
-    if (lim == 0 || lim == 2) imnmx = TMath::LocMin(4,hlim[lim]); // find lower limit
-    else imnmx = TMath::LocMax(4,hlim[lim]); // find upper limit
-
-    // Adjust histogram limit
-    for (Int_t i = 0; i < 4; i++) {
-      if (!fHist[i]) continue;
-      nExtra = TMath::Nint ((hlim[lim][imnmx]-hlim[lim][i]) / hbin[lim][i]);
-      hlim[lim][i] += nExtra * hbin[lim][i];
-    }
-  }
-    
-  // Rebuild histograms 
-  TH2F *hist = 0;
-  Int_t nx, ny;
-  Double_t x, y, cont, cmax=0;
-  char hName[4];
-  for (Int_t ihist=0; ihist<4; ihist++) {
-    if (!fHist[ihist]) continue;
-    nx = TMath::Nint((hlim[1][ihist]-hlim[0][ihist])/hbin[0][ihist]);
-    ny = TMath::Nint((hlim[3][ihist]-hlim[2][ihist])/hbin[2][ihist]);
-    cout << ihist << " " << hlim[0][ihist] << " " << hlim[1][ihist] << " " << nx;
-    cout << " " << hlim[2][ihist] << " " << hlim[3][ihist] << " " << ny << endl;
-    sprintf(hName,"hh%d",ihist);
-    hist =  new TH2F(hName,"hist",nx,hlim[0][ihist],hlim[1][ihist],ny,hlim[2][ihist],hlim[3][ihist]);
-    for (Int_t i=1; i<=fHist[ihist]->GetNbinsX(); i++) {
-      x = fHist[ihist]->GetXaxis()->GetBinCenter(i);
-      for (Int_t j=1; j<=fHist[ihist]->GetNbinsY(); j++) {
-       y = fHist[ihist]->GetYaxis()->GetBinCenter(j);
-       cont = fHist[ihist]->GetCellContent(i,j);
-       hist->Fill(x,y,cont);
-      }
-    }
-    cmax = TMath::Max (cmax,hist->GetMaximum());
-    sprintf(hName,"%s%d",fHist[ihist]->GetName(),ihist);
-    fHist[ihist]->Delete();
-    fHist[ihist] = new TH2F(*hist);
-    fHist[ihist]->SetName(hName);
-    fHist[ihist]->SetNdivisions(505,"Z");
-    hist->Delete(); 
-  }
-  if (fDebug) printf("%f \n",cmax);
-
-  for (Int_t ihist=0; ihist<4; ihist++) {
-    if (!fHist[ihist]) continue;
-    fHist[ihist]->SetMaximum(cmax);
-    fHist[ihist]->SetMinimum(0);
-  }
+  if (!fDraw || fDraw->Next()) goto next;
 }
 
 //_____________________________________________________________________________
 void AliMUONClusterFinderAZ::AddPad(Int_t cath, Int_t digit)
 {
   // Add pad to the cluster
-  //AZ AliMUONDigit *mdig = (AliMUONDigit*)fMuonDigits->UncheckedAt(digit);
   AliMUONDigit *mdig = AliMUONClusterInput::Instance()->Digit(cath,digit); //AZ
 
   Int_t charge = mdig->Signal();
@@ -669,7 +216,6 @@ void AliMUONClusterFinderAZ::AddPad(Int_t cath, Int_t digit)
     fUsed[cath][digit] = kTRUE; 
     return; 
   } 
-  
   Int_t isec = fSegmentation[cath]->Sector(fInput->DetElemId(),mdig->PadX(), mdig->PadY());
   Int_t nPads = fnPads[0] + fnPads[1];
   fXyq[0][nPads] = xpad;
@@ -689,17 +235,14 @@ void AliMUONClusterFinderAZ::AddPad(Int_t cath, Int_t digit)
   Int_t nn, ix, iy, xList[10], yList[10];
   AliMUONDigit  *mdig1;
 
-  //AZ Int_t ndigits = fMuonDigits->GetEntriesFast();
-  Int_t ndigits = AliMUONClusterInput::Instance()->NDigits(cath); //AZ
+  Int_t ndigits = AliMUONClusterInput::Instance()->NDigits(cath); 
   fSegmentation[cath]->Neighbours(fInput->DetElemId(),mdig->PadX(),mdig->PadY(),&nn,xList,yList); 
   for (Int_t in=0; in<nn; in++) {
     ix=xList[in];
     iy=yList[in];
     for (Int_t digit1 = 0; digit1 < ndigits; digit1++) {
       if (digit1 == digit) continue;
-      //AZ mdig1 = (AliMUONDigit*)fMuonDigits->UncheckedAt(digit1);
-      mdig1 = AliMUONClusterInput::Instance()->Digit(cath,digit1); //AZ
-      //AZobsolete if (mdig1->Cathode() != cath) continue;
+      mdig1 = AliMUONClusterInput::Instance()->Digit(cath,digit1); 
       if (!fUsed[cath][digit1] && mdig1->PadX() == ix && mdig1->PadY() == iy) {
        //AZ--- temporary fix on edges
        //fSegmentation[cath]->GetPadC(mdig1->PadX(), mdig1->PadY(), xpad, ypad, zpad);
@@ -758,16 +301,6 @@ Bool_t AliMUONClusterFinderAZ::Overlap(Float_t *xy1, Int_t iPad, Float_t *xy12,
   return kTRUE;
 }
 
-//_____________________________________________________________________________
-/*
-Bool_t AliMUONClusterFinderAZ::Overlap(Int_t i, Int_t j, Float_t *xy12, Int_t iSkip)
-{
-  // Check if the pads i and j overlap and return overlap area
-
-  Float_t xy1[4], xy2[4];
-  return Overlap(xy1, xy2, xy12, iSkip);
-}
-*/
 //_____________________________________________________________________________
 Bool_t AliMUONClusterFinderAZ::CheckPrecluster(Int_t *nShown)
 {
@@ -867,7 +400,7 @@ Bool_t AliMUONClusterFinderAZ::CheckPrecluster(Int_t *nShown)
     for (Int_t i=0; i<npad; i++) {
       if (flags[i]) continue;
       nFlags ++;
-      if (fDebug) cout << i << " " << fPadIJ[0][i] << endl;
+      if (fDebug) cout << i << " " << fPadIJ[0][i] << " " << fXyq[0][i] << " " << fXyq[1][i] << endl;
     }
     if (fDebug && nFlags) cout << " nFlags = " << nFlags << endl;
     //if (nFlags > 2 || (Float_t)nFlags / npad > 0.2) { // why 2 ??? - empirical choice
@@ -876,6 +409,9 @@ Bool_t AliMUONClusterFinderAZ::CheckPrecluster(Int_t *nShown)
        if (flags[i]) continue;
        digit = TMath::Nint (fXyq[5][i]);
        cath = fPadIJ[0][i];
+       // Check for edge effect (missing pads on the other cathode)
+       Int_t cath1 = TMath::Even(cath), ix, iy;
+       if (!fSegmentation[cath1]->GetPadI(fInput->DetElemId(),fXyq[0][i],fXyq[1][i],fZpad,ix,iy)) continue;
        fUsed[cath][digit] = kFALSE; // release pad
        fXyq[2][i] = -2;
        fnPads[cath]--;
@@ -1003,7 +539,7 @@ Bool_t AliMUONClusterFinderAZ::CheckPrecluster(Int_t *nShown)
   npad = fnPads[0] + fnPads[1];
   if (npad > 500) { cout << " ***** Too large cluster. Give up. " << npad << endl; return kFALSE; }
   // Back up charge value
-  for (Int_t j=0; j<npad; j++) fXyq[6][j] = fXyq[2][j];
+  for (Int_t j = 0; j < npad; j++) fXyq[6][j] = fXyq[2][j];
 
   return kTRUE;
 }
@@ -1220,7 +756,8 @@ Bool_t AliMUONClusterFinderAZ::MainLoop(Int_t iSimple)
   Double_t *coef = 0, *probi = 0; 
   AddVirtualPad(); // add virtual pads if necessary
   Int_t npadTot = fnPads[0] + fnPads[1], npadOK = 0;
-  for (Int_t i=0; i<npadTot; i++) if (fPadIJ[1][i] == 0) npadOK++;
+  for (Int_t i = 0; i < npadTot; i++) if (fPadIJ[1][i] == 0) npadOK++;
+  if (fDraw) fDraw->ResetMuon();
 
   while (1) {
 
@@ -1278,21 +815,8 @@ Bool_t AliMUONClusterFinderAZ::MainLoop(Int_t iSimple)
 
     // Adjust histogram to approximately the same limits as for the pads
     // (for good presentation)
-    //*
-    Float_t xypads[4];
-    if (fHist[0]) {
-      xypads[0] = fHist[0]->GetXaxis()->GetXmin();
-      xypads[1] = -fHist[0]->GetXaxis()->GetXmax();
-      xypads[2] = fHist[0]->GetYaxis()->GetXmin();
-      xypads[3] = -fHist[0]->GetYaxis()->GetXmax();
-      for (Int_t i=0; i<4; i++) {
-       while(1) {
-         if (xylim[i] < xypads[i]) break;
-         xylim[i] -= 2*pixPtr->Size(i/2);
-       }
-      }
-    } // if (fHist[0])
-    //*/
+    if (fDraw) fDraw->AdjustHist(xylim, pixPtr);
+
     Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
     Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
     
@@ -1301,16 +825,7 @@ Bool_t AliMUONClusterFinderAZ::MainLoop(Int_t iSimple)
       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
       mlem->Fill(pixPtr->Coord(0),pixPtr->Coord(1),pixPtr->Charge());
     }
-    //gPad->GetCanvas()->cd(3);
-    if (fDraw) {
-      ((TCanvas*)gROOT->FindObject("c2"))->cd();
-      gPad->SetTheta(55);
-      gPad->SetPhi(30);
-      //mlem->SetFillColor(19);
-      mlem->Draw("lego1Fb");
-      gPad->Update();
-      gets((char*)&ix);
-    }
+    if (fDraw) fDraw->DrawHist("c2", mlem);
 
     // Check if the total charge of pixels is too low
     Double_t qTot = 0;
@@ -1318,7 +833,7 @@ Bool_t AliMUONClusterFinderAZ::MainLoop(Int_t iSimple)
       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
       qTot += pixPtr->Charge();
     }
-    //AZif (qTot < 1.e-4 || npadOK < 3 && qTot < 50) {
+    //AZ if (qTot < 1.e-4 || npadOK < 3 && qTot < 50) {
     if (qTot < 1.e-4 || npadOK < 3 && qTot < 7) {
       delete [] coef; delete [] probi; coef = 0; probi = 0;
       fPixArray->Delete(); 
@@ -1363,7 +878,7 @@ Bool_t AliMUONClusterFinderAZ::MainLoop(Int_t iSimple)
 
     if (iSimple) {
       // Simple cluster - skip further passes thru EM-procedure
-      fxyMu[0][6] = fxyMu[1][6] = 9999;
+      //fxyMu[0][6] = fxyMu[1][6] = 9999;
       Simple();
       delete [] coef; delete [] probi; coef = 0; probi = 0;
       fPixArray->Delete(); 
@@ -1508,15 +1023,9 @@ Bool_t AliMUONClusterFinderAZ::MainLoop(Int_t iSimple)
     iy = mlem->GetYaxis()->FindBin(pixPtr->Coord(1));
     mlem->SetBinContent(ix, iy, pixPtr->Charge());
   }
-  if (fDraw) {
-    ((TCanvas*)gROOT->FindObject("c2"))->cd();
-    gPad->SetTheta(55);
-    gPad->SetPhi(30);
-    mlem->Draw("lego1Fb");
-    gPad->Update();
-  }
+  if (fDraw) fDraw->DrawHist("c2", mlem);
 
-  fxyMu[0][6] = fxyMu[1][6] = 9999;
+  //fxyMu[0][6] = fxyMu[1][6] = 9999;
   // Try to split into clusters
   Bool_t ok = kTRUE;
   if (mlem->GetSum() < 1) ok = kFALSE;
@@ -2237,6 +1746,7 @@ Int_t AliMUONClusterFinderAZ::Fit(Int_t nfit, Int_t *clustFit, TObjArray **clust
   // Get number of pads in X and Y 
   Int_t nInX = 0, nInY;
   PadsInXandY(nInX, nInY);
+  //cout << " nInX and Y: " << nInX << " " << nInY << endl;
 
   // Take cluster maxima as fitting seeds
   TObjArray *pix;
@@ -2536,42 +2046,23 @@ Int_t AliMUONClusterFinderAZ::Fit(Int_t nfit, Int_t *clustFit, TObjArray **clust
     }
   }
 
-  Int_t indx, imax;
+  Int_t indx;
   fnPads[1] -= nVirtual;
-  if (fReco) {
-    Double_t coef = 0;
-    //for (Int_t j=0; j<nfit; j++) {
-    for (Int_t j=nfit-1; j>=0; j--) {
-      indx = j<2 ? j*2 : j*2+1;  
-      if (nfit == 1) coef = 1;
-      else coef = j==nfit-1 ? parOk[indx+2] : 1-coef;
-      coef = TMath::Max (coef, 0.);
-      if (nfit == 3 && j < 2) coef = j==1 ? coef*parOk[indx+2] : coef - parOk[7];
-      coef = TMath::Max (coef, 0.);
-      AddRawCluster (parOk[indx], parOk[indx+1], coef*fQtot, errOk[indx], nfit0+10*nfit, tracks,
-                    //sigCand[maxSeed[j]][0], sigCand[maxSeed[j]][1]);
-                    //sigCand[0][0], sigCand[0][1], dist[j]);
-                    sigCand[0][0], sigCand[0][1], dist[TMath::LocMin(nfit,dist)]);
-    }
-    return nfit;
-  } 
-  for (Int_t i=0; i<fnMu; i++) {
-    cmax = fxyMu[i][6];
-    for (Int_t j=0; j<nfit; j++) {
-      indx = j<2 ? j*2 : j*2+1;  
-      rad = (fxyMu[i][0]-parOk[indx])*(fxyMu[i][0]-parOk[indx]) +
-            (fxyMu[i][1]-parOk[indx+1])*(fxyMu[i][1]-parOk[indx+1]);
-      if (rad < cmax) {
-       cmax = rad; 
-       imax = indx;
-       fxyMu[i][6] = cmax;
-       fxyMu[i][2] = parOk[imax] - fxyMu[i][0];
-       fxyMu[i][4] = parOk[imax+1] - fxyMu[i][1];
-       fxyMu[i][3] = errOk[imax];
-       fxyMu[i][5] = errOk[imax+1];
-      }
-    }      
+  Double_t coef = 0;
+  //for (Int_t j=0; j<nfit; j++) {
+  for (Int_t j=nfit-1; j>=0; j--) {
+    indx = j<2 ? j*2 : j*2+1;  
+    if (nfit == 1) coef = 1;
+    else coef = j==nfit-1 ? parOk[indx+2] : 1-coef;
+    coef = TMath::Max (coef, 0.);
+    if (nfit == 3 && j < 2) coef = j==1 ? coef*parOk[indx+2] : coef - parOk[7];
+    coef = TMath::Max (coef, 0.);
+    AddRawCluster (parOk[indx], parOk[indx+1], coef*fQtot, errOk[indx], nfit0+10*nfit, tracks,
+                  //sigCand[maxSeed[j]][0], sigCand[maxSeed[j]][1]);
+                  //sigCand[0][0], sigCand[0][1], dist[j]);
+                  sigCand[0][0], sigCand[0][1], dist[TMath::LocMin(nfit,dist)]);
   }
+  if (fDraw) fDraw->FillMuon(nfit, parOk, errOk);
   return nfit;
 }  
 
@@ -2740,15 +2231,7 @@ Int_t AliMUONClusterFinderAZ::FindLocalMaxima(Int_t *localMax, Double_t *maxVal)
     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
     hist->Fill(pixPtr->Coord(0), pixPtr->Coord(1), pixPtr->Charge());
   }
-  if (fDraw) {
-    ((TCanvas*)gROOT->FindObject("c2"))->cd();
-    gPad->SetTheta(55);
-    gPad->SetPhi(30);
-    hist->Draw("lego1Fb");
-    gPad->Update();
-    int ia;
-    cin >> ia;
-  }
+  if (fDraw) fDraw->DrawHist("c2", hist);
 
   Int_t nMax = 0, indx;
   Int_t *isLocalMax = new Int_t[ny*nx];
@@ -2934,16 +2417,17 @@ void AliMUONClusterFinderAZ::AddVirtualPad()
       }
     }
 
-    Int_t iAddX = 0, iAddY = 0, ix1 = 0, iy1 = 0, iMuon = 0, iPad = 0;
+    Int_t iAddX = 0, iAddY = 0, ix1 = 0, iy1 = 0, iPad = 0;
     if (maxpad[0][0] < 0) iPad = 1;
 
+    /*
     // This part of code to take care of edge effect (problems in MC)
     Float_t sprX = fResponse->SigmaIntegration()*fResponse->ChargeSpreadX();
     Float_t sprY = fResponse->SigmaIntegration()*fResponse->ChargeSpreadY();
     Double_t rmin = 9999, rad2;
-    Int_t border = 0, iYlow = 0;
+    Int_t border = 0, iYlow = 0, iMuon = 0;
 
-    if (!fReco) {
+    if (fDraw) {
       for (Int_t i=0; i<2; i++) {
        rad2 = (fXyq[0][maxpad[iPad][0]]-fxyMu[i][0]) * (fXyq[0][maxpad[iPad][0]]-fxyMu[i][0]);
        rad2 += (fXyq[1][maxpad[iPad][0]]-fxyMu[i][1]) * (fXyq[1][maxpad[iPad][0]]-fxyMu[i][1]);
@@ -2956,6 +2440,7 @@ void AliMUONClusterFinderAZ::AddVirtualPad()
        iYlow = fSegmentation[cath]->Iy();
       }
     }
+    */
     
     for (iPad=0; iPad<2; iPad++) {
       if (iPad && !iAddX && !iAddY) break;
@@ -2995,25 +2480,29 @@ void AliMUONClusterFinderAZ::AddVirtualPad()
          if (xList[k] != ix || yList[k] != iy) continue;
          if (!mirror) {
            if ((!cath || maxpad[0][0] < 0) && 
-               (TMath::Abs(iy-iy0) == 1 || TMath::Abs(iy*iy0) == 1)) {
+               //(TMath::Abs(iy-iy0) == 1 || TMath::Abs(iy*iy0) == 1)) {
+               (TMath::Abs(iy-iy0) == 1 || iy*iy0 == -1)) {
              xList[k] = yList[k] = 0; 
              neighb--;
              break;
            }
            if ((cath || maxpad[1][0] < 0) && 
-               (TMath::Abs(ix-ix0) == 1 || TMath::Abs(ix*ix0) == 1)) {
+               //(TMath::Abs(ix-ix0) == 1 || TMath::Abs(ix*ix0) == 1)) {
+               (TMath::Abs(ix-ix0) == 1 || ix*ix0 == -1)) {
              xList[k] = yList[k] = 0; 
              neighb--;
            }
          } else {
            if ((!cath || maxpad[0][0] < 0) && 
-               (TMath::Abs(ix-ix0) == 1 || TMath::Abs(ix*ix0) == 1)) {
+               //(TMath::Abs(ix-ix0) == 1 || TMath::Abs(ix*ix0) == 1)) {
+               (TMath::Abs(ix-ix0) == 1 || ix*ix0 == -1)) {
              xList[k] = yList[k] = 0; 
              neighb--;
              break;
            }
            if ((cath || maxpad[1][0] < 0) && 
-               (TMath::Abs(iy-iy0) == 1 || TMath::Abs(iy*iy0) == 1)) {
+               //(TMath::Abs(iy-iy0) == 1 || TMath::Abs(iy*iy0) == 1)) {
+               (TMath::Abs(iy-iy0) == 1 || iy*iy0 == -1)) {
              xList[k] = yList[k] = 0; 
              neighb--;
            }
@@ -3034,7 +2523,8 @@ void AliMUONClusterFinderAZ::AddVirtualPad()
        fPadIJ[1][npads] = 0;
        ix = xList[j]; 
        iy = yList[j];
-       if (TMath::Abs(ix-ix0) == 1 || TMath::Abs(ix*ix0) == 1) {
+       //if (TMath::Abs(ix-ix0) == 1 || TMath::Abs(ix*ix0) == 1) {
+       if (TMath::Abs(ix-ix0) == 1 || ix*ix0 == -1) {
          if (iy != iy0) continue; // new segmentation - check
          if (nInX != 2) continue; // new
          if (!mirror) {
@@ -3046,6 +2536,7 @@ void AliMUONClusterFinderAZ::AddVirtualPad()
          }
          if (iPad && !iAddX) continue;
          fSegmentation[cath]->GetPadC(fInput->DetElemId(),ix,iy,fXyq[0][npads],fXyq[1][npads],zpad);
+         if (fXyq[0][npads] > 1.e+5) continue; // temporary fix
          if (ix1 == ix0) continue;
          //if (iPad && ix1 == ix0) continue;
          //if (iPad && TMath::Abs(fXyq[0][npads]-fXyq[0][iAddX]) < fXyq[3][iAddX]) continue;
@@ -3115,6 +2606,7 @@ void AliMUONClusterFinderAZ::PadsInXandY(Int_t &nInX, Int_t &nInY)
   nXsaved = nYsaved = 0;
   //if (nInX >= 0) {nInX = nXsaved; nInY = nYsaved; return; }
   Float_t *xPad0 = NULL, *yPad0 = NULL, *xPad1 = NULL, *yPad1 = NULL;
+  Float_t wMinX[2] = {99, 99}, wMinY[2] = {99, 99};
   Int_t *nPad0 = NULL, *nPad1 = NULL;
   Int_t nPads = fnPads[0] + fnPads[1];
   if (fnPads[0]) {
@@ -3136,6 +2628,10 @@ void AliMUONClusterFinderAZ::PadsInXandY(Int_t &nInX, Int_t &nInY)
     else if (nInX > 0 && fPadIJ[1][j] != 1 && fPadIJ[1][j] != -9) continue; // exclude non-marked
     if (nInX <= 0 && fXyq[2][j] > fResponse->Saturation()-1) continue; // skip overflows
     cath = fPadIJ[0][j];
+    if (fXyq[3][j] > 0) { // exclude virtual pads
+      wMinX[cath] = TMath::Min (wMinX[cath], fXyq[3][j]);
+      wMinY[cath] = TMath::Min (wMinY[cath], fXyq[4][j]);
+    }
     if (cath) { xPad1[n1] = fXyq[0][j]; yPad1[n1++] = fXyq[1][j]; }
     else { xPad0[n0] = fXyq[0][j]; yPad0[n0++] = fXyq[1][j]; }
   }
@@ -3160,8 +2656,10 @@ void AliMUONClusterFinderAZ::PadsInXandY(Int_t &nInX, Int_t &nInY)
   }
   if (fnPads[0]) { delete [] xPad0; delete [] yPad0; delete [] nPad0; }
   if (fnPads[1]) { delete [] xPad1; delete [] yPad1; delete [] nPad1; }
-  nInY = TMath::Max (npady[0], npady[1]);
-  nInX = TMath::Max (npadx[0], npadx[1]);
+  //nInY = TMath::Max (npady[0], npady[1]);
+  //nInX = TMath::Max (npadx[0], npadx[1]);
+  nInY = wMinY[0] < wMinY[1] ? npady[0] : npady[1];
+  nInX = wMinX[0] < wMinX[1] ? npadx[0] : npadx[1];
 }
 
 //_____________________________________________________________________________
@@ -3288,7 +2786,7 @@ void AliMUONClusterFinderAZ::Errors(AliMUONRawCluster *clus)
   //errY = 0.01;
   //errX = TMath::Max (errX, 0.144);
   clus->SetX(0, xreco); clus->SetY(0, yreco);
-  clus->SetX(1, errX); clus->SetY(1, errY);
+  clus->SetErrX(errX); clus->SetErrY(errY);
 }
 
 //_____________________________________________________________________________
index d046e1f8409e50923ef91e99ec6581d85536fccc..7730d20315444687d3a9382c84bd5bb5a9cc21b9 100644 (file)
@@ -12,7 +12,6 @@
 
 #include "AliMUONClusterFinderVS.h"
 
-class TH2F;
 class TH2D;
 class TClonesArray;
 class TMinuit;
@@ -21,16 +20,21 @@ class TMatrixD;
 class AliSegmentation;
 class AliMUONResponse;
 class AliMUONPixel;
+class AliMUONClusterDrawAZ;
 
 class AliMUONClusterFinderAZ : public AliMUONClusterFinderVS 
 {
 public:
-  AliMUONClusterFinderAZ(Bool_t draw = 0, Int_t iReco = 1);// Constructor
+  AliMUONClusterFinderAZ(Bool_t draw = 0); // Constructor
   virtual ~AliMUONClusterFinderAZ(); // Destructor
 
   void     FindRawClusters(); // the same interface as for old cluster finder
-  void     EventLoop(Int_t nev, Int_t ch); // first event 
+  void     EventLoop(Int_t nev = 0, Int_t ch = 0); // first event 
   Bool_t   TestTrack(Int_t t) const; // test if track was selected
+  Int_t    GetNPads(Int_t cath) const { return fnPads[cath]; }
+  Int_t    GetIJ(Int_t indx, Int_t iPad) const { return fPadIJ[indx][iPad]; }
+  Float_t  GetXyq(Int_t indx, Int_t iPad) const { return fXyq[indx][iPad]; }
+  Float_t  GetZpad() { return fZpad; }
  
 protected:
   AliMUONClusterFinderAZ(const AliMUONClusterFinderAZ& rhs);
@@ -46,7 +50,6 @@ protected:
   Int_t      fnPads[2];        // ! number of pads in the cluster on 2 cathodes
   Float_t    fXyq[7][fgkDim];    // ! pad information
   Int_t      fPadIJ[2][fgkDim];  // ! pad information
-  //AZ AliSegmentation *fSegmentation[2]; // ! old segmentation
   AliMUONGeometrySegmentation *fSegmentation[2]; // ! new segmentation
   AliMUONResponse *fResponse;// ! response
   Float_t    fZpad;            // ! z-coordinate of the hit
@@ -56,20 +59,14 @@ protected:
 
   static     TMinuit* fgMinuit; // ! Fitter
   Bool_t     fUsed[2][fgkDim]; // ! flags for used pads
-  TH2F*      fHist[4]; // ! histograms
-  TClonesArray *fMuonDigits; // ! pointer to digits
-  Bool_t     fDraw; // ! draw flag
-  Int_t      fnMu; // ! number of muons passing thru the selected area
-  Double_t   fxyMu[2][7]; // ! muon information
+  AliMUONClusterDrawAZ *fDraw; // ! drawing object 
   TObjArray* fPixArray; // ! collection of pixels
   Int_t fnCoupled; // ! number of coupled clusters in precluster
   Int_t fDebug; // ! debug level
 
   // Functions
 
-  void   ModifyHistos(void); // modify histograms
   void   AddPad(Int_t cath, Int_t digit); // add a pad to the cluster
-  //AZ Bool_t Overlap(Int_t cath, TObject *dig); // check if the pad from one cathode overlaps with a pad in the cluster on the other cathode
   Bool_t Overlap(Int_t cath, AliMUONDigit *dig); // check if the pad from one cathode overlaps with a pad in the cluster on the other cathode
   Bool_t Overlap(Float_t *xy1, Int_t iPad, Float_t *xy12, Int_t iSkip); // check if pads xy1 and iPad overlap and return overlap area
   Bool_t CheckPrecluster(Int_t *nShown); // check precluster to simplify it (if possible)
@@ -104,8 +101,6 @@ protected:
              Double_t wy, Double_t wx, Int_t iover, 
              Double_t dyc, Double_t dxc, Double_t qtot, 
              Double_t &yrec, Double_t &xrec, Double_t &erry, Double_t &errx);
-  void DrawCluster(Int_t nev0, Int_t ch0); // draw precluster
-  Int_t Next(Int_t &nev0, Int_t &ch0); // commands for drawing
 
   // Dummy methods for overloading warnings
   void FindCluster(int, int, int, AliMUONRawCluster&) {return;}
index 68b80a943b6ef975fcccc84d5d9a086ee40bd894..8898bb50bd9dfecd96a8b32c8f6c85df7bcf6e12 100644 (file)
 /**************************************************************************
-
  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
-
  *                                                                        *
-
  * Author: The ALICE Off-line Project.                                    *
-
  * Contributors are mentioned in the code where appropriate.              *
-
  *                                                                        *
-
  * Permission to use, copy, modify and distribute this software and its   *
-
  * documentation strictly for non-commercial purposes is hereby granted   *
-
  * without fee, provided that the above copyright notice appears in all   *
-
  * copies and that both the copyright notice and this permission notice   *
-
  * appear in the supporting documentation. The authors make no claims     *
-
  * about the suitability of this software for any purpose. It is          *
-
  * provided "as is" without express or implied warranty.                  *
-
  **************************************************************************/
-
 /* $Id$ */
 
-
-
 ///////////////////////////////////////////////////////////////////////////////
-
 //                                                                           //
-
 // class for MUON reconstruction                                             //
-
 //                                                                           //
-
 ///////////////////////////////////////////////////////////////////////////////
-
 #include <TParticle.h>
-
 #include <TArrayF.h>
 
-
-
 #include "AliRunLoader.h"
-
 #include "AliHeader.h"
-
 #include "AliGenEventHeader.h"
-
 #include "AliESD.h"
-
 #include "AliMUONReconstructor.h"
-
  
-
 #include "AliMUONData.h"
-
 #include "AliMUONTrackReconstructor.h"
-
 #include "AliMUONClusterReconstructor.h"
-
 #include "AliMUONClusterFinderVS.h"
-
 #include "AliMUONClusterFinderAZ.h"
-
 #include "AliMUONTrack.h"
-
 #include "AliMUONTrackParam.h"
-
 #include "AliMUONTriggerTrack.h"
-
 #include "AliESDMuonTrack.h"
-
 #include "AliMUONRawData.h"
 
-
-
 #include "AliRawReader.h"
 
 
-
-
-
 ClassImp(AliMUONReconstructor)
-
 //_____________________________________________________________________________
-
 AliMUONReconstructor::AliMUONReconstructor()
-
   : AliReconstructor()
-
 {
-
 }
-
 //_____________________________________________________________________________
-
 AliMUONReconstructor::~AliMUONReconstructor()
-
 {
-
 }
-
 //_____________________________________________________________________________
-
 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader) const
-
 {
-
 //  AliLoader
-
   AliLoader* loader = runLoader->GetLoader("MUONLoader");
-
   Int_t nEvents = runLoader->GetNumberOfEvents();
 
-
-
 // used local container for each method
-
 // passing fLoader as argument, could be avoided ???
-
   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader);
-
   AliMUONData* dataEvent = recoEvent->GetMUONData();
-
   if (strstr(GetOption(),"Kalman")) recoEvent->SetTrackMethod(2); // Kalman
-
   else recoEvent->SetTrackMethod(1); // original
 
-
-
   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader);
-
   AliMUONData* dataCluster = recoCluster->GetMUONData();
-
   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
-
   if (strstr(GetOption(),"AZ")) {
-
-    recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ(0,1);
-
+    recModel = (AliMUONClusterFinderVS*) new AliMUONClusterFinderAZ();
     recoCluster->SetRecoModel(recModel);
-
   }
-
   recModel->SetGhostChi2Cut(10);
 
-
-
   loader->LoadDigits("READ");
-
   loader->LoadRecPoints("RECREATE");
-
   loader->LoadTracks("RECREATE");
-
   
-
   //   Loop over events              
-
   for(Int_t ievent = 0; ievent < nEvents; ievent++) {
-
     printf("Event %d\n",ievent);
-
     runLoader->GetEvent(ievent);
 
-
-
     //----------------------- digit2cluster & Trigger2Trigger -------------------
-
     if (!loader->TreeR()) loader->MakeRecPointsContainer();
-
      
-
     // tracking branch
-
     dataCluster->MakeBranch("RC");
-
     dataCluster->SetTreeAddress("D,RC");
-
     recoCluster->Digits2Clusters(); 
-
     dataCluster->Fill("RC"); 
 
-
-
     // trigger branch
-
     dataCluster->MakeBranch("TC");
-
     dataCluster->SetTreeAddress("TC");
-
     recoCluster->Trigger2Trigger(); 
-
     dataCluster->Fill("TC");
 
-
-
     loader->WriteRecPoints("OVERWRITE");
 
-
-
     //---------------------------- Track & TriggerTrack ---------------------
-
     if (!loader->TreeT()) loader->MakeTracksContainer();
 
-
-
     // trigger branch
-
     dataEvent->MakeBranch("RL"); //trigger track
-
     dataEvent->SetTreeAddress("RL");
-
     recoEvent->EventReconstructTrigger();
-
     dataEvent->Fill("RL");
 
-
-
     // tracking branch
-
     dataEvent->MakeBranch("RT"); //track
-
     dataEvent->SetTreeAddress("RT");
-
     recoEvent->EventReconstruct();
-
     dataEvent->Fill("RT");
 
-
-
     loader->WriteTracks("OVERWRITE"); 
-
   
-
     //--------------------------- Resetting branches -----------------------
-
     dataCluster->ResetDigits();
-
     dataCluster->ResetRawClusters();
-
     dataCluster->ResetTrigger();
 
-
-
     dataEvent->ResetRawClusters();
-
     dataEvent->ResetTrigger();
-
     dataEvent->ResetRecTracks();  
-
     dataEvent->ResetRecTriggerTracks();
 
-
-
   }
-
   loader->UnloadDigits();
-
   loader->UnloadRecPoints();
-
   loader->UnloadTracks();
 
-
-
   delete recoCluster;
-
   delete recoEvent;
-
 }
 
-
-
 //_____________________________________________________________________________
-
 void AliMUONReconstructor::Reconstruct(AliRunLoader* runLoader, AliRawReader* rawReader) const
-
 {
-
 //  AliLoader
-
   AliLoader* loader = runLoader->GetLoader("MUONLoader");
 
-
-
 // used local container for each method
-
 // passing fLoader as argument, could be avoided ???
-
   AliMUONTrackReconstructor* recoEvent = new AliMUONTrackReconstructor(loader);
-
   AliMUONData* dataEvent = recoEvent->GetMUONData();
 
-
-
   AliMUONRawData* rawData = new AliMUONRawData(loader);
-
   AliMUONData* dataCluster = rawData->GetMUONData();
 
-
-
   AliMUONClusterReconstructor* recoCluster = new AliMUONClusterReconstructor(loader, dataCluster);
-
   AliMUONClusterFinderVS *recModel = recoCluster->GetRecoModel();
-
   recModel->SetGhostChi2Cut(10);
 
-
-
   loader->LoadRecPoints("RECREATE");
-
   loader->LoadTracks("RECREATE");
-
   loader->LoadDigits("RECREATE");
 
 
-
-
-
   //   Loop over events  
-
   Int_t iEvent = 0;
-
             
-
   while (rawReader->NextEvent()) {
-
     printf("Event %d\n",iEvent);
-
     runLoader->GetEvent(iEvent++);
 
-
-
     //----------------------- raw2digits & raw2trigger-------------------
-
     if (!loader->TreeD()) loader->MakeDigitsContainer();
 
-
-
     // tracking branch
-
     dataCluster->MakeBranch("D");
-
     dataCluster->SetTreeAddress("D");
-
     rawData->ReadTrackerDDL(rawReader);
-
     dataCluster->Fill("D"); 
 
-
-
     // trigger branch
-
     dataCluster->MakeBranch("GLT");
-
     dataCluster->SetTreeAddress("GLT");
-
     rawData->ReadTriggerDDL(rawReader);
-
     dataCluster->Fill("GLT"); 
 
-
-
     loader->WriteDigits("OVERWRITE");
 
-
-
     //----------------------- digit2cluster & Trigger2Trigger -------------------
-
     if (!loader->TreeR()) loader->MakeRecPointsContainer();
-
      
-
     // tracking branch
-
     dataCluster->MakeBranch("RC");
-
     dataCluster->SetTreeAddress("RC");
-
     recoCluster->Digits2Clusters(); 
-
     dataCluster->Fill("RC"); 
 
-
-
     // trigger branch
-
     dataCluster->MakeBranch("TC");
-
     dataCluster->SetTreeAddress("TC");
-
     recoCluster->Trigger2Trigger(); 
-
     dataCluster->Fill("TC");
 
-
-
     loader->WriteRecPoints("OVERWRITE");
 
-
-
     //---------------------------- Track & TriggerTrack ---------------------
-
     if (!loader->TreeT()) loader->MakeTracksContainer();
 
-
-
     // trigger branch
-
     dataEvent->MakeBranch("RL"); //trigger track
-
     dataEvent->SetTreeAddress("RL");
-
     recoEvent->EventReconstructTrigger();
-
     dataEvent->Fill("RL");
 
-
-
     // tracking branch
-
     dataEvent->MakeBranch("RT"); //track
-
     dataEvent->SetTreeAddress("RT");
-
     recoEvent->EventReconstruct();
-
     dataEvent->Fill("RT");
 
-
-
     loader->WriteTracks("OVERWRITE");  
-
   
-
     //--------------------------- Resetting branches -----------------------
-
     dataCluster->ResetDigits();
-
     dataCluster->ResetRawClusters();
-
     dataCluster->ResetTrigger();
 
-
-
     dataEvent->ResetRawClusters();
-
     dataEvent->ResetTrigger();
-
     dataEvent->ResetRecTracks();
-
     dataEvent->ResetRecTriggerTracks();
-
   
-
   }
-
   loader->UnloadRecPoints();
-
   loader->UnloadTracks();
-
   loader->UnloadDigits();
 
-
-
   delete recoCluster;
-
   delete recoEvent;
-
 }
 
-
-
 //_____________________________________________________________________________
-
 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliESD* esd) const
-
 {
-
   TClonesArray* recTracksArray = 0;
-
   TClonesArray* recTrigTracksArray = 0;
-
   
-
   AliLoader* loader = runLoader->GetLoader("MUONLoader");
-
   loader->LoadTracks("READ");
-
   AliMUONData* muonData = new AliMUONData(loader,"MUON","MUON");
 
-
-
    // declaration  
-
   Int_t iEvent;// nPart;
-
   Int_t nTrackHits;// nPrimary;
-
   Double_t fitFmin;
-
   TArrayF vertex(3);
 
-
-
   Double_t bendingSlope, nonBendingSlope, inverseBendingMomentum;
-
   Double_t xRec, yRec, zRec, chi2MatchTrigger;
-
   Bool_t matchTrigger;
 
-
-
   // setting pointer for tracks, triggertracks & trackparam at vertex
-
   AliMUONTrack* recTrack = 0;
-
   AliMUONTrackParam* trackParam = 0;
-
   AliMUONTriggerTrack* recTriggerTrack = 0;
-
 //   TParticle* particle = new TParticle();
-
 //   AliGenEventHeader* header = 0;
-
   iEvent = runLoader->GetEventNumber(); 
-
   runLoader->GetEvent(iEvent);
 
-
-
   // vertex calculation (maybe it exists already somewhere else)
-
   vertex[0] = vertex[1] = vertex[2] = 0.;
-
  //  nPrimary = 0;
-
 //   if ( (header = runLoader->GetHeader()->GenEventHeader()) ) {
-
 //     header->PrimaryVertex(vertex);
-
 //   } else {
-
 //     runLoader->LoadKinematics("READ");
-
 //     runLoader->TreeK()->GetBranch("Particles")->SetAddress(&particle);
-
 //     nPart = (Int_t)runLoader->TreeK()->GetEntries();
-
 //     for(Int_t iPart = 0; iPart < nPart; iPart++) {
-
 //       runLoader->TreeK()->GetEvent(iPart);
-
 //       if (particle->GetFirstMother() == -1) {
-
 //     vertex[0] += particle->Vx();
-
 //     vertex[1] += particle->Vy();
-
 //     vertex[2] += particle->Vz();
-
 //     nPrimary++;
-
 //       }
-
 //       if (nPrimary) {
-
 //     vertex[0] /= (double)nPrimary;
-
 //     vertex[1] /= (double)nPrimary;
-
 //     vertex[2] /= (double)nPrimary;
-
 //       }
-
 //     }
-
 //   }
-
   // setting ESD MUON class
-
   AliESDMuonTrack* theESDTrack = new  AliESDMuonTrack() ;
 
-
-
   //-------------------- trigger tracks-------------
-
   Long_t trigPat = 0;
-
   muonData->SetTreeAddress("RL");
-
   muonData->GetRecTriggerTracks();
-
   recTrigTracksArray = muonData->RecTriggerTracks();
 
-
-
   // ready global trigger pattern from first track
-
   if (recTrigTracksArray) 
-
     recTriggerTrack = (AliMUONTriggerTrack*) recTrigTracksArray->First();
-
   if (recTriggerTrack) trigPat = recTriggerTrack->GetGTPattern();
 
-
-
   //printf(">>> Event %d Number of Recconstructed tracks %d \n",iEvent, nrectracks);
-
  
-
   // -------------------- tracks-------------
-
   muonData->SetTreeAddress("RT");
-
   muonData->GetRecTracks();
-
   recTracksArray = muonData->RecTracks();
-
         
-
   Int_t nRecTracks = 0;
-
   if (recTracksArray)
-
     nRecTracks = (Int_t) recTracksArray->GetEntriesFast(); //
-
   
-
   // loop over tracks
-
   for (Int_t iRecTracks = 0; iRecTracks <  nRecTracks;  iRecTracks++) {
 
-
-
     // reading info from tracks
-
     recTrack = (AliMUONTrack*) recTracksArray->At(iRecTracks);
 
-
-
     trackParam = (AliMUONTrackParam*) (recTrack->GetTrackParamAtHit())->First();
-
     trackParam->ExtrapToVertex(vertex[0],vertex[1],vertex[2]);
 
-
-
     bendingSlope            = trackParam->GetBendingSlope();
-
     nonBendingSlope         = trackParam->GetNonBendingSlope();
-
     inverseBendingMomentum = trackParam->GetInverseBendingMomentum();
-
     xRec  = trackParam->GetNonBendingCoor();
-
     yRec  = trackParam->GetBendingCoor();
-
     zRec  = trackParam->GetZ();
 
-
-
     nTrackHits       = recTrack->GetNTrackHits();
-
     fitFmin          = recTrack->GetFitFMin();
-
     matchTrigger     = recTrack->GetMatchTrigger();
-
     chi2MatchTrigger = recTrack->GetChi2MatchTrigger();
 
-
-
     // setting data member of ESD MUON
-
     theESDTrack->SetInverseBendingMomentum(inverseBendingMomentum);
-
     theESDTrack->SetThetaX(TMath::ATan(nonBendingSlope));
-
     theESDTrack->SetThetaY(TMath::ATan(bendingSlope));
-
     theESDTrack->SetZ(zRec);
-
     theESDTrack->SetBendingCoor(yRec); // calculate vertex at ESD or Tracking level ?
-
     theESDTrack->SetNonBendingCoor(xRec);
-
     theESDTrack->SetChi2(fitFmin);
-
     theESDTrack->SetNHit(nTrackHits);
-
     theESDTrack->SetMatchTrigger(matchTrigger);
-
     theESDTrack->SetChi2MatchTrigger(chi2MatchTrigger);
 
-
-
     // storing ESD MUON Track into ESD Event 
-
     if (nRecTracks != 0)  
-
       esd->AddMuonTrack(theESDTrack);
-
   } // end loop tracks
 
-
-
   // add global trigger pattern
-
   if (nRecTracks != 0)  
-
     esd->SetTrigger(trigPat);
 
-
-
   // reset muondata
-
   muonData->ResetRecTracks();
-
   muonData->ResetRecTriggerTracks();
 
-
-
   //} // end loop on event  
-
   loader->UnloadTracks(); 
-
  //  if (!header)
-
 //     runLoader->UnloadKinematics();
-
   delete theESDTrack;
-
   delete muonData;
-
   // delete particle;
-
 }//_____________________________________________________________________________
-
 void AliMUONReconstructor::FillESD(AliRunLoader* runLoader, AliRawReader* /*rawReader*/, AliESD* esd) const
-
 {
-
   // don't need rawReader ???
-
   FillESD(runLoader, esd);
-
 }
-
index ce2df793d36603f05b58c4b65a988c37f53dbaf5..ea99ff40da2733e2213168a637d10e4e5c56df18 100644 (file)
@@ -20,6 +20,7 @@
 #pragma link C++ class AliMUONHitForRec+; 
 #pragma link C++ class AliMUONTrackHit+; 
 #pragma link C++ class AliMUONSegment+;
+#pragma link C++ class AliMUONClusterDrawAZ+;
 
 // raw data
 #pragma link C++ class AliMUONDDLTrigger+;
index 3c94103da68fd1df2b414c8857fb87908194c61b..40d9d6a1fb9f78ff346471b0e677c698e8436246 100644 (file)
@@ -22,7 +22,8 @@ SRCS:= AliMUONClusterReconstructor.cxx \
        AliMUONDDLTracker.cxx \
        AliMUONSubEventTracker.cxx \
        AliMUONRawData.cxx \
-       AliMUONRawStream.cxx 
+       AliMUONRawStream.cxx \
+       AliMUONClusterDrawAZ.cxx
  
 
 HDRS:= $(SRCS:.cxx=.h)