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