X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=MUON%2FAliMUONClusterFinderMLEM.cxx;h=6596734d2fdf9fa140729e5c963498029965a391;hb=4011b280271f6006055ee9b7568f745ad583c8e3;hp=888edf48e71bd21b08031d5c224a40e971384202;hpb=72dae9fff6e7e1656e0f9069eb01b49bbd1f3040;p=u%2Fmrichter%2FAliRoot.git diff --git a/MUON/AliMUONClusterFinderMLEM.cxx b/MUON/AliMUONClusterFinderMLEM.cxx index 888edf48e71..6596734d2fd 100644 --- a/MUON/AliMUONClusterFinderMLEM.cxx +++ b/MUON/AliMUONClusterFinderMLEM.cxx @@ -51,8 +51,6 @@ ClassImp(AliMUONClusterFinderMLEM) /// \endcond -const Double_t AliMUONClusterFinderMLEM::fgkZeroSuppression = 6; // average zero suppression value -//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); @@ -81,10 +79,13 @@ fDebug(0), fPlot(plot), fSplitter(0x0), fNClusters(0), -fNAddVirtualPads(0) +fNAddVirtualPads(0), +fLowestPixelCharge(0), +fLowestPadCharge(0), +fLowestClusterCharge(0) { /// Constructor - + fkSegmentation[1] = fkSegmentation[0] = 0x0; if (fPlot) fDebug = 1; @@ -110,7 +111,7 @@ AliMUONClusterFinderMLEM::Prepare(Int_t detElemId, const AliMpVSegmentation* seg[2]) { /// Prepare for clustering -// AliCodeTimerAuto("") +// AliCodeTimerAuto("",0) for ( Int_t i = 0; i < 2; ++i ) { @@ -121,7 +122,11 @@ AliMUONClusterFinderMLEM::Prepare(Int_t detElemId, 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 @@ -129,7 +134,8 @@ AliMUONClusterFinderMLEM::Prepare(Int_t detElemId, fEventNumber = runLoader ? runLoader->GetEventNumber() : 0; fClusterNumber = -1; fClusterList.Delete(); - + fPixArray->Delete(); + AliDebug(3,Form("EVT %d DE %d",fEventNumber,fDetElemId)); if ( fPreClusterFinder->NeedSegmentation() ) @@ -147,7 +153,7 @@ 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); @@ -162,6 +168,7 @@ AliMUONClusterFinderMLEM::NextCluster() fPreCluster = fPreClusterFinder->NextCluster(); + fPixArray->Delete(); fClusterList.Delete(); // reset the list of clusters for this pre-cluster fClusterNumber = -1; //AZ @@ -196,7 +203,7 @@ AliMUONClusterFinderMLEM::WorkOnPreCluster() /// Starting from a precluster, builds a pixel array, and then /// extract clusters from this array - // AliCodeTimerAuto("") + // AliCodeTimerAuto("",0) if (fDebug) { cout << " *** Event # " << fEventNumber @@ -293,11 +300,11 @@ 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; } @@ -363,10 +370,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 = fkSegmentation[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() < 20) 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); @@ -536,7 +544,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; ); */ @@ -608,19 +616,22 @@ void AliMUONClusterFinderMLEM::BuildPixArrayOneCathode(AliMUONCluster& cluster) if (cluster.Multiplicity(0) == 0) cath0 = 1; else if (cluster.Multiplicity(1) == 0) cath1 = 0; - TVector2 leftDown = cluster.Area(cath0).LeftDownCorner(); - TVector2 rightUp = cluster.Area(cath0).RightUpCorner(); - min[0] = leftDown.X(); - min[1] = leftDown.Y(); - max[0] = rightUp.X(); - max[1] = rightUp.Y(); + + 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) { - leftDown = cluster.Area(cath1).LeftDownCorner(); - rightUp = cluster.Area(cath1).RightUpCorner(); - min[0] = TMath::Max (min[0], leftDown.X()); - min[1] = TMath::Max (min[1], leftDown.Y()); - max[0] = TMath::Min (max[0], rightUp.X()); - max[1] = TMath::Min (max[1], rightUp.Y()); + 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 @@ -802,7 +813,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; @@ -836,10 +847,11 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple { ++lc; 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]; @@ -889,6 +901,10 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple xylim[0],-xylim[1],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) @@ -904,7 +920,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; @@ -941,7 +957,7 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple 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 @@ -967,7 +983,7 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple { if (!i) { - pixPtr2->SetCharge(10); + pixPtr2->SetCharge(fLowestPadCharge); pixPtr2->SetSize(indx, pixPtr2->Size(indx)/2); width = -pixPtr2->Size(indx); pixPtr2->Shift(indx, width); @@ -1047,8 +1063,8 @@ 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 (fHistMlem->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 @@ -1092,7 +1108,7 @@ Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple // Try to split into clusters Bool_t ok = kTRUE; - if (fHistMlem->GetSum() < 1) + if (fHistMlem->GetSum() < 2.0*fLowestPixelCharge) { ok = kFALSE; } @@ -1291,7 +1307,7 @@ Int_t AliMUONClusterFinderMLEM::FindNearest(const AliMUONPad *pixPtr0) for (Int_t i = 0; i < nPix; ++i) { pixPtr = (AliMUONPad*) fPixArray->UncheckedAt(i); - if (pixPtr == pixPtr0 || 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; @@ -1396,7 +1412,7 @@ 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 (fHistAnode->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(fHistAnode, i, j, isLocalMax); @@ -1546,7 +1562,7 @@ AliMUONClusterFinderMLEM::AddBinSimple(TH2D *hist, Int_t ic, Int_t jc) for (Int_t j = TMath::Max(jc-1,1); j <= je; ++j) { cont1 = hist->GetCellContent(j,i); if (cont1 > cont) continue; - if (cont1 < 0.5) continue; + if (cont1 < fLowestPixelCharge) continue; pixPtr = new AliMUONPad (hist->GetXaxis()->GetBinCenter(j), hist->GetYaxis()->GetBinCenter(i), 0, 0, cont1); fPixArray->Add(pixPtr); @@ -1586,19 +1602,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) { @@ -1661,16 +1677,15 @@ void AliMUONClusterFinderMLEM::AddVirtualPad(AliMUONCluster& cluster) 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 = fkSegmentation[nonb[inb]]->PadByPosition(pos,kTRUE); - AliMpPad mppad = fkSegmentation[nonb[inb]]->PadByPosition(pos,kFALSE); + 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.)); + 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, 6.)); - if (muonPad.Charge() < 1.) muonPad.SetCharge(1.); + 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(), @@ -1697,10 +1712,10 @@ void AliMUONClusterFinderMLEM::PadsInXandY(AliMUONCluster& cluster, 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); } //_____________________________________________________________________________ @@ -1764,4 +1779,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; +} +