#include <TMinuit.h>
#include <TCanvas.h>
#include <TMath.h>
-#include <TROOT.h>
#include "AliCodeTimer.h"
/// \cond CLASSIMP
/// \endcond
const Double_t AliMUONClusterFinderMLEM::fgkZeroSuppression = 6; // average zero suppression value
-const Double_t AliMUONClusterFinderMLEM::fgkSaturation = 3000; // average saturation level
//const Double_t AliMUONClusterFinderMLEM::fgkDistancePrecision = 1e-6; // (cm) used to check overlaps and so on
const Double_t AliMUONClusterFinderMLEM::fgkDistancePrecision = 1e-3; // (cm) used to check overlaps and so on
const TVector2 AliMUONClusterFinderMLEM::fgkIncreaseSize(-AliMUONClusterFinderMLEM::fgkDistancePrecision,-AliMUONClusterFinderMLEM::fgkDistancePrecision);
const TVector2 AliMUONClusterFinderMLEM::fgkDecreaseSize(AliMUONClusterFinderMLEM::fgkDistancePrecision,AliMUONClusterFinderMLEM::fgkDistancePrecision);
- TMinuit* AliMUONClusterFinderMLEM::fgMinuit = 0x0;
+// Status flags for pads
+const Int_t AliMUONClusterFinderMLEM::fgkZero = 0x0; ///< pad "basic" state
+const Int_t AliMUONClusterFinderMLEM::fgkMustKeep = 0x1; ///< do not kill (for pixels)
+const Int_t AliMUONClusterFinderMLEM::fgkUseForFit = 0x10; ///< should be used for fit
+const Int_t AliMUONClusterFinderMLEM::fgkOver = 0x100; ///< processing is over
+const Int_t AliMUONClusterFinderMLEM::fgkModified = 0x1000; ///< modified pad charge
+const Int_t AliMUONClusterFinderMLEM::fgkCoupled = 0x10000; ///< coupled pad
//_____________________________________________________________________________
AliMUONClusterFinderMLEM::AliMUONClusterFinderMLEM(Bool_t plot, AliMUONVClusterFinder* clusterFinder)
fEventNumber(0),
fDetElemId(-1),
fClusterNumber(0),
-fZpad(0.0),
-fReco(1),
-fCathBeg(0),
+fHistMlem(0x0),
+fHistAnode(0x0),
fPixArray(new TObjArray(20)),
fDebug(0),
fPlot(plot),
{
/// Constructor
- fSegmentation[1] = fSegmentation[0] = 0x0;
-
- fPadBeg[0] = fPadBeg[1] = fCathBeg;
-
- if (!fgMinuit) fgMinuit = new TMinuit(8);
+ fkSegmentation[1] = fkSegmentation[0] = 0x0;
if (fPlot) fDebug = 1;
}
AliMUONClusterFinderMLEM::~AliMUONClusterFinderMLEM()
{
/// Destructor
- delete fgMinuit; fgMinuit = 0; delete fPixArray; fPixArray = 0;
+ delete fPixArray; fPixArray = 0;
// delete fDraw;
delete fPreClusterFinder;
delete fSplitter;
//_____________________________________________________________________________
Bool_t
-AliMUONClusterFinderMLEM::Prepare(const AliMpVSegmentation* segmentations[2],
- const AliMUONVDigitStore& digitStore)
+AliMUONClusterFinderMLEM::Prepare(Int_t detElemId,
+ TClonesArray* pads[2],
+ const AliMpArea& area,
+ const AliMpVSegmentation* seg[2])
{
/// Prepare for clustering
-// AliCodeTimerAuto("")
+// AliCodeTimerAuto("",)
for ( Int_t i = 0; i < 2; ++i )
{
- fSegmentation[i] = segmentations[i];
+ fkSegmentation[i] = seg[i];
}
// Find out the DetElemId
- fDetElemId = -1;
-
- TIter next(digitStore.CreateIterator());
- AliMUONVDigit* d = static_cast<AliMUONVDigit*>(next());
-
- if (d)
- {
- fDetElemId = d->DetElemId();
- }
-
- if ( fDetElemId < 0 )
- {
- AliWarning("Could not find DE. Probably no digits at all ?");
- return kFALSE;
- }
+ fDetElemId = detElemId;
delete fSplitter;
fSplitter = new AliMUONClusterSplitterMLEM(fDetElemId,fPixArray);
fSplitter->SetDebug(fDebug);
// find out current event number, and reset the cluster number
- fEventNumber = AliRunLoader::GetRunLoader()->GetEventNumber();
+ AliRunLoader *runLoader = AliRunLoader::Instance();
+ fEventNumber = runLoader ? runLoader->GetEventNumber() : 0;
fClusterNumber = -1;
fClusterList.Delete();
AliDebug(3,Form("EVT %d DE %d",fEventNumber,fDetElemId));
- return fPreClusterFinder->Prepare(segmentations,digitStore);
+ if ( fPreClusterFinder->NeedSegmentation() )
+ {
+ return fPreClusterFinder->Prepare(detElemId,pads,area,seg);
+ }
+ else
+ {
+ return fPreClusterFinder->Prepare(detElemId,pads,area);
+ }
}
//_____________________________________________________________________________
AliMUONClusterFinderMLEM::NextCluster()
{
/// Return next cluster
-// AliCodeTimerAuto("")
+// AliCodeTimerAuto("",)
// if the list of clusters is not void, pick one from there
TObject* o = fClusterList.At(++fClusterNumber);
// pre-cluster and treat it
fPreCluster = fPreClusterFinder->NextCluster();
-
+
+ fClusterList.Delete(); // reset the list of clusters for this pre-cluster
+ fClusterNumber = -1; //AZ
+
if (!fPreCluster)
{
// we are done
return 0x0;
}
- fClusterList.Delete(); // reset the list of clusters for this pre-cluster
- fClusterNumber = -1; //AZ
-
WorkOnPreCluster();
// WorkOnPreCluster may have used only part of the pads, so we check that
/// Starting from a precluster, builds a pixel array, and then
/// extract clusters from this array
- // AliCodeTimerAuto("")
+ // AliCodeTimerAuto("",)
if (fDebug) {
- cout << " *** Event # " << AliRunLoader::GetRunLoader()->GetEventNumber()
+ cout << " *** Event # " << fEventNumber
<< " det. elem.: " << fDetElemId << endl;
for (Int_t j = 0; j < fPreCluster->Multiplicity(); ++j) {
AliMUONPad* pad = fPreCluster->Pad(j);
Int_t nMax = 1, localMax[100], maxPos[100];
Double_t maxVal[100];
- if (cluster->Multiplicity() > 50)
- {
- nMax = FindLocalMaxima(fPixArray, localMax, maxVal);
- }
- //nMax = 1; // just for test
-
- if (nMax > 1)
- {
- TMath::Sort(nMax, maxVal, maxPos, kTRUE); // in decreasing order
- }
-
Int_t iSimple = 0, nInX = -1, nInY;
PadsInXandY(*cluster,nInX, nInY);
- if (nMax == 1 && nInX < 4 && nInY < 4)
+ if (nInX < 4 && nInY < 4)
{
- iSimple = 1; //1; // simple cluster
+ iSimple = 1; // simple cluster
+ }
+ else
+ {
+ nMax = FindLocalMaxima(fPixArray, localMax, maxVal); // for small clusters just to tag pixels
+ if (nMax > 1) {
+ if (cluster->Multiplicity() <= 50) nMax = 1; // for small clusters
+ if (nMax > 1) TMath::Sort(nMax, maxVal, maxPos, kTRUE); // in descending order
+ }
}
for (Int_t i = 0; i < nMax; ++i)
for (Int_t j = 0; j < mult; ++j)
{
AliMUONPad* pad = cluster->Pad(j);
- if ( pad->Status() == 0 ) continue; // pad charge was not modified
- pad->SetStatus(0);
+ //if ( pad->Status() == 0 ) continue; // pad charge was not modified
+ if ( pad->Status() != fgkOver ) continue; // pad was not used
+ //pad->SetStatus(0);
+ pad->SetStatus(fgkZero);
pad->RevertCharge(); // use backup charge value
}
}
} // for (Int_t i=0; i<nMax;
- if (nMax > 1) ((TH2D*) gROOT->FindObject("anode"))->Delete();
- //TH2D *mlem = (TH2D*) gROOT->FindObject("mlem");
- //if (mlem) mlem->Delete();
+ delete fHistMlem;
+ delete fHistAnode;
+ fHistMlem = fHistAnode = 0x0;
delete cluster;
return kTRUE;
}
/// Check precluster in order to attempt to simplify it (mostly for
/// two-cathode preclusters)
- // AliCodeTimerAuto("")
+ AliCodeTimerAuto("",)
// Disregard small clusters (leftovers from splitting or noise)
if ((origCluster.Multiplicity()==1 || origCluster.Multiplicity()==2) &&
return 0x0;
}
- AliMUONCluster* cluster = static_cast<AliMUONCluster*>(origCluster.Clone());
+ AliMUONCluster* cluster = new AliMUONCluster(origCluster);
AliDebug(2,"Start of CheckPreCluster=");
//StdoutToAliDebug(2,cluster->Print("full"));
{
/// Check two-cathode cluster
- Int_t i1 = cluster->Multiplicity(0) ? 0 : 1;
- Int_t i2 = cluster->Multiplicity(1) ? 1 : 0;
-
Int_t npad = cluster->Multiplicity();
Int_t* flags = new Int_t[npad];
for (Int_t j = 0; j < npad; ++j) flags[j] = 0;
for ( Int_t i = 0; i < npad; ++i)
{
AliMUONPad* padi = cluster->Pad(i);
- if ( padi->Cathode() != i1 ) continue;
+ if ( padi->Cathode() != 0 ) continue;
for (Int_t j = i+1; j < npad; ++j)
{
AliMUONPad* padj = cluster->Pad(j);
- if ( padj->Cathode() != i2 ) continue;
+ if ( padj->Cathode() != 1 ) continue;
if ( !AliMUONPad::AreOverlapping(*padi,*padj,fgkDecreaseSize) ) continue;
flags[i] = flags[j] = 1; // mark overlapped pads
}
Int_t cath = pad->Cathode();
Int_t cath1 = TMath::Even(cath);
// Check for edge effect (missing pads on the _other_ cathode)
- AliMpPad mpPad = fSegmentation[cath1]->PadByPosition(pad->Position(),kFALSE);
+ AliMpPad mpPad =
+ fkSegmentation[cath1]->PadByPosition(pad->Position().X(),
+ pad->Position().Y(),kFALSE);
if (!mpPad.IsValid()) continue;
- if (nFlags == 1 && pad->Charge() < fgkZeroSuppression * 3) continue;
+ //if (nFlags == 1 && pad->Charge() < fgkZeroSuppression * 3) continue;
+ if (nFlags == 1 && pad->Charge() < 20) continue;
AliDebug(2,Form("Releasing the following pad : de,cath,ix,iy %d,%d,%d,%d charge %e",
fDetElemId,pad->Cathode(),pad->Ix(),pad->Iy(),pad->Charge()));
toBeRemoved.AddLast(pad);
{
// big difference
Int_t cathode = cluster->MaxRawChargeCathode();
- Int_t imin(0);
- Int_t imax(0);
+ Int_t imin(-1);
+ Int_t imax(-1);
Double_t cmax(0);
Double_t cmin(1E9);
{
cmin = pad->Charge();
imin = i;
+ if (imax < 0) {
+ imax = imin;
+ cmax = cmin;
+ }
}
- if ( pad->Charge() > cmax )
+ else if ( pad->Charge() > cmax )
{
cmax = pad->Charge();
imax = i;
cout << " Surface = " << pixelI->Size(0)*pixelI->Size(1)*4 << endl;
pixelJ->Print();
cout << " Surface = " << pixelJ->Size(0)*pixelJ->Size(1)*4 << endl;
- cout << " Area surface = " << area.Dimensions().X()*area.Dimensions().Y()*4 << endl;
+ cout << " Area surface = " << area.GetDimensionX()*area.GetDimensionY()*4 << endl;
cout << "-------" << endl;
);
*/
}
fPixArray->Delete();
-
- if ( !cluster.Multiplicity(0) || !cluster.Multiplicity(1) )
- {
- BuildPixArrayOneCathode(cluster);
- }
- else
- {
- //BuildPixArrayTwoCathodes(cluster);
- BuildPixArrayOneCathode(cluster);
- }
-
- //fPixArray->Sort(); // FIXME : not really needed, only to compare with ClusterFinderAZ
+ BuildPixArrayOneCathode(cluster);
Int_t nPix = fPixArray->GetLast()+1;
// AliDebug(2,Form("nPix after BuildPixArray=%d",nPix));
- Double_t xPadMin(1E9);
- Double_t yPadMin(1E9);
-
- //for ( Int_t i = 0; i < cluster.Multiplicity(); ++i )
- for ( Int_t i = 0; i < npad; ++i )
- {
- AliMUONPad* pad = cluster.Pad(i);
- xPadMin = TMath::Min (xPadMin, pad->DX());
- yPadMin = TMath::Min (yPadMin, pad->DY());
- }
-
- Double_t wxmin(1E9);
- Double_t wymin(1E9);
-
- for ( Int_t i = 0; i < nPix; ++i )
- {
- AliMUONPad* pixPtr = Pixel(i);
- wxmin = TMath::Min(wxmin, pixPtr->Size(0));
- wymin = TMath::Min(wymin, pixPtr->Size(1));
- }
-
- //wxmin = TMath::Abs (wxmin - xPadMin/2) > 0.001 ? xPadMin : xPadMin / 2;
- //wymin = TMath::Abs (wymin - yPadMin/2) > 0.001 ? yPadMin : yPadMin / 2;
- wxmin = xPadMin;
- wymin = yPadMin;
-
- /*
- // Check if small pixel X-size
- AdjustPixel(cluster,wxmin, 0);
- // Check if small pixel Y-size
- AdjustPixel(cluster,wymin, 1);
- // Check if large pixel size
- AdjustPixel(wxmin, wymin);
- */
-
- // Remove discarded pixels
- for (Int_t i = 0; i < nPix; ++i)
- {
- AliMUONPad* pixPtr = Pixel(i);
- if (pixPtr->Charge() < 1)
- {
- AliDebug(2,Form("Removing pixel %d with charge<1 : ",i));
- //StdoutToAliDebug(2,pixPtr->Print());
- RemovePixel(i);
- }
- }
-
- fPixArray->Compress();
- nPix = fPixArray->GetEntriesFast();
-
-// AliDebug(2,Form("nPix after AdjustPixel=%d",nPix));
-
- //if ( nPix > cluster.Multiplicity() )
if ( nPix > npad )
{
// AliDebug(2,Form("Will trim number of pixels to number of pads"));
// Too many pixels - sort and remove pixels with the lowest signal
fPixArray->Sort();
- for ( Int_t i = cluster.Multiplicity(); i < nPix; ++i )
+ for ( Int_t i = npad; i < nPix; ++i )
{
RemovePixel(i);
}
fPixArray->Compress();
- nPix = fPixArray->GetEntriesFast();
} // if (nPix > npad)
// StdoutToAliDebug(2,cout << "End of BuildPixelArray:" << endl;
// AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
- // Find min and max cluster dimensions
- Double_t minx[2] = {9999,9999}, maxx[2] = {-9999,-9999};
- Double_t miny[2] = {9999,9999}, maxy[2] = {-9999,-9999};
-
TVector2 dim = cluster.MinPadDimensions (-1, kFALSE);
- Double_t width[2] = {dim.X(), dim.Y()}, xy0[2];
- Int_t found[2] = {0};
+ Double_t width[2] = {dim.X(), dim.Y()}, xy0[2]={99999,99999};
+ Int_t found[2] = {0,0}, mult = cluster.Multiplicity();
- for ( Int_t i = 0; i < cluster.Multiplicity(); ++i) {
+ for ( Int_t i = 0; i < mult; ++i) {
AliMUONPad* pad = cluster.Pad(i);
- Int_t cath = pad->Cathode();
- minx[cath] = TMath::Min (minx[cath], pad->Coord(0)-pad->Size(0));
- maxx[cath] = TMath::Max (maxx[cath], pad->Coord(0)+pad->Size(0));
- miny[cath] = TMath::Min (miny[cath], pad->Coord(1)-pad->Size(1));
- maxy[cath] = TMath::Max (maxy[cath], pad->Coord(1)+pad->Size(1));
for (Int_t j = 0; j < 2; ++j) {
if (found[j] == 0 && TMath::Abs(pad->Size(j)-width[j]) < fgkDistancePrecision) {
xy0[j] = pad->Coord(j);
found[j] = 1;
}
}
+ if (found[0] && found[1]) break;
}
- /*
- TVector2 leftDown = cluster.Area(0).LeftDownCorner();
- TVector2 rightUp = cluster.Area(0).RightUpCorner();
- cout << leftDown.X() << " " << leftDown.Y() << " " << rightUp.X() << " " << rightUp.Y() << endl;
- leftDown = cluster.Area(1).LeftDownCorner();
- rightUp = cluster.Area(1).RightUpCorner();
- cout << leftDown.X() << " " << leftDown.Y() << " " << rightUp.X() << " " << rightUp.Y() << endl;
- */
- //cout << minx[0] << " " << maxx[0] << " " << minx[1] << " " << maxx[1] << endl;
- //cout << miny[0] << " " << maxy[0] << " " << miny[1] << " " << maxy[1] << endl;
- //cout << width[0] << " " << width[1] << endl;
Double_t min[2], max[2];
Int_t cath0 = 0, cath1 = 1;
if (cluster.Multiplicity(0) == 0) cath0 = 1;
else if (cluster.Multiplicity(1) == 0) cath1 = 0;
- min[0] = TMath::Max (minx[cath0], minx[cath1]);
- min[1] = TMath::Max (miny[cath0], miny[cath1]);
- max[0] = TMath::Min (maxx[cath0], maxx[cath1]);
- max[1] = TMath::Min (maxy[cath0], maxy[cath1]);
+
+
+ Double_t leftDownX, leftDownY;
+ cluster.Area(cath0).LeftDownCorner(leftDownX, leftDownY);
+ Double_t rightUpX, rightUpY;
+ cluster.Area(cath0).RightUpCorner(rightUpX, rightUpY);
+ min[0] = leftDownX;
+ min[1] = leftDownY;
+ max[0] = rightUpX;
+ max[1] = rightUpY;;
+ if (cath1 != cath0) {
+ cluster.Area(cath1).LeftDownCorner(leftDownX, leftDownY);
+ cluster.Area(cath1).RightUpCorner(rightUpX, rightUpY);
+ min[0] = TMath::Max (min[0], leftDownX);
+ min[1] = TMath::Max (min[1], leftDownY);
+ max[0] = TMath::Min (max[0], rightUpX);
+ max[1] = TMath::Min (max[1], rightUpY);
+ }
// Adjust limits
//width[0] /= 2; width[1] /= 2; // just for check
- Int_t nbins[2];
+ Int_t nbins[2]={0,0};
for (Int_t i = 0; i < 2; ++i) {
Double_t dist = (min[i] - xy0[i]) / width[i] / 2;
+ if (TMath::Abs(dist) < 1.e-6) dist = -1.e-6;
min[i] = xy0[i] + (TMath::Nint(dist-TMath::Sign(1.e-6,dist))
+ TMath::Sign(0.5,dist)) * width[i] * 2;
nbins[i] = TMath::Nint ((max[i] - min[i]) / width[i] / 2);
TAxis *yaxis = hist1->GetYaxis();
// Fill histogram
- Int_t mult = cluster.Multiplicity();
for ( Int_t i = 0; i < mult; ++i) {
AliMUONPad* pad = cluster.Pad(i);
Int_t ix0 = xaxis->FindBin(pad->X());
Int_t iy0 = yaxis->FindBin(pad->Y());
- PadOverHist(0, ix0, iy0, pad);
+ PadOverHist(0, ix0, iy0, pad, hist1, hist2);
}
// Store pixels
Double_t x = xaxis->GetBinCenter(i);
for (Int_t j = 1; j <= nbins[1]; ++j) {
if (hist2->GetCellContent(i,j) < 0.1) continue;
- if (hist2->GetCellContent(i,j) < 1.1 && cluster.Multiplicity(0) &&
- cluster.Multiplicity(1)) continue;
+ //if (hist2->GetCellContent(i,j) < 1.1 && cluster.Multiplicity(0) &&
+ // cluster.Multiplicity(1)) continue;
+ if (cath0 != cath1) {
+ // Two-sided cluster
+ Double_t cont = hist2->GetCellContent(i,j);
+ if (cont < 999.) continue;
+ if (cont-Int_t(cont/1000.)*1000. < 0.5) continue;
+ }
Double_t y = yaxis->GetBinCenter(j);
Double_t charge = hist1->GetCellContent(i,j);
AliMUONPad* pixPtr = new AliMUONPad(x, y, width[0], width[1], charge);
fPixArray->Add(pixPtr);
}
}
+ //*
+ if (fPixArray->GetEntriesFast() == 1) {
+ // Split pixel into 2
+ AliMUONPad* pixPtr = static_cast<AliMUONPad*> (fPixArray->UncheckedAt(0));
+ pixPtr->SetSize(0,width[0]/2.);
+ pixPtr->Shift(0,-width[0]/4.);
+ pixPtr = new AliMUONPad(pixPtr->X()+width[0], pixPtr->Y(), width[0]/2., width[1], pixPtr->Charge());
+ fPixArray->Add(pixPtr);
+ }
+ //*/
//fPixArray->Print();
delete hist1;
delete hist2;
}
//_____________________________________________________________________________
-void AliMUONClusterFinderMLEM::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, AliMUONPad *pad)
+void AliMUONClusterFinderMLEM::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, AliMUONPad *pad,
+ TH2D *hist1, TH2D *hist2)
{
/// "Span" pad over histogram in the direction idir
- TH2D *hist1 = static_cast<TH2D*> (gROOT->FindObject("Grid"));
- TH2D *hist2 = static_cast<TH2D*> (gROOT->FindObject("Entries"));
TAxis *axis = idir == 0 ? hist1->GetXaxis() : hist1->GetYaxis();
- Int_t nbins = axis->GetNbins();
- Double_t bin = axis->GetBinWidth(1);
+ Int_t nbins = axis->GetNbins(), cath = pad->Cathode();
+ Double_t bin = axis->GetBinWidth(1), amask = TMath::Power(1000.,cath*1.);
Int_t nbinPad = (Int_t)(pad->Size(idir)/bin*2+fgkDistancePrecision) + 1; // number of bins covered by pad
if (ixy > nbins) break;
Double_t lowEdge = axis->GetBinLowEdge(ixy);
if (lowEdge + fgkDistancePrecision > pad->Coord(idir) + pad->Size(idir)) break;
- if (idir == 0) PadOverHist(1, ixy, iy0, pad); // span in the other direction
+ if (idir == 0) PadOverHist(1, ixy, iy0, pad, hist1, hist2); // span in the other direction
else {
// Fill histogram
Double_t cont = pad->Charge();
if (hist2->GetCellContent(ix0, ixy) > 0.1)
cont = TMath::Min (hist1->GetCellContent(ix0, ixy), cont);
hist1->SetCellContent(ix0, ixy, cont);
- hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+1);
+ //hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+1);
+ hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+amask);
}
}
if (ixy < 1) break;
Double_t upEdge = axis->GetBinUpEdge(ixy);
if (upEdge - fgkDistancePrecision < pad->Coord(idir) - pad->Size(idir)) break;
- if (idir == 0) PadOverHist(1, ixy, iy0, pad); // span in the other direction
+ if (idir == 0) PadOverHist(1, ixy, iy0, pad, hist1, hist2); // span in the other direction
else {
// Fill histogram
Double_t cont = pad->Charge();
if (hist2->GetCellContent(ix0, ixy) > 0.1)
cont = TMath::Min (hist1->GetCellContent(ix0, ixy), cont);
hist1->SetCellContent(ix0, ixy, cont);
- hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+1);
+ //hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+1);
+ hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+amask);
}
}
}
-/*
-//_____________________________________________________________________________
-void AliMUONClusterFinderMLEM::BuildPixArrayOneCathode(AliMUONCluster& cluster)
-{
- /// From a single-cathode cluster, build the pixel array
-
-// AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
-
- for ( Int_t j=0; j<cluster.Multiplicity(); ++j)
- {
- AliMUONPad* pad = cluster.Pad(j);
- AliMUONPad* pixPtr = new AliMUONPad(pad->Position(),pad->Dimensions(),
- pad->Charge());
- fPixArray->Add(pixPtr);
- }
-}
-*/
-
-//_____________________________________________________________________________
-void AliMUONClusterFinderMLEM::BuildPixArrayTwoCathodes(AliMUONCluster& cluster)
-{
- /// From a two-cathodes cluster, build the pixel array
-
-// AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
-
- Int_t i1 = cluster.Pad(0)->Cathode();
- Int_t i2 = TMath::Even(i1);
-
- for ( Int_t i = 0; i < cluster.Multiplicity(); ++i)
- {
- AliMUONPad* padi = cluster.Pad(i);
- if (padi->Cathode() != i1) continue;
-
- for ( Int_t j = 1; j < cluster.Multiplicity(); ++j)
- {
- AliMUONPad* padj = cluster.Pad(j);
- if (padj->Cathode() != i2) continue;
-
- AliMpArea overlap;
-
- if ( AliMUONPad::AreOverlapping(*padi,*padj,fgkDecreaseSize,overlap) )
- {
- AliMUONPad* pixPtr = new AliMUONPad(overlap.Position(),
- overlap.Dimensions(),
- TMath::Min(padi->Charge(),padj->Charge()));
- if ( ( padi->Charge() <= padj->Charge() && padi->IsSaturated() ) ||
- ( padj->Charge() < padi->Charge() && padj->IsSaturated() ) )
- {
- // if smallest charge of the 2 pads is already above saturation, then
- // the pixel is saturated...
- pixPtr->SetSaturated(kTRUE);
- }
- pixPtr->SetReal(kFALSE);
- fPixArray->Add(pixPtr);
- }
- }
- }
-}
-
-//_____________________________________________________________________________
-void AliMUONClusterFinderMLEM::AdjustPixel(AliMUONCluster& /*cluster*/,
- Float_t width, Int_t ixy)
-{
- /// Check if some pixels have smaller size than others (adjust if necessary)
-
- AliDebug(2,Form("width=%e ixy=%d",width,ixy));
-
- AliMUONPad *pixPtr, *pixPtr1 = 0;
- Int_t ixy1 = !ixy;
- Int_t nPix = fPixArray->GetEntriesFast(), iOK = 1;
-
- Double_t xy0 = 0, minmax[2] = {9999,-9999}, dist = 0;
- // First, find a "normal" pixel
- for (Int_t i = 0; i < nPix; ++i) {
- pixPtr = Pixel(i);
- if (pixPtr->Charge() < 1) continue; // discarded pixel
- minmax[0] = TMath::Min (minmax[0], pixPtr->Size(ixy));
- minmax[1] = TMath::Max (minmax[1], pixPtr->Size(ixy));
- if (pixPtr->Size(ixy) - width < -fgkDistancePrecision) iOK = 0;
- if (TMath::Abs(pixPtr->Size(ixy)-width) > fgkDistancePrecision) continue;
- xy0 = pixPtr->Coord(ixy);
- }
- if (TMath::Abs(minmax[0]-minmax[1]) < fgkDistancePrecision) iOK = 1; // the same size
- if (iOK == 1) return; // all pixels have the same size in the direction IXY
-
- //cout << " --- " << xy0 << endl; fPixArray->Print();
- for (Int_t i = 0; i < nPix; ++i)
- {
- pixPtr = Pixel(i);
- if (pixPtr->Charge() < 1) continue; // discarded pixel
- if (pixPtr->Size(ixy) - width < -fgkDistancePrecision)
- {
- // try to merge
- if (fDebug) cout << i << " Small X or Y: " << ixy << " " << pixPtr->Size(ixy) << " " << width << " " << pixPtr->Coord(0) << " " << pixPtr->Coord(1) << endl;
- for (Int_t j = i + 1; j < nPix; ++j)
- {
- pixPtr1 = Pixel(j);
- if (pixPtr1->Charge() < 1) continue; // discarded pixel
- if (TMath::Abs(pixPtr1->Size(ixy)-width) < fgkDistancePrecision) continue; // right size
- if (TMath::Abs(pixPtr1->Coord(ixy1)-pixPtr->Coord(ixy1)) > fgkDistancePrecision) continue; // different rows/columns
- dist = TMath::Abs (pixPtr1->Coord(ixy) - pixPtr->Coord(ixy));
- if (TMath::Abs(dist-pixPtr1->Size(ixy)-pixPtr->Size(ixy)) < fgkDistancePrecision) // neighbours
- {
- // merge
- //Double_t dist = (pixPtr->Coord(ixy) + pixPtr1->Coord(ixy)) / 2;
- //dist = TMath::Nint((dist-xy0)/width/2) * width * 2;
- dist = (pixPtr->Coord(ixy)-xy0) / width / 2;
- cout << j << " " << dist << endl;
- dist = TMath::Nint(dist) * width * 2;
- pixPtr->SetCoord(ixy, xy0+dist);
- pixPtr->SetSize(ixy, width);
- pixPtr->SetCharge(TMath::Min (pixPtr->Charge(),pixPtr1->Charge()));
- pixPtr1->SetCharge(0);
- pixPtr1 = 0x0;
- break;
- }
- } // for (Int_t j = i + 1;
- if (pixPtr1 || i == nPix-1) {
- // edge pixel - just increase its size
- if (fDebug) cout << " No pair ..." << endl;
- cout << (pixPtr->Coord(ixy)-xy0)/width/2 << endl;
- dist = (pixPtr->Coord(ixy) - xy0) / width / 2;
- dist = TMath::Nint(dist) * width * 2;
- pixPtr->SetCoord(ixy, xy0+dist);
- pixPtr->SetSize(ixy, width);
- }
- } // if (pixPtr->Size(ixy)-width < -fgkDistancePrecision)
- } // for (Int_t i = 0; i < nPix;
- //cout << " *** " << endl; fPixArray->Print();
-}
-
-//_____________________________________________________________________________
-void AliMUONClusterFinderMLEM::AdjustPixel(Double_t wxmin, Double_t wymin)
-{
-/// Check if some pixels have large size (adjust if necessary)
-
- AliDebug(2,Form("wxmin=%e wymin=%e",wxmin,wymin));
-
- Int_t n2[2], iOK = 1, nPix = fPixArray->GetEntriesFast();
- AliMUONPad *pixPtr, pix;
- Double_t xy0[2] = {9999, 9999}, wxy[2], dist[2] = {0};
-
- // Check if large pixel size
- for (Int_t i = 0; i < nPix; i++) {
- pixPtr = (AliMUONPad*) fPixArray->UncheckedAt(i);
- if (pixPtr->Charge() < 1) continue; // discarded pixel
- 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;
- Int_t update[2] = {0};
- //cout << " --- " << endl; fPixArray->Print();
- cout << xy0[0] << " " << xy0[1] << endl;
- for (Int_t i = 0; i < nPix; i++) {
- pixPtr = (AliMUONPad*) fPixArray->UncheckedAt(i);
- if (pixPtr->Charge() < 1) continue; // discarded pixel
- n2[0] = n2[1] = 1;
- update[0] = update[1] = 0;
- 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]; // distance to "normal" pixel
- // Go back to position of the first updated pixel
- dist[j] += (pixPtr->Size(j) - wxy[j]) * TMath::Sign(1.,-dist[j]);
- n2[j] = TMath::Nint (pixPtr->Size(j) / wxy[j]);
- update[j] = 1;
- }
- if (update[0] == 0 && update[1] == 0) continue;
- if (fDebug) cout << " Different " << pixPtr->Size(0) << " " << wxy[0] << " "
- << pixPtr->Size(1) << " " << wxy[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 (update[0]) pix.SetCoord(0, xy0[0] + dist[0] + TMath::Sign(2.,dist[0]) * ii * wxy[0]);
- for (Int_t jj = 0; jj < n2[1]; jj++) {
- if (update[1]) pix.SetCoord(1, xy0[1] + dist[1] + TMath::Sign(2.,dist[1]) * jj * wxy[1]);
- fPixArray->Add(new AliMUONPad(pix));
- //pix.Print();
- }
- }
- pixPtr->SetCharge(0);
- } // for (Int_t i = 0; i < nPix;
- cout << " *** " << endl; fPixArray->Print();
-}
-
//_____________________________________________________________________________
void
AliMUONClusterFinderMLEM::Plot(const char* basename)
{
/// Compute coefficients needed for MLEM algorithm
+ Int_t npadTot = cluster.Multiplicity();
Int_t nPix = fPixArray->GetLast()+1;
//memset(probi,0,nPix*sizeof(Double_t));
+ for (Int_t j = 0; j < npadTot*nPix; ++j) coef[j] = 0.;
for (Int_t j = 0; j < nPix; ++j) probi[j] = 0.;
Int_t mult = cluster.Multiplicity();
for ( Int_t ipix = 0; ipix < nPix; ++ipix )
{
Int_t indx1 = indx + ipix;
- if (pad->Status() < 0)
+ //if (pad->Status() < 0)
+ if (pad->Status() != fgkZero)
{
coef[indx1] = 0;
continue;
{
/// Repeat MLEM algorithm until pixel size becomes sufficiently small
- // AliCodeTimerAuto("")
+ // AliCodeTimerAuto("",)
Int_t nPix = fPixArray->GetLast()+1;
Int_t npadOK = 0;
for (Int_t i = 0; i < npadTot; ++i)
{
- if (cluster.Pad(i)->Status() == 0) ++npadOK;
+ //if (cluster.Pad(i)->Status() == 0) ++npadOK;
+ if (cluster.Pad(i)->Status() == fgkZero) ++npadOK;
}
- TH2D* mlem(0x0);
Double_t* coef(0x0);
Double_t* probi(0x0);
- Int_t lc(0); // loop counter (for debug)
+ Int_t lc(0); // loop counter
//Plot("mlem.start");
-
+ AliMUONPad* pixPtr = Pixel(0);
+ Double_t xylim[4] = {pixPtr->X(), -pixPtr->X(), pixPtr->Y(), -pixPtr->Y()};
+
while (1)
{
++lc;
- mlem = (TH2D*) gROOT->FindObject("mlem");
- delete mlem;
+ delete fHistMlem;
AliDebug(2,Form("lc %d nPix %d(%d) npadTot %d npadOK %d",lc,nPix,fPixArray->GetLast()+1,npadTot,npadOK));
AliDebug(2,Form("EVT%d PixArray=",fEventNumber));
// MLEM algorithm
Mlem(cluster,coef, probi, 15);
- Double_t xylim[4] = {999, 999, 999, 999};
- AliMUONPad* pixPtr(0x0);
-
- for ( Int_t ipix = 0; ipix < nPix; ++ipix )
- {
- pixPtr = Pixel(ipix);
- for ( Int_t i = 0; i < 4; ++i )
- {
- xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
- }
- }
+ // Find histogram limits for the 1'st pass only - for others computed below
+ if (lc == 1) {
+ for ( Int_t ipix = 1; ipix < nPix; ++ipix )
+ {
+ pixPtr = Pixel(ipix);
+ for ( Int_t i = 0; i < 2; ++i )
+ {
+ Int_t indx = i * 2;
+ if (pixPtr->Coord(i) < xylim[indx]) xylim[indx] = pixPtr->Coord(i);
+ else if (-pixPtr->Coord(i) < xylim[indx+1]) xylim[indx+1] = -pixPtr->Coord(i);
+ }
+ }
+ } else pixPtr = Pixel(0);
+
for (Int_t i = 0; i < 4; i++)
{
xylim[i] -= pixPtr->Size(i/2);
}
-
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);
xylim[0],-xylim[1],xylim[2],-xylim[3]
));
- mlem = new TH2D("mlem","mlem",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
+ fHistMlem = new TH2D("mlem","mlem",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
for (Int_t ipix = 0; ipix < nPix; ++ipix)
{
- AliMUONPad* pixPtr = Pixel(ipix);
- mlem->Fill(pixPtr->Coord(0),pixPtr->Coord(1),pixPtr->Charge());
+ AliMUONPad* pixPtr2 = Pixel(ipix);
+ fHistMlem->Fill(pixPtr2->Coord(0),pixPtr2->Coord(1),pixPtr2->Charge());
}
// Check if the total charge of pixels is too low
for ( Int_t i = 0; i < npadTot; ++i)
{
AliMUONPad* pad = cluster.Pad(i);
- if ( pad->Status() == 0) pad->SetStatus(-1);
+ //if ( pad->Status() == 0) pad->SetStatus(-1);
+ if ( pad->Status() == fgkZero) pad->SetStatus(fgkOver);
}
return kFALSE;
}
// Calculate position of the center-of-gravity around the maximum pixel
Double_t xyCOG[2];
- FindCOG(mlem, xyCOG);
+ FindCOG(xyCOG);
if (TMath::Min(pixPtr->Size(0),pixPtr->Size(1)) < 0.07 &&
pixPtr->Size(0) > pixPtr->Size(1)) break;
// Sort pixels according to the charge
+ MaskPeaks(1); // mask local maxima
fPixArray->Sort();
+ MaskPeaks(0); // unmask local maxima
Double_t pixMin = 0.01*Pixel(0)->Charge();
pixMin = TMath::Min(pixMin,50.);
nPix = 0;
for (Int_t ipix = 0; ipix < nPix1; ++ipix)
{
- AliMUONPad* pixPtr = Pixel(ipix);
+ AliMUONPad* pixPtr2 = Pixel(ipix);
if ( nPix >= npadOK // too many pixels already
||
- pixPtr->Charge() < pixMin // too low charge
+ ((pixPtr2->Charge() < pixMin) && (pixPtr2->Status() != fgkMustKeep)) // too low charge
)
{
RemovePixel(ipix);
{
if (!i)
{
- pixPtr->SetCharge(10);
- pixPtr->SetSize(indx, pixPtr->Size(indx)/2);
- width = -pixPtr->Size(indx);
- pixPtr->Shift(indx, width);
+ pixPtr2->SetCharge(10);
+ pixPtr2->SetSize(indx, pixPtr2->Size(indx)/2);
+ width = -pixPtr2->Size(indx);
+ pixPtr2->Shift(indx, width);
// Shift pixel position
if (ix)
{
ix = 0;
for (Int_t j = 0; j < 2; ++j)
{
- shift[j] = pixPtr->Coord(j) - xyCOG[j];
- shift[j] -= ((Int_t)(shift[j]/pixPtr->Size(j)/2))*pixPtr->Size(j)*2;
+ shift[j] = pixPtr2->Coord(j) - xyCOG[j];
+ shift[j] -= ((Int_t)(shift[j]/pixPtr2->Size(j)/2))*pixPtr2->Size(j)*2;
}
} // if (ix)
- pixPtr->Shift(0, -shift[0]);
- pixPtr->Shift(1, -shift[1]);
+ pixPtr2->Shift(0, -shift[0]);
+ pixPtr2->Shift(1, -shift[1]);
+ ++nPix;
}
- else
+ else if (nPix < npadOK)
{
- pixPtr = new AliMUONPad(*pixPtr);
- pixPtr->Shift(indx, -2*width);
- fPixArray->Add(pixPtr);
+ pixPtr2 = new AliMUONPad(*pixPtr2);
+ pixPtr2->Shift(indx, -2*width);
+ pixPtr2->SetStatus(fgkZero);
+ fPixArray->Add(pixPtr2);
+ ++nPix;
}
- for (Int_t i = 0; i < 4; ++i)
+ else continue; // skip adjustment of histo limits
+ for (Int_t j = 0; j < 4; ++j)
{
- xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
+ xylim[j] = TMath::Min (xylim[j], (j%2 ? -1 : 1)*pixPtr2->Coord(j/2));
}
} // for (Int_t i=0; i<2;
- nPix += 2;
} // for (Int_t ipix=0;
fPixArray->Compress();
- nPix = fPixArray->GetEntriesFast();
AliDebug(2,Form("After shift:"));
//StdoutToAliDebug(2,fPixArray->Print("","full"););
xylim[0],xylim[1],
xylim[2],xylim[3]));
- // Remove excessive pixels
- if (nPix > npadOK)
+ if (nPix < npadOK)
{
- for (Int_t ipix = npadOK; ipix < nPix; ++ipix)
- {
- RemovePixel(ipix);
- }
- }
- else
- {
- AliMUONPad* pixPtr = Pixel(0);
- // add pixels if the maximum is at the limit of pixel area
+ AliMUONPad* pixPtr2 = Pixel(0);
+ // add pixels if the maximum is at the limit of pixel area:
// start from Y-direction
Int_t j = 0;
for (Int_t i = 3; i > -1; --i)
{
if (nPix < npadOK &&
- TMath::Abs((i%2 ? -1 : 1)*xylim[i]-xyCOG[i/2]) < pixPtr->Size(i/2))
+ TMath::Abs((i%2 ? -1 : 1)*xylim[i]-xyCOG[i/2]) < pixPtr2->Size(i/2))
{
- AliMUONPad* p = static_cast<AliMUONPad*>(pixPtr->Clone());
- p->SetCoord(i/2, xyCOG[i/2]+(i%2 ? 2:-2)*pixPtr->Size(i/2));
+ //AliMUONPad* p = static_cast<AliMUONPad*>(pixPtr->Clone());
+ AliMUONPad* p = new AliMUONPad(*pixPtr2);
+ p->SetCoord(i/2, xyCOG[i/2]+(i%2 ? 2:-2)*pixPtr2->Size(i/2));
+ xylim[i] = p->Coord(i/2) * (i%2 ? -1 : 1); // update histo limits
j = TMath::Even (i/2);
p->SetCoord(j, xyCOG[j]);
AliDebug(2,Form("Adding pixel on the edge (i=%d) ",i));
}
}
}
- fPixArray->Compress();
nPix = fPixArray->GetEntriesFast();
delete [] coef;
delete [] probi;
// remove pixels with low signal or low visibility
// Cuts are empirical !!!
- Double_t thresh = TMath::Max (mlem->GetMaximum()/100.,1.);
+ Double_t thresh = TMath::Max (fHistMlem->GetMaximum()/100.,1.);
thresh = TMath::Min (thresh,50.);
Double_t charge = 0;
// Mark pixels which should be removed
for (Int_t i = 0; i < nPix; ++i)
{
- AliMUONPad* pixPtr = Pixel(i);
- charge = pixPtr->Charge();
+ AliMUONPad* pixPtr2 = Pixel(i);
+ charge = pixPtr2->Charge();
if (charge < thresh)
{
- pixPtr->SetCharge(-charge);
+ pixPtr2->SetCharge(-charge);
}
}
Int_t near = 0;
for (Int_t i = 0; i < nPix; ++i)
{
- AliMUONPad* pixPtr = Pixel(i);
- charge = pixPtr->Charge();
+ AliMUONPad* pixPtr2 = Pixel(i);
+ charge = pixPtr2->Charge();
if (charge > 0) continue;
- near = FindNearest(pixPtr);
- pixPtr->SetCharge(0);
+ near = FindNearest(pixPtr2);
+ pixPtr2->SetCharge(0);
probi[i] = 0; // make it "invisible"
AliMUONPad* pnear = Pixel(near);
pnear->SetCharge(pnear->Charge() + (-charge));
//StdoutToAliDebug(2,fPixArray->Print("","full"););
//Plot("mlem.beforesplit");
- // Update histogram
+ // Update histogram
for (Int_t i = 0; i < nPix; ++i)
{
- AliMUONPad* pixPtr = Pixel(i);
- Int_t ix = mlem->GetXaxis()->FindBin(pixPtr->Coord(0));
- Int_t iy = mlem->GetYaxis()->FindBin(pixPtr->Coord(1));
- mlem->SetBinContent(ix, iy, pixPtr->Charge());
+ AliMUONPad* pixPtr2 = Pixel(i);
+ Int_t ix = fHistMlem->GetXaxis()->FindBin(pixPtr2->Coord(0));
+ Int_t iy = fHistMlem->GetYaxis()->FindBin(pixPtr2->Coord(1));
+ fHistMlem->SetBinContent(ix, iy, pixPtr2->Charge());
}
// Try to split into clusters
Bool_t ok = kTRUE;
- if (mlem->GetSum() < 1)
+ if (fHistMlem->GetSum() < 1)
{
ok = kFALSE;
}
else
{
- fSplitter->Split(cluster,mlem,coef,fClusterList);
+ fSplitter->Split(cluster,fHistMlem,coef,fClusterList);
}
delete [] coef;
return ok;
}
+//_____________________________________________________________________________
+void AliMUONClusterFinderMLEM::MaskPeaks(Int_t mask)
+{
+ /// Mask/unmask pixels corresponding to local maxima (add/subtract 10000 to their charge
+ /// - to avoid loosing low charge pixels after sorting)
+
+ for (Int_t i = 0; i < fPixArray->GetEntriesFast(); ++i) {
+ AliMUONPad* pix = Pixel(i);
+ if (pix->Status() == fgkMustKeep) {
+ if (mask == 1) pix->SetCharge(pix->Charge()+10000.);
+ else pix->SetCharge(pix->Charge()-10000.);
+ }
+ }
+}
+
//_____________________________________________________________________________
void AliMUONClusterFinderMLEM::Mlem(AliMUONCluster& cluster,
- Double_t* coef, Double_t* probi,
+ const Double_t* coef, Double_t* probi,
Int_t nIter)
{
/// Use MLEM to find pixel charges
{
Pixel(ipix)->SetChargeBackup(0);
// Correct each pixel
+ probi1[ipix] = 0;
if (probi[ipix] < 0.01) continue; // skip "invisible" pixel
Double_t sum = 0;
probi1[ipix] = probMax;
for (Int_t j = 0; j < npad; ++j)
{
AliMUONPad* pad = cluster.Pad(j);
- if (pad->Status() < 0) continue;
+ //if (pad->Status() < 0) continue;
+ if (pad->Status() != fgkZero) continue;
Double_t sum1 = 0;
Int_t indx1 = j*nPix;
Int_t indx = indx1 + ipix;
}
if (sum1 > 1.e-6) sum += pad->Charge()*coef[indx]/sum1;
- //if (coef[indx] > 1.e-6) sum += pad->Charge()*coef[indx]/sum1;
} // for (Int_t j=0;
AliMUONPad* pixPtr = Pixel(ipix);
if (probi1[ipix] > 1.e-6)
}
//_____________________________________________________________________________
-void AliMUONClusterFinderMLEM::FindCOG(TH2D *mlem, Double_t *xyc)
+void AliMUONClusterFinderMLEM::FindCOG(Double_t *xyc)
{
/// Calculate position of the center-of-gravity around the maximum pixel
Int_t ixmax, iymax, ix, nsumx=0, nsumy=0, nsum=0;
Int_t i1 = -9, j1 = -9;
- mlem->GetMaximumBin(ixmax,iymax,ix);
- Int_t nx = mlem->GetNbinsX();
- Int_t ny = mlem->GetNbinsY();
- Double_t thresh = mlem->GetMaximum()/10;
+ fHistMlem->GetMaximumBin(ixmax,iymax,ix);
+ Int_t nx = fHistMlem->GetNbinsX();
+ Int_t ny = fHistMlem->GetNbinsY();
+ Double_t thresh = fHistMlem->GetMaximum()/10;
Double_t x, y, cont, xq=0, yq=0, qq=0;
Int_t ie = TMath::Min(ny,iymax+1), je = TMath::Min(nx,ixmax+1);
for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
- y = mlem->GetYaxis()->GetBinCenter(i);
+ y = fHistMlem->GetYaxis()->GetBinCenter(i);
for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
- cont = mlem->GetCellContent(j,i);
+ cont = fHistMlem->GetCellContent(j,i);
if (cont < thresh) continue;
if (i != i1) {i1 = i; nsumy++;}
if (j != j1) {j1 = j; nsumx++;}
- x = mlem->GetXaxis()->GetBinCenter(j);
+ x = fHistMlem->GetXaxis()->GetBinCenter(j);
xq += x*cont;
yq += y*cont;
qq += cont;
for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
if (i == iymax) continue;
for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
- cont = mlem->GetCellContent(j,i);
+ cont = fHistMlem->GetCellContent(j,i);
if (cont > cmax) {
cmax = cont;
- x = mlem->GetXaxis()->GetBinCenter(j);
- y = mlem->GetYaxis()->GetBinCenter(i);
+ x = fHistMlem->GetXaxis()->GetBinCenter(j);
+ y = fHistMlem->GetYaxis()->GetBinCenter(i);
i2 = i;
j2 = j;
}
for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
if (j == ixmax) continue;
for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
- cont = mlem->GetCellContent(j,i);
+ cont = fHistMlem->GetCellContent(j,i);
if (cont > cmax) {
cmax = cont;
- x = mlem->GetXaxis()->GetBinCenter(j);
- y = mlem->GetYaxis()->GetBinCenter(i);
+ x = fHistMlem->GetXaxis()->GetBinCenter(j);
+ y = fHistMlem->GetYaxis()->GetBinCenter(i);
i2 = i;
j2 = j;
}
}
//_____________________________________________________________________________
-Int_t AliMUONClusterFinderMLEM::FindNearest(AliMUONPad *pixPtr0)
+Int_t AliMUONClusterFinderMLEM::FindNearest(const AliMUONPad *pixPtr0)
{
/// Find the pixel nearest to the given one
/// (algorithm may be not very efficient)
for (Int_t i = 0; i < nPix; ++i) {
pixPtr = (AliMUONPad*) fPixArray->UncheckedAt(i);
- if (pixPtr == pixPtr0) continue;
- if (pixPtr->Charge() < 0.5) continue;
+ if (pixPtr == pixPtr0 || pixPtr->Charge() < 0.5) continue;
dx = (xc - pixPtr->Coord(0)) / pixPtr->Size(0);
dy = (yc - pixPtr->Coord(1)) / pixPtr->Size(1);
r = dx *dx + dy * dy;
AliDebug(1,Form("nPix=%d",pixArray->GetLast()+1));
- TH2D *hist = NULL;
- //if (pixArray == fPixArray) hist = (TH2D*) gROOT->FindObject("anode");
- //else { hist = (TH2D*) gROOT->FindObject("anode1"); cout << hist << endl; }
- //if (hist) hist->Delete();
-
Double_t xylim[4] = {999, 999, 999, 999};
Int_t nPix = pixArray->GetEntriesFast();
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);
- if (pixArray == fPixArray) hist = new TH2D("anode","anode",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
- else hist = new TH2D("anode1","anode1",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
+ if (pixArray == fPixArray) fHistAnode = new TH2D("anode","anode",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
+ else fHistAnode = new TH2D("anode1","anode1",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
for (Int_t ipix = 0; ipix < nPix; ++ipix) {
pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix);
- hist->Fill(pixPtr->Coord(0), pixPtr->Coord(1), pixPtr->Charge());
+ fHistAnode->Fill(pixPtr->Coord(0), pixPtr->Coord(1), pixPtr->Charge());
}
// if (fDraw && pixArray == fPixArray) fDraw->DrawHist("c2", hist);
for (Int_t i = 1; i <= ny; ++i) {
indx = (i-1) * nx;
for (Int_t j = 1; j <= nx; ++j) {
- if (hist->GetCellContent(j,i) < 0.5) continue;
+ if (fHistAnode->GetCellContent(j,i) < 0.5) continue;
//if (isLocalMax[indx+j-1] < 0) continue;
if (isLocalMax[indx+j-1] != 0) continue;
- FlagLocalMax(hist, i, j, isLocalMax);
+ FlagLocalMax(fHistAnode, i, j, isLocalMax);
}
}
for (Int_t j = 1; j <= nx; ++j) {
if (isLocalMax[indx+j-1] > 0) {
localMax[nMax] = indx + j - 1;
- maxVal[nMax++] = hist->GetCellContent(j,i);
- if (nMax > 99) AliFatal(" Too many local maxima !!!");
+ maxVal[nMax++] = fHistAnode->GetCellContent(j,i);
+ ((AliMUONPad*)fSplitter->BinToPix(fHistAnode, j, i))->SetStatus(fgkMustKeep);
+ if (nMax > 99) break;
}
}
+ if (nMax > 99) {
+ AliError(" Too many local maxima !!!");
+ break;
+ }
}
if (fDebug) cout << " Local max: " << nMax << endl;
delete [] isLocalMax;
- if (nMax == 1) hist->Delete();
return nMax;
}
//_____________________________________________________________________________
void AliMUONClusterFinderMLEM::FindCluster(AliMUONCluster& cluster,
- Int_t *localMax, Int_t iMax)
+ const Int_t *localMax, Int_t iMax)
{
/// Find pixel cluster around local maximum \a iMax and pick up pads
/// overlapping with it
- TH2D *hist = (TH2D*) gROOT->FindObject("anode");
/* Just for check
TCanvas* c = new TCanvas("Anode","Anode",800,600);
c->cd();
Int_t tmp;
cin >> tmp;
*/
- Int_t nx = hist->GetNbinsX();
- Int_t ny = hist->GetNbinsY();
+ Int_t nx = fHistAnode->GetNbinsX();
+ Int_t ny = fHistAnode->GetNbinsY();
Int_t ic = localMax[iMax] / nx + 1;
Int_t jc = localMax[iMax] % nx + 1;
Int_t nxy = ny * nx;
// Drop all pixels from the array - pick up only the ones from the cluster
fPixArray->Delete();
- Double_t wx = hist->GetXaxis()->GetBinWidth(1)/2;
- Double_t wy = hist->GetYaxis()->GetBinWidth(1)/2;
- Double_t yc = hist->GetYaxis()->GetBinCenter(ic);
- Double_t xc = hist->GetXaxis()->GetBinCenter(jc);
- Double_t cont = hist->GetCellContent(jc,ic);
+ Double_t wx = fHistAnode->GetXaxis()->GetBinWidth(1)/2;
+ Double_t wy = fHistAnode->GetYaxis()->GetBinWidth(1)/2;
+ Double_t yc = fHistAnode->GetYaxis()->GetBinCenter(ic);
+ Double_t xc = fHistAnode->GetXaxis()->GetBinCenter(jc);
+ Double_t cont = fHistAnode->GetCellContent(jc,ic);
fPixArray->Add(new AliMUONPad (xc, yc, wx, wy, cont));
used[(ic-1)*nx+jc-1] = kTRUE;
- AddBinSimple(hist, ic, jc);
+ AddBinSimple(fHistAnode, ic, jc);
//fSplitter->AddBin(hist, ic, jc, 1, used, (TObjArray*)0); // recursive call
Int_t nPix = fPixArray->GetEntriesFast();
// Pick up pads which overlap with found pixels
for (Int_t i = 0; i < npad; ++i)
{
- cluster.Pad(i)->SetStatus(-1);
+ //cluster.Pad(i)->SetStatus(-1);
+ cluster.Pad(i)->SetStatus(fgkOver); // just the dirty trick
}
for (Int_t i = 0; i < nPix; ++i)
for (Int_t j = 0; j < npad; ++j)
{
AliMUONPad* pad = cluster.Pad(j);
- if (pad->Status() == 0) continue;
+ //if (pad->Status() == 0) continue;
+ if (pad->Status() == fgkZero) continue;
if ( Overlap(*pad,*pixPtr) )
{
- pad->SetStatus(0);
+ //pad->SetStatus(0);
+ pad->SetStatus(fgkZero);
if (fDebug) { cout << j << " "; pad->Print("full"); }
}
}
//_____________________________________________________________________________
void
-AliMUONClusterFinderMLEM::AddBinSimple(TH2D *mlem, Int_t ic, Int_t jc)
+AliMUONClusterFinderMLEM::AddBinSimple(TH2D *hist, Int_t ic, Int_t jc)
{
/// Add adjacent bins (+-1 in X and Y) to the cluster
- Int_t nx = mlem->GetNbinsX();
- Int_t ny = mlem->GetNbinsY();
- Double_t cont1, cont = mlem->GetCellContent(jc,ic);
+ Int_t nx = hist->GetNbinsX();
+ Int_t ny = hist->GetNbinsY();
+ Double_t cont1, cont = hist->GetCellContent(jc,ic);
AliMUONPad *pixPtr = 0;
Int_t ie = TMath::Min(ic+1,ny), je = TMath::Min(jc+1,nx);
for (Int_t i = TMath::Max(ic-1,1); i <= ie; ++i) {
for (Int_t j = TMath::Max(jc-1,1); j <= je; ++j) {
- cont1 = mlem->GetCellContent(j,i);
+ cont1 = hist->GetCellContent(j,i);
if (cont1 > cont) continue;
if (cont1 < 0.5) continue;
- pixPtr = new AliMUONPad (mlem->GetXaxis()->GetBinCenter(j),
- mlem->GetYaxis()->GetBinCenter(i), 0, 0, cont1);
+ pixPtr = new AliMUONPad (hist->GetXaxis()->GetBinCenter(j),
+ hist->GetYaxis()->GetBinCenter(i), 0, 0, cont1);
fPixArray->Add(pixPtr);
}
}
return *this;
}
-//_____________________________________________________________________________
-void
-AliMUONClusterFinderMLEM::Neighbours(Int_t cathode, Int_t ix, Int_t iy,
- Int_t& n, Int_t* xList, Int_t* yList)
-{
- /// Get the list of neighbours of pad at (cathode,ix,iy)
- n = 0;
-
- const AliMpVSegmentation* seg = fSegmentation[cathode];
-
- AliMpPad pad = seg->PadByIndices(AliMpIntPair(ix,iy),kTRUE);
-
- // Define the region to look into : a region slightly bigger
- // than the pad itself (5% bigger), in order to catch first neighbours.
-
- AliMpArea area(pad.Position(),pad.Dimensions()*1.05);
-
- AliMpVPadIterator* it = seg->CreateIterator(area);
- it->First();
- while ( !it->IsDone() && n < 10 )
- {
- AliMpPad p = it->CurrentItem();
- if ( p != pad ) // skip self
- {
- xList[n] = p.GetIndices().GetFirst();
- yList[n] = p.GetIndices().GetSecond();
- ++n;
- }
- it->Next();
- }
- delete it;
-}
-
//_____________________________________________________________________________
void AliMUONClusterFinderMLEM::AddVirtualPad(AliMUONCluster& cluster)
{
Bool_t same = kFALSE;
if (TMath::Abs(dim0.Y()-dim1.Y()) < fgkDistancePrecision) same = kTRUE; // the same pad size on both planes
- AliMpIntPair cn;
+ Long_t cn;
Bool_t check[2] = {kFALSE, kFALSE};
Int_t nxy[2];
nxy[0] = nxy[1] = 0;
for (Int_t inb = 0; inb < 2; ++inb) {
cn = cluster.NofPads(nonb[inb], 0, kTRUE);
- if (inb == 0 && cn.GetFirst() == 2) check[inb] = kTRUE; // check non-bending plane
- else if (inb == 1 && cn.GetSecond() == 2) check[inb] = kTRUE; // check bending plane
+ if (inb == 0 && AliMp::PairFirst(cn) == 2) check[inb] = kTRUE; // check non-bending plane
+ else if (inb == 1 && AliMp::PairSecond(cn) == 2) check[inb] = kTRUE; // check bending plane
if (same) {
- nxy[0] = TMath::Max (nxy[0], cn.GetFirst());
- nxy[1] = TMath::Max (nxy[1], cn.GetSecond());
+ nxy[0] = TMath::Max (nxy[0], AliMp::PairFirst(cn));
+ nxy[1] = TMath::Max (nxy[1], AliMp::PairSecond(cn));
if (inb == 0 && nxy[0] < 2) nonb[inb] = !nonb[inb];
- else if (inb == 1 && cn.GetSecond() < 2) nonb[inb] = !nonb[inb];
+ else if (inb == 1 && AliMp::PairSecond(cn) < 2) nonb[inb] = !nonb[inb];
}
}
if (same) {
TVector2 pos;
if (inb == 0) pos.Set (pad->X() + idir * (pad->DX()+fgkDistancePrecision), pad->Y());
else pos.Set (pad->X(), pad->Y() + idir * (pad->DY()+fgkDistancePrecision));
- //AliMpPad mppad = fSegmentation[nonb[inb]]->PadByPosition(pos,kTRUE);
- AliMpPad mppad = fSegmentation[nonb[inb]]->PadByPosition(pos,kFALSE);
+ //AliMpPad mppad = fkSegmentation[nonb[inb]]->PadByPosition(pos,kTRUE);
+ AliMpPad mppad = fkSegmentation[nonb[inb]]->PadByPosition(pos.X(), pos.Y(),kFALSE);
if (!mppad.IsValid()) continue; // non-existing pad
- cn = mppad.GetIndices();
- AliMUONPad muonPad(fDetElemId, nonb[inb], cn.GetFirst(), cn.GetSecond(),
- mppad.Position().X(), mppad.Position().Y(),
- mppad.Dimensions().X(), mppad.Dimensions().Y(), 0);
+ AliMUONPad muonPad(fDetElemId, nonb[inb], mppad.GetIx(), mppad.GetIy(),
+ mppad.GetPositionX(), mppad.GetPositionY(),
+ mppad.GetDimensionX(), mppad.GetDimensionY(), 0);
if (inb == 0) muonPad.SetCharge(TMath::Min (amax[j]/100, 5.));
- else muonPad.SetCharge(TMath::Min (amax[j]/15, fgkZeroSuppression));
+ //else muonPad.SetCharge(TMath::Min (amax[j]/15, fgkZeroSuppression));
+ else muonPad.SetCharge(TMath::Min (amax[j]/15, 6.));
if (muonPad.Charge() < 1.) muonPad.SetCharge(1.);
muonPad.SetReal(kFALSE);
if (fDebug) printf(" ***** Add virtual pad in %d direction ***** %f %f %f %3d %3d %f %f \n",
/// Find number of pads in X and Y-directions (excluding virtual ones and
/// overflows)
- Int_t statusToTest = 1;
+ //Int_t statusToTest = 1;
+ Int_t statusToTest = fgkUseForFit;
- if ( nInX < 0 ) statusToTest = 0;
+ //if ( nInX < 0 ) statusToTest = 0;
+ if ( nInX < 0 ) statusToTest = fgkZero;
Bool_t mustMatch(kTRUE);
- AliMpIntPair cn = cluster.NofPads(statusToTest,mustMatch);
+ Long_t cn = cluster.NofPads(statusToTest,mustMatch);
- nInX = cn.GetFirst();
- nInY = cn.GetSecond();
+ nInX = AliMp::PairFirst(cn);
+ nInY = AliMp::PairSecond(cn);
}
//_____________________________________________________________________________
for (Int_t i = 0; i < mult; ++i)
{
AliMUONPad* pad = cluster.Pad(i);
+ /*
if ( pad->IsSaturated())
{
pad->SetStatus(-9);
{
pad->SetStatus(1);
}
+ */
+ if (!pad->IsSaturated()) pad->SetStatus(fgkUseForFit);
}
- nfit = fSplitter->Fit(cluster,1, nForFit, clustFit, clusters, parOk, fClusterList);
+ nfit = fSplitter->Fit(cluster,1, nForFit, clustFit, clusters, parOk, fClusterList, fHistMlem);
}
//_____________________________________________________________________________
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff