AlgShowerSR Class Reference

#include <AlgShowerSR.h>

Inheritance diagram for AlgShowerSR:

List of all members.

Public Member Functions
	AlgShowerSR ()
virtual	~AlgShowerSR ()
virtual void	RunAlg (AlgConfig &ac, CandHandle &ch, CandContext &cx)
virtual void	Trace (const char *c) const
Protected Member Functions
void	SetUVT (CandShowerHandle *)
void	SetUV (CandShowerHandle *)
void	SetT (CandShowerHandle *)
Int_t	FindTimingDirection (CandShowerHandle &shwh, AlgConfig &ac, Double_t *fitparm, Double_t &timefitchi2)

---

Constructor & Destructor Documentation

AlgShowerSR::AlgShowerSR (   )

Within AlgShowerSR::RunAlg, the vertex location and shower direction is determined. First, the most upstream U and V planes are found. Next, a loop over all shower clusters is used to obtain the total shower charge. The charge-weighted U/V positions of the CandClusters in the plane interval from the most upstream planes to VtxPlaneSpan planes downstream is identified as the shower vertex location. To obtain the shower direction, a linear fit is performed to the charge-weighted transverse position in a given plane vs the Z position of that plane. A weighted linear fit is performed in each viewed, weighted by the total cluster charge for this shower in a given plane. The du/dz and dv/dz slopes obtained in these fits are then used to obtain the direction cosines of the shower axis. Finally, the shower energy is set. Definition at line 42 of file AlgShowerSR.cxx. { }

AlgShowerSR::~AlgShowerSR (   )  [virtual]

Definition at line 47 of file AlgShowerSR.cxx. { }

---

Member Function Documentation

Int_t AlgShowerSR::FindTimingDirection (  CandShowerHandle &  shwh, AlgConfig &  ac, Double_t *  fitparm, Double_t &  timefitchi2 )  [protected]

Definition at line 801 of file AlgShowerSR.cxx. References Detector::AsString(), UgliStripHandle::ClearFiber(), digit(), CandDigitHandle::GetCharge(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), Registry::GetDouble(), PlexSEIdAltL::GetEnd(), UgliStripHandle::GetHalfLength(), Registry::GetInt(), CandHandle::GetNDaughters(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandDigitHandle::GetPlexSEIdAltL(), CandStripHandle::GetStripEndId(), UgliGeomHandle::GetStripHandle(), CandStripHandle::GetTime(), CandShowerHandle::GetU(), CandShowerHandle::GetV(), CandHandle::GetVldContext(), CandStripHandle::GetZPos(), MSG, PropagationVelocity::Velocity(), LinearFit::Weighted(), and UgliStripHandle::WlsPigtail(). Referenced by RunAlg(). { Double_t parmMaxTimeFitRes = ac.GetDouble("MaxTimeFitRes"); Double_t parmMinTimeFitSize = ac.GetInt("MinTimeFitSize"); Double_t parmMaxTimeFitOutlier = ac.GetDouble("MaxTimeFitOutlier"); MSG("AlgShowerSR", Msg::kVerbose) << "in find timing direction" << endl; fitparm[0] = 0; fitparm[1] = 0; CandStripHandle *strip = 0; CandDigitHandle *digit = 0; StripEnd::StripEnd_t stripEnd = StripEnd::kUnknown; //get VldContext and UgliGeomHandle CandStripHandleItr stripItr(shwh.GetDaughterIterator()); CandStripHandle *firststrip = stripItr.Ptr(); assert(firststrip); const VldContext *vldcontext = firststrip->GetVldContext(); assert(vldcontext); UgliGeomHandle ugh(*vldcontext); //arrays for determining the timing direction - limit to //1000 entries - really shouldnt need more than that to //be able to figure out if the shower is going north or //south Double_t time[1000]; Double_t pathLength[1000]; Double_t weight[1000]; Double_t adcsum[1000]; for(Int_t i = 0; i<1000; ++i){ time[i] = 0.; pathLength[i] = 0.; weight[i] = 0.; adcsum[i] = 0; } //need to put in a check to make sure you dont have more planes in one //view than the other Int_t uBeg = 1000; Int_t uEnd = 0; Int_t vBeg = 1000; Int_t vEnd = 0; while( (strip = stripItr())){ Int_t plane = strip->GetPlane(); if(strip->GetPlaneView() == PlaneView::kV && plane>vEnd) vEnd=plane; if(strip->GetPlaneView() == PlaneView::kV && plane<vBeg) vBeg=plane; if(strip->GetPlaneView() == PlaneView::kU && plane>uEnd) uEnd=plane; if(strip->GetPlaneView() == PlaneView::kU && plane<uBeg) uBeg=plane; } //get the endpoints Int_t showerBeg = TMath::Min(uBeg,vBeg); Int_t showerEnd = TMath::Max(uEnd,vEnd); if(TMath::Abs(uBeg-vBeg)>3.){ showerBeg = TMath::Max(uBeg,vBeg) - 1; } if(TMath::Abs(uEnd-vEnd)>3.){ showerEnd = TMath::Min(uEnd,vEnd) + 1; } MSG("AlgShowerSR", Msg::kDebug) << "u end points = " << uBeg << " " << uEnd << " v end points = " << vBeg << " " << vEnd << endl << " shower end points = " << showerBeg << " " << showerEnd << endl; Double_t halfLength = 0.; Double_t apos = 0.; Double_t distFromCenter = 0.; Double_t fiberDist = 0.; Int_t nctr = 0; Double_t maxPathLength = 0; for(Int_t iplane = showerBeg;iplane <= showerEnd;iplane++){ adcsum[nctr] = 0; time[nctr] = 0; pathLength[nctr] = 0; CandStripHandleItr shwstripItr(shwh.GetDaughterIterator()); while( (strip = shwstripItr()) && nctr<1000){ if(strip->GetPlane() == iplane){ MSG("AlgShowerSR", Msg::kDebug) << "strip from plane " << iplane << " " << Detector::AsString(vldcontext->GetDetector()) << endl; //only look at double ended strips if in the far detector //and strips from planes that have orthogonal measurements if( (strip->GetNDaughters() == 2 || vldcontext->GetDetector() == Detector::kNear)){ MSG("AlgShowerSR", Msg::kVerbose) << "strip good for timing" << endl; //get the iterator over the digits in the strip CandDigitHandleItr digitItr(strip->GetDaughterIterator()); while( (digit = digitItr())){ stripEnd = digit->GetPlexSEIdAltL().GetEnd(); pathLength[nctr] = strip->GetZPos(); maxPathLength = TMath::Max(maxPathLength,pathLength[nctr] ); UgliStripHandle stripHandle = ugh.GetStripHandle(strip->GetStripEndId()); halfLength = stripHandle.GetHalfLength(); apos = 0.; fiberDist = 0.; if(strip->GetPlaneView() == PlaneView::kV) apos = shwh.GetU(iplane); else if(strip->GetPlaneView() == PlaneView::kU) apos = shwh.GetV(iplane); distFromCenter = 0.; if(strip->GetPlaneView() == PlaneView::kV) distFromCenter = (stripEnd == StripEnd::kNegative) ? apos : -apos; else if(strip->GetPlaneView() == PlaneView::kU) distFromCenter = (stripEnd == StripEnd::kNegative) ? -apos : apos; fiberDist = (halfLength + distFromCenter + stripHandle.ClearFiber(stripEnd) + stripHandle.WlsPigtail(stripEnd)); time[nctr] += digit->GetCharge() * (strip->GetTime(stripEnd) - fiberDist/PropagationVelocity::Velocity()); weight[nctr] = 1.0; adcsum[nctr] += digit->GetCharge(); Double_t temptime = time[nctr]; if(adcsum[nctr]!=0)temptime /= adcsum[nctr]; MSG("AlgShowerSR", Msg::kDebug) << nctr << " " << pathLength[nctr] << " " << temptime << " " << strip->GetTime(stripEnd) << " " << weight[nctr] << endl; }//end loop over digits in strip }//end if double ended strip }//end loop over strips in plane } if(adcsum[nctr] > 0){ time[nctr] /= adcsum[nctr]; nctr++; } } for(Int_t i=0;i<nctr;i++){ MSG("AlgShowerSR",Msg::kDebug) << "TIMEFIT " << i << " " <<pathLength[i] << " " << (time[i]-time[0])*1.e9 << " " << weight[i] << endl; } //loop over all strips in the shower Double_t efitparm[2]; Double_t maxresid = parmMaxTimeFitRes + 1.; Double_t resid = 0.; Int_t ctr = nctr; Int_t imaxresid = 0; Int_t nremove = 0; Bool_t dofit = false; while(maxresid > parmMaxTimeFitRes && nctr-nremove > parmMinTimeFitSize && (1.*nremove) < parmMaxTimeFitOutlier*nctr ){ dofit = true; LinearFit::Weighted(ctr,pathLength,time,weight,fitparm,efitparm); maxresid = 0.; imaxresid = 0; for(int i=0; i<ctr; ++i){ resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i])*1.e9; if(weight[i]>0. && TMath::Abs(resid)>maxresid){ maxresid = TMath::Abs(resid); imaxresid = i; } } MSG("AlgShowerSR",Msg::kDebug) << "constrained fit, dt/ds slope, offset = " << fitparm[1] << " " << fitparm[0] << endl; MSG("AlgShowerSR",Msg::kDebug) << "maximum residual " << maxresid << " " << pathLength[imaxresid] << " " << time[imaxresid] << " " << weight[imaxresid] << endl; if(maxresid>parmMaxTimeFitRes){ MSG("AlgShowerSR",Msg::kDebug) << " removing" << endl; weight[imaxresid] = 0.; nremove++; } timefitchi2 = 0.; for(int i = 0; i < ctr; ++i){ resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i])*1e9; timefitchi2 += weight[i]*resid*resid; } MSG("AlgShowerSR",Msg::kDebug) << "chi2 = " << timefitchi2 \ << " n = " << nctr-nremove << endl; }//end loop to find timing fit and remove outliers timefitchi2 = 0.; if(dofit){ for(int i = 0; i<ctr; ++i){ if(weight[i] > 0.){ resid = (fitparm[1]*pathLength[i]+fitparm[0]-time[i])*1.e9; timefitchi2 += weight[i]*resid*resid; } } } MSG("AlgShowerSR",Msg::kDebug) << " TIMEFIT " << fitparm[1] << " chi^2/ndf = " << timefitchi2/(1.*(nctr-nremove)) << " " << TMath::Abs(uBeg-uEnd) << " " << TMath::Abs(vBeg-vEnd) << endl; //check the chi^2 for the fit - if it is really bad dont change the //initial guess at the direction for the event, ie just make sure //that the slope in time vs pathlength is positive; if(fitparm[1] < 0. && timefitchi2 >= 10.*(nctr-nremove)) fitparm[1] *= -1.; return nctr-nremove; }

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

Implements AlgBase. Reimplemented in AlgShowerSS. Definition at line 52 of file AlgShowerSR.cxx. References CandShowerHandle::AddCluster(), CandHandle::AddDaughterLink(), AlgReco::Calibrate(), CandShowerHandle::CalibrateEnergy(), FindTimingDirection(), CandClusterHandle::GetBegPlane(), CandStripHandle::GetBegTime(), CandStripHandle::GetCharge(), CandContext::GetDataIn(), CandHandle::GetDaughterIterator(), CandRecoHandle::GetEndPlane(), CandClusterHandle::GetEndPlane(), CandRecoHandle::GetEndT(), CandRecoHandle::GetEndZ(), Registry::GetInt(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), CandClusterHandle::GetPlaneView(), UgliGeomHandle::GetScintPlnHandle(), CandStripHandle::GetStripEndId(), CandStripHandle::GetTPos(), CandShowerHandle::GetU(), CandShowerHandle::GetV(), CandHandle::GetVldContext(), CandRecoHandle::GetVtxPlane(), CandRecoHandle::GetVtxT(), CandRecoHandle::GetVtxZ(), UgliPlnHandle::GetZ0(), MSG, CandRecoHandle::SetDirCosU(), CandRecoHandle::SetDirCosV(), CandRecoHandle::SetDirCosZ(), CandRecoHandle::SetEndDirCosU(), CandRecoHandle::SetEndDirCosV(), CandRecoHandle::SetEndDirCosZ(), CandRecoHandle::SetEndPlane(), CandRecoHandle::SetEndT(), CandRecoHandle::SetEndU(), CandRecoHandle::SetEndV(), CandRecoHandle::SetEndZ(), CandShowerHandle::SetEnergy(), SetUV(), CandRecoHandle::SetVtxPlane(), CandRecoHandle::SetVtxT(), CandRecoHandle::SetVtxU(), CandRecoHandle::SetVtxV(), CandRecoHandle::SetVtxZ(), and LinearFit::Weighted(). { assert(ch.InheritsFrom("CandShowerHandle")); CandShowerHandle &csh = dynamic_cast<CandShowerHandle &>(ch); assert(cx.GetDataIn()); assert(cx.GetDataIn()->InheritsFrom("TObjArray")); // grab parameters // Double_t MIPperGeV = ac.GetDouble("MIPperGeV"); Int_t vtxplanespan = ac.GetInt("VtxPlaneSpan"); Int_t isCosmic = ac.GetInt("IsCosmic"); const TObjArray *tary = dynamic_cast<const TObjArray*>(cx.GetDataIn()); Int_t ubegplane=0; Int_t vbegplane=0; Int_t uendplane=0; Int_t vendplane=0; CandStripHandle *begstrip=0; CandStripHandle *endstrip=0; for (Int_t i=0; i<=tary->GetLast(); i++) { TObject *tobj = tary->At(i); assert(tobj->InheritsFrom("CandClusterHandle")); CandClusterHandle *clusterhandle = dynamic_cast<CandClusterHandle*>(tobj); csh.AddCluster(clusterhandle); switch (clusterhandle->GetPlaneView()) { case PlaneView::kU: if (!ubegplane || clusterhandle->GetBegPlane()<ubegplane) { ubegplane = clusterhandle->GetBegPlane(); } if(!uendplane || clusterhandle->GetEndPlane()>uendplane) { uendplane = clusterhandle->GetEndPlane(); } break; case PlaneView::kV: if (!vbegplane || clusterhandle->GetBegPlane()<vbegplane) { vbegplane = clusterhandle->GetBegPlane(); } if(!vendplane || clusterhandle->GetEndPlane()>vendplane) { vendplane = clusterhandle->GetEndPlane(); } break; default: break; } if (!begstrip) { CandStripHandleItr stripItr(clusterhandle->GetDaughterIterator()); begstrip = stripItr(); } } // bail if these shower has an empty daughter list if(!begstrip){ MSG("ShowerSR", Msg::kError) << "Found empty daughter list in CandShowerSR" << endl; csh.SetDirCosU(0.); csh.SetDirCosV(0.); csh.SetDirCosZ(1.); csh.SetVtxPlane(-1); csh.SetVtxZ(-1.); csh.SetVtxU(0.); csh.SetVtxV(0.); csh.SetVtxT(0.); csh.SetEndZ(-1.); csh.SetEndU(0.); csh.SetEndV(0.); csh.SetEndT(0.); csh.SetEnergy(0); return; } const VldContext * vld = begstrip->GetVldContext(); UgliGeomHandle ugh(*vld); begstrip = 0; Double_t uvtx=0.; Double_t vvtx=0.; Double_t uvtxph=0.; Double_t vvtxph=0.; Double_t tvtx=0.; Double_t uend=0.; Double_t vend=0.; Double_t uendph=0.; Double_t vendph=0.; Double_t tend=0.; Bool_t firstvtx(1); Bool_t firstend(1); Double_t *uvtxfit = new Double_t[vtxplanespan+1]; Double_t *uvtxphfit = new Double_t[vtxplanespan+1]; Double_t *uvtxzfit = new Double_t[vtxplanespan+1]; Double_t *vvtxfit = new Double_t[vtxplanespan+1]; Double_t *vvtxphfit = new Double_t[vtxplanespan+1]; Double_t *vvtxzfit = new Double_t[vtxplanespan+1]; Double_t *uendfit = new Double_t[vtxplanespan+1]; Double_t *uendphfit = new Double_t[vtxplanespan+1]; Double_t *uendzfit = new Double_t[vtxplanespan+1]; Double_t *vendfit = new Double_t[vtxplanespan+1]; Double_t *vendphfit = new Double_t[vtxplanespan+1]; Double_t *vendzfit = new Double_t[vtxplanespan+1]; for (Int_t i=0; i<=vtxplanespan; i++) { uvtxfit[i] = 0.; uvtxphfit[i] = 0.; uvtxzfit[i] = 0.; vvtxfit[i] = 0.; vvtxphfit[i] = 0.; vvtxzfit[i] = 0.; uendfit[i] = 0.; uendphfit[i] = 0.; uendzfit[i] = 0.; vendfit[i] = 0.; vendphfit[i] = 0.; vendzfit[i] = 0.; } // loop over cluster daughters, calculating a charge sum, and // constructing arrays of charge-weighted u vs z and v vs z for the // planes within VtxPlaneSpan of beginning and end of shower Double_t charge=0.; for (Int_t i=0; i<=tary->GetLast(); i++) { TObject *tobj = tary->At(i); assert(tobj->InheritsFrom("CandClusterHandle")); CandClusterHandle *clusterhandle = dynamic_cast<CandClusterHandle*>(tobj); CandStripHandleItr stripItr(clusterhandle->GetDaughterIterator()); while (CandStripHandle *strip = stripItr()) { csh.AddDaughterLink(*strip); charge += strip->GetCharge(); if (!begstrip || strip->GetPlane()<begstrip->GetPlane()) { begstrip = strip; } if (!endstrip || strip->GetPlane()>endstrip->GetPlane()) { endstrip = strip; } switch (strip->GetPlaneView()) { case PlaneView::kU: if (strip->GetPlane()<=ubegplane+vtxplanespan) { Double_t charge = strip->GetCharge(); Double_t tpos = strip->GetTPos(); uvtx += tpos*charge; uvtxph += charge; if (firstvtx || strip->GetBegTime()<tvtx) { firstvtx = 0; tvtx = strip->GetBegTime(); } Int_t dplane = strip->GetPlane()-ubegplane; if(dplane<=vtxplanespan && dplane>+0){ uvtxfit[dplane] += tpos*charge; uvtxphfit[dplane] += charge; UgliScintPlnHandle upid = ugh.GetScintPlnHandle(strip->GetStripEndId()); uvtxzfit[dplane] = upid.GetZ0(); } } if(strip->GetPlane()>=uendplane-vtxplanespan) { Double_t charge = strip->GetCharge(); Double_t tpos = strip->GetTPos(); uend += tpos*charge; uendph += charge; if (firstend || strip->GetBegTime()<tend) { firstend = 0; tend = strip->GetBegTime(); } Int_t dplane = uendplane-strip->GetPlane(); if(dplane<=vtxplanespan && dplane>=0){ uendfit[dplane] += tpos*charge; uendphfit[dplane] += charge; UgliScintPlnHandle upid = ugh.GetScintPlnHandle(strip->GetStripEndId()); uendzfit[dplane] = upid.GetZ0(); } } break; case PlaneView::kV: if (strip->GetPlane()<=vbegplane+vtxplanespan) { Double_t charge = strip->GetCharge(); Double_t tpos = strip->GetTPos(); vvtx += tpos*charge; vvtxph += charge; if (firstvtx || strip->GetBegTime()<tvtx) { firstvtx = 0; tvtx = strip->GetBegTime(); } Int_t dplane = strip->GetPlane()-vbegplane; if(dplane<=vtxplanespan && dplane>=0){ vvtxfit[dplane] += tpos*charge; vvtxphfit[dplane] += charge; UgliScintPlnHandle upid = ugh.GetScintPlnHandle(strip->GetStripEndId()); vvtxzfit[dplane] = upid.GetZ0(); } } if(strip->GetPlane()>=uendplane-vtxplanespan) { Double_t charge = strip->GetCharge(); Double_t tpos = strip->GetTPos(); vend += tpos*charge; vendph += charge; if (firstend || strip->GetBegTime()<tend) { firstend = 0; tend = strip->GetBegTime(); } Int_t dplane = vendplane-strip->GetPlane(); if(dplane<=vtxplanespan && dplane>=0){ vendfit[dplane] += tpos*charge; vendphfit[dplane] += charge; UgliScintPlnHandle upid = ugh.GetScintPlnHandle(strip->GetStripEndId()); vendzfit[dplane] = upid.GetZ0(); } } break; default: break; } } } // normalize charge weighted fit points, and remove array entries // with zero pulse height Int_t iuvtxfit=0; Int_t ivvtxfit=0; Int_t iuendfit=0; Int_t ivendfit=0; for (Int_t i=0; i<=vtxplanespan; i++) { if (uvtxphfit[i]>0.) { uvtxfit[iuvtxfit] = uvtxfit[i]/uvtxphfit[i]; uvtxzfit[iuvtxfit]=uvtxzfit[i]; uvtxphfit[iuvtxfit]=uvtxphfit[i]; iuvtxfit++; } if (vvtxphfit[i]>0.) { vvtxfit[ivvtxfit] = vvtxfit[i]/vvtxphfit[i]; vvtxzfit[ivvtxfit]=vvtxzfit[i]; vvtxphfit[ivvtxfit]=vvtxphfit[i]; ivvtxfit++; } if (uendphfit[i]>0.) { uendfit[iuendfit] = uendfit[i]/uendphfit[i]; uendzfit[iuendfit]=uendzfit[i]; uendphfit[iuendfit]=uendphfit[i]; iuendfit++; } if (vendphfit[i]>0.) { vendfit[ivendfit] = vendfit[i]/vendphfit[i]; vendzfit[ivendfit]=vendzfit[i]; vendphfit[ivendfit]=vendphfit[i]; ivendfit++; } } //perform weighted fit of u, v vs z positions Double_t uvtxparm[2]; Double_t vvtxparm[2]; Double_t dudz=0.; Double_t dvdz=0.; if (iuvtxfit>1) { LinearFit::Weighted(iuvtxfit,uvtxzfit,uvtxfit,uvtxphfit,uvtxparm); dudz = uvtxparm[1]; } if (ivvtxfit>1) { LinearFit::Weighted(ivvtxfit,vvtxzfit,vvtxfit,vvtxphfit,vvtxparm); dvdz = vvtxparm[1]; } Double_t uendparm[2]; Double_t vendparm[2]; if (iuendfit>1) { LinearFit::Weighted(iuendfit,uendzfit,uendfit,uendphfit,uendparm); } if (ivendfit>1) { LinearFit::Weighted(ivendfit,vendzfit,vendfit,vendphfit,vendparm); } // calculate direction cosines from linear fits above Double_t cosz=1./sqrt(1.+dudz*dudz+dvdz*dvdz); Double_t cosu=dudz*cosz; Double_t cosv=dvdz*cosz; csh.SetDirCosU(cosu); csh.SetDirCosV(cosv); csh.SetDirCosZ(cosz); csh.SetEndDirCosU(cosu); csh.SetEndDirCosV(cosv); csh.SetEndDirCosZ(cosz); // set the vertex planes and z position. The vertex Z position // is the position of the most upstream strip. This is not // correct for cosmics! PlexStripEndId stripid = begstrip->GetStripEndId(); UgliScintPlnHandle planehandle = ugh.GetScintPlnHandle(stripid); csh.SetVtxPlane(begstrip->GetPlane()); csh.SetVtxZ(planehandle.GetZ0()); stripid = endstrip->GetStripEndId(); planehandle = ugh.GetScintPlnHandle(stripid); csh.SetEndPlane(endstrip->GetPlane()); csh.SetEndZ(planehandle.GetZ0()); // the vertex time is defined to be the earliest strip begin time in // the vertex region defined by VtxPlaneSpan csh.SetVtxT(tvtx); csh.SetEndT(tend); SetUV(&csh); // if this is a cosmic event, base direction on timing if(isCosmic){ Double_t fitparm[2] = {0,0}; Double_t timefitchi2; FindTimingDirection(csh, ac, fitparm, timefitchi2); if(fitparm[1]<0){ MSG("AlgShowerSR",Msg::kDebug) << " swapping vertex " << endl; Int_t newvtx = csh.GetEndPlane(); Int_t newend = csh.GetVtxPlane(); Double_t newvtxz = csh.GetEndZ(); Double_t newendz = csh.GetVtxZ(); Double_t newvtxt = csh.GetEndT(); Double_t newendt = csh.GetVtxT(); csh.SetVtxZ(newvtxz); csh.SetEndZ(newendz); csh.SetVtxPlane(newvtx); csh.SetEndPlane(newend); csh.SetVtxT(newvtxt); csh.SetEndT(newendt); } } // the vertex and end U/V positions are based on fits if the fit was done csh.SetVtxU(csh.GetU(csh.GetVtxPlane())); csh.SetVtxV(csh.GetV(csh.GetVtxPlane())); csh.SetEndU(csh.GetU(csh.GetEndPlane())); csh.SetEndV(csh.GetV(csh.GetEndPlane())); Calibrate(&csh); CandTrackHandle * associatedtrk=0; csh.CalibrateEnergy(associatedtrk,ac); delete [] uvtxfit; delete [] uvtxphfit; delete [] uvtxzfit; delete [] vvtxfit; delete [] vvtxphfit; delete [] vvtxzfit; delete [] uendfit; delete [] uendphfit; delete [] uendzfit; delete [] vendfit; delete [] vendphfit; delete [] vendzfit; }

void AlgShowerSR::SetT (  CandShowerHandle *   )  [protected]

Definition at line 653 of file AlgShowerSR.cxx. References UgliStripHandle::ClearFiber(), CandStripHandle::GetBegTime(), CandStripHandle::GetCalTimeType(), CandStripHandle::GetCharge(), CandHandle::GetDaughterIterator(), UgliStripHandle::GetHalfLength(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), VldContext::GetSimFlag(), CandStripHandle::GetStripEndId(), UgliGeomHandle::GetStripHandle(), CandShowerHandle::GetU(), CandShowerHandle::GetV(), CandHandle::GetVldContext(), max, min, MSG, CandShowerHandle::SetT(), PropagationVelocity::Velocity(), ArrivalTime::Weight(), and UgliStripHandle::WlsPigtail(). Referenced by SetUVT(). { // we take a weighted average of hits in the same // plane. the proper way to do this would be to consider hits coming // from the same PMT and use the earliest time. MSG("Shower",Msg::kDebug) << "starting SetT" << endl; TIter stripItr(candshower->GetDaughterIterator()); CandStripHandle *firststrip = 0; firststrip = dynamic_cast<CandStripHandle*>(stripItr()); UgliGeomHandle *ugh = 0; SimFlag::SimFlag_t simFlag = SimFlag::kUnknown; if (firststrip) { ugh = new UgliGeomHandle(*firststrip->GetVldContext()); simFlag = firststrip->GetVldContext()->GetSimFlag(); } Double_t timewgt[2] = {0.,0.}; Double_t time[2] = {0.,0.}; StripEnd::StripEnd_t stripend[2] = {StripEnd::kNegative,StripEnd::kPositive}; Double_t apos; ArrivalTime arrtime; CalTimeType::CalTimeType_t caltimetype = CalTimeType::kNone; Double_t propagation_velocity = PropagationVelocity::Velocity(simFlag)*Munits::ns; CandStripHandleItr orderedstripItr(candshower->GetDaughterIterator()); CandStripHandleKeyFunc *stripKf = orderedstripItr.CreateKeyFunc(); stripKf->SetFun(CandStripHandle::KeyFromPlane); orderedstripItr.GetSet()->AdoptSortKeyFunc(stripKf); stripKf = 0; Int_t oldplane = -1; PlaneView::PlaneView_t oldview = PlaneView::kUnknown; Bool_t first(1); while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(orderedstripItr())) { if (strip->GetPlaneView()==PlaneView::kU || strip->GetPlaneView()==PlaneView::kV) { if (strip->GetPlaneView()==PlaneView::kU) { apos = candshower->GetV(strip->GetPlane()); } else { apos = candshower->GetU(strip->GetPlane()); } apos=min(apos,4.0); apos=max(apos,-4.0); // sanity check on apos bounds UgliStripHandle striphandle = ugh->GetStripHandle(strip->GetStripEndId()); caltimetype = strip->GetCalTimeType(); if (first || strip->GetPlane()==oldplane) { if (first) { first = 0; oldplane = strip->GetPlane(); oldview = strip->GetPlaneView(); } Double_t thisph[2] = { strip->GetCharge(CalDigitType::kPE,stripend[0]), strip->GetCharge(CalDigitType::kPE,stripend[1])}; Double_t thiswgt[2] = { min(0.4,arrtime.Weight(max(1.,thisph[0]))), min(0.4,arrtime.Weight(max(1.,thisph[1])))}; for (int i=0; i<2; i++) { if (thisph[i]>0.) { Double_t plusminus = 1.; if ((strip->GetPlaneView()==PlaneView::kV && i==1) || (strip->GetPlaneView()==PlaneView::kU && i==0)) { plusminus = -1.; } Double_t thistime = strip->GetBegTime(stripend[i]) - 1.e-9*(striphandle.ClearFiber(stripend[i]) + striphandle.WlsPigtail(stripend[i]) + (striphandle.GetHalfLength() + plusminus*apos))/propagation_velocity; MSG("Shower",Msg::kDebug) << "plane " << strip->GetPlane() << " end " << i << " raw time " << strip->GetBegTime(stripend[i])*1e9 << " pulse height " << thisph[i] << " weight " << thiswgt[i] << " clearfiber " << striphandle.ClearFiber(stripend[i]) << " wlspigtail " << striphandle.WlsPigtail(stripend[i]) << " halflength " << striphandle.GetHalfLength() << " apos " << apos << " 3D time " << thistime*1e9 << endl; time[i] += thiswgt[i]*thistime; timewgt[i] += thiswgt[i]; MSG("Shower",Msg::kDebug) << "summed 3D time " << time[i]/timewgt[i] << endl; } } } else { for (int i=0; i<2; i++) { if (timewgt[i]>0.) { time[i] /= timewgt[i]; candshower->SetT(oldplane,stripend[i],time[i]); MSG("Shower",Msg::kDebug) << "setting time " << oldplane << " " << stripend[i] << " " << time[i]*1e9 << endl; } } oldplane = strip->GetPlane(); oldview = strip->GetPlaneView(); Double_t thisph[2] = { strip->GetCharge(CalDigitType::kPE,stripend[0]), strip->GetCharge(CalDigitType::kPE,stripend[1])}; Double_t thiswgt[2] = { min(0.4,arrtime.Weight(max(1.,thisph[0]))), min(0.4,arrtime.Weight(max(1.,thisph[1])))}; for (int i=0; i<2; i++) { time[i] = 0.; timewgt[i] = 0.; if (thisph[i]>0.) { Double_t plusminus = 1.; if ((strip->GetPlaneView()==PlaneView::kV && i==1) || (strip->GetPlaneView()==PlaneView::kU && i==0)) { plusminus = -1.; } Double_t thistime = strip->GetBegTime(stripend[i]) - 1.e-9*(striphandle.ClearFiber(stripend[i]) + striphandle.WlsPigtail(stripend[i]) + (striphandle.GetHalfLength() + plusminus*apos))/propagation_velocity; MSG("Shower",Msg::kDebug) << "plane " << strip->GetPlane() << " end " << i << " raw time " << strip->GetBegTime(stripend[i])*1e9 << " pulse height " << thisph[i] << " weight " << thiswgt[i] << " clearfiber " << striphandle.ClearFiber(stripend[i]) << " wlspigtail " << striphandle.WlsPigtail(stripend[i]) << " halflength " << striphandle.GetHalfLength() << " apos " << apos << " 3D time " << thistime*1e9 << endl; time[i] = thiswgt[i]*thistime; timewgt[i] = thiswgt[i]; MSG("Shower",Msg::kDebug) << "summed 3D time " << time[i]/timewgt[i] << endl; } } } } } for (int i=0; i<2; i++) { if (timewgt[i]>0.) { time[i] /= timewgt[i]; candshower->SetT(oldplane,stripend[i],time[i]); MSG("Shower",Msg::kDebug) << "setting time " << oldplane << " " << stripend[i] << " " << time[i]*1e9 << endl; } } if (ugh) delete ugh; }

void AlgShowerSR::SetUV (  CandShowerHandle *   )  [protected]

Definition at line 443 of file AlgShowerSR.cxx. References PlexPlaneId::GetAdjoinScint(), CandRecoHandle::GetBegPlane(), CandStripHandle::GetCharge(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), CandRecoHandle::GetEndPlane(), CandStripHandle::GetNDigit(), PlexPlaneId::GetPlane(), CandStripHandle::GetPlane(), CandStripHandle::GetPlaneView(), UgliGeomHandle::GetScintPlnHandle(), CandStripHandle::GetTPos(), CandShowerHandle::GetU(), CandShowerHandle::GetV(), CandHandle::GetVldContext(), CandShowerHandle::GetZ(), UgliPlnHandle::GetZ0(), UgliScintPlnHandle::IsValid(), PlexPlaneId::IsValid(), max, min, MSG, CandShowerHandle::SetU(), and CandShowerHandle::SetV(). Referenced by RunAlg(), and SetUVT(). { TIter stripItr(candshower->GetDaughterIterator()); CandStripHandle *firststrip = 0; firststrip = dynamic_cast<CandStripHandle*>(stripItr()); if (!firststrip) return; // if no strips, do nothing const VldContext *vldc = firststrip->GetVldContext(); Detector::Detector_t det = vldc->GetDetector(); UgliGeomHandle ugh = UgliGeomHandle(*vldc); CandStripHandleItr orderedstripItr(candshower->GetDaughterIterator()); CandStripHandleKeyFunc *stripKf = orderedstripItr.CreateKeyFunc(); stripKf->SetFun(CandStripHandle::KeyFromPlane); orderedstripItr.GetSet()->AdoptSortKeyFunc(stripKf); stripKf = 0; Bool_t first(1); Int_t oldplane = -1; PlaneView::PlaneView_t oldview = PlaneView::kUnknown; Float_t tpos = 0.; Float_t ph = 0.; Int_t idir=1; Int_t vtxplane = candshower->GetBegPlane(); Int_t endplane = candshower->GetEndPlane(); if (vtxplane>endplane) idir=-1; Int_t begplane = min(vtxplane,endplane); endplane = max(vtxplane,endplane); begplane = -1; endplane = -1; // find begin and end planes while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(orderedstripItr())) { if (strip->GetPlaneView()==PlaneView::kU || strip->GetPlaneView()==PlaneView::kV) { if (begplane<0 || strip->GetPlane()<begplane) { begplane = strip->GetPlane(); } if (endplane<0 || strip->GetPlane()>endplane) { endplane = strip->GetPlane(); } } } Int_t *planelist = new Int_t[endplane+1]; Int_t *stripendlist = new Int_t[endplane+1]; if(begplane>-1 && endplane>-1){ for (int i=0; i<=endplane; i++) { planelist[i] = 0; stripendlist[i] = 0; } // transverse positions orderedstripItr.Reset(); while (CandStripHandle *strip = dynamic_cast<CandStripHandle*>(orderedstripItr())) { if (strip->GetPlaneView()==PlaneView::kU || strip->GetPlaneView()==PlaneView::kV) { planelist[strip->GetPlane()] = strip->GetPlaneView(); if (strip->GetNDigit(StripEnd::kNegative)) { stripendlist[strip->GetPlane()] += (1<<StripEnd::kNegative); } if (strip->GetNDigit(StripEnd::kPositive)) { stripendlist[strip->GetPlane()] += (1<<StripEnd::kPositive); } if (first || strip->GetPlane()==oldplane) { if (first) { first = 0; oldplane = strip->GetPlane(); oldview = strip->GetPlaneView(); } Float_t stripph = strip->GetCharge(); tpos += strip->GetTPos()*stripph; ph += stripph; } else { tpos /= ph; if (oldview==PlaneView::kU) { candshower->SetU(oldplane,tpos); } else if (oldview==PlaneView::kV) { candshower->SetV(oldplane,tpos); } oldplane = strip->GetPlane(); oldview = strip->GetPlaneView(); Float_t stripph = strip->GetCharge(); tpos = strip->GetTPos()*stripph; ph = stripph; } } } if (ph>0.) { tpos /= ph; if (oldview==PlaneView::kU) { candshower->SetU(oldplane,tpos); } else if (oldview==PlaneView::kV) { candshower->SetV(oldplane,tpos); } } // create a starting and stopping PlexPlaneId Int_t idirAdjoin = (begplane>endplane) ? -1 : 1; // direction PlexPlaneId planeid(det,begplane,false); PlexPlaneId planeidend(det,endplane,false); PlexPlaneId planeid_toofar = planeidend.GetAdjoinScint(idirAdjoin); // after each step, if not reached end, move to adjoining *scint* plane for ( ; planeid != planeid_toofar ; planeid = planeid.GetAdjoinScint(idirAdjoin) ) { // quit if hit something meaningless if ( ! planeid.IsValid() ) break; UgliScintPlnHandle scintpln = ugh.GetScintPlnHandle(planeid); if ( ! scintpln.IsValid() ) continue; Int_t iplane = planeid.GetPlane(); Float_t zpos = scintpln.GetZ0(); for (Int_t planeview=PlaneView::kU; planeview<=PlaneView::kV; planeview++) { if (planelist[iplane]!=planeview) { Int_t found = 0; Int_t nfound = 0; Int_t jplane[2]; // search upstream for (Int_t iplane1=iplane-1; iplane1>=begplane && !found; iplane1--) { if (planelist[iplane1]==planeview) { found++; jplane[nfound] = iplane1; nfound++; } } Int_t ndown = 1; if (!found) { ndown++; } found = 0; // search downstream for (Int_t iplane1=iplane+1; iplane1<=endplane && found<ndown; iplane1++) { if (planelist[iplane1]==planeview) { found++; jplane[nfound] = iplane1; nfound++; } } if (!found) { for (Int_t iplane1=jplane[0]-1; iplane1>=begplane && !found; iplane1--) { if (planelist[iplane1]==planeview) { found++; jplane[nfound] = iplane1; nfound++; } } } if (nfound==1) { if (planeview==PlaneView::kU) { candshower->SetU(iplane,candshower->GetU(jplane[0])); } else { candshower->SetV(iplane,candshower->GetV(jplane[0])); } } else { Float_t z[2]; z[0] = candshower->GetZ(jplane[0]); z[1] = candshower->GetZ(jplane[1]); Float_t x[2]; if (planeview==PlaneView::kU) { x[0] = candshower->GetU(jplane[0]); x[1] = candshower->GetU(jplane[1]); Float_t denom = (z[1]-z[0]); if (denom) candshower->SetU(iplane, x[0] + (zpos-z[0]) * ((x[1]-x[0])/denom)); else { MSG("Shower",Msg::kWarning) << endl << "Denominator (z[1]-z[0]) is zero." << " SetU not called." << endl; } } else { x[0] = candshower->GetV(jplane[0]); x[1] = candshower->GetV(jplane[1]); Float_t denom = (z[1]-z[0]); if (denom) candshower->SetV(iplane, x[0] + (zpos-z[0]) * ((x[1]-x[0])/denom)); else { MSG("Shower",Msg::kWarning) << endl << "Denominator is (z[1]-z[0]) zero." << " SetV not called." << endl; } } } } } } } delete [] planelist; delete [] stripendlist; }

void AlgShowerSR::SetUVT (  CandShowerHandle *   )  [protected]

Definition at line 435 of file AlgShowerSR.cxx. References SetT(), and SetUV(). Referenced by AlgShowerSS::RunAlg(). { this->SetUV(candshower); this->SetT(candshower); }

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

Reimplemented from AlgBase. Reimplemented in AlgShowerSS. Definition at line 1021 of file AlgShowerSR.cxx. { }

---

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

• AlgShowerSR.h
• AlgShowerSR.cxx

---