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