AlgSubShowerSRList Class Reference

#include <AlgSubShowerSRList.h>

Inheritance diagram for AlgSubShowerSRList:

List of all members.

Public Member Functions
	AlgSubShowerSRList ()
virtual	~AlgSubShowerSRList ()
virtual void	RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void	Trace (const char *c) const
Static Public Member Functions
NavKey	StripKeyFromPlane (const CandStripHandle *)
Private Member Functions
Bool_t	FindCluster (TObjArray *, TObjArray *, Int_t)
Bool_t	TestOverLap (TObjArray *, CandHandle &, const CandSliceHandle *)
std::vector< window >	TransCluster (std::vector< window >, Int_t)
std::vector< window >	HoughTransCluster (std::vector< window >, Int_t)
Bool_t	MergeCluster (TObjArray *, CandHandle &, Int_t, Bool_t test=false)
Bool_t	CleanUp (TObjArray *, CandHandle &, Int_t, const CandSliceHandle *cslice=NULL)
Bool_t	FormHalo (TObjArray *, CandHandle &, CandContext &, AlgHandle &, const CandSliceListHandle *, Int_t)
Bool_t	BestHough (std::vector< window >, Double_t *, Bool_t)
Private Attributes
Double_t	fMergeFracEnergyLimit
Double_t	fMergeBestFracBoundary
Double_t	fMergeBestFracWithin
Double_t	fMergeBestFracCross
Double_t	fLongWindowMipCut
Int_t	fMaxPlaneGap

---

Constructor & Destructor Documentation

AlgSubShowerSRList::AlgSubShowerSRList (   )

Definition at line 57 of file AlgSubShowerSRList.cxx. : fMergeFracEnergyLimit(0.1),fMergeBestFracBoundary(0.8), fMergeBestFracWithin(0.6),fMergeBestFracCross(0.9), fLongWindowMipCut(0.0),fMaxPlaneGap(2) { }

AlgSubShowerSRList::~AlgSubShowerSRList (   )  [virtual]

Definition at line 65 of file AlgSubShowerSRList.cxx. { }

---

Member Function Documentation

Bool_t AlgSubShowerSRList::BestHough (  std::vector< window >  , Double_t *  , Bool_t  )  [private]

Definition at line 3109 of file AlgSubShowerSRList.cxx. References MSG, name, CandHandle::Print(), and STRIPWIDTHINMETRES. Referenced by HoughTransCluster(), and TransCluster(). { std::vector<window>::iterator winIterBeg = winVec.begin(); std::vector<window>::iterator winIterEnd = winVec.end(); //get extrema of windows: Int_t MaxPlane = 0; Int_t MinPlane = 500; Double_t MaxStrip = -5.; Double_t MinStrip = 5.; //total PH: Double_t totPH = 0; //PH of window with max PH: Double_t maxWinPH = 0; //plane with maximum PH, and value of max PH: Int_t maxPlanePos = -1; Double_t maxPlanePosPH = 0; //plane by plane sum PH: Double_t PlanePH[500] = {0}; while(winIterBeg!=winIterEnd){ if(MaxPlane<winIterBeg->plane) MaxPlane = winIterBeg->plane; if(MinPlane>winIterBeg->plane) MinPlane = winIterBeg->plane; if(MaxStrip<winIterBeg->upper_tpos) MaxStrip = winIterBeg->upper_tpos; if(MinStrip>winIterBeg->lower_tpos) MinStrip = winIterBeg->lower_tpos; if(maxWinPH<winIterBeg->ph) maxWinPH = winIterBeg->ph; totPH+=winIterBeg->ph; PlanePH[winIterBeg->plane] += winIterBeg->ph; if(PlanePH[winIterBeg->plane]>maxPlanePosPH) { maxPlanePos = winIterBeg->plane; maxPlanePosPH = PlanePH[maxPlanePos]; } winIterBeg++; } //window tpos position with max PH, and value of max PH Double_t maxWindowPos = 0.; Double_t maxWindowPH = 0.; winIterBeg = winVec.begin(); while(winIterBeg!=winIterEnd){ if(winIterBeg->plane==maxPlanePos){ if(winIterBeg->ph>maxWindowPH){ maxWindowPos = winIterBeg->phpos; maxWindowPH = winIterBeg->ph; } } winIterBeg++; } if(winVec.size()<=1 || MinPlane==MaxPlane){ MCV[0] = 0; MCV[1] = maxWindowPos; MCV[2] = maxPlanePos; MCV[3] = 0; MCV[4] = 0; MCV[5] = 0; MCV[6] = maxWindowPos; MCV[7] = maxPlanePos; MCV[8] = 0; MCV[9] = 0; MCV[10] = 0; MCV[11] = maxWindowPos; MCV[12] = 0; MCV[13] = maxWindowPos; return true; } Double_t TPosWin = MaxStrip-MinStrip+0.041; MSG("SubShowerSR",Msg::kDebug) << "Hough extrema: plane " << MinPlane << " to " << MaxPlane << " strip tpos " << MinStrip << " to " << MaxStrip << endl; //define "x=0" to be at MinPlane-2, //so y-intercept is strip value at MinPlane-2 //gradient is in units of Delta(tpos) per plane TH2F *houghHist = new TH2F("houghHist","Hough Histogram", Int_t(TPosWin/STRIPWIDTHINMETRES) + MaxPlane-MinPlane+2, -2.*TPosWin/Float_t(MaxPlane-MinPlane+2), 2.*TPosWin/Float_t(MaxPlane-MinPlane+2), Int_t(2.*TPosWin/STRIPWIDTHINMETRES), MinStrip-TPosWin, MaxStrip+TPosWin); TH2F *phHoughHist = new TH2F("phHoughHist","PH Weighted Hough Histogram", Int_t(TPosWin/STRIPWIDTHINMETRES) + MaxPlane-MinPlane+2, -2.*TPosWin/Float_t(MaxPlane-MinPlane+2), 2.*TPosWin/Float_t(MaxPlane-MinPlane+2), Int_t(2.*TPosWin/STRIPWIDTHINMETRES), MinStrip-TPosWin, MaxStrip+TPosWin); winIterBeg = winVec.begin(); while(winIterBeg!=winIterEnd){ Int_t plane = winIterBeg->plane; Double_t strip = winIterBeg->phpos; Double_t ph = winIterBeg->ph; Int_t lastBin = -1; for(Double_t c=MinStrip-TPosWin; c<=MaxStrip+TPosWin; c+=STRIPWIDTHINMETRES){ Double_t m = (strip - c)/Float_t(plane-MinPlane+2); MSG("SubShowerSR",Msg::kVerbose) << "Hough entries: " << m << " " << c << " " << ph << endl; if(houghHist->FindBin(m,c)!=lastBin){ houghHist->Fill(m,c,1); phHoughHist->Fill(m,c,TMath::Sqrt(ph)); lastBin = houghHist->FindBin(m,c); } } winIterBeg++; } Int_t maxBin = houghHist->GetMaximumBin(); Float_t maxVal = houghHist->GetMaximum(); Int_t xbin = maxBin%(houghHist->GetNbinsX()+2); Int_t ybin = maxBin/(houghHist->GetNbinsX()+2); Int_t phMaxBin = phHoughHist->GetMaximumBin(); Float_t phMaxVal = phHoughHist->GetMaximum(); Int_t phXbin = phMaxBin%(phHoughHist->GetNbinsX()+2); Int_t phYbin = phMaxBin/(phHoughHist->GetNbinsX()+2); Double_t totx = 0; Double_t sumx=0,sumx2=0; Double_t ptotx = 0; Double_t psumx=0,psumx2=0; Double_t toty=0; Double_t sumy=0,sumy2=0; Double_t ptoty = 0; Double_t psumy=0,psumy2=0; for(int xx = 1;xx<=houghHist->GetNbinsX();xx++){ for(int yy = 1;yy<=houghHist->GetNbinsY();yy++){ if(houghHist->GetBinContent(xx,yy)>maxVal*0.75) { totx += houghHist->GetBinContent(xx,yy); sumx += ( houghHist->GetBinContent(xx,yy) * houghHist->GetXaxis()->GetBinCenter(xx) ); sumx2 += ( houghHist->GetBinContent(xx,yy) * houghHist->GetXaxis()->GetBinCenter(xx) * houghHist->GetXaxis()->GetBinCenter(xx) ); toty += houghHist->GetBinContent(xx,yy); sumy += ( houghHist->GetBinContent(xx,yy) * houghHist->GetYaxis()->GetBinCenter(yy) ); sumy2 += ( houghHist->GetBinContent(xx,yy) * houghHist->GetYaxis()->GetBinCenter(yy) * houghHist->GetYaxis()->GetBinCenter(yy) ); } if(phHoughHist->GetBinContent(xx,yy)>phMaxVal*0.75) { ptotx += phHoughHist->GetBinContent(xx,yy); psumx += ( phHoughHist->GetBinContent(xx,yy) * phHoughHist->GetXaxis()->GetBinCenter(xx) ); psumx2 += ( phHoughHist->GetBinContent(xx,yy) * phHoughHist->GetXaxis()->GetBinCenter(xx) * phHoughHist->GetXaxis()->GetBinCenter(xx) ); ptoty += phHoughHist->GetBinContent(xx,yy); psumy += ( phHoughHist->GetBinContent(xx,yy) * phHoughHist->GetYaxis()->GetBinCenter(yy) ); psumy2 += ( phHoughHist->GetBinContent(xx,yy) * phHoughHist->GetYaxis()->GetBinCenter(yy) * phHoughHist->GetYaxis()->GetBinCenter(yy) ); } } } MCV[2] = MinPlane - 2; if(totx>0) { MCV[0] = sumx/totx; MCV[3] = (sumx2/totx - TMath::Power(MCV[0],2))/totx; if(MCV[3]>0) MCV[3] = TMath::Sqrt(MCV[3]); else MCV[3] = 0; } if(toty>0) { MCV[1] = sumy/toty; MCV[4] = (sumy2/toty - TMath::Power(MCV[1],2))/toty; if(MCV[4]>0) MCV[4] = TMath::Sqrt(MCV[4]); else MCV[4] = 0; } if(MCV[3]<houghHist->GetXaxis()->GetBinWidth(1)) MCV[3] = houghHist->GetXaxis()->GetBinWidth(1); if(MCV[4]<houghHist->GetYaxis()->GetBinWidth(1)) MCV[4] = houghHist->GetYaxis()->GetBinWidth(1); MCV[7] = MinPlane - 2; if(ptotx>0) { MCV[5] = psumx/ptotx; MCV[8] = (psumx2/ptotx - TMath::Power(MCV[5],2))/ptotx; if(MCV[8]>0) MCV[8] = TMath::Sqrt(MCV[8]); else MCV[8] = 0; } if(ptoty>0) { MCV[6] = psumy/ptoty; MCV[9] = (psumy2/ptoty - TMath::Power(MCV[6],2))/ptoty; if(MCV[9]>0) MCV[9] = TMath::Sqrt(MCV[9]); else MCV[9] = 0; } if(MCV[8]<phHoughHist->GetXaxis()->GetBinWidth(1)) MCV[8] = phHoughHist->GetXaxis()->GetBinWidth(1); if(MCV[9]<phHoughHist->GetYaxis()->GetBinWidth(1)) MCV[9] = phHoughHist->GetYaxis()->GetBinWidth(1); MCV[10] = houghHist->GetXaxis()->GetBinCenter(xbin); MCV[11] = houghHist->GetYaxis()->GetBinCenter(ybin); MCV[12] = phHoughHist->GetXaxis()->GetBinCenter(phXbin); MCV[13] = phHoughHist->GetYaxis()->GetBinCenter(phYbin); if(makeEPS) { gStyle->SetOptStat(0); TCanvas *can = new TCanvas(); can->Divide(2,2); can->cd(1); houghHist->Draw("COLZ"); can->cd(2); phHoughHist->Draw("COLZ"); can->cd(3); houghHist->Draw("LEGO2"); can->cd(4); phHoughHist->Draw("LEGO2"); Int_t countr = 0; Bool_t noPlot = true; while(noPlot){ char name[256]; sprintf(name,"best_%i.root",countr); std::ifstream Test(name); if(!Test) { //if files does not exist: can->Print(name); noPlot = false; } countr+=1; } delete can; } delete houghHist; delete phHoughHist; MSG("SubShowerSR",Msg::kDebug) << "Before Cuts in BestHough:" << "\nHough gradient, intercept, plane vertex: " << MCV[0] << "+/-" << MCV[3] << " " << MCV[1] << "+/-" << MCV[4] << " " << MCV[2] << "\nPH Hough gradient, intercept, plane vertex: " << MCV[5] << "+/-" << MCV[8] << " " << MCV[6] << "+/-" << MCV[9] << " " << MCV[7] << "\n Hough Max Bin Gradient and Intercept Values: " << MCV[10] << " " << MCV[11] << " " << MCV[12] << " " << MCV[13] << endl; //check that there is a coincidence of 3 or more windows //and gradient is bigger than error if(maxVal < 3 || (TMath::Abs(MCV[0])-MCV[3]<0 && TMath::Abs(MCV[10])-MCV[3]<0)) { MCV[0] = 0; MCV[1] = maxWindowPos; MCV[2] = maxPlanePos; MCV[3] = 0; MCV[4] = 0; MCV[10] = 0; MCV[11] = maxWindowPos; } if((phMaxVal - maxWindowPH)/(totPH - maxWindowPH) < 0.1 || (TMath::Abs(MCV[5])-MCV[8]<0 ) && (TMath::Abs(MCV[12])-MCV[8]<0 )) { MCV[5] = 0; MCV[6] = maxWindowPos; MCV[7] = maxPlanePos; MCV[8] = 0; MCV[9] = 0; MCV[12] = 0; MCV[13] = maxWindowPos; } return true; }

Bool_t AlgSubShowerSRList::CleanUp (  TObjArray *  , CandHandle &  , Int_t  , const CandSliceHandle *  cslice = NULL )  [private]

Definition at line 2568 of file AlgSubShowerSRList.cxx. References CandHandle::AddDaughterLink(), CandSubShowerSRHandle::GetAveTime(), CandSubShowerSRHandle::GetAvgDev(), CandRecoHandle::GetBegPlane(), CandRecoHandle::GetCandSlice(), CandRecoHandle::GetCharge(), CandStripHandle::GetCharge(), CandHandle::GetDaughterIterator(), CandRecoHandle::GetEndPlane(), CandSubShowerSRHandle::GetEnergy(), Calibrator::GetMIP(), CandHandle::GetNDaughters(), CandRecoHandle::GetNStrip(), CandStripHandle::GetPlane(), CandSubShowerSRHandle::GetPlaneView(), CandStripHandle::GetPlaneView(), CandSubShowerSRHandle::GetSlope(), CandStripHandle::GetStrip(), CandStripHandle::GetTime(), CandStripHandle::GetTPos(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), CandStripHandle::GetZPos(), Calibrator::Instance(), MSG, and CandSubShowerSRHandle::SetEnergy(). Referenced by FormHalo(), and RunAlg(). { MSG("SubShowerSR", Msg::kDebug) << "Running CleanUp()" << endl; Int_t totUnassigned = allStrips->GetLast()+1; if(totUnassigned==0) return false; Calibrator& calibrator=Calibrator::Instance(); assert(ch.InheritsFrom("CandSubShowerSRListHandle")); CandSubShowerSRListHandle &csslh = dynamic_cast<CandSubShowerSRListHandle &>(ch); //Get planeview for this cleanup TObject *firstObj = allStrips->At(0); assert(firstObj->InheritsFrom("CandStripHandle")); CandStripHandle *firstStp = dynamic_cast<CandStripHandle*>(firstObj); PlaneView::PlaneView_t pv = firstStp->GetPlaneView(); MSG("SubShowerSR", Msg::kDebug) << "PlaneView = " << pv << endl; //Any changes made? Bool_t anyChanges = false; //once strips have been assigned do not try to assign again: Bool_t *useStrip = new Bool_t[totUnassigned]; Bool_t *assignedStrip = new Bool_t[totUnassigned]; for(int i=0;i<totUnassigned;i++) { assignedStrip[i] = false; //if cslice defined, check that strip is part of slice: if(cslice){ useStrip[i] = false; TObject *tobj = allStrips->At(i); CandStripHandle *unassigned = dynamic_cast<CandStripHandle*>(tobj); CandStripHandleItr stpSliceItr(cslice->GetDaughterIterator()); while(CandStripHandle *stp = stpSliceItr()) { if(*stp==*unassigned) { useStrip[i] = true; break; } } } else useStrip[i] = true; } MSG("SubShowerSR", Msg::kDebug) << "Checking " << totUnassigned << " strips" << endl; //going to check if any isolated strips are nearby within +/- 2 strip for (Int_t i=0; i<totUnassigned; i++) { MSG("SubShowerSR", Msg::kVerbose) << "Checking unassigned strip " << i << endl; TObject *tobj = allStrips->At(i); if(!useStrip[i]) continue; MSG("SubShowerSR", Msg::kVerbose) << "Can still use unassigned strip " << i << endl; CandStripHandle *unassigned = dynamic_cast<CandStripHandle*>(tobj); Float_t *StripDistance = new Float_t[csslh.GetNDaughters()]; Float_t theCharge = calibrator.GetMIP(unassigned->GetCharge(CalDigitType::kSigCorr)); //loop through existing subshowers Int_t ssCount = 0; CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator()); while (CandSubShowerSRHandle *subshower = subshowerItr()) { StripDistance[ssCount] = 99999; //set to 999999 as a flag //increment early because of continue statements, later use ssCount-1 ssCount++; //if a slice is given in argument list //demand that slices are the same if(cslice){ if(*subshower->GetCandSlice()!=*cslice) continue; } else if(TMath::Abs(subshower->GetAveTime() - unassigned->GetTime())*1.0e9 > 100.0) continue; //only look at the right planeview subshowers if(subshower->GetPlaneView()!=pv) continue; MSG("SubShowerSR", Msg::kVerbose) << "SubShower " << ssCount-1 << " in the same planeview" << endl; // if plane of strip is more than 1 plane (of same type away) // don't try to add it to this subshower // (also beware of SM boundary in FarDet) Int_t thePlane = unassigned->GetPlane(); if(thePlane<subshower->GetBegPlane()-2 || thePlane>subshower->GetEndPlane()+2 || (det==1 && ( (subshower->GetBegPlane()>249 && thePlane<249) || (subshower->GetEndPlane()<249 && thePlane>249) ))) continue; MSG("SubShowerSR", Msg::kVerbose) << "SubShower " << ssCount-1 << " in the same plane range" << endl; //loop through strips in subshower //check for proximity to other strips in subshower //(on a plane-by-plane basis) Bool_t checkProximity = false; Float_t theStrip_tpos = unassigned->GetTPos(); CandStripHandleItr stripItr(subshower->GetDaughterIterator()); while(CandStripHandle *stp = stripItr()){ //check we're on the same plane //or on the begin plane if the unassigned plane < begin plane //or on the end plane if the unassigned plane > end plane if((stp->GetPlane() == thePlane) || (thePlane < subshower->GetBegPlane() && stp->GetPlane() == subshower->GetBegPlane()) || (thePlane > subshower->GetEndPlane() && stp->GetPlane() == subshower->GetEndPlane())) { //allow +/- 2.5 strips: if(TMath::Abs(theStrip_tpos-stp->GetTPos())<=2.5*STRIPWIDTHINMETRES) { MSG("SubShowerSR",Msg::kVerbose) << "TMath::Abs(theStrip_tpos-stp->GetTPos())=" << TMath::Abs(theStrip_tpos-stp->GetTPos()) << " stp->Plane()=" << stp->GetPlane() << " stp->Strip()=" << stp->GetStrip() << endl; checkProximity = true; } } } if(!checkProximity) continue; //not close enough to this subshower MSG("SubShowerSR", Msg::kVerbose) << "SubShower " << ssCount-1 << " has strips in the same tpos range" << endl; //calculate pH weighted distance from cluster centre //get subshower vertex and angle: Double_t tposVtx = 0; if(pv==PlaneView::kU || pv==PlaneView::kX) tposVtx = subshower->GetVtxU(); else if(pv==PlaneView::kV || pv==PlaneView::kY) tposVtx = subshower->GetVtxV(); Double_t zposVtx = subshower->GetVtxZ(); Float_t slope = subshower->GetSlope(); Float_t avgdev = subshower->GetAvgDev(); //get strip position and charge: Float_t tpos = unassigned->GetTPos(); Float_t zpos = unassigned->GetZPos(); if(avgdev!=0.) { StripDistance[ssCount-1] = (TMath::Abs(tpos - (tposVtx + slope*(zpos - zposVtx))) * TMath::Cos(TMath::ATan(slope)))*theCharge/avgdev; } else { MSG("SubShowerSR",Msg::kVerbose) << "SubShower " << ssCount-1 << " has AvgDev = 0!" << endl; StripDistance[ssCount-1] = 99998; } MSG("SubShowerSR", Msg::kVerbose) << "Unassigned strip " << i << " and SubShower " << ssCount-1 << " are a good match with dev(strip)/avgdev(subshower) = " << StripDistance[ssCount-1] << endl; } //make decision here about closest subshower Float_t minSD = 99999; Int_t bestSS = -1; for(int j=0;j<csslh.GetNDaughters();j++){ if(StripDistance[j]<minSD) { minSD = StripDistance[j]; bestSS = j; } } if(bestSS!=-1){ MSG("SubShowerSR", Msg::kVerbose) << "Best match for strip " << i << ": SubShower " << bestSS << " Adding plane: " << unassigned->GetPlane() << " strip: " << unassigned->GetStrip() << endl; CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator()); ssCount = 0; while(CandSubShowerSRHandle *cssh = subshowerItr()){ if(ssCount==bestSS){ MSG("SubShowerSR",Msg::kVerbose) << "Before: cssh->GetNStrip() = " << cssh->GetNStrip() << " cssh->GetEnergy() = " << cssh->GetEnergy() << endl; cssh->AddDaughterLink(*unassigned); Double_t mip2gev = cssh->GetEnergy() / calibrator.GetMIP(cssh->GetCharge(CalStripType::kSigCorr)); cssh->SetEnergy(cssh->GetEnergy()+mip2gev*theCharge); MSG("SubShowerSR",Msg::kVerbose) << "After: cssh->GetNStrip() = " << cssh->GetNStrip() << " cssh->GetEnergy() = " << cssh->GetEnergy() << endl; assignedStrip[i] = true; anyChanges = true; break; } ssCount+=1; } } else { MSG("SubShowerSR", Msg::kVerbose) << "No matching subshower for unassigned strip " << i << endl; } delete [] StripDistance; } //remove assigned strips Int_t assigned = 0; for(int i=0;i<totUnassigned;i++) { if(assignedStrip[i]) { allStrips->RemoveAt(i); assigned+=1; } } allStrips->Compress(); if(anyChanges) { MSG("SubShowerSR", Msg::kDebug) << assigned << "/" << totUnassigned << " strips assigned to SubShowers" << endl; } else { MSG("SubShowerSR", Msg::kDebug) << "No Changes made to SubShowers" << endl; } return anyChanges; }

Bool_t AlgSubShowerSRList::FindCluster (  TObjArray *  allStrips, TObjArray *  newSubShower, Int_t  det )  [private]

Read in info: find maximum PH plane within selected longitudinal window: Definition at line 510 of file AlgSubShowerSRList.cxx. References window::centroid, CandStripHandle::GetCharge(), Calibrator::GetMIP(), CandStripHandle::GetPlane(), CandStripHandle::GetStrip(), CandStripHandle::GetTPos(), CandHandle::GetVldContext(), CandStripHandle::GetZPos(), Calibrator::Instance(), window::lower, window::lower_tpos, MSG, window::ph, window::phpos, window::phwidth, window::plane, s(), TransCluster(), window::type, window::upper, window::upper_tpos, and win. Referenced by RunAlg(). { Int_t nstp = allStrips->GetLast()+1; if(nstp==0){ MSG("SubShower", Msg::kError) << "No strips present!" << endl; return true; } //Get useful info from candidates for clustering Int_t *plane = new Int_t[nstp]; Int_t *strip = new Int_t[nstp]; Double_t *ph = new Double_t[nstp]; Double_t *pherr = new Double_t[nstp]; Double_t *z = new Double_t[nstp]; Double_t *tpos = new Double_t[nstp]; //neighbouring strip info Double_t *TransGradientPlus = new Double_t[nstp]; Double_t *TransGradientMinus = new Double_t[nstp]; Double_t *TransGradientErrorPlus = new Double_t[nstp]; Double_t *TransGradientErrorMinus = new Double_t[nstp]; CandStripHandle *usefulStrip; int begPlane = 500; int endPlane = 0; Int_t doubleCounter = 0; Int_t alreadyGot[500][192] = {{0}}; Calibrator& calibrator=Calibrator::Instance(); for (Int_t i=0; i<=allStrips->GetLast(); i++) { TObject *tobj = allStrips->At(i); assert(tobj->InheritsFrom("CandStripHandle")); CandStripHandle *stp = dynamic_cast<CandStripHandle*>(tobj); if(alreadyGot[stp->GetPlane()][stp->GetStrip()]==1) { for(int j=0;j<i-doubleCounter;j++){ if(plane[j]==stp->GetPlane() && strip[j]==stp->GetStrip()){ ph[j] += calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr)); TransGradientPlus[j] = ph[j]; TransGradientMinus[j] = ph[j]; break; } } doubleCounter+=1; nstp -= 1; } else { alreadyGot[stp->GetPlane()][stp->GetStrip()] = 1; plane[i-doubleCounter] = stp->GetPlane(); if (plane[i-doubleCounter]<begPlane) begPlane = plane[i-doubleCounter]; if (plane[i-doubleCounter]>endPlane) endPlane = plane[i-doubleCounter]; strip[i-doubleCounter] = stp->GetStrip(); ph[i-doubleCounter] = calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr)); pherr[i-doubleCounter] = ph[i-doubleCounter] / TMath::Sqrt(stp->GetCharge(CalDigitType::kPE)); //pe stats err z[i-doubleCounter] = stp->GetZPos(); tpos[i-doubleCounter] = stp->GetTPos(); TransGradientPlus[i-doubleCounter] = ph[i-doubleCounter]; TransGradientErrorPlus[i-doubleCounter] = TMath::Power(pherr[i-doubleCounter],2); TransGradientMinus[i-doubleCounter] = ph[i-doubleCounter]; TransGradientErrorMinus[i-doubleCounter] = TMath::Power(pherr[i-doubleCounter],2); usefulStrip = stp; } } //find transverse neighbour strip index and PH gradient for each strip for(Int_t i=0;i<nstp;i++){ for(Int_t j=0;j<nstp;j++){ if (plane[j]==plane[i]){ if(strip[j]==strip[i]-1) { TransGradientMinus[i] = TransGradientMinus[i] - ph[j]; TransGradientErrorMinus[i] = TMath::Sqrt(TransGradientErrorMinus[i]+ TMath::Power(pherr[j],2)); } if(strip[j]==strip[i]+1) { TransGradientPlus[i] = TransGradientPlus[i] - ph[j]; TransGradientErrorPlus[i] = TMath::Sqrt(TransGradientErrorPlus[i]+ TMath::Power(pherr[j],2)); TransGradientErrorMinus[i] = 0; TransGradientErrorPlus[i] = 0; } } } } //find per plane PH int npls = (endPlane-begPlane)/2+1; if(det==1&&(begPlane-249)*(endPlane-249)<0) npls = (endPlane-begPlane-1)/2+1; Int_t *allpl = new Int_t[npls]; Double_t *allplph = new Double_t[npls]; Int_t *allstp = new Int_t[npls]; for(Int_t i=0;i<npls;i++){ allpl[i] = 2*i+begPlane; if(det==1&&(begPlane-249)*(endPlane-249)<0&&(2*i+begPlane)>249) allpl[i] = 2*i+begPlane+1; allplph[i] = 0.; allstp[i] = 0; for(Int_t j=0;j<nstp;j++){ if (plane[j]==allpl[i]){ allplph[i] += ph[j]; allstp[i] += 1; } } } UgliGeomHandle ugh(*usefulStrip->GetVldContext()); //Start Clustering: //Find contiguous longitudinal windows MSG("SubShowerSR",Msg::kDebug) <<"start longitudinal clustering"<<endl; std::vector<int> zeroPlanes; //vector of planes with PH = 0 int start = begPlane-2; if(det==1&&(begPlane-2-249)*(begPlane-249)<0) start = begPlane-3; zeroPlanes.push_back(start-2); zeroPlanes.push_back(start); for(Int_t i=0;i<npls;i++){ if (allplph[i]<=fLongWindowMipCut) zeroPlanes.push_back(allpl[i]); } int over = endPlane+2; if(det==1&&(endPlane+2-249)*(endPlane-249)<0) over = endPlane+3; zeroPlanes.push_back(over); zeroPlanes.push_back(over+2); MSG("SubShowerSR",Msg::kVerbose) << "Number of Zero Planes = " << zeroPlanes.size() << " First Zero Plane = " << start << " Last Zero Plane = "<< over << endl; //---------------------------------------------- //Find "Zero PH" windows //position of first plane in a longitudinal gap window: std::vector<int> firstGapPlane; //position of last plane in a longitudinal gap window: std::vector<int> lastGapPlane; int begGap = 0; int endGap = -5; //filling gap window vectors: for (UInt_t i = 0; i<zeroPlanes.size(); i++){ bool foundgap = false; if(zeroPlanes.at(i)>endGap) { begGap = zeroPlanes.at(i); endGap = zeroPlanes.at(i); int j = 0; while ((i+j+1)<zeroPlanes.size() && zeroPlanes.at(i+j+1)-zeroPlanes.at(i+j)<4) { foundgap = true; endGap = zeroPlanes.at(i+j+1); j++; } if (foundgap){ firstGapPlane.push_back(begGap); lastGapPlane.push_back(endGap); } } } //Find highest strip density longitudinal window Int_t maxWinBegPlane = 0; //first boundary plane of max window Int_t maxWinEndPlane = 0; Float_t maxDensity = 0; for (UInt_t i = 0; i<firstGapPlane.size(); i++){ if((i+1)<firstGapPlane.size()){ Float_t density = 0; Int_t bpln = lastGapPlane.at(i); Int_t epln = firstGapPlane.at(i+1); for(int j=0;j<npls;j++){ if(allpl[j]>=bpln && allpl[j]<=epln){ density += Float_t(allstp[j]); } } density/=Float_t(epln-bpln+2); MSG("SubShowerSR",Msg::kDebug) << "Proto-long window [" << bpln << "," << epln << "] density = " << density << endl; if(density > maxDensity) { maxWinBegPlane = bpln; maxWinEndPlane = epln; maxDensity = density; } } } MSG("SubShowerSR",Msg::kDebug) << "Using Longitudinal Window: " << maxWinBegPlane << "-" << maxWinEndPlane << endl; Int_t maxPlane = -1; Double_t maxPlanePH = 0; for(int i=0;i<npls;i++){ if(allpl[i]<maxWinEndPlane && allpl[i]>maxWinBegPlane) { if(allplph[i]>maxPlanePH) { maxPlane = allpl[i]; maxPlanePH = allplph[i]; } } } MSG("SubShowerSR",Msg::kVerbose) << "Maximum PH plane within Long. window = " << maxPlane << endl; Int_t numPreMaxPlanes = maxPlane - maxWinBegPlane; Int_t numPostMaxPlanes = maxWinEndPlane - maxPlane; Int_t numLongPlanes = numPostMaxPlanes + numPreMaxPlanes - 1; if(numLongPlanes/2<2){ for (Int_t i=0; i<=allStrips->GetLast(); i++) { TObject *tobj = allStrips->At(i); assert(tobj->InheritsFrom("CandStripHandle")); CandStripHandle *stp = dynamic_cast<CandStripHandle*>(tobj); if(stp->GetPlane()>=maxWinBegPlane && stp->GetPlane()<=maxWinEndPlane){ newSubShower->Add(stp); } } return true; } //Start Transverse Clustering //window info: std::vector<window> winvec; //vector containing all transverse windows //Loop over all planes in long. window //Find transverse windows per plane for (int i=0;i<numLongPlanes;i++){ //find next appropriate plane to look at starting from max plane: int pl = -1; if(numPreMaxPlanes>numPostMaxPlanes){ if (i<numPreMaxPlanes) pl = maxPlane-i; else pl = maxPlane + i - numPreMaxPlanes + 1; //add 1 to prevent re-doing max plane } else{ if (i<numPostMaxPlanes) pl = maxPlane+i; else pl = maxPlane - i + numPostMaxPlanes - 1; //subtract 1 to prevent re-doing max plane } std::vector<Int_t> winu; // vector of upper dip std::vector<Int_t> wind; // vector of lower dip std::vector<Double_t> winu_tpos; // vector of tpos of upper dip std::vector<Double_t> wind_tpos; // vector of tpos of lower dip std::vector<Int_t> wintu; // vector of upper dip type std::vector<Int_t> wintd; // vector of lower dip type // //find dips and categorize them MSG("SubShowerSR",Msg::kVerbose) << "Find Upper/Lower boundaries in plane = " << pl << endl; for(Int_t i=0;i<nstp;i++){ if(plane[i]==pl){ //using neighbouring strip info from above: //(only called dips if delta PH is significant given pe stats) if(TransGradientPlus[i]<0 /*-(TransGradientErrorPlus[i])*/ && TransGradientMinus[i]<0 /*-(TransGradientErrorMinus[i])*/){ winu.push_back(strip[i]); winu_tpos.push_back(tpos[i]); wintu.push_back(0); // type-0 (this is a dip) MSG("SubShowerSR",Msg::kVerbose) << "up dip strip = " << strip[i] << " ph = " << ph[i] << " tgp = " << TransGradientPlus[i] << " tgm = " << TransGradientMinus[i] << endl; } else if(TransGradientPlus[i]==ph[i] && TransGradientMinus[i]<=ph[i]){ winu.push_back(strip[i]); winu_tpos.push_back(tpos[i]); MSG("SubShowerSR",Msg::kVerbose) << "up edge/iso strip = " << strip[i] << " ph = " << ph[i] << " tgp = " << TransGradientPlus[i] << " tgm = " << TransGradientMinus[i] << endl; if (TransGradientMinus[i]<ph[i]) wintu.push_back(1); // type-1 (this is an edge strip) else wintu.push_back(2); // type-2 (this is an isolated strip) } if(TransGradientPlus[i]<0 /*-(TransGradientErrorPlus[i])*/ && TransGradientMinus[i]<0 /*-(TransGradientErrorMinus[i])*/ ){ wind.push_back(strip[i]); wind_tpos.push_back(tpos[i]); wintd.push_back(0); // type-0 MSG("SubShowerSR",Msg::kVerbose) << "down dip strip = " << strip[i] << " ph = " << ph[i] << " tgp = " << TransGradientPlus[i] << " tgm = " << TransGradientMinus[i] << endl; } else if(TransGradientMinus[i]==ph[i] && TransGradientPlus[i]<=ph[i]){ wind.push_back(strip[i]); wind_tpos.push_back(tpos[i]); MSG("SubShowerSR",Msg::kVerbose) << "down edge/iso strip = " << strip[i] << " ph = " << ph[i] << " tgp = " << TransGradientPlus[i] << " tgm = " << TransGradientMinus[i] << endl; if (TransGradientPlus[i]<ph[i]) wintd.push_back(1); //type-1 else wintd.push_back(2); //type-2 } } } //use dips to form windows UInt_t wins = winu.size(); if(wins!=wind.size()) { MSG("SubShowerSR",Msg::kError) << "Number of Upper Transverse Window Boundaries does not equal " << "number of Lower Transverse Window Boundaries... " << "Leaving FindCluster()" << endl; MSG("SubShowerSR",Msg::kError) << "winu.size() = " << winu.size() << endl; MSG("SubShowerSR",Msg::kError) << "wind.size() = " << wind.size() << endl; MSG("SubShowerSR",Msg::kError) << "wins = " << wins << endl; return false; } //temporarily hold window info Int_t *win0 = new Int_t[wins]; // strip of upper dip Int_t *win1 = new Int_t[wins]; // strip of lower dip Float_t *win0_tpos = new Float_t[wins]; // tpos of strip of upper dip Float_t *win1_tpos = new Float_t[wins]; // tpos of strip of lower dip Float_t *win2 = new Float_t[wins]; // PH of window Int_t *win3 = new Int_t[wins]; // ID of window: // 0-both dips; // 1 or 10-one dip,one gap; // 11-both gaps; // 22-isolated strip Float_t *win4 = new Float_t[wins]; // Centroid of window in tpos Float_t *win5 = new Float_t[wins]; // PH weighted tpos of window Float_t *win6 = new Float_t[wins]; // PH weighted width of window //match up dips to form windows for (UInt_t w = 0; w<wins; w++){ win0[w] = -1; win1[w] = -1; win0_tpos[w] = 0; win1_tpos[w] = 0; win2[w] = 0; win3[w] = -1; win4[w] = 0; win5[w] = 0; win6[w] = 0; int winw = 200; for (UInt_t s = 0; s<wind.size(); s++){ int diff = winu.at(w)-wind.at(s); if ( ( (diff>0 && wintu.at(w)<2 && wintd.at(s)<2) || (diff==0 && wintu.at(w)==2 && wintd.at(s)==2) ) && diff<winw) { winw = diff; win0[w] = winu.at(w); win1[w] = wind.at(s); win0_tpos[w] = winu_tpos.at(w); win1_tpos[w] = wind_tpos.at(s); win3[w] = wintu.at(w)*10+wintd.at(s); } } } //now get PH of windows //and push_back windows into vectors for (UInt_t w = 0; w<wins; w++){ Double_t biggestStripPH = 0; Int_t n = 0; for(Int_t i=0;i<nstp;i++){ if (plane[i]==pl && strip[i]<=win0[w] && strip[i]>=win1[w]) { win2[w] += ph[i]; if(ph[i]>biggestStripPH) { win4[w] = tpos[i]; biggestStripPH = ph[i]; } win5[w] += ph[i]*tpos[i]; win6[w] += ph[i]*tpos[i]*tpos[i]; n++; } } win5[w] /= win2[w]; if(n!=1){ win6[w] /= win2[w]; win6[w] -= win5[w]*win5[w]; if(win6[w]>0) win6[w] = TMath::Sqrt(win6[w]); else win6[w] = 0.041/TMath::Sqrt(12.); } else win6[w] = 0.041/TMath::Sqrt(12.); window win; win.plane = pl; win.upper = win0[w]; win.lower = win1[w]; win.upper_tpos = win0_tpos[w]; win.lower_tpos = win1_tpos[w]; win.ph = win2[w]; win.type = win3[w]; win.centroid = win4[w]; win.phpos = win5[w]; win.phwidth = win6[w]; winvec.push_back(win); } delete [] win0; delete [] win1; delete [] win0_tpos; delete [] win1_tpos; delete [] win2; delete [] win3; delete [] win4; delete [] win5; delete [] win6; } //Perform clustering on windows: std::vector<window> clusterWinVec = TransCluster(winvec,det); //no windows added to cluster: if(clusterWinVec.size()==0) return false; //Check for longitudinal gaps within the proto-cluster //Remove windows after first gap on both sides of shower max plane maxPlane = -1; //recalculate maxPlane maxPlanePH = 0.; Double_t PlanePH[500] = {0.}; Int_t non0Win = 0; std::vector<window>::iterator winIter = clusterWinVec.begin(); std::vector<window>::iterator winEnd = clusterWinVec.end(); while(winIter!=winEnd){ if(winIter->plane>0){ PlanePH[winIter->plane] += winIter->ph; if(PlanePH[winIter->plane]>0) non0Win++; if(PlanePH[winIter->plane]>maxPlanePH) { maxPlane = winIter->plane; maxPlanePH = PlanePH[maxPlane]; } } winIter++; } if(non0Win==0){ MSG("SubShowerSR",Msg::kDebug) <<"No window solutions found"<<endl; return false; } Int_t maxLongPlane = -1; Int_t minLongPlane = -1; Int_t counter = 0; Int_t delta = 2; while(PlanePH[maxPlane+counter]>0||(maxPlane+counter)==249) { maxLongPlane = maxPlane+counter; delta = 2; if(det==1&&(maxLongPlane-249)*(maxLongPlane+delta-249)<0) delta = 3; counter+=delta; } counter = 0; while(PlanePH[maxPlane+counter]>0||(maxPlane+counter)==249) { minLongPlane = maxPlane+counter; delta = 2; if(det==1&&(minLongPlane-249)*(minLongPlane-delta-249)<0) delta = 3; counter-=delta; } MSG("SubShowerSR",Msg::kDebug) <<"New Longitudinal Window: "<< minLongPlane << "-" << maxLongPlane << endl; for (Int_t i=0; i<=allStrips->GetLast(); i++) { TObject *tobj = allStrips->At(i); assert(tobj->InheritsFrom("CandStripHandle")); CandStripHandle *stp = dynamic_cast<CandStripHandle*>(tobj); if(stp->GetPlane()>=minLongPlane && stp->GetPlane()<=maxLongPlane){ winIter = clusterWinVec.begin(); while(winIter!=winEnd){ if(winIter->plane==stp->GetPlane()){ Int_t strp = stp->GetStrip(); if(strp<=winIter->upper && strp>=winIter->lower) { newSubShower->Add(stp); break; } } winIter++; } } } delete [] plane; delete [] strip; delete [] ph; delete [] z; delete [] tpos; delete [] TransGradientPlus; delete [] TransGradientMinus; delete [] TransGradientErrorPlus; delete [] TransGradientErrorMinus; delete [] allplph; delete [] allpl; delete [] allstp; return true; }

Bool_t AlgSubShowerSRList::FormHalo (  TObjArray *  , CandHandle &  , CandContext &  , AlgHandle &  , const CandSliceListHandle *  , Int_t  )  [private]

Definition at line 2808 of file AlgSubShowerSRList.cxx. References CandHandle::AddDaughterLink(), CleanUp(), CandSubShowerSRHandle::GetAvgDev(), CandRecoHandle::GetBegPlane(), CandRecoHandle::GetCandSlice(), CandSubShowerSRHandle::GetClusterID(), CandHandle::GetDaughterIterator(), CandRecoHandle::GetEndPlane(), CandSubShowerSRHandle::GetEnergy(), CandHandle::GetNDaughters(), CandRecoHandle::GetNStrip(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandSubShowerSRHandle::GetPlaneView(), CandSubShowerSRHandle::GetSlope(), CandStripHandle::GetTPos(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandSubShowerSR::MakeCandidate(), MSG, PITCHINMETRES, CandHandle::Print(), CandRecoHandle::SetCandSlice(), CandSubShowerSRHandle::SetClusterID(), and CandContext::SetDataIn(). Referenced by RunAlg(). { MSG("SubShowerSR", Msg::kDebug) << "Running FormHalo()" << endl; assert(ch.InheritsFrom("CandSubShowerSRListHandle")); Int_t totUnassigned = allStrips->GetLast()+1; if(totUnassigned==0) return false; MSG("SubShowerSR", Msg::kDebug) << "Forming slice-wise, " << "plane-wise Halo clusters from " << totUnassigned << " unassigned strips" << endl; //loop over slices Int_t nslice = 0; CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator()); while (CandSliceHandle *slice = sliceItr()) { MSG("SubShowerSR",Msg::kDebug) << "Slice " << nslice << endl; //get valid plane ranges for this slice: CandSubShowerSRListHandle &subshowerlist = dynamic_cast<CandSubShowerSRListHandle&>(ch); CandSubShowerSRHandleItr subshowerItr(subshowerlist.GetDaughterIterator()); Int_t plane[500] = {0}; Int_t firstPlane = 500; Int_t lastPlane = 0; while (CandSubShowerSRHandle *subshower = subshowerItr()) { if (*subshower->GetCandSlice()==*slice) { for(int i=subshower->GetBegPlane();i<=subshower->GetEndPlane();i++){ plane[i] = 1; if(i<firstPlane) firstPlane = i; if(i>lastPlane) lastPlane = i; } } } //loop over valid plane ranges and look for gaps > fMaxPlaneGap while(firstPlane<=lastPlane){ Int_t begPlane = 500; Int_t endPlane = 0; Int_t nGaps = 0; for(int i=firstPlane;i<=lastPlane;i++){ if(plane[i]==0) { if(begPlane!=500) nGaps+=1; } else { if(begPlane==500) begPlane = i; endPlane = i; } plane[i] = 0; if(nGaps>fMaxPlaneGap) { break; } } firstPlane = endPlane+1; //use plane range between begPlane and endPlane to form subshowers //loop over views: for(Int_t pv=0;pv<2;pv++){ MSG("SubShowerSR",Msg::kDebug) << "View : " << pv << endl; //first loop through subshower list again and histogram tposvtx //in order to look for peaks TH1F *vtxHist = new TH1F("vtxHist","vtxHist",66,-4.125,4.125); subshowerItr.ResetFirst(); while (CandSubShowerSRHandle *subshower = subshowerItr()) { if (*subshower->GetCandSlice()==*slice && subshower->GetClusterID()!=ClusterType::kXTalk) { if(pv==0 && (subshower->GetPlaneView()==PlaneView::kV || subshower->GetPlaneView()==PlaneView::kY)) continue; else if(pv==1 && (subshower->GetPlaneView()==PlaneView::kU || subshower->GetPlaneView()==PlaneView::kX)) continue; if(subshower->GetBegPlane()>=begPlane && subshower->GetBegPlane()<=endPlane){ if(pv==0) { Float_t weightedCentre = subshower->GetVtxU() + subshower->GetSlope() * PITCHINMETRES * Float_t(subshower->GetEndPlane() - subshower->GetBegPlane() + 1)/2.; if(weightedCentre<-4. || weightedCentre>4.) weightedCentre = subshower->GetVtxU(); vtxHist->Fill(weightedCentre); MSG("SubShowerSR",Msg::kDebug) << "Filling vtxHist with Subshower " << " with tposvtv = " << subshower->GetVtxU() << " and energy = " << subshower->GetEnergy() << " and slope = " << subshower->GetSlope() << " and AvgDev = " << subshower->GetAvgDev() << " and weightedCentre = " << weightedCentre << endl; } else if(pv==1) { Float_t weightedCentre = subshower->GetVtxV() + subshower->GetSlope() * PITCHINMETRES * Float_t(subshower->GetEndPlane() - subshower->GetBegPlane() + 1)/2.; if(weightedCentre<-4. || weightedCentre>4.) weightedCentre = subshower->GetVtxV(); vtxHist->Fill(weightedCentre); MSG("SubShowerSR",Msg::kDebug) << "Filling vtxHist with Subshower " << " with tposvtv = " << subshower->GetVtxV() << " and energy = " << subshower->GetEnergy() << " and slope = " << subshower->GetSlope() << " and AvgDev = " << subshower->GetAvgDev() << " and weightedCentre = " << weightedCentre << endl; } } } } if(false){ TCanvas *can = new TCanvas(); can->cd(1); vtxHist->Draw(); char nomus[256]; sprintf(nomus,"peakPlot_FH_view%i_slice%i.eps",pv,nslice); can->Print(nomus); delete can; } TSpectrum *spec = new TSpectrum(10); spec->Search(vtxHist,1); Int_t nPeaks = spec->GetNPeaks(); Float_t *peakPos = spec->GetPositionX(); Float_t *peakVal = spec->GetPositionY(); //if no peaks are found //then there's no reconstructed subshower //in this slice + plane range + view if(nPeaks==0) { MSG("SubShowerSR",Msg::kDebug) << "No peaks found" << endl; //so take everything and put it in a subshower CandStripHandleItr stripslcItr(slice->GetDaughterIterator()); TObjArray newSubShower; while (CandStripHandle *stp = stripslcItr()) { if(pv==0 && (stp->GetPlaneView()==PlaneView::kV || stp->GetPlaneView()==PlaneView::kY)) continue; else if(pv==1 && (stp->GetPlaneView()==PlaneView::kU || stp->GetPlaneView()==PlaneView::kX)) continue; //is this strip in the subshower plane range? //allow +/- 2 planes in this window if(stp->GetPlane()>endPlane+2 || stp->GetPlane()<begPlane-2) continue; //loop over unassigned strips, and see if it's there for (Int_t i=0; i<allStrips->GetLast()+1; i++) { TObject *tobj = allStrips->At(i); CandStripHandle *unassigned = dynamic_cast<CandStripHandle*>(tobj); if(*unassigned == *stp) { newSubShower.Add(unassigned); allStrips->RemoveAt(i); } } allStrips->Compress(); } //make a new subshower with the unassigned strips in this slice if(newSubShower.GetLast()+1>0){ cxx.SetDataIn(&newSubShower); MSG("SubShowerSR",Msg::kDebug) << "forming SubShowerSR\n"; CandSubShowerSRHandle subshowersrhandle = CandSubShowerSR::MakeCandidate(ah,cxx); if(subshowersrhandle.GetNStrip()!=0) { MSG("SubShowerSR",Msg::kDebug) << "Formed SubShowerSR " << ch.GetNDaughters() << " with " << subshowersrhandle.GetNStrip() << " strips" << endl; subshowersrhandle.SetClusterID(ClusterType::kHalo); subshowersrhandle.SetCandSlice(slice); ch.AddDaughterLink(subshowersrhandle); } } delete spec; delete vtxHist; continue; } //have one or more peak => maybe one or //more shower in this plane range //first find acceptable limits in tpos for a halo //subshower to be formed: Float_t *peakPosLow = new Float_t[nPeaks]; Float_t *peakPosUpp = new Float_t[nPeaks]; MSG("SubShowerSR",Msg::kDebug) << "Number of Peaks found = " << nPeaks << " with positions, heights " << "and bounds: " << endl; for (int i=0;i<nPeaks;i++){ if(i==0) peakPosLow[i] = -4; else peakPosLow[i] = peakPosUpp[i-1]; if(i==nPeaks-1) peakPosUpp[i] = 4; else { if(peakVal[i]+peakVal[i+1] == 0) { MSG("SubShowerSR",Msg::kError) << "peakVal[" << i << "]=" << peakVal[i] << " peakVal[" << i+1 << "]=" << peakVal[i+1] << " - Setting nPeaks=1!" << endl; nPeaks=1; peakPosUpp[0] = 4; break; } peakPosUpp[i] = peakPos[i] + (peakVal[i])*(peakPos[i+1] - peakPos[i]) / ((peakVal[i]) + (peakVal[i+1])); } MSG("SubShowerSR",Msg::kDebug) << peakPos[i] << " " << peakVal[i] << " [" << peakPosLow[i] << "," << peakPosUpp[i] << "]" << endl; } //loop through peaks and assign strips to form halo subshowers //according to acceptable tpos and plane ranges for(int i=0;i<nPeaks;i++){ //loop through slice strips CandStripHandleItr stripslcItr(slice->GetDaughterIterator()); TObjArray newSubShower; //loop over strips in slices in this view: // (have to loop through slice strips and not just // unassigned strips because you cannot get slice number // from strip handle) while (CandStripHandle *stp = stripslcItr()) { //only look at one view at a time: if(pv==0 && (stp->GetPlaneView()==PlaneView::kV || stp->GetPlaneView()==PlaneView::kY)) continue; else if(pv==1 && (stp->GetPlaneView()==PlaneView::kU || stp->GetPlaneView()==PlaneView::kX)) continue; //is this strip in the subshower plane range? //allow +/- 2 planes in this window if(stp->GetPlane()>endPlane+2 || stp->GetPlane()<begPlane-2) continue; //is this strip in the tpos range if(stp->GetTPos()>peakPosUpp[i] || stp->GetTPos()<=peakPosLow[i]) continue; //loop over unassigned strips, and see if it's there for (Int_t i=0; i<allStrips->GetLast()+1; i++) { TObject *tobj = allStrips->At(i); CandStripHandle *unassigned = dynamic_cast<CandStripHandle*>(tobj); if(*unassigned == *stp) { newSubShower.Add(unassigned); allStrips->RemoveAt(i); } } allStrips->Compress(); } //make a new subshower with the unassigned strips in this slice if(newSubShower.GetLast()+1>0){ cxx.SetDataIn(&newSubShower); MSG("SubShowerSR",Msg::kDebug) << "forming SubShowerSR\n"; CandSubShowerSRHandle subshowersrhandle = CandSubShowerSR::MakeCandidate(ah,cxx); if(subshowersrhandle.GetNStrip()!=0) { MSG("SubShowerSR",Msg::kDebug) << "Formed SubShowerSR with " << subshowersrhandle.GetNStrip() << " strips" << endl; subshowersrhandle.SetClusterID(ClusterType::kHalo); subshowersrhandle.SetCandSlice(slice); ch.AddDaughterLink(subshowersrhandle); } } } delete spec; delete vtxHist; delete [] peakPosLow; delete [] peakPosUpp; } } //slice-wise clean up Int_t numberOfCleanUps = 0; while(CleanUp(allStrips,ch,det,slice)) { numberOfCleanUps+=1; //prevent against infinite loop! (which should NEVER happen) if(numberOfCleanUps>1000) break; } nslice++; } //slice-wise clean up Int_t numberOfCleanUps = 0; while(CleanUp(allStrips,ch,det)) { numberOfCleanUps+=1; //prevent against infinite loop! (which should NEVER happen) if(numberOfCleanUps>1000) break; } return true; }

std::vector< window > AlgSubShowerSRList::HoughTransCluster (  std::vector< window >  , Int_t  )  [private]

Definition at line 1742 of file AlgSubShowerSRList.cxx. References BestHough(), MSG, and win. Referenced by TransCluster(). { MSG("SubShowerSR",Msg::kDebug) << "In HoughTransCluster()" << endl; std::vector<window> cluster; std::vector<window> nonCluster; std::vector<window> emptyCluster; //copy all windows into nonCluster: //also get best vertex estimate in first plane Int_t vertexPlane = 500; Double_t vertexPlanePH = 0; Double_t vertexTPos = 0; Double_t vertexTPosWidth = 0; //use these to calculate the density of hits in the long. window //this can be used to decide which Hough method to use Double_t hitDensity = 0; Double_t totStrips = 0; Double_t maxTpos = -5.0; Double_t minTpos = 5.0; Int_t lastPlane = 0; std::vector<window>::iterator wIt = win.begin(); std::vector<window>::iterator wEn = win.end(); while(wIt!=wEn){ nonCluster.push_back(*wIt); if(wIt->plane<vertexPlane) { if(wIt->ph>1.0) { vertexPlane = wIt->plane; vertexPlanePH = wIt->ph; vertexTPos = wIt->phpos; vertexTPosWidth = wIt->phwidth; if(vertexTPosWidth<0.041) vertexTPosWidth=0.041; } else if(vertexPlanePH<1.0){ vertexPlane = wIt->plane; vertexPlanePH = wIt->ph; vertexTPos = wIt->phpos; vertexTPosWidth = wIt->phwidth; if(vertexTPosWidth<0.041) vertexTPosWidth=0.041; } } else if(wIt->plane==vertexPlane){ if(wIt->ph>vertexPlanePH){ vertexPlanePH = wIt->ph; vertexTPos = wIt->phpos; vertexTPosWidth = wIt->phwidth; if(vertexTPosWidth<0.041) vertexTPosWidth=0.041; } } if(wIt->ph>1.0) { if(wIt->plane>lastPlane) lastPlane = wIt->plane; if(wIt->upper_tpos>maxTpos) maxTpos = wIt->upper_tpos; if(wIt->lower_tpos<minTpos) minTpos = wIt->lower_tpos; totStrips += Double_t(wIt->upper - wIt->lower + 1); } wIt++; } //calculate hitDensity: if(lastPlane-vertexPlane+1 != 0) { //multiple planes hitDensity = totStrips/(((maxTpos-minTpos+0.041)/0.041) * Double_t(lastPlane-vertexPlane+1)); } else hitDensity = 0; //single plane MSG("SubShowerSR",Msg::kDebug) << "hitDensity = " << hitDensity << endl; //keep track of some quantities for quality checking at the end: Int_t lowPlane = 500; Int_t uppPlane = 0; Double_t nStrips = 0; Double_t avePH = 0; Double_t aveWidth = 0; //keep looping until we have all windows Bool_t keepGoing = true; Int_t nLoops = 0; while(keepGoing && nLoops<5) { MSG("SubShowerSR",Msg::kDebug) << "Running BestHough(): Loop number " << nLoops << endl; //remember number of strips in cluster after last loop //once this value does not change, then stop looping Double_t oldNStrips = nStrips; //Get Hough m,c: Double_t *MCV = new Double_t[14]; MCV[0] = 0; MCV[1] = 0; MCV[2] = 0; MCV[3] = 0; MCV[4] = 0; MCV[5] = 0; MCV[6] = 0; MCV[7] = 0; MCV[8] = 0; MCV[9] = 0; MCV[10] = 0; MCV[11] = 0; MCV[12] = 0; MCV[13] = 0; Int_t badGradient[4] = {0}; //To decide which Hough to use first Int_t offset = 5; //by default PH first if(hitDensity<0.05) offset = 0; //unless hit density <0.05 of max possible //If we are on the first loop, use all windows to get angle //If we are not on first loop, just use previously //selected windows to get better determination of angle if(cluster.size()==0){ MSG("SubShowerSR",Msg::kDebug) << "cluster.size()=0" << endl; if(BestHough(win,MCV,false)) { if(true){ //verify that best hough passes reasonably close to vertex: MSG("SubShowerSR",Msg::kDebug) << "vertexPlane = " << vertexPlane << endl; MSG("SubShowerSR",Msg::kDebug) << "vertexTPos = " << vertexTPos << endl; MSG("SubShowerSR",Msg::kDebug) << "vertexTPosWidth = " << vertexTPosWidth << endl; double vertexY_simple = (vertexPlane - MCV[2])*MCV[0] + MCV[1]; double vertexY_simple_peak = (vertexPlane - MCV[2])*MCV[10] + MCV[11]; MSG("SubShowerSR",Msg::kDebug) << "vertexY_simple = " << vertexY_simple << endl; if( vertexY_simple < vertexTPos - vertexTPosWidth || vertexY_simple > vertexTPos + vertexTPosWidth || MCV[0]==0 ){ offset = 5; badGradient[0] = 1; } if( (vertexY_simple_peak < vertexTPos - vertexTPosWidth) || (vertexY_simple_peak > vertexTPos + vertexTPosWidth) || MCV[10]==0 ){ badGradient[1] = 1; } double vertexY_ph = (vertexPlane - MCV[7])*MCV[5] + MCV[6]; double vertexY_ph_peak = (vertexPlane - MCV[7])*MCV[12] + MCV[13]; MSG("SubShowerSR",Msg::kDebug) << "vertexY_ph = " << vertexY_ph << endl; if( vertexY_ph < vertexTPos - vertexTPosWidth || vertexY_ph > vertexTPos + vertexTPosWidth || MCV[5]==0 ) { offset = 0; //use simple Hough if this is no good badGradient[2] = 1; } if( vertexY_ph_peak < vertexTPos - vertexTPosWidth || vertexY_ph_peak > vertexTPos + vertexTPosWidth || MCV[12]==0 ) { badGradient[3] = 1; } } } } else { MSG("SubShowerSR",Msg::kDebug) << "cluster.size()=" << cluster.size() << endl; BestHough(cluster,MCV,false); } MSG("SubShowerSR",Msg::kDebug) << "\nHough gradient, intercept, plane vertex: " << MCV[0] << "+/-" << MCV[3] << " " << MCV[1] << "+/-" << MCV[4] << " " << MCV[2] << "\nPH Hough gradient, intercept, plane vertex: " << MCV[5] << "+/-" << MCV[8] << " " << MCV[6] << "+/-" << MCV[9] << " " << MCV[7] << "\n Hough Max Bin Gradient and Intercept Values: " << MCV[10] << " " << MCV[11] << " " << MCV[12] << " " << MCV[13] << endl; //if we've learnt nothing from first run of BestHough //stop now and try other clustering algorithm if(badGradient[0]==1 && badGradient[1]==1 && badGradient[2]==1 && badGradient[3]==1 && cluster.size()==0) { delete [] MCV; return emptyCluster; } if(true){ cluster.clear(); avePH = 0; aveWidth = 0; nStrips = 0; lowPlane = 500; uppPlane = 0; std::vector<int> tempcluster[4]; std::vector<window>::iterator winIter = nonCluster.begin(); std::vector<window>::iterator winEnd = nonCluster.end(); Int_t counter = 0; while(winIter!=winEnd){ double y[4] = {}; double y_upp[4] = {}; double y_low[4] = {}; y[0] = (winIter->plane - MCV[2])*MCV[0] + MCV[1]; y[1] = (winIter->plane - MCV[2])*MCV[10] + MCV[11]; y[2] = (winIter->plane - MCV[7])*MCV[5] + MCV[6]; y[3] = (winIter->plane - MCV[7])*MCV[12] + MCV[13]; if(cluster.size()==0) { for(int ii=0;ii<4;ii++) { y_upp[ii] = y[ii] + STRIPWIDTHINMETRES; y_low[ii] = y[ii] - STRIPWIDTHINMETRES; } } else { for(int ii=0;ii<4;ii++) { y_upp[ii] = y[ii] + MCV[4+int(ii/2)*5]; y_low[ii] = y[ii] - MCV[4+int(ii/2)*5]; } } for(int ii=0;ii<4;ii++){ if( ( winIter->upper_tpos>=y[ii] && winIter->lower_tpos<=y[ii] ) || ( winIter->phpos<=y_upp[ii] && winIter->phpos>=y_low[ii] ) ) { if(badGradient[ii]==0) { tempcluster[ii].push_back(counter); } } } winIter++; counter++; } Int_t best_cluster = -1; UInt_t max_size = 0; for(int ii=0;ii<4;ii++){ if(tempcluster[ii].size()>max_size) { max_size = tempcluster[ii].size(); best_cluster = ii; } } if(best_cluster>=0 && best_cluster<=3) { std::vector<window>::iterator winIter = nonCluster.begin(); std::vector<window>::iterator winEnd = nonCluster.end(); std::vector<int>::iterator tempIter = tempcluster[best_cluster].begin(); std::vector<int>::iterator tempEnd = tempcluster[best_cluster].end(); counter = 0; while(winIter!=winEnd){ if(tempIter!=tempEnd && (*tempIter) == counter) { cluster.push_back(*winIter); avePH += winIter->ph; aveWidth += winIter->upper_tpos - winIter->lower_tpos; nStrips += Double_t(winIter->upper - winIter->lower + 1); if(winIter->plane<lowPlane) lowPlane = winIter->plane; if(winIter->plane>uppPlane) uppPlane = winIter->plane; tempIter++; } winIter++; counter++; } } } else { cluster.clear(); avePH = 0; aveWidth = 0; nStrips = 0; lowPlane = 500; uppPlane = 0; Bool_t stop = false; while(!stop && offset>=0){ std::vector<window>::iterator winIter = nonCluster.begin(); std::vector<window>::iterator winEnd = nonCluster.end(); MSG("SubShowerSR",Msg::kVerbose) << "Adding to Proto SubShower: " << endl; while(winIter!=winEnd){ double y = (winIter->plane - MCV[2+offset])*MCV[offset] + MCV[1+offset]; double y_upp = (winIter->plane - MCV[2+offset])*MCV[offset] + MCV[1+offset]; double y_low = (winIter->plane - MCV[2+offset])*MCV[offset] + MCV[1+offset]; if(cluster.size()==0) { y_upp += 0.041; y_low -= 0.041; } else { y_upp += MCV[4+offset]; y_low -= MCV[4+offset]; } MSG("SubShowerSR",Msg::kVerbose) << "Plane: " << winIter->plane << endl; MSG("SubShowerSR",Msg::kVerbose) << "y, upp_tpos, low_tpos: " << y << " " << winIter->upper_tpos << " " << winIter->lower_tpos << endl; MSG("SubShowerSR",Msg::kVerbose) << "phpos, y_upp, y_low: " << winIter->phpos << " " << y_upp << " " << y_low << endl; if( ( winIter->upper_tpos>=y && winIter->lower_tpos<=y ) || ( winIter->phpos<=y_upp && winIter->phpos>=y_low ) ) { MSG("SubShowerSR",Msg::kVerbose) << "Passed!" << endl; cluster.push_back(*winIter); avePH += winIter->ph; aveWidth += winIter->upper_tpos - winIter->lower_tpos; nStrips += Double_t(winIter->upper - winIter->lower + 1); if(winIter->plane<lowPlane) lowPlane = winIter->plane; if(winIter->plane>uppPlane) uppPlane = winIter->plane; } winIter++; } offset-=5; if(cluster.size()!=0) stop = true; } } delete [] MCV; MSG("SubShowerSR",Msg::kDebug) << "Proto-SubShower currently has " << (uppPlane-lowPlane+2)/2 << " planes and " << nStrips << " strips" << endl; if(nStrips == oldNStrips) keepGoing = false; nLoops+=1; } if(nStrips<=0) return emptyCluster; MSG("SubShowerSR",Msg::kDebug) << "Returning Proto-SubShower from HoughTransCluster()" << endl; if( det==1 && uppPlane>249 && lowPlane<249 ) return cluster; avePH/=nStrips; if( ((uppPlane - lowPlane + 2)/2)<3 && avePH<1. ) return emptyCluster; return cluster; }

Bool_t AlgSubShowerSRList::MergeCluster (  TObjArray *  , CandHandle &  , Int_t  , Bool_t  test = false )  [private]

Definition at line 2252 of file AlgSubShowerSRList.cxx. References CandStripHandle::DupHandle(), CandSubShowerSRHandle::GetAvgDev(), CandRecoHandle::GetBegPlane(), CandRecoHandle::GetCharge(), CandStripHandle::GetCharge(), CandSubShowerSRHandle::GetClusterID(), CandHandle::GetDaughterIterator(), CandRecoHandle::GetEndPlane(), CandSubShowerSRHandle::GetEnergy(), Calibrator::GetMIP(), CandHandle::GetNDaughters(), CandStripHandle::GetPlane(), CandSubShowerSRHandle::GetPlaneView(), CandStripHandle::GetPlaneView(), CandSubShowerSRHandle::GetSlope(), CandStripHandle::GetTPos(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), CandStripHandle::GetZPos(), Calibrator::Instance(), MSG, PITCHINMETRES, CandHandle::RemoveDaughter(), and STRIPWIDTHINMETRES. Referenced by RunAlg(). { MSG("SubShowerSR", Msg::kDebug) << "Running MergeCluster()" << endl; Int_t totStp = newSubShower->GetLast()+1; if(totStp==0) return false; Calibrator& calibrator=Calibrator::Instance(); assert(ch.InheritsFrom("CandSubShowerSRListHandle")); CandSubShowerSRListHandle &csslh = dynamic_cast<CandSubShowerSRListHandle &>(ch); if(csslh.GetNDaughters()==0) return false; //Get planeview for this merge test TObject *firstObj = newSubShower->At(0); assert(firstObj->InheritsFrom("CandStripHandle")); CandStripHandle *firstStp = dynamic_cast<CandStripHandle*>(firstObj); PlaneView::PlaneView_t pv = firstStp->GetPlaneView(); MSG("SubShowerSR", Msg::kDebug) << "PlaneView = " << pv << endl; //check that PH of new subshower is less than some //fraction of the primary subshower in the view Float_t PrimarySubShowerEnergy = 0; Float_t NewSubShowerEnergy = 0; Int_t firstPlaneInNewSubShower = 500; Int_t lastPlaneInNewSubShower = 0; for (Int_t i=0; i<totStp; i++) { TObject *tobj = newSubShower->At(i); CandStripHandle *nssStrip = dynamic_cast<CandStripHandle*>(tobj); NewSubShowerEnergy += calibrator.GetMIP(nssStrip->GetCharge(CalDigitType::kSigCorr)); Int_t thePlane = nssStrip->GetPlane(); if(thePlane<firstPlaneInNewSubShower) firstPlaneInNewSubShower = thePlane; if(thePlane>lastPlaneInNewSubShower) lastPlaneInNewSubShower = thePlane; } //calculate the number of planes in this newsubshower: Float_t nPlanesInNewSubShower = Float_t(lastPlaneInNewSubShower-firstPlaneInNewSubShower+2)/2.; MSG("SubShowerSR", Msg::kDebug) << "First Plane in NewSubShower = " << firstPlaneInNewSubShower << " Last Plane in NewSubShower = " << lastPlaneInNewSubShower << " Number of Planes in New SubShower " << nPlanesInNewSubShower << endl; Double_t mip2gev; CandSubShowerSRHandleItr PHsubshowerItr(csslh.GetDaughterIterator()); while (CandSubShowerSRHandle *subshower = PHsubshowerItr()) { if(subshower->GetPlaneView()!=pv) continue; mip2gev = subshower->GetEnergy() / calibrator.GetMIP(subshower->GetCharge(CalStripType::kSigCorr)); if(subshower->GetEnergy()>PrimarySubShowerEnergy) PrimarySubShowerEnergy=subshower->GetEnergy(); } NewSubShowerEnergy*=mip2gev; if(PrimarySubShowerEnergy<=0) return false; if(NewSubShowerEnergy/PrimarySubShowerEnergy>fMergeFracEnergyLimit) return false; //keep track of the following to decide whether or not to merge: //1) number of shared boundary planes: //2) number of planes which have strips nearby to another subshower core: //3) charge asymmetry of newsubshower around axis of other subshowers: // (if this is not high => newsubshower is being formed // from outliers of another subshower) //4) summed distance of each hit in newsubshower to the closest hit in another subshower: // (only used when newsubshower has <=2 planes) Float_t *FracBoundaryPlanes = new Float_t[csslh.GetNDaughters()]; Float_t *FracWithinPlanes = new Float_t[csslh.GetNDaughters()]; Float_t *FracDiffCross = new Float_t[csslh.GetNDaughters()]; Float_t *SumDistanceFromNewSubShower = new Float_t[csslh.GetNDaughters()]; for(int i=0;i<csslh.GetNDaughters();i++) { FracBoundaryPlanes[i] = 0; FracWithinPlanes[i] = 0; FracDiffCross[i] = 1.; SumDistanceFromNewSubShower[i] = 0.; } //loop through existing subshowers Int_t ssCount = 0; CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator()); while (CandSubShowerSRHandle *subshower = subshowerItr()) { ssCount++; //increment early because of continue statements, later use ssCount-1 //only look at the right planeview subshowers if(subshower->GetPlaneView()!=pv) { SumDistanceFromNewSubShower[ssCount-1] = 99999.; continue; } MSG("SubShowerSR", Msg::kVerbose) << "SubShower " << ssCount-1 << " in the same planeview" << endl; //get subshower vertex and angle: Double_t tposVtx = 0; if(pv==PlaneView::kU || pv==PlaneView::kX) tposVtx = subshower->GetVtxU(); else if(pv==PlaneView::kV || pv==PlaneView::kY) tposVtx = subshower->GetVtxV(); Double_t zposVtx = subshower->GetVtxZ(); Float_t slope = subshower->GetSlope(); Float_t avgdev = subshower->GetAvgDev(); ClusterType::ClusterType_t id = subshower->GetClusterID(); Bool_t stripVeto = false; Int_t PlaneChecker[500] = {0}; Double_t PlaneTot[500] = {0}; Double_t PlaneWithin[500] = {0}; Double_t PlaneUp = 0; Double_t PlaneDown = 0; //shortest distance between any strip in newsubshower and any strip in existing subshower: Float_t minimumDistanceFromNewSubShower = 999.; //perform checks to see if new potential cluster can be merged with an existing cluster for (Int_t i=0; i<totStp; i++) { TObject *tobj = newSubShower->At(i); CandStripHandle *nssStrip = dynamic_cast<CandStripHandle*>(tobj); // if plane of any strip is more than 1 plane (of same type) away // don't try to merge with this subshower unless new subshower contains <=2 planes // (also beware of SM boundary in FarDet) Int_t thePlane = nssStrip->GetPlane(); Float_t theCharge = calibrator.GetMIP(nssStrip->GetCharge(CalDigitType::kSigCorr)); if(thePlane<subshower->GetBegPlane()-2 || thePlane>subshower->GetEndPlane()+2 || (det==1 && ( (subshower->GetBegPlane()>249 && thePlane<249) || (subshower->GetEndPlane()<249 && thePlane>249) ) ) ) { stripVeto = true; } else MSG("SubShowerSR", Msg::kVerbose) << "New Strip " << i << " in the same plane range" << endl; //if we got a stripVeto above, //and this NewSubShower is longer than 2 planes //or the existing SubShower is not of type XTalk //don't merge, break here: if(stripVeto) { if(nPlanesInNewSubShower>2) break; if(id!=ClusterType::kXTalk) { SumDistanceFromNewSubShower[ssCount-1] = 99999.; break; } } //loop through strips in subshower //check for proximity to other strips in subshower (on a plane-by-plane basis) Float_t theStrip_tpos = nssStrip->GetTPos(); //remember closest distance from nssStrip to any strip in subshower: Float_t closestDistance = 999.; CandStripHandleItr stripItr(subshower->GetDaughterIterator()); while(CandStripHandle *stp = stripItr()){ //check we're on the same plane //or on the begin plane if the nssStrip plane < begin plane //or on the end plane if the nssStrip plane > end plane if((stp->GetPlane() == thePlane) || (thePlane < subshower->GetBegPlane() && stp->GetPlane() == subshower->GetBegPlane()) || (thePlane > subshower->GetEndPlane() && stp->GetPlane() == subshower->GetEndPlane())) { if(TMath::Abs(theStrip_tpos-stp->GetTPos())<=2.5*STRIPWIDTHINMETRES) { //allow +/- 2.5 strips PlaneChecker[thePlane] = 1; //this plane has a shared boundary } if(TMath::Abs(theStrip_tpos-stp->GetTPos())<closestDistance) { closestDistance = TMath::Abs(theStrip_tpos-stp->GetTPos()); } } if(stripVeto) { //newsubshower must be short (<2 planes) //this strip is not adjacent to any existing subshower //calculate the closestDistance anyway, but factor in delta(plane): Float_t testDistance = TMath::Abs(theStrip_tpos-stp->GetTPos()) + PITCHINMETRES*TMath::Abs(Float_t(thePlane-stp->GetPlane())); if(testDistance < closestDistance) closestDistance = testDistance; } } if(closestDistance<minimumDistanceFromNewSubShower) minimumDistanceFromNewSubShower = closestDistance; SumDistanceFromNewSubShower[ssCount-1] += closestDistance; //calculate pH weighted distance from cluster centre //get strip position and charge: Float_t tpos = nssStrip->GetTPos(); Float_t zpos = nssStrip->GetZPos(); Double_t StripDistance = 99999; if(avgdev!=0.) { StripDistance = (tpos-(tposVtx+slope*(zpos-zposVtx))* TMath::Cos(TMath::ATan(slope)))*theCharge/avgdev; } else MSG("SubShowerSR",Msg::kWarning) << "SubShower " << ssCount-1 << " has AvgDev = 0!" << endl; if(TMath::Abs(StripDistance)<=1.) PlaneWithin[thePlane]+=theCharge; if(StripDistance>0) PlaneUp+=theCharge; if(StripDistance<0) PlaneDown+=theCharge; PlaneTot[thePlane]+=theCharge; } if(stripVeto) continue; //If this is true, Must be comparing a subshower with itself if(SumDistanceFromNewSubShower[ssCount-1] == 0.000000) { SumDistanceFromNewSubShower[ssCount-1] = 99999; continue; } MSG("SubShowerSR",Msg::kVerbose) << "SumDistanceFromNewSubShower for SubShower " << ssCount-1 << " = " << SumDistanceFromNewSubShower[ssCount-1] << endl; //fill the arrays for this subshower: Float_t n0planes = 0.; for(int i = firstPlaneInNewSubShower;i<=lastPlaneInNewSubShower;i+=2){ FracBoundaryPlanes[ssCount-1] += PlaneChecker[i]; if (PlaneTot[i]>0){ FracWithinPlanes[ssCount-1] += PlaneWithin[i]/PlaneTot[i]; } else n0planes++; } FracBoundaryPlanes[ssCount-1] /= nPlanesInNewSubShower-n0planes; FracWithinPlanes[ssCount-1] /= nPlanesInNewSubShower-n0planes; if (PlaneUp*PlaneDown>0) { FracDiffCross[ssCount-1] = TMath::Abs(PlaneUp-PlaneDown)/(PlaneUp+PlaneDown); } //lastly we need to worry about a special case: //if newSubShower is in the same plane range, and the number of planes is <=2 //but the closest strip to the current subshower is >2.5 strips away and the //cluster id is not kXTalk, then don't merge with this subshower if(!stripVeto && nPlanesInNewSubShower<=2 && minimumDistanceFromNewSubShower>2.5*STRIPWIDTHINMETRES && id!=ClusterType::kXTalk) SumDistanceFromNewSubShower[ssCount-1] = 99999.; } Float_t BestSumDistance = 9999.; Int_t closestSubShowerSumDistance = -1; Float_t BestFracBoundary = fMergeBestFracBoundary; Int_t closestSubShowerBoundary = -1; Float_t BestFracWithin = fMergeBestFracWithin; Int_t closestSubShowerWithin = -1; Float_t BestFracCross = fMergeBestFracCross; Int_t closestSubShowerCross = -1; //Float_t BestFracBoundaryHalf = BestFracBoundary/2.; //Float_t BestFracWithinHalf = BestFracWithin/2.; //Int_t closestSubShowerBoundWithin = -1; Int_t closestSubShower = -1; for(int i=0;i<csslh.GetNDaughters();i++){ if(false&&nPlanesInNewSubShower<=2){ if(SumDistanceFromNewSubShower[i]<BestSumDistance) { BestSumDistance = SumDistanceFromNewSubShower[i]; closestSubShowerSumDistance = i; } } else { if(FracBoundaryPlanes[i]>BestFracBoundary){ BestFracBoundary = FracBoundaryPlanes[i]; closestSubShowerBoundary = i; } if(FracWithinPlanes[i]>BestFracWithin) { BestFracWithin = FracWithinPlanes[i]; closestSubShowerWithin = i; } //if(FracBoundaryPlanes[i]>BestFracBoundaryHalf && // FracWithinPlanes[i]>BestFracWithinHalf) { //BestFracBoundaryHalf = FracBoundaryPlanes[i]; //BestFracWithinHalf = FracWithinPlanes[i]; //closestSubShowerBoundWithin = i; //} if(FracDiffCross[i]<BestFracCross) { BestFracCross = FracDiffCross[i]; closestSubShowerCross = i; } } } delete [] FracBoundaryPlanes; delete [] FracWithinPlanes; delete [] FracDiffCross; delete [] SumDistanceFromNewSubShower; if(closestSubShowerCross>=0) closestSubShower = closestSubShowerCross; else if(closestSubShowerSumDistance>=0) closestSubShower = closestSubShowerSumDistance; else if(closestSubShowerBoundary>=0) closestSubShower = closestSubShowerBoundary; else if(closestSubShowerWithin>=0) closestSubShower = closestSubShowerWithin; //else if(closestSubShowerBoundWithin>=0) closestSubShower = closestSubShowerBoundWithin; if(closestSubShower!=-1){ CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator()); ssCount = 0; while(CandSubShowerSRHandle *cssh = subshowerItr()){ if(ssCount==closestSubShower){ MSG("SubShowerSR", Msg::kDebug) << "Found Sub Shower to merge "<<closestSubShower<<endl; if(test) return true; CandStripHandleItr stripItr(cssh->GetDaughterIterator()); while(CandStripHandle *stp = stripItr()) newSubShower->Add(stp->DupHandle()); csslh.RemoveDaughter(cssh); MSG("SubShowerSR", Msg::kDebug) << "Old Sub Shower removed"<<endl; return true; } ssCount+=1; } } return false; //couldn't merge }

void AlgSubShowerSRList::RunAlg (  AlgConfig &  ac, CandHandle &  ch, CandContext &  cx )  [virtual]

Implements AlgBase. Definition at line 70 of file AlgSubShowerSRList.cxx. References CandHandle::AddDaughterLink(), CleanUp(), CandStripHandle::DupHandle(), FindCluster(), fLongWindowMipCut, fMaxPlaneGap, fMergeBestFracBoundary, fMergeBestFracCross, fMergeBestFracWithin, fMergeFracEnergyLimit, FormHalo(), Registry::Get(), AlgFactory::GetAlgHandle(), CandClusterHandle::GetBegPlane(), CandRecoHandle::GetBegPlane(), CandContext::GetCandRecord(), CandClusterHandle::GetCandSlice(), CandRecoHandle::GetCandSlice(), CandStripHandle::GetCharge(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), Registry::GetDouble(), CandClusterHandle::GetEndPlane(), CandRecoHandle::GetEndPlane(), AlgFactory::GetInstance(), Registry::GetInt(), Calibrator::GetMIP(), CandContext::GetMom(), CandHandle::GetNDaughters(), CandRecoHandle::GetNStrip(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), RecMinos::GetVldContext(), Calibrator::Instance(), CandSubShowerSR::MakeCandidate(), MergeCluster(), MSG, CandHandle::RemoveDaughter(), CandRecoHandle::SetCandSlice(), and TestOverLap(). { MSG("SubShowerSR", Msg::kDebug) << "Starting AlgSubShowerSRList::RunAlg()" << endl; assert(cx.GetDataIn()); if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) { return; } const CandSliceListHandle *slicelist = 0; const CandClusterListHandle *clusterlist = 0; const CandTrackListHandle *tracklist = 0; const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn()); for (Int_t i=0; i<=cxin->GetLast(); i++) { TObject *tobj = cxin->At(i); MSG("SubShowerSR", Msg::kDebug) << "cxin->At(" << i << ")"<< endl; if (tobj->InheritsFrom("CandSliceListHandle")) { slicelist = dynamic_cast<CandSliceListHandle*>(tobj); MSG("SubShowerSR", Msg::kDebug) << "Got SliceList" << endl; } if (tobj->InheritsFrom("CandClusterListHandle")) { clusterlist = dynamic_cast<CandClusterListHandle*>(tobj); MSG("SubShowerSR", Msg::kDebug) << "Got ClusterList" << endl; } if (tobj->InheritsFrom("CandTrackListHandle")) { tracklist = dynamic_cast<CandTrackListHandle*>(tobj); MSG("SubShowerSR", Msg::kDebug) << "Got TrackList" << endl; } } if (!slicelist || !clusterlist || !tracklist) { MSG("SubShowerSR",Msg::kInfo) << "CandSliceListHandle or " << "CandClusterListHandle or " << "CandTrackListHandle missing\n"; return; } MSG("SubShowerSR", Msg::kDebug) << "Creating CandContext" << endl; // Create Candcontext CandContext cxx(this,cx.GetMom()); //get config for AlgSubShowerSRList const char *charSubShowerSRAlgConfig = 0; ac.Get("SubShowerSRAlgConfig",charSubShowerSRAlgConfig); fMergeFracEnergyLimit = ac.GetDouble("MergeFracEnergyLimit"); fMergeBestFracBoundary = ac.GetDouble("MergeBestFracBoundary"); fMergeBestFracWithin = ac.GetDouble("MergeBestFracWithin"); fMergeBestFracCross = ac.GetDouble("MergeBestFracCross"); fLongWindowMipCut = ac.GetDouble("LongWindowMipCut"); fMaxPlaneGap = ac.GetInt("MaxPlaneGap"); Int_t minNTrkOnlyPlanesToUseTrackStrips = ac.GetInt("MinNTrkOnlyPlanesToUseTrackStrips"); Int_t maxPlanesFromVtxToUseTrackStrips = ac.GetInt("MaxPlanesFromVertexToUseTrackStrips"); Double_t minMIPsToUseTrackStrips = ac.GetDouble("MinMIPsToUseTrackStrips"); //Get singleton instance of AlgFactory AlgFactory &af = AlgFactory::GetInstance(); AlgHandle ah = af.GetAlgHandle("AlgSubShowerSR",charSubShowerSRAlgConfig); const CandRecord *candrec = cx.GetCandRecord(); assert(candrec); const VldContext *vldcptr = candrec->GetVldContext(); assert(vldcptr); VldContext vldc = *vldcptr; UgliGeomHandle ugh(vldc); Int_t detector = 0; if(vldc.GetDetector()==DetectorType::kFar) detector = 1; Calibrator& calibrator = Calibrator::Instance(); //keep track of all unassigned strips in all slices and try to add them later TObjArray unassignedStrips; //loop over all slices Int_t sliceCounter = 0; CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator()); while (CandSliceHandle *slice = sliceItr()) { MSG("SubShowerSR",Msg::kVerbose) << "Slice " << sliceCounter << endl; //also keep track of cluster planes and track planes Int_t cluPln[500] = {0}; Int_t trkPln[500] = {0}; Int_t nTrkOnlyPlanes = 0; if(tracklist!=0) { //loop through all tracks in track list CandTrackHandleItr trackItr(tracklist->GetDaughterIterator()); while (CandTrackHandle *track = trackItr()){ //if this track is in the current slice then... if (*track->GetCandSlice()==*slice) { for(Int_t i=track->GetBegPlane();i<=track->GetEndPlane();i++){ trkPln[i] += 1; //if trkPln > 1 => multiple tracks in plane nTrkOnlyPlanes += 1; } } } } if(clusterlist!=0) { //loop through all clusters in cluster list CandClusterHandleItr clusterItr(clusterlist->GetDaughterIterator()); while (CandClusterHandle *cluster = clusterItr()){ //if this cluster is in the current slice then... if (*cluster->GetCandSlice()==*slice) { for(Int_t i=cluster->GetBegPlane();i<=cluster->GetEndPlane();i+=2){ cluPln[i] += 1; //if cluPln > 1 => multiple clusters in plane if(trkPln[i]>0) nTrkOnlyPlanes-=1; if(detector==1){//account for SM gap if(i==247) i = 250; if(i==248) i = 251; } } } } } //for each slice, make an Object Array to hold all the //non track strips in the slice TObjArray allUStrips; TObjArray allVStrips; //loop through the strips in the slice CandStripHandleItr stripslcItr(slice->GetDaughterIterator()); MSG("SubShowerSR",Msg::kVerbose) << "About to loop through slice strips" << endl; while (CandStripHandle *stp =stripslcItr()) { //while over slice stp bool gotStrip = false; MSG("SubShowerSR",Msg::kVerbose) << "PlaneView " << stp->GetPlaneView() << endl; CandClusterHandleItr clusterItr(clusterlist->GetDaughterIterator()); Int_t clusterCounter = 0; while(CandClusterHandle *cluster = clusterItr()){ MSG("SubShowerSR",Msg::kVerbose) << "Looking for stp match in cluster " << clusterCounter << endl; if(*cluster->GetCandSlice()==*slice){ CandStripHandleItr stripcluItr(cluster->GetDaughterIterator()); while (CandStripHandle *stpclu =stripcluItr()){ if(*stp==*stpclu){ gotStrip = true; //add strip to object array then break out of loop if(stp->GetPlaneView()==PlaneView::kU || stp->GetPlaneView()==PlaneView::kX) { allUStrips.Add(stp); MSG("SubShowerSR",Msg::kVerbose) << "Found match! U strip" << endl; } else if(stp->GetPlaneView()==PlaneView::kV || stp->GetPlaneView()==PlaneView::kY) { allVStrips.Add(stp); MSG("SubShowerSR",Msg::kVerbose) << "Found match! V strip" << endl; } else MSG("SubShowerSR",Msg::kError) << "Unknown PlaneView! Not adding strip to initial ObjArray" << endl; break; } } } if(gotStrip) break; clusterCounter++; } //if this is an ND strip in spectrometer, don't add to array if(!gotStrip && detector==0 && stp->GetPlane()>=121) gotStrip=true; if(!gotStrip){ //if !gotStrip bool inTrack = false; Int_t nPlanesFromBegPlane = 0; if(tracklist!=0) { MSG("SubShowerSR",Msg::kVerbose) << "No matched strip found in clusters, looking in tracks" << endl; //loop through all tracks in track list CandTrackHandleItr trackItr(tracklist->GetDaughterIterator()); Int_t trackCounter = 0; while (CandTrackHandle *track = trackItr()){ MSG("SubShowerSR",Msg::kVerbose) << "Looking for match of slc stp with trk stp in track " << trackCounter << endl; //if this track is in the current slice then... if (*track->GetCandSlice()==*slice) { //...loop over all strips in this track CandStripHandleItr striptrkItr(track->GetDaughterIterator()); while (CandStripHandle *stptrk =striptrkItr()){ //if this strip is in this track then... if (*stp == *stptrk) { MSG("SubShowerSR",Msg::kVerbose) << "Found match! This strip is in a track!" << endl; //...break and move on to the next strip in the slice inTrack = true; nPlanesFromBegPlane = TMath::Abs(stp->GetPlane()-track->GetBegPlane()); break; } } } if(inTrack) break; trackCounter++; } } MSG("SubShowerSR",Msg::kVerbose) << "nTrkOnlyPlanes = " << nTrkOnlyPlanes << endl; if(inTrack && nTrkOnlyPlanes<minNTrkOnlyPlanesToUseTrackStrips) { MSG("SubShowerSR",Msg::kDebug) << "nTrkOnlyPlanes<" << minNTrkOnlyPlanesToUseTrackStrips << " - going to include this strip " << "even though it's in a track." << endl; inTrack = false; //ignore tracks that are buried in showers } else if(inTrack && nPlanesFromBegPlane < maxPlanesFromVtxToUseTrackStrips) { //near vertex if(calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr)) > minMIPsToUseTrackStrips) { MSG("SubShowerSR",Msg::kDebug) << "nPlanesFromBegPlane < " << maxPlanesFromVtxToUseTrackStrips << " for this strip and ph = " << calibrator.GetMIP(stp->GetCharge(CalDigitType::kSigCorr)) << " mips. " << "Going to include it even though it's in a track." << endl; inTrack = false; } } //then: add it to our array if(!inTrack){ MSG("SubShowerSR",Msg::kVerbose) << "Strip not found in track! "; if(stp->GetPlaneView()==PlaneView::kU || stp->GetPlaneView()==PlaneView::kX) { allUStrips.Add(stp); MSG("SubShowerSR",Msg::kVerbose) << "Adding U strip" << endl; } else if(stp->GetPlaneView()==PlaneView::kV || stp->GetPlaneView()==PlaneView::kY) { allVStrips.Add(stp); MSG("SubShowerSR",Msg::kVerbose) << "Adding V strip" << endl; } else MSG("SubShowerSR",Msg::kError) << "Unknown PlaneView! Not adding strip to initial ObjArray" << endl; } } } MSG("SubShowerSR",Msg::kDebug) << "Found total of " << allUStrips.GetLast()+1 << " U strips" << " and " << allVStrips.GetLast()+1 << " V strips" << endl; //In here put code to pull out the strips belonging to the SubShowers: bool keepGoingU = true; bool keepGoingV = true; while(keepGoingU || keepGoingV){ if(keepGoingU) MSG("SubShowerSR",Msg::kDebug) << "Still have " << allUStrips.GetLast()+1 << " U strips to process" << endl; else if(keepGoingV) MSG("SubShowerSR",Msg::kDebug) << "Still have " << allVStrips.GetLast()+1 << " V strips to process" << endl; if(allUStrips.GetLast()+1<=0) keepGoingU = false; if(allVStrips.GetLast()+1<=0) keepGoingV = false; TObjArray newSubShower; if(keepGoingU){ //do U first this->FindCluster(&allUStrips,&newSubShower,detector); while(this->TestOverLap(&newSubShower,ch,slice)); MSG("SubShowerSR",Msg::kDebug) << "Added " << newSubShower.GetLast()+1 << " U strips to newSubShower" << endl; } else if(keepGoingV) { //then V strips TObjArray protoSubShower; this->FindCluster(&allVStrips,&newSubShower,detector); while(this->TestOverLap(&newSubShower,ch,slice)); MSG("SubShowerSR",Msg::kDebug) << "Added " << newSubShower.GetLast()+1 << " V strips to newSubShower" << endl; } //Everytime you have a newSubShower with all the strips //you want to add do this cxx.SetDataIn(&newSubShower); MSG("SubShowerSR",Msg::kDebug) << "forming SubShowerSR\n"; CandSubShowerSRHandle subshowersrhandle = CandSubShowerSR::MakeCandidate(ah,cxx); if(subshowersrhandle.GetNStrip()!=0) { MSG("SubShowerSR",Msg::kDebug) << "Formed SubShowerSR with " << subshowersrhandle.GetNStrip() << " strips" << endl; subshowersrhandle.SetCandSlice(slice); ch.AddDaughterLink(subshowersrhandle); CandStripHandleItr stripItr(subshowersrhandle.GetDaughterIterator()); int *tmpUArray = new int[allUStrips.GetLast()+1]; for(int i=0;i<allUStrips.GetLast()+1;i++) tmpUArray[i] = 0; int *tmpVArray = new int[allVStrips.GetLast()+1]; for(int i=0;i<allVStrips.GetLast()+1;i++) tmpVArray[i] = 0; while(CandStripHandle *stp = stripItr()){ if(stp->GetPlaneView()==PlaneView::kU || stp->GetPlaneView()==PlaneView::kX){ for (int i=0; i<=allUStrips.GetLast(); i++) { TObject *tobj = allUStrips.At(i); assert(tobj->InheritsFrom("CandStripHandle")); CandStripHandle *striphandle = dynamic_cast<CandStripHandle*>(tobj); if(*stp==*striphandle) { tmpUArray[i] = 1; break; } } } else if(stp->GetPlaneView()==PlaneView::kV || stp->GetPlaneView()==PlaneView::kY){ for (int i=0; i<=allVStrips.GetLast(); i++) { TObject *tobj = allVStrips.At(i); assert(tobj->InheritsFrom("CandStripHandle")); CandStripHandle *striphandle = dynamic_cast<CandStripHandle*>(tobj); if(*stp==*striphandle) { tmpVArray[i] = 1; break; } } } } //remove used strips from arrays for(int i=0;i<allUStrips.GetLast()+1;i++) { if(tmpUArray[i]==1) allUStrips.RemoveAt(i); } allUStrips.Compress(); if(allUStrips.GetLast()==-1) keepGoingU = false; delete [] tmpUArray; for(int i=0;i<allVStrips.GetLast()+1;i++) { if(tmpVArray[i]==1) allVStrips.RemoveAt(i); } allVStrips.Compress(); if(allVStrips.GetLast()==-1) keepGoingV = false; delete [] tmpVArray; } else { MSG("SubShowerSR",Msg::kDebug) << "SubShowerSR formed with 0 strips - not adding to list" << endl; if(keepGoingU) { MSG("SubShowerSR",Msg::kDebug) << "Running CleanUp() on U Strips" << endl; keepGoingU = CleanUp(&allUStrips,ch,detector,slice); if(!keepGoingU) { for(int i=0;i<allUStrips.GetLast()+1;i++) { unassignedStrips.Add(allUStrips.At(i)); } } } else { MSG("SubShowerSR",Msg::kDebug) << "Running CleanUp() on V Strips" << endl; keepGoingV = CleanUp(&allVStrips,ch,detector,slice); if(!keepGoingV) { for(int i=0;i<allVStrips.GetLast()+1;i++) { unassignedStrips.Add(allVStrips.At(i)); } } } } } sliceCounter++; continue; //Try to merge some clusters: Bool_t stillMergeChanges = true; while(stillMergeChanges){ stillMergeChanges = false; CandSubShowerSRListHandle &csslh = dynamic_cast<CandSubShowerSRListHandle &>(ch); MSG("SubShowerSR",Msg::kDebug) << "In Merge Loop: NDaughters = " << csslh.GetNDaughters() << endl; CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator()); CandSubShowerSRHandleKeyFunc *engKF = subshowerItr.CreateKeyFunc(); engKF->SetFun(CandSubShowerSRHandle::KeyFromReverseViewEnergy); subshowerItr.GetSet()->AdoptSortKeyFunc(engKF); while(CandSubShowerSRHandle *cssh = subshowerItr()){ if(*cssh->GetCandSlice()==*slice){ CandStripHandleItr stripItr(cssh->GetDaughterIterator()); TObjArray mergeSubShower; while(CandStripHandle *stp = stripItr()) mergeSubShower.Add(stp->DupHandle()); if(this->MergeCluster(&mergeSubShower,csslh,detector,true)){ MSG("SubShowerSR",Msg::kDebug) << "Going to Merge: NDaughters = " << csslh.GetNDaughters() << endl; csslh.RemoveDaughter(cssh); MSG("SubShowerSR",Msg::kDebug) << "Removed Initial Daughter: NDaughters = " << csslh.GetNDaughters() << endl; this->MergeCluster(&mergeSubShower,csslh,detector); MSG("SubShowerSR",Msg::kDebug) << "Merged: NDaughters = " << csslh.GetNDaughters() << endl; cxx.SetDataIn(&mergeSubShower); CandSubShowerSRHandle subshowersrhandle = CandSubShowerSR::MakeCandidate(ah,cxx); subshowersrhandle.SetCandSlice(slice); csslh.AddDaughterLink(subshowersrhandle); MSG("SubShowerSR",Msg::kDebug) << "Made New Daughter: NDaughters = " << csslh.GetNDaughters() << endl; stillMergeChanges = true; break; } if(stillMergeChanges) break; } if(stillMergeChanges) break; } } } //slice-wise cleanup if(FormHalo(&unassignedStrips,ch,cxx,ah,slicelist,detector)) MSG("SubShowerSR",Msg::kDebug) << "Halo SubShowers formed" << endl; }

NavKey AlgSubShowerSRList::StripKeyFromPlane (  const CandStripHandle *   )  [static]

Definition at line 3397 of file AlgSubShowerSRList.cxx. { return const_cast<CandStripHandle *>(strip)->GetPlane(); }

Bool_t AlgSubShowerSRList::TestOverLap (  TObjArray *  , CandHandle &  , const CandSliceHandle *  )  [private]

Definition at line 2078 of file AlgSubShowerSRList.cxx. References CandRecoHandle::GetBegPlane(), CandRecoHandle::GetCandSlice(), CandStripHandle::GetCharge(), CandSubShowerSRHandle::GetClusterID(), CandHandle::GetDaughterIterator(), CandRecoHandle::GetEndPlane(), CandHandle::GetNDaughters(), CandStripHandle::GetPlane(), CandSubShowerSRHandle::GetPlaneView(), CandStripHandle::GetPlaneView(), CandSubShowerSRHandle::GetSlope(), CandStripHandle::GetTPos(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandRecoHandle::GetVtxZ(), CandStripHandle::GetZPos(), and MSG. Referenced by RunAlg(). { //Going to check that proto-subshower does not go through an existing subshower //if it does, take largest, plane-wise contiguous set of strips MSG("SubShowerSR", Msg::kDebug) << "Running TestOverLap()" << endl; //if no strips return false. Empty subshower will be caught later. Int_t totStp = newSubShower->GetLast()+1; if(totStp<=1) return false; assert(ch.InheritsFrom("CandSubShowerSRListHandle")); CandSubShowerSRListHandle &csslh = dynamic_cast<CandSubShowerSRListHandle &>(ch); //if this is the first subshower, just return false: if(csslh.GetNDaughters()==0) return false; //Get plane view: TObject *firstObj = newSubShower->At(0); assert(firstObj->InheritsFrom("CandStripHandle")); CandStripHandle *firstStp = dynamic_cast<CandStripHandle*>(firstObj); PlaneView::PlaneView_t pv = firstStp->GetPlaneView(); MSG("SubShowerSR", Msg::kDebug) << "PlaneView = " << pv << endl; //loop through subshowers in view + energy order //so that newSubShower is tested against largest subshower first CandSubShowerSRHandleItr subshowerItr(csslh.GetDaughterIterator()); CandSubShowerSRHandleKeyFunc *engKF = subshowerItr.CreateKeyFunc(); engKF->SetFun(CandSubShowerSRHandle::KeyFromViewEnergy); subshowerItr.GetSet()->AdoptSortKeyFunc(engKF); while (CandSubShowerSRHandle *subshower = subshowerItr()) { if((*subshower->GetCandSlice())!=(*slice)) continue; if(subshower->GetPlaneView()!=pv) continue; if(subshower->GetClusterID()==ClusterType::kXTalk || subshower->GetClusterID()==ClusterType::kHalo) continue; //Get vertex and angle of subshower: Double_t tpos_vtx = 0; Double_t z_vtx = subshower->GetVtxZ(); if((subshower->GetPlaneView()==PlaneView::kU || subshower->GetPlaneView()==PlaneView::kX)) tpos_vtx = subshower->GetVtxU(); if((subshower->GetPlaneView()==PlaneView::kV || subshower->GetPlaneView()==PlaneView::kY)) { tpos_vtx = subshower->GetVtxV(); } Double_t slope = subshower->GetSlope(); Int_t beg_plane = subshower->GetBegPlane(); Int_t end_plane = subshower->GetEndPlane(); //for calculating asymmetry: Double_t n_upp = 0; Double_t n_low = 0; Double_t ph_upp = 0; Double_t ph_low = 0; Double_t n_plane_upp[500] = {0}; Double_t n_plane_low[500] = {0}; Double_t ph_plane_upp[500] = {0}; Double_t ph_plane_low[500] = {0}; Int_t firstPlane = 500; Int_t lastPlane = 0; //calculate number of strips above and below axis of subshower for(int i=0;i<totStp;i++){ CandStripHandle *stp = dynamic_cast<CandStripHandle*>(newSubShower->At(i)); if(firstPlane>stp->GetPlane()) firstPlane = stp->GetPlane(); if(lastPlane<stp->GetPlane()) lastPlane = stp->GetPlane(); if(stp->GetTPos() > (slope*(stp->GetZPos()-z_vtx) + tpos_vtx)) { n_plane_upp[stp->GetPlane()] += 1; ph_plane_upp[stp->GetPlane()] += stp->GetCharge(CalDigitType::kSigCorr); if(stp->GetPlane()>=beg_plane && stp->GetPlane()<=end_plane) { n_upp += 1; ph_upp += stp->GetCharge(CalDigitType::kSigCorr); } } else { n_plane_low[stp->GetPlane()] += 1; ph_plane_low[stp->GetPlane()] += stp->GetCharge(CalDigitType::kSigCorr); if(stp->GetPlane()>=beg_plane && stp->GetPlane()<=end_plane) { n_low += 1; ph_low += stp->GetCharge(CalDigitType::kSigCorr); } } } //if asymmetry < 95% (should be 100%!) in the plane //range of the current subshower then only keep the //larger part of the new subshower: Double_t asymmetry = 1; Double_t asymmetryPH = 1; if( n_upp + n_low > 1 ) asymmetry = TMath::Abs(n_upp-n_low)/(n_upp+n_low); if( ph_upp + ph_low > 1 ) asymmetryPH = TMath::Abs(ph_upp-ph_low)/(ph_upp+ph_low); MSG("SubShowerSR", Msg::kDebug) << "Asymmetry = " << asymmetry << " AsymmetryPH = " << asymmetryPH << " TotStp = " << totStp << " n_upp = " << n_upp << " n_low = " << n_low << " ph_upp = " << ph_upp << " ph_low = " << ph_low << endl; if(asymmetry < 0.95) { //we know now that this proto-subshower is cutting through another subshower //find the greatest number of contiguous planes above or below slope //and take those strips only Int_t bestAsymFirstPlane = firstPlane; Int_t bestAsymLastPlane = firstPlane; Int_t AsymFirstPlane = firstPlane; Int_t currentSense = 0; Int_t bestSense = 0; //loop through all planes and find longest segments for(int i=firstPlane;i<=lastPlane;i+=2){ //find longest set of contiguous planes with same asym sign: if(n_plane_upp[i]>n_plane_low[i] || ( n_plane_upp[i]==n_plane_low[i] && ph_plane_upp[i]>ph_plane_low[i]) ) { if(currentSense==0) { AsymFirstPlane = i; currentSense = 1; } else if(currentSense>0) currentSense += 1; else if(currentSense<0) { //sense has changed: if(TMath::Abs(currentSense)>TMath::Abs(bestSense)){ bestSense = currentSense; bestAsymFirstPlane = AsymFirstPlane; bestAsymLastPlane = i-1; } currentSense = 1; AsymFirstPlane = i; } } else if(n_plane_low[i]>n_plane_upp[i] || ( n_plane_upp[i]==n_plane_low[i] && ph_plane_low[i]>ph_plane_upp[i]) ){ if(currentSense==0) { AsymFirstPlane = i; currentSense = -1; } else if(currentSense<0) currentSense -= 1; else if(currentSense>0) { if(TMath::Abs(currentSense)>TMath::Abs(bestSense)){ bestSense = currentSense; bestAsymFirstPlane = AsymFirstPlane; bestAsymLastPlane = i-1; } currentSense = -1; AsymFirstPlane = i; } } } if(TMath::Abs(currentSense)>TMath::Abs(bestSense)){ bestSense = currentSense; bestAsymFirstPlane = AsymFirstPlane; bestAsymLastPlane = lastPlane; } for(int i=0;i<totStp;i++){ CandStripHandle *stp = dynamic_cast<CandStripHandle*>(newSubShower->At(i)); if(stp->GetPlane()<bestAsymFirstPlane || stp->GetPlane()>bestAsymLastPlane){ newSubShower->RemoveAt(i); } } newSubShower->Compress(); if(newSubShower->GetLast()+1<totStp) { return true; } } } return false; }

void AlgSubShowerSRList::Trace (  const char *  c  )  const [virtual]

Reimplemented from AlgBase. Definition at line 3403 of file AlgSubShowerSRList.cxx. { }

std::vector< window > AlgSubShowerSRList::TransCluster (  std::vector< window >  , Int_t  )  [private]

Definition at line 1053 of file AlgSubShowerSRList.cxx. References BestHough(), window::centroid, HoughTransCluster(), window::lower, window::lower_tpos, MSG, window::ph, window::phpos, window::phwidth, window::plane, window::type, window::upper, window::upper_tpos, and win. Referenced by FindCluster(). { std::vector<window> cluster; Int_t maxPlane = -1; //calculate maxPlane Double_t maxPlanePH = 0; Double_t PlanePH[500] = {0}; std::vector<window>::iterator winIter = win.begin(); std::vector<window>::iterator winEnd = win.end(); while(winIter!=winEnd){ PlanePH[winIter->plane] += winIter->ph; if(PlanePH[winIter->plane]>maxPlanePH) { maxPlane = winIter->plane; maxPlanePH = PlanePH[maxPlane]; } winIter++; } Int_t maxLongPlane = -1; Int_t minLongPlane = -1; Int_t counter = 0; Int_t delta = 2; while(PlanePH[maxPlane+counter]>0||(maxPlane+counter)==249) { maxLongPlane = maxPlane+counter; delta = 2; if(det==1&&(maxLongPlane-249)*(maxLongPlane+delta-249)<0) delta = 3; counter+=delta; } counter = 0; while(PlanePH[maxPlane+counter]<=0||(maxPlane+counter)==249) { minLongPlane = maxPlane+counter; delta = 2; if(det==1&&(minLongPlane-249)*(minLongPlane-delta-249)<0) delta = 3; counter-=delta; } Int_t numPreMaxPlanes = maxPlane-minLongPlane+1; Int_t numPostMaxPlanes = maxLongPlane-maxPlane+1; Int_t numLongPlanes = numPostMaxPlanes + numPreMaxPlanes - 1; //match up windows plane by plane Double_t shwmid = 0; // centroid of cluster transversely Double_t shwwid = 0; // trans width of widest window in proto-cluster Double_t maxwid = 0; Double_t clusterph = 0; // total PH of cluster Double_t refp = 0; // reference plane } essentially pivot point for Double_t refs = 0; // reference strip } cluster //long. window info: Int_t *winpl = new Int_t[numLongPlanes]; // plane number Int_t *winpl0 = new Int_t[numLongPlanes]; // upper bound Int_t *winpl1 = new Int_t[numLongPlanes]; // lower bound Float_t *winpl0_tpos = new Float_t[numLongPlanes]; // upper bound in tpos Float_t *winpl1_tpos = new Float_t[numLongPlanes]; // lower bound in tpos Double_t *winpl2 = new Double_t[numLongPlanes]; // PH Int_t *winpl3 = new Int_t[numLongPlanes]; // window type Double_t *winpl4 = new Double_t[numLongPlanes]; // width Double_t *winpl5 = new Double_t[numLongPlanes]; // centroid Double_t slopen = 0; Double_t eslopen = 0; Bool_t trustslopen = false; std::vector<window> houghcluster = HoughTransCluster(win,det); if(houghcluster.size()==0) { //houghcluster.size()=0 try usual tracking for (int d = 0;d<numLongPlanes;d++){ int pl = -1; if(numPreMaxPlanes>numPostMaxPlanes){ if (d<numPreMaxPlanes) pl = maxPlane-d; else pl = maxPlane+d-numPreMaxPlanes+1; } else{ if (d<numPostMaxPlanes) pl = maxPlane+d; else pl = maxPlane-d+numPostMaxPlanes-1; } winpl[d] = 0; winpl0[d] = -1; winpl1[d] = -1; winpl0_tpos[d] = -10; winpl1_tpos[d] = -10; winpl2[d] = 0; winpl3[d] = 0; winpl4[d] = -1; winpl5[d] = -1; // PH for current best matched window: double weight = 0.; //deviation of best matched window from prediction of //centroid for this plane: double weight1 = 10000.; //loop keeps track of # windows already checked on this plane: int loop = 0; // counts number of windows with a good match: int within = 0; //prediction of centroid for this plane (not using slope in info): double mid0 = 0; //prediction of centroid for this plane using slope info where available: double mid = 0; //loop over all in windows winIter = win.begin(); while (winIter!=winEnd){ if (winIter->plane==pl){ //if it's in current plane: loop++; double dev = shwwid; if (dev<(winIter->upper_tpos-winIter->lower_tpos+STRIPWIDTHINMETRES)) dev = winIter->upper_tpos-winIter->lower_tpos+STRIPWIDTHINMETRES; if (dev<2.*STRIPWIDTHINMETRES) dev = 2.*STRIPWIDTHINMETRES; if (shwmid!=0) { mid0 = shwmid; mid = shwmid; } else { mid = (winIter->upper+winIter->lower)/2.; mid0 = mid; } // check if this window is consistent with previously selected // windows and has a higher PH than current best match for this // plane or else if it's shower max plane if((shwmid==0 || (TMath::Abs(winIter->upper_tpos-mid)<dev && TMath::Abs(winIter->lower_tpos-mid)<dev)) && winIter->ph>weight){ //loop through planes in long. window again to do a local //scan to ensure best match is found. //step forward and backward once each if possible int loopp = 0; int loopm = 0; int withinp = 0; int withinm = 0; //double midp0 = 0; //double midm0 = 0; double midp = 0; double midm = 0; // step forward: std::vector<window>::iterator winIter1 = win.begin(); while(winIter1!=winEnd){ if ((((det==1&&(pl-249)*(pl+2-249)>0)||det==0) &&winIter1->plane==pl+2) ||(det==1&&(pl-249)*(pl+2-249)<0 &&winIter1->plane==pl+3)){ //go to next same view plane loopp++; if(shwmid!=0){ //not shower max plane midp = shwmid; if(TMath::Abs(winIter1->upper_tpos-midp)<dev/2. ||TMath::Abs(winIter1->lower_tpos-midp)<dev/2.){ //good match withinp++; } } else if(shwmid==0){ //in shower max plane midp = (winIter->upper_tpos+winIter->lower_tpos)/2.; if(TMath::Abs(winIter1->upper_tpos-midp)<dev/2. ||TMath::Abs(winIter1->lower_tpos-midp)<dev/2.){ //good match withinp++; } } } //step backwards: if ((((det==1&&(pl-249)*(pl-2-249)>0)||det==0) &&winIter1->plane==pl-2) ||(det==1&&(pl-249)*(pl-2-249)<0 &&winIter1->plane==pl-3)){ //go to previous same view plane loopm++; if(shwmid!=0){ // not shower max plane midm = shwmid; if(TMath::Abs(winIter1->upper_tpos-midm)<dev ||TMath::Abs(winIter1->lower_tpos-midm)<dev){ //good match withinm++; } } else if(shwmid==0){ //shower max plane midm = (winIter->upper_tpos+winIter->lower_tpos)/2.; if(TMath::Abs(winIter1->upper_tpos-midm)<dev ||TMath::Abs(winIter1->lower_tpos-midm)<dev){ //good match withinm++; } } } winIter1++; } if((withinp>0 && withinm>0) || (loopp==0 && withinm>0) || (loopm==0 && withinp>0)){ //success! current window is a good solution within++; // write out solution: winpl[d] = winIter->plane; winpl0[d] = winIter->upper; winpl1[d] = winIter->lower; winpl0_tpos[d] = winIter->upper_tpos; winpl1_tpos[d] = winIter->lower_tpos; winpl2[d] = winIter->ph; winpl3[d] = winIter->type; weight = winIter->ph; weight1 = (winIter->upper_tpos+winIter->lower_tpos)/2.-mid0; //these will be overwritten if there is another window //which is also a good match, but has more PH, //ensuring we get best possible match } } } winIter++; } if(loop*within>0){ //got a matched window if (shwmid==0) { // shower max plane refp = winpl[d]; // set reference point refs = (winpl0_tpos[d]+winpl1_tpos[d])/2.; } //determine widest window: if(shwwid<(winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES)) shwwid = winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES; // average centroid of event so far: shwmid = (shwmid*clusterph+(winpl0_tpos[d]+winpl1_tpos[d]) * winpl2[d]/2.)/(clusterph+winpl2[d]); clusterph+=winpl2[d]; winpl4[d] = shwwid; winpl5[d] = shwmid; } } //Fitting the boundaries of windows matched Double_t slopetn = 0; Double_t slopebn = 0; Double_t eslopetn = 0; Double_t eslopebn = 0; Double_t chi2tn = 0; Double_t chi2bn = 0; int nt = 0; int nb = 0; for (int d = 0;d<numLongPlanes;d++){ int pl = -1; if(numPreMaxPlanes>numPostMaxPlanes){ if (d<numPreMaxPlanes) pl = maxPlane-d; else pl = maxPlane+d-numPreMaxPlanes+1; } else{ if (d<numPostMaxPlanes) pl = maxPlane+d; else pl = maxPlane-d+numPostMaxPlanes-1; } if(winpl2[d]>0) {nt++; nb++;} if (maxwid<winpl4[d]) maxwid = winpl4[d]; if (maxwid<2.*STRIPWIDTHINMETRES) maxwid = 2.*STRIPWIDTHINMETRES; } if (nt>2&&nb>2) { TGraphErrors *grat = new TGraphErrors(nt); TGraphErrors *grab = new TGraphErrors(nb); int t = 0; int b = 0; for (int d = 0;d<numLongPlanes;d++){ int pl = -1; if(numPreMaxPlanes>numPostMaxPlanes){ if (d<numPreMaxPlanes) pl = maxPlane-d; else pl = maxPlane+d-numPreMaxPlanes+1; } else{ if (d<numPostMaxPlanes) pl = maxPlane+d; else pl = maxPlane-d+numPostMaxPlanes-1; } if(winpl2[d]>0) { grat->SetPoint(t,winpl[d],winpl0_tpos[d]); grat->SetPointError(t,0,1./TMath::Sqrt(winpl2[d])); t++; grab->SetPoint(b,winpl[d],winpl1_tpos[d]); grab->SetPointError(b,0,1./TMath::Sqrt(winpl2[d])); b++; } } TF1 *p1 = new TF1("p1","pol1",0,500); grat->Fit("p1","NNLQ"); slopetn = p1->GetParameter(1); eslopetn = p1->GetParError(1); chi2tn = p1->GetChisquare(); grab->Fit("p1","NNLQ"); slopebn = p1->GetParameter(1); eslopebn = p1->GetParError(1); chi2bn = p1->GetChisquare(); delete grat; delete grab; delete p1; if ((slopetn+slopebn)!=0. && slopetn*slopebn!=0.){ if(chi2tn/nt<100. && chi2bn/nb<100. && ((((eslopetn/TMath::Abs(slopetn)<0.1 && eslopebn/TMath::Abs(slopebn)<0.1 &&chi2tn/nt>0.1 && chi2bn/nb>0.1)) && (TMath::Abs(slopetn-slopebn)/TMath::Abs(slopetn+slopebn)<0.25)) || TMath::Abs(slopetn-slopebn)==0)) { slopen = (slopetn+slopebn)/2.; eslopen = TMath::Sqrt(TMath::Power(eslopetn,2)+TMath::Power(eslopebn,2))/2.; if(TMath::Abs(slopetn-slopebn)>eslopen) eslopen = TMath::Abs(slopetn-slopebn); trustslopen = true; } } MSG("SubShowerSR",Msg::kDebug) <<"fit window "<<slopetn<<" "<<slopebn <<" "<<slopen<<" "<<eslopetn<<" " <<eslopebn<<" "<<eslopen<<" " <<chi2tn/nt<<" "<<chi2bn/nb<<" " <<nt<<" "<<nb<<endl; } //check if solution is good: double usewid = maxwid; for (Int_t d = 0;d<numLongPlanes;d++){ int pl = -1; if(numPreMaxPlanes>numPostMaxPlanes){ if (d<numPreMaxPlanes) pl = maxPlane-d; else pl = maxPlane+d-numPreMaxPlanes+1; } else{ if (d<numPostMaxPlanes) pl = maxPlane+d; else pl = maxPlane-d+numPostMaxPlanes-1; } if(d>0 && (winpl0[d]>=0||winpl1[d]>=0) ){ //if there is a solution if(!trustslopen && //don't trust new slope //window is too far away from centroid (winpl0_tpos[d]-shwmid)*(winpl1_tpos[d]-shwmid) > (maxwid/2.)*(maxwid/2.) ){ MSG("SubShowerSR",Msg::kDebug) << "removed plane " << winpl[d] << " " << winpl0[d] << " " << winpl1[d] << " " << winpl0_tpos[d] << " " << winpl1_tpos[d] << endl; shwmid = (shwmid*clusterph-(winpl0_tpos[d]+winpl1_tpos[d]) * winpl2[d]/2.)/(clusterph-winpl2[d]); winpl0[d] = -1; winpl1[d] = -1; winpl0_tpos[d] = -10; winpl1_tpos[d] = -10; winpl2[d] = 0; winpl3[d] = 0; winpl4[d] = -1; winpl5[d] = -1; } else if(trustslopen){ //this time account for slope on window boundaries: //try to find another window as a reference point //to compare the window in this plane with int i = 0; while(winpl2[i]==0&&i<=d) i++; if(i==d){ while(winpl2[i]==0&&i<numLongPlanes) i++; } if(usewid<TMath::Abs(winpl[d]-winpl[i])*eslopen) //find appropriate cut value on window width usewid = TMath::Abs(winpl[d]-winpl[i])*eslopen; //expected centroid for plane: double expmid = (winpl0_tpos[i] + winpl1_tpos[i])/2.+(winpl[d]-winpl[i])*slopen; if((winpl0_tpos[d]-expmid)*(winpl1_tpos[d]-expmid) > (usewid/2.)*(usewid/2.)){ MSG("SubShowerSR",Msg::kDebug) <<"removed plane using slope " <<winpl[d]<<" "<<winpl0[d]<<" " <<winpl1[d]<<" with pl " <<maxPlane<<" "<<winpl[i] <<" "<<expmid<<" "<<usewid<<endl; shwmid = (shwmid*clusterph-(winpl0_tpos[d]+winpl1_tpos[d]) * winpl2[d]/2.)/(clusterph-winpl2[d]); clusterph-=winpl2[d]; winpl0[d] = -1; winpl1[d] = -1; winpl0_tpos[d] = -10; winpl1_tpos[d] = -10; winpl2[d] = 0; winpl3[d] = 0; winpl4[d] = -1; winpl5[d] = -1; } } } } } else { //use hough solution as a starting point: for (int d = 0;d<numLongPlanes;d++){ int pl = -1; if(numPreMaxPlanes>numPostMaxPlanes){ if (d<numPreMaxPlanes) pl = maxPlane-d; else pl = maxPlane+d-numPreMaxPlanes+1; } else{ if (d<numPostMaxPlanes) pl = maxPlane+d; else pl = maxPlane-d+numPostMaxPlanes-1; } winpl[d] = 0; winpl0[d] = -1; winpl1[d] = -1; winpl0_tpos[d] = -10; winpl1_tpos[d] = -10; winpl2[d] = 0; winpl3[d] = 0; winpl4[d] = -1; winpl5[d] = -1; //loop over all in windows std::vector<window>::iterator hcIter = houghcluster.begin(); std::vector<window>::iterator hcEnd = houghcluster.end(); while (hcIter!=hcEnd){ if (hcIter->plane==pl){ //if it's in current plane: winpl[d] = hcIter->plane; winpl0[d] = hcIter->upper; winpl1[d] = hcIter->lower; winpl0_tpos[d] = hcIter->upper_tpos; winpl1_tpos[d] = hcIter->lower_tpos; winpl2[d] = hcIter->ph; winpl3[d] = hcIter->type; if (shwmid==0) { // shower max plane refp = winpl[d]; // set reference point refs = (winpl0_tpos[d]+winpl1_tpos[d])/2.; } //determine widest window: if(shwwid<(winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES)) shwwid = winpl0_tpos[d]-winpl1_tpos[d]+STRIPWIDTHINMETRES; // average centroid of event so far: shwmid = (shwmid*clusterph+(winpl0_tpos[d]+winpl1_tpos[d]) * winpl2[d]/2.)/(clusterph+winpl2[d]); // calculate average distance to centroid: clusterph += winpl2[d]; winpl4[d] = shwwid; winpl5[d] = shwmid; if (maxwid<winpl4[d]) maxwid = winpl4[d]; if (maxwid<2.*STRIPWIDTHINMETRES) maxwid = 2.*STRIPWIDTHINMETRES; } hcIter++; } } Double_t *MCV = new Double_t[14]; MCV[0] = 0; MCV[1] = 0; MCV[2] = 0; MCV[3] = 0; MCV[4] = 0; MCV[5] = 0; MCV[6] = 0; MCV[7] = 0; MCV[8] = 0; MCV[9] = 0; MCV[10] = 0; MCV[11] = 0; MCV[12] = 0; MCV[13] = 0; BestHough(houghcluster,MCV,false); if(MCV[5]!=0){ slopen = MCV[5]; eslopen = MCV[8]; trustslopen = true; } else if(MCV[0]!=0){ slopen = MCV[0]; eslopen = MCV[3]; trustslopen = true; } else { slopen = 0; eslopen = 0; trustslopen = false; } delete [] MCV; } //Try to fill planes with no windows //First find centroid and width for +/- 1 planes for (Int_t d = 0;d<numLongPlanes;d++){ int pl = -1; if(numPreMaxPlanes>numPostMaxPlanes){ if (d<numPreMaxPlanes) pl = maxPlane-d; else pl = maxPlane+d-numPreMaxPlanes+1; } else{ if (d<numPostMaxPlanes) pl = maxPlane+d; else pl = maxPlane-d+numPostMaxPlanes-1; } double minwei = 200.; //minimum distance to centroid (weight) double minwei_tpos = 10.; //minimum distance to centroid (weight) int newwinpl0 = -1; //new window upper and lower strip bounds int newwinpl1 = -1; Float_t newwinpl0_tpos = -10; //new window upp and low strip tpos bounds Float_t newwinpl1_tpos = -10; double newwinpl2 = 0.; int newwinpl3 = 0; if(winpl0[d]==-1&&winpl1[d]==-1){ //does this plane have a solution? int p = 0; int m = 0; if (det==0||(det==1&&(pl-249)*(pl+2-249)>0)) p = 2; else p = 3; if (det==0||(det==1&&(pl-249)*(pl-2-249)>0)) m = 2; else m = 3; double pmid = -1; //centroid for next plane (plus) double mmid = -1; //centroid for prev plane (minus) double pwid = -1; //width double mwid = -1; int pwt = -1; // window type plus int mwt = -1; // minus bool useslopep = false; bool useslopem = false; for (Int_t c = 0;c<numLongPlanes;c++){ if (winpl[c]==pl+p) { pmid = winpl5[c]; if (trustslopen) useslopep = true; pwid = winpl4[c]; pwt = winpl3[c]; } if (winpl[c]==pl-m) { mmid = winpl5[c]; if (trustslopen) useslopem = true; mwid = winpl4[c]; mwt = winpl3[c]; } } //Find compatible windows double usewid = maxwid; if(maxwid>0){ if((mmid==-1&&pmid==-1)|| (!useslopem&&mmid>=0&&TMath::Abs(shwmid-mmid)>maxwid)|| (!useslopep&&pmid>=0&&TMath::Abs(shwmid-pmid)>maxwid)|| (useslopem&&mmid>=0&& (TMath::Abs(shwmid-mmid-(maxPlane-pl+m)*slopen)>maxwid)|| (TMath::Abs(shwmid-mmid-(maxPlane-pl+m)*slopen) > TMath::Abs(maxPlane-pl+m)*eslopen))|| (useslopep&&pmid>=0&& (TMath::Abs(shwmid-pmid-(maxPlane-pl-p)*slopen)>maxwid)|| (TMath::Abs(shwmid-pmid-(maxPlane-pl-p)*slopen) > TMath::Abs(maxPlane-pl-p)*eslopen))){ minwei = 200.; minwei_tpos = 10.; newwinpl0 = -1; newwinpl1 = -1; newwinpl0_tpos = -10; newwinpl1_tpos = -10; newwinpl2 = 0.; winIter = win.begin(); while(winIter!=winEnd){ if(winIter->plane==pl){ if (//shwmid>=0&& ((winIter->upper_tpos-shwmid) * (winIter->lower_tpos-shwmid)<= (maxwid/2.)*(maxwid/2.)) && TMath::Abs((winIter->upper_tpos + winIter->lower_tpos)/2.-shwmid)<minwei_tpos){ minwei = TMath::Abs((winIter->upper_tpos + winIter->lower_tpos)/2.-shwmid); newwinpl0 = winIter->upper; newwinpl1 = winIter->lower; newwinpl0_tpos = winIter->upper_tpos; newwinpl1_tpos = winIter->lower_tpos; newwinpl2 = winIter->ph; newwinpl3 = winIter->type; } } winIter++; } } else{ minwei = 200.; minwei_tpos = 10.; newwinpl0 = -1; newwinpl1 = -1; newwinpl0_tpos = -10.; newwinpl1_tpos = -10.; newwinpl2 = 0.; winIter = win.begin(); while(winIter!=winEnd){ if(winIter->plane==pl){ if (!useslopem&&mmid>=0 &&((winIter->upper_tpos-mmid) * (winIter->lower_tpos-mmid) <= (maxwid/2.)*(maxwid/2.)) && TMath::Abs((winIter->upper_tpos + winIter->lower_tpos)/2.-mmid)<minwei_tpos){ minwei = TMath::Abs((winIter->upper_tpos + winIter->lower_tpos)/2.-mmid); newwinpl0 = winIter->upper; newwinpl1 = winIter->lower; newwinpl0_tpos = winIter->upper_tpos; newwinpl1_tpos = winIter->lower_tpos; newwinpl2 = winIter->ph; newwinpl3 = winIter->type; } else if (!useslopep&&pmid>=0 &&((winIter->upper_tpos-pmid) * (winIter->lower_tpos-pmid) <= (maxwid/2.)*(maxwid/2.)) && TMath::Abs((winIter->upper_tpos + winIter->lower_tpos)/2.-pmid)<minwei){ minwei = TMath::Abs((winIter->upper_tpos + winIter->lower_tpos)/2.-pmid); newwinpl0 = winIter->upper; newwinpl1 = winIter->lower; newwinpl0_tpos = winIter->upper_tpos; newwinpl1_tpos = winIter->lower_tpos; newwinpl2 = winIter->ph; newwinpl3 = winIter->type; } else if (useslopem&&mmid>=0){ if(usewid<m*eslopen) usewid = m*eslopen; else usewid = maxwid; if((winIter->upper_tpos-mmid-m*slopen) * (winIter->lower_tpos-mmid-m*slopen) <= (usewid/2.)*(usewid/2.) &&TMath::Abs((winIter->upper_tpos + winIter->lower_tpos)/2.-mmid-m*slopen)<minwei){ minwei = TMath::Abs((winIter->upper_tpos + winIter->lower_tpos)/2.-mmid-m*slopen); newwinpl0 = winIter->upper; newwinpl1 = winIter->lower; newwinpl0_tpos = winIter->upper_tpos; newwinpl1_tpos = winIter->lower_tpos; newwinpl2 = winIter->ph; newwinpl3 = winIter->type; } } else if (useslopep&&pmid>=0){ if(usewid<p*eslopen) usewid = p*eslopen; else usewid = maxwid; if((winIter->upper_tpos-pmid+p*slopen) * (winIter->lower_tpos-pmid+p*slopen) <= (usewid/2.)*(usewid/2.) &&TMath::Abs((winIter->upper_tpos + winIter->lower_tpos)/2.-pmid+p*slopen)<minwei){ minwei = TMath::Abs((winIter->upper_tpos + winIter->lower_tpos)/2.-pmid+p*slopen); newwinpl0 = winIter->upper; newwinpl1 = winIter->lower; newwinpl0_tpos = winIter->upper_tpos; newwinpl1_tpos = winIter->lower_tpos; newwinpl2 = winIter->ph; newwinpl3 = winIter->type; } } } winIter++; } } } if(newwinpl0>=0) { winpl0[d] = newwinpl0; winpl1[d] = newwinpl1; winpl0_tpos[d] = newwinpl0_tpos; winpl1_tpos[d] = newwinpl1_tpos; winpl2[d] = newwinpl2; winpl3[d] = newwinpl3; shwmid = (shwmid*clusterph + (winpl0_tpos[d]+winpl1_tpos[d]) * winpl2[d]/2.)/(clusterph+winpl2[d]); clusterph+=winpl2[d]; MSG("SubShowerSR",Msg::kDebug) <<"find missed windows "<<pl <<" "<<newwinpl0<<" " <<newwinpl1<<" "<<minwei<<endl; } } window winSolution; winSolution.plane = pl; winSolution.upper = winpl0[d]; winSolution.lower = winpl1[d]; winSolution.upper_tpos = winpl0_tpos[d]; winSolution.lower_tpos = winpl1_tpos[d]; winSolution.ph = winpl2[d]; winSolution.type = winpl3[d]; winSolution.centroid = 0.; winSolution.phpos = 0.; winSolution.phwidth = 0.; cluster.push_back(winSolution); } delete [] winpl; delete [] winpl0; delete [] winpl1; delete [] winpl0_tpos; delete [] winpl1_tpos; delete [] winpl2; delete [] winpl3; delete [] winpl4; delete [] winpl5; return cluster; }

---

Member Data Documentation

Double_t AlgSubShowerSRList::fLongWindowMipCut [private]

Definition at line 59 of file AlgSubShowerSRList.h. Referenced by RunAlg().

Int_t AlgSubShowerSRList::fMaxPlaneGap [private]

Definition at line 61 of file AlgSubShowerSRList.h. Referenced by RunAlg().

Double_t AlgSubShowerSRList::fMergeBestFracBoundary [private]

Definition at line 56 of file AlgSubShowerSRList.h. Referenced by RunAlg().

Double_t AlgSubShowerSRList::fMergeBestFracCross [private]

Definition at line 58 of file AlgSubShowerSRList.h. Referenced by RunAlg().

Double_t AlgSubShowerSRList::fMergeBestFracWithin [private]

Definition at line 57 of file AlgSubShowerSRList.h. Referenced by RunAlg().

Double_t AlgSubShowerSRList::fMergeFracEnergyLimit [private]

Definition at line 55 of file AlgSubShowerSRList.h. Referenced by RunAlg().

---

The documentation for this class was generated from the following files:

• AlgSubShowerSRList.h
• AlgSubShowerSRList.cxx

---