clusterss::MakeClusterSS Class Reference

#include <MakeClusterSS.h>

Inheritance diagram for clusterss::MakeClusterSS:

List of all members.

Public Member Functions
	MakeClusterSS (const char *version)
void	Update (const cfg::Config &c)
	~MakeClusterSS ()
jobc::Result	Reco (edm::EventHandle &evt)
jobc::Result	Ana (const edm::EventHandle &evt)
void	SplitIntoViews (std::vector< const recobase::CellHit * > &hitlist)
bool	DoCluster (geo::View_t view)
void	MakePlaneClusters (geo::View_t view)
void	FindLongWindow (geo::View_t view)
void	FindAllTransverseWindows (geo::View_t view)
bool	ClusterWindows (geo::View_t view)
void	CleanUpWindows ()
std::vector< clusterss::WindowSS * >	HoughTransCluster ()
bool	StupidHoughNoWork (geo::View_t view, ClustVars &cv)
void	CheckCurrSoln (ClustVars &cv)
bool	PruneCellList (geo::View_t view)
void	MergeClusters (geo::View_t view)
houghResult	BestHough (std::vector< clusterss::WindowSS * > hvec)
void	CleanUpContainers ()
void	CleanUpAllContainers ()
int	FindPlaneFromShwMax (geo::View_t view, int away)
Static Public Attributes
const double	ADCTOMIP = 1.
Private Attributes
int	printlevel
geo::Geometry *	geom
int	evtcntr
int	run
int	event
int	begPlane
int	endPlane
int	begWinPlane
int	endWinPlane
int	shwMaxPlane
int	numPreMaxPlanes
int	numPostMaxPlanes
Bool_t	pMergeClusters
Double_t	pMaxOvlpLink
Int_t	pMaxPDiff
Double_t	pMaxZDiff
Double_t	pMaxTDiff
Double_t	pMaxAngleDiff
Float_t	pMaxWidToLenP2
Double_t	pLongWindowMipCut
Int_t	pLongWindowGapPlaneCut
Int_t	pMinHoughPlanes
Double_t	pMinHoughPH
Double_t	pBestHoughMaxPHFrac
int	pMAXHOUGHITER
std::vector< clusterss::CellHitGrad * >	xchglist
std::vector< clusterss::CellHitGrad * >	ychglist
std::map< int, recobase::Cluster * >	plnclusters
std::vector< clusterss::WindowSS * >	plnwindows
std::vector< clusterss::WindowSS * >	windows
std::vector< recobase::Cluster * >	xclusterlist
std::vector< recobase::Cluster * >	yclusterlist
TH2D *	houghHist
TH2D *	phHoughHist

---

Constructor & Destructor Documentation

MakeClusterSS::MakeClusterSS (  const char *  version  )

Definition at line 51 of file MakeClusterSS.cxx. References houghHist, phHoughHist, and jobc::Module::SetCfgVersion(). : jobc::Module("MakeClusterSS"), printlevel(0), geom(), evtcntr(0), run(-1), event(-1), begPlane(-1), endPlane(-1), begWinPlane(-1), endWinPlane(-1), shwMaxPlane(-1), numPreMaxPlanes(-1), numPostMaxPlanes(-1), pMergeClusters(true), pMaxOvlpLink(0), pMaxPDiff(3), pMaxZDiff(18.), pMaxTDiff(12.), pMaxAngleDiff(1.), pMaxWidToLenP2(5.), pLongWindowMipCut(0.), pLongWindowGapPlaneCut(3), pMinHoughPlanes(5), pMinHoughPH(1.0), pBestHoughMaxPHFrac(0.75), pMAXHOUGHITER(5), xclusterlist(), yclusterlist() { houghHist = new TH2D("houghHist", "Hough Histogram",1,-1,1,1,-1,1); phHoughHist = new TH2D("phHoughHist", "PH Weighted Hough Histogram",1,-1,1,1,-1,1); // Set the configuration version tag this->SetCfgVersion(version); // Required! }

MakeClusterSS::~MakeClusterSS (   )

Definition at line 111 of file MakeClusterSS.cxx. References CleanUpAllContainers(). { CleanUpAllContainers(); }

---

Member Function Documentation

jobc::Result MakeClusterSS::Ana (  const edm::EventHandle &  evt  )  [virtual]

Reimplemented from jobc::Module. Definition at line 208 of file MakeClusterSS.cxx. References edm::EventHandle::Reco(), xclusterlist, and yclusterlist. { std::vector<const recobase::Cluster *>xclusterlist(0); bool foundx=true; try{evt.Reco().Get("./SSClusterX",xclusterlist);} catch(edm::Exception e){ cerr<<"Couldn't find an xcluster"<<endl; foundx=false; } std::vector<const recobase::Cluster *>yclusterlist(0); bool foundy=true; try{evt.Reco().Get("./SSClusterY",yclusterlist);} catch(edm::Exception e){ cerr<<"Couldn't find an ycluster"<<endl; foundy=false; } std::vector<const recobase::Cluster *> &clist=xclusterlist; string sv[2]={"X","Y"}; for(int i=0;i<2;i++){ bool fc=false; geo::View_t view; if(i==0){ clist=xclusterlist; fc = foundx; view=geo::kX; } else{ clist=yclusterlist; fc = foundy; view=geo::kY; } cout<<"******************************MakeClusterPrinting "<<sv[i]<<" clusters. Found "<<clist.size()<<endl; if(fc){ for(unsigned int i=0;i<clist.size();i++){ int ncell = clist[i]->NCell(view); cout<<""<<endl; cout<<"**********Printing cluster #"<<i<<" with NCELLS="<<ncell<<endl; for(int j=0;j<ncell;j++){ cout<<"Cell: "<<j<<endl; clist[i]->Cell(view,j)->Print(""); } } } } return jobc::kPassed; // kFailed if you want to fail the event }

houghResult MakeClusterSS::BestHough (  std::vector< clusterss::WindowSS * >  hvec  )

Definition at line 2173 of file MakeClusterSS.cxx. References geo::PlaneGeo::Cell(), clusterss::houghResult::flat, geom, geo::CellGeo::HalfW(), houghHist, clusterss::houghStats::inter, clusterss::houghStats::iZ, clusterss::stats::mean, clusterss::stats::peak, clusterss::houghResult::ph, phHoughHist, geo::Geometry::Plane(), clusterss::stats::rms, and clusterss::houghStats::slope. Referenced by ClusterWindows(), and HoughTransCluster(). { //TV: this code can be streamlined quite a bit. //There are a few loops over things we've already figured out //Come back to this and make it cleaner soon std::vector<WindowSS*>::iterator winIterBeg = winVec.begin(); std::vector<WindowSS*>::iterator winIterEnd = winVec.end(); houghResult MCV; Double_t STRIPWIDTHINMETRES=0; //get extrema of windows: Int_t MaxPlane = -1; Int_t MinPlane = 1985; Double_t MaxStrip = -8.; Double_t MinStrip = 8.; //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[2000] = {0}; //loop over vector of windoes while(winIterBeg!=winIterEnd){ //keep track of max/mins // cout<<"on plane "<<(*winIterBeg)->plane<<" upper "<<(*winIterBeg)->upper<<" lower "<<(*winIterBeg)->lower<<endl; // cout<<" upper tpos "<<(*winIterBeg)->upper_tpos<<" lower tpos "<<(*winIterBeg)->lower_tpos<<" ph "<<(*winIterBeg)->ph<<endl; 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; //keep track of total totPH+=(*winIterBeg)->ph; //keep track of plane PH PlanePH[(*winIterBeg)->plane]+= (*winIterBeg)->ph; if(PlanePH[(*winIterBeg)->plane]>maxPlanePosPH) { maxPlanePos = (*winIterBeg)->plane;//this is plane of shw max. maxPlanePosPH = PlanePH[maxPlanePos]; } // const geo::PlaneGeo gplane = geom->Plane((*winIterBeg)->plane); const geo::PlaneGeo gplane = geom->Plane(10); const geo::CellGeo gstrip = gplane.Cell(10); STRIPWIDTHINMETRES=gstrip.HalfW()*2.; winIterBeg++; } //window tpos position with max PH, and value of max PH Double_t maxWindowPos = 0.; Double_t maxWindowPH = 0.; winIterBeg = winVec.begin(); //loop over windows again //to find transverse position in plane of shower max while(winIterBeg!=winIterEnd){ if((*winIterBeg)->plane==maxPlanePos){ if((*winIterBeg)->ph>maxWindowPH){ maxWindowPos = (*winIterBeg)->phpos; maxWindowPH = (*winIterBeg)->ph; } } winIterBeg++; } // cout<<"Best hough: max plane "<<MaxPlane<<" min "<<MinPlane<<" max strip "<<MaxStrip<<" min "<<MinStrip<<endl; // cout<<" tot pH "<<totPH<<" maxWinPH "<<maxWinPH<<" maxplane pos "<<maxPlanePos<<" maxPlanePosPH "<<maxPlanePosPH<<endl; //if we didn't have any windows, its just a one planer if(winVec.size()<=1 || MinPlane==MaxPlane){ MCV.flat.slope.mean = 0; MCV.flat.inter.mean = maxWindowPos; MCV.flat.iZ = maxPlanePos; MCV.flat.slope.rms = 0; MCV.flat.inter.rms = 0; MCV.ph.slope.mean = 0; MCV.ph.inter.mean = maxWindowPos; MCV.ph.iZ = maxPlanePos; MCV.ph.slope.rms = 0; MCV.ph.inter.rms = 0; MCV.flat.slope.peak = 0; MCV.flat.inter.peak = maxWindowPos; MCV.ph.slope.peak = 0; MCV.ph.inter.peak = maxWindowPos; return MCV; } Double_t TPosWin = MaxStrip-MinStrip+STRIPWIDTHINMETRES; //MSG("SubShowerSR",Msg::kDebug) << "Hough extrema: plane " << MinPlane // 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 houghHist->Reset(); Int_t param1 = Int_t(TPosWin/STRIPWIDTHINMETRES) + MaxPlane-MinPlane+2 ; Double_t param2 = -2.*TPosWin/Float_t(MaxPlane-MinPlane+2); Double_t param3 = 2.*TPosWin/Float_t(MaxPlane-MinPlane+2); Int_t param4 = Int_t(2.*TPosWin/STRIPWIDTHINMETRES); Double_t param5 = MinStrip-TPosWin ; Double_t param6 = MaxStrip+TPosWin; houghHist->SetBins(param1,param2,param3,param4,param5,param6); phHoughHist->Reset(); phHoughHist->SetBins(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; // loop over all possible intercepts for(Int_t ic=1; ic<=houghHist->GetYaxis()->GetNbins(); ic++) { Double_t c = houghHist->GetYaxis()->GetBinCenter(ic); Double_t m = (strip - c)/Float_t(plane-MinPlane+2); //MSG("SubShowerSR",Msg::kVerbose) << "Hough entries: " << m << " " // << c << " " << ph << endl; 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*pBestHoughMaxPHFrac){ 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*pBestHoughMaxPHFrac) { 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.flat.iZ = MinPlane - 2; if(totx>0) { MCV.flat.slope.mean = sumx/totx; MCV.flat.slope.rms = (sumx2/totx - TMath::Power(MCV.flat.slope.mean,2))/totx; if(MCV.flat.slope.rms>0) MCV.flat.slope.rms = TMath::Sqrt(MCV.flat.slope.rms); else MCV.flat.slope.rms = 0; } if(toty>0) { MCV.flat.inter.mean = sumy/toty; MCV.flat.inter.rms = (sumy2/toty - TMath::Power(MCV.flat.inter.mean,2))/toty; if(MCV.flat.inter.rms>0) MCV.flat.inter.rms = TMath::Sqrt(MCV.flat.inter.rms); else MCV.flat.inter.rms = 0; } if(MCV.flat.slope.rms<houghHist->GetXaxis()->GetBinWidth(1)) MCV.flat.slope.rms = houghHist->GetXaxis()->GetBinWidth(1); if(MCV.flat.inter.rms<houghHist->GetYaxis()->GetBinWidth(1)) MCV.flat.inter.rms = houghHist->GetYaxis()->GetBinWidth(1); MCV.ph.iZ = MinPlane - 2; if(ptotx>0) { MCV.ph.slope.mean = psumx/ptotx; MCV.ph.slope.rms = (psumx2/ptotx - TMath::Power(MCV.ph.slope.mean,2))/ptotx; if(MCV.ph.slope.rms>0) MCV.ph.slope.rms = TMath::Sqrt(MCV.ph.slope.rms); else MCV.ph.slope.rms = 0; } if(ptoty>0) { MCV.ph.inter.mean = psumy/ptoty; MCV.ph.inter.rms = (psumy2/ptoty - TMath::Power(MCV.ph.inter.mean,2))/ptoty; if(MCV.ph.inter.rms>0) MCV.ph.inter.rms = TMath::Sqrt(MCV.ph.inter.rms); else MCV.ph.inter.rms = 0; } if(MCV.ph.slope.rms<phHoughHist->GetXaxis()->GetBinWidth(1)) MCV.ph.slope.rms = phHoughHist->GetXaxis()->GetBinWidth(1); if(MCV.ph.inter.rms<phHoughHist->GetYaxis()->GetBinWidth(1)) MCV.ph.inter.rms = phHoughHist->GetYaxis()->GetBinWidth(1); MCV.flat.slope.peak = houghHist->GetXaxis()->GetBinCenter(xbin); MCV.flat.inter.peak = houghHist->GetYaxis()->GetBinCenter(ybin); MCV.ph.slope.peak = phHoughHist->GetXaxis()->GetBinCenter(phXbin); MCV.ph.inter.peak = phHoughHist->GetYaxis()->GetBinCenter(phYbin); if(maxVal < 3 || (TMath::Abs(MCV.flat.slope.mean)-MCV.flat.slope.rms<0 && TMath::Abs(MCV.flat.slope.peak)-MCV.flat.slope.rms<0)) { MCV.flat.slope.mean = 0; MCV.flat.inter.mean = maxWindowPos; MCV.flat.iZ = maxPlanePos; MCV.flat.slope.rms = 0; MCV.flat.inter.rms = 0; MCV.flat.slope.peak = 0; MCV.flat.inter.peak = maxWindowPos; } if( (TMath::Abs(MCV.ph.slope.mean)-MCV.ph.slope.rms<0 ) && (TMath::Abs(MCV.ph.slope.peak)-MCV.ph.slope.rms<0 )) { MCV.ph.slope.mean = 0; MCV.ph.inter.mean = maxWindowPos; MCV.ph.iZ = maxPlanePos; MCV.ph.slope.rms = 0; MCV.ph.inter.rms = 0; MCV.ph.slope.peak = 0; MCV.ph.inter.peak = maxWindowPos; } return MCV; }

void MakeClusterSS::CheckCurrSoln (  ClustVars &  cv  )

Definition at line 1677 of file MakeClusterSS.cxx. References clusterss::ClustVars::eslopen, clusterss::ClustVars::maxwid, printlevel, clusterss::ClustVars::RemoveWindow(), clusterss::ClustVars::shwmid, clusterss::ClustVars::slopen, clusterss::ClustVars::trustslopen, and windows. Referenced by ClusterWindows(). { if(printlevel) cout<<"CheckCurrSoln "<<endl; //first decide which variables to use based on the view // CHECK IF SOLUTION IS GOOD double usewid = cv.maxwid; if(printlevel) cout<<"Checkking solutions. use wid is "<<usewid<<" maxwid "<<cv.maxwid<<endl; for(unsigned int d = 0;d<windows.size();d++){ // cout<<"Looping over windows "<<d<<" out of "<<windows.size()<<endl; if(d>0 && (windows[d]->upper>0||windows[d]->lower>0)){ //there is a solution if(printlevel){ cout<<"Found a soln to check: "<<d<<" trust slope? "<<cv.trustslopen<<endl <<" cv.shwmid "<<cv.shwmid<<" upper tpos "<<windows[d]->upper_tpos <<" lower tpos "<<windows[d]->lower_tpos <<" usewid: "<<usewid<<endl <<" window width sq: "<<(windows[d]->upper_tpos-cv.shwmid)*(windows[d]->lower_tpos-cv.shwmid) <<" usewid sq: "<<(usewid/2.)*(usewid/2.)<<endl; } if(!cv.trustslopen && (((windows[d]->upper_tpos-cv.shwmid)*(windows[d]->lower_tpos-cv.shwmid) > (usewid/2.)*(usewid/2.)))){ //don't trust new slope //and window is too far away from centroid //remove these windows // cout<<"removing a soln"<<endl; cv.RemoveWindow(*windows[d]); if(printlevel==1) cout<<"removing a window at "<<windows[d]->plane<<endl; windows[d]->upper=-1; windows[d]->lower=-1; windows[d]->upper_tpos=-1000;// x 100 to make cm windows[d]->lower_tpos=-1000;// x 100 to make cm windows[d]->ph=0; windows[d]->type=0; windows[d]->phwidth=-100; windows[d]->centroid=-1000;// x 100 to make cm } else if(cv.trustslopen){ //slope from fit was ok //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 unsigned int i = 0; while(windows[i]->ph==0.&&i<=d){ if(printlevel==1) cout<<"Looping over i "<<i<<" d "<<d<<endl; i++; } if(i==d){ //TV: 7-30-09 this doesnt' make sense while(windows[i]->ph==0.&&i<windows.size()){ if(printlevel==1) cout<<"i==d"<<" i "<<i<<" d "<<d<<endl; i++; } } // cout<<" i is "<<i<<endl; int diff = windows[d]->plane-windows[i]->plane; if(printlevel==1) cout<<"diff is "<<diff<<endl; if(usewid<TMath::Abs(diff)*cv.eslopen){ //find appropriate cut value on window width usewid = TMath::Abs(diff)*cv.eslopen; // cout<<"in if usewid "<<usewid<<endl; } if(printlevel==1) cout<<"usewid "<<usewid<<endl; //expected centroid for plane: double expmid = (windows[i]->upper_tpos+windows[i]->lower_tpos)/2.+ (windows[d]->plane-windows[i]->plane)*cv.slopen; if(printlevel==1){ cout<<"slopen "<<cv.slopen<<" i upper tpos "<<windows[i]->upper_tpos<<" lower "<<windows[i]->lower_tpos<<endl; cout<<" planes "<<windows[d]->plane<<" "<<windows[i]->plane<<endl; cout<<"expmid "<<expmid<<endl; } if((windows[d]->upper_tpos-expmid)*(windows[d]->lower_tpos-expmid)>(usewid/2.)*(usewid/2.)){ if(printlevel==1) cout<<"too far away from expected middle"<<endl; //too far away from expected middle //get rid of this window //remove these windows cv.RemoveWindow(*windows[d]); if(printlevel==1) cout<<"removing a window at "<<windows[d]->plane<<endl; windows[d]->upper=-1; windows[d]->lower=-1; windows[d]->upper_tpos=-1000;// x 100 to make cm windows[d]->lower_tpos=-1000;// x 100 to make cm windows[d]->ph=0; windows[d]->type=0; windows[d]->phwidth=-100; windows[d]->centroid=-1000;// x 100 to make cm } }// trustslopen } //if there is a solution (d>0 etc. // cout<<"done checking this solution"<<endl; }// for d (plane index) // cout<<"all done with CheckingSolns"<<endl; //now vectors of windows are filled return; }

void MakeClusterSS::CleanUpAllContainers (   )

Definition at line 2617 of file MakeClusterSS.cxx. References CleanUpContainers(), xchglist, xclusterlist, ychglist, and yclusterlist. Referenced by Reco(), and ~MakeClusterSS(). { CleanUpContainers(); // cout<<"in clean up"<<endl; for(unsigned int i=0;i<xchglist.size();i++){ if(xchglist[i]){ delete xchglist[i]; xchglist[i]=0; } } // cout<<"Cleans xchglist "<<xchglist.size()<<endl; for(unsigned int i=0;i<ychglist.size();i++){ if(ychglist[i]){ delete ychglist[i]; ychglist[i]=0; } } // cout<<"Cleans ychglist "<<ychglist.size()<<endl; for(unsigned int i=0;i<xclusterlist.size();i++){ if(xclusterlist[i]){ delete xclusterlist[i]; xclusterlist[i]=0; } } // cout<<"Cleans xclusterlist "<<xclusterlist.size()<<endl; for(unsigned int i=0;i<yclusterlist.size();i++){ if(yclusterlist[i]){ delete yclusterlist[i]; yclusterlist[i]=0; } } // cout<<"Cleans yclusterlist "<<yclusterlist.size()<<endl; xchglist.clear(); ychglist.clear(); xclusterlist.clear(); yclusterlist.clear(); }

void MakeClusterSS::CleanUpContainers (   )

Definition at line 2663 of file MakeClusterSS.cxx. References plnclusters, plnwindows, and windows. Referenced by CleanUpAllContainers(), and Reco(). { for(unsigned int i=0;i<plnclusters.size();i++){ if(plnclusters[i]){ delete plnclusters[i]; plnclusters[i]=0; } } // cout<<"Cleans plnclusters "<<plnclusters.size()<<endl; for(unsigned int i=0;i<plnwindows.size();i++){ if(plnwindows[i]){ delete plnwindows[i]; plnwindows[i]=0; } } // cout<<"Cleans plnwindows "<<plnwindows.size()<<endl; // cout<<"about to Cleans windows "<<windows.size()<<endl; for(unsigned int i=0;i<windows.size();i++){ if(windows[i]){ delete windows[i]; windows[i]=0; } } // cout<<"Cleans windows "<<windows.size()<<endl; // cout<<"done with clean up"<<endl; plnclusters.clear(); plnwindows.clear(); windows.clear(); }

void MakeClusterSS::CleanUpWindows (   )

Definition at line 2509 of file MakeClusterSS.cxx. References numPreMaxPlanes, printlevel, and windows. Referenced by DoCluster(). { // Clean up the clustered window vector. // Check for longitudinal gaps within the proto-subshower // find new max plane based on the window ph // Remove windows after first gap on both sides of shower max plane // update begWin and endWin plane if(printlevel==1){ cout<<"Cleaning up Windows"<<endl; } bool foundnonzerowin=false; int newshwmax = -1; double maxph=0.; // sort(windows.begin(),windows.end()); std::vector<WindowSS *>::iterator winit = windows.begin(); while(winit!=windows.end()){ if(((*winit)->upper!=-1&&(*winit)->lower!=-1)&&!foundnonzerowin){ //we've found a real window; foundnonzerowin=true; } // if(!foundnonzerowin){ //get rid of any leading zero windows // windows.erase(winit); // continue; // } if((*winit)->ph>maxph){ //keep track of new maxph plane maxph = (*winit)->ph; newshwmax = (*winit)->plane; } winit++; } if(!foundnonzerowin) return; int newbeg=begWinPlane; int newend=endWinPlane; bool foundpregap=false; bool foundpostgap=false; winit = windows.begin(); while(winit!=windows.end()){ if((*winit)->upper==-1&&(*winit)->lower==-1){ //we've found a gap; if(numPreMaxPlanes>numPostMaxPlanes){ if((*winit)->plane<=newshwmax){ //its a pregap foundpregap=true; if((*winit)->plane>newbeg){ //find the biggest plane with a gap less than the newshw max newbeg = (*winit)->plane; } } else{ //its a postgap foundpostgap=true; if((*winit)->plane<newend){ //find the smallest plane with a gap bigger than the newshw max newend = (*winit)->plane; } } } else{ if((*winit)->plane>newshwmax){ //its a postgap foundpostgap=true; if((*winit)->plane<newend){ //find the smallest plane with a gap bigger than the newshw max newend = (*winit)->plane; } } else{ //its a pregap foundpregap=true; if((*winit)->plane>newbeg){ //find the biggest plane with a gap less than the newshw max newbeg = (*winit)->plane; } } } } winit++; } if(printlevel==1) cout<<"In clean up: newshwmax "<<newshwmax<<" new beg "<<newbeg<<" new end "<<newend<<endl; // now loop one last time to get rid of everything not between newbeg and newend winit = windows.begin(); while(winit!=windows.end()){ if(((*winit)->plane<=newbeg&&foundpregap)||((*winit)->plane>=newend&&foundpostgap)){ windows.erase(winit); continue; } winit++; } if(printlevel==1){ cout<<"Printing windows after Window Cleanup"<<endl; for(unsigned int i=0;i<windows.size();i++){ cout<<" i "<<i<<" plane "<<windows[i]->plane<<" upper "<<windows[i]->upper<<" lower "<<windows[i]->lower<<endl; } } return; }

bool MakeClusterSS::ClusterWindows (  geo::View_t  view  )

Definition at line 988 of file MakeClusterSS.cxx. References clusterss::ClustVars::AddNewWindow(), begWinPlane, BestHough(), geo::PlaneGeo::Cell(), clusterss::WindowSS::centroid, CheckCurrSoln(), endWinPlane, clusterss::ClustVars::eslopen, FindPlaneFromShwMax(), clusterss::houghResult::flat, geom, geo::Geometry::GetPlanesByView(), geo::CellGeo::HalfW(), HoughTransCluster(), clusterss::WindowSS::lower, clusterss::WindowSS::lower_tpos, clusterss::ClustVars::maxwid, clusterss::stats::mean, geo::Geometry::NextPlaneInView(), numPreMaxPlanes, clusterss::houghResult::ph, clusterss::WindowSS::ph, clusterss::WindowSS::phwidth, geo::Geometry::Plane(), clusterss::WindowSS::plane, plnwindows, clusterss::ClustVars::Print(), printlevel, clusterss::stats::rms, shwMaxPlane, clusterss::ClustVars::shwmid, clusterss::ClustVars::shwwid, clusterss::houghStats::slope, clusterss::ClustVars::slopen, StupidHoughNoWork(), clusterss::ClustVars::trustslopen, clusterss::WindowSS::type, clusterss::WindowSS::upper, clusterss::WindowSS::upper_tpos, and windows. Referenced by DoCluster(). { //now, cluster the windows that we've found // cout<<"In Cluster Windows "<<view<<endl; windows.clear(); if(printlevel==1){ cout<<"Size of plnwindows "<<plnwindows.size()<<endl; } if(plnwindows.size()<1){ return false; } Int_t numLongPlanes = numPreMaxPlanes+numPostMaxPlanes; // cout<<"In Cluster windows: numLongPlanes: "<<numLongPlanes<<endl; // if(numLongPlanes/2<2){ // if(begWinPlane==endWinPlane){ //if there's only one plane in the Long Window, //we need to find the window that contains the max pulse height //then fill the windows vector up //with this one plane window //TV: actually the original code does this if //the count of the two end planes, and all the planes in between (possibly in the other view) //is less than 3 (or more precisely (int)(N/2)<2) //this will do weird things at module boundaries, but I will code it //so it behaves the way it did in RecoSubShower if(endWinPlane-begWinPlane<3){ std::vector<WindowSS *>::iterator iter = plnwindows.begin(); float maxph=0.; WindowSS* bestwin=0; while(iter!=plnwindows.end()){ if((*iter)->ph>maxph){ maxph=(*iter)->ph; bestwin = (*iter); } iter++; } if(endWinPlane==begWinPlane){ WindowSS *bw = new WindowSS(*bestwin); windows.push_back(bw); } else{ const std::set<UInt_t>& plv=geom->GetPlanesByView(view); std::set<UInt_t>::iterator fp = plv.find(begWinPlane); std::set<UInt_t>::iterator ep = plv.find(endWinPlane); ep++; if(printlevel==1){ cout<<"begWInPlane is "<<begWinPlane<<" end "<<endWinPlane<<endl; cout<<"fp is "<<(*fp)<<" ep is "<<(*ep)<<endl; } while(fp!=ep){ WindowSS *bw = new WindowSS(*bestwin); bw->plane=(*fp); windows.push_back(bw); fp++; } } if(printlevel==1){ cout<<"About to return from special case of ClusterWindows"<<endl; cout<<"**And the Windows are:"<<endl; for(unsigned int i=0;i<windows.size();i++){ cout<<" plane "<<windows[i]->plane<<" upper "<<windows[i]->upper<<" lower "<<windows[i]->lower<<endl; } } return false; //we don't want to do the clean up step in this case } const geo::PlaneGeo gplane = geom->Plane(10); const geo::CellGeo gstrip = gplane.Cell(10); // IS THAT OK ? TRA LA RI LA (this must mean something in greek) Double_t STRIPWIDTHINMETRES=gstrip.HalfW()*2.; //otherwise cluster the transverse clusters ("windows") //to make a new vector of windows // Int_t numPreMaxPlanes = shwMaxPlane-begWinPlane+1;//plus 1 to get the end points right // Int_t numPostMaxPlanes = endWinPlane-shwMaxPlane+1; // cout<<"numPreMaxPlanes "<<numPreMaxPlanes<<" numPostMaxPlanes "<<numPostMaxPlanes<<endl; //match up windows plane by plane ClustVars cv; std::vector<clusterss::WindowSS*> houghcluster = HoughTransCluster(); if(houghcluster.size()==0){ // cout<<"Stupid Hough found nothing. Using alternate"<<endl; StupidHoughNoWork(view, cv); CheckCurrSoln(cv); if(printlevel==1) cout<<"Finished removing solns from Stupid"<<endl; } else { // cout<<"Hough Solution found"<<endl; //use hough solution as a starting point for (int d = 0;d<numLongPlanes;d++){ //the way this code was originally written, //null windows were added in the planes in the other view //rewriting so it doesn't worry about the other view //will be more efficient int pl = FindPlaneFromShwMax(view,d); // cout<<"on d: "<<d<<" pl is "<<pl<<endl; WindowSS *win = new WindowSS(); //set equal to -1, will only get set to a real plane if //there was a hough window selected that failed the slope check //in which case we want to try to find a new window for that plane //if there was no hough solution to begin with, we dont want to //try to find a new window later win->plane=-1; //loop over hough windows std::vector<WindowSS *>::iterator hcIter = houghcluster.begin(); std::vector<WindowSS *>::iterator hcEnd = houghcluster.end(); while (hcIter!=hcEnd){ if ((*hcIter)->plane==pl){ //if it's in current plane: win->plane = (*hcIter)->plane; win->upper = (*hcIter)->upper; win->lower = (*hcIter)->lower; win->upper_tpos = (*hcIter)->upper_tpos; win->lower_tpos = (*hcIter)->lower_tpos; win->ph = (*hcIter)->ph; win->type = (*hcIter)->type; // cout<<"there's a Hough "<<win->plane<<" upper "<<win->upper<<" lower "<<win->lower<<endl; cv.AddNewWindow(*win,STRIPWIDTHINMETRES); win->phwidth = cv.shwwid; win->centroid = cv.shwmid; } hcIter++; } windows.push_back(win); }//loop over d houghResult hough = BestHough(houghcluster); if(hough.ph.slope.mean!=0){ cv.slopen = hough.ph.slope.mean; cv.eslopen = hough.ph.slope.rms; cv.trustslopen = true; } else if(hough.flat.slope.mean!=0){ cv.slopen = hough.flat.slope.mean; cv.eslopen = hough.flat.slope.rms; cv.trustslopen = true; } else { cv.slopen = 0; cv.eslopen = 0; cv.trustslopen = false; } //check if solution is good: CheckCurrSoln(cv); if(printlevel==1) cout<<"Finished removing solns from Hough"<<endl; } if(printlevel==1){ cout<<"*****After first pass clustering CV holds: "<<endl; cv.Print(); cout<<"**And the Windows are:"<<endl; for(unsigned int i=0;i<windows.size();i++){ cout<<" plane "<<windows[i]->plane<<" upper "<<windows[i]->upper<<" lower "<<windows[i]->lower<<endl; } } if(cv.maxwid<2*STRIPWIDTHINMETRES) cv.maxwid = 2*STRIPWIDTHINMETRES; if(printlevel==1) cout<<"STRIPWIDTHINMETRES "<<STRIPWIDTHINMETRES<<endl; //Try to fill planes with no windows //First find centroid and width for +/- 1 planes if(printlevel==1){ cout<<"Trying to fill planes with no windows"<<endl; } for (int d = 0;d<numLongPlanes;d++){ // cout<<"in d loop again "<<d<<" numlongplanes "<<numLongPlanes<<endl; int pl = FindPlaneFromShwMax(view,d); if(printlevel==1) cout<<"d: "<<d<<" on pl "<<pl<<endl; double minwei = 20000; // x 100 to make cm //minimum distance to centroid (weight) if(windows[d]->upper==-1&&windows[d]->lower==-1&&windows[d]->plane!=-1){ if(printlevel==1) cout<<"No solution in entry d "<<d<<" plane "<<windows[d]->plane<<" pl "<<pl<<endl; //does this plane have a solution? //no double pmid = -1000; // x 100 to make cm //centroid for next plane (plus) double mmid = -1000; // x 100 to make cm //centroid for prev plane (minus) double pwid = -100; // x 100 to make cm //width double mwid = -100; // x 100 to make cm int pwt = -1; // window type plus int mwt = -1; // minus bool useslopep = false; bool useslopem = false; bool foundp = false; bool foundm = false; for (int c = 0;c<numLongPlanes;c++){ //look one plane ahead if(windows[c]->plane==(int)geom->NextPlaneInView(pl,1)) { if(windows[c]->upper>=0){ //if upper is not -1, then we have a real entry pmid = windows[c]->centroid; if (cv.trustslopen) useslopep = true; pwid = windows[c]->phwidth; pwt = windows[c]->type; foundp=true; } } //and one plane back if(windows[c]->plane==(int)geom->NextPlaneInView(pl,-1)) { if(windows[c]->upper>=0){ //if upper is not -1, then we have a real entry mmid = windows[c]->centroid; if (cv.trustslopen) useslopem = true; mwid = windows[c]->phwidth; mwt = windows[c]->type; foundm=true; } } } //Find compatible windows double usewid = cv.maxwid; // 10 became x 100 to turn in cm if(cv.maxwid>0){ if((!foundp&&!foundm)|| (!useslopem&&TMath::Abs(cv.shwmid-mmid)>2*cv.maxwid)|| (!useslopep&&TMath::Abs(cv.shwmid-pmid)>2*cv.maxwid)|| (useslopem&&(TMath::Abs(cv.shwmid-mmid-(geom->NextPlaneInView(pl,-1))*cv.slopen)>3*cv.maxwid))|| (useslopep&&(TMath::Abs(cv.shwmid-pmid-(shwMaxPlane-geom->NextPlaneInView(pl,1))*cv.slopen)>3*cv.maxwid))){ //do I really need this if? //either I didn't find a next plane or a previous plane //or I don't trust the minus slope and the centroid of the next plane is further way from //current shower middle than 2*the maximum width //or I don't trust he positive slope and the centroid is further away from current middle //than 2*maximum width //or I do trust the slope (in either view) and the predicted middle is further away than 3 times the maxwie minwei = 20000;// x 100 to make cm while(usewid<=2*cv.maxwid){ std::vector<clusterss::WindowSS *>::iterator winIter = plnwindows.begin(); int tempcount=0; //loop over the transverse plane windows while(winIter!=plnwindows.end()){ tempcount++; if((*winIter)->plane==pl){ double justmid = cv.shwmid; if(printlevel==1){ cout<<"justmid "<<justmid<<" pmid "<<pmid<<" mmid "<<mmid<<endl; } //10 became x 100 to turn in cm if(foundp) justmid = pmid; if(foundm&&(pl>shwMaxPlane||justmid==cv.shwmid)) justmid = mmid; if(printlevel==1){ cout<<"Will I add this window? upper tpos "<<(*winIter)->upper_tpos<<" justmid "<<justmid <<" lower tpos "<<(*winIter)->lower_tpos<<" usewid "<<usewid<<" minwei "<< minwei<<endl; cout<<"first condition: "<<((*winIter)->upper_tpos-justmid)*((*winIter)->lower_tpos-justmid) <<" usewid/2^2"<<(usewid/2.)*(usewid/2.) <<" midpt "<<TMath::Abs(((*winIter)->upper_tpos+(*winIter)->lower_tpos)/2.-justmid)<<endl; } if ((((*winIter)->upper_tpos-justmid)*((*winIter)->lower_tpos-justmid)<=(usewid/2.)*(usewid/2.)) && TMath::Abs(((*winIter)->upper_tpos + (*winIter)->lower_tpos)/2.-justmid)<minwei){ if(printlevel==1) cout<<"YES!! Case 1"<<endl; //this window is the right size around justmid //and average distance from justmid is below minimum //reset minimum //this is an acceptible window for this plane, so put it in the real window vector //TV: we don't increment the info in the ClustVars holder, though minwei = TMath::Abs(((*winIter)->upper_tpos + (*winIter)->lower_tpos)/2.-justmid); windows[d]->upper = (*winIter)->upper; windows[d]->lower = (*winIter)->lower; windows[d]->upper_tpos = (*winIter)->upper_tpos; windows[d]->lower_tpos = (*winIter)->lower_tpos; windows[d]->ph = (*winIter)->ph; windows[d]->type = (*winIter)->type; windows[d]->centroid = (*winIter)->centroid; windows[d]->phpos = (*winIter)->phpos; windows[d]->phwidth = (*winIter)->phwidth; } } winIter++; }//end loop over plane windows usewid+=cv.maxwid; } } else{ usewid = cv.maxwid; minwei = 20000;// x 100 to make cm while(usewid<=3*cv.maxwid){ std::vector<clusterss::WindowSS *>::iterator winIter = plnwindows.begin(); while(winIter!=plnwindows.end()){ if((*winIter)->plane==pl){ // 10 became x 100 to turn in cm if (!useslopem&&foundm&&(((*winIter)->upper_tpos-mmid)*((*winIter)->lower_tpos-mmid)<= (usewid/2.)*(usewid/2.)) && TMath::Abs(((*winIter)->upper_tpos+(*winIter)->lower_tpos)/2.-mmid)<minwei){ //this window is the right size around justmid //and average distance from justmid is below minimum //reset minimum //this is an acceptible window for this plane, so put it in the real window vector //TV, again we didn't increment CV info. if(printlevel==1) cout<<"YES!! Case 2"<<endl; minwei = TMath::Abs(((*winIter)->upper_tpos + (*winIter)->lower_tpos)/2.-mmid); windows[d]->upper = (*winIter)->upper; windows[d]->lower = (*winIter)->lower; windows[d]->upper_tpos = (*winIter)->upper_tpos; windows[d]->lower_tpos = (*winIter)->lower_tpos; windows[d]->ph = (*winIter)->ph; windows[d]->type = (*winIter)->type; windows[d]->centroid = (*winIter)->centroid; windows[d]->phpos = (*winIter)->phpos; windows[d]->phwidth = (*winIter)->phwidth; } // 10 became x 100 to turn in cm else if (!useslopep&&foundp&&(((*winIter)->upper_tpos-pmid)*((*winIter)->lower_tpos-pmid) <=(usewid/2.)*(usewid/2.)) && TMath::Abs(((*winIter)->upper_tpos + (*winIter)->lower_tpos)/2.-pmid)<minwei){ //this window is the right size around justmid //and average distance from justmid is below minimum //reset minimum //this is an acceptible window for this plane, so put it in the real window vector if(printlevel==1) cout<<"YES!! Case 3"<<endl; minwei = TMath::Abs(((*winIter)->upper_tpos + (*winIter)->lower_tpos)/2.-pmid); windows[d]->upper = (*winIter)->upper; windows[d]->lower = (*winIter)->lower; windows[d]->upper_tpos = (*winIter)->upper_tpos; windows[d]->lower_tpos = (*winIter)->lower_tpos; windows[d]->ph = (*winIter)->ph; windows[d]->type = (*winIter)->type; windows[d]->centroid = (*winIter)->centroid; windows[d]->phpos = (*winIter)->phpos; windows[d]->phwidth = (*winIter)->phwidth; } // 10 became x 100 to turn in cm else if (useslopem&&foundm){ Int_t m=pl-geom->NextPlaneInView(pl,-1); if(usewid<m*cv.eslopen) usewid = m*cv.eslopen; if(((*winIter)->upper_tpos-mmid-m*cv.slopen) * ((*winIter)->lower_tpos-mmid-m*cv.slopen) <= (usewid/2.)*(usewid/2.) &&TMath::Abs(((*winIter)->upper_tpos + (*winIter)->lower_tpos)/2.-mmid-m*cv.slopen)<minwei){ if(printlevel==1) cout<<"YES!! Case 4"<<endl; minwei = TMath::Abs(((*winIter)->upper_tpos + (*winIter)->lower_tpos)/2.-mmid-m*cv.slopen); windows[d]->upper = (*winIter)->upper; windows[d]->lower = (*winIter)->lower; windows[d]->upper_tpos = (*winIter)->upper_tpos; windows[d]->lower_tpos = (*winIter)->lower_tpos; windows[d]->ph = (*winIter)->ph; windows[d]->type = (*winIter)->type; windows[d]->centroid = (*winIter)->centroid; windows[d]->phpos = (*winIter)->phpos; windows[d]->phwidth = (*winIter)->phwidth; } } // 10 became x 100 to turn in cm else if (useslopep&&foundp){ Int_t p=geom->NextPlaneInView(pl,1)-pl; if(usewid<p*cv.eslopen) usewid = p*cv.eslopen; if(((*winIter)->upper_tpos-pmid+p*cv.slopen) * ((*winIter)->lower_tpos-pmid+p*cv.slopen) <= (usewid/2.)*(usewid/2.) &&TMath::Abs(((*winIter)->upper_tpos + (*winIter)->lower_tpos)/2.-pmid+p*cv.slopen)<minwei){ if(printlevel==1) cout<<"YES!! Case 5"<<endl; minwei = TMath::Abs(((*winIter)->upper_tpos + (*winIter)->lower_tpos)/2.-pmid+p*cv.slopen); windows[d]->upper = (*winIter)->upper; windows[d]->lower = (*winIter)->lower; windows[d]->upper_tpos = (*winIter)->upper_tpos; windows[d]->lower_tpos = (*winIter)->lower_tpos; windows[d]->ph = (*winIter)->ph; windows[d]->type = (*winIter)->type; windows[d]->centroid = (*winIter)->centroid; windows[d]->phpos = (*winIter)->phpos; windows[d]->phwidth = (*winIter)->phwidth; } } } winIter++; } usewid+=cv.maxwid; } } } cv.AddNewWindow(*windows[d],STRIPWIDTHINMETRES); } //if(newwinpl0>=0) windows[d]->plane=pl; //set the plane so we can sort later } // solution or not??? // cout<<"End of ClusterWindows. We have "<<windows.size()<<endl; // cout<<"Printing windows found by ClusterWindows"<<endl; // for(unsigned int i=0;i<windows.size();i++){ // cout<<" i "<<i<<" plane "<<windows[i]->plane<<" upper "<<windows[i]->upper<<" lower "<<windows[i]->lower<<endl; // } return true; //if we get this far, we do want to clean up the windows found in the next function }

bool MakeClusterSS::DoCluster (  geo::View_t  view  )

Definition at line 294 of file MakeClusterSS.cxx. References CleanUpWindows(), ClusterWindows(), FindAllTransverseWindows(), FindLongWindow(), MakePlaneClusters(), printlevel, PruneCellList(), and windows. Referenced by Reco(). { //group the strips into clusters in individial planes MakePlaneClusters(view); //find an appropriate window for longitudnal clustering FindLongWindow(view); //find transverse windows in all planes FindAllTransverseWindows(view); //group the transverse windows into actual clusters bool docleanup=ClusterWindows(view); if(docleanup){ if(printlevel==1){ cout<<"About to clean up windows"<<endl; } //only clean up the windows in certain circumstances //that are determined in cluster windows CleanUpWindows(); } if(printlevel==1){ cout<<"About to start filling clusters and PRUNING"<<endl; cout<<"I have "<<windows.size()<<" windows"<<endl; for(unsigned int i=0;i<windows.size();i++){ cout<<"Plane: "<<windows[i]->plane<<" upper "<<windows[i]->upper<<" lower "<<windows[i]->lower<<endl; } } //get rid of the Cells that we've clustered bool foundnewcluster = PruneCellList(view); // if(foundnewcluster){ // cout<<" Found a new cluster "<<endl; // } if(foundnewcluster) return true; return false; }

void MakeClusterSS::FindAllTransverseWindows (  geo::View_t  view  )

Definition at line 636 of file MakeClusterSS.cxx. References rawdata::RawDigit::ADC(), begWinPlane, geo::PlaneGeo::Cell(), recobase::CellHit::Cell(), recobase::Cluster::Cell(), clusterss::WindowSS::centroid, endWinPlane, geom, geo::CellGeo::HalfW(), clusterss::WindowSS::lower, clusterss::WindowSS::lower_tpos, recobase::Cluster::NCell(), numPostMaxPlanes, numPreMaxPlanes, clusterss::WindowSS::ph, clusterss::WindowSS::phpos, clusterss::WindowSS::phwidth, clusterss::WindowSS::plane, geo::Geometry::Plane(), recobase::CellHit::Plane(), plnclusters, plnwindows, clusterss::WindowSS::Print(), printlevel, shwMaxPlane, clusterss::CellHitGrad::TGM(), clusterss::CellHitGrad::TGP(), clusterss::WindowSS::type, clusterss::WindowSS::upper, and clusterss::WindowSS::upper_tpos. Referenced by DoCluster(). { //a function that first looks for clusters of hits in each plane //records the upper and lower edges of possible windows in each plane //upper(lower) edges are defined either by a Cell that has no signal in the Cell above(below) //or by a cell that has less pulseheight in it than either of its neighbors //(in which case, the current cell becomes both upper and lower edge) //after finding transverse windows in a plane, it tries to find the smallest window //final result is a vector of WindowSS classes, which holds all the smallest windows in every plane //code cut and pasted (more or less) from RecoSubShower.cxx // cout<<"In FindAllTransverseWindows "<<view<<endl; plnwindows.clear(); //Start Transverse Clustering //Loop over all planes in long. window //Find transverse windows per plane if(plnclusters.size()==0) return; std::map<int,recobase::Cluster *>::iterator pBegMax=plnclusters.find(begWinPlane); std::map<int,recobase::Cluster *>::iterator pEndMax=plnclusters.find(endWinPlane); std::map<int,recobase::Cluster *>::iterator pbeg=plnclusters.find(shwMaxPlane); std::map<int,recobase::Cluster *>::iterator pend; std::map<int,recobase::Cluster *>::iterator piter; std::map<int,recobase::Cluster *>::iterator pitersav; // start at max plane. First, go upstream or downstream from there, // depending on which end of the window has more planes. Then, go // back and do the other end of the window. //figure out where to start and which way to go //if shwmax plane is 0 or 1, we can only go forwards //there's probably an end constraint too that I haven't taken the time to figure out //this works if shwMaxPlane isn't the begining or the end of the plane window vec //if numPreMaxPlanes==1 (i.e. we only have one plane) we want to go forwards at first //I think this will still fail if we end up with just one entry in the plncluster vector //but that isn't really a cluster anymore, is it? if(numPreMaxPlanes>numPostMaxPlanes&&!(shwMaxPlane==0||shwMaxPlane==1)&&numPreMaxPlanes>1){ pend=pBegMax; pend--; // cout<<"Going backwards at first"<<endl; } else{ pend=pEndMax; pend++; // cout<<"Going forwards at first"<<endl; } piter = pbeg; pitersav = pend; Int_t pass=1; // keep track of which pass we're on //keep track of upper and lower limits of strip windows 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 // a jagged loop over planes if(printlevel==1) cout<<"Starting jagged loop: starting on "<<piter->first<<" ending on "<<pend->first<<endl; while(piter != pend){ // find strip windows. //For this, we care either about: edges (type 1), isolated strips (type 2), pulseheight dips (type 0) winu.clear(); wind.clear(); winu_tpos.clear(); wind_tpos.clear(); wintu.clear(); wintd.clear(); // find dips and categorize them Int_t plane=piter->first; if(printlevel==1){ cout<<"Finding plane clusters: on plane "<<plane<<endl; } //get the cluster of strips in this plane recobase::Cluster *c=piter->second; //loop over CellHits in c int NCELL = c->NCell(view); for(int icell = 0;icell<NCELL;icell++){ const CellHitGrad *thiscell=0; bool foundcellup=false; bool foundcelldown=false; thiscell=dynamic_cast<const CellHitGrad *>(c->Cell(view,icell)); int cn = thiscell->Cell(); int pln = thiscell->Plane(); const CellHitGrad *cellup=0; cellup=dynamic_cast<const CellHitGrad *>(c->Cell(view,cn+1,pln)); if(cellup!=0&&cellup->Cell()==cn+1){ //we found a cell and it has the right cell number foundcellup=true; } const CellHitGrad *celldown=0; celldown=dynamic_cast<const CellHitGrad *>(c->Cell(view,cn-1,pln)); if(celldown!=0&&celldown->Cell()==cn-1){ //we found a cell and it has the right cell number foundcelldown=true; } // PH significance is not currently used //(only called dips if delta PH is significant given pe stats) if(!foundcellup&&!foundcelldown){ // cout<<"Found an isolated strip "<<plane<<" "<<cn<<endl; // isolated strip //this cell is both an upper and lower edge winu.push_back(cn); wind.push_back(cn); winu_tpos.push_back(thiscell->TCPos()); wind_tpos.push_back(thiscell->TCPos()); wintu.push_back(2); // isolated wintd.push_back(2); // isolated } else if(thiscell->TGP()<0&&thiscell->TGM()<0){ // cout<<"Found a dip strip "<<plane<<" "<<cn<<endl; //wow. this is the only plane I've needed the gradient info so far //and I have the up and down strips, so I could have just computed //the gradients here. //It would simplify things to get rid of the CellHitGrad objects all together // dip - higher PH on either side //this cell is both an upper and lower edge winu.push_back(cn); winu_tpos.push_back(thiscell->TCPos()); wintu.push_back(0); // type-0 (this is a dip) wind.push_back(cn); wind_tpos.push_back(thiscell->TCPos()); wintd.push_back(0); // type-0 (this is a dip) } else if(foundcelldown&&!foundcellup){ // cout<<"Found an upper edge"<<endl; // edge: there is a strip below, but not above //this cell is an upper edge winu.push_back(cn); winu_tpos.push_back(thiscell->TCPos()); wintu.push_back(1); // type-1 (this is an edge strip) } else if(foundcellup&&!foundcelldown){ // cout<<"Found a lower edge "<<plane<<" "<<cn<<endl; // edge: there is a strip above, but not below //this cell is a lower edge wind.push_back(cn); wind_tpos.push_back(thiscell->TCPos()); wintd.push_back(1); //type-1 } } //end loop over cells in this plane //use dips to form windows //if we didn't find the same number of upper edges as lower edges something isn't right //TV: 7-30-09 //this code will do what RecoSubShower was doing when I took it over, //but I suspect there might be a problem with possible windows at edges of the detector... //once I get this working and validated, I will have to come back to this and check it out //if this comment is still here, it means I haven't done this yet. //if you're looking for a bug that affects clustering on the edges of the detector, its probably here. unsigned int wins = winu.size(); if(wins!=wind.size()) { cerr<<"UP edges and down edges don't match. Returning"<<endl; return; } //temporarily hold window info std::vector<Int_t> win0(wins,0); // strip of upper dip std::vector<Int_t> win1(wins,0); // strip of lower dip std::vector<Float_t> win0_tpos(wins,0); // tpos of strip of upper dip std::vector<Float_t> win1_tpos(wins,0); // tpos of strip of lower dip std::vector<Float_t> win2(wins,0); // PH of window std::vector<Int_t> win3(wins,0); // ID of window: // 0-both dips; // 1 or 10-one dip,one gap; // 11-both gaps; // 22-isolated strip std::vector<Float_t> win4(wins,0); // Centroid of window in tpos std::vector<Float_t> win5(wins,0); // PH weighted tpos of window std::vector<Float_t> win6(wins,0); // PH weighted width of window //match up dips to form windows //loop over upper edges for(unsigned int 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 = geom->Plane(plane).Ncells(); //loop over lower edges for (unsigned int s = 0; s<wind.size(); s++){ //compute size of distance between upper edge and lower edge 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){ //either we have two edges and the upper and lower cells are not the same cell //or we have an isolated hit //and this diff is smaller than the any window diff we've already seen //TV 7-30-09: this will save the smallest difference between any upper and any down //but if we have an isolated hit, won't that be the one we always save? //perhaps we get the isolated ones first, then get the others in subsquent clustering? //not sure I completely understand what this code does //will investigate once I get it working... 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); } } //for s } // for w //now get PH of windows //and push_back windows into vectors // const geo::PlaneGeo gplane = geom->Plane(plane); const geo::PlaneGeo gplane = geom->Plane(10); const geo::CellGeo gstrip = gplane.Cell(10); // IS THAT OK ? TRA LA RI LA (this must mean something in greek) Double_t STRIPWIDTHINMETRES=gstrip.HalfW()*2.; //loop over windows for(unsigned int w = 0; w<wins; w++){ Double_t biggestStripPH = 0; Int_t n = 0; //loop over CellHits in c int NCELL = c->NCell(view); for(int icell = 0;icell<NCELL;icell++){ const CellHitGrad *thiscell=0; thiscell=dynamic_cast<const CellHitGrad *>(c->Cell(view,icell)); int celln=thiscell->Cell(); float mip; bool calibrated = thiscell->MIP(mip); if(!calibrated){ unsigned short adc; thiscell->ADC(adc); mip = adc/ADCTOMIP; } if(celln<=win0[w] && celln>=win1[w]){ //the current cell number is in this window //add this cell's pulseheight to the appropriate vector win2[w]+=mip; if(mip>biggestStripPH) { win4[w]=thiscell->TCPos(); biggestStripPH=mip; } win5[w]+= mip*thiscell->TCPos(); win6[w]+= mip*thiscell->TCPos()*thiscell->TCPos(); n++; } } //divide by total pulseheight to find ph weighted TPOS of this window win5[w] /= win2[w]; if(n!=1){ //more than one cell in this window //find phweighted TPOS^2 of the window win6[w] /= win2[w]; //subtract off the mean tpos to find the width squared win6[w] -= win5[w]*win5[w]; //then take the square root if(win6[w]>0) win6[w] = sqrt(win6[w]); else win6[w] = STRIPWIDTHINMETRES/sqrt(12.); } else{ win6[w] = STRIPWIDTHINMETRES/sqrt(12.); } //now fill WindowSS Class WindowSS *win = new WindowSS(); win->plane = plane; 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]; //and put in our vector of windows // cout<<"Pushing back a window: "<<endl; // win->Print(); plnwindows.push_back(win); } //end loop over window edges //figure out where we started and which way we're going in the jagged loop //increment iterators if(numPreMaxPlanes>numPostMaxPlanes){ // we started out in the beginning of the LongWindow if(piter == pBegMax) { // cout<<"Reached the end of the upstream bit"<<endl; // reached the upstream end, need to reverse pbeg++; // don't repeat max plane pend=pEndMax; pend++; // set end iterator just beyond window piter=pbeg; // restart pass++; // change direction } else{ //if its the first pass keep going backward, otherwise go forwards if(pass == 1) piter--; else piter++; } } else{ // we started out in the end of the LongWindow if(piter == pEndMax){ // cout<<"Reached the end of the downstream bit"<<endl; // reached the downstream end, need to reverse pbeg--; // don't repeat max plane pend=pBegMax; pend--; piter=pbeg; pass++; } else{ //if its the first pass keep going forwards, otherwise go backwards if(pass == 1) piter++; else piter--; } } pitersav = piter; //TV: 7-30-09, don't know why we need the next line, will comment out for now //if(winvec.size() > 3000 ) break; // !!!!!!!!!!!! FIX THIS!!!!!!!!!!!! } if(printlevel==1){ cout<<"all done finding transverse windows. Found "<<plnwindows.size()<<endl; cout<<"Those windows are: "<<endl; for(unsigned int i=0;i<plnwindows.size();i++){ cout<<"plnwindows plane "<<i<<endl; plnwindows[i]->Print(); } } }

void MakeClusterSS::FindLongWindow (  geo::View_t  view  )

Definition at line 456 of file MakeClusterSS.cxx. References begPlane, begWinPlane, endPlane, endWinPlane, event, geom, geo::Geometry::GetPlanesByView(), numPostMaxPlanes, numPreMaxPlanes, plnclusters, pLongWindowGapPlaneCut, printlevel, run, and shwMaxPlane. Referenced by DoCluster(). { //find the set of planes with the largest pulseheight per plane //fill begWinPlane and endWinPlane with the start and end of this range of planes //code more or less cut and pasted from RecoSubShower::FindMaxWindow // cout<<"In find long window "<<view<<endl; //TV: this code can not be doing the right thing for 1 plane windows //modifying it just go on if begPlane and endPlane are the same if(begPlane==endPlane){ begWinPlane = begPlane; endWinPlane = endPlane; shwMaxPlane = begPlane; numPreMaxPlanes=1; numPostMaxPlanes=1; return; } // cout<<"Doing work over plnclusters of size "<<plnclusters.size()<<endl; // cout<<"beginning plane: "<<begPlane<<" ending plane "<<endPlane<<endl; const std::set<UInt_t>& plv=geom->GetPlanesByView(view); // std::set<UInt_t>::iterator fp = plv.begin(); // std::set<UInt_t>::iterator ep = plv.end(); // cout<<"Printing all planes in this view:"<<endl; // int cc=0; // while(fp!=ep){ // cout<<"This plane "<<cc<<" is "<<*fp<<endl; // cc++; // fp++; // } // ep--; // cout<<"Got planes from geom: first plane this view: "<<*fp<<" last: "<<*ep<<endl; // cout<<" tried to get planes from geom. size: "<<plv.size()<<endl; if (plv.size() == 0) assert(0); std::set<UInt_t>::iterator itPlvBeg=plv.find(begPlane); std::set<UInt_t>::iterator itPlvEnd=plv.find(endPlane); itPlvEnd++; // pull end back if it's gone past THE end if(itPlvEnd==plv.end()) itPlvEnd--; //Find highest strip density longitudinal window int maxbeg = 0; //first boundary plane of max window int maxend = 0; Float_t maxDensity = 0; Float_t density = 0; Bool_t inAWindow=false; Bool_t currZero=true; Int_t contigZero = 0; Int_t totalPlaneIndex = 0; // keep track of number of planes, since the // plane ID number is not a smooth function of total plane Int_t startPlane = -1; Int_t finPlane = 0; Int_t lastSeenPlane = 0; Int_t firstSeenTotalIndex = 0; Int_t lastSeenTotalIndex = 0; for(std::set<UInt_t>::iterator iter=itPlvBeg; iter!=itPlvEnd; iter++){ totalPlaneIndex++; const Int_t plane=*iter; //does this plane have activity above threshold? currZero=true; if(plnclusters.find(plane) != plnclusters.end()){ float mip=(plnclusters.find(plane))->second->Q(); if(mip>pLongWindowMipCut){ // cout<<"Not in a gap. On Plane "<<plane<<endl; //TV 8-3-09: Changed the following blocks of code //pretty sure that original RecoSubShower had a logical flaw here //that would mean we miss some windows of 1 plane long //is this the first plane of a new window? if (!inAWindow){ //we weren't in a window before //reset the start and the density startPlane=plane; firstSeenTotalIndex=totalPlaneIndex; inAWindow=true; density=0; } lastSeenPlane=plane; lastSeenTotalIndex=totalPlaneIndex; currZero=false; // keep track of number of contiguous empty planes contigZero=0; density+=mip; // cout<<"Current density "<<density<<" lastSeenTotalIndex "<<lastSeenTotalIndex //<<" firstSeenTotalIndex "<<firstSeenTotalIndex<<" in a window "<<inAWindow //<<" startPlane "<<startPlane<<" finPlane "<<finPlane<<endl; } } if(currZero){ contigZero++; //cout<<"In a gap of length "<<contigZero<<" Plane is "<<plane<<endl; } // if there have been more than a certain number of empty planes, // close off previous window if(contigZero>=pLongWindowGapPlaneCut && finPlane<startPlane){ // cout<<"This gap has gotten longer than the threshold: "<<contigZero<<endl; finPlane=lastSeenPlane; Int_t npln=lastSeenTotalIndex-firstSeenTotalIndex+1; inAWindow=false; density/=float(npln); //cout<<"Window this gap ends had "<<npln<<" planes and a density of "<<density<<endl; if(density>maxDensity){ // check for max ph density window maxbeg=startPlane; maxend=finPlane; maxDensity=density; } } } // all planes in view (within this event) // if (!currZero && finPlane<startPlane){ // there was an unfinished window at the end //TV: 8-3-09 think the above is not sufficient to end all started windows //instead, just check whether or not we think we're in a window if(inAWindow){ // there was an unfinished window at the end //cout<<"Out of loop, checking unfinished window"<<endl; finPlane=lastSeenPlane; Int_t npln=lastSeenTotalIndex-firstSeenTotalIndex+1; inAWindow=false; density/=float(npln); //cout<<"Window this gap ends had "<<npln<<" planes and a density of "<<density<<endl; if (density>maxDensity) // check for max ph density window { maxbeg=startPlane; maxend=finPlane; maxDensity=density; } } // cout<<"Max density window: "<<maxDensity<<" start "<<maxbeg<<" end "<<maxend<<endl; //now find the plane of maximum pulse height in our highest density window int maxpl=-1; Double_t maxPlanePH = 0; itPlvBeg=plv.find(maxbeg); itPlvEnd=plv.find(maxend); itPlvEnd++; // pull end back if it's gone past THE end if(itPlvEnd==plv.end()) itPlvEnd--; for(std::set<UInt_t>::iterator iter=itPlvBeg; iter!=itPlvEnd; iter++){ const Int_t plane=*iter; if(plnclusters.find(plane) != plnclusters.end()){ float mip=(plnclusters.find(plane))->second->Q(); if(mip>maxPlanePH){ maxpl = plane; maxPlanePH = mip; } } } begWinPlane=maxbeg; endWinPlane=maxend; shwMaxPlane=maxpl; if(printlevel==1){ cout<<"Found a LongWindow "<<view<<" run "<<run<<" event "<<event<<endl; cout<<"begWinPlane "<<maxbeg<<" endWinPlane "<<maxend<<" shwMaxPlane "<<maxpl<<endl; } numPreMaxPlanes=0; numPostMaxPlanes=0; for(std::set<UInt_t>::iterator iter=itPlvBeg; iter!=itPlvEnd; iter++){ const int thisplane = *iter; if(thisplane<=maxpl) numPreMaxPlanes++; else numPostMaxPlanes++; } if(printlevel==1){ cout<<"numPreMaxPlanes "<<numPreMaxPlanes<<" Post "<<numPostMaxPlanes<<endl; } }

int MakeClusterSS::FindPlaneFromShwMax (  geo::View_t  view, int  away )

Definition at line 2451 of file MakeClusterSS.cxx. References geom, geo::Geometry::GetPlanesByView(), numPostMaxPlanes, numPreMaxPlanes, and shwMaxPlane. Referenced by ClusterWindows(), and StupidHoughNoWork(). { //TV: I've changed this code so it //will no longer put zeros in for planes in the view we're not working on //this code may be slow, in which case, I'll have to rethink this //goal is to order planes starting with shwMax, //if pre is larger than post, step backwards through the planes, till you get to beg plane //then step forwards after shwMax until you get to end Plane //if post is larger than pre, step forwards till you get to end //then step backwards until you get to beg Plane unsigned int pl=0; const std::set<UInt_t>& plv=geom->GetPlanesByView(view); std::set<UInt_t>::iterator fp = plv.find(shwMaxPlane); int d =away; if(numPreMaxPlanes>numPostMaxPlanes){ if(away<numPreMaxPlanes){ // cout<<"away<=NumPreMaxPlanes"<<endl; while(d){ fp--; d--; } } else{ // cout<<"away is bigger than Premaxplanes"<<endl; while(d>=numPreMaxPlanes){ fp++; d--; } } } else{ // cout<<"Post is bigger than pre"<<endl; if(away<=numPostMaxPlanes){ // cout<<"away is less than numpost"<<endl; while(d){ fp++; d--; } } else{ // cout<<"away is bigger than numpost"<<endl; while(d>=numPostMaxPlanes+1){ fp--; d--; } } } pl=*fp; // cout<<"Returning "<<pl<<endl; return pl; }

std::vector< clusterss::WindowSS * > MakeClusterSS::HoughTransCluster (   )

Definition at line 1863 of file MakeClusterSS.cxx. References BestHough(), geo::PlaneGeo::Cell(), clusterss::houghResult::flat, geom, geo::CellGeo::HalfW(), clusterss::houghStats::inter, clusterss::houghStats::iZ, clusterss::stats::mean, clusterss::stats::peak, clusterss::houghResult::ph, geo::Geometry::Plane(), plnwindows, pMinHoughPlanes, clusterss::stats::rms, and clusterss::houghStats::slope. Referenced by ClusterWindows(). { // cout<<"In HoughTransCluster: view: "<<view<<endl; std::vector<WindowSS *> cluster; std::vector<WindowSS *> emptyCluster; //copy all windows into nonCluster: //also get best vertex estimate in first plane Int_t vertexPlane = 1985; Double_t vertexPlanePH = 0; Double_t vertexTPos = 0; Double_t vertexTPosWidth = 0; Int_t vertexPlane1 = 1985; Double_t vertexPlanePH1 = 0; Double_t vertexTPos1 = 0; Double_t vertexTPosWidth1 = 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 totStrips = 0; Double_t maxTpos = -800.0;// x 100 to make cm Double_t minTpos = 800.0; Int_t lastPlane = 0; Double_t PHavgTpos = -800.0;// x 100 to make cm Double_t totPH1 = 0.; std::vector<WindowSS *>::iterator wIt = plnwindows.begin(); std::vector<WindowSS *>::iterator wEn = plnwindows.end(); //loop over the transverse windows in each plane while(wIt!=wEn){ //save a copy of this window if((*wIt)->plane<vertexPlane){ //save info on the most upstream plane vertexPlane1 = vertexPlane; vertexPlanePH1 = vertexPlanePH; vertexTPos1 = vertexTPos; vertexTPosWidth1 = vertexTPosWidth; if((*wIt)->ph>1.0) { //save info on the most upstream plane with pulseheight > 1 mip //TV: this will depend on MIP calibration... vertexPlane = (*wIt)->plane; vertexPlanePH = (*wIt)->ph; vertexTPos = (*wIt)->phpos; vertexTPosWidth = (*wIt)->phwidth; // const geo::PlaneGeo gplane = geom->Plane((*wIt)->plane); const geo::PlaneGeo gplane = geom->Plane(10); const geo::CellGeo gstrip = gplane.Cell(10); Double_t STRIPWIDTHINMETRES=gstrip.HalfW()*2.; if(vertexTPosWidth<STRIPWIDTHINMETRES) vertexTPosWidth=STRIPWIDTHINMETRES; } else if(vertexPlanePH<1.0){ //pulse height less than one mip vertexPlane = (*wIt)->plane; vertexPlanePH = (*wIt)->ph; vertexTPos = (*wIt)->phpos; vertexTPosWidth = (*wIt)->phwidth; // const geo::PlaneGeo gplane = geom->Plane((*wIt)->plane); const geo::PlaneGeo gplane = geom->Plane(10); const geo::CellGeo gstrip = gplane.Cell(10); Double_t STRIPWIDTHINMETRES=gstrip.HalfW()*2.; if(vertexTPosWidth<STRIPWIDTHINMETRES) vertexTPosWidth=STRIPWIDTHINMETRES; } } else if((*wIt)->plane==vertexPlane){ if((*wIt)->ph>vertexPlanePH){ //if we find another window in the vertex plane //save the one with the max pulse height vertexPlanePH = (*wIt)->ph; vertexTPos = (*wIt)->phpos; vertexTPosWidth = (*wIt)->phwidth; // const geo::PlaneGeo gplane = geom->Plane((*wIt)->plane); const geo::PlaneGeo gplane = geom->Plane(10); const geo::CellGeo gstrip = gplane.Cell(10); Double_t STRIPWIDTHINMETRES=gstrip.HalfW()*2.; if(vertexTPosWidth<STRIPWIDTHINMETRES) vertexTPosWidth=STRIPWIDTHINMETRES; } } if((*wIt)->ph>1.0) { //save last plane if((*wIt)->plane>lastPlane) lastPlane = (*wIt)->plane; //save maximum upper edge if((*wIt)->upper_tpos>maxTpos) maxTpos = (*wIt)->upper_tpos; //save minimum lower edge if((*wIt)->lower_tpos<minTpos) minTpos = (*wIt)->lower_tpos; //compute pulseheight weighted position PHavgTpos+= (*wIt)->phpos * (*wIt)->ph; //add up the number of contained strips totStrips+= Double_t((*wIt)->upper - (*wIt)->lower + 1); //add up total pulseheight totPH1+= (*wIt)->ph; } wIt++; } //compute pulseheight weighted position PHavgTpos /= totPH1; if(TMath::Abs(vertexTPos-vertexTPos1)>2*(vertexTPosWidth+vertexTPosWidth1)&& TMath::Abs(vertexTPos-PHavgTpos)>TMath::Abs(vertexTPos1-PHavgTpos)){ //vertex plane transverse position too far away from average trans position //reset the values vertexPlane = vertexPlane1; vertexPlanePH = vertexPlanePH1; vertexTPos = vertexTPos1; vertexTPosWidth = vertexTPosWidth1; } //keep track of some quantities for quality checking at the end: Int_t lowPlane = 2000; 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<pMAXHOUGHITER) { // 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: houghResult MCV; Int_t badGradient[4] = {0}; //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){ // cout << "cluster.size()=0" << endl; MCV = BestHough(plnwindows); //verify that best hough passes reasonably close to vertex: //TV: huh? Does it ever do this? if(false){ // cout<<"HOW THE FUCK DID I GET IN HERE?"<<endl; double vertexY_simple = (vertexPlane - MCV.flat.iZ)* MCV.flat.slope.mean + MCV.flat.inter.mean; double vertexY_simple_peak = (vertexPlane - MCV.flat.iZ)* MCV.flat.slope.peak + MCV.flat.inter.peak; // MSG("SubShowerSR",Msg::kDebug) // cout<< "vertexY_simple = " // << vertexY_simple << endl; if( vertexY_simple < vertexTPos - 2*vertexTPosWidth || vertexY_simple > vertexTPos + 2*vertexTPosWidth || MCV.flat.slope.mean==0 ) badGradient[0] = 1; if( vertexY_simple_peak < vertexTPos - 2*vertexTPosWidth || vertexY_simple_peak > vertexTPos + 2*vertexTPosWidth || MCV.flat.slope.peak==0 ) badGradient[1] = 1; double vertexY_ph = (vertexPlane - MCV.ph.iZ)*MCV.ph.slope.mean + MCV.ph.inter.mean; double vertexY_ph_peak = (vertexPlane - MCV.ph.iZ)*MCV.ph.slope.peak + MCV.ph.inter.peak; // MSG("SubShowerSR",Msg::kDebug) // cout << "vertexY_ph = " // << vertexY_ph << endl; if( vertexY_ph < vertexTPos - 2*vertexTPosWidth || vertexY_ph > vertexTPos + 2*vertexTPosWidth || MCV.ph.slope.mean==0 ) badGradient[2] = 1; if( vertexY_ph_peak < vertexTPos - 2*vertexTPosWidth || vertexY_ph_peak > vertexTPos + 2*vertexTPosWidth || MCV.ph.slope.peak==0 ) badGradient[3] = 1; }//end of the statement that does nothing. } // cluster.size()==0, i.e., fisrt loop else { // cout << "cluster.size()="<< cluster.size() << endl; MCV = BestHough(cluster); } /* cout<<"MCV: flat slope "<<MCV.flat.slope.mean<<" inter "<<MCV.flat.inter.mean<<" iZ "<<MCV.flat.iZ<<endl; cout<<"MCV: flat slope "<<MCV.flat.slope.rms<<" inter "<<MCV.flat.inter.rms<<" iZ "<<MCV.flat.iZ<<endl; cout<<"MCV: flat slope "<<MCV.flat.slope.peak<<" inter "<<MCV.flat.inter.peak<<" iZ "<<MCV.flat.iZ<<endl; cout<<"MCV: ph slope "<<MCV.ph.slope.mean<<" inter "<<MCV.ph.inter.mean<<" iZ "<<MCV.ph.iZ<<endl; cout<<"MCV: ph slope "<<MCV.ph.slope.rms<<" inter "<<MCV.ph.inter.rms<<" iZ "<<MCV.ph.iZ<<endl; cout<<"MCV: ph slope "<<MCV.ph.slope.peak<<" inter "<<MCV.ph.inter.peak<<" iZ "<<MCV.ph.iZ<<endl; */ cluster.clear(); avePH = 0; aveWidth = 0; nStrips = 0; lowPlane = 2000; uppPlane = 0; std::vector<int> tempcluster[4]; std::vector<WindowSS*>::iterator winIter = plnwindows.begin(); std::vector<WindowSS*>::iterator winEnd = plnwindows.end(); Int_t counter = 0; //loop over original windows while(winIter!=winEnd){ double y[4] = {}; double y_upp[4] = {}; double y_low[4] = {}; //compute predicted t. pos from the flat (mean and peak) and ph (mean and peak) solns. y[0] = ((*winIter)->plane - MCV.flat.iZ)*MCV.flat.slope.mean + MCV.flat.inter.mean; y[1] = ((*winIter)->plane - MCV.flat.iZ)*MCV.flat.slope.peak + MCV.flat.inter.peak; y[2] = ((*winIter)->plane - MCV.ph.iZ)*MCV.ph.slope.mean + MCV.ph.inter.mean; y[3] = ((*winIter)->plane - MCV.ph.iZ)*MCV.ph.slope.peak + MCV.ph.inter.peak; // const geo::PlaneGeo gplane = geom->Plane((*winIter)->plane); const geo::PlaneGeo gplane = geom->Plane(10); const geo::CellGeo gstrip = gplane.Cell(10); Double_t STRIPWIDTHINMETRES=gstrip.HalfW()*2.; if(cluster.size()==0) { for(int ii=0;ii<4;ii++) { //compute 1 strip tolerances around predicted tpos's y_upp[ii] = y[ii] + 1.*STRIPWIDTHINMETRES; y_low[ii] = y[ii] - 1.*STRIPWIDTHINMETRES; } } else { for(int ii=0;ii<4;ii++) { //compute tolerances based on rms around predicted tpos's Double_t rms = ii<2 ? MCV.flat.inter.rms : MCV.ph.inter.rms; y_upp[ii] = y[ii] + rms; y_low[ii] = y[ii] - rms; } } 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])) { //predicted tpos is in the window, or within the tolerances //bad gradient is always zero, unless false sometimes equals true. //so keep track of this window if(badGradient[ii]==0) tempcluster[ii].push_back(counter); } } winIter++; counter++; } Int_t best_cluster = -1; Int_t max_size = 0; for(unsigned int ii=0;ii<4;ii++){ if((int)tempcluster[ii].size()>max_size) { //figure out which of the types of clusters has the most entries max_size = tempcluster[ii].size(); best_cluster = ii; } } if(best_cluster>=0 && best_cluster<=3) { std::vector<WindowSS*>::iterator winIter = plnwindows.begin(); std::vector<WindowSS*>::iterator winEnd = plnwindows.end(); std::vector<int>::iterator tempIter = tempcluster[best_cluster].begin(); std::vector<int>::iterator tempEnd = tempcluster[best_cluster].end(); counter = 0; //loop over plane windows while(winIter!=winEnd){ //loop over tmpcluster // cout<<"Testing a hough to push back in plane "<<(*winIter)->plane<<" upper "<<(*winIter)->upper<<" lower "<<(*winIter)->lower<<" counter "<<counter<<" (*tempiter )"<<(*tempIter)<<endl; if(tempIter!=tempEnd && (*tempIter) == counter) { //push back this window into the cluster // cout<<"Found a hough to push back in plane "<<(*winIter)->plane<<" upper "<<(*winIter)->upper<<" lower "<<(*winIter)->lower<<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; tempIter++; } winIter++; counter++; } } if(nStrips == oldNStrips) keepGoing = false; nLoops+=1; } if(nStrips<=0) { return emptyCluster; } // MSG("SubShowerSR",Msg::kDebug) // << "Returning Proto-SubShower from HoughTransCluster()" << endl; avePH/=nStrips; if(((uppPlane-lowPlane + 2)/2)<pMinHoughPlanes&&avePH<pMinHoughPH ) { return emptyCluster; } return cluster; }

void MakeClusterSS::MakePlaneClusters (  geo::View_t  view  )

Definition at line 340 of file MakeClusterSS.cxx. References rawdata::RawDigit::ADC(), begPlane, recobase::CellHit::Cell(), endPlane, clusterss::CellHitGrad::PHErr(), recobase::CellHit::Plane(), plnclusters, printlevel, clusterss::CellHitGrad::SetTGM(), clusterss::CellHitGrad::SetTGMe(), clusterss::CellHitGrad::TGM(), and clusterss::CellHitGrad::TGMe(). Referenced by DoCluster(). { // cout<<" In Make PlaneClusters "<<view<<endl; //first decide which variables to use based on the view std::vector<CellHitGrad*> *chglist; if(view==geo::kX){ //xview chglist=&xchglist; } else{ //yview chglist=&ychglist; } plnclusters.clear(); begPlane=-1; endPlane=-1; if(printlevel==1){ // cout<<" IN MakePlaneClusters "<<view<<endl; // cout<<"Working with strips: "<<endl; // for(unsigned int i=0;i<chglist->size();i++){ // (*chglist)[i]->Print(""); // } } //loop over the list of CellHits and put into map of clusters, indexed by plane number for(unsigned int ih=0;ih<chglist->size();ih++){ CellHitGrad *chg = (*chglist)[ih]; //look for the beginning plane if(begPlane==-1){ begPlane=chg->Plane(); } if(chg->Plane()<begPlane){ begPlane = chg->Plane(); } //look for the end plane if(endPlane==-1){ endPlane=chg->Plane(); } if(chg->Plane()>endPlane){ endPlane = chg->Plane(); } //loop over CellHits again to fill the extra information required by a CellHitGrad //look for adjacent cells //first look for one higher int cellp = chg->Cell()+1; //then one lower int cellm = chg->Cell()-1; for(unsigned int icell = 0;icell<chglist->size();icell++){ CellHitGrad *nxtchg = (*chglist)[icell]; if(nxtchg->Cell()==cellp){ //found one up float mip; bool calibrated = nxtchg->MIP(mip); if(!calibrated){ unsigned short adc; nxtchg->ADC(adc); mip = adc/ADCTOMIP; } chg->SetTGM(chg->TGM()+mip); chg->SetTGMe(chg->TGMe()+nxtchg->PHErr()*nxtchg->PHErr()); } //now one lower if(nxtchg->Cell()==cellm){ float mip; bool calibrated = nxtchg->MIP(mip); if(!calibrated){ unsigned short adc; nxtchg->ADC(adc); mip = adc/ADCTOMIP; } chg->SetTGM(chg->TGM()+mip); chg->SetTGMe(chg->TGMe()+nxtchg->PHErr()*nxtchg->PHErr()); } } //add this to a cluster for the given plane //if we havent seen this plane before, setup a cluster for it if(plnclusters.find(chg->Plane())==plnclusters.end()){ recobase::Cluster *c = new recobase::Cluster(); plnclusters.insert(pair<int,recobase::Cluster*>(chg->Plane(),c)); } //add this Cell to the plane cluster plnclusters[chg->Plane()]->Add(chg); // cout<<" *****Adding a Cell to a plane cluster! plane: "<<chg->Plane()<<endl; // chg->Print(); // cout<<" *****"<<endl; //add pulseheight to cluster Q double Q = plnclusters[chg->Plane()]->Q(); float mip; bool calibrated = chg->MIP(mip); if(!calibrated){ unsigned short adc; chg->ADC(adc); mip = adc/ADCTOMIP; } short unsigned int adc; chg->ADC(adc); plnclusters[chg->Plane()]->SetQ(Q+mip); } // cout<<" Put all Cells into Plane Clusters found "<<plnclusters.size()<<" of them"<<endl; if(printlevel==1){ cout<<"Leaving MakePlaneClusters "<<view<<endl; cout<<"Found "<<plnclusters.size()<<" planes"<<endl; cout<<"begPlane "<<begPlane<<" endPlane "<<endPlane<<endl; } }

void MakeClusterSS::MergeClusters (  geo::View_t  view  )

Definition at line 2444 of file MakeClusterSS.cxx. Referenced by Reco(). { if(view==geo::kX){} }

bool MakeClusterSS::PruneCellList (  geo::View_t  view  )

Definition at line 1778 of file MakeClusterSS.cxx. References rawdata::RawDigit::ADC(), recobase::Cluster::Add(), begWinPlane, recobase::CellHit::Cell(), endWinPlane, event, recobase::Cluster::NCell(), recobase::CellHit::Plane(), printlevel, recobase::Cluster::Q(), run, recobase::Cluster::SetQ(), and windows. Referenced by DoCluster(). { // cout<<"In PruneCellList "<<view<<endl; //first decide which variables to use based on the view std::vector<CellHitGrad*> *chglist; std::vector<recobase::Cluster *> *finalcluster; if(view==geo::kX){ //xview chglist=&xchglist; finalcluster = &xclusterlist; } else{ //yview chglist=&ychglist; finalcluster = &yclusterlist; } if(printlevel==1){ cout<<"Printing windows in PruneStep"<<endl; for(unsigned int i=0;i<windows.size();i++){ cout<<" i "<<i<<" plane "<<windows[i]->plane<<" upper "<<windows[i]->upper<<" lower "<<windows[i]->lower<<endl; } cout<<"at begining of prune, have "<<chglist->size()<<endl; cout<<"begWinPlane "<<begWinPlane<<" endwinPlane "<<endWinPlane<<endl; } recobase::Cluster *c = new recobase::Cluster(); std::vector<CellHitGrad *>::iterator iter = chglist->begin(); while(iter!=chglist->end()){ recobase::CellHit *stp = *iter; if(printlevel==1){ // cout<<"Checking if this cell is in a window"<<endl; // stp->Print(); } if(stp->Plane()>=begWinPlane && stp->Plane()<=endWinPlane){ std::vector<WindowSS *>::iterator winIter = windows.begin(); while(winIter!=windows.end()){ // cout<<"Checking if "<<(*winIter)->plane<<" equals "<<stp->Plane()<<endl; if((*winIter)->plane==stp->Plane()){ Int_t strp = stp->Cell(); // cout<<"IT DOES! Checkign if "<<strp<<" is between "<<(*winIter)->upper<<" and "<<(*winIter)->lower<<endl; if(strp<=(*winIter)->upper && strp>=(*winIter)->lower) { c->Add(stp); if(printlevel==1){ // cout<<"Added a cell to cluster "<<stp->Plane()<<" "<<strp<<endl; } double Q = c->Q(); float mip; bool calibrated = stp->MIP(mip); if(!calibrated){ unsigned short adc; stp->ADC(adc); mip = adc/ADCTOMIP; } c->SetQ(Q+mip); chglist->erase(iter); iter--; break; } } winIter++; } } iter++; } bool foundcluster = false; if(c->NCell()>0){ foundcluster=true; finalcluster->push_back(c); } if(printlevel==1){ cout<<"Filled another cluster for Run: "<<run<<" event "<<event<<endl; cout<<"Now we have "<<finalcluster->size()<<" clusters for view: "<<view<<endl; cout<<"New cluster has "<<c->NCell()<<" strips "<<endl; cout<<"strips left over "<<chglist->size()<<endl; } return foundcluster; }

jobc::Result MakeClusterSS::Reco (  edm::EventHandle &  evt  )  [virtual]

Reimplemented from jobc::Module. Definition at line 118 of file MakeClusterSS.cxx. References edm::EventHandle::Cal(), CleanUpAllContainers(), CleanUpContainers(), edm::EventHandle::DAQ(), DoCluster(), edm::EventHeader::Event(), event, evtcntr, geom, edm::EventHandle::Header(), geo::Geometry::Instance(), MergeClusters(), edm::EventHandle::Reco(), edm::EventHeader::Run(), run, SplitIntoViews(), xclusterlist, and yclusterlist. { evtcntr++; // printlevel=1; // cout<<"PrintLevel in Make Cluster "<<printlevel<<endl; //get the list of CellHits std::vector<const recobase::CellHit*> hitlist(0); try{ evt.Cal().Get("./MergeHits",hitlist); } catch(edm::Exception e){ cerr<<"Couldn't find the MergedCellHit list, trying for the original"<<endl; try{ evt.Cal().Get("./Hits",hitlist); } catch(edm::Exception e){ cerr<<"Couldn't find the Original CellHit list either, failing."<<endl; return jobc::kFailed; } } //now get a header std::vector<const hdr::Header*> head(0); try{ evt.DAQ().Get(".",head); } catch(edm::Exception e){ cerr<<"Couldn't get a head"<<endl; return false; } geom = &geo::Geometry::Instance(1,head[0]->DetectorID()); edm::EventHeader &eh = evt.Header(); run = eh.Run(); event = eh.Event(); // cout<<"*************************************************************************"<<endl; // cout<<"***IN RECO: On Run "<<run<<" "<<event<<" which is the "<<evtcntr<<" event we've seen"<<endl; //split initial CellHit list into separate views //fill initial values of CellHitGrad class SplitIntoViews(hitlist); //loop over the views (x or y) const int NVIEWS=2; for(int iv=0;iv<NVIEWS;iv++){ std::vector<CellHitGrad*> *chglist; geo::View_t view; if(iv==0){ //xview chglist=&xchglist; view=geo::kX; } else{ //yview chglist=&ychglist; view=geo::kY; } bool foundclus=true; //keep trying to find clusters until we run out of Cells //or we didn'y find another cluster while(chglist->size()!=0&&foundclus){ // cout<<"TRYING TO FIND CLUSTERS. "<<chglist->size()<<" strips left"<<endl; foundclus = DoCluster(view); // if(foundclus){ // cout<<"FOUND A CLUSTER, going around again"<<endl; // } CleanUpContainers(); } //if so configured, try to merge small clusters into bigger clusters if(pMergeClusters){ MergeClusters(view); } } //write the 2d clusters out if(evt.Reco().GetFolder("./SSClusterX")==0){ evt.Reco().MakeFolder("./SSClusterX"); } evt.Reco().PutVector(xclusterlist,"SSClusterX"); if(evt.Reco().GetFolder("./SSClusterY")==0){ evt.Reco().MakeFolder("./SSClusterY"); } evt.Reco().PutVector(yclusterlist,"SSClusterY"); // cout<<"About to clean up"<<endl; CleanUpAllContainers(); return jobc::kPassed; // kFailed if you want to fail the event }

void MakeClusterSS::SplitIntoViews (  std::vector< const recobase::CellHit * > &  hitlist  )

Definition at line 263 of file MakeClusterSS.cxx. References rawdata::RawDigit::ADC(), clusterss::CellHitGrad::SetPHErr(), clusterss::CellHitGrad::SetTGM(), clusterss::CellHitGrad::SetTGMe(), clusterss::CellHitGrad::SetTGP(), clusterss::CellHitGrad::SetTGPe(), recobase::CellHit::View(), xchglist, and ychglist. Referenced by Reco(). { for(unsigned int ih=0;ih<hitlist.size();ih++){ const recobase::CellHit* stp = hitlist[ih]; CellHitGrad *chg = new CellHitGrad(*stp); float mip; bool calibrated = stp->MIP(mip); if(!calibrated){ unsigned short adc; stp->ADC(adc); mip = adc/ADCTOMIP; } chg->SetTGP(mip); double pherr = 0.; if(mip>0) pherr = mip/sqrt(mip); chg->SetPHErr(pherr); chg->SetTGP(mip); chg->SetTGM(mip); chg->SetTGPe(pherr*pherr); chg->SetTGMe(pherr*pherr); if(stp->View()==geo::kX) xchglist.push_back(chg); else if(stp->View()==geo::kY) ychglist.push_back(chg); } return; }

bool MakeClusterSS::StupidHoughNoWork (  geo::View_t  view, ClustVars &  cv )

Definition at line 1405 of file MakeClusterSS.cxx. References clusterss::ClustVars::AddNewWindow(), geo::PlaneGeo::Cell(), clusterss::WindowSS::centroid, clusterss::ClustVars::chi2bn, clusterss::ClustVars::chi2tn, clusterss::ClustVars::eslopebn, clusterss::ClustVars::eslopen, clusterss::ClustVars::eslopetn, FindPlaneFromShwMax(), geom, geo::CellGeo::HalfW(), clusterss::WindowSS::lower, clusterss::WindowSS::lower_tpos, geo::Geometry::NextPlaneInView(), numPreMaxPlanes, clusterss::WindowSS::ph, clusterss::WindowSS::phwidth, clusterss::WindowSS::plane, geo::Geometry::Plane(), plnwindows, clusterss::ClustVars::shwmid, clusterss::ClustVars::shwwid, clusterss::ClustVars::slopebn, clusterss::ClustVars::slopen, clusterss::ClustVars::slopetn, clusterss::ClustVars::trustslopen, clusterss::WindowSS::type, clusterss::WindowSS::upper, clusterss::WindowSS::upper_tpos, and windows. Referenced by ClusterWindows(). { //TV: FIxes needed here. all of the ClustVars variables are doubles, so == comparisons won't work //found no houghcluster //so try another approach to clustering // cout<<"StupidHoughNoWork "<<view<<endl; //first decide which variables to use based on the view windows.clear(); const geo::PlaneGeo gplane = geom->Plane(10); const geo::CellGeo gstrip = gplane.Cell(10); // IS THAT OK ? TRA LA RI LA (this must mean something in greek) Double_t STRIPWIDTHINMETRES=gstrip.HalfW()*2.; Int_t numLongPlanes = numPreMaxPlanes+numPostMaxPlanes; // !!! change uint int // cout<<"Numlongplanes in Stupid "<<numLongPlanes<<endl; //loop over planes int nt = 0; //count of number of windows found int nb = 0; //also a count of the number of winows found for (int d = 0;d<numLongPlanes;d++){ //find the dth plane from shower max //before if there are more planes in the prewindow //after if there are more planes in the postwindow int pl = FindPlaneFromShwMax(view,d); // cout<<"In stupid, d: "<<d<<" pl "<<pl<<endl; WindowSS *win = new WindowSS(); //here, set the plane to be something real //if later we found no windows in this plane //we will want to try again, so we set the plane number win->plane=pl; // PH for current best matched window: double weight = 0.; //deviation of best matched window from prediction of //centroid for this plane: double weight1 = 1000000.; // 100 to make cm //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 std::vector<clusterss::WindowSS *>::iterator winIter = plnwindows.begin(); while(winIter!=plnwindows.end()){ if((*winIter)->plane==pl){ //we have a window in our plane of interest // cout<<"Found a plnwindow for plane "<<pl<<" upper: "<<(*winIter)->upper<<" lower "<<(*winIter)->lower<<endl; // cout<<"what is cv.shwmid? "<<cv.shwmid<<endl; loop++; double dev = cv.shwwid; //find the widest window in this plane if(dev<((*winIter)->upper_tpos-(*winIter)->lower_tpos+STRIPWIDTHINMETRES)) dev =(*winIter)->upper_tpos-(*winIter)->lower_tpos+STRIPWIDTHINMETRES; //if its too small, make it bigger (I don't know why) if(dev<2.*STRIPWIDTHINMETRES) dev = 2.*STRIPWIDTHINMETRES; // cout<<"dev is "<<dev<<endl; //set the prediction of the centriod for this plane //both with and without slope info if(cv.shwmid!=0){ mid0 = cv.shwmid; mid = cv.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((cv.shwmid==0 || (TMath::Abs((*winIter)->upper_tpos-mid)<dev && TMath::Abs((*winIter)->lower_tpos-mid)<dev)) && ((*winIter)->upper_tpos+(*winIter)->lower_tpos)/2.-mid0<weight1){ //shower centroid is currently zero //or plane half width is smaller than dev (upper minus lower) //and half width minus best guess at half width is less than allowable deviation //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 midp = 0; double midm = 0; //find windows in next plane and in previous plane //to check to see if this window matches windows adjacent //matches if a window in next(prev) plane //has upper an lower edges that are less far away from the middle of the long. cluster //than the biggest window in the current plane seen so far std::vector<WindowSS *>::iterator winIter1 = plnwindows.begin(); while(winIter1!=plnwindows.end()){ //first look at next plane if((*winIter1)->plane==(int)geom->NextPlaneInView(pl,1)){ //go to next plane in same view loopp++; if(cv.shwmid!=0){ //not shower max plane midp = cv.shwmid; if(TMath::Abs((*winIter1)->upper_tpos-midp)<dev ||TMath::Abs((*winIter1)->lower_tpos-midp)<dev){ //good match withinp++; } } else if(cv.shwmid==0){ //in shower max plane midp = ((*winIter)->upper_tpos+(*winIter)->lower_tpos)/2.; if(TMath::Abs((*winIter1)->upper_tpos-midp)<dev ||TMath::Abs((*winIter1)->lower_tpos-midp)<dev){ //good match withinp++; } } } // if (*winIter1)? //now previous plane if((*winIter1)->plane==(int)geom->NextPlaneInView(pl,-1)){ //go to previous same view plane loopm++; if(cv.shwmid!=0){ // not shower max plane midm = cv.shwmid; if(TMath::Abs((*winIter1)->upper_tpos-midm)<dev ||TMath::Abs((*winIter1)->lower_tpos-midm)<dev){ //good match withinm++; } } else if(cv.shwmid==0){ //shower max plane midm = ((*winIter)->upper_tpos+(*winIter)->lower_tpos)/2.; // cout<<"midm "<<midm<<endl; // cout<<" prev plane upper: "<<(*winIter1)->upper_tpos<<" lower "<<(*winIter1)->lower_tpos<<endl; if(TMath::Abs((*winIter1)->upper_tpos-midm)<dev ||TMath::Abs((*winIter1)->lower_tpos-midm)<dev){ //good match withinm++; } }// if shwmid!=0 }//if winIter1 winIter1++; } //done looking in planes ahead and behind the current one // cout<<"Did we find a good window? within p "<<withinp<<" m "<<withinm<<" loopp "<<loopp<<" loopm "<<loopm<<endl; if((withinp>0 && withinm>0)||(loopp==0 && withinm>0)||(loopm==0 && withinp>0)){ // cout<<"Found a window inplane "<<pl<<endl; //success! current window is a good solution within++; //save solution: win->plane = (*winIter)->plane; win->upper = (*winIter)->upper; win->lower = (*winIter)->lower; win->upper_tpos = (*winIter)->upper_tpos; win->lower_tpos = (*winIter)->lower_tpos; win->ph = (*winIter)->ph; win->type = (*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 //TV: 7-30-09 where is it checking the PH? Won't this just store the last window //found in the current plane? //TV: weight1 can be negative, which means we won't add a better cluster much of the time //TV: I also think the criterion is not more PH, but is supposed to be distance (maybe comment is an old one) } // if withinp } // if shmid=0, etc. } //if (winIter->plane==pl) winIter++; } //while (winIter!=winEnd) if(loop*within>0){ //got a matched window //we do this bit of the loop again if we find a hough solution //determine widest window: cv.AddNewWindow(*win,STRIPWIDTHINMETRES); win->phwidth = cv.shwwid; win->centroid = cv.shwmid; } // if (loop*within...) windows.push_back(win); if(windows[d]->ph>0){nt++; nb++;} } // for d (plane index) //if we found at least two windows 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++){ //fill some graphs of upper edge tpos vs. plane number //fill some graphs of lower edge tpos vs. plane number if(windows[d]->ph>0) { grat->SetPoint(t,windows[d]->plane,windows[d]->upper_tpos); grat->SetPointError(t,0,1./sqrt(windows[d]->ph)); t++; grab->SetPoint(b,windows[d]->plane,windows[d]->lower_tpos); grab->SetPointError(b,0,1./sqrt(windows[d]->ph)); b++; } } // for d (plane) //fit second order polys to these graphs (why 2nd order?) //find the slopes of the upper limits/plane grat->Fit("pol2","Q", 0, 2000); TF1 *ssPol1 = grat->GetFunction("pol2"); // if(!ssPol1 ){ // cout << "**********************************************************" <<endl // << "ECK NO POL! " <<endl; // } cv.slopetn = ssPol1->GetParameter(1); cv.eslopetn = ssPol1->GetParError(1); cv.chi2tn = ssPol1->GetChisquare(); grab->Fit("pol2","Q", 0, 2000); ssPol1 = grab->GetFunction("pol2"); // if(!ssPol1 ){ // cout << "**********************************************************" <<endl // << "ECK NO POL! " <<endl; // } //find the slopes of the lower limits/plane cv.slopebn = ssPol1->GetParameter(1); cv.eslopebn = ssPol1->GetParError(1); cv.chi2bn = ssPol1->GetChisquare(); ssPol1 = NULL; cv.slopen = (cv.slopetn+cv.slopebn)/2.; cv.eslopen = sqrt(cv.eslopetn*cv.eslopetn+cv.eslopebn*cv.eslopebn)/2.; // cout<<"fit window "<<cv.slopetn<<" "<<cv.slopebn // <<" " <<cv.slopen<<" "<<cv.eslopetn<<" " // <<cv.eslopebn<<" "<<cv.eslopen<<" " // <<cv.chi2tn/nt<<" "<<cv.chi2bn/nb<<" " // <<nt<<" "<<nb<<endl; if((cv.slopetn+cv.slopebn)!=0. && cv.slopetn*cv.slopebn!=0.){ //non zero slopes if(cv.chi2tn/nt<100. && cv.chi2bn/nb<100. && ((((cv.eslopetn/TMath::Abs(cv.slopetn)<0.1 && cv.eslopebn/TMath::Abs(cv.slopebn)<0.1 && cv.chi2tn/nt>0.1 && cv.chi2bn/nb>0.1)) && (TMath::Abs(cv.slopetn-cv.slopebn)/TMath::Abs(cv.slopetn+cv.slopebn)<0.25)) || TMath::Abs(cv.slopetn-cv.slopebn)==0)){ //not crazy chi2/ndf, errors on slopes<10% //rel. difference between up and down slopes<0.25% //or we got the same slopes //use average slope and error, trust them if(TMath::Abs(cv.slopetn-cv.slopebn)>cv.eslopen) cv.eslopen = TMath::Abs(cv.slopetn-cv.slopebn); cv.trustslopen = true; }// if chi2 } // if slope blah //get rid of the graphs delete grat; grat=NULL; delete grab; grab=NULL; } return true; //done if no hough solution }

void MakeClusterSS::Update (  const cfg::Config &  c  )  [virtual]

Implements cfg::Observer. Definition at line 91 of file MakeClusterSS.cxx. References pBestHoughMaxPHFrac, pLongWindowGapPlaneCut, pLongWindowMipCut, pMaxAngleDiff, pMAXHOUGHITER, pMaxOvlpLink, pMaxPDiff, pMaxTDiff, pMaxWidToLenP2, pMaxZDiff, pMergeClusters, pMinHoughPH, pMinHoughPlanes, and printlevel. { c("PrintLevel").Get(printlevel); c("MergeClusters").Get(pMergeClusters); c("MaxOvlpLink").Get(pMaxOvlpLink); c("MaxPDiff").Get(pMaxPDiff); c("MaxZDiff").Get(pMaxZDiff); c("MaxTDiff").Get(pMaxTDiff); c("MaxAngleDiff").Get(pMaxAngleDiff); c("MaxWidToLenP2").Get(pMaxWidToLenP2); c("LongWindowMipCut").Get(pLongWindowMipCut); c("LongWindowGapPlaneCut").Get(pLongWindowGapPlaneCut); c("MinHoughPlanes").Get(pMinHoughPlanes); c("MinHoughPH").Get(pMinHoughPH); c("BestHoughMaxPHFrac").Get(pBestHoughMaxPHFrac); c("MaxHoughIter").Get(pMAXHOUGHITER); }

---

Member Data Documentation

const double MakeClusterSS::ADCTOMIP = 1. [static]

Definition at line 49 of file MakeClusterSS.cxx.

int clusterss::MakeClusterSS::begPlane [private]

Definition at line 79 of file MakeClusterSS.h. Referenced by FindLongWindow(), and MakePlaneClusters().

int clusterss::MakeClusterSS::begWinPlane [private]

Definition at line 82 of file MakeClusterSS.h. Referenced by ClusterWindows(), FindAllTransverseWindows(), FindLongWindow(), and PruneCellList().

int clusterss::MakeClusterSS::endPlane [private]

Definition at line 80 of file MakeClusterSS.h. Referenced by FindLongWindow(), and MakePlaneClusters().

int clusterss::MakeClusterSS::endWinPlane [private]

Definition at line 83 of file MakeClusterSS.h. Referenced by ClusterWindows(), FindAllTransverseWindows(), FindLongWindow(), and PruneCellList().

int clusterss::MakeClusterSS::event [private]

Definition at line 78 of file MakeClusterSS.h. Referenced by FindLongWindow(), PruneCellList(), and Reco().

int clusterss::MakeClusterSS::evtcntr [private]

Definition at line 76 of file MakeClusterSS.h. Referenced by Reco().

geo::Geometry* clusterss::MakeClusterSS::geom [private]

Definition at line 75 of file MakeClusterSS.h. Referenced by BestHough(), ClusterWindows(), FindAllTransverseWindows(), FindLongWindow(), FindPlaneFromShwMax(), HoughTransCluster(), Reco(), and StupidHoughNoWork().

TH2D* clusterss::MakeClusterSS::houghHist [private]

Definition at line 115 of file MakeClusterSS.h. Referenced by BestHough(), and MakeClusterSS().

int clusterss::MakeClusterSS::numPostMaxPlanes [private]

Definition at line 86 of file MakeClusterSS.h. Referenced by FindAllTransverseWindows(), FindLongWindow(), and FindPlaneFromShwMax().

int clusterss::MakeClusterSS::numPreMaxPlanes [private]

Definition at line 85 of file MakeClusterSS.h. Referenced by CleanUpWindows(), ClusterWindows(), FindAllTransverseWindows(), FindLongWindow(), FindPlaneFromShwMax(), and StupidHoughNoWork().

Double_t clusterss::MakeClusterSS::pBestHoughMaxPHFrac [private]

Definition at line 101 of file MakeClusterSS.h. Referenced by Update().

TH2D* clusterss::MakeClusterSS::phHoughHist [private]

Definition at line 116 of file MakeClusterSS.h. Referenced by BestHough(), and MakeClusterSS().

std::map<int,recobase::Cluster *> clusterss::MakeClusterSS::plnclusters [private]

Definition at line 107 of file MakeClusterSS.h. Referenced by CleanUpContainers(), FindAllTransverseWindows(), FindLongWindow(), and MakePlaneClusters().

std::vector<clusterss::WindowSS *> clusterss::MakeClusterSS::plnwindows [private]

Definition at line 108 of file MakeClusterSS.h. Referenced by CleanUpContainers(), ClusterWindows(), FindAllTransverseWindows(), HoughTransCluster(), and StupidHoughNoWork().

Int_t clusterss::MakeClusterSS::pLongWindowGapPlaneCut [private]

Definition at line 98 of file MakeClusterSS.h. Referenced by FindLongWindow(), and Update().

Double_t clusterss::MakeClusterSS::pLongWindowMipCut [private]

Definition at line 97 of file MakeClusterSS.h. Referenced by Update().

Double_t clusterss::MakeClusterSS::pMaxAngleDiff [private]

Definition at line 94 of file MakeClusterSS.h. Referenced by Update().

int clusterss::MakeClusterSS::pMAXHOUGHITER [private]

Definition at line 102 of file MakeClusterSS.h. Referenced by Update().

Double_t clusterss::MakeClusterSS::pMaxOvlpLink [private]

Definition at line 90 of file MakeClusterSS.h. Referenced by Update().

Int_t clusterss::MakeClusterSS::pMaxPDiff [private]

Definition at line 91 of file MakeClusterSS.h. Referenced by Update().

Double_t clusterss::MakeClusterSS::pMaxTDiff [private]

Definition at line 93 of file MakeClusterSS.h. Referenced by Update().

Float_t clusterss::MakeClusterSS::pMaxWidToLenP2 [private]

Definition at line 95 of file MakeClusterSS.h. Referenced by Update().

Double_t clusterss::MakeClusterSS::pMaxZDiff [private]

Definition at line 92 of file MakeClusterSS.h. Referenced by Update().

Bool_t clusterss::MakeClusterSS::pMergeClusters [private]

Definition at line 89 of file MakeClusterSS.h. Referenced by Update().

Double_t clusterss::MakeClusterSS::pMinHoughPH [private]

Definition at line 100 of file MakeClusterSS.h. Referenced by Update().

Int_t clusterss::MakeClusterSS::pMinHoughPlanes [private]

Definition at line 99 of file MakeClusterSS.h. Referenced by HoughTransCluster(), and Update().

int clusterss::MakeClusterSS::printlevel [private]

Definition at line 74 of file MakeClusterSS.h. Referenced by CheckCurrSoln(), CleanUpWindows(), ClusterWindows(), DoCluster(), FindAllTransverseWindows(), FindLongWindow(), MakePlaneClusters(), PruneCellList(), and Update().

int clusterss::MakeClusterSS::run [private]

Definition at line 77 of file MakeClusterSS.h. Referenced by FindLongWindow(), PruneCellList(), and Reco().

int clusterss::MakeClusterSS::shwMaxPlane [private]

Definition at line 84 of file MakeClusterSS.h. Referenced by ClusterWindows(), FindAllTransverseWindows(), FindLongWindow(), and FindPlaneFromShwMax().

std::vector<clusterss::WindowSS *> clusterss::MakeClusterSS::windows [private]

Definition at line 109 of file MakeClusterSS.h. Referenced by CheckCurrSoln(), CleanUpContainers(), CleanUpWindows(), ClusterWindows(), DoCluster(), PruneCellList(), and StupidHoughNoWork().

std::vector<clusterss::CellHitGrad*> clusterss::MakeClusterSS::xchglist [private]

Definition at line 104 of file MakeClusterSS.h. Referenced by CleanUpAllContainers(), and SplitIntoViews().

std::vector<recobase::Cluster*> clusterss::MakeClusterSS::xclusterlist [private]

Definition at line 112 of file MakeClusterSS.h. Referenced by Ana(), CleanUpAllContainers(), and Reco().

std::vector<clusterss::CellHitGrad*> clusterss::MakeClusterSS::ychglist [private]

Definition at line 105 of file MakeClusterSS.h. Referenced by CleanUpAllContainers(), and SplitIntoViews().

std::vector<recobase::Cluster*> clusterss::MakeClusterSS::yclusterlist [private]

Definition at line 113 of file MakeClusterSS.h. Referenced by Ana(), CleanUpAllContainers(), and Reco().

---

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

• MakeClusterSS.h
• MakeClusterSS.cxx

---