X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONClusterFinderMLEM.cxx;h=55b4799033a70fb80509733a05f408fea81b6242;hb=716ccba730076136b23f8c72e8be85ee1f3766e5;hp=4904413ec1af4469d1c2890c69b8662bb0b21f14;hpb=2abdae6e20a296abde062cc695230f47624f0451;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONClusterFinderMLEM.cxx b/MUON/AliMUONClusterFinderMLEM.cxx index 4904413ec1a..55b4799033a 100644 --- a/MUON/AliMUONClusterFinderMLEM.cxx +++ b/MUON/AliMUONClusterFinderMLEM.cxx @@ -45,21 +45,25 @@ #include #include #include -#include #include "AliCodeTimer.h" +using std::endl; +using std::cout; /// \cond CLASSIMP ClassImp(AliMUONClusterFinderMLEM) /// \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) @@ -70,23 +74,21 @@ fClusterList(), 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), fSplitter(0x0), fNClusters(0), -fNAddVirtualPads(0) +fNAddVirtualPads(0), +fLowestPixelCharge(0), +fLowestPadCharge(0), +fLowestClusterCharge(0) { /// 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; } @@ -95,7 +97,7 @@ fNAddVirtualPads(0) AliMUONClusterFinderMLEM::~AliMUONClusterFinderMLEM() { /// Destructor - delete fgMinuit; fgMinuit = 0; delete fPixArray; fPixArray = 0; + delete fPixArray; fPixArray = 0; // delete fDraw; delete fPreClusterFinder; delete fSplitter; @@ -105,46 +107,47 @@ AliMUONClusterFinderMLEM::~AliMUONClusterFinderMLEM() //_____________________________________________________________________________ Bool_t -AliMUONClusterFinderMLEM::Prepare(const AliMpVSegmentation* segmentations[2], - const AliMUONVDigitStore& digitStore) +AliMUONClusterFinderMLEM::Prepare(Int_t detElemId, + TObjArray* pads[2], + const AliMpArea& area, + const AliMpVSegmentation* seg[2]) { /// Prepare for clustering -// AliCodeTimerAuto("") +// AliCodeTimerAuto("",0) 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(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 = new AliMUONClusterSplitterMLEM(fDetElemId, + fPixArray, + fLowestPixelCharge, + fLowestPadCharge, + fLowestClusterCharge); 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(); - + fPixArray->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); + } } //_____________________________________________________________________________ @@ -152,10 +155,17 @@ AliMUONCluster* AliMUONClusterFinderMLEM::NextCluster() { /// Return next cluster -// AliCodeTimerAuto("") +// AliCodeTimerAuto("",0) // if the list of clusters is not void, pick one from there - TObject* o = fClusterList.At(++fClusterNumber); + TObject* o(0x0); + + // do we have clusters in our list ? + if ( fClusterNumber < fClusterList.GetLast() ) + { + o = fClusterList.At(++fClusterNumber); + } + if ( o != 0x0 ) return static_cast(o); //FIXME : at this point, must check whether we've used all the digits @@ -166,16 +176,17 @@ AliMUONClusterFinderMLEM::NextCluster() // pre-cluster and treat it fPreCluster = fPreClusterFinder->NextCluster(); - + + fPixArray->Delete(); + 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 @@ -201,10 +212,10 @@ AliMUONClusterFinderMLEM::WorkOnPreCluster() /// Starting from a precluster, builds a pixel array, and then /// extract clusters from this array - // AliCodeTimerAuto("") + // AliCodeTimerAuto("",0) 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); @@ -230,24 +241,21 @@ AliMUONClusterFinderMLEM::WorkOnPreCluster() 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) @@ -265,15 +273,17 @@ AliMUONClusterFinderMLEM::WorkOnPreCluster() 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 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; } @@ -299,16 +309,16 @@ AliMUONClusterFinderMLEM::CheckPrecluster(const AliMUONCluster& origCluster) /// Check precluster in order to attempt to simplify it (mostly for /// two-cathode preclusters) - // AliCodeTimerAuto("") + AliCodeTimerAuto("",0) // Disregard small clusters (leftovers from splitting or noise) if ((origCluster.Multiplicity()==1 || origCluster.Multiplicity()==2) && - origCluster.Charge(0)+origCluster.Charge(1) < 10) + origCluster.Charge(0)+origCluster.Charge(1) < fLowestClusterCharge ) { return 0x0; } - AliMUONCluster* cluster = static_cast(origCluster.Clone()); + AliMUONCluster* cluster = new AliMUONCluster(origCluster); AliDebug(2,"Start of CheckPreCluster="); //StdoutToAliDebug(2,cluster->Print("full")); @@ -332,9 +342,6 @@ AliMUONClusterFinderMLEM::CheckPreclusterTwoCathodes(AliMUONCluster* cluster) { /// 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; @@ -343,11 +350,11 @@ AliMUONClusterFinderMLEM::CheckPreclusterTwoCathodes(AliMUONCluster* cluster) 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 } @@ -372,9 +379,11 @@ AliMUONClusterFinderMLEM::CheckPreclusterTwoCathodes(AliMUONCluster* cluster) 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() < fLowestPadCharge) 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); @@ -392,8 +401,8 @@ AliMUONClusterFinderMLEM::CheckPreclusterTwoCathodes(AliMUONCluster* cluster) { // 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); @@ -414,8 +423,12 @@ AliMUONClusterFinderMLEM::CheckPreclusterTwoCathodes(AliMUONCluster* cluster) { 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; @@ -540,7 +553,7 @@ AliMUONClusterFinderMLEM::CheckOverlaps() 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; ); */ @@ -561,88 +574,23 @@ void AliMUONClusterFinderMLEM::BuildPixArray(AliMUONCluster& cluster) } 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; @@ -657,54 +605,50 @@ void AliMUONClusterFinderMLEM::BuildPixArrayOneCathode(AliMUONCluster& cluster) // 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); @@ -720,12 +664,11 @@ void AliMUONClusterFinderMLEM::BuildPixArrayOneCathode(AliMUONCluster& cluster) 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 @@ -733,29 +676,44 @@ void AliMUONClusterFinderMLEM::BuildPixArrayOneCathode(AliMUONCluster& cluster) 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 (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 (gROOT->FindObject("Grid")); - TH2D *hist2 = static_cast (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 @@ -764,14 +722,15 @@ void AliMUONClusterFinderMLEM::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, Ali 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); } } @@ -780,232 +739,35 @@ void AliMUONClusterFinderMLEM::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, Ali 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; jPosition(),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] <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) +AliMUONClusterFinderMLEM::Plot(const char* /*basename*/) { /// Make a plot and save it as png return; //AZ - if (!fPlot) return; - - TCanvas* c = new TCanvas("MLEM","MLEM",800,600); - c->Draw(); - Draw(); - c->Modified(); - c->Update(); - c->Print(Form("%s.EVT%d.DE%d.CLU%d.png",basename,fEventNumber, - fDetElemId,fClusterNumber)); +// if (!fPlot) return; +// +// TCanvas* c = new TCanvas("MLEM","MLEM",800,600); +// c->Draw(); +// Draw(); +// c->Modified(); +// c->Update(); +// c->Print(Form("%s.EVT%d.DE%d.CLU%d.png",basename,fEventNumber, +// fDetElemId,fClusterNumber)); } //_____________________________________________________________________________ @@ -1016,9 +778,11 @@ AliMUONClusterFinderMLEM::ComputeCoefficients(AliMUONCluster& cluster, { /// 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(); @@ -1030,7 +794,8 @@ AliMUONClusterFinderMLEM::ComputeCoefficients(AliMUONCluster& cluster, 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; @@ -1057,7 +822,7 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple { /// Repeat MLEM algorithm until pixel size becomes sufficiently small - // AliCodeTimerAuto("") + // AliCodeTimerAuto("",0) Int_t nPix = fPixArray->GetLast()+1; @@ -1075,25 +840,27 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple 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; + fHistMlem = 0x0; 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)); - //StdoutToAliDebug(2,fPixArray->Print("","full")); + //StdoutToAliDebug(2,fPixArray->Print("full")); coef = new Double_t [npadTot*nPix]; probi = new Double_t [nPix]; @@ -1115,22 +882,24 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple // 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); @@ -1141,12 +910,16 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple xylim[0],-xylim[1],xylim[2],-xylim[3] )); - mlem = new TH2D("mlem","mlem",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]); + AliDebug(2,Form("LowestPadCharge=%e",fLowestPadCharge)); + + delete fHistMlem; + + 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 @@ -1156,7 +929,7 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple qTot += Pixel(i)->Charge(); } - if ( qTot < 1.e-4 || ( npadOK < 3 && qTot < 7 ) ) + if ( qTot < 1.e-4 || ( npadOK < 3 && qTot < fLowestClusterCharge ) ) { AliDebug(1,Form("Deleting the above cluster (charge %e too low, npadOK=%d)",qTot,npadOK)); delete [] coef; @@ -1165,7 +938,8 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple 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; } @@ -1182,15 +956,17 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple // 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.); + pixMin = TMath::Min(pixMin,100*fLowestPixelCharge); // Decrease pixel size and shift pixels to make them centered at // the maximum one @@ -1203,10 +979,10 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple 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); @@ -1216,39 +992,41 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple { if (!i) { - pixPtr->SetCharge(10); - pixPtr->SetSize(indx, pixPtr->Size(indx)/2); - width = -pixPtr->Size(indx); - pixPtr->Shift(indx, width); + pixPtr2->SetCharge(fLowestPadCharge); + 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");); @@ -1259,27 +1037,21 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple 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(pixPtr->Clone()); - p->SetCoord(i/2, xyCOG[i/2]+(i%2 ? 2:-2)*pixPtr->Size(i/2)); + //AliMUONPad* p = static_cast(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)); @@ -1290,7 +1062,6 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple } } } - fPixArray->Compress(); nPix = fPixArray->GetEntriesFast(); delete [] coef; delete [] probi; @@ -1301,18 +1072,18 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple // remove pixels with low signal or low visibility // Cuts are empirical !!! - Double_t thresh = TMath::Max (mlem->GetMaximum()/100.,1.); - thresh = TMath::Min (thresh,50.); + Double_t thresh = TMath::Max (fHistMlem->GetMaximum()/100.,2.0*fLowestPixelCharge); + thresh = TMath::Min (thresh,100.0*fLowestPixelCharge); 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); } } @@ -1320,11 +1091,11 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple 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)); @@ -1335,24 +1106,24 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple //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() < 2.0*fLowestPixelCharge) { ok = kFALSE; } else { - fSplitter->Split(cluster,mlem,coef,fClusterList); + fSplitter->Split(cluster,fHistMlem,coef,fClusterList); } delete [] coef; @@ -1362,9 +1133,24 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple 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 @@ -1383,13 +1169,15 @@ void AliMUONClusterFinderMLEM::Mlem(AliMUONCluster& cluster, { 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; @@ -1408,7 +1196,6 @@ void AliMUONClusterFinderMLEM::Mlem(AliMUONCluster& cluster, } 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) @@ -1434,27 +1221,27 @@ void AliMUONClusterFinderMLEM::Mlem(AliMUONCluster& cluster, } //_____________________________________________________________________________ -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; @@ -1470,11 +1257,11 @@ void AliMUONClusterFinderMLEM::FindCOG(TH2D *mlem, Double_t *xyc) 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; } @@ -1494,11 +1281,11 @@ void AliMUONClusterFinderMLEM::FindCOG(TH2D *mlem, Double_t *xyc) 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; } @@ -1517,7 +1304,7 @@ void AliMUONClusterFinderMLEM::FindCOG(TH2D *mlem, Double_t *xyc) } //_____________________________________________________________________________ -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) @@ -1529,8 +1316,7 @@ Int_t AliMUONClusterFinderMLEM::FindNearest(AliMUONPad *pixPtr0) 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() < fLowestPixelCharge) continue; dx = (xc - pixPtr->Coord(0)) / pixPtr->Size(0); dy = (yc - pixPtr->Coord(1)) / pixPtr->Size(1); r = dx *dx + dy * dy; @@ -1607,14 +1393,12 @@ Int_t AliMUONClusterFinderMLEM::FindLocalMaxima(TObjArray *pixArray, Int_t *loca 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(); + + if ( nPix <= 0 ) return 0; + AliMUONPad *pixPtr = 0; for (Int_t ipix = 0; ipix < nPix; ++ipix) { pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix); @@ -1625,11 +1409,11 @@ Int_t AliMUONClusterFinderMLEM::FindLocalMaxima(TObjArray *pixArray, Int_t *loca 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); @@ -1640,10 +1424,10 @@ Int_t AliMUONClusterFinderMLEM::FindLocalMaxima(TObjArray *pixArray, Int_t *loca 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) < fLowestPixelCharge) 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); } } @@ -1652,14 +1436,18 @@ Int_t AliMUONClusterFinderMLEM::FindLocalMaxima(TObjArray *pixArray, Int_t *loca 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; } @@ -1699,12 +1487,11 @@ void AliMUONClusterFinderMLEM::FlagLocalMax(TH2D *hist, Int_t i, Int_t j, Int_t //_____________________________________________________________________________ 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(); @@ -1713,8 +1500,8 @@ void AliMUONClusterFinderMLEM::FindCluster(AliMUONCluster& cluster, 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; @@ -1724,14 +1511,14 @@ void AliMUONClusterFinderMLEM::FindCluster(AliMUONCluster& cluster, // 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(); @@ -1747,7 +1534,8 @@ void AliMUONClusterFinderMLEM::FindCluster(AliMUONCluster& cluster, // 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) @@ -1756,10 +1544,12 @@ void AliMUONClusterFinderMLEM::FindCluster(AliMUONCluster& cluster, 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"); } } } @@ -1770,23 +1560,23 @@ void AliMUONClusterFinderMLEM::FindCluster(AliMUONCluster& cluster, //_____________________________________________________________________________ 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); + if (cont1 < fLowestPixelCharge) continue; + pixPtr = new AliMUONPad (hist->GetXaxis()->GetBinCenter(j), + hist->GetYaxis()->GetBinCenter(i), 0, 0, cont1); fPixArray->Add(pixPtr); } } @@ -1805,39 +1595,6 @@ AliMUONClusterFinderMLEM::operator=(const AliMUONClusterFinderMLEM& rhs) 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) { @@ -1857,19 +1614,19 @@ 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) { @@ -1931,16 +1688,16 @@ void AliMUONClusterFinderMLEM::AddVirtualPad(AliMUONCluster& cluster) 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); - if (inb == 0) muonPad.SetCharge(TMath::Min (amax[j]/100, 5.)); - else muonPad.SetCharge(TMath::Min (amax[j]/15, fgkZeroSuppression)); - if (muonPad.Charge() < 1.) muonPad.SetCharge(1.); + 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, fLowestPadCharge)); + //else muonPad.SetCharge(TMath::Min (amax[j]/15, fgkZeroSuppression)); + else muonPad.SetCharge(TMath::Min (amax[j]/15, fLowestPadCharge)); + if (muonPad.Charge() < 2.0*fLowestPixelCharge) muonPad.SetCharge(2.0*fLowestPixelCharge); muonPad.SetReal(kFALSE); if (fDebug) printf(" ***** Add virtual pad in %d direction ***** %f %f %f %3d %3d %f %f \n", inb, muonPad.Charge(), muonPad.X(), muonPad.Y(), muonPad.Ix(), @@ -1959,16 +1716,18 @@ void AliMUONClusterFinderMLEM::PadsInXandY(AliMUONCluster& cluster, /// 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); } //_____________________________________________________________________________ @@ -1996,6 +1755,7 @@ void AliMUONClusterFinderMLEM::Simple(AliMUONCluster& cluster) for (Int_t i = 0; i < mult; ++i) { AliMUONPad* pad = cluster.Pad(i); + /* if ( pad->IsSaturated()) { pad->SetStatus(-9); @@ -2004,8 +1764,10 @@ void AliMUONClusterFinderMLEM::Simple(AliMUONCluster& cluster) { 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); } //_____________________________________________________________________________ @@ -2029,4 +1791,14 @@ AliMUONClusterFinderMLEM::Print(Option_t* what) const } } +//_____________________________________________________________________________ +void +AliMUONClusterFinderMLEM::SetChargeHints(Double_t lowestPadCharge, Double_t lowestClusterCharge) +{ + /// Set some thresholds we use later on in the algorithm + fLowestPadCharge=lowestPadCharge; + fLowestClusterCharge=lowestClusterCharge; + fLowestPixelCharge=fLowestPadCharge/12.0; +} +