AlgShowerSSList Class Reference

#include <AlgShowerSSList.h>

Inheritance diagram for AlgShowerSSList:

List of all members.

Public Member Functions
	AlgShowerSSList ()
virtual	~AlgShowerSSList ()
virtual void	RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void	Trace (const char *c) const
Private Attributes
TH1F *	vtxUHistAll
TH1F *	vtxUHist
TH1F *	vtxVHistAll
TH1F *	vtxVHist

---

Constructor & Destructor Documentation

AlgShowerSSList::AlgShowerSSList (   )

Definition at line 50 of file AlgShowerSSList.cxx. { vtxUHistAll = new TH1F("vtxUHistAll","vtxUHistAll",1,-4.,+4); vtxUHist = new TH1F("vtxUHist","vtxUHist",1,-4.,+4); vtxVHistAll = new TH1F("vtxVHistAll","vtxVHistAll",1,-4.,+4); vtxVHist = new TH1F("vtxVHist","vtxVHist",1,-4.,+4); }

AlgShowerSSList::~AlgShowerSSList (   )  [virtual]

Definition at line 61 of file AlgShowerSSList.cxx. { delete vtxUHistAll; delete vtxUHist; delete vtxVHistAll; delete vtxVHist; }

---

Member Function Documentation

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

We begin with a nested loop over all CandSlices, and over all CandClusters within each CandSlice, generating list of all clusters in the U and V planes. Next, we execute a nested loop over U-view and V-view CandClusters, checking for matching 2D clusters. The following criteria are used: Beginning planes for the two 2D 'long' CandClusters (defined by PlaneScale) must be within DiffViewPlaneMatch., or DiffViewPlaneMatchShort if one or more clusters is short. Beginning times for the two 2D CandClusters must be within DiffViewTimeMatch. The fractional overlap of U and V clusters must exceed DiffViewPlaneCoverage, if both clusters are long. The total pulse heights of the U and V clusters must agree to within DiffViewPulseHeightCut if at least one of the clusters is long. (Note: This should probably be changed to a charge ratio) If a U/V cluster set is found which satisfies these requirements, a second loop over U clusters is performed, searching for alternative U/V pairs involving the original V cluster satisfying the requirements above, but higher in total energy. If no better alternative is found a CandShower is constructed. In this way, a given CandCluster will only be used once, in the shower with the highest total energy. Reimplemented from AlgShowerSRList. Definition at line 73 of file AlgShowerSSList.cxx. References CandHandle::AddDaughterLink(), Registry::Get(), AlgFactory::GetAlgHandle(), CandSubShowerSRHandle::GetAveTime(), CandSubShowerSRHandle::GetAvgDev(), CandRecoHandle::GetBegPlane(), CandContext::GetCandRecord(), CandRecoHandle::GetCandSlice(), CandStripHandle::GetCharge(), CandSubShowerSRHandle::GetClusterID(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), CandRecoHandle::GetEndPlane(), CandShowerHandle::GetEnergy(), CandSubShowerSRHandle::GetEnergy(), AlgFactory::GetInstance(), CandSubShowerSRHandle::GetMaxStpTime(), CandSubShowerSRHandle::GetMaxU(), CandSubShowerSRHandle::GetMaxV(), CandSubShowerSRHandle::GetMinStpTime(), CandSubShowerSRHandle::GetMinU(), CandSubShowerSRHandle::GetMinV(), CandContext::GetMom(), CandSubShowerSRHandle::GetPlaneView(), CandSubShowerSRHandle::GetSlope(), CandStripHandle::GetTPos(), CandHandle::GetUidInt(), RecMinos::GetVldContext(), CandRecoHandle::GetVtxU(), CandRecoHandle::GetVtxV(), CandShowerSR::MakeCandidate(), MSG, PITCHINMETRES, CandHandle::Print(), CandRecoHandle::SetCandSlice(), vtxUHist, vtxUHistAll, vtxVHist, and vtxVHistAll. { assert(cx.GetDataIn()); if (!(cx.GetDataIn()->InheritsFrom("TObjArray"))) { return; } const CandSliceListHandle *slicelist = 0; const CandClusterListHandle *clusterlist = 0; const CandSubShowerSRListHandle *subshowerlist = 0; const TObjArray *cxin = dynamic_cast<const TObjArray *>(cx.GetDataIn()); for (Int_t i=0; i<=cxin->GetLast(); i++) { TObject *tobj = cxin->At(i); if (tobj->InheritsFrom("CandSliceListHandle")) { slicelist = dynamic_cast<CandSliceListHandle*>(tobj); } if (tobj->InheritsFrom("CandClusterListHandle")) { clusterlist = dynamic_cast<CandClusterListHandle*>(tobj); } if (tobj->InheritsFrom("CandSubShowerSRListHandle")) { subshowerlist = dynamic_cast<CandSubShowerSRListHandle*>(tobj); } } if (!slicelist || !clusterlist || !subshowerlist) { MSG("AlgShowerSS",Msg::kWarning) << "CandSliceListHandle, CandClusterListHandle or CandSubShowerSRListHandle missing\n"; return; } CandContext cxx(this,cx.GetMom()); Double_t MinNonPhysicsEnergy = 0.2; //in GeV Int_t MaxPlaneGap = 2; Int_t NPeakFindingBins = 32; Double_t SubShowerTimingCut = 50.0; //in ns Double_t EnergyAsymmetryCut = 0.3; Double_t PeakMergeTimeCut = 10.0; //in ns //get config for AlgSubShowerSS const char *charShowerSSAlgConfig = 0; ac.Get("ShowerSSAlgConfig",charShowerSSAlgConfig); ac.Get("MaxPlaneGap",MaxPlaneGap); ac.Get("NPeakFindingBins",NPeakFindingBins); ac.Get("MinNonPhysicsEnergy",MinNonPhysicsEnergy); ac.Get("SubShowerTimingCut",SubShowerTimingCut); ac.Get("EnergyAsymmetryCut",EnergyAsymmetryCut); ac.Get("PeakMergeTimeCut",PeakMergeTimeCut); //reset all hists and set bins // (add an extra bin at each end outside of expected range, // in order to be able to detect peaks at (far) detector edges) vtxUHistAll->Reset(); vtxUHistAll->SetBins(NPeakFindingBins+2, -4. - 8./Float_t(NPeakFindingBins), +4. + 8./Float_t(NPeakFindingBins)); vtxVHistAll->Reset(); vtxVHistAll->SetBins(NPeakFindingBins+2, -4. - 8./Float_t(NPeakFindingBins), +4. + 8./Float_t(NPeakFindingBins)); vtxUHist->Reset(); vtxUHist->SetBins(NPeakFindingBins+2, -4. - 8./Float_t(NPeakFindingBins), +4. + 8./Float_t(NPeakFindingBins)); vtxVHist->Reset(); vtxVHist->SetBins(NPeakFindingBins+2, -4. - 8./Float_t(NPeakFindingBins), +4. + 8./Float_t(NPeakFindingBins)); //Get singleton instance of AlgFactory AlgFactory &af = AlgFactory::GetInstance(); AlgHandle ah = af.GetAlgHandle("AlgShowerSS",charShowerSSAlgConfig); const CandRecord *candrec = cx.GetCandRecord(); assert(candrec); const VldContext *vldcptr = candrec->GetVldContext(); assert(vldcptr); VldContext vldc = *vldcptr; UgliGeomHandle ugh(vldc); CandSliceHandleItr sliceItr(slicelist->GetDaughterIterator()); Int_t nslice = 0; Int_t nsubshower = 0; while (CandSliceHandle *slice = sliceItr()) { //slice loop MSG("AlgShowerSS",Msg::kDebug) << "Slice: " << nslice << endl; // Select SubShowers within this slice CandSubShowerSRHandleItr subshowerItr(subshowerlist->GetDaughterIterator()); // iterate over subshowers, and make a shower out of // all subshowers in each slice, in contiguous planes int plane[500] = {0}; int firstPlane = 500; int lastPlane = 0; nsubshower = 0; while (CandSubShowerSRHandle *subshower = subshowerItr()) { if (*subshower->GetCandSlice()==*slice) { MSG("AlgShowerSS",Msg::kDebug) << "SubShower: " << nsubshower << " in slice " << nslice << endl; for(int i=subshower->GetBegPlane();i<=subshower->GetEndPlane();i++){ MSG("AlgShowerSS",Msg::kDebug) << "Plane: " << i << endl; plane[i] = 1; if(i<firstPlane) firstPlane = i; if(i>lastPlane) lastPlane = i; } } nsubshower++; } MSG("AlgShowerSS",Msg::kDebug) << "First Plane = " << firstPlane << endl; MSG("AlgShowerSS",Msg::kDebug) << "Last Plane = " << lastPlane << endl; //separate up into plane chunks according to MaxPlaneGap Int_t nshower = 0; while(firstPlane<=lastPlane){ //while loop over planes Int_t begPlane = 500; Int_t endPlane = 0; Int_t nGaps = 0; for(int i=firstPlane;i<=lastPlane;i++){ if(plane[i]==0) { if(begPlane!=500) nGaps+=1; } else { if(begPlane==500) begPlane = i; endPlane = i; } plane[i] = 0; if(nGaps>MaxPlaneGap) { MSG("AlgShowerSS",Msg::kDebug) << "begPlane = " << begPlane << endl; MSG("AlgShowerSS",Msg::kDebug) << "endPlane = " << endPlane << endl; break; } } firstPlane = endPlane+1; //could have multiple showers per plane chunk //pull out all subshowers in plane chunk and look for TPos peaks: //loop through subshower list again and histogram tposvtx //in order to look for peaks vtxUHistAll->Reset(); vtxVHistAll->Reset(); vtxUHist->Reset(); vtxVHist->Reset(); subshowerItr.ResetFirst(); nsubshower = 0; while (CandSubShowerSRHandle *subshower = subshowerItr()) { if (*subshower->GetCandSlice()==*slice) { if(subshower->GetBegPlane()>=begPlane && subshower->GetBegPlane()<=endPlane){ if(subshower->GetPlaneView()==PlaneView::kU || subshower->GetPlaneView()==PlaneView::kX) { MSG("AlgShowerSS",Msg::kDebug) << "Filling vtxUHist with Subshower " << nsubshower << " with vtvU = " << subshower->GetVtxU() << " and energy = " << subshower->GetEnergy() << " and slope = " << subshower->GetSlope() << " and AvgDev = " << subshower->GetAvgDev() << endl; CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); while (CandStripHandle *stp = stripssItr()) { //fill with max charge for this tpos Float_t val = stp->GetCharge(CalDigitType::kSigCorr) - vtxUHistAll->GetBinContent(vtxUHistAll->FindBin(stp->GetTPos())); if(val>0) vtxUHistAll->Fill(stp->GetTPos(),val); if( subshower->GetClusterID()==ClusterType::kEMLike || subshower->GetClusterID()==ClusterType::kHadLike || subshower->GetClusterID()==ClusterType::kTrkLike){ val = stp->GetCharge(CalDigitType::kSigCorr) - vtxUHist->GetBinContent(vtxUHist->FindBin(stp->GetTPos())); if(val>0) vtxUHist->Fill(stp->GetTPos(),val); } } } else if(subshower->GetPlaneView()==PlaneView::kV || subshower->GetPlaneView()==PlaneView::kY){ MSG("AlgShowerSS",Msg::kDebug) << "Filling vtxVHist with Subshower " << nsubshower << " with vtvV = " << subshower->GetVtxV() << " and energy = " << subshower->GetEnergy() << " and slope = " << subshower->GetSlope() << " and AvgDev = " << subshower->GetAvgDev() << endl; CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); while (CandStripHandle *stp = stripssItr()) { //fill with max charge for this tpos Float_t val = stp->GetCharge(CalDigitType::kSigCorr) - vtxVHistAll->GetBinContent(vtxVHistAll->FindBin(stp->GetTPos())); if(val>0) vtxVHistAll->Fill(stp->GetTPos(),val); if( subshower->GetClusterID()==ClusterType::kEMLike || subshower->GetClusterID()==ClusterType::kHadLike || subshower->GetClusterID()==ClusterType::kTrkLike){ val = stp->GetCharge(CalDigitType::kSigCorr) - vtxVHist->GetBinContent(vtxVHist->FindBin(stp->GetTPos())); if(val>0) vtxVHist->Fill(stp->GetTPos(),val); } } } } } nsubshower++; } if(false){ TCanvas *can = new TCanvas(); can->Divide(2,2); can->cd(1); vtxUHistAll->Draw(); can->cd(2); vtxVHistAll->Draw(); can->cd(3); vtxUHist->Draw(); can->cd(4); vtxVHist->Draw(); char nomus[256]; sprintf(nomus,"peakPlot_slice%i.eps",nslice); can->Print(nomus); delete can; } //width of shower in units of number of bins //assume shower typically has a width of ~4 strips Double_t sigma = (4.*0.0417/8.)*(NPeakFindingBins); //don't let width be <=1 bin to prevent //against clipping error messages if(sigma<=1.0) sigma = 1.01; Int_t nUPeaks = 0; Int_t nUPeaks_save = 0; Int_t nUPeaksAll = 0; Float_t peakUPos[10] = {0}; Float_t peakUVal[10] = {0}; Float_t tmp_peakUPos = vtxUHist->GetBinCenter(vtxUHist->GetMaximumBin()); Float_t tmp_peakUVal = vtxUHist->GetMaximum(); MSG("AlgShowerSS",Msg::kDebug) << "U Hist Integrals: " << vtxUHist->Integral() << " " << vtxUHistAll->Integral() << endl; for(int bins = 1;bins<=vtxUHist->GetNbinsX();bins++){ if(vtxUHist->GetBinContent(bins)>0) nUPeaks++; if(vtxUHistAll->GetBinContent(bins)>0) nUPeaksAll++; } if(nUPeaks==1) { peakUPos[0] = tmp_peakUPos; peakUVal[0] = tmp_peakUVal; } else if(nUPeaks>0) { TSpectrum spec(10); Float_t *source = new Float_t[vtxUHist->GetNbinsX()]; Float_t *dest = new Float_t[vtxUHist->GetNbinsX()]; for(int h_loop=0;h_loop<vtxUHist->GetNbinsX();h_loop++){ source[h_loop] = vtxUHist->GetBinContent(h_loop+1); dest[h_loop] = 0; } #if ROOT_VERSION_CODE >= ROOT_VERSION(5,9,1) nUPeaks = spec.SearchHighRes(source,dest, vtxUHist->GetNbinsX(), sigma,5,0,1,0,1); #else nUPeaks = spec.Search1HighRes(source,dest, vtxUHist->GetNbinsX(), sigma,5,0,1,0,1); #endif Float_t *posX = spec.GetPositionX(); Float_t *posX_copy = new Float_t[nUPeaks]; for(int p_loop=0;p_loop<nUPeaks;p_loop++) posX_copy[p_loop] = posX[p_loop]; //make sure peaks are in ascending order for(int p_loop=0;p_loop<nUPeaks-1;p_loop++) { if(posX_copy[p_loop+1]<posX_copy[p_loop]) { Float_t temp = posX_copy[p_loop]; posX_copy[p_loop] = posX_copy[p_loop+1]; posX_copy[p_loop+1] = temp; p_loop = -1; } } for(int p_loop=0;p_loop<nUPeaks;p_loop++){ Int_t bin = 1+Int_t(posX_copy[p_loop]+0.5); peakUPos[p_loop] = vtxUHist->GetBinCenter(bin); peakUVal[p_loop] = vtxUHist->GetBinContent(bin); } delete [] posX_copy; delete [] source; delete [] dest; } if(nUPeaksAll>1) { TSpectrum spec(10); Float_t *source = new Float_t[vtxUHistAll->GetNbinsX()]; Float_t *dest = new Float_t[vtxUHistAll->GetNbinsX()]; for(int h_loop=0;h_loop<vtxUHistAll->GetNbinsX();h_loop++){ source[h_loop] = vtxUHistAll->GetBinContent(h_loop+1); dest[h_loop] = 0; } #if ROOT_VERSION_CODE >= ROOT_VERSION(5,9,1) nUPeaksAll = spec.SearchHighRes(source,dest, vtxUHistAll->GetNbinsX(), sigma,5,0,1,0,1); #else nUPeaksAll = spec.Search1HighRes(source,dest, vtxUHistAll->GetNbinsX(), sigma,5,0,1,0,1); #endif delete [] source; delete [] dest; } nUPeaks_save = nUPeaks; if(nUPeaksAll==0 && nUPeaks>0) nUPeaksAll = nUPeaks; Int_t nVPeaks = 0; Int_t nVPeaks_save = 0; Int_t nVPeaksAll = 0; Float_t peakVPos[10] = {0}; Float_t peakVVal[10] = {0}; Float_t tmp_peakVPos = vtxVHist->GetBinCenter(vtxVHist->GetMaximumBin()); Float_t tmp_peakVVal = vtxVHist->GetMaximum(); MSG("AlgShowerSS",Msg::kDebug) << "V Hist Integrals: " << vtxVHist->Integral() << " " << vtxVHistAll->Integral() << endl; for(int bins = 1;bins<=vtxVHist->GetNbinsX();bins++){ if(vtxVHist->GetBinContent(bins)>0) nVPeaks++; if(vtxVHistAll->GetBinContent(bins)>0) nVPeaksAll++; } if(nVPeaks==1) { peakVPos[0] = tmp_peakVPos; peakVVal[0] = tmp_peakVVal; } else if(nVPeaks>0) { TSpectrum spec(10); Float_t *source = new Float_t[vtxVHist->GetNbinsX()]; Float_t *dest = new Float_t[vtxVHist->GetNbinsX()]; for(int h_loop=0;h_loop<vtxVHist->GetNbinsX();h_loop++){ source[h_loop] = vtxVHist->GetBinContent(h_loop+1); dest[h_loop] = 0; } #if ROOT_VERSION_CODE >= ROOT_VERSION(5,9,1) nVPeaks = spec.SearchHighRes(source,dest, vtxVHist->GetNbinsX(), sigma,5,0,1,0,1); #else nVPeaks = spec.Search1HighRes(source,dest, vtxVHist->GetNbinsX(), sigma,5,0,1,0,1); #endif Float_t *posX = spec.GetPositionX(); Float_t *posX_copy = new Float_t[nVPeaks]; for(int p_loop=0;p_loop<nVPeaks;p_loop++) posX_copy[p_loop] = posX[p_loop]; //make sure peaks are in ascending order for(int p_loop=0;p_loop<nVPeaks-1;p_loop++) { if(posX_copy[p_loop+1]<posX_copy[p_loop]) { Float_t temp = posX_copy[p_loop]; posX_copy[p_loop] = posX_copy[p_loop+1]; posX_copy[p_loop+1] = temp; p_loop = -1; } } for(int p_loop=0;p_loop<nVPeaks;p_loop++){ Int_t bin = 1+Int_t(posX_copy[p_loop]+0.5); peakVPos[p_loop] = vtxVHist->GetBinCenter(bin); peakVVal[p_loop] = vtxVHist->GetBinContent(bin); } delete [] posX_copy; delete [] source; delete [] dest; } if(nVPeaksAll>1) { TSpectrum spec(10); Float_t *source = new Float_t[vtxVHistAll->GetNbinsX()]; Float_t *dest = new Float_t[vtxVHistAll->GetNbinsX()]; for(int h_loop=0;h_loop<vtxVHistAll->GetNbinsX();h_loop++){ source[h_loop] = vtxVHistAll->GetBinContent(h_loop+1); dest[h_loop] = 0; } #if ROOT_VERSION_CODE >= ROOT_VERSION(5,9,1) nVPeaksAll = spec.SearchHighRes(source,dest, vtxVHistAll->GetNbinsX(), sigma,5,0,1,0,1); #else nVPeaksAll = spec.Search1HighRes(source,dest, vtxVHistAll->GetNbinsX(), sigma,5,0,1,0,1); #endif delete [] source; delete [] dest; } nVPeaks_save = nVPeaks; if(nVPeaksAll==0 && nVPeaks>0) nVPeaksAll = nVPeaks; //if no peaks in either view if(nUPeaksAll==0 || nVPeaksAll==0) { //can't form a 3D shower in these conditions! MSG("AlgShowerSS",Msg::kDebug) << "nUPeaks = " << nUPeaks << " nVPeaks = " << nVPeaks << " Can't form a 3D shower" << endl; continue; } //one or both views have no good "physics" subshower: if(nUPeaks==0 || nVPeaks==0){ MSG("AlgShowerSS",Msg::kDebug) << "Have something in both views but " << "nUPeaks = " << nUPeaks << "and nVPeaks = " << nVPeaks << endl; //try to form anyway with whatever's there if(nUPeaks==0) nUPeaks = 1; if(nVPeaks==0) nVPeaks = 1; } //Try to form some 3D showers, yeah! //have one or more peak => maybe one or //more shower in this plane range //get boundaries that define each shower in U view: std::vector<Float_t> peakUPosLow(nUPeaks,0); std::vector<Float_t> peakUPosUpp(nUPeaks,0); MSG("AlgShowerSS",Msg::kDebug) << "Number of U Peaks found = " << nUPeaks << " with positions, heights " << "and bounds: " << endl; for (int i=0;i<nUPeaks;i++){ if(i==0) peakUPosLow[i] = -5; else peakUPosLow[i] = peakUPosUpp[i-1]; if(i==nUPeaks-1) peakUPosUpp[i] = 5; else { Int_t bin = vtxUHist->GetXaxis()->FindBin(peakUPos[i]); Int_t bestBin = vtxUHist->GetNbinsX(); Float_t threeBinAve = 1000000; //sigcor! while(bin<vtxUHist->GetNbinsX() && vtxUHist->GetBinCenter(bin)<peakUPos[i+1]) { Float_t ave = (vtxUHist->GetBinContent(bin-1) + vtxUHist->GetBinContent(bin) + vtxUHist->GetBinContent(bin+1)); if(ave<threeBinAve) { threeBinAve = ave; bestBin = bin; } bin++; } peakUPosUpp[i] = vtxUHist->GetBinCenter(bestBin); } MSG("AlgShowerSS",Msg::kDebug) << peakUPos[i] << " " << peakUVal[i] << " [" << peakUPosLow[i] << "," << peakUPosUpp[i] << "]" << endl; } //get boundaries that define each shower in V view: std::vector<Float_t> peakVPosLow(nVPeaks,0); std::vector<Float_t> peakVPosUpp(nVPeaks,0); MSG("AlgShowerSS",Msg::kDebug) << "Number of V Peaks found = " << nVPeaks << " with positions, heights " << "and bounds: " << endl; for (int i=0;i<nVPeaks;i++){ if(i==0) peakVPosLow[i] = -5; else peakVPosLow[i] = peakVPosUpp[i-1]; if(i==nVPeaks-1) peakVPosUpp[i] = 5; else { Int_t bin = vtxVHist->GetXaxis()->FindBin(peakVPos[i]); Int_t bestBin = vtxVHist->GetNbinsX(); Float_t threeBinAve = 1000000; //gev! while(bin<vtxVHist->GetNbinsX() && vtxVHist->GetBinCenter(bin)<peakVPos[i+1]) { Float_t ave = (vtxVHist->GetBinContent(bin-1) + vtxVHist->GetBinContent(bin) + vtxVHist->GetBinContent(bin+1)); if(ave<threeBinAve) { threeBinAve = ave; bestBin = bin; } bin++; } peakVPosUpp[i] = vtxVHist->GetBinCenter(bestBin); } MSG("AlgShowerSS",Msg::kDebug) << peakVPos[i] << " " << peakVVal[i] << " [" << peakVPosLow[i] << "," << peakVPosUpp[i] << "]" << endl; } //for matching up U/V views of showers: std::vector<Float_t> totUEnergy(nUPeaks,0); std::vector<Double_t> aveUTime(nUPeaks,0); std::vector<Float_t> approxVfromTime(nUPeaks,0); std::vector<Float_t> maxPhysU(nUPeaks,0); std::vector<Float_t> minPhysU(nUPeaks,0); std::vector<Int_t> maxPlaneU(nUPeaks,0); std::vector<Int_t> minPlaneU(nUPeaks,0); std::vector<Int_t> nSubShowersU(nUPeaks,0); for(int i=0;i<nUPeaks;i++) { totUEnergy[i] = 0; aveUTime[i] = 0; approxVfromTime[i] = 0; maxPhysU[i] = -1; maxPlaneU[i] = -1; minPhysU[i] = 999; minPlaneU[i] = -1; nSubShowersU[i] = 0; } std::vector<Float_t> totVEnergy(nVPeaks,0); std::vector<Double_t> aveVTime(nVPeaks,0); std::vector<Float_t> approxUfromTime(nVPeaks,0); std::vector<Float_t> maxPhysV(nVPeaks,0); std::vector<Float_t> minPhysV(nVPeaks,0); std::vector<Int_t> maxPlaneV(nVPeaks,0); std::vector<Int_t> minPlaneV(nVPeaks,0); std::vector<Int_t> nSubShowersV(nVPeaks,0); for(int i=0;i<nVPeaks;i++) { totVEnergy[i] = 0; aveVTime[i] = 0; approxUfromTime[i] = 0; maxPhysV[i] = -1; maxPlaneV[i] = -1; minPhysV[i] = 999; minPlaneV[i] = -1; nSubShowersV[i] = 0; } //loop through subshowerlist and add up energy of showers subshowerItr.ResetFirst(); while (CandSubShowerSRHandle *subshower = subshowerItr()) { if (*subshower->GetCandSlice()==*slice) { if(subshower->GetBegPlane()>=begPlane && subshower->GetEndPlane()<=endPlane) { if(subshower->GetPlaneView()==PlaneView::kU || subshower->GetPlaneView()==PlaneView::kX){ for(int i=0;i<nUPeaks;i++){ Float_t fracPHInRange = 0; Float_t totPHInSS = 0; CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); while (CandStripHandle *stp = stripssItr()) { totPHInSS += stp->GetCharge(CalDigitType::kSigCorr); if(stp->GetTPos()>peakUPosLow[i] && stp->GetTPos()<=peakUPosUpp[i]) fracPHInRange+=stp->GetCharge(CalDigitType::kSigCorr); } if(totPHInSS>=0 && fracPHInRange/totPHInSS >= 0.5) { totUEnergy[i] += subshower->GetEnergy(); aveUTime[i] += 1e9*subshower->GetAveTime()*subshower->GetEnergy(); if(nUPeaks>1 && subshower->GetClusterID()==ClusterType::kEMLike || subshower->GetClusterID()==ClusterType::kHadLike || subshower->GetClusterID()==ClusterType::kTrkLike){ approxVfromTime[i] += subshower->GetVtxV(); nSubShowersU[i] += 1; for(int j=subshower->GetBegPlane();j<=subshower->GetEndPlane();j++){ Float_t testMinU = subshower->GetMinU(j); Float_t testMaxU = subshower->GetMaxU(j); if(testMaxU>maxPhysU[i]) { maxPhysU[i] = testMaxU; maxPlaneU[i] = j; } if(testMinU<minPhysU[i]) { minPhysU[i] = testMinU; minPlaneU[i] = j; } } } } } } else if(subshower->GetPlaneView()==PlaneView::kV || subshower->GetPlaneView()==PlaneView::kY){ for(int i=0;i<nVPeaks;i++){ Float_t fracPHInRange = 0; Float_t totPHInSS = 0; CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); while (CandStripHandle *stp = stripssItr()) { totPHInSS += stp->GetCharge(CalDigitType::kSigCorr); if(stp->GetTPos()>peakVPosLow[i] && stp->GetTPos()<=peakVPosUpp[i]) fracPHInRange+=stp->GetCharge(CalDigitType::kSigCorr); } if(totPHInSS>=0 && fracPHInRange/totPHInSS >= 0.5) { totVEnergy[i] += subshower->GetEnergy(); aveVTime[i] += 1e9*subshower->GetAveTime()*subshower->GetEnergy(); if(nVPeaks>1 && subshower->GetClusterID()==ClusterType::kEMLike || subshower->GetClusterID()==ClusterType::kHadLike || subshower->GetClusterID()==ClusterType::kTrkLike){ approxUfromTime[i] += subshower->GetVtxU(); nSubShowersV[i] += 1; for(int j=subshower->GetBegPlane();j<=subshower->GetEndPlane();j++){ Float_t testMinV = subshower->GetMinV(j); Float_t testMaxV = subshower->GetMaxV(j); if(testMaxV>maxPhysV[i]) { maxPhysV[i] = testMaxV; maxPlaneV[i] = j; } if(testMinV<minPhysV[i]) { minPhysV[i] = testMinV; minPlaneV[i] = j; } } } } } } } } } //define number of showers there will be in this plane range Int_t nShowers = nUPeaks; if(nShowers>nVPeaks) nShowers = nVPeaks; MSG("AlgShowerSS",Msg::kDebug) << "nShowers = " << nShowers << endl; //Combine peaks if there is a U/V mismatch in number of peaks //U planes: Int_t sanity_counter = 0; //just in case it decides to keep on looping... Int_t nUPeaks_tmp = nUPeaks; while(nUPeaks_tmp!=nShowers && sanity_counter<3){ Int_t merge_peak_1 = -1; //lowest merge peak Int_t merge_peak_2 = -1; //highest merge peak Float_t best_merge_metric = 999999; for(int i=0;i<nUPeaks-1;i++){ //don't care about empty peak regions: if(totUEnergy[i]<=0) continue; MSG("AlgShowerSS",Msg::kDebug) << "Peak " << i << " is not empty" << endl; //find next peak region that isn't empty: Int_t comp_index = 1; for (; (i+comp_index)<nUPeaks; ++comp_index) if (totUEnergy[i+comp_index]>0) break; /* while( totUEnergy[i+comp_index]<=0 && i+comp_index<nUPeaks) { comp_index += 1; } */ if(i+comp_index>=nUPeaks) break; //no more comparisons available MSG("AlgShowerSS",Msg::kDebug) << "Compare with peak " << i+comp_index << endl; Float_t merge_metric = 999999; Float_t time_diff = 999999; if(totUEnergy[i+comp_index]>0 && minPlaneU[i+comp_index]!=-1) { //get distance of closest approach in metres: merge_metric = TMath::Power(minPhysU[i+comp_index] - maxPhysU[i],2); merge_metric += PITCHINMETRES*TMath::Power(minPlaneU[i+comp_index] - maxPlaneU[i],2); merge_metric = TMath::Sqrt(merge_metric); MSG("AlgShowerSS",Msg::kDebug) << "Distance metric = " << merge_metric << endl; //add in effective distance based on timing: merge_metric += Mphysical::c_light*1e-9*TMath::Abs(aveUTime[i+comp_index] / totUEnergy[i+comp_index] - aveUTime[i]/totUEnergy[i]); MSG("AlgShowerSS",Msg::kDebug) << "Time metric = " << Mphysical::c_light*1e-9* TMath::Abs(aveUTime[i+comp_index] / totUEnergy[i+comp_index] - aveUTime[i] / totUEnergy[i]) << endl; time_diff = ( TMath::Abs(aveUTime[i+comp_index] / totUEnergy[i+comp_index] - aveUTime[i] / totUEnergy[i] ) ); } else if(totUEnergy[i+comp_index]>0) { // just check timing since this is a non-physics subshower cluster // this will lead to a smaller merge_metric value in general // but it is probably better to merge the non-physics clusters in any case merge_metric = Mphysical::c_light*1e-9*TMath::Abs(aveUTime[i+comp_index] / totUEnergy[i+comp_index] - aveUTime[i]/totUEnergy[i]); time_diff = ( TMath::Abs(aveUTime[i+comp_index] / totUEnergy[i+comp_index] - aveUTime[i] / totUEnergy[i] ) ); } MSG("AlgShowerSS",Msg::kDebug) << "Final metric = " << merge_metric << endl; if(merge_metric<best_merge_metric && time_diff<PeakMergeTimeCut){ best_merge_metric = merge_metric; merge_peak_1 = i; merge_peak_2 = i+comp_index; MSG("AlgShowerSS",Msg::kDebug) << "New best merge peak!\n" << "Merge peak " << merge_peak_1 << " with " << merge_peak_2 << " based on metric value " << best_merge_metric << endl; } } if(best_merge_metric<999999 && merge_peak_1!=-1 && merge_peak_2!=-1){ //add energy to second peak and //remove energy from first peak totUEnergy[merge_peak_2] += totUEnergy[merge_peak_1]; totUEnergy[merge_peak_1] = 0; aveUTime[merge_peak_2] += aveUTime[merge_peak_1]; aveUTime[merge_peak_1] = 0; //add other-view vertex estimate sum approxVfromTime[merge_peak_2] += approxVfromTime[merge_peak_1]; approxVfromTime[merge_peak_1] = 0; if(nSubShowersU[merge_peak_1]>=0 && nSubShowersU[merge_peak_2]>=0) { //shift min U of second peak down and //set min U plane of second peak minPhysU[merge_peak_2] = minPhysU[merge_peak_1]; minPlaneU[merge_peak_2] = minPlaneU[merge_peak_1]; //set min/max limits of first peak to be the same maxPhysU[merge_peak_1] = minPhysU[merge_peak_1]; maxPlaneU[merge_peak_1] = minPlaneU[merge_peak_1]; //set peak lower limit of second peak to that of first peakUPosLow[merge_peak_2] = peakUPosLow[merge_peak_1]; //set peak lower and upper limit of first peak to be the same peakUPosUpp[merge_peak_1] = peakUPosLow[merge_peak_1]; } else if(nSubShowersU[merge_peak_1]==-1) { //first peak is non-physical so don't change anything for merge_peak_2 //set min/max limits of first peak to be the same maxPhysU[merge_peak_1] = minPhysU[merge_peak_1]; maxPlaneU[merge_peak_1] = minPlaneU[merge_peak_1]; //set peak lower and upper limit of first peak to be the same peakUPosUpp[merge_peak_1] = peakUPosLow[merge_peak_1]; } else if(nSubShowersU[merge_peak_2]==-1) { //set min/max limits of second peak to be the same as the first peak minPhysU[merge_peak_2] = minPhysU[merge_peak_1]; minPlaneU[merge_peak_2] = minPlaneU[merge_peak_1]; maxPhysU[merge_peak_2] = maxPhysU[merge_peak_1]; maxPlaneU[merge_peak_2] = maxPlaneU[merge_peak_1]; //set min/max limits of first peak to be the same maxPhysU[merge_peak_1] = minPhysU[merge_peak_1]; maxPlaneU[merge_peak_1] = minPlaneU[merge_peak_1]; //set peak lower limit of second peak to that of first peakUPosLow[merge_peak_2] = peakUPosLow[merge_peak_1]; //set peak lower and upper limit of first peak to be the same peakUPosUpp[merge_peak_1] = peakUPosLow[merge_peak_1]; } nSubShowersU[merge_peak_2] += nSubShowersU[merge_peak_1]; nSubShowersU[merge_peak_1] = 0; nUPeaks_tmp -= 1; } else { MSG("AlgShowerSS",Msg::kDebug) << "Incrementing U view sanity_counter:\n" << " nUPeaks_tmp = " << nUPeaks_tmp << " best_merge_metric = " << best_merge_metric << " merge_peak_1 = " << merge_peak_1 << " merge_peak_2 = " << merge_peak_2 << endl; sanity_counter+=1; } } //V planes: Int_t nVPeaks_tmp = nVPeaks; while(nVPeaks_tmp!=nShowers && sanity_counter<5){ Int_t merge_peak_1 = -1; //lowest merge peak Int_t merge_peak_2 = -1; //highest merge peak Float_t best_merge_metric = 999999; for(int i=0;i<nVPeaks-1;i++){ //don't care about empty peak regions: if(totVEnergy[i]<=0) continue; MSG("AlgShowerSS",Msg::kDebug) << "Peak " << i << " is not empty" << endl; //find next peak region that isn't empty: Int_t comp_index = 1; for (; (i+comp_index)<nVPeaks; ++comp_index) if (totVEnergy[i+comp_index]>0) break; /* while( totVEnergy[i+comp_index]<=0 && i+comp_index<nVPeaks) { comp_index += 1; } */ if(i+comp_index>=nVPeaks) break; //no more comparisons available MSG("AlgShowerSS",Msg::kDebug) << "Compare with peak " << i+comp_index << endl; Float_t merge_metric = 999999; if(totVEnergy[i+comp_index]>0 && minPlaneV[i+comp_index]!=-1) { //get distance of closest approach in metres: merge_metric = TMath::Power(minPhysV[i+comp_index] - maxPhysV[i],2); merge_metric += PITCHINMETRES*TMath::Power(minPlaneV[i+comp_index] - maxPlaneV[i],2); merge_metric = TMath::Sqrt(merge_metric); MSG("AlgShowerSS",Msg::kDebug) << "Distance metric = " << merge_metric << endl; //add in effective distance based on timing: merge_metric += Mphysical::c_light*1e-9*TMath::Abs(aveVTime[i+comp_index] / totVEnergy[i+comp_index] - aveVTime[i]/totVEnergy[i]); MSG("AlgShowerSS",Msg::kDebug) << "Time metric = " << Mphysical::c_light*1e-9* TMath::Abs(aveVTime[i+comp_index] / totVEnergy[i+comp_index] - aveVTime[i] / totVEnergy[i]) << endl; } else if(totVEnergy[i+comp_index]>0) { // just check timing since this is a non-physics subshower cluster // this will lead to a smaller merge_metric value in general // but it is probably better to merge the non-physics clusters in any case merge_metric = Mphysical::c_light*1e-9*TMath::Abs(aveVTime[i+comp_index] / totVEnergy[i+comp_index] - aveVTime[i]/totVEnergy[i]); } MSG("AlgShowerSS",Msg::kDebug) << "Final metric = " << merge_metric << endl; if(merge_metric<best_merge_metric){ best_merge_metric = merge_metric; merge_peak_1 = i; merge_peak_2 = i+comp_index; MSG("AlgShowerSS",Msg::kDebug) << "New best merge peak!\n" << "Merge peak " << merge_peak_1 << " with " << merge_peak_2 << " based on metric value " << best_merge_metric << endl; } } if(best_merge_metric<999999 && merge_peak_1!=-1 && merge_peak_2!=-1){ //add energy to second peak and //remove energy from first peak totVEnergy[merge_peak_2] += totVEnergy[merge_peak_1]; totVEnergy[merge_peak_1] = 0; aveVTime[merge_peak_2] += aveVTime[merge_peak_1]; aveVTime[merge_peak_1] = 0; //add other-view vertex estimate sum approxUfromTime[merge_peak_2] += approxUfromTime[merge_peak_1]; approxUfromTime[merge_peak_1] = 0; if(nSubShowersV[merge_peak_1]>=0 && nSubShowersV[merge_peak_2]>=0) { //shift min V of second peak down and //set min V plane of second peak minPhysV[merge_peak_2] = minPhysV[merge_peak_1]; minPlaneV[merge_peak_2] = minPlaneV[merge_peak_1]; //set min/max limits of first peak to be the same maxPhysV[merge_peak_1] = minPhysV[merge_peak_1]; maxPlaneV[merge_peak_1] = minPlaneV[merge_peak_1]; //set peak lower limit of second peak to that of first peakVPosLow[merge_peak_2] = peakVPosLow[merge_peak_1]; //set peak lower and upper limit of first peak to be the same peakVPosUpp[merge_peak_1] = peakVPosLow[merge_peak_1]; } else if(nSubShowersV[merge_peak_1]==-1) { //first peak is non-physical so don't change anything for merge_peak_2 //set min/max limits of first peak to be the same maxPhysV[merge_peak_1] = minPhysV[merge_peak_1]; maxPlaneV[merge_peak_1] = minPlaneV[merge_peak_1]; //set peak lower and upper limit of first peak to be the same peakVPosUpp[merge_peak_1] = peakVPosLow[merge_peak_1]; } else if(nSubShowersV[merge_peak_2]==-1) { //set min/max limits of second peak to be the same as the first peak minPhysV[merge_peak_2] = minPhysV[merge_peak_1]; minPlaneV[merge_peak_2] = minPlaneV[merge_peak_1]; maxPhysV[merge_peak_2] = maxPhysV[merge_peak_1]; maxPlaneV[merge_peak_2] = maxPlaneV[merge_peak_1]; //set min/max limits of first peak to be the same maxPhysV[merge_peak_1] = minPhysV[merge_peak_1]; maxPlaneV[merge_peak_1] = minPlaneV[merge_peak_1]; //set peak lower limit of second peak to that of first peakVPosLow[merge_peak_2] = peakVPosLow[merge_peak_1]; //set peak lower and upper limit of first peak to be the same peakVPosUpp[merge_peak_1] = peakVPosLow[merge_peak_1]; } nSubShowersV[merge_peak_2] += nSubShowersV[merge_peak_1]; nSubShowersV[merge_peak_1] = 0; nVPeaks_tmp -= 1; } else { MSG("AlgShowerSS",Msg::kDebug) << "Incrementing V view sanity_counter:\n" << " nVPeaks_tmp = " << nVPeaks_tmp << " best_merge_metric = " << best_merge_metric << " merge_peak_1 = " << merge_peak_1 << " merge_peak_2 = " << merge_peak_2 << endl; sanity_counter+=1; } } for(int i=0;i<nShowers;i++){ MSG("AlgShowerSS",Msg::kDebug) << "Positions, peaks and limits after merging:" << endl; MSG("AlgShowerSS",Msg::kDebug) << peakUPos[i] << " " << peakUVal[i] << " [" << peakUPosLow[i] << "," << peakUPosUpp[i] << "]" << endl; MSG("AlgShowerSS",Msg::kDebug) << peakVPos[i] << " " << peakVVal[i] << " [" << peakVPosLow[i] << "," << peakVPosUpp[i] << "]" << endl; } //get best order for matching subshower groups in U/V std::vector<Int_t> bestOrderU(nUPeaks,0); std::vector<Int_t> bestOrderV(nVPeaks,0); for(int i=0;i<nUPeaks;i++){ if(totUEnergy[i]>0){ Int_t nbigger = 0; for(int j=0;j<nUPeaks;j++){ if(totUEnergy[i]<totUEnergy[j]) nbigger+=1; } bestOrderU[nbigger] = i; } else { // if there is no energy in this "peak" Int_t nsmaller = 0; // put zero peaks in reverse order for(int j=0;j<i;j++){ // purely to ensure that all indices if(totUEnergy[j]==0) nsmaller+=1; // of bestOrderX are filled } // with legitimate values (prevents bestOrderU[nUPeaks-1-nsmaller] = i; // against a rare case seg fault) } } for(int i=0;i<nVPeaks;i++){ if(totVEnergy[i]>0){ Int_t nbigger = 0; for(int j=0;j<nVPeaks;j++){ if(totVEnergy[i]<totVEnergy[j]) nbigger+=1; } bestOrderV[nbigger] = i; } else { Int_t nsmaller = 0; for(int j=0;j<i;j++){ if(totVEnergy[j]==0) nsmaller+=1; } bestOrderV[nVPeaks-1-nsmaller] = i; } } //check energy ordering is ok: for(int i=0;i<nShowers-1;i++){ //loop over showers Float_t EUi = totUEnergy[bestOrderU[i]]; Float_t EVi = totVEnergy[bestOrderV[i]]; Float_t EUii = totUEnergy[bestOrderU[i+1]]; Float_t EVii = totVEnergy[bestOrderV[i+1]]; if(EUi==0 || EUii==0 || EVi==0 || EVii==0) continue; //can we swap the order and still get two good energy matches? Float_t asym = TMath::Abs(EUi - EVi)/(EUi + EVi) + TMath::Abs(EUii - EVii)/(EUii + EVii); Float_t asymSwap = TMath::Abs(EUi - EVii)/(EUi + EVii) + TMath::Abs(EUii - EVi)/(EUii + EVi); //if difference between two asymmetry estimates is <30% //see if vertex estimate from timing //or if plane range information helps if( TMath::Abs(asym - asymSwap)/asym < EnergyAsymmetryCut || TMath::Abs(aveUTime[bestOrderU[i]]/EUi - aveVTime[bestOrderV[i]]/EVi) > SubShowerTimingCut ) { if(nSubShowersU[bestOrderU[i]]==0 || nSubShowersV[bestOrderV[i]]==0 || nSubShowersU[bestOrderU[i+1]]==0 || nSubShowersV[bestOrderV[i+1]]==0) continue; Double_t timediff = TMath::Abs( (aveUTime[bestOrderU[i]]/EUi) - (aveVTime[bestOrderV[i]]/EVi) ) + TMath::Abs( (aveUTime[bestOrderU[i+1]]/EUii) - (aveVTime[bestOrderV[i+1]]/EVii) ); Double_t timediff_swap = TMath::Abs( (aveUTime[bestOrderU[i]]/EUi) - (aveVTime[bestOrderV[i+1]]/EVii) ) + TMath::Abs( (aveUTime[bestOrderU[i+1]]/EUii) - (aveVTime[bestOrderV[i]]/EVi) ); if(timediff_swap<timediff && 1.0e9*timediff_swap<SubShowerTimingCut) { //check that timing diff is reasonable //swap! //arbitrary choice which view to switch: Int_t tempIndex = bestOrderU[i]; bestOrderU[i] = bestOrderU[i+1]; bestOrderU[i+1] = tempIndex; //start checks again so that clusters //can propagate if better matches are found i=0; continue; } //now check plane range agreement: Int_t planediff = TMath::Abs(minPlaneU[bestOrderU[i]] - minPlaneV[bestOrderV[i]]) + TMath::Abs(maxPlaneU[bestOrderU[i]] - maxPlaneV[bestOrderV[i]]) + TMath::Abs(minPlaneU[bestOrderU[i+1]] - minPlaneV[bestOrderV[i+1]]) + TMath::Abs(maxPlaneU[bestOrderU[i+1]] - maxPlaneV[bestOrderV[i+1]]); Int_t planediff_swap = TMath::Abs(minPlaneU[bestOrderU[i]] - minPlaneV[bestOrderV[i+1]]) + TMath::Abs(maxPlaneU[bestOrderU[i]] - maxPlaneV[bestOrderV[i+1]]) + TMath::Abs(minPlaneU[bestOrderU[i+1]] - minPlaneV[bestOrderV[i]]) + TMath::Abs(maxPlaneU[bestOrderU[i+1]] - maxPlaneV[bestOrderV[i]]); if(planediff_swap<planediff && 1.0e9*timediff_swap<SubShowerTimingCut) { //swap! //arbitrary choice which view to switch: Int_t tempIndex = bestOrderU[i]; bestOrderU[i] = bestOrderU[i+1]; bestOrderU[i+1] = tempIndex; //start checks again so that clusters //can propagate if better matches are found i=0; continue; } if(false){ //space + time - not implemented yet approxVfromTime[bestOrderU[i]] /= nSubShowersU[bestOrderU[i]]; approxUfromTime[bestOrderV[i]] /= nSubShowersV[bestOrderV[i]]; approxVfromTime[bestOrderU[i+1]] /= nSubShowersU[bestOrderU[i+1]]; approxUfromTime[bestOrderV[i+1]] /= nSubShowersV[bestOrderV[i+1]]; //sum distance from both showers Float_t vertexDistance = TMath::Sqrt(TMath::Power(approxUfromTime[bestOrderV[i]] - peakUPos[bestOrderU[i]],2.) + TMath::Power(approxVfromTime[bestOrderU[i]] - peakVPos[bestOrderV[i]],2)) + TMath::Sqrt(TMath::Power(approxUfromTime[bestOrderV[i+1]] - peakUPos[bestOrderU[i+1]],2.) + TMath::Power(approxVfromTime[bestOrderU[i+1]] - peakVPos[bestOrderV[i+1]],2)); Float_t vertexDistanceSwap = TMath::Sqrt(TMath::Power(approxUfromTime[bestOrderV[i+1]] - peakUPos[bestOrderU[i]],2.) + TMath::Power(approxVfromTime[bestOrderU[i]] - peakVPos[bestOrderV[i+1]],2)) + TMath::Sqrt(TMath::Power(approxUfromTime[bestOrderV[i]] - peakUPos[bestOrderU[i+1]],2.) + TMath::Power(approxVfromTime[bestOrderU[i+1]] - peakVPos[bestOrderV[i]],2)); if(vertexDistanceSwap<vertexDistance && vertexDistanceSwap<0.5) { //check that best vertex is reasonable //swap! //arbitrary choice which view to switch: Int_t tempIndex = bestOrderU[i]; bestOrderU[i] = bestOrderU[i+1]; bestOrderU[i+1] = tempIndex; //start checks again so that clusters //can propagate if better matches are found i=0; } } } } //loop over number of expected showers: for(int i=0;i<nShowers;i++){ TObjArray newshower; Int_t nsubshower = 0; subshowerItr.ResetFirst(); while (CandSubShowerSRHandle *subshower = subshowerItr()) { if (*subshower->GetCandSlice()==*slice) { if(subshower->GetBegPlane()>=begPlane && subshower->GetEndPlane()<=endPlane) { if(subshower->GetPlaneView()==PlaneView::kU || subshower->GetPlaneView()==PlaneView::kX){ Float_t fracPHInRange = 0; Float_t totPHInSS = 0; CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); while (CandStripHandle *stp = stripssItr()) { totPHInSS += stp->GetCharge(CalDigitType::kSigCorr); if(stp->GetTPos()>peakUPosLow[bestOrderU[i]] && stp->GetTPos()<=peakUPosUpp[bestOrderU[i]]) fracPHInRange+=stp->GetCharge(CalDigitType::kSigCorr); } if(totPHInSS>=0 && fracPHInRange/totPHInSS >= 0.5) { newshower.Add(subshower); MSG("AlgShowerSS",Msg::kDebug) << "Adding SubShower: " << nsubshower << " in slice " << nslice << " to newshower " << i << endl; } } if(subshower->GetPlaneView()==PlaneView::kV || subshower->GetPlaneView()==PlaneView::kY){ Float_t fracPHInRange = 0; Float_t totPHInSS = 0; CandStripHandleItr stripssItr(subshower->GetDaughterIterator()); while (CandStripHandle *stp = stripssItr()) { totPHInSS += stp->GetCharge(CalDigitType::kSigCorr); if(stp->GetTPos()>peakVPosLow[bestOrderV[i]] && stp->GetTPos()<=peakVPosUpp[bestOrderV[i]]) fracPHInRange+=stp->GetCharge(CalDigitType::kSigCorr); } if(totPHInSS>=0 && fracPHInRange/totPHInSS >= 0.5) { newshower.Add(subshower); MSG("AlgShowerSS",Msg::kDebug) << "Adding SubShower: " << nsubshower << " in slice " << nslice << " to newshower " << i << endl; } } } } nsubshower+=1; } //now check that max/min strip times from subshowers do not have gaps //greater than 50ns //make vectors to hold min/max times per subshower //then loop through other subshowers looking for gaps>100ns //find "contiguous" time range and add up the energy //will end up with duplicate entries in the vectors since most subshowers //should be contiguous Int_t nAddedSS = newshower.GetLast()+1; vector<Double_t> contiguousMin(nAddedSS,0); vector<Double_t> contiguousMax(nAddedSS,0); vector<Double_t> totalEnergyInRange(nAddedSS,0); for(int j=0;j<nAddedSS;j++){ CandSubShowerSRHandle *subshower = dynamic_cast<CandSubShowerSRHandle*>(newshower.At(j)); contiguousMin[j] = 1e9*subshower->GetMinStpTime(); contiguousMax[j] = 1e9*subshower->GetMaxStpTime(); MSG("AlgShowerSS",Msg::kDebug) << "Initial " << j << " " << contiguousMin[j] << " " << contiguousMax[j] << endl; totalEnergyInRange[j] = subshower->GetEnergy(); for(int k=0;k<nAddedSS;k++){ if(j==k) continue; CandSubShowerSRHandle *subshower2 = dynamic_cast<CandSubShowerSRHandle*>(newshower.At(k)); Double_t max2 = 1e9*subshower2->GetMaxStpTime(); Double_t min2 = 1e9*subshower2->GetMinStpTime(); if(min2 < contiguousMin[j] && max2+SubShowerTimingCut < contiguousMin[j]) continue; if(min2-SubShowerTimingCut > contiguousMax[j] && max2 > contiguousMax[j]) continue; if(min2 > contiguousMin[j] && max2 < contiguousMax[j]) continue; if(min2 < contiguousMin[j] && max2 > contiguousMax[j]){ contiguousMin[j] = min2; contiguousMax[j] = max2; totalEnergyInRange[j] += subshower2->GetEnergy(); continue; } if(max2 <= contiguousMax[j] && max2+SubShowerTimingCut >= contiguousMin[j] && min2 < contiguousMin[j]) { contiguousMin[j] = min2; totalEnergyInRange[j] += subshower2->GetEnergy(); continue; } if(min2 >= contiguousMin[j] && min2-SubShowerTimingCut <= contiguousMax[j] && max2 > contiguousMax[j]) { contiguousMax[j] = max2; totalEnergyInRange[j] += subshower2->GetEnergy(); continue; } } MSG("AlgShowerSS",Msg::kDebug) << "Final " << j << " " << contiguousMin[j] << " " << contiguousMax[j] << endl; } //find largest energy range: Int_t rangeIndex = -1; Double_t largestEnergy = 0; for(int j=0;j<nAddedSS;j++){ if(totalEnergyInRange[j]>largestEnergy){ rangeIndex = j; largestEnergy = totalEnergyInRange[j]; } } if(nAddedSS>0) MSG("AlgShowerSS",Msg::kDebug) << "Largest " << rangeIndex << " " << largestEnergy << " " << contiguousMin[rangeIndex] << " " << contiguousMax[rangeIndex] << endl; //remove subshowers out of this range: for(int j=0;j<nAddedSS;j++){ CandSubShowerSRHandle *subshower = dynamic_cast<CandSubShowerSRHandle*>(newshower.At(j)); if(1e9*subshower->GetMinStpTime()<contiguousMin[rangeIndex] || 1e9*subshower->GetMaxStpTime()>contiguousMax[rangeIndex]) { newshower.RemoveAt(j); MSG("AlgShowerSS",Msg::kDebug) << "Removing SubShower: " << j << " in slice " << nslice << " from newshower " << i << endl; } } newshower.Compress(); cxx.SetDataIn(&newshower); CandShowerSRHandle showerhandle = CandShowerSR::MakeCandidate(ah,cxx); if((showerhandle.GetEnergy()>0 && nUPeaks_save>0 && nVPeaks_save>0) || showerhandle.GetEnergy()>MinNonPhysicsEnergy) { showerhandle.SetCandSlice(slice); ch.AddDaughterLink(showerhandle); MSG("AlgShowerSS",Msg::kDebug) << "Slice " << nslice << " Shower "<< nshower << " uid = " << showerhandle.GetUidInt() << endl; } else { MSG("AlgShowerSS",Msg::kDebug) << "Shower " << nshower << " has zero energy." << " - Not adding shower to list." << endl; } } } nslice++; } }

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

Reimplemented from AlgShowerSRList. Definition at line 1281 of file AlgShowerSSList.cxx. { }

---

Member Data Documentation

TH1F* AlgShowerSSList::vtxUHist [private]

Definition at line 27 of file AlgShowerSSList.h. Referenced by RunAlg().

TH1F* AlgShowerSSList::vtxUHistAll [private]

Definition at line 26 of file AlgShowerSSList.h. Referenced by RunAlg().

TH1F* AlgShowerSSList::vtxVHist [private]

Definition at line 29 of file AlgShowerSSList.h. Referenced by RunAlg().

TH1F* AlgShowerSSList::vtxVHistAll [private]

Definition at line 28 of file AlgShowerSSList.h. Referenced by RunAlg().

---

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

• AlgShowerSSList.h
• AlgShowerSSList.cxx

---