AlgDeMuxCosmics Class Reference

#include <AlgDeMuxCosmics.h>

Inheritance diagram for AlgDeMuxCosmics:

List of all members.

Public Member Functions
	AlgDeMuxCosmics ()
virtual	~AlgDeMuxCosmics ()
virtual void	Trace (const char *c) const
virtual void	RunAlg (AlgConfig &acd, CandHandle &ch, CandContext &cx)
Private Member Functions
void	FindFit (DmxPlaneItr &planeItr, Double_t &prevChiSq, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &a4, Int_t *hypotheses, Int_t direction, Int_t leverArm, Float_t offset)
Bool_t	FindPlanesToDropFromFit (DmxPlaneItr &planeItr, Double_t a1, Double_t a2, Double_t a3, Double_t a4, Int_t direction, Int_t leverArm, Float_t offset)
void	SetPlanesToDeterminedFit (DmxPlaneItr &planeItr, Double_t a1, Double_t a2, Double_t a3, Double_t a4, Float_t offset)
void	SetPlanesToDeterminedFit (DmxPlaneItr &planeItr, Double_t a, Double_t b)
void	FindCosmicSolution (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, DmxStatus *status, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &a4, Double_t &chiSq, Float_t offset)
void	FindWindowCosmicSolution (DmxPlaneItr &validPlaneItr, DmxStatus *status, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &a4, Double_t &chiSq, Float_t offset)
void	UseSingleFit (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, DmxStatus *status)
void	UseSlidingWindow (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, DmxPlaneItr &windowPlaneItr, DmxStatus *status, Int_t vertex, Int_t endPlane)
Private Attributes
Double_t	fRequiredMatedSignalRatio
Double_t	fSlopeRMS
Int_t	fStrayPlanes
UInt_t	fPlanesInSet
Int_t	fStrayCut

---

Constructor & Destructor Documentation

AlgDeMuxCosmics::AlgDeMuxCosmics (   )

Definition at line 73 of file AlgDeMuxCosmics.cxx. : fRequiredMatedSignalRatio(0.33), fSlopeRMS(0.), fStrayPlanes(0), fPlanesInSet(6), fStrayCut(6) { //default constructor }

AlgDeMuxCosmics::~AlgDeMuxCosmics (   )  [virtual]

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

---

Member Function Documentation

void AlgDeMuxCosmics::FindCosmicSolution (  DmxPlaneItr &  planeItr, DmxPlaneItr &  validPlaneItr, DmxStatus *  status, Double_t &  a1, Double_t &  a2, Double_t &  a3, Double_t &  a4, Double_t &  chiSq, Float_t  offset )  [private]

Definition at line 558 of file AlgDeMuxCosmics.cxx. References FindFit(), FindPlanesToDropFromFit(), DmxStatus::GetEventDirection(), and DmxUtilities::IsValidFit(). Referenced by UseSingleFit(), and UseSlidingWindow(). { //variables to keep track of the best reconstruction Double_t fitA1 = 0.; Double_t fitA2 = 0.; Double_t fitA3 = 0.; Double_t fitA4 = 0.; Double_t useA1 = -10000.; Double_t useA2 = -10000.; Double_t useA3 = -10000.; Double_t useA4 = -10000.; Double_t bestChi2 = 1.e20; Double_t chi2 = 1.000001e20; Int_t direction = status->GetEventDirection(); DmxUtilities util; //make the lever arm for the demuxer the # of valid planes or 6, whichever //is smaller UInt_t leverArm = validPlaneItr.SizeSelect(); if(leverArm > fPlanesInSet){ leverArm = fPlanesInSet;} //MSG("DMXX", Msg::kDebug) << "in FindCosmicSolution" << endl; //MSG("DMX", Msg::kDebug) << "\tlever arm = " << leverArm << endl; //declare an array to tell the fitter which hypotheses to use Int_t hypsToUse[] = {0,0,0,0,0,0}; for(Int_t plane0Hyp = 1; plane0Hyp < 4; plane0Hyp++){ hypsToUse[0] = plane0Hyp; for(Int_t plane1Hyp = 1; plane1Hyp < 4; plane1Hyp++){ hypsToUse[1] = plane1Hyp; if(leverArm == 2){ FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset); if(chi2 < bestChi2){ bestChi2 = chi2; useA1 = fitA1; useA2 = fitA2; useA3 = fitA3; useA4 = fitA4; }//end if its a better straight line fit } else if(leverArm >= 3){ for(Int_t plane2Hyp = 1; plane2Hyp < 4; plane2Hyp++){ hypsToUse[2] = plane2Hyp; if(leverArm == 3){ FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset); if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){ bestChi2 = chi2; useA1 = fitA1; useA2 = fitA2; useA3 = fitA3; useA4 = fitA4; //see if you can drop some planes and improve the fit if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){ FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset); if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){ bestChi2 = chi2; useA1 = fitA1; useA2 = fitA2; useA3 = fitA3; useA4 = 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 } else if(leverArm >= 4){ for(Int_t plane3Hyp = 1; plane3Hyp < 4; plane3Hyp++){ hypsToUse[3] = plane3Hyp; if(leverArm == 4){ FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset); if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){ bestChi2 = chi2; useA1 = fitA1; useA2 = fitA2; useA3 = fitA3; useA4 = fitA4; //see if you can drop some planes and improve the fit if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){ FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset); if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){ bestChi2 = chi2; useA1 = fitA1; useA2 = fitA2; useA3 = fitA3; useA4 = 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 } else if(leverArm >= 5){ for(Int_t plane4Hyp = 1; plane4Hyp < 4; plane4Hyp++){ hypsToUse[4] = plane4Hyp; if(leverArm == 5){ FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset); if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){ bestChi2 = chi2; useA1 = fitA1; useA2 = fitA2; useA3 = fitA3; useA4 = fitA4; //see if you can drop some planes and improve the fit if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){ FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset); if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){ bestChi2 = chi2; useA1 = fitA1; useA2 = fitA2; useA3 = fitA3; useA4 = 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 } else if(leverArm == 6){ for(Int_t plane5Hyp = 1; plane5Hyp < 4; plane5Hyp++){ hypsToUse[5] = plane5Hyp; FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset); if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){ bestChi2 = chi2; useA1 = fitA1; useA2 = fitA2; useA3 = fitA3; useA4 = fitA4; //see if you can drop some planes and improve the fit if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){ FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset); if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){ bestChi2 = chi2; useA1 = fitA1; useA2 = fitA2; useA3 = fitA3; useA4 = 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 }//end for loop over 6th plane's hyps }//end if using 6 planes }//end for loop over 5th plane's hyps }//end if at least 5 planes } //loop over 4th plane } //end if at least 4 planes } //loop over 3rd plane } //end if at least 3 planes } //loop over 2nd plane } //loop over 1st plane //force the fit to the one found if(useA1 != -10000. && useA2 != -10000. && useA3 != 10000. && useA4 != 10000.){ a1 = useA1; a2 = useA2; a3 = useA3; a4 = useA4; chiSq = bestChi2; } //MSG("DMX", Msg::kDebug) << "finished in FindCosmicsSolution" << endl; return; }

void AlgDeMuxCosmics::FindFit (  DmxPlaneItr &  planeItr, Double_t &  prevChiSq, Double_t &  a1, Double_t &  a2, Double_t &  a3, Double_t &  a4, Int_t *  hypotheses, Int_t  direction, Int_t  leverArm, Float_t  offset )  [private]

Definition at line 1476 of file AlgDeMuxCosmics.cxx. References DmxUtilities::FindLinearFit(), DmxPlane::GetCoG(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneType(), DmxPlane::GetZPosition(), and DmxPlane::IsGolden(). Referenced by FindCosmicSolution(). { //MSG("DMXX", Msg::kDebug) << "in FindFit" << endl; //cout << "in find fit for first n planes" << endl; Double_t a1L = 0.; Double_t a2L = 0.; Double_t a1Q = 0.; Double_t a2Q = 0.; Double_t a3Q = 0.; Double_t chiSqL = 1.e20; Double_t chiSqQ = 1.e20; //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.}; DmxUtilities util; //cout << "got util object" << endl; if(direction == -1){ planeItr.ResetLast(); } else if( direction == 1){ planeItr.ResetFirst(); } //planeItr.ResetFirst(); Int_t planeCtr = 0; Int_t arrayCtr = 0; DmxPlane *plane = 0; //MSG("DMX", Msg::kDebug) << "\t" << hypotheses[0] << "\t" << hypotheses[1] // << "\t" << hypotheses[2] << "\t" << hypotheses[3] // << "\t" << hypotheses[4] << "\t" << hypotheses[5] << endl; //loop over the plane iterator to fill the variables while( (plane = planeItr()) && planeCtr < leverArm){ if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){ x[arrayCtr] = plane->GetZPosition() - offset; //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){ //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetCoG(); 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){ //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetCoG(); y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG(); weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio()); } else if(hypotheses[planeCtr] == 3){ //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetCoG(); y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG(); weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio()); } } ++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; Double_t chiSq = 0.; //use the DmxUtilities function to find a linear fit util.FindLinearFit(x, y, weight, arrayCtr, a1L, a2L, chiSqL); //util.FindQuadraticFit(x, y, weight, arrayCtr, a1Q, a2Q, a3Q, chiSqQ); //util.FindCubicFit(x, y, weight, arrayCtr, a1, a2, a3, a4, chiSq); //MSG("DMXX", Msg::kDebug) << "\tlinear fit " << chiSqL << "\t" << a1L << "\t" << a2L << endl; //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q // << "\t" << a2Q << "\t" << a3Q << endl; //pick the better fit if(chiSqL<=chiSqQ){ prevChiSq = chiSqL; a1 = a1L; a2 = a2L; a3 = 0.; a4 = 0.; } else{ prevChiSq = chiSqQ; a1 = a1Q; a2 = a2Q; a3 = a3Q; a4 = 0.; } //MSG("DMX", Msg::kDebug) << "\tinitial fit " << a1 << "\t" << a2 << "\t" << a3 << "\t" << a4 << 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){ a1L = 0.; a2L = 0.; a1Q = 0.; a2Q = 0.; a3Q = 0.; 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){ x[arrayCtr] = plane->GetZPosition() - offset; y[arrayCtr] = plane->GetCoG(); //cout << "filling arrays " << arrayCtr << " " << x[arrayCtr] << " " << y[arrayCtr] // << " " << weightSq[arrayCtr] << endl; 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("DMX", 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, a1L, a2L, chiSqL); //util.FindQuadraticFit(x, y, weight, arrayCtr, a1Q, a2Q, a3Q, chiSqQ); //util.FindCubicFit(x, y, weight, arrayCtr, fitA1, fitA2, fitA3, fitA4, chiSq); //MSG("DMXX", Msg::kDebug) << "\tlinear fit " << chiSqL << "\t" << a1L << "\t" << a2L << endl; //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q // << "\t" << a2Q << "\t" << a3Q << endl; //pick the better fit if(chiSqL<=chiSqQ){ chiSq = chiSqL; fitA1 = a1L; fitA2 = a2L; } else{ chiSq = chiSqQ; fitA1 = a1Q; fitA2 = a2Q; fitA3 = a3Q; } planeItr.Reset(); planeCtr = 0; if(chiSq < prevChiSq){ prevChiSq = chiSq; a1 = fitA1; a2 = fitA2; a3 = fitA3; a4 = fitA4; //MSG("DMX", Msg::kDebug) << "\trefined fit, drop plane " << i // << "\t" << prevChiSq << "\t" << a << "\t" << b << endl; } }//end for loop to improve fit by dropping a plane }//end if enough planes to improve fit planeItr.ResetFirst(); //MSG("DMX", Msg::kDebug) << "\tfinal fit" << "\t" << a << "\t" << b << "\t" << prevChiSq << endl; //cout << "finished find fit for first n planes" << endl; return; }

Bool_t AlgDeMuxCosmics::FindPlanesToDropFromFit (  DmxPlaneItr &  planeItr, Double_t  a1, Double_t  a2, Double_t  a3, Double_t  a4, Int_t  direction, Int_t  leverArm, Float_t  offset )  [private]

Definition at line 1789 of file AlgDeMuxCosmics.cxx. References DmxPlane::GetCoG(), DmxPlane::GetPlaneType(), DmxPlane::GetZPosition(), DmxPlane::IsGolden(), and DmxPlane::SetGolden(). Referenced by FindCosmicSolution(). { //MSG("DMXX", Msg::kDebug) << "in FindPlanesToDropFromFit" << endl; Double_t diff1 = 0.; Double_t diff2 = 0.; Double_t diff3 = 0.; bool planesDropped = false; Int_t dropPlanes = 0; Int_t planeCtr = 0; if( direction == 1 ){ planeItr.ResetFirst(); } else if( direction == -1 ){ planeItr.ResetLast(); } //MSG("DMX", 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){ Double_t fitCoG = (a1 + (plane->GetZPosition()-offset)*a2 + TMath::Power(plane->GetZPosition()-offset,2)*a3 + TMath::Power(plane->GetZPosition()-offset,3)*a4); //MSG("DMX", 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("DMX", 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("DMX", 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; }

void AlgDeMuxCosmics::FindWindowCosmicSolution (  DmxPlaneItr &  validPlaneItr, DmxStatus *  status, Double_t &  a1, Double_t &  a2, Double_t &  a3, Double_t &  a4, Double_t &  chiSq, Float_t  offset )  [private]

Definition at line 746 of file AlgDeMuxCosmics.cxx. References DmxUtilities::FindLinearFit(), DmxPlane::GetCoG(), DmxStatus::GetEventDirection(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlane::GetZPosition(), DmxPlane::IsGolden(), MSG, and DmxPlane::SetGolden(). Referenced by UseSlidingWindow(). { //MSG("DMXX", Msg::kDebug) << "in FindWindowCosmicSolution" << endl; //variables to keep track of the best reconstruction Double_t fitA1L = 0.; Double_t fitA2L = 0.; Double_t chi2L = 1.0e20; // Double_t fitA1Q = 0.; // Double_t fitA2Q = 0.; // Double_t fitA3Q = 0.; // Double_t chi2Q = 1.0e20; Double_t fitA1 = 0.; Double_t fitA2 = 0.; Double_t fitA3 = 0.; Double_t fitA4 = 0.; Double_t chi2 = 1.0e20; Double_t chi2Fit = 1.0e20; Int_t direction = status->GetEventDirection(); //make the lever arm for the demuxer the # of valid planes or 6, whichever //is smaller Int_t leverArm = fPlanesInSet; //offset is the lowest z pos in the window. subtract it from each x array value //- this has the effect of changing your y axis location //relative to the z = 0 location so that when doing higher order fits, you //dont get really stupid numbers //declare the variables to do a least squares fit to a straight line. //use three chiSq values to keep track of the 3 different hyps in the //unset plane - also have 3 different y values for those hyps and 3 //different weights 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 weight2 = 1.; Double_t weight3 = 1.; Double_t y2 = 0.; Double_t y3 = 0.; //variable to let you know if the unset plane is a golden plane or not bool lastGolden = false; Double_t a1L = 0.; Double_t a2L = 0.; Double_t a1LB = 0.; Double_t a2LB = 0.; Double_t a1LS = 0.; Double_t a2LS = 0.; Double_t a1LT = 0.; Double_t a2LT = 0.; // Double_t a1QB = 0.; // Double_t a2QB = 0.; // Double_t a3QB = 0.; // Double_t a1QS = 0.; // Double_t a2QS = 0.; // Double_t a3QS = 0.; // Double_t a1QT = 0.; // Double_t a2QT = 0.; // Double_t a3QT = 0.; // Double_t a1Q = 0.; // Double_t a2Q = 0.; // Double_t a3Q = 0.; Double_t chiSqL = 1.e20; Double_t chiSqLB = 1.e20; Double_t chiSqLS = 1.e20; Double_t chiSqLT = 1.e20; // Double_t chiSqQB = 1.e20; // Double_t chiSqQS = 1.e20; // Double_t chiSqQT = 1.e20; // Double_t chiSqQ = 1.e20; Double_t a1Fit = 0.; Double_t a2Fit = 0.; Double_t a3Fit = 0.; Double_t a4Fit = 0.; DmxUtilities util; if(direction == -1) validPlaneItr.ResetLast(); else if( direction == 1) validPlaneItr.ResetFirst(); Int_t planeCtr = 0; DmxPlane *plane = 0; while( (plane = validPlaneItr()) && planeCtr < leverArm){ //cout << "plane = " << plane->GetPlaneNumber() << endl; x[planeCtr] = plane->GetZPosition() - offset; y[planeCtr] = plane->GetCoG(); weight[planeCtr] = 1. / TMath::Sqrt(plane->GetPlaneCharge()); if(plane->IsGolden()){ //golden plane, so make the weight smaller to give it more //clout in the fit weight[planeCtr] /= TMath::Sqrt(10.); if(planeCtr == leverArm-1){ lastGolden = true; y2 = y[planeCtr]; y3 = y[planeCtr]; weight2 = weight[planeCtr]; weight3 = weight[planeCtr]; } } //only make changes if not a golden plane if(plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden() ){ //make sure you dont divide by zero when getting the mated signal ratio of //the set hypothesis - could be that you set the plane to a hypothesis //with no mated signal, even though it is a valid plane if(planeCtr < leverArm-1 && dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio() > 0.){ weight[planeCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio()); // MSG("Dmx", Msg::kDebug) << dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio() << endl; } else if( planeCtr == leverArm-1 ){ y2 = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG(); y3 = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG(); //dont have to worry about dividing by zero here because you havent //set the plane yet and you know that the best hypotheses in the planes //must have at least 0.3 mated signal weight2 = weight[planeCtr]/TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio()); weight3 = weight[planeCtr]/TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio()); weight[planeCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()); } } ++planeCtr; } //end loop over planes to find variables for fit //do fit using function in DmxUtilities //the best hyp fit first util.FindLinearFit(x, y, weight, leverArm, a1LB, a2LB, chiSqLB); //util.FindQuadraticFit(x, y, weight, leverArm, a1QB, a2QB, a3QB, chiSqQB); //find fit for using the last plane's second best hyp y[leverArm-1] = y2; weight[leverArm-1] = weight2; util.FindLinearFit(x, y, weight, leverArm, a1LS, a2LS, chiSqLS); //util.FindQuadraticFit(x, y, weight, leverArm, a1QS, a2QS, a3QS, chiSqQS); //find fit using the last plane's third best hyp y[leverArm-1] = y3; weight[leverArm-1] = weight3; util.FindLinearFit(x, y, weight, leverArm, a1LT, a2LT, chiSqLT); //util.FindQuadraticFit(x, y, weight, leverArm, a1QT, a2QT, a3QT, chiSqQT); //util.FindCubicFit(x, y, weight, leverArm, a1, a2, a3, a4, chiSq); // MSG("DMXX", Msg::kDebug) << "\tlinear fit1 " << chiSqLB << "\t" << a1LB << "\t" << a2LB << endl; // MSG("DMXX", Msg::kDebug) << "\tlinear fit2 " << chiSqLS << "\t" << a1LS << "\t" << a2LS << endl; // MSG("DMXX", Msg::kDebug) << "\tlinear fit3 " << chiSqLT << "\t" << a1LT << "\t" << a2LT << endl; //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q // << "\t" << a2Q << "\t" << a3Q << endl; //pick the better linear fit if(chiSqLB<=chiSqLS && chiSqLB<=chiSqLT){ chiSqL = chiSqLB; a1L = a1LB; a2L = a2LB; } else if(chiSqLS<chiSqLB && chiSqLS<=chiSqLT){ chiSqL = chiSqLS; a1L = a1LS; a2L = a2LS; } else if(chiSqLT<chiSqLB && chiSqLT<chiSqLS){ chiSqL = chiSqLT; a1L = a1LT; a2L = a2LT; } //pick the better quadratic fit // if(chiSqQB<=chiSqQS && chiSqQB<=chiSqQT){ // chiSqQ = chiSqQB; // a1Q = a1QB; // a2Q = a2QB; // a3Q = a3QB; // } // else if(chiSqQS<chiSqQB && chiSqQS<=chiSqQT){ // chiSqQ = chiSqQS; // a1Q = a1QS; // a2Q = a2QS; // a3Q = a3QS; // } // else if(chiSqQT<chiSqQB && chiSqQT<chiSqQS){ // chiSqQ = chiSqQT; // a1Q = a1QT; // a2Q = a2QT; // a3Q = a3QT; // } //pick the better fit - linear or quadratic //if(chiSqL<=chiSqQ){ chi2Fit = chiSqL; a1Fit = a1L; a2Fit = a2L; a3Fit = 0.; a4Fit = 0.; //} // else{ // chi2Fit = chiSqQ; // a1Fit = a1Q; // a2Fit = a2Q; // a3Fit = a3Q; // a4Fit = 0.; // } //MSG("DMXX", Msg::kDebug) << "\tinitial fit " << chi2Fit << "\t" << a1Fit << "\t" << a2Fit // << "\t" << a3Fit << "\t" << a4Fit << endl; //redo the fit using only the previously set planes and only if you have more than 3 planes in the set if(leverArm > 3){ //do fit using function in DmxUtilities util.FindLinearFit(x, y, weight, leverArm-1, fitA1L, fitA2L, chi2L); //util.FindQuadraticFit(x, y, weight, planeCtr, fitA1Q, fitA2Q, fitA3Q, chi2Q); //util.FindCubicFit(x, y, weight, planeCtr, fitA1, fitA2, fitA3, fitA4, chiSq); //MSG("DMXX", Msg::kDebug) << "\trefined linear fit " << chi2L << "\t" << fitA1L << "\t" << fitA2L << endl; //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q // << "\t" << a2Q << "\t" << a3Q << endl; //pick the better fit - linear or quadratic //if(chi2L<=chi2Q){ chi2 = chi2L; fitA1 = fitA1L; fitA2 = fitA2L; fitA3 = 0.; fitA4 = 0.; //} // else{ // chi2 = chi2Q; // fitA1 = fitA1Q; // fitA2 = fitA2Q; // fitA3 = fitA3Q; // fitA4 = 0.; // } if( lastGolden ){ if(direction == 1) validPlaneItr.ResetLast(); else if( direction == -1) validPlaneItr.ResetFirst(); plane = validPlaneItr(); MSG("AlgDmx", Msg::kDebug) << "last plane number = " << plane->GetPlaneNumber() << " this is a golden plane " << plane->IsGolden() << endl; //unset the IsGolden flag for the last plane if the difference between the //fit center of gravity and the plane's center of gravity is > 0.5m and //any of the following //1. the golden CoG is different from the previous set valid plane's // CoG by more than 0.25m/plane (plane to plane spacing is 0.1188m) //2. chi^2 of fit without the golden plane is an order of magnitude // smaller than the chi^2 of the fit with it //3. the average difference between the set CoG's and the golden CoG and // the fit CoG's for all planes in the window is greater than 0.5m //cout << x[leverArm-1] << " " << x[leverArm-2] << " " << a1Fit << " " << a2Fit << endl; Double_t fitCoG = (a1Fit + a2Fit*(x[leverArm-1]) + a3Fit*TMath::Power(x[leverArm-1],2) + a4Fit*TMath::Power(x[leverArm-1],3)); Double_t fitDiff = 0.; for(Int_t ii = 0; ii<leverArm; ii++){ fitDiff += TMath::Abs((a1Fit + a2Fit*(x[ii]) + a3Fit*TMath::Power(x[ii],2) + a4Fit*TMath::Power(x[ii],3) - y[ii])); } Double_t planeDiff = TMath::Abs(x[leverArm-1]-x[leverArm-2])/0.1188; MSG("DmxAlg", Msg::kDebug) << "fit is " << fitCoG << " golden CoG = " << plane->GetCoG() << endl <<" y diff = " << TMath::Abs(y[leverArm-2]-plane->GetCoG()) << "/" << planeDiff*.25 << endl <<" fitDiff = " << fitDiff <<"/" << 0.5*leverArm << endl << "chi2Fit = " << chi2Fit << "/" << 10.*chi2 << endl; if(TMath::Abs(plane->GetCoG()-fitCoG)>0.5 && (fitDiff>0.5*leverArm || TMath::Abs(y[leverArm-2]-plane->GetCoG())>(0.25*planeDiff) || chi2Fit>10.*chi2)){ plane->SetGolden(false); MSG("DmxAlg", Msg::kDebug) << plane->GetPlaneNumber() << " was golden " << plane->IsGolden() << endl; lastGolden = false; } }//end check of last plane is golden //if the fit found by dropping the last plane has a better chi^2 than the //fit using all planes, and that last plane was not golden, then pass the //fit from the first planes in the window out. if the last plane was initially //golden, this only happens for the last time through. // if(chi2 < chi2Fit && !lastGolden){ chi2Fit = chi2; a1Fit = fitA1; a2Fit = fitA2; a3Fit = fitA3; a4Fit = fitA4; } }//end for loop to improve fit by dropping a plane chiSq = chi2Fit; a1 = a1Fit; a2 = a2Fit; a3 = a3Fit; a4 = a4Fit; //MSG("DMXX", Msg::kDebug) << "\tchosen linear fit " << chiSq << "\t" << a1 << "\t" << a2 << endl; validPlaneItr.Reset(); return; }

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

Implements AlgBase. Definition at line 93 of file AlgDeMuxCosmics.cxx. References RawTriggerCodes::AsString(), DmxUtilities::CheckFitWithTiming(), DmxUtilities::FillPlaneArray(), DmxUtilities::FindEndPlane(), DmxUtilities::FindPlanesOffFit(), DmxUtilities::FindVertexPlane(), fPlanesInSet, fRequiredMatedSignalRatio, fSlopeRMS, fStrayCut, fStrayPlanes, DmxStatus::GetAverageTimingOffset(), CandRecord::GetCandHeader(), CandHandle::GetCandRecord(), VldContext::GetDetector(), Registry::GetDouble(), DmxStatus::GetEndPlaneNumber(), DmxStatus::GetEventDeMuxed(), DmxStatus::GetEventNumber(), DmxStatus::GetFigureOfMeritFailure(), Registry::GetInt(), CandHandle::GetNDaughters(), DmxStatus::GetNonPhysicalFailure(), DmxStatus::GetPlaneArray(), CandHeader::GetSnarl(), DmxStatus::GetVertexPlaneNumber(), RecMinosHdr::GetVldContext(), DmxUtilities::IsOverlappingMultiple(), KeyOnPlane(), KeyOnUView(), KeyOnValidU(), KeyOnValidV(), KeyOnVView(), MSG, DmxStatus::ResetStatus(), DmxStatus::SetAverageTimingOffset(), CandDeMuxDigitListHandle::SetAvgTimeOffset(), CandDeMuxDigitListHandle::SetDeMuxDigitListFlagBit(), DmxStatus::SetEndPlaneNumber(), DmxStatus::SetEventNumber(), DmxStatus::SetFigureOfMerit(), DmxStatus::SetMultipleMuon(), DmxStatus::SetNoVertexFailure(), DmxStatus::SetNumberOfPlanes(), CandDeMuxDigitListHandle::SetNumStrayPlanesU(), CandDeMuxDigitListHandle::SetNumStrayPlanesV(), CandDeMuxDigitListHandle::SetNumValidPlanesU(), CandDeMuxDigitListHandle::SetNumValidPlanesV(), SetPlanesToDeterminedFit(), DmxStatus::SetUOverlappingMultiple(), DmxStatus::SetUStrayPlanes(), DmxStatus::SetUValidPlanes(), DmxStatus::SetValidPlanesFailure(), DmxStatus::SetVertexPlaneNumber(), DmxStatus::SetVertexPlaneZPosition(), DmxStatus::SetVOverlappingMultiple(), DmxStatus::SetVStrayPlanes(), DmxStatus::SetVValidPlanes(), UseSingleFit(), and UseSlidingWindow(). { assert( ch.InheritsFrom("CandDigitListHandle") ); CandDeMuxDigitListHandle &cdlh = dynamic_cast<CandDeMuxDigitListHandle&>(ch); MSG("AlgDeMuxCosmics", Msg::kDebug) << cdlh.GetNDaughters() << " digits in AlgDeMuxCosmics" << endl; VldContext vldc = ch.GetCandRecord()->GetCandHeader()->GetVldContext(); if (vldc.GetDetector() != Detector::kFar) { static int msglimit = 20; // no more than 20 warnings about the detector if (msglimit) { MSG("DMX",Msg::kDebug) << "AlgDeMuxCosmics::RunAlg() can not DeMux " << Detector::AsString(vldc.GetDetector()) << " detector. Skip futher DeMux processing." << endl; if (--msglimit == 0) MSG("DMX",Msg::kDebug) << " ... last message of this type." << endl; } return; } //get the AlgConfigDeMux object fPlanesInSet = acd.GetInt("PlanesInSet"); fRequiredMatedSignalRatio = acd.GetDouble("RatioMatedSignalForValid"); fStrayCut = acd.GetInt("StrayDeltaStripLimit"); //MSG("DMX", Msg::kDebug) << "starting cosmics algorithm" << 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("DMXX", 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()); //instantiate the DmxInitialize object and then fill DmxStatus with the event information DmxUtilities util; util.FillPlaneArray(status, cdlh, acd); const TObjArray *planeArray = status->GetPlaneArray(); //create the TObjectItr over planes and program it to sort the planes DmxPlaneItr planeItr(planeArray); DmxPlaneKeyFunc *pnKF = planeItr.CreateKeyFunc(); pnKF->SetFun(KeyOnPlane); planeItr.GetSet()->AdoptSortKeyFunc(pnKF); pnKF = 0; planeItr.ResetFirst(); // while(planeItr.IsValid()){ // cout << dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber() << endl; // planeItr.Next(); // } // 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(); //make the iterators for the window, valid planes, and vertex DmxPlaneItr vertexPlaneItr(planeArray); DmxPlaneItr validPlaneItr(planeArray); DmxPlaneItr windowPlaneItr(planeArray); //create a KeyFunc to sort planes by number DmxPlaneKeyFunc *vertexPlaneKF = vertexPlaneItr.CreateKeyFunc(); DmxPlaneKeyFunc *validPlaneKF = validPlaneItr.CreateKeyFunc(); DmxPlaneKeyFunc *validWindowPlaneKF = windowPlaneItr.CreateKeyFunc(); //program the KeyFunc with the sort function vertexPlaneKF->SetFun(KeyOnPlane); validPlaneKF->SetFun(KeyOnPlane); validWindowPlaneKF->SetFun(KeyOnPlane); //get the NavSet from the iterator and pass the KeyFunc to it vertexPlaneItr.GetSet()->AdoptSortKeyFunc(vertexPlaneKF); validPlaneItr.GetSet()->AdoptSortKeyFunc(validPlaneKF); windowPlaneItr.GetSet()->AdoptSortKeyFunc(validWindowPlaneKF); //clear the KF pointer because we no longer own the KeyFunc vertexPlaneKF = 0; validPlaneKF = 0; validWindowPlaneKF = 0; planeItr.ResetFirst(); validPlaneItr.ResetFirst(); windowPlaneItr.ResetFirst(); //now program a key function to select on plane orientation - U first DmxPlaneKeyFunc *orientValidUKF = validPlaneItr.CreateKeyFunc(); DmxPlaneKeyFunc *orientValidWindowUKF = windowPlaneItr.CreateKeyFunc(); DmxPlaneKeyFunc *orientUKF = planeItr.CreateKeyFunc(); DmxPlaneKeyFunc *orientValidVKF = validPlaneItr.CreateKeyFunc(); DmxPlaneKeyFunc *orientValidWindowVKF = windowPlaneItr.CreateKeyFunc(); DmxPlaneKeyFunc *orientVKF = planeItr.CreateKeyFunc(); Int_t viewCtr = 0; while( viewCtr < 2 ){ fSlopeRMS = 0.; fStrayPlanes = 0; if( viewCtr == 0){ //program this function with the select function orientUKF->SetFun(KeyOnUView); orientValidUKF->SetFun(KeyOnValidU); orientValidWindowUKF->SetFun(KeyOnValidU); //adopt it as a selection function planeItr.GetSet()->AdoptSelectKeyFunc(orientUKF); orientUKF = 0; planeItr.Reset(); validPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidUKF); orientValidUKF = 0; validPlaneItr.Reset(); windowPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidWindowUKF); orientValidWindowUKF = 0; windowPlaneItr.Reset(); } else if(viewCtr == 1){ //program this function with the select function orientVKF->SetFun(KeyOnVView); orientValidVKF->SetFun(KeyOnValidV); orientValidWindowVKF->SetFun(KeyOnValidV); //adopt it as a selection function planeItr.GetSet()->AdoptSelectKeyFunc(orientVKF); orientVKF = 0; planeItr.Reset(); validPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidVKF); orientValidVKF = 0; validPlaneItr.Reset(); windowPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidWindowVKF); orientValidWindowVKF = 0; windowPlaneItr.Reset(); } // validPlaneItr.ResetLast(); // MSG("AlgDeMuxCosmics", Msg::kDebug)<<"event " << status->GetEventNumber() << endl; // while(validPlaneItr.IsValid()){ // MSG("DMX1", Msg::kDebug)<<"\tvalid plane " // << dynamic_cast<DmxPlane* >(validPlaneItr.Ptr())->GetPlaneNumber() << endl; // validPlaneItr.Prev(); // } // validPlaneItr.Reset(); // DmxPlane *plane = 0; // while( (plane = dynamic_cast<DmxPlane *>(validPlaneItr())) ){ // MSG("DMX1", Msg::kDebug) << "\t" << plane->GetPlaneNumber() <<endl; // } // validPlaneItr.Reset(); //slice to those planes from the vertex on planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber()); planeItr.GetSet()->Update(); validPlaneItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber()); validPlaneItr.GetSet()->Update(); windowPlaneItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber()); windowPlaneItr.GetSet()->Update(); //if only 2 planes set them to their best hyps if(validPlaneItr.SizeSelect() == 2){ validPlaneItr.ResetFirst(); //variables for straight line fit Float_t cogBeg = -1., zBeg = -1., cogEnd = -1., zEnd = -1.; while( validPlaneItr.IsValid() ){ //only set it if it is a shower plane, muon planes are set in the constructor if(validPlaneItr.Ptr()->GetPlaneType() == DmxPlaneTypes::kShower && !validPlaneItr.Ptr()->IsGolden()) dynamic_cast<DmxShowerPlane *>(validPlaneItr.Ptr())->SetStrips("best"); //cout << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetCoG() << "\t" // << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetPlaneNumber() // << "\t" << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetZPosition() << endl; //set the variables for the straighe line fit between the two planes if(cogBeg == -1.){ cogBeg = validPlaneItr.Ptr()->GetCoG(); zBeg = validPlaneItr.Ptr()->GetZPosition(); } else{ cogEnd = validPlaneItr.Ptr()->GetCoG(); zEnd = validPlaneItr.Ptr()->GetZPosition(); } validPlaneItr.Next(); }//end loop over planes //set the nonvalid planes if they exist if(planeItr.SizeSelect()>2 && cogEnd-cogBeg != 0. && cogBeg != -1. && cogEnd != -1.){ Double_t slope = (cogEnd-cogBeg)/(zEnd - zBeg); Double_t intercept = cogBeg - slope*zBeg; //cout << slope << "\t" << intercept << endl; SetPlanesToDeterminedFit(planeItr, intercept, slope); fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut); } //fStrayPlanes -= util.CheckFit(validPlaneItr); validPlaneItr.ResetFirst(); }//end if only 2 valid planes else if(validPlaneItr.SizeSelect() >2 ){ //if there are only enough planes for 1 instances of the sliding window + initial fit, //change fPlanesInSet to give you 3 + the initial fit, but at least 4 always if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) && fPlanesInSet==6) fPlanesInSet -= 2; else if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) && fPlanesInSet==5) fPlanesInSet -= 1; if(validPlaneItr.SizeSelect() > fPlanesInSet ){ UseSlidingWindow(planeItr, validPlaneItr, windowPlaneItr, status, status->GetVertexPlaneNumber(), status->GetEndPlaneNumber()); } else if(validPlaneItr.SizeSelect() <= fPlanesInSet){ UseSingleFit(planeItr, validPlaneItr, status); } //check the end points of the event - send the validPlaneItr to the function //fStrayPlanes -= util.CheckFit(validPlaneItr); //what about multiple muons? look for another vertex after the end plane planeItr.GetSet()->ClearSlice(false); validPlaneItr.GetSet()->ClearSlice(false); windowPlaneItr.GetSet()->ClearSlice(false); vertexPlaneItr.GetSet()->Slice(status->GetEndPlaneNumber()+1, 500); vertexPlaneItr.GetSet()->Update(); //MSG("DMX1", Msg::kDebug)<< "event = " << status->GetEventNumber() // << "\tend plane = " << status->GetEndPlaneNumber() << endl; Int_t nextVertex = -1; if(vertexPlaneItr.SizeSelect() > 3){ nextVertex = util.FindVertexPlane(vertexPlaneItr);} while( nextVertex != -1){ status->SetMultipleMuon(true); cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kMultipleMuonEvent); vertexPlaneItr.GetSet()->ClearSlice(false); vertexPlaneItr.GetSet()->Slice(nextVertex, 500); vertexPlaneItr.GetSet()->Update(); //get the next end plane Int_t endPlane = util.FindEndPlane(vertexPlaneItr); //MSG("DMX1", Msg::kDebug)<< "\tin double muon loop, vertex at plane " << nextVertex // << "\tend plane = " << endPlane << "\tplanes in slice = " // << vertexPlaneItr.SizeSelect() << endl; planeItr.GetSet()->Slice(nextVertex, endPlane); planeItr.GetSet()->Update(); validPlaneItr.GetSet()->Slice(nextVertex, endPlane); validPlaneItr.GetSet()->Update(); windowPlaneItr.GetSet()->Slice(nextVertex, endPlane); windowPlaneItr.GetSet()->Update(); //if only 2 planes set them to their best hyps if(validPlaneItr.SizeSelect() == 2){ validPlaneItr.ResetFirst(); //variables for straight line fit Float_t cogBeg = -1., zBeg = -1., cogEnd = -1., zEnd = -1.; while( validPlaneItr.IsValid() ){ //only set it if it is a shower plane, muon planes are set in the constructor if(validPlaneItr.Ptr()->GetPlaneType() == DmxPlaneTypes::kShower && !validPlaneItr.Ptr()->IsGolden()) dynamic_cast<DmxShowerPlane *>(validPlaneItr.Ptr())->SetStrips("best"); //set the variables for the straighe line fit between the two planes if(cogBeg == -1.){ cogBeg = validPlaneItr.Ptr()->GetCoG(); zBeg = validPlaneItr.Ptr()->GetZPosition(); } else{ cogEnd = validPlaneItr.Ptr()->GetCoG(); zEnd = validPlaneItr.Ptr()->GetZPosition(); } validPlaneItr.Next(); }//end loop over planes if(planeItr.SizeSelect()>2 && cogEnd-cogBeg != 0. && cogBeg != -1. && cogEnd != -1.){ Double_t slope = (cogEnd-cogBeg)/(zEnd - zBeg); Double_t intercept = cogBeg - slope*zBeg; SetPlanesToDeterminedFit(planeItr, intercept, slope); fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut); } //fStrayPlanes -= util.CheckFit(validPlaneItr); validPlaneItr.ResetFirst(); }//end if only 2 valid planes else if(validPlaneItr.SizeSelect() > 2){ //if there are only enough planes for 1 instances of the sliding window + initial fit, //change fPlanesInSet to give you 3 + the initial fit if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) && fPlanesInSet==6) fPlanesInSet -= 2; else if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) && fPlanesInSet==5) fPlanesInSet -= 1; if(validPlaneItr.SizeSelect() > fPlanesInSet ){ UseSlidingWindow(planeItr, validPlaneItr, windowPlaneItr, status, nextVertex, endPlane); } else if(validPlaneItr.SizeSelect() <= fPlanesInSet ){ UseSingleFit(planeItr, validPlaneItr, status); } }//end if more than 2 valid planes in event //check the end points of the event - send the validPlaneItr to the function //fStrayPlanes -= util.CheckFit(validPlaneItr); //look for another vertex after the end plane planeItr.GetSet()->ClearSlice(false); validPlaneItr.GetSet()->ClearSlice(false); windowPlaneItr.GetSet()->ClearSlice(false); vertexPlaneItr.GetSet()->ClearSlice(false); vertexPlaneItr.GetSet()->Slice(endPlane+1, 500); vertexPlaneItr.GetSet()->Update(); //MSG("DMX1", 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 }//end if more than 2 valid planes in set else if(validPlaneItr.SizeSelect() < 2){ MSG("DMX", Msg::kDebug)<< "Event " << status->GetEventNumber() << " not demuxed; only " << validPlaneItr.SizeSelect() << " valid planes in view" << endl; status->SetValidPlanesFailure(true); cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kTooFewValidPlanes); } //not enough planes to demux //done using the first view, so clear the selection function //MSG("Dmx", Msg::kDebug) <<"Done with View " << viewCtr << 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( status->GetFigureOfMeritFailure() ){ //see if the event failed the figure of Merit test, if so, see if it is an //overlapping muon cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kEventFailedFilter); if(util.IsOverlappingMultiple(validPlaneItr, 1.*status->GetVertexPlaneNumber(), acd.GetDouble("HoughInterceptRMSLimit"), acd.GetDouble("HoughSlopeRMSLimit"), acd.GetDouble("HoughPeakLimit"))){ status->SetUOverlappingMultiple(true); cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kMultipleMuonEvent); } } } if(viewCtr == 1 && status->GetEventDeMuxed() ){ status->SetVStrayPlanes(fStrayPlanes); status->SetVValidPlanes(validPlaneItr.SizeSelect()); cdlh.SetNumValidPlanesV(validPlaneItr.SizeSelect()); cdlh.SetNumStrayPlanesV(fStrayPlanes); //see if the event failed the figure of Merit test, if so, see if it is an //overlapping muon if( status->GetFigureOfMeritFailure() ){ cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kEventFailedFilter); if(util.IsOverlappingMultiple(validPlaneItr, 1.*status->GetVertexPlaneNumber(), acd.GetDouble("HoughInterceptRMSLimit"), acd.GetDouble("HoughSlopeRMSLimit"), acd.GetDouble("HoughPeakLimit"))){ status->SetVOverlappingMultiple(true); cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kMultipleMuonEvent); } } } vertexPlaneItr.GetSet()->ClearSlice(false); planeItr.GetSet()->ClearSlice(false); planeItr.GetSet()->AdoptSelectKeyFunc(0); planeItr.Reset(); validPlaneItr.GetSet()->ClearSlice(false); validPlaneItr.GetSet()->AdoptSelectKeyFunc(0); validPlaneItr.Reset(); windowPlaneItr.GetSet()->ClearSlice(false); windowPlaneItr.GetSet()->AdoptSelectKeyFunc(0); windowPlaneItr.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("DMX", 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 AlgDeMuxCosmics::SetPlanesToDeterminedFit (  DmxPlaneItr &  planeItr, Double_t  a, Double_t  b )  [private]

Definition at line 1751 of file AlgDeMuxCosmics.cxx. References DmxPlane::GetZPosition(), DmxPlane::IsGolden(), and DmxPlane::SetStrips(). { //MSG("DMXX", Msg::kDebug) << "in SetPlanesToDeterminedFit" << endl; planeItr.ResetFirst(); DmxPlane *plane = 0; while( planeItr.IsValid() ){ plane = planeItr.Ptr(); //if(plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){ //MSG("DMXX", Msg::kDebug)<<"\tin SetPlanesToDeterminedFit " << plane->GetPlaneNumber() // <<" has golden flag " << (Int_t)plane->IsGolden() << endl; Float_t fitCoG = a + (plane->GetZPosition() * b); //MSG("DMXX", Msg::kDebug)<< a << "\t" << b << "\t" << fitCoG << endl; //golden planes have already been set in their constructors if( !plane->IsGolden() ){ //MSG("DMXX", Msg::kDebug)<<"\tsetting plane " << plane->GetPlaneNumber() << endl; plane->SetStrips(fitCoG); } //}//end if plane is in the window bounds planeItr.Next(); } planeItr.Reset(); return; }

void AlgDeMuxCosmics::SetPlanesToDeterminedFit (  DmxPlaneItr &  planeItr, Double_t  a1, Double_t  a2, Double_t  a3, Double_t  a4, Float_t  offset )  [private]

Definition at line 1708 of file AlgDeMuxCosmics.cxx. References DmxPlane::GetZPosition(), DmxPlane::IsGolden(), and DmxPlane::SetStrips(). Referenced by RunAlg(), UseSingleFit(), and UseSlidingWindow(). { //MSG("DMXX", Msg::kDebug) << "in SetPlanesToDeterminedFit" << endl; planeItr.ResetFirst(); DmxPlane *plane = 0; while( planeItr.IsValid() ){ plane = planeItr.Ptr(); //MSG("DMXX", Msg::kDebug)<<"\tin SetPlanesToDeterminedFit " << plane->GetPlaneNumber() // <<" has golden flag " << (Int_t)plane->IsGolden() << endl; Float_t fitCoG = (a1 + (plane->GetZPosition()-offset)*a2 + TMath::Power(plane->GetZPosition()-offset,2)*a3 + TMath::Power(plane->GetZPosition()-offset,3)*a4); //MSG("DMXX", Msg::kDebug)<< a1 << "\t" << a2 << "\t" << a3 // << "\t" << a4 << "\t" << fitCoG // << "\t" << plane->GetZPosition() << endl; //golden planes have already been set in their constructors if( !plane->IsGolden() ){ //MSG("DMXX", Msg::kDebug)<<"\tsetting plane " << plane->GetPlaneNumber() // << " to " << fitCoG << endl; plane->SetStrips(fitCoG); } planeItr.Next(); } planeItr.Reset(); //cout << "finished set planes for fit" << endl; return; }

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

Reimplemented from AlgBase. Definition at line 551 of file AlgDeMuxCosmics.cxx. { }

void AlgDeMuxCosmics::UseSingleFit (  DmxPlaneItr &  planeItr, DmxPlaneItr &  validPlaneItr, DmxStatus *  status )  [private]

Definition at line 1087 of file AlgDeMuxCosmics.cxx. References FindCosmicSolution(), DmxUtilities::FindPlanesOffFit(), fStrayCut, fStrayPlanes, DmxPlane::GetCoG(), DmxStatus::GetEventNumber(), DmxPlane::GetZPosition(), DmxPlane::IsGolden(), MSG, DmxStatus::SetEventDirection(), DmxPlane::SetGolden(), DmxStatus::SetNonPhysicalFailure(), and SetPlanesToDeterminedFit(). Referenced by RunAlg(). { //MSG("DMXX", Msg::kDebug) << "in UseSingleFit" << endl; DmxUtilities util = DmxUtilities(); Double_t a1F = -10000.; Double_t a2F = -10000.; Double_t a3F = -10000.; Double_t a4F = -10000.; Double_t chiSqF = 1.1e20; Double_t a1B = -10000.; Double_t a2B = -10000.; Double_t a3B = -10000.; Double_t a4B = -10000.; Double_t chiSqB = 1.1e20; validPlaneItr.ResetFirst(); Float_t offset = validPlaneItr.Ptr()->GetZPosition(); validPlaneItr.ResetLast(); //set the event direction forwards status->SetEventDirection(1); //cout << "calling FindCosmicSolution" << endl; FindCosmicSolution(planeItr, validPlaneItr, status, a1F, a2F, a3F, a4F, chiSqF, offset); //set the event direction to backwards status->SetEventDirection(-1); FindCosmicSolution(planeItr, validPlaneItr, status, a1B, a2B, a3B, a4B, chiSqB, offset); DmxPlane *plane = 0; if( chiSqF <= chiSqB && a1F != -10000. && a2F != -10000. && a3F != -10000. && a4F != 10000.){ //loop over the planes and if a golden plane cog is more than 0.25m off from the fit //set it as non-golden. validPlaneItr.Reset(); while( (plane = planeItr()) ){ Double_t fitCoG = (a1F + a2F*plane->GetZPosition() + a3F*TMath::Power(plane->GetZPosition(),2) + a4F*TMath::Power(plane->GetZPosition(),3)); if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false); } SetPlanesToDeterminedFit(planeItr, a1F, a2F, a3F, a4F, offset); fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut); } else if(chiSqB < chiSqF && a1B != -10000. && a2B != -10000. && a3B != 10000. && a4B != 10000.){ //loop over the planes and if a golden plane cog is more than 0.25m off from the fit //set it as non-golden. validPlaneItr.Reset(); while( (plane = planeItr()) ){ Double_t fitCoG = (a1B + a2B*plane->GetZPosition() + a3B*TMath::Power(plane->GetZPosition(),2) + a4B*TMath::Power(plane->GetZPosition(),3)); if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false); } SetPlanesToDeterminedFit(planeItr, a1B, a2B, a3B, a4B, offset); fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut); } else if( a1B == -10000. && a2B == -10000. && a3B != 10000. && a4B != 10000. && a1F == -10000. && a2F == -10000. && a3F != -10000. && a4F != 10000.){ MSG("DMX", Msg::kDebug)<< "Event " << status->GetEventNumber() << " not demuxed; can't find a physical solution" << endl; status->SetNonPhysicalFailure(true); } return; }

void AlgDeMuxCosmics::UseSlidingWindow (  DmxPlaneItr &  planeItr, DmxPlaneItr &  validPlaneItr, DmxPlaneItr &  windowPlaneItr, DmxStatus *  status, Int_t  vertex, Int_t  endPlane )  [private]

Definition at line 1171 of file AlgDeMuxCosmics.cxx. References FindCosmicSolution(), DmxUtilities::FindPlanesOffFit(), FindWindowCosmicSolution(), fPlanesInSet, fStrayCut, fStrayPlanes, DmxPlane::GetCoG(), DmxStatus::GetEventNumber(), DmxPlane::GetPlaneNumber(), DmxStatus::GetVertexZPosition(), DmxPlane::GetZPosition(), DmxPlane::IsGolden(), MSG, DmxStatus::SetEventDirection(), DmxPlane::SetGolden(), DmxStatus::SetNonPhysicalFailure(), and SetPlanesToDeterminedFit(). Referenced by RunAlg(). { //MSG("DMXX", Msg::kDebug) << "in UseSlidingWindow" << endl; //start the sliding window - keep the window to 6 valid planes //and use the fit parameters you found for the previous set. if the //new plane is a shower plane, set to the hypothesis nearest the fit. //if none of the 3 best hypotheses is within 12 strips of the fit, then //just use the same fit parameters. redo the fit, and then set all //non-set planes in the window. move the the window down another plane //and do it again. DmxUtilities util; UInt_t validPlaneCnt = validPlaneItr.SizeSelect(); UInt_t loopCtr = 0; UInt_t ctr = 0; Int_t firstPlane = 0; Int_t almostLastPlane = 0; Int_t lastPlane = 0; Double_t a1 = -10000.; Double_t a2 = -10000.; Double_t a3 = -10000.; Double_t a4 = -10000.; Double_t chiSq = 1.1e20; Float_t offset = 0.; DmxPlane *plane = 0; Int_t direction = 1; //status->GetEventDirection(); // MSG("dmx", Msg::kDebug) << "Event = " << status->GetEventNumber() << "\tvertex = " << vertex // << "\tend = " << endPlane << endl; // while(validPlaneItr.IsValid() ){ // MSG("dmx", Msg::kDebug) << "\t" << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetPlaneNumber() // << "\t" << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->IsGolden() // <<endl; // validPlaneItr.Next(); // } validPlaneItr.ResetFirst(); windowPlaneItr.ResetFirst(); //now demux the event the <= makes sure that you find an initial fit for the first //set of planes, then go from there while(loopCtr <= (validPlaneCnt - fPlanesInSet) ){ ctr = 0; //MSG("dmx", Msg::kDebug) << "find new window" << endl; //MSG("dmx", Msg::kDebug) << "\tdirection = " << direction; while( validPlaneItr.IsValid() && ctr < fPlanesInSet){ plane = validPlaneItr.Ptr(); //MSG("dmx", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber() // << endl; if(direction == 1){ if( ctr == 0){ firstPlane = plane->GetPlaneNumber(); offset = plane->GetZPosition(); } if( ctr == fPlanesInSet-2){ almostLastPlane = plane->GetPlaneNumber(); } if( ctr == fPlanesInSet-1){ lastPlane = plane->GetPlaneNumber(); } validPlaneItr.Next(); } else if(direction == -1){ if( ctr == 0){ lastPlane = plane->GetPlaneNumber(); } if( ctr == fPlanesInSet-2){ almostLastPlane = plane->GetPlaneNumber(); } if( ctr == fPlanesInSet-1){ firstPlane = plane->GetPlaneNumber(); offset = plane->GetZPosition(); } validPlaneItr.Prev(); } ++ctr; } //MSG("dmx", Msg::kDebug) << "end find new window" << endl; ctr = 0; //MSG("dmx", Msg::kDebug) << "\tfirst plane = " << firstPlane << "\talmost last plane = " // << almostLastPlane << "\tlast plane = " << lastPlane << endl; //move the window windowPlaneItr.GetSet()->ClearSlice(false); planeItr.GetSet()->ClearSlice(false); windowPlaneItr.GetSet()->Slice(firstPlane, lastPlane); windowPlaneItr.GetSet()->Update(); planeItr.GetSet()->Slice(firstPlane, lastPlane); planeItr.GetSet()->Update(); if( loopCtr == 0){ //set the initial window - try it both ways and take the best one as //determined by the chiSq Double_t a1B = -10000.; Double_t a2B = -10000.; Double_t a3B = -10000.; Double_t a4B = -10000.; Double_t chiSqB = 1.1e20; //set the event direction forwards status->SetEventDirection(1); //you are looking at the first planes of the event, so your offset should be the //z position of the vertex Float_t offsetF = status->GetVertexZPosition(); FindCosmicSolution(planeItr, windowPlaneItr, status, a1, a2, a3, a4, chiSq, offsetF); //MSG("dmx", Msg::kDebug) << "\tforward direction" << a1 << "\t" << a2 // << "\t" << chiSq << "\t" << offset << endl; //set the event direction to backwards Int_t forwardLastPlane = lastPlane; status->SetEventDirection(-1); direction = -1; validPlaneItr.ResetLast(); windowPlaneItr.GetSet()->ClearSlice(false); planeItr.GetSet()->ClearSlice(false); windowPlaneItr.ResetLast(); ctr = 0; //MSG("dmx", Msg::kDebug) << "check reverse direction" << endl; while( validPlaneItr.IsValid() && ctr < fPlanesInSet){ plane = validPlaneItr.Ptr(); // MSG("dmx", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber() // << endl; if( ctr == 0){ lastPlane = plane->GetPlaneNumber(); } if( ctr == fPlanesInSet-2){ almostLastPlane = plane->GetPlaneNumber(); } if( ctr == fPlanesInSet-1){ firstPlane = plane->GetPlaneNumber(); offset = plane->GetZPosition(); } ++ctr; validPlaneItr.Prev(); } //MSG("dmx", Msg::kDebug) << "end check reverse direction" << endl; ctr = 0; //MSG("dmx", Msg::kDebug) << "first plane = " << firstPlane << "\talmost last plane = " // << almostLastPlane << "\tlast plane = " << lastPlane << endl; windowPlaneItr.GetSet()->Slice(firstPlane, lastPlane); windowPlaneItr.GetSet()->Update(); windowPlaneItr.ResetLast(); planeItr.GetSet()->Slice(firstPlane, lastPlane); planeItr.GetSet()->Update(); FindCosmicSolution(planeItr, windowPlaneItr, status, a1B, a2B, a3B, a4B, chiSqB, offset); //MSG("dmx", Msg::kDebug) << "\tbackward direction" << a1B << "\t" << a2B // << "\t" << chiSqB << "\t" << offset << endl; if( chiSq <= chiSqB && a1 != -10000. && a2 != -10000. && a3 != 10000. && a4 != 10000.){ planeItr.GetSet()->ClearSlice(false); //slice from the vertex to the last plane in the forward direction to get them //all planeItr.GetSet()->Slice(vertex, forwardLastPlane); planeItr.GetSet()->Update(); //loop over the planes and if a golden plane cog is more than 0.25m off from the fit //set it as non-golden. windowPlaneItr.GetSet()->ClearSlice(false); windowPlaneItr.GetSet()->Slice(vertex, forwardLastPlane); windowPlaneItr.GetSet()->Update(); while( (plane = windowPlaneItr()) ){ Double_t fitCoG = (a1 + a2*plane->GetZPosition() + a3*TMath::Power(plane->GetZPosition(),2) + a4*TMath::Power(plane->GetZPosition(),3)); if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false); } //MSG("dmx", Msg::kDebug) << "\tloop = " << loopCtr << "\tsetting planes " << vertex // << "\tto " << forwardLastPlane << endl; SetPlanesToDeterminedFit(planeItr, a1, a2, a3, a4, offsetF); fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut); status->SetEventDirection(1); direction = 1; validPlaneItr.ResetFirst(); ctr = 0; //get the validPlaneItr back to where it should be for this direction while(ctr<fPlanesInSet){ validPlaneItr.Next(); ++ctr; } ctr = 0; } else if(chiSqB < chiSq && a1B != -10000. && a2B != -10000. && a3B != 10000. && a4B != 10000.){ //set a, b, and chiSq to the backwards values a1 = a1B; a2 = a2B; a3 = a3B; a4 = a4B; chiSq = chiSqB; //slice from the end plane to the first plane planeItr.GetSet()->ClearSlice(false); planeItr.GetSet()->Slice(firstPlane, endPlane); planeItr.GetSet()->Update(); //loop over the planes and if a golden plane cog is more than 0.25m off from the fit //set it as non-golden. windowPlaneItr.ResetFirst(); while( (plane = windowPlaneItr()) ){ Double_t fitCoG = (a1 + a2*plane->GetZPosition() + a3*TMath::Power(plane->GetZPosition(),2) + a4*TMath::Power(plane->GetZPosition(),3)); if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false); } //MSG("dmx", Msg::kDebug) << "\tloop = " << loopCtr << "\tsetting planes " << firstPlane // << "\tto " << endPlane << endl; SetPlanesToDeterminedFit(planeItr, a1, a2, a3, a4, offset); fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut); //event direction is already backwards and the validPlaneItr is where it should be } else if( a1B == -10000. && a2B == -10000. && a1 == -10000. && a2 == -10000. && a3 == -10000. && a4 == -10000. && a3B == -10000. && a4B == -10000.){ MSG("dmx", Msg::kDebug)<< "Event " << status->GetEventNumber() << " not demuxed;" << " can't find a physical solution" << endl; status->SetNonPhysicalFailure(true); } //MSG("dmx", Msg::kDebug) << "\tloop 0 fit - " << a1 << "\t" << a2 << "\t" << a3 << "\tchi^2 = " // << chiSq << "\tstray planes = " << fStrayPlanes << endl; }//end if initial set else{ //set the window FindWindowCosmicSolution(windowPlaneItr, status, a1, a2, a3, a4, chiSq, offset); //MSG("dmx", Msg::kDebug) << "\tloop " << loopCtr << " " << firstPlane << " " << lastPlane // << "\t" << a1 << "\t" << a2 << "\t" << a3 << "\tchi^2 = " // << chiSq << "\tstray planes = " << fStrayPlanes << endl; planeItr.GetSet()->ClearSlice(false); if(direction == 1 ){ firstPlane = almostLastPlane + 1; //if its the last loop, select the rest of the unset planes to set, otherwise //stay in the loop bounds if(loopCtr == (validPlaneItr.SizeSelect()-fPlanesInSet)) lastPlane = endPlane; } else if(direction == -1 ){ lastPlane = almostLastPlane-1; if(loopCtr == (validPlaneItr.SizeSelect()-fPlanesInSet)) firstPlane = vertex; } planeItr.GetSet()->Slice(firstPlane, lastPlane); planeItr.GetSet()->Update(); //MSG("dmx", Msg::kDebug) << "\tloop = " << loopCtr << "\tsetting planes " << firstPlane // << "\tto " << lastPlane << endl; SetPlanesToDeterminedFit(planeItr, a1, a2, a3, a4, offset); fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut); //status->SetHistogramEntries(b, a, plane->GetPlaneView()); } ctr = 0; //set the validPlaneItr to the opposite direction and back up 5 planes //MSG("dmx", Msg::kDebug) << "move iterator to first plane in new window" << endl; while(ctr<fPlanesInSet-1 && validPlaneItr.IsValid()){ //MSG("dmx", Msg::kDebug) << "plane number = " // << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetPlaneNumber() // << endl; ++ctr; if(direction == 1){ validPlaneItr.Prev(); } else if(direction == -1){ validPlaneItr.Next(); } } //MSG("dmx", Msg::kDebug) << "end move iterator to first plane in new window" << endl; ctr = 0; //should go through the loop and get the new window ++loopCtr; } return; }

---

Member Data Documentation

UInt_t AlgDeMuxCosmics::fPlanesInSet [private]

Definition at line 40 of file AlgDeMuxCosmics.h. Referenced by RunAlg(), and UseSlidingWindow().

Double_t AlgDeMuxCosmics::fRequiredMatedSignalRatio [private]

Definition at line 37 of file AlgDeMuxCosmics.h. Referenced by RunAlg().

Double_t AlgDeMuxCosmics::fSlopeRMS [private]

Definition at line 38 of file AlgDeMuxCosmics.h. Referenced by RunAlg().

Int_t AlgDeMuxCosmics::fStrayCut [private]

Definition at line 41 of file AlgDeMuxCosmics.h. Referenced by RunAlg(), UseSingleFit(), and UseSlidingWindow().

Int_t AlgDeMuxCosmics::fStrayPlanes [private]

Definition at line 39 of file AlgDeMuxCosmics.h. Referenced by RunAlg(), UseSingleFit(), and UseSlidingWindow().

---

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

• AlgDeMuxCosmics.h
• AlgDeMuxCosmics.cxx

---