/* $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;
+ fnPads[0]=fnPads[1]=0;
+
+ for (Int_t i=0; i<7; i++)
+ for (Int_t j=0; j<fgkDim; j++)
+ fXyq[i][j]= 9999.;
+
+ for (Int_t i=0; i<2; i++)
+ for (Int_t j=0; j<fgkDim; j++) {
+ fPadIJ[i][j]=-1;
+ fUsed[i][j] = 0;
+ }
+
+ fSegmentation[1] = fSegmentation[0] = 0;
+ fResponse = 0x0;
+
+ fZpad = 100000;
+ fNpar = 0;
+ fQtot = 0;
+ fReco = 1;
+
+ fCathBeg = 0;
+ fPadBeg[0] = fPadBeg[1] = 0;
if (!fgMinuit) fgMinuit = new TMinuit(8);
- fDraw = draw;
- fReco = iReco;
+
+ if (!fgClusterFinder) fgClusterFinder = this;
+ fDraw = 0;
fPixArray = new TObjArray(20);
+ fnCoupled = 0;
fDebug = 0; //0;
- if (draw) fDebug = 1;
- AliWarning("*** Running AZ cluster finder ***");
+
+ if (draw) {
+ fDebug = 1;
+ fReco = 0;
+ fDraw = new AliMUONClusterDrawAZ(this);
+ }
+ cout << " *** Running AZ cluster finder *** " << endl;
}
//_____________________________________________________________________________
{
// 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};
+ if (fReco != 2) { // skip initialization for the combined cluster / track
+ fCathBeg = fPadBeg[0] = fPadBeg[1] = 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 (fReco == 2 && (nShown[0] || nShown[1])) return; // only one precluster for the combined finder
+ 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 = fCathBeg; 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;
Int_t digit;
Bool_t eEOC = kTRUE; // end-of-cluster
- for (digit = 0; digit < ndigits[cath]; digit++) {
+ for (digit = fPadBeg[cath]; digit < ndigits[cath]; digit++) {
mdig = AliMUONClusterInput::Instance()->Digit(cath,digit);
if (first) {
// Find first unused pad
} // 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 isec = fSegmentation[cath]->Sector(fInput->DetElemId(), mdig->PadX(), mdig->PadY());
Int_t nPads = fnPads[0] + fnPads[1];
fXyq[0][nPads] = xpad;
fXyq[1][nPads] = ypad;
fPadIJ[0][nPads] = cath;
fPadIJ[1][nPads] = 0;
fUsed[cath][digit] = kTRUE;
- if (fDebug) printf(" bbb %d %d %f %f %f %f %f %d\n", nPads, cath, xpad, ypad, zpad0, fXyq[3][nPads]*2, fXyq[4][nPads]*2, charge);
+ if (fDebug) printf(" bbb %d %d %f %f %f %f %f %4d %3d %3d\n", nPads, cath, xpad, ypad, zpad0, fXyq[3][nPads]*2, fXyq[4][nPads]*2, charge, mdig->PadX(), mdig->PadY());
fnPads[cath]++;
// Check neighbours
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]--;
}
+ if (fDraw) fDraw->UpdateCluster(npad);
} // if (nFlags > 2)
// Check correlations of cathode charges
if ((over[0] || over[1]) && TMath::Abs(sum[0]-sum[1])/(sum[0]+sum[1])*2 > 1) { // 3 times difference
if (fDebug) cout << " Release " << endl;
// Big difference
- cath = sum[0]>sum[1] ? 0 : 1;
- Int_t imax = 0;
- Double_t cmax=-1;
+ cath = sum[0] > sum[1] ? 0 : 1;
+ Int_t imax = 0, imin = 0;
+ Double_t cmax = -1, cmin = 9999, dxMin = 0, dyMin = 0;
Double_t *dist = new Double_t[npad];
- for (Int_t i=0; i<npad; i++) {
- if (fPadIJ[0][i] != cath) continue;
+ for (Int_t i = 0; i < npad; i++) {
+ if (fPadIJ[0][i] != cath || fXyq[2][i] < 0) continue;
+ if (fXyq[2][i] < cmin) {
+ cmin = fXyq[2][i];
+ imin = i;
+ }
if (fXyq[2][i] < cmax) continue;
cmax = fXyq[2][i];
imax = i;
}
// Arrange pads according to their distance to the max,
// normalized to the pad size
- for (Int_t i=0; i<npad; i++) {
+ for (Int_t i = 0; i < npad; i++) {
dist[i] = 0;
- if (fPadIJ[0][i] != cath) continue;
+ if (fPadIJ[0][i] != cath || fXyq[2][i] < 0) continue;
if (i == imax) continue;
- if (fXyq[2][i] < 0) continue;
- dist[i] = (fXyq[0][i]-fXyq[0][imax])*(fXyq[0][i]-fXyq[0][imax])/
- fXyq[3][imax]/fXyq[3][imax]/4;
- dist[i] += (fXyq[1][i]-fXyq[1][imax])*(fXyq[1][i]-fXyq[1][imax])/
- fXyq[4][imax]/fXyq[4][imax]/4;
- dist[i] = TMath::Sqrt (dist[i]);
+ Double_t dx = (fXyq[0][i] - fXyq[0][imax]) / fXyq[3][imax] / 2;
+ Double_t dy = (fXyq[1][i] - fXyq[1][imax]) / fXyq[4][imax] / 2;
+ dist[i] = TMath::Sqrt (dx * dx + dy * dy);
+ if (i == imin) {
+ cmin = dist[i] + 0.001; // distance to the pad with minimum charge
+ dxMin = dx;
+ dyMin = dy;
+ }
}
TMath::Sort(npad, dist, flags, kFALSE); // in increasing order
Int_t indx;
Double_t xmax = -1;
- for (Int_t i=0; i<npad; i++) {
+ for (Int_t i = 0; i < npad; i++) {
indx = flags[i];
- if (fPadIJ[0][indx] != cath) continue;
- if (fXyq[2][indx] < 0) continue;
- if (fXyq[2][indx] <= cmax || TMath::Abs(dist[indx]-xmax)<1.e-3) {
+ if (fPadIJ[0][indx] != cath || fXyq[2][indx] < 0) continue;
+ if (dist[indx] > cmin) {
+ // Farther than the minimum pad
+ Double_t dx = (fXyq[0][indx] - fXyq[0][imax]) / fXyq[3][imax] / 2;
+ Double_t dy = (fXyq[1][indx] - fXyq[1][imax]) / fXyq[4][imax] / 2;
+ dx *= dxMin;
+ dy *= dyMin;
+ if (dx >= 0 && dy >= 0) continue;
+ if (TMath::Abs(dx) > TMath::Abs(dy) && dx >= 0) continue;
+ if (TMath::Abs(dy) > TMath::Abs(dx) && dy >= 0) continue;
+ }
+ if (fXyq[2][indx] <= cmax || TMath::Abs(dist[indx]-xmax) < 1.e-3) {
// Release pads
- if (TMath::Abs(dist[indx]-xmax)<1.e-3)
+ if (TMath::Abs(dist[indx]-xmax) < 1.e-3)
cmax = TMath::Max((Double_t)(fXyq[2][indx]),cmax);
else cmax = fXyq[2][indx];
xmax = dist[indx];
fUsed[cath][digit] = kFALSE;
fXyq[2][indx] = -2;
fnPads[cath]--;
- // xmax = dist[i]; // Bug?
- } else {
- // Check pad overlaps once more
- for (Int_t i=0; i<npad; i++) flags[i] = 0;
- for (Int_t i=0; i<npad; i++) {
- if (fXyq[2][i] < 0) continue;
- if (fPadIJ[0][i] != i1) continue;
- xy1[0] = fXyq[0][i] - fXyq[3][i];
- xy1[1] = fXyq[0][i] + fXyq[3][i];
- xy1[2] = fXyq[1][i] - fXyq[4][i];
- xy1[3] = fXyq[1][i] + fXyq[4][i];
- for (Int_t j=0; j<npad; j++) {
- if (fXyq[2][j] < 0) continue;
- if (fPadIJ[0][j] != i2) continue;
- if (!Overlap(xy1, j, xy12, 0)) continue;
- flags[i] = flags[j] = 1; // mark overlapped pads
- } // for (Int_t j=0;
- } // for (Int_t i=0;
- nFlags=0;
- for (Int_t i=0; i<npad; i++) {
- if (fXyq[2][i] < 0 || flags[i]) continue;
- nFlags ++;
- }
- if (nFlags == fnPads[0] + fnPads[1]) {
- // No overlap
- for (Int_t i=0; i<npad; i++) {
- if (fXyq[2][i] < 0 || fPadIJ[0][i] != cath) continue;
- fXyq[2][i] = -2;
- fnPads[cath]--;
- }
- }
- break;
+ }
+ } // for (Int_t i = 0; i < npad;
+
+ // Check pad overlaps once more
+ for (Int_t j = 0; j < npad; j++) flags[j] = 0;
+ for (Int_t k = 0; k < npad; k++) {
+ if (fXyq[2][k] < 0 || fPadIJ[0][k] != i1) continue;
+ xy1[0] = fXyq[0][k] - fXyq[3][k];
+ xy1[1] = fXyq[0][k] + fXyq[3][k];
+ xy1[2] = fXyq[1][k] - fXyq[4][k];
+ xy1[3] = fXyq[1][k] + fXyq[4][k];
+ for (Int_t j = 0; j < npad; j++) {
+ if (fXyq[2][j] < 0) continue;
+ if (fPadIJ[0][j] != i2) continue;
+ if (!Overlap(xy1, j, xy12, 0)) continue;
+ flags[k] = flags[j] = 1; // mark overlapped pads
+ } // for (Int_t j = 0;
+ } // for (Int_t k = 0;
+ nFlags = 0;
+ for (Int_t j = 0; j < npad; j++) {
+ if (fXyq[2][j] < 0 || flags[j]) continue;
+ nFlags ++;
+ }
+ if (nFlags == fnPads[0] + fnPads[1]) {
+ // No overlap
+ for (Int_t j = 0; j < npad; j++) {
+ if (fXyq[2][j] < 0 || fPadIJ[0][j] != cath) continue;
+ fXyq[2][j] = -2;
+ fnPads[cath]--;
}
- }
+ }
delete [] dist; dist = 0;
+ if (fDraw) fDraw->UpdateCluster(npad);
} // TMath::Abs(sum[0]-sum[1])...
} // if (fnPads[0] && fnPads[1])
delete [] flags; flags = 0;
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;
}
{
// Check if some pixels have large size (adjust if necessary)
- Int_t nx, ny;
- Int_t nPix = fPixArray->GetEntriesFast();
- AliMUONPixel *pixPtr, *pixPtr1, pix;
+ Int_t n1[2], n2[2], iOK = 1, nPix = fPixArray->GetEntriesFast();
+ AliMUONPixel *pixPtr, pix;
+ Double_t xy0[2] = {9999, 9999}, wxy[2], dist[2];
// Check if large pixel size
- for (Int_t i=0; i<nPix; i++) {
+ for (Int_t i = 0; i < nPix; i++) {
pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
if (pixPtr->Charge() < 1) continue; // discarded pixel
- if (pixPtr->Size(0)-wxmin > 1.e-4 || pixPtr->Size(1)-wymin > 1.e-4) {
- if (fDebug) cout << " Different " << pixPtr->Size(0) << " " << wxmin << " " << pixPtr->Size(1) << " " << wymin << endl;
- pix = *pixPtr;
- nx = TMath::Nint (pix.Size(0)/wxmin);
- ny = TMath::Nint (pix.Size(1)/wymin);
- pix.Shift(0, -pix.Size(0)-wxmin);
- pix.Shift(1, -pix.Size(1)-wymin);
- pix.SetSize(0, wxmin);
- pix.SetSize(1, wymin);
- for (Int_t ii=0; ii<nx; ii++) {
- pix.Shift(0, wxmin*2);
- for (Int_t jj=0; jj<ny; jj++) {
- pix.Shift(1, wymin*2);
- pixPtr1 = new AliMUONPixel(pix);
- fPixArray->Add((TObject*)pixPtr1);
- }
+ if (pixPtr->Size(0) - wxmin < 1.e-4) {
+ if (xy0[0] > 9998) xy0[0] = pixPtr->Coord(0); // position of a "normal" pixel
+ if (pixPtr->Size(1) - wymin < 1.e-4) {
+ if (xy0[1] > 9998) xy0[1] = pixPtr->Coord(1); // position of a "normal" pixel
+ continue;
+ } else iOK = 0; // large pixel
+ } else {
+ iOK = 0; // large pixel
+ if (xy0[1] > 9998 && pixPtr->Size(1) - wymin < 1.e-4) xy0[1] = pixPtr->Coord(1); // "normal" pixel
+ }
+ if (xy0[0] < 9998 && xy0[1] < 9998) break;
+ }
+ if (iOK) return;
+
+ wxy[0] = wxmin;
+ wxy[1] = wymin;
+ //cout << xy0[0] << " " << xy0[1] << endl;
+ for (Int_t i = 0; i < nPix; i++) {
+ pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
+ if (pixPtr->Charge() < 1) continue; // discarded pixel
+ n1[0] = n1[1] = 999;
+ n2[0] = n2[1] = 1;
+ for (Int_t j = 0; j < 2; j++) {
+ if (pixPtr->Size(j) - wxy[j] < 1.e-4) continue;
+ dist[j] = (pixPtr->Coord(j) - xy0[j]) / wxy[j] / 2; // normalized distance to "normal" pixel
+ n2[j] = TMath::Nint (pixPtr->Size(j) / wxy[j]);
+ n1[j] = n2[j] == 1 ? TMath::Nint(dist[j]) : (Int_t)dist[j];
+ }
+ if (n1[0] > 998 && n1[1] > 998) continue;
+ if (fDebug) cout << " Different " << pixPtr->Size(0) << " " << wxy[0] << " "
+ << pixPtr->Size(1) << " " << wxy[1] <<endl;
+
+ if (n2[0] > 2 || n2[1] > 2) { cout << n2[0] << " " << n2[1] << endl; AliFatal("Too large pixel."); }
+ //cout << n1[0] << " " << n2[0] << " " << n1[1] << " " << n2[1] << endl;
+ pix = *pixPtr;
+ pix.SetSize(0, wxy[0]); pix.SetSize(1, wxy[1]);
+ //pixPtr->Print();
+ for (Int_t ii = 0; ii < n2[0]; ii++) {
+ if (n1[0] < 999) pix.SetCoord(0, xy0[0] + (n1[0] + TMath::Sign(1.,dist[0]) * ii) * 2 * wxy[0]);
+ for (Int_t jj = 0; jj < n2[1]; jj++) {
+ if (n1[1] < 999) pix.SetCoord(1, xy0[1] + (n1[1] + TMath::Sign(1.,dist[1]) * jj) * 2 * wxy[1]);
+ fPixArray->Add(new AliMUONPixel(pix));
+ //pix.Print();
}
- pixPtr->SetCharge(0);
}
- } // for (Int_t i=0; i<nPix;
- return;
+ pixPtr->SetCharge(0);
+ } // for (Int_t i = 0; i < nPix;
}
//_____________________________________________________________________________
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) {
Double_t xylim[4] = {999, 999, 999, 999};
for (Int_t ipix=0; ipix<nPix; ipix++) {
pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
+ //cout << ipix+1; pixPtr->Print();
for (Int_t i=0; i<4; i++)
xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
- //cout << ipix+1; pixPtr->Print();
}
for (Int_t i=0; i<4; i++) {
xylim[i] -= pixPtr->Size(i/2); if (fDebug) cout << (i%2 ? -1 : 1)*xylim[i] << " "; }
// 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) {
+ if (!fDraw) {
Double_t coef = 0;
//for (Int_t j=0; j<nfit; j++) {
for (Int_t j=nfit-1; j>=0; j--) {
//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];
- }
- }
- }
+ } else 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;
fXyq[3][npads] = -2; // flag
fnPads[1]++;
iAddX = npads;
- if (fDebug) cout << " ***** Add virtual pad in X ***** " << fXyq[2][npads]
- << " " << fXyq[0][npads] << " " << fXyq[1][npads] << endl;
+ if (fDebug) printf(" ***** Add virtual pad in X ***** %f %f %f %3d %3d \n", fXyq[2][npads],
+ fXyq[0][npads], fXyq[1][npads], ix, iy);
ix1 = ix0;
continue;
}
//if (iPad && TMath::Abs(fXyq[1][npads]-fXyq[1][iAddY]) < fXyq[4][iAddY]) continue;
//if (TMath::Abs(fXyq[0][npads]) < 1 && TMath::Abs(fXyq[1][npads]) < 1) continue; // strange case (something with pad mapping)
if (maxpad[0][0] < 0 || mirror && maxpad[1][0] >= 0) {
- if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[1]/20, 5.);
- else fXyq[2][npads] = TMath::Min (aamax[1]/20, 5.);
+ //if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[1]/20, 5.);
+ //else fXyq[2][npads] = TMath::Min (aamax[1]/20, 5.);
+ if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[1]/15, (double)fResponse->ZeroSuppression());
+ else fXyq[2][npads] = TMath::Min (aamax[1]/15, (double)fResponse->ZeroSuppression());
}
else {
- if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[0]/20, 5.);
- else fXyq[2][npads] = TMath::Min (aamax[0]/20, 5.);
+ //if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[0]/20, 5.);
+ //else fXyq[2][npads] = TMath::Min (aamax[0]/20, 5.);
+ if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[0]/15, (double)fResponse->ZeroSuppression());
+ else fXyq[2][npads] = TMath::Min (aamax[0]/15, (double)fResponse->ZeroSuppression());
}
fXyq[2][npads] = TMath::Max (fXyq[2][npads], (float)1);
//isec = fSegmentation[cath]->Sector(fInput->DetElemId(),ix, iy);
fXyq[3][npads] = -2; // flag
fnPads[1]++;
iAddY = npads;
- if (fDebug) cout << " ***** Add virtual pad in Y ***** " << fXyq[2][npads]
- << " " << fXyq[0][npads] << " " << fXyq[1][npads] << endl;
+ if (fDebug) printf(" ***** Add virtual pad in Y ***** %f %f %f %3d %3d \n", fXyq[2][npads],
+ fXyq[0][npads], fXyq[1][npads], ix, iy);
iy1 = iy0;
}
} // for (Int_t j=0; j<nn;
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]);
+ if (TMath::Abs (wMinY[0] - wMinY[1]) < 1.e-3) nInY = TMath::Max (npady[0], npady[1]);
+ else nInY = wMinY[0] < wMinY[1] ? npady[0] : npady[1];
+ if (TMath::Abs (wMinX[0] - wMinX[1]) < 1.e-3) nInX = TMath::Max (npadx[0], npadx[1]);
+ else 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);
}
//_____________________________________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff