AlgDeMuxBeam Class Reference

#include <AlgDeMuxBeam.h>

Inheritance diagram for AlgDeMuxBeam:

List of all members.

Public Member Functions
	AlgDeMuxBeam ()
virtual	~AlgDeMuxBeam ()
virtual void	Trace (const char *c) const
virtual void	RunAlg (AlgConfig &acd, CandHandle &ch, CandContext &cx)
Private Member Functions
void	DeMuxFirstNPlanesTest (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr)
void	FindShowerWindowSolution (DmxPlaneItr &windowPlaneItr, Float_t cog, Int_t firstPlane, Int_t lastPlane)
void	SetFirstNPlanes (DmxPlaneItr &planeItr, Float_t slope, Float_t intercept)
void	SetMuonRegionWindow (DmxPlaneItr &planeItr, Float_t slope, Float_t intercept, bool setAll)
void	SetShowerRegionWindow (DmxPlaneItr &planeItr, Float_t slope, Float_t intercept)
void	UseShowerSlidingWindow (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, DmxStatus *status)
void	UseMuonSlidingWindow (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, DmxStatus *status)
void	UseReverseMuonSlidingWindow (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, DmxStatus *status)
void	RegionLinearFit (Float_t &slope, Float_t &intercept, Float_t &chiSqr, Float_t *x, Float_t *y, Float_t *sigma, Float_t *residual, Int_t numPoints)
void	FindFitOverNPlanes (Float_t *x, Float_t *y, Float_t *sigma, Float_t *residual, Float_t &slope, Float_t &intercept, Int_t numPoints, Int_t *planeNum)
Float_t	FindRegionMatedChargeFraction (DmxPlaneItr &validPlaneItr)
void	ReconcileShowerAndMuonRegions (DmxPlaneItr &planeItr, CandDeMuxDigitHandleItr &cdhit)
void	DeMuxFirstNPlanes (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, CandDeMuxDigitHandleItr &cdhit)
Float_t	FindDigitCoG (DmxPlaneItr &validPlaneItr, CandDeMuxDigitHandleItr &cdhit)
void	FindFit (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &a4, Double_t &bestChi2, Int_t leverArm, Int_t *hypotheses)
void	FindFitForHypothesisSet (DmxPlaneItr &planeItr, Double_t &prevChiSq, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &a4, Int_t *hypotheses, Int_t leverArm)
Bool_t	FindPlanesToDropFromFit (DmxPlaneItr &planeItr, Double_t a1, Double_t a2, Double_t a3, Double_t a4, Int_t leverArm)
void	FindRMSOfPlanes (DmxPlaneItr &planeItr, Int_t leverArm, Int_t *hypsToUse, Float_t &cog, Float_t &rms)
Private Attributes
Float_t	fFinalShowerRegionSlope
Float_t	fFinalShowerRegionIntercept
Float_t	fInitialMuonRegionSlope
Float_t	fInitialMuonRegionIntercept
Float_t	fDigitResetCut
Int_t	fMaxBadResidual
Float_t	fMaxResidual
Float_t	fMaxMuonRemovalFraction
Int_t	fMaxStripAdjustment
Int_t	fMinMuonPlanesRequired
Float_t	fMinMatedChargeFrac
Int_t	fMuonPlanesInSet
Float_t	fMuonSlopeLimit
Int_t	fMuonStartPlaneNumber
Float_t	fMuonStartPlaneZPos
Int_t	fPlanesInSet
Int_t	fStrayCut
Int_t	fStrayPlanes
Int_t	fNumberOfHypotheses
Float_t	fShowerMuonTransverseDifference
VldContext	fVldContext

---

Constructor & Destructor Documentation

AlgDeMuxBeam::AlgDeMuxBeam (   )

Definition at line 72 of file AlgDeMuxBeam.cxx. : fMaxResidual(0.1), fMaxMuonRemovalFraction(0.2), fMuonPlanesInSet(5), fMuonStartPlaneNumber(500), fStrayCut(12), fStrayPlanes(0) { //default constructor }

AlgDeMuxBeam::~AlgDeMuxBeam (   )  [virtual]

Definition at line 85 of file AlgDeMuxBeam.cxx. { }

---

Member Function Documentation

void AlgDeMuxBeam::DeMuxFirstNPlanes (  DmxPlaneItr &  planeItr, DmxPlaneItr &  validPlaneItr, CandDeMuxDigitHandleItr &  cdhit )  [private]

Definition at line 2008 of file AlgDeMuxBeam.cxx. References FindFit(), DmxPlane::GetCoG(), DmxPlane::GetPlaneType(), DmxPlane::GetZPosition(), MSG, SetFirstNPlanes(), and DmxPlane::SetStrips(). { Int_t hypsToUse[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0} ; Int_t leverArm = (Int_t)validPlaneItr.SizeSelect(); if(leverArm > fPlanesInSet) leverArm = fPlanesInSet; //find the bounds for the first n planes MSG("AlgDeMuxBeam", Msg::kDebug) << "In DeMuxFirstNPlanes with lever arm = " << leverArm << endl; //useA2 is the slope, useA1 is the intercept //variables to keep track of the best reconstruction Double_t useA1 = -10000.; Double_t useA2 = -10000.; Double_t useA3 = -10000.; Double_t useA4 = -10000.; Double_t bestChi2 = 1.e20; cdhit.ResetFirst(); validPlaneItr.ResetFirst(); planeItr.ResetFirst(); if(leverArm == 2){ //only two planes in the set. if they are shower planes, just set it //each plane to the best hypothesis. valid muon planes are already set. //get the slope and intercept for a straight line through the two and //then set all remaining planes based on that Float_t z1 = 0., z2 = 0.; Float_t tPos1 = 0., tPos2 = 0.; Int_t ctr = 0; DmxPlane *plane = 0; while( (plane = validPlaneItr()) ){ if(plane->GetPlaneType() == DmxPlaneTypes::kShower) plane->SetStrips(); if(ctr == 0){ z1 = plane->GetZPosition(); tPos1 = plane->GetCoG(); } else{ z2 = plane->GetZPosition(); tPos2 = plane->GetCoG(); } ++ctr; }//end loop to find the first and last plane tpos and zpos if(z1 != z2) useA2 = (tPos2 - tPos1)/(z2 - z1); useA1 = tPos2 - useA2*z2; }//end if only 2 valid planes else if(leverArm==3){ //evaluate the rms of the centers of gravity - take the tightest //distribution for(Int_t ia = 1; ia < 4; ++ia){ hypsToUse[0] = ia; for(Int_t ib = 1; ib < 4; ++ib){ hypsToUse[1] = ib; for(Int_t ic = 1; ic < 4; ++ic){ hypsToUse[2] = ic; FindFit(planeItr, validPlaneItr, useA1, useA2, useA3, useA4, bestChi2, leverArm, hypsToUse); } } } }//end if lever arm = 3 else if(leverArm == 4){ for(Int_t ia = 1; ia < 4; ++ia){ hypsToUse[0] = ia; for(Int_t ib = 1; ib < 4; ++ib){ hypsToUse[1] = ib; for(Int_t ic = 1; ic < 4; ++ic){ hypsToUse[2] = ic; for(Int_t id = 1; id < 4; ++id){ hypsToUse[3] = id; FindFit(planeItr, validPlaneItr, useA1, useA2, useA3, useA4, bestChi2, leverArm, hypsToUse); } } } } }//end if 4 plane lever arm else if(leverArm == 5){ for(Int_t ia = 1; ia < 4; ++ia){ hypsToUse[0] = ia; for(Int_t ib = 1; ib < 4; ++ib){ hypsToUse[1] = ib; for(Int_t ic = 1; ic < 4; ++ic){ hypsToUse[2] = ic; for(Int_t id = 1; id < 4; ++id){ hypsToUse[3] = id; for(Int_t ie = 1; ie < 4; ++ie){ hypsToUse[4] = ie; FindFit(planeItr, validPlaneItr, useA1, useA2, useA3, useA4, bestChi2, leverArm, hypsToUse); } } } } } } else if(leverArm == 6){ for(Int_t ia = 1; ia < 4; ++ia){ hypsToUse[0] = ia; for(Int_t ib = 1; ib < 4; ++ib){ hypsToUse[1] = ib; for(Int_t ic = 1; ic < 4; ++ic){ hypsToUse[2] = ic; for(Int_t id = 1; id < 4; ++id){ hypsToUse[3] = id; for(Int_t ie = 1; ie < 4; ++ie){ hypsToUse[4] = ie; for(Int_t ig = 1; ig < 4; ++ig){ hypsToUse[5] = ig; FindFit(planeItr, validPlaneItr, useA1, useA2, useA3, useA4, bestChi2, leverArm, hypsToUse); } } } } } } }//end if lever arm = 6 planes MSG("AlgDeMuxBeam", Msg::kDebug) << "\tfinal fit" << "\t" << useA1 << "\t" << useA2 << "\t" << bestChi2 << endl; SetFirstNPlanes(planeItr, useA2, useA1); return; }

void AlgDeMuxBeam::DeMuxFirstNPlanesTest (  DmxPlaneItr &  planeItr, DmxPlaneItr &  validPlaneItr )  [private]

Definition at line 646 of file AlgDeMuxBeam.cxx. References fFinalShowerRegionIntercept, fFinalShowerRegionSlope, FindFitOverNPlanes(), fVldContext, DmxMuonPlane::GetCoG(), DmxShowerPlane::GetHypothesis(), DmxShowerPlane::GetHypothesisCoG(), DmxMuonPlane::GetInitialCoG(), DmxMuonPlane::GetInitialStripCoG(), DmxHypothesis::GetMatedSignalRatio(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), UgliGeomHandle::GetStripHandle(), UgliStripHandle::GetTPos(), DmxPlane::GetZPosition(), MSG, and SetFirstNPlanes(). Referenced by RunAlg(). { Int_t leverArm = (Int_t)validPlaneItr.SizeSelect(); if(leverArm > fPlanesInSet) leverArm = fPlanesInSet; Int_t planeCtr = 0; //loop over the valid planes and take the first leverArm planes. DmxPlane *plane = 0; DmxShowerPlane *showerPlane = 0; DmxMuonPlane *muonPlane = 0; validPlaneItr.ResetFirst(); Float_t matedCharge = 0.; Float_t planeCharge = 0.; Float_t maxMatedCharge = 0.; Int_t maxMatedChargeHyp = 0; Float_t zPos[] = {0., 0., 0., 0., 0., 0.}; Float_t tPos[] = {0., 0., 0., 0., 0., 0.}; Float_t weight[] = {1., 1., 1., 1., 1., 1.}; Float_t residual[] = {0., 0., 0., 0., 0., 0.}; Int_t planeNum[] = {0,0,0,0,0,0}; Float_t slope = 0.; Float_t intercept = 0.; Float_t averageRMS = 0.; Int_t rmsCtr = 0; //now loop over all hypotheses in turn for the planes and see which set gives you //the most mated signal over the leverArm planes for(Int_t i = 0; i<fNumberOfHypotheses; ++i){ validPlaneItr.ResetFirst(); planeCtr = 0; matedCharge = 0.; rmsCtr = 0; averageRMS = 0.; while( (plane = validPlaneItr()) && planeCtr<leverArm){ planeCharge = plane->GetPlaneCharge(); //check out the type of plane you have here if(plane->GetPlaneType() == DmxPlaneTypes::kShower){ matedCharge += dynamic_cast<DmxShowerPlane *>(plane)->GetHypothesis(i)->GetMatedSignalRatio()*planeCharge; averageRMS += dynamic_cast<DmxShowerPlane *>(plane)->GetHypothesis(i)->GetTieBreakerStat(); ++rmsCtr; MSG("AlgDeMuxBeam", Msg::kDebug) << "plane = " << plane->GetPlaneNumber() << " i = " << i << " mated frac = " << dynamic_cast<DmxShowerPlane *>(plane)->GetHypothesis(i)->GetMatedSignalRatio() << " rms = " << dynamic_cast<DmxShowerPlane *>(plane)->GetHypothesis(i)->GetTieBreakerStat() << " charge = " << planeCharge << endl; } //its a valid muon plane so check to see if the current hypothesis contains the initial strip else{ if(dynamic_cast<DmxMuonPlane *>(plane)->GetInitialStripCoG()>=i && dynamic_cast<DmxMuonPlane *>(plane)->GetInitialStripCoG()<=i+23) matedCharge += planeCharge; } ++planeCtr; }//end loop over first leverArm valid planes //does the current hypothesis represent the most mated signal? if(matedCharge>maxMatedCharge){ maxMatedCharge = matedCharge; //minAverageRMS = averageRMS/(1.*rmsCtr); maxMatedChargeHyp = i; } }//end loop over hypotheses MSG("AlgDeMuxBeam", Msg::kDebug) << "best initial hypothesis = " << maxMatedChargeHyp << endl; validPlaneItr.ResetFirst(); planeCtr = 0; //make an array holding the best hypothesis numbers for these planes Int_t bestHyps[] = {maxMatedChargeHyp, maxMatedChargeHyp, maxMatedChargeHyp, maxMatedChargeHyp, maxMatedChargeHyp, maxMatedChargeHyp}; //get the transverse position for the best hypothesis UgliGeomHandle ugh(fVldContext); //check the number of valid planes - if only 2, then just set these suckers to the best hyp and be done with it. if(leverArm == 2){ //get the strip end id for the center of gravity for the center of the best hypothesis PlexStripEndId seid(Detector::kFar,dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetPlaneNumber(),maxMatedChargeHyp+12); UgliStripHandle ush = ugh.GetStripHandle(seid); intercept = ush.GetTPos(); //set the first n planes - slope is 0 because you are just setting all the planes to one //hypothesis SetFirstNPlanes(planeItr, 0., intercept); return; } //loop over the valid planes again and see if by adjusting the hypothesis number you can get more //mated signal for the plane while ((plane = validPlaneItr()) && planeCtr<leverArm){ maxMatedCharge = 0.; maxMatedChargeHyp = bestHyps[planeCtr]; zPos[planeCtr] = plane->GetZPosition(); //get the strip end id for the center of gravity for the center of the best hypothesis PlexStripEndId seid(Detector::kFar,plane->GetPlaneNumber(),bestHyps[planeCtr]+12); UgliStripHandle ush = ugh.GetStripHandle(seid); tPos[planeCtr] = ush.GetTPos(); planeNum[planeCtr] = plane->GetPlaneNumber(); planeCharge = plane->GetPlaneCharge(); showerPlane = 0; muonPlane = 0; if( plane->GetPlaneType() == DmxPlaneTypes::kShower) showerPlane = dynamic_cast<DmxShowerPlane *>(plane); else muonPlane = dynamic_cast<DmxMuonPlane *>(plane); if(showerPlane){ tPos[planeCtr] = showerPlane->GetHypothesisCoG(maxMatedChargeHyp); maxMatedCharge = showerPlane->GetHypothesis(bestHyps[planeCtr])->GetMatedSignalRatio()*planeCharge; } else{ if(muonPlane->GetInitialStripCoG()>=maxMatedChargeHyp && muonPlane->GetInitialStripCoG()<=maxMatedChargeHyp+23){ tPos[planeCtr] = muonPlane->GetCoG(); maxMatedCharge = planeCharge; } } //look at hypotheses n strips on either side of the chosen hyp for the plane for(Int_t i = 1; i<=fMaxStripAdjustment; ++i){ if(showerPlane){ if(bestHyps[planeCtr]-i>=0){ matedCharge = showerPlane->GetHypothesis(bestHyps[planeCtr]-i)->GetMatedSignalRatio()*planeCharge; if(matedCharge>maxMatedCharge){ MSG("AlgDeMuxBeam", Msg::kDebug) << "initial best hyp = " << bestHyps[planeCtr] << " new best hyp = " << bestHyps[planeCtr]-i << " new tpos = " << showerPlane->GetHypothesisCoG(maxMatedChargeHyp) << endl; maxMatedCharge = matedCharge; maxMatedChargeHyp = bestHyps[planeCtr]-i; tPos[planeCtr] = showerPlane->GetHypothesisCoG(maxMatedChargeHyp); } }//end if the new hypothesis is allowed if(bestHyps[planeCtr]+i<fNumberOfHypotheses){ matedCharge = showerPlane->GetHypothesis(bestHyps[planeCtr]+i)->GetMatedSignalRatio()*planeCharge; if(matedCharge>maxMatedCharge){ MSG("AlgDeMuxBeam", Msg::kDebug) << "initial best hyp = " << bestHyps[planeCtr] << " new best hyp = " << bestHyps[planeCtr]+i << " new tpos = " << showerPlane->GetHypothesisCoG(maxMatedChargeHyp)<< endl; maxMatedCharge = matedCharge; maxMatedChargeHyp = bestHyps[planeCtr]+i; tPos[planeCtr] = showerPlane->GetHypothesisCoG(maxMatedChargeHyp); } }//end if the other new hypothesis is allowed }//end if it is a shower plane if(muonPlane){ if(bestHyps[planeCtr]+i<fNumberOfHypotheses){ if(muonPlane->GetInitialStripCoG()>=bestHyps[planeCtr]+i && muonPlane->GetInitialStripCoG()<=bestHyps[planeCtr]+i+23){ //this hypothesis contains the original decoding for the valid muon plane //so the mated charge is all charge on the plane and the transverse position is //the original cog MSG("AlgDeMuxBeam", Msg::kDebug) << "initial best hyp = " << bestHyps[planeCtr] << " new best hyp = " << bestHyps[planeCtr]+i << " new tpos = " << muonPlane->GetCoG()<< endl; maxMatedChargeHyp = bestHyps[planeCtr] + i; maxMatedCharge = planeCharge; tPos[planeCtr] = muonPlane->GetInitialCoG(); } }//end if the new hypothesis is allowed if(bestHyps[planeCtr]-i>=0){ if(muonPlane->GetInitialStripCoG()>=bestHyps[planeCtr]-i && muonPlane->GetInitialStripCoG()<=bestHyps[planeCtr]-i+23){ //this hypothesis contains the original decoding for the valid muon plane //so the mated charge is all charge on the plane and the transverse position is //the original cog MSG("AlgDeMuxBeam", Msg::kDebug) << "initial best hyp = " << bestHyps[planeCtr] << " new best hyp = " << bestHyps[planeCtr]-i << " new tpos = " << muonPlane->GetCoG() << endl; maxMatedChargeHyp = bestHyps[planeCtr] - i; maxMatedCharge = planeCharge; tPos[planeCtr] = muonPlane->GetInitialCoG(); } }//end if the other new hypothesis is allowed }//end if it is a muon plane }//end loop over other hypotheses bestHyps[planeCtr] = maxMatedChargeHyp; MSG("AlgDeMuxBeam", Msg::kDebug) << plane->GetPlaneNumber() << " zpos = " << zPos[planeCtr] << " tpos = " << tPos[planeCtr] << " best hyp = " << bestHyps[planeCtr] << endl; ++planeCtr; }//end loop over valid planes to adjust the fit FindFitOverNPlanes(zPos,tPos,weight,residual,slope, intercept, leverArm, planeNum); MSG("AlgDeMuxBeam", Msg::kDebug) << "set to slope = " << slope << " intercept = " << intercept << endl; fFinalShowerRegionSlope = slope; fFinalShowerRegionIntercept = intercept; //set the first n planes SetFirstNPlanes(planeItr, slope, intercept); return; }

Float_t AlgDeMuxBeam::FindDigitCoG (  DmxPlaneItr &  validPlaneItr, CandDeMuxDigitHandleItr &  cdhit )  [private]

Definition at line 1717 of file AlgDeMuxBeam.cxx. References fVldContext, UgliGeomHandle::GetStripHandle(), and UgliStripHandle::GetTPos(). { // MSG("AlgDeMuxBeam", Msg::kDebug) << "In FindDigitCoG" << endl; Float_t cog = -5.; Float_t sigTPosSum = 0.; Float_t sigSum = 0.; UgliGeomHandle ugh(fVldContext); UgliStripHandle ush; validPlaneItr.ResetFirst(); //loop over the valid planes and get the digits from those planes while(validPlaneItr.IsValid()){ cdhit.GetSet()->Slice(validPlaneItr.Ptr()->GetPlaneNumber()); //got to tell the iterator to sort the slice cdhit.GetSet()->Update(); // MSG("AlgDeMuxBeam", Msg::kDebug) << "getting CoG with plane " // << validPlaneItr.Ptr()->GetPlaneNumber() // << " digits" << endl; while(cdhit.IsValid()){ if(cdhit.Ptr()->GetDeMuxDigitFlagWord() == 0){ ush = ugh.GetStripHandle(cdhit.Ptr()->GetPlexSEIdAltL().GetBestSEId()); sigTPosSum += ush.GetTPos()*cdhit.Ptr()->GetCharge(); sigSum += cdhit.Ptr()->GetCharge(); } cdhit.Next(); }//end loop over digits //clear the digit slice cdhit.GetSet()->ClearSlice(false); //advance the plane iterator validPlaneItr.Next(); }//end loop over planes validPlaneItr.ResetFirst(); cdhit.ResetFirst(); //find the cog if(sigSum > 0.) cog = sigTPosSum/sigSum; // MSG("AlgDeMuxBeam", Msg::kDebug) << "digit cog = " << cog << endl; return cog; }

void AlgDeMuxBeam::FindFit (  DmxPlaneItr &  planeItr, DmxPlaneItr &  validPlaneItr, Double_t &  a1, Double_t &  a2, Double_t &  a3, Double_t &  a4, Double_t &  bestChi2, Int_t  leverArm, Int_t *  hypotheses )  [private]

Definition at line 2149 of file AlgDeMuxBeam.cxx. References FindFitForHypothesisSet(), FindPlanesToDropFromFit(), and DmxUtilities::IsValidFit(). Referenced by DeMuxFirstNPlanes(). { Double_t fitA1 = 0.; Double_t fitA2 = 0.; Double_t fitA3 = 0.; Double_t fitA4 = 0.; Double_t chi2 = 0.; DmxUtilities util; FindFitForHypothesisSet(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypotheses, leverArm); if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4) ){ bestChi2 = chi2; a1 = fitA1; a2 = fitA2; a3 = fitA3; a4 = fitA4; //see if you can drop some planes and improve the fit if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, leverArm) ){ FindFitForHypothesisSet(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypotheses, leverArm); if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4) ){ bestChi2 = chi2; a1 = fitA1; a2 = fitA2; a3 = fitA3; a4 = fitA4; } validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE); }//end if planes should be dropped because they are way off }//end if its a better straight line fit return; }

void AlgDeMuxBeam::FindFitForHypothesisSet (  DmxPlaneItr &  planeItr, Double_t &  prevChiSq, Double_t &  a1, Double_t &  a2, Double_t &  a3, Double_t &  a4, Int_t *  hypotheses, Int_t  leverArm )  [private]

Definition at line 2197 of file AlgDeMuxBeam.cxx. References DmxUtilities::FindLinearFit(), DmxPlane::GetCoG(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlaneTypes::GetPlaneTypeString(), DmxPlane::GetZPosition(), DmxPlane::IsGolden(), DmxPlane::IsValid(), and MSG. Referenced by FindFit(). { MSG("AlgDeMuxBeam", Msg::kDebug) << "in FindFitForHypothesisSet" << endl; planeItr.ResetFirst(); //cout << "in find fit for first n planes" << endl; //declare the variables to do a least squares fit to a straight line. //weight is the inverse fraction of the total charge on the plane that //is on opposite ends of common strips. multiply the inverse by a //sensible number to take account for digits with 1000+ adc's //the arrays are the number of planes in the set, ie the lever arm Double_t x[] = {0.,0.,0.,0.,0.,0.}; Double_t y[] = {0.,0.,0.,0.,0.,0.}; Double_t weight[] = {1.,1.,1.,1.,1.,1.}; Double_t chiSq = 0.; DmxUtilities util; Int_t planeCtr = 0; Int_t arrayCtr = 0; DmxPlane *plane = 0; //loop over the plane iterator to fill the variables while( (plane = planeItr()) && planeCtr < leverArm){ //get some generic info about the plane MSG("AlgDeMuxBeam", Msg::kDebug) << "plane " << plane->GetPlaneNumber() << " is valid " << (Int_t)plane->IsValid() << " is golden " << (Int_t)plane->IsGolden() << " masked " << (Int_t)planeItr.GetSet()->GetMasks().GetMask(plane) << " " << DmxPlaneTypes::GetPlaneTypeString(plane->GetPlaneType()) << endl; if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){ x[arrayCtr] = plane->GetZPosition(); //get the maximum possible weight first, then for each hypothesis //multiply by the square of the fraction weight[arrayCtr] = 1./ TMath::Sqrt(plane->GetPlaneCharge()); if(plane->IsGolden()) weight[arrayCtr] /= TMath::Sqrt(10.); y[arrayCtr] = plane->GetCoG(); //only do this if the shower plane is not golden - if it is, //you know where it should be reconstructed to if( plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){ if(hypotheses[planeCtr] == 1){ weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()); } //plane->GetCoG() returns the best hypothesis CoG so only change the value //of y if you are wanting the 2nd or 3rd best hypotheses else if(hypotheses[planeCtr] == 2){ y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG(); weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio()); } else if(hypotheses[planeCtr] == 3){ y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG(); weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio()); } } MSG("AlgDeMuxBeam", Msg::kDebug) << "zpos = " << x[arrayCtr] << " tpos = " << y[arrayCtr] << " hypothesis = " << hypotheses[arrayCtr] << endl; ++arrayCtr; }//end if the plane isnt marked as sucking && it is in the window bounds ++planeCtr; } //end loop over planes to find variables for fit //cout << "filled arrays" << endl; //use the DmxUtilities function to find a linear fit util.FindLinearFit(x, y, weight, arrayCtr, a1, a2, chiSq); prevChiSq = (double)chiSq; //a1 = intercept; //a2 = slope; MSG("AlgDeMuxBeam", Msg::kDebug) << "\tinitial fit " << a1 << "\t" << a2 << "\t" << a3 << "\t" << a4 << "\t" << chiSq << endl; planeItr.Reset(); planeCtr = 0; arrayCtr = 0; //cout << "reset the iterator" << endl; //redo the fit as many times as there are planes, dropping each plane in //turn to see if it produces a better chi^2 value - only do it if you have //more than 3 planes - ie you can drop one and still do a reasonable fit //2 planes = great straight line every time. if(leverArm > 3){ Double_t fitA1 = 0.; Double_t fitA2 = 0.; Double_t fitA3 = 0.; Double_t fitA4 = 0.; for(Int_t i = 0; i < leverArm; i++){ arrayCtr = 0; //cout <<"in loop " << i << endl; //loop over the plane iterator to fill the variables while( (plane = planeItr()) && planeCtr < leverArm){ //use those planes not marked as kTRUE - a mark of kTRUE means to drop the plane //from the set when finding the fit. ie if its true, the plane truely sucks if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){ if(i != planeCtr || plane->IsGolden() ){ x[arrayCtr] = plane->GetZPosition(); y[arrayCtr] = plane->GetCoG(); weight[arrayCtr] = 1. / TMath::Sqrt(plane->GetPlaneCharge()); if(plane->IsGolden()) weight[arrayCtr] /= TMath::Sqrt(10.); //only make changes if not a golden plane if( plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){ //only change y if the hypothesis isnt the best one - the unset plane //returns the best cog from the GetCoG() method if(hypotheses[planeCtr] == 1){ weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()); } else if(hypotheses[planeCtr] == 2){ y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG(); weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio()); } else if(hypotheses[planeCtr] == 3){ y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG(); weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio()); } } ++arrayCtr; //else MSG("AlgDeMuxBeam", Msg::kDebug) << "\tdrop plane " << plane->GetPlaneNumber() //<< endl; } //end if plane is not the dropped one }//end if plane is in window bounds ++planeCtr; } //end loop over planes to find variables for fit //use the DmxUtilities function to find a linear fit util.FindLinearFit(x, y, weight, arrayCtr, fitA1, fitA2, chiSq); MSG("AlgDeMuxBeam", Msg::kDebug) << "\tlinear fit drop plane " << chiSq << "\t" << fitA1 << "\t" << fitA2 << endl; planeItr.Reset(); planeCtr = 0; if(chiSq < prevChiSq){ prevChiSq = chiSq; a1 = fitA1; a2 = fitA2; a3 = fitA3; a4 = fitA4; MSG("AlgDeMuxBeam", Msg::kDebug) << "\trefined fit, drop plane " << i << "\t" << prevChiSq << "\t" << a1 << "\t" << a2 << endl; } }//end for loop to improve fit by dropping a plane }//end if enough planes to improve fit planeItr.ResetFirst(); return; }

void AlgDeMuxBeam::FindFitOverNPlanes (  Float_t *  x, Float_t *  y, Float_t *  sigma, Float_t *  residual, Float_t &  slope, Float_t &  intercept, Int_t  numPoints, Int_t *  planeNum )  [private]

Definition at line 896 of file AlgDeMuxBeam.cxx. References fMaxMuonRemovalFraction, fMaxResidual, MSG, and RegionLinearFit(). Referenced by DeMuxFirstNPlanesTest(), UseMuonSlidingWindow(), and UseShowerSlidingWindow(). { Float_t chiSqr; Float_t averageResidual = 0.; Float_t maxResidual = 0.; Float_t maxChiSqr = 0.; Float_t prevWeight = 0.; Int_t maxResidualLoc = 0; Int_t numBadResidual = 0; //now fit a straight line RegionLinearFit(slope, intercept, chiSqr, x, y, sigma, residual, numPoints); //find the fit for the first n planes and look for //planes with large residuals. repeat until planes with //large residuals are removed Int_t nRemoved = 0; bool pointRemoved = true; while(nRemoved<=fMaxMuonRemovalFraction*numPoints && pointRemoved){ pointRemoved = false; RegionLinearFit(slope, intercept, chiSqr, x, y, sigma, residual, numPoints); averageResidual = 0.; maxResidual = 0.; numBadResidual = 0; for(Int_t i = 0; i<numPoints; ++i){ MSG("AlgDeMuxBeam", Msg::kDebug) << "zpos = " << x[i] << " tpos = " << y[i] << " plane = " << planeNum[i] << " residual = " << residual[i] << endl; if(sigma[i]>0.){ averageResidual += TMath::Abs(residual[i]); if(TMath::Abs(residual[i])>fMaxResidual) ++numBadResidual; } //find the point with the largest residual more than the allowed value if(TMath::Abs(residual[i])>maxResidual){ maxResidual = TMath::Abs(residual[i]); maxResidualLoc = i; MSG("AlgDeMuxBeam", Msg::kDebug) << "biggest residual for plane " << planeNum[i] << endl; } }//end loop to find average residual and entry with the largest residual averageResidual /= (1.*(numPoints-nRemoved)); MSG("AlgDeMuxBeam", Msg::kDebug) << "average residual for first n planes = " << averageResidual << "/" << fMaxResidual << endl; //is the average residual too big? then remove the point with the //largest residual and try again. if(averageResidual>fMaxResidual || numBadResidual>fMaxBadResidual){ ++nRemoved; pointRemoved = true; maxChiSqr = chiSqr; //try removing one plane at a time from the fit and see which one causes the bad fit for(Int_t j = 0; j<numPoints; ++j){ prevWeight = sigma[j]; sigma[j] = 0.; RegionLinearFit(slope, intercept, chiSqr, x, y, sigma, residual, numPoints); sigma[j] = prevWeight; if(chiSqr<maxChiSqr){ maxChiSqr = chiSqr; maxResidualLoc = j; } }//end loop over planes to get rid of the bad one MSG("AlgDeMuxBeam", Msg::kDebug) << "!!!! removed plane " << planeNum[maxResidualLoc] << " from fit !!!!" << endl; sigma[maxResidualLoc] = 0.; } }//end loop to check fit return; }

Bool_t AlgDeMuxBeam::FindPlanesToDropFromFit (  DmxPlaneItr &  planeItr, Double_t  a1, Double_t  a2, Double_t  a3, Double_t  a4, Int_t  leverArm )  [private]

Definition at line 1773 of file AlgDeMuxBeam.cxx. References DmxPlane::GetCoG(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlane::GetZPosition(), DmxPlane::IsGolden(), MSG, and DmxPlane::SetGolden(). Referenced by FindFit(). { MSG("AlgDeMuxBeam", Msg::kDebug) << "in FindPlanesToDropFromFit" << endl; planeItr.ResetFirst(); Double_t diff1 = 0.; Double_t diff2 = 0.; Double_t diff3 = 0.; Double_t fitCoG = 0.; bool planesDropped = false; Int_t dropPlanes = 0; Int_t planeCtr = 0; //MSG("AlgDeMuxBeam", Msg::kDebug)<<"\tin FindPlanesToDropFromFit"; DmxPlane *plane = 0; //find the plane farthest off from the fit, keep going until you have just 3 planes left while( (plane = planeItr()) && dropPlanes<(Int_t)(0.1*leverArm) && planeCtr<leverArm){ fitCoG = (a1 + (plane->GetZPosition())*a2 + TMath::Power(plane->GetZPosition(),2)*a3 + TMath::Power(plane->GetZPosition(),3)*a4); //MSG("AlgDeMuxBeam", Msg::kDebug)<<"\t" << plane->GetPlaneNumber() << endl; if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){ if( plane->GetPlaneType() == DmxPlaneTypes::kShower ){ diff1 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG()); diff2 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG()); diff3 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG()); } else if( plane->GetPlaneType() == DmxPlaneTypes::kMuon ){ diff1 = TMath::Abs(fitCoG - plane->GetCoG()); diff2 = TMath::Abs(fitCoG - plane->GetCoG()); diff3 = TMath::Abs(fitCoG - plane->GetCoG()); } //if all the 3 best hypotheses are way off from the fit, that plane is most //likely messing up the fit. so mark it as kTRUE - ie it, it truely sucks if(diff1 > .504 && diff2 > .504 && diff3 > .504 && !plane->IsGolden() ){ planeItr.GetSet()->GetMasks().SetMask(kTRUE, planeItr.Ptr()); planesDropped = true; MSG("AlgDeMuxBeam", Msg::kDebug)<<"\tmarking plane " << plane->GetPlaneNumber() << endl; //decrement the number of planes now used for a fit. ++dropPlanes; } else if(diff1 > 0.504 && plane->IsGolden()){ //if this was a golden plane but it just doesnt work, make a non-golden plane planeItr.GetSet()->GetMasks().SetMask(kTRUE, planeItr.Ptr()); plane->SetGolden(false); planesDropped = true; MSG("AlgDeMuxBeam", Msg::kDebug)<<"\tmarking plane " << plane->GetPlaneNumber() << endl; //decrement the number of planes now used for a fit. ++dropPlanes; }//end if dropping a golden plane } //end if the mask hasnt already been set for this plane ++planeCtr; } planeItr.ResetFirst(); //cout << "finished drop planes from fit" << endl; return planesDropped; }

Float_t AlgDeMuxBeam::FindRegionMatedChargeFraction (  DmxPlaneItr &  validPlaneItr  )  [private]

Definition at line 1541 of file AlgDeMuxBeam.cxx. References DmxPlane::GetCoG(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlane::GetStripCoG(), and MSG. Referenced by RunAlg(). { MSG("AlgDeMuxBeam", Msg::kDebug) << "In FindRegionMatedChargeFraction" << endl; //reset the validPlaneItr validPlaneItr.ResetFirst(); DmxPlane *plane = 0; Float_t totalCharge = 0.; Float_t matedCharge = 0.; Float_t charge = 0.; while( (plane = validPlaneItr()) ){ charge = plane->GetPlaneCharge(); totalCharge += charge; MSG("AlgDeMuxBeam", Msg::kDebug) << "plane = " << plane->GetPlaneNumber() << " charge = " << charge << endl; if(plane->GetPlaneType() == DmxPlaneTypes::kShower) matedCharge += dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio()*charge; else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){ MSG("AlgDeMuxBeam", Msg::kDebug) << "muon plane cog = " << plane->GetCoG() << "/" << dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG() << "/" << dynamic_cast<DmxMuonPlane *>(plane)->GetStripCoG() << "/" << dynamic_cast<DmxMuonPlane *>(plane)->GetInitialStripCoG() << endl; //since this is a valid muon plane, all the charge is mated, assuming that the initial //center of gravity is the same as the set center of gravity if(plane->GetCoG() == dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG()) matedCharge += charge; }//end if muon plane MSG("AlgDeMuxBeam", Msg::kDebug) << " mated charge = " << matedCharge << endl; }//end loop over valid planes //get the fraction of mated charge in this region. check it against the //passed in parameter Float_t matedChargeFrac = matedCharge; if(totalCharge > 0.) matedChargeFrac /= totalCharge; else matedChargeFrac = -1.; MSG("AlgDeMuxBeam", Msg::kDebug) << "mated charge fraction in region = " << matedChargeFrac << endl; return matedChargeFrac; }

void AlgDeMuxBeam::FindRMSOfPlanes (  DmxPlaneItr &  planeItr, Int_t  leverArm, Int_t *  hypsToUse, Float_t &  cog, Float_t &  rms )  [private]

Definition at line 1847 of file AlgDeMuxBeam.cxx. References DmxPlane::GetCoG(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneType(), and MSG. { // MSG("AlgDeMuxBeam", Msg::kDebug) << "In FindRMSOfPlanes" << endl; //initialize the variables for finding the RMS Double_t sigTPos2Sum = 0.; Double_t sigSum = 0.; Double_t sigTPosSum = 0.; Double_t signal = 0.; Double_t planeCoG = 0.; Int_t ctr = 0; planeItr.ResetFirst(); //loop over the planes in the set and use only those that are //between firstPlane and lastPlane DmxPlane *plane = 0; while( planeItr.IsValid() && ctr < leverArm){ plane = planeItr.Ptr(); signal = plane->GetPlaneCharge(); planeCoG = plane->GetCoG(); if( plane->GetPlaneType() == DmxPlaneTypes::kShower ){ if( hypsToUse[ctr] == 1){ signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio(); planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG(); } else if( hypsToUse[ctr] == 2){ signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio(); planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG(); } else if( hypsToUse[ctr] == 3){ signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio(); planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG(); } } sigTPos2Sum += signal * planeCoG * planeCoG; sigTPosSum += signal * planeCoG; sigSum += signal; //increment the ctr ++ctr; //move to the next plane in the list planeItr.Next(); } planeItr.Reset(); //calculate the RMS using the formula from Bevington //[Sum(S_i * x_i^2) / Sum(S_i) - (av_x)^2] * (N / N -1) where //S_i is the signal at a distance x_i and av_x is the signal //weighted average distance from the line for a digit and N is the //number of measures in the set Double_t weighter = (leverArm * 1.0) / (1.0 * (leverArm - 1)); Double_t avTPos = sigTPosSum / sigSum; Double_t rmsSq = ((sigTPos2Sum / sigSum) - (avTPos * avTPos)) * weighter; cog = avTPos; //stupid round off error - make sure it doesnt make the rms^2 negative if(TMath::Abs(rmsSq)<1.e-3) rmsSq = 0.; rms = TMath::Sqrt(rmsSq); MSG("AlgDeMuxBeam", Msg::kDebug) << "cog = " << cog << " rms = " << rms << endl; //see what happens if you drop a plane from the set ctr = 0; //redo the fit using only the previously set planes and only if you have //more than 3 planes in the set if(leverArm > 3){ Double_t rmsRef = 0.; Double_t cogRef = 0.; for(Int_t i = 0; i < leverArm; i++){ ctr = 0; sigTPosSum = 0; sigTPos2Sum = 0; sigSum = 0; //loop over the plane iterator to fill the variables while( planeItr.IsValid() && ctr < leverArm){ plane = planeItr.Ptr(); if(i != ctr){ signal = plane->GetPlaneCharge(); planeCoG = plane->GetCoG(); if( plane->GetPlaneType() == DmxPlaneTypes::kShower ){ if( hypsToUse[ctr] == 1){ signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio(); planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG(); } else if( hypsToUse[ctr] == 2){ signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio(); planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG(); } else if( hypsToUse[ctr] == 3){ signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio(); planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG(); } } sigTPos2Sum += signal * planeCoG * planeCoG; sigTPosSum += signal * planeCoG; sigSum += signal; // MSG("AlgDeMuxBeam", Msg::kDebug) << "inloop " << sigTPos2Sum << " " << sigTPosSum // << " " << sigSum << endl; }//end if not the plane to drop ++ctr; planeItr.Next(); } //end loop over planes to find variables for fit planeItr.Reset(); Double_t weighterRef = (leverArm * 1.0) / (1.0 * (leverArm - 1)); Double_t avTPosRef = sigTPosSum / sigSum; Double_t rmsSqRef = ((sigTPos2Sum / sigSum) - (avTPosRef * avTPosRef)) * weighterRef; // MSG("AlgDeMuxBeam", Msg::kDebug) << avTPosRef << " " << sigTPos2Sum/sigSum << " " // << weighterRef << endl; //stupid round off error - make sure it doesnt make the rms^2 negative if(TMath::Abs(rmsSqRef)<1.e-3) rmsSqRef = 0.; cogRef = avTPosRef; rmsRef = TMath::Sqrt(rmsSqRef); if(rmsRef < rms){ rms = rmsRef; cog = cogRef; // MSG("AlgDeMuxBeam", Msg::kDebug) << "new cog = " << cog << " rms = " << rms << endl; } }//end loop over planes to drop one at a time }//end if at least 3 planes in set return; }

void AlgDeMuxBeam::FindShowerWindowSolution (  DmxPlaneItr &  windowPlaneItr, Float_t  cog, Int_t  firstPlane, Int_t  lastPlane )  [private]

void AlgDeMuxBeam::ReconcileShowerAndMuonRegions (  DmxPlaneItr &  planeItr, CandDeMuxDigitHandleItr &  cdhit )  [private]

Definition at line 1600 of file AlgDeMuxBeam.cxx. References PlexSEIdAltL::ClearWeights(), digit(), fDigitResetCut, fFinalShowerRegionSlope, fInitialMuonRegionSlope, fMuonStartPlaneNumber, fShowerMuonTransverseDifference, fVldContext, PlexSEIdAltL::GetBestSEId(), DmxPlane::GetCoG(), PlexSEIdAltL::GetCurrentSEId(), CandDeMuxDigitHandle::GetDeMuxDigitFlagWord(), DmxPlane::GetPlaneNumber(), CandDigitHandle::GetPlexSEIdAltL(), CandDigitHandle::GetPlexSEIdAltLWritable(), UgliGeomHandle::GetStripHandle(), UgliStripHandle::GetTPos(), DmxPlane::GetZPosition(), PlexSEIdAltL::IsValid(), MSG, PlexSEIdAltL::Next(), PlexSEIdAltL::SetCurrentWeight(), and PlexSEIdAltL::SetFirst(). Referenced by RunAlg(). { //get the last slope and intercept for the shower and muon regions and //compare the predicted tpos at the muon start plane number. //also find the extent of the shower at that point. if the muon track //is more than the allowed value away from the shower, then reconcile //the shower to the track. Float_t showerPredictedTPos = fFinalShowerRegionSlope*fMuonStartPlaneZPos; showerPredictedTPos += fFinalShowerRegionIntercept; Float_t muonPredictedTPos = fInitialMuonRegionSlope*fMuonStartPlaneZPos; muonPredictedTPos += fInitialMuonRegionIntercept; MSG("AlgDeMuxBeam", Msg::kDebug) << "shower prediction = " << showerPredictedTPos << " muon prediction = " << muonPredictedTPos << " muon start zpos = " << fMuonStartPlaneZPos << endl; CandDeMuxDigitHandle *digit = 0; DmxPlane *plane = 0; UgliGeomHandle ugh(fVldContext); UgliStripHandle ush; Float_t maxPlaneTPos = -5.; Float_t minPlaneTPos = 5.; Float_t digitTPos = 0.; Float_t setDigitTPos = 0.; if(TMath::Abs(muonPredictedTPos-showerPredictedTPos) > fShowerMuonTransverseDifference){ //it looks like the muon and shower don't match up. loop over the //CandDeMuxDigitList and see what the extent of the shower region is while( (plane = planeItr()) ){ maxPlaneTPos = -5.; minPlaneTPos = 5.; cdhit.GetSet()->Slice(plane->GetPlaneNumber()); //got to tell the iterator to sort the slice cdhit.GetSet()->Update(); while( (digit = cdhit()) ){ if(digit->GetDeMuxDigitFlagWord() == 0){ ush = ugh.GetStripHandle(digit->GetPlexSEIdAltL().GetBestSEId()); digitTPos = ush.GetTPos(); if(digitTPos > maxPlaneTPos) maxPlaneTPos = digitTPos; else if(digitTPos < minPlaneTPos) minPlaneTPos = digitTPos; } }//end loop over digits in the plane //reset the digit iterator cdhit.ResetFirst(); //check the demuxed cog of each plane in the shower region against //the predicted position based on the initial muon linear fit. also //check that the predicted position is outside the minimum and maximum //transverse positions of the reconstructed strips in the the plane. //if the difference is more than the allowed value, or the predicted //location is outside the reconstructed strip bounds, then redo the reconstruction muonPredictedTPos = fInitialMuonRegionSlope*plane->GetZPosition(); muonPredictedTPos += fInitialMuonRegionIntercept; if(plane->GetPlaneNumber() < fMuonStartPlaneNumber && TMath::Abs(plane->GetCoG()-muonPredictedTPos) > fShowerMuonTransverseDifference && (muonPredictedTPos<minPlaneTPos || muonPredictedTPos>maxPlaneTPos)){ //plane->SetStrips(muonPredictedTPos); //loop over the digits again and look at the SEIdAltLists - if there is //an alternative within the accepted bound of the predicted position, //take it. while( (digit = cdhit()) ){ ush = ugh.GetStripHandle(digit->GetPlexSEIdAltL().GetBestSEId()); setDigitTPos = ush.GetTPos(); //set the alternatives list to the first one digit->GetPlexSEIdAltL().SetFirst(); while( digit->GetPlexSEIdAltL().IsValid() ){ ush = ugh.GetStripHandle(digit->GetPlexSEIdAltL().GetCurrentSEId()); digitTPos = ush.GetTPos(); MSG("AlgDeMuxBeam", Msg::kDebug) << "current tpos = " << digitTPos << "/" << muonPredictedTPos << endl; if(TMath::Abs(digitTPos-muonPredictedTPos) < fDigitResetCut || (TMath::Abs(setDigitTPos-muonPredictedTPos) > fShowerMuonTransverseDifference && TMath::Abs(digitTPos-muonPredictedTPos) < fShowerMuonTransverseDifference) ){ //clear all weights for the digit digit->GetPlexSEIdAltLWritable().ClearWeights(); MSG("AlgDeMuxBeam", Msg::kDebug) << "setting" << endl; digit->GetPlexSEIdAltLWritable().SetCurrentWeight(1.0); break; } digit->GetPlexSEIdAltL().Next(); }//end loop over alt list }//end loop over digits to find the ones near the predicted tpos }//end if this plane should have some digits reset //clear the digit slice cdhit.GetSet()->ClearSlice(false); }//end loop over planeItr }//end if the shower and muon regions have different predicted transverse positions return; }

void AlgDeMuxBeam::RegionLinearFit (  Float_t &  slope, Float_t &  intercept, Float_t &  chiSqr, Float_t *  x, Float_t *  y, Float_t *  sigma, Float_t *  residual, Int_t  numPoints )  [private]

Definition at line 1476 of file AlgDeMuxBeam.cxx. References MSG. Referenced by FindFitOverNPlanes(). { //define the necessary sums Double_t sumXSqrDivSigmaSqr = 0.; Double_t sumYDivSigmaSqr = 0.; Double_t sumXDivSigmaSqr = 0.; Double_t sumXYDivSigmaSqr = 0.; Double_t sum1DivSigmaSqr = 0.; Double_t sigmaSqr = 0.; chiSqr = 0.; for(Int_t i=0; i<numPoints; ++i){ if(TMath::Abs(sigma[i])>0.){ sigmaSqr = sigma[i]*sigma[i]; sum1DivSigmaSqr += 1./sigmaSqr; sumXDivSigmaSqr += x[i]/sigmaSqr; sumYDivSigmaSqr += y[i]/sigmaSqr; sumXYDivSigmaSqr += (x[i]*y[i])/sigmaSqr; sumXSqrDivSigmaSqr += (x[i]*x[i])/sigmaSqr; } MSG("AlgDeMuxBeam", Msg::kDebug) << "x = " << x[i] << " y = " << y[i] << " sigma = " << sigma[i] << " sum1DivSigmaSqr = " << sum1DivSigmaSqr << " sumXDivSigmaSqr = " << sumXDivSigmaSqr << " sumYDivSigmaSqr = " << sumYDivSigmaSqr << " sumXYDivSigmaSqr = " << sumXYDivSigmaSqr << " sumXSqrDivSigmaSqr = " << sumXSqrDivSigmaSqr << endl; } Float_t delta = sum1DivSigmaSqr*sumXSqrDivSigmaSqr - sumXDivSigmaSqr*sumXDivSigmaSqr; if(delta==0.){ MSG("AlgDeMuxBeam", Msg::kWarning) << "delta = 0 in linear fit - kludge an answer and bail" << endl; intercept = 0.; slope = 0.; chiSqr = 10000.; return; } intercept = (sumXSqrDivSigmaSqr*sumYDivSigmaSqr - sumXDivSigmaSqr*sumXYDivSigmaSqr)/delta; slope = (sum1DivSigmaSqr*sumXYDivSigmaSqr - sumXDivSigmaSqr*sumYDivSigmaSqr)/delta; //now get the chi^2 and residual for the fit - zero chiSqr to be safe chiSqr = 0.; Double_t numerator = 0.; Double_t denominator = 0.; for(Int_t i=0; i<numPoints; ++i){ residuals[i] = (y[i]-intercept-slope*x[i]); numerator = TMath::Power(y[i]-intercept-slope*x[i],2.); denominator = sigma[i]*sigma[i]; if(denominator>0.) chiSqr += numerator/denominator; } //fit to a straight line has 2 degrees of freedom if(numPoints>2) chiSqr /= 1.*(numPoints-2); return; }

void AlgDeMuxBeam::RunAlg (  AlgConfig &  acd, CandHandle &  ch, CandContext &  cx )  [virtual]

Implements AlgBase. Definition at line 91 of file AlgDeMuxBeam.cxx. References Detector::AsString(), DmxUtilities::CheckFitWithTiming(), DeMuxFirstNPlanesTest(), fDigitResetCut, DmxUtilities::FillPlaneArray(), DmxUtilities::FindBeamMuonStartPlane(), DmxUtilities::FindEndPlane(), DmxUtilities::FindPlanesOffFit(), FindRegionMatedChargeFraction(), DmxUtilities::FindVertexPlane(), fMaxBadResidual, fMaxMuonRemovalFraction, fMaxResidual, fMaxStripAdjustment, fMinMatedChargeFrac, fMinMuonPlanesRequired, fMuonPlanesInSet, fMuonSlopeLimit, fMuonStartPlaneNumber, fMuonStartPlaneZPos, fNumberOfHypotheses, fPlanesInSet, fShowerMuonTransverseDifference, fStrayCut, fStrayPlanes, fVldContext, DmxStatus::GetAverageTimingOffset(), CandRecord::GetCandHeader(), CandHandle::GetCandRecord(), CandHandle::GetDaughterIterator(), VldContext::GetDetector(), Registry::GetDouble(), DmxStatus::GetEndPlaneNumber(), DmxStatus::GetEventDeMuxed(), DmxStatus::GetEventNumber(), Registry::GetInt(), DmxStatus::GetMuonStartPlaneNumber(), DmxStatus::GetNonPhysicalFailure(), DmxStatus::GetPlaneArray(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlaneTypes::GetPlaneTypeString(), DmxPlane::GetPlaneView(), CandHeader::GetSnarl(), DmxStatus::GetVertexPlaneNumber(), CandHandle::GetVldContext(), DmxPlane::GetZPosition(), DmxPlane::IsValid(), KeyFromPlane1(), KeyOnPlane1(), KeyOnUView1(), KeyOnValidU1(), KeyOnValidV1(), KeyOnVView1(), MSG, ReconcileShowerAndMuonRegions(), DmxStatus::ResetStatus(), DmxStatus::SetAverageTimingOffset(), CandDeMuxDigitListHandle::SetAvgTimeOffset(), CandDeMuxDigitListHandle::SetDeMuxDigitListFlagBit(), DmxStatus::SetEndPlaneNumber(), DmxStatus::SetEventNumber(), DmxStatus::SetFigureOfMerit(), DmxStatus::SetMultipleMuon(), DmxStatus::SetMuonStartPlaneNumber(), DmxStatus::SetNoVertexFailure(), DmxStatus::SetNumberOfPlanes(), CandDeMuxDigitListHandle::SetNumStrayPlanesU(), CandDeMuxDigitListHandle::SetNumStrayPlanesV(), CandDeMuxDigitListHandle::SetNumValidPlanesU(), CandDeMuxDigitListHandle::SetNumValidPlanesV(), DmxStatus::SetUMuonMatedFraction(), DmxStatus::SetUShowerMatedFraction(), DmxStatus::SetUStrayPlanes(), DmxStatus::SetUValidPlanes(), DmxStatus::SetValidPlanesFailure(), DmxStatus::SetVertexPlaneNumber(), DmxStatus::SetVertexPlaneZPosition(), DmxStatus::SetVMuonMatedFraction(), DmxStatus::SetVShowerMatedFraction(), DmxStatus::SetVStrayPlanes(), DmxStatus::SetVValidPlanes(), UseMuonSlidingWindow(), and UseShowerSlidingWindow(). { assert( ch.InheritsFrom("CandDigitListHandle") ); CandDeMuxDigitListHandle &cdlh = dynamic_cast<CandDeMuxDigitListHandle&>(ch); CandDeMuxDigitHandleItr cdhit = cdlh.GetDaughterIterator(); fVldContext = *(ch.GetVldContext()); if (fVldContext.GetDetector() != Detector::kFar) { static int msglimit = 20; // no more than 20 warnings about the detector if (msglimit) { MSG("AlgDeMuxBeam",Msg::kDebug) << "AlgDeMuxBeam::RunAlg() can not DeMux " << Detector::AsString(fVldContext.GetDetector()) << " detector. Skip futher DeMux processing." << endl; if (--msglimit == 0) MSG("AlgDeMuxBeam",Msg::kDebug) << " ... last message of this type." << endl; } return; } //get the AlgConfigDeMux object fPlanesInSet = acd.GetInt("PlanesInSet"); fMinMuonPlanesRequired = acd.GetInt("MinMuonPlanesRequired"); fMinMatedChargeFrac = acd.GetDouble("MinMatedChargeFraction"); fMuonPlanesInSet = acd.GetInt("MuonPlanesInSet"); fMuonSlopeLimit = acd.GetDouble("MuonTrackSlopeLimit"); fMaxResidual = acd.GetDouble("MaxResidual"); fMaxBadResidual = acd.GetInt("MaxBadResidual"); fMaxMuonRemovalFraction = acd.GetDouble("MaxMuonRemovalFraction"); fMaxStripAdjustment = acd.GetInt("MaxStripAdjustment"); fStrayCut = acd.GetInt("StrayDeltaStripLimit"); fNumberOfHypotheses = acd.GetInt("NumberOfHypotheses"); fShowerMuonTransverseDifference = acd.GetDouble("ShowerMuonTransverseDifference"); fDigitResetCut = acd.GetDouble("ResetDigitCut"); //MSG("AlgDeMuxBeam", Msg::kDebug) << "window size = " << fPlanesInSet << endl; //get the DmxStatus object // Find the DmxStatus object or create one for needed scratch space DmxStatus *status = dynamic_cast<DmxStatus *>(gROOT->GetRootFolder()->FindObject("Loon/DeMux/DmxStatus")); bool tempStatus = false; if (status == 0) { MSG("AlgDeMuxBeam", Msg::kDebug) << "//root/Loon/DeMux/DmxStatus not found." << " Create a temporary DmxStatus." << endl; status = new DmxStatus; // Make a temporary DmxStatus if needed tempStatus = true; } else status->ResetStatus(); status->SetEventNumber(ch.GetCandRecord()->GetCandHeader()->GetSnarl()); MSG("AlgDeMuxBeam", Msg::kDebug) << "starting beam demuxing for event " << status->GetEventNumber() << endl; DmxUtilities util; util.FillPlaneArray(status, cdlh, acd); const TObjArray *planeArray = status->GetPlaneArray(); //DmxPlane *plane = dynamic_cast<DmxPlane *>(planeArray->First()); //MSG("AlgDeMuxBeam", Msg::kDebug) <<"Event = " << status->GetEventNumber() << endl; //create the DmxPlaneItr over planes and program it to sort the planes DmxPlaneItr planeItr(planeArray); //create a KeyFunc to sort planes by number DmxPlaneKeyFunc *pnKF = planeItr.CreateKeyFunc(); pnKF->SetFun(KeyOnPlane1); planeItr.GetSet()->AdoptSortKeyFunc(pnKF); pnKF = 0; planeItr.ResetFirst(); status->SetNumberOfPlanes(planeItr.SizeSelect()); status->SetVertexPlaneNumber(util.FindVertexPlane(planeItr)); //demux the event if there is an identified vertex if( status->GetVertexPlaneNumber() != -1){ planeItr.ResetFirst(); planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), 500); planeItr.GetSet()->Update(); //set the end plane number since there was a vertex plane status->SetEndPlaneNumber(util.FindEndPlane(planeItr)); planeItr.GetSet()->ClearSlice(false); planeItr.ResetFirst(); //set the z position of the vertex plane planeItr.GetSet()->Slice(status->GetVertexPlaneNumber()); planeItr.GetSet()->Update(); planeItr.ResetFirst(); status->SetVertexPlaneZPosition(dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetZPosition()); planeItr.GetSet()->ClearSlice(false); planeItr.ResetFirst(); //set the plane number for the muon start plane planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber()); planeItr.GetSet()->Update(); //set the muon start plane number status->SetMuonStartPlaneNumber(util.FindBeamMuonStartPlane(planeItr, acd.GetDouble("AveragePEGainConversion"))); fMuonStartPlaneNumber = status->GetMuonStartPlaneNumber(); planeItr.GetSet()->ClearSlice(false); planeItr.ResetFirst(); //create a KeyFunc to sort digits by plane CandDeMuxDigitHandleKeyFunc *planeKF = cdhit.CreateKeyFunc(); //program the KeyFunc with the sort function planeKF->SetFun(KeyFromPlane1); //get the NavSet from the iterator and pass the KeyFunc to it cdhit.GetSet()->AdoptSortKeyFunc(planeKF); //clear the KF pointer because we no longer own the KeyFunc planeKF = 0; //make the iterators for the window, valid planes, and vertex DmxPlaneItr vertexPlaneItr(planeArray); DmxPlaneItr validPlaneItr(planeArray); //create a KeyFunc to sort planes by number DmxPlaneKeyFunc *vertexPlaneKF = vertexPlaneItr.CreateKeyFunc(); DmxPlaneKeyFunc *validPlaneKF = validPlaneItr.CreateKeyFunc(); //program the KeyFunc with the sort function vertexPlaneKF->SetFun(KeyOnPlane1); validPlaneKF->SetFun(KeyOnPlane1); //get the NavSet from the iterator and pass the KeyFunc to it vertexPlaneItr.GetSet()->AdoptSortKeyFunc(vertexPlaneKF); validPlaneItr.GetSet()->AdoptSortKeyFunc(validPlaneKF); //clear the KF pointer because we no longer own the KeyFunc vertexPlaneKF = 0; validPlaneKF = 0; planeItr.ResetFirst(); validPlaneItr.ResetFirst(); //now program a key function to select on plane orientation - U first DmxPlaneKeyFunc *orientValidUKF = validPlaneItr.CreateKeyFunc(); DmxPlaneKeyFunc *orientUKF = planeItr.CreateKeyFunc(); DmxPlaneKeyFunc *orientValidVKF = validPlaneItr.CreateKeyFunc(); DmxPlaneKeyFunc *orientVKF = planeItr.CreateKeyFunc(); //initialize the muon and shower region mated fraction variables in the DmxStatus object status->SetUShowerMatedFraction(-1.); status->SetUMuonMatedFraction(-1.); status->SetVShowerMatedFraction(-1.); status->SetVMuonMatedFraction(-1.); DmxPlane *plane = 0; bool endPlaneValid = false; planeItr.GetSet()->Update(); planeItr.ResetFirst(); //see if it worked while( (plane = planeItr()) ){ MSG("AlgDeMuxBeam", Msg::kDebug) << "number = " << plane->GetPlaneNumber() << " orientation = " << (Int_t)plane->GetPlaneView() << " type = " << DmxPlaneTypes::GetPlaneTypeString(plane->GetPlaneType()) << " IsValid() = " << (Int_t)plane->IsValid() << endl; if(plane->GetPlaneNumber() == fMuonStartPlaneNumber) fMuonStartPlaneZPos = plane->GetZPosition(); if(plane->GetPlaneNumber() == status->GetEndPlaneNumber() && plane->IsValid()) endPlaneValid = true; } planeItr.ResetFirst(); Int_t viewCtr = 0; Int_t numValidMuonPlanes = 0; Float_t muonMatedFraction = -1.; Float_t showerMatedFraction = -1.; while( viewCtr < 2 ){ //set the muon start plane number at the start of the loop for each //view as one view may have reset it to -1 if there were not enough valid //muon planes in that view after the muon start plane number fMuonStartPlaneNumber = status->GetMuonStartPlaneNumber(); muonMatedFraction = -1.; showerMatedFraction = -1.; MSG("AlgDeMuxBeam", Msg::kDebug) << "vertex plane = " << status->GetVertexPlaneNumber() << " muon start plane = " << fMuonStartPlaneNumber << " end plane = " << status->GetEndPlaneNumber() << endl; if( viewCtr == 0){ //program this function with the select function orientUKF->SetFun(KeyOnUView1); orientValidUKF->SetFun(KeyOnValidU1); //adopt it as a selection function planeItr.GetSet()->AdoptSelectKeyFunc(orientUKF); orientUKF = 0; planeItr.Reset(); validPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidUKF); orientValidUKF = 0; validPlaneItr.Reset(); } else if(viewCtr == 1){ //program this function with the select function orientVKF->SetFun(KeyOnVView1); orientValidVKF->SetFun(KeyOnValidV1); //adopt it as a selection function planeItr.GetSet()->AdoptSelectKeyFunc(orientVKF); orientVKF = 0; planeItr.Reset(); validPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidVKF); orientValidVKF = 0; validPlaneItr.Reset(); } planeItr.ResetFirst(); //there must be at least 2 valid planes for the demuxing to work. //slice the validPlaneItr to all planes in the event validPlaneItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber()); validPlaneItr.GetSet()->Update(); if(validPlaneItr.SizeSelect() >= 2){ planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber()+1); planeItr.GetSet()->Update(); //demux the first n planes in the detector. the method adjusts the size of n //such that if there are < fPlanesInSet planes hit, n becomes the number of planes //hit. Otherwise, n is fPlanesInSet if((Int_t)validPlaneItr.SizeSelect() < fPlanesInSet) DeMuxFirstNPlanesTest(planeItr, validPlaneItr); else if( (Int_t)validPlaneItr.SizeSelect() >= fPlanesInSet){ //still have to demux the first n planes of the event the same for //muon or shower like events. loop over the planes to find the last //plane in the first set of n Int_t ctr = 0; Int_t lastPlane = 0; Int_t firstPlane = 0; while( validPlaneItr.IsValid() && ctr < fPlanesInSet){ plane = validPlaneItr.Ptr(); //MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber() // << endl; if( ctr == 0 ) firstPlane = plane->GetPlaneNumber(); lastPlane = plane->GetPlaneNumber(); validPlaneItr.Next(); ++ctr; } MSG("AlgDeMuxBeam", Msg::kDebug) << "\tfirst plane = " << firstPlane << "\tlast plane = " << lastPlane << endl; MSG("AlgDeMuxBeam", Msg::kDebug) << "demuxing first N planes" << endl; //clear the slice of the validPlaneItr validPlaneItr.GetSet()->ClearSlice(false); planeItr.GetSet()->ClearSlice(false); //set the initial window validPlaneItr.GetSet()->Slice(status->GetVertexPlaneNumber(), lastPlane+1); validPlaneItr.GetSet()->Update(); //check to see if the number of valid planes == fPlanesInSet. if so, then //it means that there wont be any sliding windows, so set all planes in //the event if(status->GetEndPlaneNumber()-lastPlane <= 2 && !endPlaneValid) lastPlane = status->GetEndPlaneNumber(); //set the plane iterator to be from the vertex to the last plane in the window, //since this is the first window done planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), lastPlane+1); planeItr.GetSet()->Update(); DeMuxFirstNPlanesTest(planeItr, validPlaneItr); validPlaneItr.GetSet()->ClearSlice(false); planeItr.GetSet()->ClearSlice(false); //check to make sure that there are valid planes after the //muon start plane validPlaneItr.GetSet()->Slice(firstPlane, status->GetEndPlaneNumber()); validPlaneItr.GetSet()->Update(); numValidMuonPlanes = 0; while( (plane = validPlaneItr() ) ){ if(plane->GetPlaneNumber()>fMuonStartPlaneNumber && fMuonStartPlaneNumber > 0 && plane->GetPlaneType()==DmxPlaneTypes::kMuon) ++numValidMuonPlanes; } validPlaneItr.GetSet()->ClearSlice(false); MSG("AlgDeMuxBeam", Msg::kDebug) << "number valid muon planes in view " << viewCtr << " = " << numValidMuonPlanes << endl; //if the number of valid muon planes is less than 2, then set the //fMuonStartPlaneNumber to -1 since you will just demux the muon stuff //with the shower region if(numValidMuonPlanes<fMinMuonPlanesRequired) fMuonStartPlaneNumber = -1; //demux the shower region if it exists and it extends beyond the //first set of planes demuxed if(status->GetVertexPlaneNumber() != fMuonStartPlaneNumber && fMuonStartPlaneNumber > 0 && TMath::Abs(fMuonStartPlaneNumber-lastPlane)>2){ //slice the valid plane iterator to do the shower region MSG("AlgDeMuxBeam", Msg::kDebug) << "slice shower region with track from " << firstPlane+1 << " to " << fMuonStartPlaneNumber << endl; validPlaneItr.GetSet()->Slice(firstPlane+1, fMuonStartPlaneNumber); validPlaneItr.GetSet()->Update(); UseShowerSlidingWindow(planeItr, validPlaneItr, status); showerMatedFraction = FindRegionMatedChargeFraction(validPlaneItr); if(viewCtr == 0) status->SetUShowerMatedFraction(showerMatedFraction); if(viewCtr == 1) status->SetVShowerMatedFraction(showerMatedFraction); validPlaneItr.GetSet()->ClearSlice(false); MSG("AlgDeMuxBeam", Msg::kDebug) <<"Done Setting Shower Region for snarl " << status->GetEventNumber() << endl; } else if(fMuonStartPlaneNumber < 0 && TMath::Abs(status->GetEndPlaneNumber()-lastPlane)>2){ //slice the valid plane iterator to do the shower region - there is //no track to speak of, so just slice the planes out to the end of //the detector validPlaneItr.GetSet()->Slice(firstPlane+1, 486); validPlaneItr.GetSet()->Update(); UseShowerSlidingWindow(planeItr, validPlaneItr, status); showerMatedFraction = FindRegionMatedChargeFraction(validPlaneItr); if(viewCtr == 0) status->SetUShowerMatedFraction(showerMatedFraction); if(viewCtr == 1) status->SetVShowerMatedFraction(showerMatedFraction); validPlaneItr.GetSet()->ClearSlice(false); MSG("AlgDeMuxBeam", Msg::kDebug) <<"Done Setting Shower Region for no track or small track " << status->GetEventNumber() << endl; } //now do the muon region if it exists //if( fMuonStartPlaneNumber == status->GetVertexPlaneNumber() // && numValidMuonPlanes>1){ // validPlaneItr.GetSet()->Slice(firstPlane+1, // status->GetEndPlaneNumber()); // validPlaneItr.GetSet()->Update(); // UseMuonSlidingWindow(planeItr,validPlaneItr, status); // validPlaneItr.GetSet()->ClearSlice(false); // MSG("AlgDeMuxBeam", Msg::kDebug) << "done setting the muon for snarl " // << status->GetEventNumber() << endl; //} if(fMuonStartPlaneNumber > 0 && numValidMuonPlanes>1){ validPlaneItr.GetSet()->Slice(fMuonStartPlaneNumber, status->GetEndPlaneNumber()); validPlaneItr.GetSet()->Update(); UseMuonSlidingWindow(planeItr,validPlaneItr, status); muonMatedFraction = FindRegionMatedChargeFraction(validPlaneItr); if(viewCtr == 0) status->SetUMuonMatedFraction(muonMatedFraction); if(viewCtr == 1) status->SetVMuonMatedFraction(muonMatedFraction); validPlaneItr.GetSet()->ClearSlice(false); MSG("AlgDeMuxBeam", Msg::kDebug) << "done setting the muon for snarl " << status->GetEventNumber() << endl; } //check and see that both tracklike and showerlike regions, //that the fraction of the charge in the muon region that is mated //is above the threshold and that there are a reasonable number of //planes in the muon region for the fit //if so, see that the track and shower solutions are consistent planeItr.GetSet()->Update(); planeItr.ResetFirst(); if(status->GetVertexPlaneNumber() != fMuonStartPlaneNumber && muonMatedFraction > fMinMatedChargeFrac && numValidMuonPlanes >= fPlanesInSet) ReconcileShowerAndMuonRegions(planeItr, cdhit); //MSG("AlgDeMuxBeam", Msg::kDebug)<<"Done with enough planes to demux"<< endl; }//end if at least fPlanesInSet valid planes to reconstruct } //end if >=2 valid planes to reconstruct the event else{ MSG("AlgDeMuxBeam", Msg::kDebug)<< "Event " << status->GetEventNumber() << " not demuxed; not enough valid planes" << endl; status->SetValidPlanesFailure(true); } //what if there is more to the track? perhaps the muon went through the //coil hole and disappeared for several planes? treat it in a way similar //to multiple cosmics planeItr.GetSet()->ClearSlice(false); validPlaneItr.GetSet()->ClearSlice(false); vertexPlaneItr.GetSet()->ClearSlice(false); //MSG("AlgDeMuxBeam", Msg::kDebug) << "slice vertex itr" << endl; vertexPlaneItr.GetSet()->Slice(status->GetEndPlaneNumber()+1, 500); vertexPlaneItr.GetSet()->Update(); //MSG("AlgDeMuxBeam", Msg::kDebug)<< "event = " << status->GetEventNumber() // << "\tend plane = " << status->GetEndPlaneNumber() << endl; //MSG("AlgDeMuxBeam", Msg::kDebug) << "check for new vertex" << endl; Int_t nextVertex = -1; if(vertexPlaneItr.SizeSelect() > 3) nextVertex = util.FindVertexPlane(vertexPlaneItr); vertexPlaneItr.GetSet()->ClearSlice(false); while( nextVertex != -1){ status->SetMultipleMuon(true); vertexPlaneItr.GetSet()->ClearSlice(false); vertexPlaneItr.GetSet()->Slice(nextVertex, 500); vertexPlaneItr.GetSet()->Update(); //get the next end plane Int_t endPlane = util.FindEndPlane(vertexPlaneItr); status->SetEndPlaneNumber(endPlane); MSG("AlgDeMuxBeam", Msg::kDebug)<< "\tin gap loop, vertex at plane " << nextVertex << "\tend plane = " << endPlane << "\tplanes in slice = " << vertexPlaneItr.SizeSelect() << endl; //if we are in this loop, then it is most likely because a muon went through //the coil hole. so the thing to do is just treat all the remaining planes as //being in the muon track and only demux that way validPlaneItr.GetSet()->Slice(nextVertex, endPlane); validPlaneItr.GetSet()->Update(); validPlaneItr.ResetFirst(); //loop over the valid planes and count the number of muon planes numValidMuonPlanes = 0; while( (plane = validPlaneItr()) ){ MSG("AlgDeMuxBeam", Msg::kDebug) << "number = " << plane->GetPlaneNumber() << " orientation = " << (Int_t)plane->GetPlaneView() << " type = " << DmxPlaneTypes::GetPlaneTypeString(plane->GetPlaneType()) << " IsValid() = " << (Int_t)plane->IsValid() << endl; if(plane->GetPlaneType()==DmxPlaneTypes::kMuon) ++numValidMuonPlanes; } validPlaneItr.ResetFirst(); if(numValidMuonPlanes<2){ nextVertex = -1; continue; } UseMuonSlidingWindow(planeItr,validPlaneItr, status); validPlaneItr.GetSet()->ClearSlice(false); MSG("AlgDeMuxBeam", Msg::kDebug) << "done setting the muon for snarl " << status->GetEventNumber() << endl; //look for another vertex after the end plane planeItr.GetSet()->ClearSlice(false); validPlaneItr.GetSet()->ClearSlice(false); vertexPlaneItr.GetSet()->ClearSlice(false); vertexPlaneItr.GetSet()->Slice(endPlane+1, 500); vertexPlaneItr.GetSet()->Update(); //MSG("AlgDeMuxBeam", Msg::kDebug)<< "\tlook for new vertex, slice from " << endPlane+1 // << "\tto plane 500, planes in slice = " // << vertexPlaneItr.SizeSelect() << endl; vertexPlaneItr.Reset(); if(vertexPlaneItr.SizeSelect() > 0) nextVertex = util.FindVertexPlane(vertexPlaneItr); else nextVertex = -1; }//end loop to find next vertex for spread out multiples fStrayPlanes = util.FindPlanesOffFit(planeItr, fStrayCut); MSG("AlgDeMuxBeam", Msg::kDebug) << status->GetEventNumber() << " " << fStrayPlanes << endl; status->SetFigureOfMerit(validPlaneItr.SizeSelect(), fStrayPlanes); if(viewCtr == 0 && status->GetEventDeMuxed() ){ status->SetUStrayPlanes(fStrayPlanes); status->SetUValidPlanes(validPlaneItr.SizeSelect()); cdlh.SetNumValidPlanesU(validPlaneItr.SizeSelect()); cdlh.SetNumStrayPlanesU(fStrayPlanes); } if(viewCtr == 1 && status->GetEventDeMuxed() ){ status->SetVStrayPlanes(fStrayPlanes); status->SetVValidPlanes(validPlaneItr.SizeSelect()); cdlh.SetNumValidPlanesV(validPlaneItr.SizeSelect()); cdlh.SetNumStrayPlanesV(fStrayPlanes); } //done using the first view, so clear the selection function //MSG("AlgDeMuxBeam", Msg::kDebug) <<"Done with View " << viewCtr << endl; planeItr.GetSet()->ClearSlice(false); planeItr.GetSet()->AdoptSelectKeyFunc(0); planeItr.ResetFirst(); validPlaneItr.GetSet()->ClearSlice(false); validPlaneItr.GetSet()->AdoptSelectKeyFunc(0); validPlaneItr.Reset(); vertexPlaneItr.GetSet()->ClearSlice(false); vertexPlaneItr.Reset(); viewCtr++; } //end loop over views //if the event was demuxed, check to see how the demuxing and timing jive //otherwise see if you need to set the nonphysical solution flag if( status->GetEventDeMuxed() ){ status->SetAverageTimingOffset(util.CheckFitWithTiming(validPlaneItr)); cdlh.SetAvgTimeOffset(status->GetAverageTimingOffset()); } else if(status->GetNonPhysicalFailure() ) cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kNonPhysicalStripSolution); }//end if there was a vertex else{ status->SetNoVertexFailure(true); cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kNoVertex); MSG("AlgDeMuxBeam", Msg::kDebug) << "no vertex found for event " << status->GetEventNumber() << endl; } //get rid of the temporary status object if you made it. if(tempStatus) delete status; return; }

void AlgDeMuxBeam::SetFirstNPlanes (  DmxPlaneItr &  planeItr, Float_t  slope, Float_t  intercept )  [private]

Definition at line 861 of file AlgDeMuxBeam.cxx. References DmxPlane::GetCoG(), DmxPlane::GetPlaneNumber(), DmxPlane::GetZPosition(), MSG, and DmxPlane::SetStrips(). Referenced by DeMuxFirstNPlanes(), and DeMuxFirstNPlanesTest(). { MSG("AlgDeMuxBeam", Msg::kDebug) << "In SetFirstNPlanes " << slope << " " << intercept << endl; Float_t cog = 0.; planeItr.ResetFirst(); DmxPlane *currentPlane = 0; while( (currentPlane = planeItr()) ){ MSG("AlgDeMuxBeam", Msg::kDebug) << "on plane " << currentPlane->GetPlaneNumber() << endl; cog = slope*currentPlane->GetZPosition() + intercept; //set those strips if the plane isnt a valid muon plane //if(currentPlane->GetPlaneType() == DmxPlaneTypes::kMuon // || !currentPlane->IsValid() ){ currentPlane->SetStrips(cog); MSG("AlgDeMuxBeam", Msg::kDebug) << "setting first n planes " << currentPlane->GetPlaneNumber() << " to " << cog << "/" << currentPlane->GetCoG() << endl; //} }//end loop over planes planeItr.ResetFirst(); return; }

void AlgDeMuxBeam::SetMuonRegionWindow (  DmxPlaneItr &  planeItr, Float_t  slope, Float_t  intercept, bool  setAll )  [private]

Definition at line 1435 of file AlgDeMuxBeam.cxx. References fMuonStartPlaneNumber, DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlane::GetZPosition(), DmxPlane::IsValid(), MSG, and DmxPlane::SetStrips(). Referenced by UseMuonSlidingWindow(). { MSG("AlgDeMuxBeam", Msg::kDebug) << "In SetMuonRegionWindow" << endl; Float_t cog = 0.; planeItr.ResetFirst(); DmxPlane *currentPlane = 0; while( planeItr.IsValid() ){ //get the current plane currentPlane = planeItr.Ptr(); MSG("AlgDeMuxBeam", Msg::kDebug) << "on plane " << currentPlane->GetPlaneNumber() << endl; cog = slope*currentPlane->GetZPosition() + intercept; //set those strips if the plane isnt a valid muon plane in the //muon track or it has been determined that the plane shouldnt be a //valid muon plane if(currentPlane->GetPlaneNumber() >= fMuonStartPlaneNumber && (!currentPlane->IsValid() || currentPlane->GetPlaneType() == DmxPlaneTypes::kShower || setAll) ){ currentPlane->SetStrips(cog); MSG("AlgDeMuxBeam", Msg::kDebug) << "setting muon region plane " << currentPlane->GetPlaneNumber() << " to cog = " << cog << endl; }//end if this plane to be set planeItr.Next(); }//end loop over planes planeItr.ResetFirst(); return; }

void AlgDeMuxBeam::SetShowerRegionWindow (  DmxPlaneItr &  planeItr, Float_t  slope, Float_t  intercept )  [private]

Definition at line 1206 of file AlgDeMuxBeam.cxx. References fMuonStartPlaneNumber, DmxPlane::GetPlaneNumber(), DmxPlane::GetZPosition(), MSG, and DmxPlane::SetStrips(). Referenced by UseShowerSlidingWindow(). { MSG("AlgDeMuxBeam", Msg::kDebug) << "In SetShowerRegionWindow" << endl; Float_t cog = 0.; planeItr.ResetFirst(); DmxPlane *currentPlane = 0; while( (currentPlane = planeItr()) ){ //get the current plane MSG("AlgDeMuxBeam", Msg::kDebug) << "on plane " << currentPlane->GetPlaneNumber() << endl; cog = slope*currentPlane->GetZPosition() + intercept; if(fMuonStartPlaneNumber>0&&currentPlane->GetPlaneNumber()<fMuonStartPlaneNumber) currentPlane->SetStrips(cog); else if(fMuonStartPlaneNumber<0) currentPlane->SetStrips(cog); MSG("AlgDeMuxBeam", Msg::kDebug) << "setting shower region plane " << currentPlane->GetPlaneNumber() << " to " << cog << endl; }//end loop over planes planeItr.ResetFirst(); return; }

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

Reimplemented from AlgBase. Definition at line 637 of file AlgDeMuxBeam.cxx. { }

void AlgDeMuxBeam::UseMuonSlidingWindow (  DmxPlaneItr &  planeItr, DmxPlaneItr &  validPlaneItr, DmxStatus *  status )  [private]

Definition at line 1243 of file AlgDeMuxBeam.cxx. References fFinalShowerRegionSlope, FindFitOverNPlanes(), fInitialMuonRegionIntercept, fInitialMuonRegionSlope, fMuonSlopeLimit, DmxPlane::GetCoG(), DmxStatus::GetEndPlaneNumber(), DmxStatus::GetMultipleMuon(), DmxStatus::GetMuonStartPlaneNumber(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlaneTypes::GetPlaneTypeString(), DmxPlane::GetZPosition(), MSG, and SetMuonRegionWindow(). Referenced by RunAlg(). { MSG("AlgDeMuxBeam", Msg::kDebug) << "In UseMuonSlidingWindowTest" << endl; DmxPlane *plane = 0; validPlaneItr.ResetFirst(); Int_t loopCtr = 0; Int_t ctr = 0; Float_t slope = 0.; Float_t intercept = 0.; Float_t prevSlope = 0.; Float_t prevIntercept = 0.; //arrays to hold the information for fitting the straight lines //weight each plane by 1/charge Float_t tpos[10]; Float_t zpos[10]; Float_t weight[10]; Float_t residual[10]; Int_t planeNum[10]; Int_t muonPlanesInSet = fMuonPlanesInSet; //find the number of valid muon planes Int_t validPlaneCnt = 0; while( (plane = validPlaneItr()) ){ MSG("AlgDeMuxBeam", Msg::kDebug) << " on plane " << plane->GetPlaneNumber() << " type " << DmxPlaneTypes::GetPlaneTypeString(plane->GetPlaneType()) << endl; if(plane->GetPlaneType()==DmxPlaneTypes::kMuon) ++validPlaneCnt; } MSG("AlgDeMuxBeam", Msg::kDebug) << "num valid muon planes = " << validPlaneCnt << endl; //make sure you dont have fewer valid planes than the default expected if(validPlaneCnt<muonPlanesInSet) muonPlanesInSet = validPlaneCnt; validPlaneItr.ResetFirst(); //flag to reset valid muon planes if they are determined to be misreconstructions bool setAll = false; //loop over the valid planes and get the first n valid muon planes while( validPlaneItr.IsValid() && ctr<muonPlanesInSet){ plane = validPlaneItr.Ptr(); if(plane->GetPlaneType()==DmxPlaneTypes::kMuon){ planeNum[ctr] = plane->GetPlaneNumber(); zpos[ctr] = plane->GetZPosition(); tpos[ctr] = plane->GetCoG(); weight[ctr] = 1./plane->GetPlaneCharge(); ++ctr; } validPlaneItr.Next(); }//end loop to fill array for first n planes //find the fit for the first n planes and look for //planes with large residuals. repeat until planes with //large residuals are removed FindFitOverNPlanes(zpos,tpos,weight,residual,slope,intercept,ctr,planeNum); //loop over the weights and if any are 0, then set the setAll flag for(Int_t i = 0; i<ctr; ++i){ if(weight[i]==0.) setAll = true; } if(validPlaneCnt == muonPlanesInSet) planeNum[ctr-1] = status->GetEndPlaneNumber(); if(status->GetMuonStartPlaneNumber() < planeNum[0] && !status->GetMultipleMuon() ) planeNum[0] = status->GetMuonStartPlaneNumber(); //check that the slope seems reasonable - if not use the last good //slope and intercept if(TMath::Abs(slope) > fMuonSlopeLimit){ MSG("AlgDeMuxBeam", Msg::kDebug) << "slope = " << slope << "/" << fMuonSlopeLimit << "reset slope to " << fFinalShowerRegionSlope << endl; slope = fFinalShowerRegionSlope; intercept = fFinalShowerRegionIntercept; } //set the initial muon region slope and intercept variables fInitialMuonRegionSlope = slope; fInitialMuonRegionIntercept = intercept; prevSlope = slope; prevIntercept = intercept; //set the first planes to the slope and intercept found planeItr.GetSet()->Slice(planeNum[0], planeNum[ctr-1]+1); planeItr.GetSet()->Update(); SetMuonRegionWindow(planeItr, slope, intercept, setAll); planeItr.GetSet()->ClearSlice(false); //back the iterator up fMuonPlanesInSet-1 places, unless the are only //fMuonPlanesInSet. in that case, just return because the work here is //done if(validPlaneCnt > muonPlanesInSet){ validPlaneItr.Prev(); MSG("AlgDeMuxBeam", Msg::kDebug) << validPlaneItr.Ptr()->GetPlaneNumber() << "/" << planeNum[1] << endl; while(validPlaneItr.Ptr()->GetPlaneNumber() > planeNum[1]){ validPlaneItr.Prev(); MSG("AlgDeMuxBeam", Msg::kDebug) << validPlaneItr.Ptr()->GetPlaneNumber() << "/" << planeNum[1] << endl; } } else return; //now demux the event while(loopCtr < (validPlaneCnt - muonPlanesInSet) ){ setAll = false; ctr = 0; MSG("AlgDeMuxBeam", Msg::kDebug) << "find new window" << endl; while( validPlaneItr.IsValid() && ctr < muonPlanesInSet){ plane = validPlaneItr.Ptr(); MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber() << endl; if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){ planeNum[ctr] = plane->GetPlaneNumber(); tpos[ctr] = plane->GetCoG(); zpos[ctr] = plane->GetZPosition(); weight[ctr] = 1./plane->GetPlaneCharge(); ++ctr; MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber() << " z = " << zpos[ctr-1] << " cog = " << tpos[ctr-1] << endl; } validPlaneItr.Next(); }//end loop over valid planes //find the straight line to valid muon planes FindFitOverNPlanes(zpos, tpos, weight, residual, slope, intercept, ctr, planeNum); MSG("AlgDeMuxBeam", Msg::kDebug) << "slope = " << slope << " intercept = " << intercept << endl; //loop over the weights and if any are 0, then set the setAll flag for(Int_t i = 0; i<ctr; ++i){ if(weight[i]==0.) setAll = true; } //check that the slope seems reasonable - if not use the last good //slope and intercept if(TMath::Abs(slope) > fMuonSlopeLimit){ MSG("AlgDeMuxBeam", Msg::kDebug) << "slope = " << slope << "/" << fMuonSlopeLimit << "reset slope to " << prevSlope << endl; slope = prevSlope; intercept = prevIntercept; } //check to see if this is the last loop. if it is, make sure to set all unset planes //just use the fit found for the last window of valid planes if(loopCtr+1 == validPlaneCnt - muonPlanesInSet) planeNum[ctr-1] = status->GetEndPlaneNumber(); MSG("AlgDeMuxBeam", Msg::kDebug) << "loopCtr = " << loopCtr << " valid planes = " << validPlaneCnt << " last plane in set = " << planeNum[ctr-1] << " last plane in event " << status->GetEndPlaneNumber() << endl; //set the planes from the last set plane to the current one planeItr.GetSet()->Slice(planeNum[ctr-2]+1, planeNum[ctr-1]); planeItr.GetSet()->Update(); SetMuonRegionWindow(planeItr, slope, intercept, setAll); planeItr.GetSet()->ClearSlice(false); ++loopCtr; prevSlope = slope; prevIntercept = intercept; //back the iterator up validPlaneItr.Prev(); while( validPlaneItr.Ptr()->GetPlaneNumber() > planeNum[1]) validPlaneItr.Prev(); }//end loop over windows planeItr.GetSet()->ClearSlice(false); return; }

void AlgDeMuxBeam::UseReverseMuonSlidingWindow (  DmxPlaneItr &  planeItr, DmxPlaneItr &  validPlaneItr, DmxStatus *  status )  [private]

void AlgDeMuxBeam::UseShowerSlidingWindow (  DmxPlaneItr &  planeItr, DmxPlaneItr &  validPlaneItr, DmxStatus *  status )  [private]

Definition at line 975 of file AlgDeMuxBeam.cxx. References fFinalShowerRegionIntercept, fFinalShowerRegionSlope, FindFitOverNPlanes(), fMuonSlopeLimit, fMuonStartPlaneNumber, DmxPlane::GetCoG(), DmxStatus::GetEndPlaneNumber(), DmxStatus::GetMultipleMuon(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxStatus::GetVertexPlaneNumber(), DmxPlane::GetZPosition(), MSG, and SetShowerRegionWindow(). Referenced by RunAlg(). { MSG("AlgDeMuxBeam", Msg::kDebug) << "In UseShowerSlidingWindow" << endl; //figure out how many windows your are going to need for the shower region. //you need to find your lever arm //for the window too - it is either fPlanesInSet or the size of the validPlaneItr set. //do like in the cosmic algorithm, start from the vertex and go on. //only use the window until you get to the start of the muon track. validPlaneItr.ResetFirst(); DmxPlane *plane = 0; validPlaneItr.ResetFirst(); Int_t loopCtr = 0; Int_t ctr = 0; Float_t slope = 0.; Float_t intercept = 0.; //arrays to hold the information for fitting the straight lines //weight each plane by 1/charge Float_t tpos[10]; Float_t zpos[10]; Float_t weight[10]; Float_t residual[10]; Int_t planeNum[10]; Int_t planesInSet = fPlanesInSet; Float_t minResidual = 1.e20; Float_t minResidualTPos = 0.; //find the number of valid planes - in the shower region both muon and shower type planes //are ok to use Int_t validPlaneCnt = validPlaneItr.SizeSelect(); MSG("AlgDeMuxBeam", Msg::kDebug) << "num valid planes = " << validPlaneCnt << endl; if(validPlaneCnt==0){ MSG("AlgDeMuxBeam", Msg::kDebug) << "no valid planes in shower window" << endl; return; } //make sure you dont have fewer valid planes than the default expected if(validPlaneCnt<planesInSet) planesInSet = validPlaneCnt; validPlaneItr.ResetFirst(); //flag to reset valid muon planes if they are determined to be misreconstructions bool setAll = false; //loop over the valid planes and get the first n valid muon planes while( validPlaneItr.IsValid() && ctr<planesInSet){ plane = validPlaneItr.Ptr(); planeNum[ctr] = plane->GetPlaneNumber(); zpos[ctr] = plane->GetZPosition(); tpos[ctr] = plane->GetCoG(); weight[ctr] = 1./plane->GetPlaneCharge(); ++ctr; validPlaneItr.Next(); }//end loop to fill array for first n planes //find the fit for the first window of n planes. check to see if the //average residual is too big, if so then remove the last plane in the //window and redo the fit. FindFitOverNPlanes(zpos, tpos, weight, residual, slope, intercept, ctr, planeNum); //check that the slope seems reasonable - if not use the last good //slope and intercept if(TMath::Abs(slope) > fMuonSlopeLimit){ MSG("AlgDeMuxBeam", Msg::kDebug) << "slope = " << slope << "/" << fMuonSlopeLimit << "reset slope to " << fFinalShowerRegionSlope << endl; slope = fFinalShowerRegionSlope; intercept = fFinalShowerRegionIntercept; } //if there are no more planes to do the sliding window, check and see if there //looks to be a track. if no track, set the last plane to be set at the end of the //event if(validPlaneCnt == planesInSet && fMuonStartPlaneNumber<0) planeNum[ctr-1] = status->GetEndPlaneNumber(); else if(validPlaneCnt == planesInSet && fMuonStartPlaneNumber>0) planeNum[ctr-1] = fMuonStartPlaneNumber-1; if(status->GetVertexPlaneNumber() < planeNum[0] && !status->GetMultipleMuon() ) planeNum[0] = status->GetVertexPlaneNumber(); //set the first planes to the slope and intercept found planeItr.GetSet()->Slice(planeNum[ctr-2], planeNum[ctr-1]+1); planeItr.GetSet()->Update(); SetShowerRegionWindow(planeItr, slope, intercept); planeItr.GetSet()->ClearSlice(false); //back the iterator up fMuonPlanesInSet-1 places, unless the are only //fMuonPlanesInSet. in that case, just return because the work here is //done if(validPlaneCnt > planesInSet){ validPlaneItr.Prev(); MSG("AlgDeMuxBeam", Msg::kDebug) << validPlaneItr.Ptr()->GetPlaneNumber() << "/" << planeNum[1] << endl; while(validPlaneItr.Ptr()->GetPlaneNumber() > planeNum[1]){ validPlaneItr.Prev(); MSG("AlgDeMuxBeam", Msg::kDebug) << validPlaneItr.Ptr()->GetPlaneNumber() << "/" << planeNum[1] << endl; } } else return; //now demux the event while(loopCtr < (validPlaneCnt - planesInSet) ){ setAll = false; ctr = 0; MSG("AlgDeMuxBeam", Msg::kDebug) << "find new window" << endl; //get the tpos for each of the previously set planes while( validPlaneItr.IsValid() && ctr < planesInSet-1){ plane = validPlaneItr.Ptr(); MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber() << endl; tpos[ctr] = plane->GetCoG(); zpos[ctr] = plane->GetZPosition(); weight[ctr] = 1./plane->GetPlaneCharge(); planeNum[ctr] = plane->GetPlaneNumber(); ++ctr; MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber() << " z = " << zpos[ctr-1] << " cog = " << tpos[ctr-1] << endl; validPlaneItr.Next(); }//end loop over valid planes //fit the previously set planes FindFitOverNPlanes(zpos, tpos, weight, residual, slope, intercept, ctr, planeNum); //now check out the last plane - if it is a shower plane, then loop over the 3 best hypotheses //to see which is closest to the line fit to the previously set planes plane = validPlaneItr.Ptr(); MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber() << endl; //get the plane info that is independant of plane type zpos[ctr] = plane->GetZPosition(); weight[ctr] = 1./plane->GetPlaneCharge(); planeNum[ctr] = plane->GetPlaneNumber(); minResidual = 1.e20; //loop over the best 3 hyps to see if any of them fit well with the track for(Int_t i = 0; i<3; ++i){ tpos[ctr] = plane->GetCoG(); minResidualTPos = tpos[ctr]; if(plane->GetPlaneType() == DmxPlaneTypes::kShower){ if(i<1) tpos[ctr] = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG(); else if(i<2) tpos[ctr] = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG(); else if(i<3) tpos[ctr] = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG(); } //muon planes only have one possible cog else break; residual[ctr] = TMath::Abs(tpos[ctr] - slope*zpos[ctr]+intercept); //check the residual for the last plane if(residual[ctr]<minResidual){ minResidual = residual[ctr]; minResidualTPos = tpos[ctr]; } }//end loop over best 3 hypotheses for this plane //use the transverse position with the minimum residual tpos[ctr] = minResidualTPos; //increment the ctr for later ++ctr; //found the best fit of the hypotheses, so now fit the entire set FindFitOverNPlanes(zpos, tpos, weight, residual, slope, intercept, ctr, planeNum); MSG("AlgDeMuxBeam", Msg::kDebug) << "fit over all planes in window slope = " << slope << " intercept = " << intercept << endl; //check that the slope seems reasonable - if not use the last good //slope and intercept if(TMath::Abs(slope) > fMuonSlopeLimit){ MSG("AlgDeMuxBeam", Msg::kDebug) << "slope = " << slope << "/" << fMuonSlopeLimit << "reset slope to " << fFinalShowerRegionSlope << endl; slope = fFinalShowerRegionSlope; intercept = fFinalShowerRegionIntercept; } //check to see if this is the last loop. if it is, make sure to set all unset planes //just use the fit found for the last window of valid planes - if no track region //is present, just set all the remaining planes in the detector - what can it hurt? if(loopCtr+1 == validPlaneCnt - planesInSet && fMuonStartPlaneNumber<0.) planeNum[ctr-1] = 486; else if(loopCtr+1 == validPlaneCnt - planesInSet && fMuonStartPlaneNumber>0.) planeNum[ctr-1] = fMuonStartPlaneNumber-1; MSG("AlgDeMuxBeam", Msg::kDebug) << "loopCtr = " << loopCtr << " valid planes = " << validPlaneCnt << " last plane in set = " << planeNum[ctr-1] << " last plane in event " << status->GetEndPlaneNumber() << endl; planeItr.GetSet()->Slice(planeNum[ctr-2]+1, planeNum[ctr-1]+1); planeItr.GetSet()->Update(); SetShowerRegionWindow(planeItr, slope, intercept); planeItr.GetSet()->ClearSlice(false); ++loopCtr; //set the final shower region slope and intercept to the current values fFinalShowerRegionSlope = slope; fFinalShowerRegionIntercept = intercept; //back the iterator up validPlaneItr.Prev(); while( validPlaneItr.Ptr()->GetPlaneNumber() > planeNum[1]) validPlaneItr.Prev(); }//end loop over windows planeItr.GetSet()->ClearSlice(false); return; }

---

Member Data Documentation

Float_t AlgDeMuxBeam::fDigitResetCut [private]

Definition at line 44 of file AlgDeMuxBeam.h. Referenced by ReconcileShowerAndMuonRegions(), and RunAlg().

Float_t AlgDeMuxBeam::fFinalShowerRegionIntercept [private]

Definition at line 41 of file AlgDeMuxBeam.h. Referenced by DeMuxFirstNPlanesTest(), and UseShowerSlidingWindow().

Float_t AlgDeMuxBeam::fFinalShowerRegionSlope [private]

Definition at line 40 of file AlgDeMuxBeam.h. Referenced by DeMuxFirstNPlanesTest(), ReconcileShowerAndMuonRegions(), UseMuonSlidingWindow(), and UseShowerSlidingWindow().

Float_t AlgDeMuxBeam::fInitialMuonRegionIntercept [private]

Definition at line 43 of file AlgDeMuxBeam.h. Referenced by UseMuonSlidingWindow().

Float_t AlgDeMuxBeam::fInitialMuonRegionSlope [private]

Definition at line 42 of file AlgDeMuxBeam.h. Referenced by ReconcileShowerAndMuonRegions(), and UseMuonSlidingWindow().

Int_t AlgDeMuxBeam::fMaxBadResidual [private]

Definition at line 45 of file AlgDeMuxBeam.h. Referenced by RunAlg().

Float_t AlgDeMuxBeam::fMaxMuonRemovalFraction [private]

Definition at line 47 of file AlgDeMuxBeam.h. Referenced by FindFitOverNPlanes(), and RunAlg().

Float_t AlgDeMuxBeam::fMaxResidual [private]

Definition at line 46 of file AlgDeMuxBeam.h. Referenced by FindFitOverNPlanes(), and RunAlg().

Int_t AlgDeMuxBeam::fMaxStripAdjustment [private]

Definition at line 48 of file AlgDeMuxBeam.h. Referenced by RunAlg().

Float_t AlgDeMuxBeam::fMinMatedChargeFrac [private]

Definition at line 50 of file AlgDeMuxBeam.h. Referenced by RunAlg().

Int_t AlgDeMuxBeam::fMinMuonPlanesRequired [private]

Definition at line 49 of file AlgDeMuxBeam.h. Referenced by RunAlg().

Int_t AlgDeMuxBeam::fMuonPlanesInSet [private]

Definition at line 51 of file AlgDeMuxBeam.h. Referenced by RunAlg().

Float_t AlgDeMuxBeam::fMuonSlopeLimit [private]

Definition at line 52 of file AlgDeMuxBeam.h. Referenced by RunAlg(), UseMuonSlidingWindow(), and UseShowerSlidingWindow().

Int_t AlgDeMuxBeam::fMuonStartPlaneNumber [private]

Definition at line 53 of file AlgDeMuxBeam.h. Referenced by ReconcileShowerAndMuonRegions(), RunAlg(), SetMuonRegionWindow(), SetShowerRegionWindow(), and UseShowerSlidingWindow().

Float_t AlgDeMuxBeam::fMuonStartPlaneZPos [private]

Definition at line 54 of file AlgDeMuxBeam.h. Referenced by RunAlg().

Int_t AlgDeMuxBeam::fNumberOfHypotheses [private]

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

Int_t AlgDeMuxBeam::fPlanesInSet [private]

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

Float_t AlgDeMuxBeam::fShowerMuonTransverseDifference [private]

Definition at line 59 of file AlgDeMuxBeam.h. Referenced by ReconcileShowerAndMuonRegions(), and RunAlg().

Int_t AlgDeMuxBeam::fStrayCut [private]

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

Int_t AlgDeMuxBeam::fStrayPlanes [private]

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

VldContext AlgDeMuxBeam::fVldContext [private]

Definition at line 60 of file AlgDeMuxBeam.h. Referenced by DeMuxFirstNPlanesTest(), FindDigitCoG(), ReconcileShowerAndMuonRegions(), and RunAlg().

---

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

• AlgDeMuxBeam.h
• AlgDeMuxBeam.cxx

---