MUON/AliMUONClusterDrawAZ.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // Cluster drawing object for AZ cluster finder
  19
  20 #include <stdlib.h>
  21 #include <Riostream.h>
  22 //#include <TROOT.h>
  23 #include <TCanvas.h>
  24 #include <TLine.h>
  25 //#include <TTree.h>
  26 #include <TH2.h>
  27 #include <TView.h>
  28 #include <TStyle.h>
  29
  30 #include "AliMUONClusterDrawAZ.h"
  31 #include "AliMUONClusterFinderAZ.h"
  32 #include "AliHeader.h"
  33 #include "AliRun.h"
  34 #include "AliMUON.h"
  35 //#include "AliMUONChamber.h"
  36 #include "AliMUONDigit.h"
  37 #include "AliMUONHit.h"
  38 #include "AliMUONRawCluster.h"
  39 //#include "AliMUONClusterInput.h"
  40 #include "AliMUONPixel.h"
  41 //#include "AliMC.h"
  42 #include "AliMUONLoader.h"
  43 #include "AliLog.h"
  44
  45 ClassImp(AliMUONClusterDrawAZ)
  46
  47 //_____________________________________________________________________________
  48 AliMUONClusterDrawAZ::AliMUONClusterDrawAZ()
  49   : TObject()
  50 {
  51 // Default constructor
  52   fFind = NULL; fData = NULL;
  53   for (Int_t i=0; i<4; i++) fHist[i] = NULL;
  54 }
  55
  56 //_____________________________________________________________________________
  57 AliMUONClusterDrawAZ::AliMUONClusterDrawAZ(AliMUONClusterFinderAZ *clusFinder)
  58   : TObject()
  59 {
  60 // Constructor
  61   fFind = clusFinder;
  62   for (Int_t i=0; i<4; i++) fHist[i] = NULL;
  63   fDebug = 1;
  64   fEvent = fChamber = 0;
  65   fModif = 0;
  66   Init();
  67 }
  68
  69 //_____________________________________________________________________________
  70 AliMUONClusterDrawAZ::~AliMUONClusterDrawAZ()
  71 {
  72   // Destructor
  73 }
  74
  75 //_____________________________________________________________________________
  76 AliMUONClusterDrawAZ::AliMUONClusterDrawAZ(const AliMUONClusterDrawAZ& rhs)
  77   : TObject(rhs)
  78 {
  79 // Protected copy constructor
  80
  81   AliFatal("Not implemented.");
  82 }
  83
  84
  85 //_____________________________________________________________________________
  86 AliMUONClusterDrawAZ&
  87 AliMUONClusterDrawAZ::operator=(const AliMUONClusterDrawAZ& rhs)
  88 {
  89 // Protected assignement operator
  90
  91   if (this == &rhs) return *this;
  92
  93   AliFatal("Not implemented.");
  94   return *this;
  95 }
  96
  97 //_____________________________________________________________________________
  98 void AliMUONClusterDrawAZ::Init()
  99 {
 100   // Initialization
 101
 102   TCanvas *c1 = new TCanvas("c1","Clusters",0,0,600,700);
 103   //c1->SetFillColor(10);
 104   c1->Divide(1,2);
 105   new TCanvas("c2","Mlem",700,0,600,350);
 106
 107   // Get pointer to Alice detectors
 108   //AliMUON *muon  = (AliMUON*) gAlice->GetModule("MUON");
 109   //if (!muon) return;
 110   //Loaders
 111   AliRunLoader *rl = AliRunLoader::GetRunLoader();
 112   AliLoader *gime = rl->GetLoader("MUONLoader");
 113   fData = ((AliMUONLoader*)gime)->GetMUONData();
 114
 115   gime->LoadHits("READ");
 116   gime->LoadRecPoints("READ");
 117 }
 118
 119 //_____________________________________________________________________________
 120 Bool_t AliMUONClusterDrawAZ::FindEvCh(Int_t nev, Int_t ch)
 121 {
 122   // Find requested event and chamber (skip the ones before the selected)
 123
 124   if (nev < fEvent) return kFALSE;
 125   else if (nev == fEvent && ch < fChamber) return kFALSE;
 126   fEvent = nev;
 127   fChamber = ch;
 128   return kTRUE;
 129 }
 130
 131 //_____________________________________________________________________________
 132 void AliMUONClusterDrawAZ::DrawCluster()
 133 {
 134   // Draw preclusters
 135
 136   TCanvas *c1 = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c1");
 137
 138   cout << " nev         " << fEvent << endl;
 139
 140   char hName[4];
 141   for (Int_t cath = 0; cath < 2; cath++) {
 142     // Build histograms
 143     if (fHist[cath*2]) {fHist[cath*2]->Delete(); fHist[cath*2] = 0;}
 144     if (fHist[cath*2+1]) {fHist[cath*2+1]->Delete(); fHist[cath*2+1] = 0;}
 145     if (fFind->GetNPads(cath) == 0) continue; // cluster on one cathode only
 146     Float_t wxMin = 999, wxMax = 0, wyMin = 999, wyMax = 0;
 147     Int_t minDx = 0, maxDx = 0, minDy = 0, maxDy = 0;
 148     for (Int_t i = 0; i < fFind->GetNPads(0)+fFind->GetNPads(1); i++) {
 149       if (fFind->GetIJ(0,i) != cath) continue;
 150       if (fFind->GetXyq(3,i) < wxMin) { wxMin = fFind->GetXyq(3,i); minDx = i; }
 151       if (fFind->GetXyq(3,i) > wxMax) { wxMax = fFind->GetXyq(3,i); maxDx = i; }
 152       if (fFind->GetXyq(4,i) < wyMin) { wyMin = fFind->GetXyq(4,i); minDy = i; }
 153       if (fFind->GetXyq(4,i) > wyMax) { wyMax = fFind->GetXyq(4,i); maxDy = i; }
 154     }
 155     cout << minDx << " " << maxDx << " " << minDy << " " << maxDy << endl;
 156     Int_t nx, ny, padSize;
 157     Float_t xmin = 9999, xmax = -9999, ymin = 9999, ymax = -9999;
 158     if (TMath::Nint(fFind->GetXyq(3,minDx)*1000) == TMath::Nint(fFind->GetXyq(3,maxDx)*1000) &&
 159         TMath::Nint(fFind->GetXyq(4,minDy)*1000) == TMath::Nint(fFind->GetXyq(4,maxDy)*1000)) {
 160       // the same segmentation
 161       cout << " Same" << endl;
 162       cout << fFind->GetXyq(3,minDx) << " " << fFind->GetXyq(3,maxDx) << " "
 163            << fFind->GetXyq(4,minDy) << " " << fFind->GetXyq(4,maxDy) << endl;
 164       for (Int_t i = 0; i < fFind->GetNPads(0)+fFind->GetNPads(1); i++) {
 165         if (fFind->GetIJ(0,i) != cath) continue;
 166         if (fFind->GetXyq(0,i) < xmin) xmin = fFind->GetXyq(0,i);
 167         if (fFind->GetXyq(0,i) > xmax) xmax = fFind->GetXyq(0,i);
 168         if (fFind->GetXyq(1,i) < ymin) ymin = fFind->GetXyq(1,i);
 169         if (fFind->GetXyq(1,i) > ymax) ymax = fFind->GetXyq(1,i);
 170       }
 171       xmin -= fFind->GetXyq(3,minDx); xmax += fFind->GetXyq(3,minDx);
 172       ymin -= fFind->GetXyq(4,minDy); ymax += fFind->GetXyq(4,minDy);
 173       nx = TMath::Nint ((xmax-xmin)/wxMin/2);
 174       ny = TMath::Nint ((ymax-ymin)/wyMin/2);
 175       cout << xmin << " " << xmax << " " << nx << " " << ymin << " " << ymax << " " << ny << endl;
 176       sprintf(hName,"h%d",cath*2);
 177       fHist[cath*2] = new TH2D(hName,"cluster",nx,xmin,xmax,ny,ymin,ymax);
 178       for (Int_t i = 0; i < fFind->GetNPads(0)+fFind->GetNPads(1); i++) {
 179         if (fFind->GetIJ(0,i) != cath) continue;
 180         fHist[cath*2]->Fill(fFind->GetXyq(0,i),fFind->GetXyq(1,i),fFind->GetXyq(2,i));
 181         //cout << fXyq[0][i] << fXyq[1][i] << fXyq[2][i] << endl;
 182       }
 183     } else {
 184       // different segmentation in the cluster
 185       cout << " Different" << endl;
 186       cout << fFind->GetXyq(3,minDx) << " " << fFind->GetXyq(3,maxDx) << " "
 187            << fFind->GetXyq(4,minDy) << " " << fFind->GetXyq(4,maxDy) << endl;
 188       Int_t nOK = 0;
 189       Int_t indx, locMin, locMax;
 190       if (TMath::Nint(fFind->GetXyq(3,minDx)*1000) != TMath::Nint(fFind->GetXyq(3,maxDx)*1000)) {
 191         // different segmentation along x
 192         indx = 0;
 193         locMin = minDx;
 194         locMax = maxDx;
 195       } else {
 196         // different segmentation along y
 197         indx = 1;
 198         locMin = minDy;
 199         locMax = maxDy;
 200       }
 201       Int_t loc = locMin;
 202       for (Int_t i = 0; i < 2; i++) {
 203         // loop over different pad sizes
 204         if (i > 0) loc = locMax;
 205         padSize = TMath::Nint(fFind->GetXyq(indx+3,loc)*1000);
 206         xmin = 9999; xmax = -9999; ymin = 9999; ymax = -9999;
 207         for (Int_t j = 0; j < fFind->GetNPads(0)+fFind->GetNPads(1); j++) {
 208           if (fFind->GetIJ(0,j) != cath) continue;
 209           if (TMath::Nint(fFind->GetXyq(indx+3,j)*1000) != padSize) continue;
 210           nOK++;
 211           xmin = TMath::Min (xmin,fFind->GetXyq(0,j));
 212           xmax = TMath::Max (xmax,fFind->GetXyq(0,j));
 213           ymin = TMath::Min (ymin,fFind->GetXyq(1,j));
 214           ymax = TMath::Max (ymax,fFind->GetXyq(1,j));
 215         }
 216         xmin -= fFind->GetXyq(3,loc); xmax += fFind->GetXyq(3,loc);
 217         ymin -= fFind->GetXyq(4,loc); ymax += fFind->GetXyq(4,loc);
 218         nx = TMath::Nint ((xmax-xmin)/fFind->GetXyq(3,loc)/2);
 219         ny = TMath::Nint ((ymax-ymin)/fFind->GetXyq(4,loc)/2);
 220         sprintf(hName,"h%d",cath*2+i);
 221         fHist[cath*2+i] = new TH2D(hName,"cluster",nx,xmin,xmax,ny,ymin,ymax);
 222         for (Int_t j = 0; j < fFind->GetNPads(0)+fFind->GetNPads(1); j++) {
 223           if (fFind->GetIJ(0,j) != cath) continue;
 224           if (TMath::Nint(fFind->GetXyq(indx+3,j)*1000) != padSize) continue;
 225           fHist[cath*2+i]->Fill(fFind->GetXyq(0,j),fFind->GetXyq(1,j),fFind->GetXyq(2,j));
 226         }
 227       } // for (Int_t i=0;
 228       if (nOK != fFind->GetNPads(cath)) cout << " *** Too many segmentations: nPads, nOK "
 229                                              << fFind->GetNPads(cath) << " " << nOK << endl;
 230     } // if (TMath::Nint(fFind->GetXyq(3,minDx)*1000)
 231   } // for (Int_t cath = 0;
 232
 233   // Draw histograms and coordinates
 234   TH2D *hFake = 0x0;
 235   if (fModif) {
 236     // This part works only after a modification of THistPainter::PaintLego(Option_t *)
 237     // in ROOT
 238     Double_t xmin = 9999, ymin = 9999, xmax = -9999, ymax = -9999, aMax = -1;
 239     for (Int_t i = 0; i < 4; i++) {
 240       if (!fHist[i]) continue;
 241       xmin = TMath::Min (xmin, fHist[i]->GetXaxis()->GetXmin());
 242       xmax = TMath::Max (xmax, fHist[i]->GetXaxis()->GetXmax());
 243       ymin = TMath::Min (ymin, fHist[i]->GetYaxis()->GetXmin());
 244       ymax = TMath::Max (ymax, fHist[i]->GetYaxis()->GetXmax());
 245       aMax = TMath:: Max (aMax, fHist[i]->GetMaximum());
 246     }
 247     if (c1->FindObject("hFake")) delete c1->FindObject("hFake");
 248     hFake = new TH2D ("hFake", "", 1, xmin, xmax, 1, ymin, ymax);
 249     hFake->SetMaximum(aMax);
 250     hFake->SetNdivisions(505,"Z");
 251     hFake->SetStats(kFALSE);
 252   }
 253
 254   TObject *stats = 0x0;
 255   for (Int_t cath = 0; cath < 2; cath++) {
 256     if (cath == 0) ModifyHistos();
 257     if (fFind->GetNPads(cath) == 0) continue; // cluster on one cathode only
 258     c1->cd(cath+1);
 259     gPad->SetTheta(55);
 260     gPad->SetPhi(30);
 261     if (fModif) {
 262       if (fHist[cath*2]) {
 263         fHist[cath*2]->Draw();
 264         gPad->Update();
 265         stats = fHist[cath*2]->GetListOfFunctions()->FindObject("stats");
 266       }
 267       hFake->Draw("legoFb");
 268     }
 269
 270     Double_t x, y, x0, y0, r1 = 999, r2 = 0;
 271     if (fHist[cath*2+1]) {
 272       //
 273       x0 = fHist[cath*2]->GetXaxis()->GetXmin() - 1000*TMath::Cos(30*TMath::Pi()/180);
 274       y0 = fHist[cath*2]->GetYaxis()->GetXmin() - 1000*TMath::Sin(30*TMath::Pi()/180);
 275       r1 = 0;
 276       Int_t ihist=cath*2;
 277       for (Int_t iy = 1; iy <= fHist[ihist]->GetNbinsY(); iy++) {
 278         y = fHist[ihist]->GetYaxis()->GetBinCenter(iy)
 279           + fHist[ihist]->GetYaxis()->GetBinWidth(iy);
 280         for (Int_t ix = 1; ix <= fHist[ihist]->GetNbinsX(); ix++) {
 281           if (fHist[ihist]->GetCellContent(ix,iy) > 0.1) {
 282             x = fHist[ihist]->GetXaxis()->GetBinCenter(ix)
 283               + fHist[ihist]->GetXaxis()->GetBinWidth(ix);
 284             r1 = TMath::Max (r1,TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)));
 285           }
 286         }
 287       }
 288       ihist = cath*2 + 1 ;
 289       for (Int_t iy = 1; iy <= fHist[ihist]->GetNbinsY(); iy++) {
 290         y = fHist[ihist]->GetYaxis()->GetBinCenter(iy)
 291           + fHist[ihist]->GetYaxis()->GetBinWidth(iy);
 292         for (Int_t ix = 1; ix <= fHist[ihist]->GetNbinsX(); ix++) {
 293           if (fHist[ihist]->GetCellContent(ix,iy) > 0.1) {
 294             x = fHist[ihist]->GetXaxis()->GetBinCenter(ix)
 295               + fHist[ihist]->GetXaxis()->GetBinWidth(ix);
 296             r2 = TMath::Max (r2,TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)));
 297           }
 298         }
 299       }
 300       cout << r1 << " " << r2 << endl;
 301     } // if (fHist[cath*2+1])
 302     if (r1 > r2) {
 303       if (fModif) fHist[cath*2]->Draw("lego1FbSame");
 304       else fHist[cath*2]->Draw("lego1Fb");
 305       if (fHist[cath*2+1]) fHist[cath*2+1]->Draw("lego1SameAxisBbFb");
 306     } else {
 307       if (fModif) fHist[cath*2+1]->Draw("lego1FbSame");
 308       else fHist[cath*2+1]->Draw("lego1Fb");
 309       fHist[cath*2]->Draw("lego1SameAxisFbBb");
 310     }
 311     c1->Update();
 312     if (fModif) stats->Draw();
 313   } // for (Int_t cath = 0;
 314
 315   // Draw simulated and reconstructed hits
 316   DrawHits();
 317 }
 318
 319 //_____________________________________________________________________________
 320 void AliMUONClusterDrawAZ::DrawHits()
 321 {
 322   // Draw simulated and reconstructed hits
 323
 324   TView *view[2] = { 0x0, 0x0 };
 325   Double_t p1[3]={0}, p2[3], xNDC[6];
 326   TLine *line[99] = {0};
 327   TCanvas *c1 = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c1");
 328   if (c1) {
 329     c1->cd(1);
 330     view[0] = c1->Pad()->GetView();
 331     c1->cd(2);
 332     view[1] = c1->Pad()->GetView();
 333   }
 334   fData->SetTreeAddress("H RC");
 335
 336   TH2D *hist = fHist[0] ? fHist[0] : fHist[2];
 337   p2[2] = hist->GetMaximum();
 338
 339   // Draw simulated hits
 340   cout << " *** Simulated hits *** " << endl;
 341   Int_t ntracks = (Int_t) fData->GetNtracks();
 342   fnMu = 0;
 343   Int_t ix, iy, iok, nLine = 0;
 344   TClonesArray *hits = NULL;
 345   for (Int_t i = 0; i < ntracks; i++) {
 346     fData->GetTrack(i);
 347     hits = fData->Hits();
 348     Int_t nhits = (Int_t) hits->GetEntriesFast();
 349     AliMUONHit* mHit;
 350     for (Int_t ihit = 0; ihit < nhits; ihit++) {
 351       mHit = (AliMUONHit*) hits->UncheckedAt(ihit);
 352       if (mHit->Chamber() != fChamber+1) continue;  // chamber number
 353       if (TMath::Abs(mHit->Z()-fFind->GetZpad()) > 1) continue; // different slat
 354       p2[0] = p1[0] = mHit->X();        // x-pos of hit
 355       p2[1] = p1[1] = mHit->Y();        // y-pos
 356       if (p1[0] < hist->GetXaxis()->GetXmin() ||
 357           p1[0] > hist->GetXaxis()->GetXmax()) continue;
 358       if (p1[1] < hist->GetYaxis()->GetXmin() ||
 359           p1[1] > hist->GetYaxis()->GetXmax()) continue;
 360       // Check if track comes thru pads with signal
 361       iok = 0;
 362       for (Int_t ihist = 0; ihist < 4; ihist++) {
 363         if (!fHist[ihist]) continue;
 364         ix = fHist[ihist]->GetXaxis()->FindBin(p1[0]);
 365         iy = fHist[ihist]->GetYaxis()->FindBin(p1[1]);
 366         if (fHist[ihist]->GetCellContent(ix,iy) > 0.5) {iok = 1; break;}
 367       }
 368       if (!iok) continue;
 369       gStyle->SetLineColor(1);
 370       if (TMath::Abs((Int_t)mHit->Particle()) == 13) {
 371         gStyle->SetLineColor(4);
 372         if (fnMu < 2) {
 373           fxyMu[fnMu][0] = p1[0];
 374           fxyMu[fnMu++][1] = p1[1];
 375         }
 376       }
 377       if (fDebug) printf(" X=%10.4f, Y=%10.4f, Z=%10.4f\n",p1[0],p1[1],mHit->Z());
 378       if (view[0] || view[1]) {
 379         // Take into account track angles
 380         p2[0] += mHit->Tlength() * TMath::Sin(mHit->Theta()/180*TMath::Pi())
 381                                  * TMath::Cos(mHit->Phi()/180*TMath::Pi()) / 2;
 382         p2[1] += mHit->Tlength() * TMath::Sin(mHit->Theta()/180*TMath::Pi())
 383                                  * TMath::Sin(mHit->Phi()/180*TMath::Pi()) / 2;
 384         for (Int_t ipad = 1; ipad < 3; ipad++) {
 385           c1->cd(ipad);
 386           if (!view[ipad-1]) continue;
 387           view[ipad-1]->WCtoNDC(p1, &xNDC[0]);
 388           view[ipad-1]->WCtoNDC(p2, &xNDC[3]);
 389           //c1->DrawLine(xpad[0],xpad[1],xpad[3],xpad[4]);
 390           line[nLine] = new TLine(xNDC[0],xNDC[1],xNDC[3],xNDC[4]);
 391           line[nLine++]->Draw();
 392         }
 393       }
 394     } // for (Int_t ihit = 0; ihit < nhits;
 395   } // for (Int_t i = 0; i < ntracks;
 396
 397   // Draw reconstructed coordinates
 398   if (fData->TreeR()) fData->GetRawClusters();
 399   TClonesArray *rawclust = fData->RawClusters(fChamber);
 400   AliMUONRawCluster *mRaw;
 401   gStyle->SetLineColor(3);
 402   cout << " *** Reconstructed hits *** " << endl;
 403   if (rawclust) {
 404     for (Int_t i = 0; i < rawclust ->GetEntriesFast(); i++) {
 405       mRaw = (AliMUONRawCluster*)rawclust ->UncheckedAt(i);
 406       if (TMath::Abs(mRaw->GetZ(0)-fFind->GetZpad()) > 1) continue; // different slat
 407       p2[0] = p1[0] = mRaw->GetX(0);        // x-pos of hit
 408       p2[1] = p1[1] = mRaw->GetY(0);        // y-pos
 409       if (p1[0] < hist->GetXaxis()->GetXmin() ||
 410           p1[0] > hist->GetXaxis()->GetXmax()) continue;
 411       if (p1[1] < hist->GetYaxis()->GetXmin() ||
 412           p1[1] > hist->GetYaxis()->GetXmax()) continue;
 413       /*
 414         treeD->GetEvent(cath);
 415         cout << mRaw->fMultiplicity[0] << mRaw->fMultiplicity[1] << endl;
 416         for (Int_t j=0; j<mRaw->fMultiplicity[cath]; j++) {
 417         Int_t digit = mRaw->fIndexMap[j][cath];
 418         cout << ((AliMUONDigit*)fMuonDigits->UncheckedAt(digit))->Signal() << endl;
 419         }
 420       */
 421       // Check if track comes thru pads with signal
 422       iok = 0;
 423       for (Int_t ihist = 0; ihist < 4; ihist++) {
 424         if (!fHist[ihist]) continue;
 425         ix = fHist[ihist]->GetXaxis()->FindBin(p1[0]);
 426         iy = fHist[ihist]->GetYaxis()->FindBin(p1[1]);
 427         if (fHist[ihist]->GetCellContent(ix,iy) > 0.5) {iok = 1; break;}
 428       }
 429       if (!iok) continue;
 430       if (fDebug) printf(" X=%10.4f, Y=%10.4f, Z=%10.4f\n",p1[0],p1[1],mRaw->GetZ(0));
 431       if (view[0] || view[1]) {
 432         for (Int_t ipad = 1; ipad < 3; ipad++) {
 433           c1->cd(ipad);
 434           if (!view[ipad-1]) continue;
 435           view[ipad-1]->WCtoNDC(p1, &xNDC[0]);
 436           view[ipad-1]->WCtoNDC(p2, &xNDC[3]);
 437           line[nLine] = new TLine(xNDC[0],xNDC[1],xNDC[3],xNDC[4]);
 438           line[nLine++]->Draw();
 439         }
 440       }
 441     } // for (Int_t i = 0; i < rawclust ->GetEntries();
 442   } // if (rawclust)
 443   c1->Update();
 444 }
 445
 446 //_____________________________________________________________________________
 447 Int_t AliMUONClusterDrawAZ::Next()
 448 {
 449   // What to do next?
 450   // File
 451   FILE *lun = 0;
 452   //lun = fopen("pull.dat","w");
 453
 454   for (Int_t i = 0; i < fnMu; i++) {
 455     // Check again if muon comes thru the used pads (due to extra splitting)
 456     for (Int_t j = 0; j < fFind->GetNPads(0)+fFind->GetNPads(1); j++) {
 457       if (TMath::Abs(fxyMu[i][0]-fFind->GetXyq(0,j))<fFind->GetXyq(3,j) &&
 458           TMath::Abs(fxyMu[i][1]-fFind->GetXyq(1,j))<fFind->GetXyq(4,j)) {
 459         if (fDebug) printf("%12.3e %12.3e %12.3e %12.3e\n",fxyMu[i][2],fxyMu[i][3],fxyMu[i][4],fxyMu[i][5]);
 460         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]);
 461         break;
 462       }
 463     }
 464   } // for (Int_t i=0; i<fnMu;
 465
 466   // What's next?
 467   char command[8];
 468   cout << " What is next? " << endl;
 469   command[0] = ' ';
 470   gets(command);
 471   if (command[0] == 'n' || command[0] == 'N') { fEvent++; fChamber = 0; } // next event
 472   else if (command[0] == 'q' || command[0] == 'Q') { if (lun) fclose(lun); } // exit display
 473   else if (command[0] == 'c' || command[0] == 'C') sscanf(command+1,"%d",&fChamber); // new chamber
 474   else if (command[0] == 'e' || command[0] == 'E') { sscanf(command+1,"%d",&fEvent); fChamber = 0; } // new event
 475   else return 1; // Next precluster
 476   return 0;
 477 }
 478
 479
 480 //_____________________________________________________________________________
 481 void AliMUONClusterDrawAZ::ModifyHistos(void)
 482 {
 483   // Modify histograms to bring them to (approximately) the same size
 484   Int_t nhist = 0;
 485   Float_t hlim[4][4], hbin[4][4]; // first index - xmin, xmax, ymin, ymax
 486
 487   Float_t binMin[4] = {999,999,999,999};
 488
 489   for (Int_t i = 0; i < 4; i++) {
 490     if (!fHist[i]) {
 491       hlim[0][i] = hlim[2][i] = 999;
 492       hlim[1][i] = hlim[3][i] = -999;
 493       continue;
 494     }
 495     hlim[0][i] = fHist[i]->GetXaxis()->GetXmin(); // xmin
 496     hlim[1][i] = fHist[i]->GetXaxis()->GetXmax(); // xmax
 497     hlim[2][i] = fHist[i]->GetYaxis()->GetXmin(); // ymin
 498     hlim[3][i] = fHist[i]->GetYaxis()->GetXmax(); // ymax
 499     hbin[0][i] = hbin[1][i] = fHist[i]->GetXaxis()->GetBinWidth(1);
 500     hbin[2][i] = hbin[3][i] = fHist[i]->GetYaxis()->GetBinWidth(1);
 501     binMin[0] = TMath::Min(binMin[0],hbin[0][i]);
 502     binMin[2] = TMath::Min(binMin[2],hbin[2][i]);
 503     nhist++;
 504   }
 505   binMin[1] = binMin[0];
 506   binMin[3] = binMin[2];
 507   cout << " Nhist: " << nhist << endl;
 508
 509   // Adjust histo limits for cathode with different segmentation
 510   for (Int_t i = 0; i < 4; i+=2) {
 511     if (!fHist[i+1]) continue;
 512     Int_t imin, imax, i1 = i + 1;
 513     for (Int_t lim = 0; lim < 4; lim++) {
 514       while (1) {
 515         if (hlim[lim][i] < hlim[lim][i1]) {
 516           imin = i;
 517           imax = i1;
 518         } else {
 519           imin = i1;
 520           imax = i;
 521         }
 522         if (TMath::Abs(hlim[lim][imin]-hlim[lim][imax])<0.01*binMin[lim]) break;
 523         if (lim == 0 || lim == 2) {
 524           // find lower limit
 525           hlim[lim][imax] -= hbin[lim][imax];
 526         } else {
 527           // find upper limit
 528           hlim[lim][imin] += hbin[lim][imin];
 529         }
 530       } // while (1)
 531     }
 532   }
 533
 534
 535   Int_t imnmx = 0, nExtra = 0;
 536   for (Int_t lim = 0; lim < 4; lim++) {
 537     if (lim == 0 || lim == 2) imnmx = TMath::LocMin(4,hlim[lim]); // find lower limit
 538     else imnmx = TMath::LocMax(4,hlim[lim]); // find upper limit
 539
 540     // Adjust histogram limit
 541     for (Int_t i = 0; i < 4; i++) {
 542       if (!fHist[i]) continue;
 543       nExtra = TMath::Nint ((hlim[lim][imnmx]-hlim[lim][i]) / hbin[lim][i]);
 544       hlim[lim][i] += nExtra * hbin[lim][i];
 545     }
 546   }
 547
 548   // Rebuild histograms
 549   TH2D *hist = 0;
 550   Int_t nx, ny;
 551   Double_t x, y, cont, cmax=0;
 552   char hName[4];
 553   for (Int_t ihist = 0; ihist < 4; ihist++) {
 554     if (!fHist[ihist]) continue;
 555     nx = TMath::Nint((hlim[1][ihist]-hlim[0][ihist])/hbin[0][ihist]);
 556     ny = TMath::Nint((hlim[3][ihist]-hlim[2][ihist])/hbin[2][ihist]);
 557     cout << ihist << " " << hlim[0][ihist] << " " << hlim[1][ihist] << " " << nx;
 558     cout << " " << hlim[2][ihist] << " " << hlim[3][ihist] << " " << ny << endl;
 559     sprintf(hName,"hh%d",ihist);
 560     hist =  new TH2D(hName,"hist",nx,hlim[0][ihist],hlim[1][ihist],ny,hlim[2][ihist],hlim[3][ihist]);
 561     for (Int_t i=1; i<=fHist[ihist]->GetNbinsX(); i++) {
 562       x = fHist[ihist]->GetXaxis()->GetBinCenter(i);
 563       for (Int_t j=1; j<=fHist[ihist]->GetNbinsY(); j++) {
 564         y = fHist[ihist]->GetYaxis()->GetBinCenter(j);
 565         cont = fHist[ihist]->GetCellContent(i,j);
 566         hist->Fill(x,y,cont);
 567       }
 568     }
 569     cmax = TMath::Max (cmax,hist->GetMaximum());
 570     sprintf(hName,"%s%d",fHist[ihist]->GetName(),ihist);
 571     fHist[ihist]->Delete();
 572     fHist[ihist] = new TH2D(*hist);
 573     fHist[ihist]->SetName(hName);
 574     fHist[ihist]->SetNdivisions(505,"Z");
 575     hist->Delete();
 576   }
 577   if (fDebug) printf("%f \n",cmax);
 578
 579   for (Int_t ihist = 0; ihist < 4; ihist++) {
 580     if (!fHist[ihist]) continue;
 581     fHist[ihist]->SetMaximum(cmax);
 582     fHist[ihist]->SetMinimum(0);
 583   }
 584 }
 585
 586 //_____________________________________________________________________________
 587 void AliMUONClusterDrawAZ::AdjustHist(Double_t *xylim, const AliMUONPixel *pixPtr)
 588 {
 589   // Adjust histogram limits for pixel drawing
 590
 591   Float_t xypads[4];
 592   if (fHist[0]) {
 593     xypads[0] = fHist[0]->GetXaxis()->GetXmin();
 594     xypads[1] = -fHist[0]->GetXaxis()->GetXmax();
 595     xypads[2] = fHist[0]->GetYaxis()->GetXmin();
 596     xypads[3] = -fHist[0]->GetYaxis()->GetXmax();
 597     for (Int_t i = 0; i < 4; i++) {
 598       while(1) {
 599         if (xylim[i] < xypads[i]) break;
 600         xylim[i] -= 2*pixPtr->Size(i/2);
 601       }
 602     }
 603   }
 604 }
 605
 606 //_____________________________________________________________________________
 607 void AliMUONClusterDrawAZ::DrawHist(const char* canvas, TH2D *hist)
 608 {
 609   // Draw histogram in given canvas
 610
 611   Int_t ix = 0;
 612   //((TCanvas*)gROOT->FindObject("c2"))->cd();
 613   ((TCanvas*)gROOT->FindObject(canvas))->cd();
 614   gPad->SetTheta(55);
 615   gPad->SetPhi(30);
 616   hist->Draw("lego1Fb");
 617   gPad->Update();
 618   gets((char*)&ix);
 619 }
 620
 621 //_____________________________________________________________________________
 622 void AliMUONClusterDrawAZ::FillMuon(Int_t nfit, const Double_t *parOk, const Double_t *errOk)
 623 {
 624   // Fill muon information
 625
 626   Int_t indx, imax;
 627   Double_t cmax, rad;
 628   for (Int_t i = 0; i < fnMu; i++) {
 629     cmax = fxyMu[i][6];
 630     for (Int_t j = 0; j < nfit; j++) {
 631       indx = j<2 ? j*2 : j*2+1;
 632       rad = (fxyMu[i][0]-parOk[indx])*(fxyMu[i][0]-parOk[indx]) +
 633             (fxyMu[i][1]-parOk[indx+1])*(fxyMu[i][1]-parOk[indx+1]);
 634       if (rad < cmax) {
 635         cmax = rad;
 636         imax = indx;
 637         fxyMu[i][6] = cmax;
 638         fxyMu[i][2] = parOk[imax] - fxyMu[i][0];
 639         fxyMu[i][4] = parOk[imax+1] - fxyMu[i][1];
 640         fxyMu[i][3] = errOk[imax];
 641         fxyMu[i][5] = errOk[imax+1];
 642       }
 643     }
 644   }
 645 }
 646
 647 //_____________________________________________________________________________
 648 void AliMUONClusterDrawAZ::UpdateCluster(Int_t npad)
 649 {
 650   // Update cluster after removing non-overlapped pads
 651
 652   Int_t cath = 0, ix = 0, iy = 0;
 653   cout << " Update cluster " << endl;
 654   gets((char*)&ix);
 655   for (Int_t i = 0; i < npad; i++) {
 656     if (TMath::Nint (fFind->GetXyq(2,i)) != -2) continue;
 657     cath = fFind->GetIJ(0,i);
 658     for (Int_t j = 0; j < 2; j++) {
 659       Int_t ihist = cath * 2 + j;
 660       if (!fHist[ihist]) continue;
 661       ix = fHist[ihist]->GetXaxis()->FindBin(fFind->GetXyq(0,i));
 662       iy = fHist[ihist]->GetYaxis()->FindBin(fFind->GetXyq(1,i));
 663       Double_t cont = fHist[ihist]->GetCellContent(ix, iy);
 664       if (cont < 0.1) continue;
 665       fHist[ihist]->Fill(fFind->GetXyq(0,i), fFind->GetXyq(1,i), -cont);
 666     }
 667   }
 668   TCanvas *c1 = (TCanvas*) gROOT->GetListOfCanvases()->FindObject("c1");
 669   if (c1) {
 670     c1->cd(1);
 671     gPad->Modified();
 672     gPad->Update();
 673     c1->cd(2);
 674     gPad->Modified();
 675     gPad->Update();
 676   }
 677 }