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