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