AlgDeMuxBeam Class Reference

#include <AlgDeMuxBeam.h>

Inheritance diagram for AlgDeMuxBeam:
AlgBase

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

Detailed Description

Definition at line 23 of file AlgDeMuxBeam.h.


Constructor & Destructor Documentation

AlgDeMuxBeam::AlgDeMuxBeam (  ) 

Definition at line 72 of file AlgDeMuxBeam.cxx.

00072                            :
00073         fMaxResidual(0.1),
00074         fMaxMuonRemovalFraction(0.2),
00075         fMuonPlanesInSet(5),
00076         fMuonStartPlaneNumber(500),
00077         fStrayCut(12),
00078         fStrayPlanes(0)
00079 {
00080   //default constructor
00081 }

AlgDeMuxBeam::~AlgDeMuxBeam (  )  [virtual]

Definition at line 85 of file AlgDeMuxBeam.cxx.

00086 {
00087 }


Member Function Documentation

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

Definition at line 2008 of file AlgDeMuxBeam.cxx.

References FindFit(), fPlanesInSet, DmxPlane::GetCoG(), DmxPlane::GetPlaneType(), DmxPlane::GetZPosition(), id, Msg::kDebug, DmxPlaneTypes::kShower, MSG, SetFirstNPlanes(), and DmxPlane::SetStrips().

02010 {
02011         
02012         Int_t hypsToUse[] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0} ;
02013         Int_t leverArm = (Int_t)validPlaneItr.SizeSelect();
02014         if(leverArm > fPlanesInSet) leverArm = fPlanesInSet;
02015         //find the bounds for the first n planes
02016         
02017         MSG("AlgDeMuxBeam", Msg::kDebug) << "In DeMuxFirstNPlanes with lever arm = " 
02018                                                                 << leverArm << endl;
02019         
02020         //useA2 is the slope, useA1 is the intercept
02021         
02022         //variables to keep track of the best reconstruction
02023         Double_t useA1 = -10000.;
02024         Double_t useA2 = -10000.;
02025         Double_t useA3 = -10000.;
02026         Double_t useA4 = -10000.;
02027         Double_t bestChi2 = 1.e20;
02028         
02029         cdhit.ResetFirst();
02030         validPlaneItr.ResetFirst();
02031         planeItr.ResetFirst();
02032         
02033         if(leverArm  == 2){
02034                 
02035                 //only two planes in the set.  if they are shower planes, just set it 
02036                 //each plane to the best hypothesis.  valid muon planes are already set.
02037                 //get the slope and intercept for a straight line through the two and 
02038                 //then set all remaining planes based on that
02039                 
02040                 Float_t z1 = 0., z2 = 0.;
02041                 Float_t tPos1 = 0., tPos2 = 0.;
02042                 Int_t ctr = 0;
02043                 DmxPlane *plane = 0;
02044                 while( (plane = validPlaneItr()) ){
02045                         
02046                         if(plane->GetPlaneType() == DmxPlaneTypes::kShower) plane->SetStrips();
02047                         if(ctr == 0){
02048                                 z1 = plane->GetZPosition();
02049                                 tPos1 = plane->GetCoG();
02050                         }
02051                         else{
02052                                 z2 = plane->GetZPosition();
02053                                 tPos2 = plane->GetCoG();
02054                         }
02055                         ++ctr;
02056                 }//end loop to find the first and last plane tpos and zpos
02057                 
02058                 if(z1 != z2) useA2 = (tPos2 - tPos1)/(z2 - z1);
02059                 useA1 = tPos2 - useA2*z2;
02060         }//end if only 2 valid planes
02061         else if(leverArm==3){
02062                 //evaluate the rms of the centers of gravity - take the tightest
02063                 //distribution
02064                 for(Int_t ia = 1; ia < 4; ++ia){
02065                         hypsToUse[0] = ia;
02066                         for(Int_t ib = 1; ib < 4; ++ib){
02067                                 hypsToUse[1] = ib;
02068                                 for(Int_t ic = 1; ic < 4; ++ic){
02069                                         hypsToUse[2] = ic;
02070                                         
02071                                         FindFit(planeItr, validPlaneItr, useA1, useA2, useA3, useA4, bestChi2,
02072                                                         leverArm, hypsToUse);
02073                                 }
02074                         }
02075                 }
02076         }//end if lever arm = 3 
02077         else if(leverArm == 4){
02078                 for(Int_t ia = 1; ia < 4; ++ia){
02079                         hypsToUse[0] = ia;
02080                         for(Int_t ib = 1; ib < 4; ++ib){
02081                                 hypsToUse[1] = ib;
02082                                 for(Int_t ic = 1; ic < 4; ++ic){
02083                                         hypsToUse[2] = ic;
02084                                         for(Int_t id = 1; id < 4; ++id){
02085                                                 hypsToUse[3] = id;
02086                                                 
02087                                                 FindFit(planeItr, validPlaneItr, useA1, useA2, useA3, useA4, bestChi2,
02088                                                                 leverArm, hypsToUse);
02089                                                 
02090                                         }
02091                                 }
02092                         }
02093                 }
02094         }//end if 4 plane lever arm
02095         else if(leverArm == 5){
02096                 for(Int_t ia = 1; ia < 4; ++ia){
02097                         hypsToUse[0] = ia;
02098                         for(Int_t ib = 1; ib < 4; ++ib){
02099                                 hypsToUse[1] = ib;
02100                                 for(Int_t ic = 1; ic < 4; ++ic){
02101                                         hypsToUse[2] = ic;
02102                                         for(Int_t id = 1; id < 4; ++id){
02103                                                 hypsToUse[3] = id;
02104                                                 for(Int_t ie = 1; ie < 4; ++ie){
02105                                                         hypsToUse[4] = ie;
02106                                                         
02107                                                         FindFit(planeItr, validPlaneItr, useA1, useA2, useA3, useA4, bestChi2,
02108                                                                         leverArm, hypsToUse);
02109                                                         
02110                                                 }
02111                                         }
02112                                 }
02113                         }
02114                 }
02115         }
02116         else if(leverArm == 6){
02117                 for(Int_t ia = 1; ia < 4; ++ia){
02118                         hypsToUse[0] = ia;
02119                         for(Int_t ib = 1; ib < 4; ++ib){
02120                                 hypsToUse[1] = ib;
02121                                 for(Int_t ic = 1; ic < 4; ++ic){
02122                                         hypsToUse[2] = ic;
02123                                         for(Int_t id = 1; id < 4; ++id){
02124                                                 hypsToUse[3] = id;
02125                                                 for(Int_t ie = 1; ie < 4; ++ie){
02126                                                         hypsToUse[4] = ie;
02127                                                         for(Int_t ig = 1; ig < 4; ++ig){
02128                                                                 hypsToUse[5] = ig;
02129                                                                 
02130                                                                 FindFit(planeItr, validPlaneItr, useA1, useA2, useA3, useA4, bestChi2,
02131                                                                                 leverArm, hypsToUse);
02132                                                         }
02133                                                 }
02134                                         }
02135                                 }
02136                         }
02137                 }
02138         }//end if lever arm = 6 planes
02139         
02140         MSG("AlgDeMuxBeam", Msg::kDebug) << "\tfinal fit" << "\t" << useA1 << "\t" << useA2 
02141                                                                 << "\t" << bestChi2 << endl; 
02142         
02143         SetFirstNPlanes(planeItr, useA2, useA1);
02144         
02145         return;
02146 }

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

Definition at line 646 of file AlgDeMuxBeam.cxx.

References fFinalShowerRegionIntercept, fFinalShowerRegionSlope, FindFitOverNPlanes(), fMaxStripAdjustment, fNumberOfHypotheses, fPlanesInSet, 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::kDebug, Detector::kFar, DmxPlaneTypes::kShower, MSG, and SetFirstNPlanes().

Referenced by RunAlg().

00648 {
00649         Int_t leverArm = (Int_t)validPlaneItr.SizeSelect();
00650         if(leverArm > fPlanesInSet) leverArm = fPlanesInSet;
00651         Int_t planeCtr = 0;
00652         
00653         //loop over the valid planes and take the first leverArm planes.
00654         DmxPlane *plane = 0;
00655         DmxShowerPlane *showerPlane = 0;
00656         DmxMuonPlane *muonPlane = 0;
00657         
00658         validPlaneItr.ResetFirst();
00659         
00660         Float_t matedCharge = 0.;
00661         Float_t planeCharge = 0.;
00662         Float_t maxMatedCharge = 0.;
00663         Int_t maxMatedChargeHyp = 0;
00664         Float_t zPos[] = {0., 0., 0., 0., 0., 0.};
00665         Float_t tPos[] = {0., 0., 0., 0., 0., 0.};
00666         Float_t weight[] = {1., 1., 1., 1., 1., 1.};
00667         Float_t residual[] = {0., 0., 0., 0., 0., 0.};
00668         Int_t planeNum[] = {0,0,0,0,0,0};
00669         Float_t slope = 0.;
00670         Float_t intercept = 0.;
00671         Float_t averageRMS = 0.;
00672         Int_t rmsCtr = 0;
00673         
00674         //now loop over all hypotheses in turn for the planes and see which set gives you 
00675         //the most mated signal over the leverArm planes
00676         for(Int_t i = 0; i<fNumberOfHypotheses; ++i){
00677                 
00678                 validPlaneItr.ResetFirst();
00679                 planeCtr = 0;
00680                 matedCharge = 0.;
00681                 rmsCtr = 0;
00682                 averageRMS = 0.;
00683                 
00684                 while( (plane = validPlaneItr()) && planeCtr<leverArm){
00685                         
00686                         planeCharge = plane->GetPlaneCharge();
00687 
00688                         //check out the type of plane you have here
00689                         if(plane->GetPlaneType() == DmxPlaneTypes::kShower){                            
00690                                 matedCharge += dynamic_cast<DmxShowerPlane *>(plane)->GetHypothesis(i)->GetMatedSignalRatio()*planeCharge;
00691                                 averageRMS += dynamic_cast<DmxShowerPlane *>(plane)->GetHypothesis(i)->GetTieBreakerStat();
00692                                 ++rmsCtr;
00693                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "plane = " << plane->GetPlaneNumber() 
00694                                                                                                 << " i = " << i 
00695                                                                                                 << " mated frac = " << dynamic_cast<DmxShowerPlane *>(plane)->GetHypothesis(i)->GetMatedSignalRatio()
00696                                                                                                 << " rms = " << dynamic_cast<DmxShowerPlane *>(plane)->GetHypothesis(i)->GetTieBreakerStat()
00697                                                                                                 << " charge = " << planeCharge << endl;
00698                         }               
00699                                 
00700                         //its a valid muon plane so check to see if the current hypothesis contains the initial strip
00701                         else{
00702                                 if(dynamic_cast<DmxMuonPlane *>(plane)->GetInitialStripCoG()>=i
00703                                    && dynamic_cast<DmxMuonPlane *>(plane)->GetInitialStripCoG()<=i+23)
00704                                         matedCharge += planeCharge;                             
00705                         }
00706                         
00707                         ++planeCtr;
00708                 }//end loop over first leverArm valid planes
00709                 
00710                 //does the current hypothesis represent the most mated signal?
00711                 if(matedCharge>maxMatedCharge){
00712                         maxMatedCharge = matedCharge;
00713                         //minAverageRMS = averageRMS/(1.*rmsCtr);
00714                         maxMatedChargeHyp = i;
00715                 }
00716                 
00717         }//end loop over hypotheses
00718         
00719         MSG("AlgDeMuxBeam", Msg::kDebug) << "best initial hypothesis = " << maxMatedChargeHyp << endl;
00720         
00721         validPlaneItr.ResetFirst();
00722         planeCtr = 0;
00723 
00724         //make an array holding the best hypothesis numbers for these planes
00725         Int_t bestHyps[] = {maxMatedChargeHyp, maxMatedChargeHyp, maxMatedChargeHyp, maxMatedChargeHyp, maxMatedChargeHyp, maxMatedChargeHyp};
00726 
00727         //get the transverse position for the best hypothesis
00728         UgliGeomHandle ugh(fVldContext);
00729         
00730         //check the number of valid planes - if only 2, then just set these suckers to the best hyp and be done with it.
00731         if(leverArm == 2){
00732                 //get the strip end id for the center of gravity for the center of the best hypothesis
00733                 PlexStripEndId seid(Detector::kFar,dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetPlaneNumber(),maxMatedChargeHyp+12);
00734                 UgliStripHandle ush = ugh.GetStripHandle(seid);
00735                 
00736                 intercept = ush.GetTPos();
00737                                 
00738                 //set the first n planes - slope is 0 because you are just setting all the planes to one
00739                 //hypothesis
00740                 SetFirstNPlanes(planeItr, 0., intercept);
00741 
00742                 return;
00743         }
00744         
00745         //loop over the valid planes again and see if by adjusting the hypothesis number you can get more 
00746         //mated signal for the plane
00747         while ((plane = validPlaneItr()) && planeCtr<leverArm){
00748                 maxMatedCharge = 0.;
00749                 maxMatedChargeHyp = bestHyps[planeCtr];
00750                 zPos[planeCtr] = plane->GetZPosition();
00751                 //get the strip end id for the center of gravity for the center of the best hypothesis
00752                 PlexStripEndId seid(Detector::kFar,plane->GetPlaneNumber(),bestHyps[planeCtr]+12);
00753                 UgliStripHandle ush = ugh.GetStripHandle(seid);
00754                 
00755                 tPos[planeCtr] = ush.GetTPos();
00756                 planeNum[planeCtr] = plane->GetPlaneNumber();
00757                 
00758                 planeCharge = plane->GetPlaneCharge();
00759                 showerPlane = 0;
00760                 muonPlane = 0;
00761                 if( plane->GetPlaneType() == DmxPlaneTypes::kShower) showerPlane = dynamic_cast<DmxShowerPlane *>(plane);
00762                 else muonPlane = dynamic_cast<DmxMuonPlane *>(plane);
00763                 
00764                 if(showerPlane){
00765                         tPos[planeCtr] = showerPlane->GetHypothesisCoG(maxMatedChargeHyp);
00766                         maxMatedCharge = showerPlane->GetHypothesis(bestHyps[planeCtr])->GetMatedSignalRatio()*planeCharge;
00767                 }
00768                 else{
00769                         if(muonPlane->GetInitialStripCoG()>=maxMatedChargeHyp 
00770                            && muonPlane->GetInitialStripCoG()<=maxMatedChargeHyp+23){
00771                                 tPos[planeCtr] = muonPlane->GetCoG();
00772                                 maxMatedCharge = planeCharge;                           
00773                         }
00774                 }
00775 
00776                 //look at hypotheses n strips on either side of the chosen hyp for the plane
00777                 for(Int_t i = 1; i<=fMaxStripAdjustment; ++i){
00778                         
00779                         if(showerPlane){
00780                                 if(bestHyps[planeCtr]-i>=0){
00781                                         matedCharge = showerPlane->GetHypothesis(bestHyps[planeCtr]-i)->GetMatedSignalRatio()*planeCharge;      
00782                                         if(matedCharge>maxMatedCharge){
00783                                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "initial best hyp = " << bestHyps[planeCtr]
00784                                                                                                                 << " new best hyp = " << bestHyps[planeCtr]-i 
00785                                                                                                                 << " new tpos = " << showerPlane->GetHypothesisCoG(maxMatedChargeHyp)
00786                                                                                                                 << endl;
00787                                                 maxMatedCharge = matedCharge;
00788                                                 maxMatedChargeHyp = bestHyps[planeCtr]-i;
00789                                                 tPos[planeCtr] = showerPlane->GetHypothesisCoG(maxMatedChargeHyp);
00790                                         }                                       
00791                                 }//end if the new hypothesis is allowed
00792                                 if(bestHyps[planeCtr]+i<fNumberOfHypotheses){
00793                                         matedCharge = showerPlane->GetHypothesis(bestHyps[planeCtr]+i)->GetMatedSignalRatio()*planeCharge;      
00794                                         if(matedCharge>maxMatedCharge){
00795                                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "initial best hyp = " << bestHyps[planeCtr]
00796                                                                                                                 << " new best hyp = " << bestHyps[planeCtr]+i   
00797                                                                                                                 << " new tpos = " << showerPlane->GetHypothesisCoG(maxMatedChargeHyp)<< endl;
00798                                                 maxMatedCharge = matedCharge;
00799                                                 maxMatedChargeHyp = bestHyps[planeCtr]+i;
00800                                                 tPos[planeCtr] = showerPlane->GetHypothesisCoG(maxMatedChargeHyp);
00801                                         }                                       
00802                                 }//end if the other new hypothesis is allowed
00803                         }//end if it is a shower plane
00804                         if(muonPlane){
00805                                 if(bestHyps[planeCtr]+i<fNumberOfHypotheses){
00806                                         if(muonPlane->GetInitialStripCoG()>=bestHyps[planeCtr]+i && muonPlane->GetInitialStripCoG()<=bestHyps[planeCtr]+i+23){
00807                                                 //this hypothesis contains the original decoding for the valid muon plane
00808                                                 //so the mated charge is all charge on the plane and the transverse position is 
00809                                                 //the original cog
00810                                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "initial best hyp = " << bestHyps[planeCtr]
00811                                                                                                                 << " new best hyp = " << bestHyps[planeCtr]+i  
00812                                                                                                                 << " new tpos = " << muonPlane->GetCoG()<< endl;
00813                                                 maxMatedChargeHyp = bestHyps[planeCtr] + i;
00814                                                 maxMatedCharge = planeCharge;                           
00815                                                 tPos[planeCtr] = muonPlane->GetInitialCoG();
00816                                         }
00817                                 }//end if the new hypothesis is allowed
00818                                 if(bestHyps[planeCtr]-i>=0){
00819                                         if(muonPlane->GetInitialStripCoG()>=bestHyps[planeCtr]-i && muonPlane->GetInitialStripCoG()<=bestHyps[planeCtr]-i+23){
00820                                                 //this hypothesis contains the original decoding for the valid muon plane
00821                                                 //so the mated charge is all charge on the plane and the transverse position is 
00822                                                 //the original cog
00823                                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "initial best hyp = " << bestHyps[planeCtr]
00824                                                                                                                 << " new best hyp = " << bestHyps[planeCtr]-i  
00825                                                                                                                 << " new tpos = " << muonPlane->GetCoG() << endl;
00826                                                 maxMatedChargeHyp = bestHyps[planeCtr] - i;
00827                                                 maxMatedCharge = planeCharge;   
00828                                                 tPos[planeCtr] = muonPlane->GetInitialCoG();
00829                                         }
00830                                 }//end if the other new hypothesis is allowed
00831                         }//end if it is a muon plane
00832                 }//end loop over other hypotheses
00833                 
00834                 bestHyps[planeCtr] = maxMatedChargeHyp;
00835                 
00836                 MSG("AlgDeMuxBeam", Msg::kDebug) << plane->GetPlaneNumber() << " zpos = " << zPos[planeCtr]
00837                                                                                 << " tpos = " << tPos[planeCtr] << " best hyp = "
00838                                                                                 << bestHyps[planeCtr] << endl;
00839                 
00840                 ++planeCtr;
00841         }//end loop over valid planes to adjust the fit
00842                 
00843         FindFitOverNPlanes(zPos,tPos,weight,residual,slope, intercept, leverArm, planeNum);
00844         
00845         MSG("AlgDeMuxBeam", Msg::kDebug) << "set to slope = " << slope << " intercept = " 
00846                                                                         << intercept << endl;
00847         
00848         fFinalShowerRegionSlope = slope;
00849         fFinalShowerRegionIntercept = intercept;
00850         
00851         //set the first n planes
00852         SetFirstNPlanes(planeItr, slope, intercept);
00853         
00854         return;
00855 }

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

Definition at line 1717 of file AlgDeMuxBeam.cxx.

References fVldContext, UgliGeomHandle::GetStripHandle(), and UgliStripHandle::GetTPos().

01719 {
01720         
01721         //   MSG("AlgDeMuxBeam", Msg::kDebug) << "In FindDigitCoG" << endl;
01722         
01723         Float_t cog = -5.;
01724         Float_t sigTPosSum = 0.;
01725         Float_t sigSum = 0.;
01726         
01727         UgliGeomHandle  ugh(fVldContext);
01728         UgliStripHandle ush;
01729         
01730         validPlaneItr.ResetFirst();
01731         
01732         //loop over the valid planes and get the digits from those planes
01733         while(validPlaneItr.IsValid()){
01734                 cdhit.GetSet()->Slice(validPlaneItr.Ptr()->GetPlaneNumber());
01735                 
01736                 //got to tell the iterator to sort the slice
01737                 cdhit.GetSet()->Update();
01738                 
01739                 //     MSG("AlgDeMuxBeam", Msg::kDebug) << "getting CoG with plane  " 
01740                 //                             <<  validPlaneItr.Ptr()->GetPlaneNumber()
01741                 //                             << " digits" << endl;
01742                 
01743                 while(cdhit.IsValid()){
01744                         if(cdhit.Ptr()->GetDeMuxDigitFlagWord() == 0){
01745                                 ush = ugh.GetStripHandle(cdhit.Ptr()->GetPlexSEIdAltL().GetBestSEId());
01746                                 sigTPosSum += ush.GetTPos()*cdhit.Ptr()->GetCharge();
01747                                 sigSum += cdhit.Ptr()->GetCharge();
01748                         }
01749                         cdhit.Next();
01750                 }//end loop over digits
01751                 
01752                 //clear the digit slice
01753                 cdhit.GetSet()->ClearSlice(false);
01754                 
01755                 //advance the plane iterator
01756                 validPlaneItr.Next();
01757         }//end loop over planes
01758         
01759         validPlaneItr.ResetFirst();
01760         cdhit.ResetFirst();
01761         
01762         //find the cog
01763         if(sigSum > 0.) cog = sigTPosSum/sigSum;
01764         
01765         //   MSG("AlgDeMuxBeam", Msg::kDebug) << "digit cog = " << cog << endl;
01766         return cog;
01767 }

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().

02153 {
02154         
02155         Double_t fitA1 = 0.;
02156         Double_t fitA2 = 0.;
02157         Double_t fitA3 = 0.;
02158         Double_t fitA4 = 0.;
02159         Double_t chi2 = 0.;
02160         DmxUtilities util;
02161         
02162         FindFitForHypothesisSet(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypotheses, 
02163                                                         leverArm);
02164         
02165         if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, 
02166                                                                                   fitA3, fitA4) ){
02167                 bestChi2 = chi2;
02168                 a1 = fitA1;
02169                 a2 = fitA2;
02170                 a3 = fitA3;
02171                 a4 = fitA4;
02172                 
02173                 //see if you can drop some planes and improve the fit
02174                 if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, 
02175                                                                         fitA4, leverArm) ){
02176                         
02177                         FindFitForHypothesisSet(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, 
02178                                                                         hypotheses, leverArm);
02179                         
02180                         if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, 
02181                                                                                                   fitA3, fitA4) ){
02182                                 bestChi2 = chi2;
02183                                 a1 = fitA1;
02184                                 a2 = fitA2;
02185                                 a3 = fitA3;
02186                                 a4 = fitA4;
02187                         }
02188                         validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
02189                 }//end if planes should be dropped because they are way off
02190         }//end if its a better straight line fit
02191         
02192         return;
02193 }

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(), Msg::kDebug, DmxPlaneTypes::kShower, and MSG.

Referenced by FindFit().

02201 {
02202         
02203         MSG("AlgDeMuxBeam", Msg::kDebug) << "in FindFitForHypothesisSet" << endl;
02204         
02205         planeItr.ResetFirst();
02206         
02207         //cout << "in find fit for first n planes" << endl;
02208         
02209         //declare the variables to do a least squares fit to a straight line.  
02210         //weight is the inverse fraction of the total charge on the plane that 
02211         //is on opposite ends of common strips. multiply the inverse by a 
02212         //sensible number to take account for digits with 1000+ adc's
02213         
02214         //the arrays are the number of planes in the set, ie the lever arm
02215         Double_t x[] = {0.,0.,0.,0.,0.,0.};
02216         Double_t y[] = {0.,0.,0.,0.,0.,0.};
02217         Double_t weight[] = {1.,1.,1.,1.,1.,1.};
02218         Double_t chiSq = 0.;
02219         
02220         DmxUtilities util;
02221         
02222         Int_t planeCtr = 0;
02223         Int_t arrayCtr = 0;
02224         DmxPlane *plane = 0;
02225         
02226         //loop over the plane iterator to fill the variables
02227         
02228         while( (plane = planeItr()) && planeCtr < leverArm){
02229                 
02230                 //get some generic info about the plane
02231                 MSG("AlgDeMuxBeam", Msg::kDebug) << "plane " << plane->GetPlaneNumber() 
02232                 << " is valid " << (Int_t)plane->IsValid()
02233                 << " is golden " << (Int_t)plane->IsGolden() 
02234                 << " masked " << (Int_t)planeItr.GetSet()->GetMasks().GetMask(plane)
02235                 << " " << DmxPlaneTypes::GetPlaneTypeString(plane->GetPlaneType())
02236                 << endl;
02237                 
02238                 if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){
02239                         x[arrayCtr] = plane->GetZPosition();
02240                         
02241                         //get the maximum possible weight first, then for each hypothesis 
02242                         //multiply by the square of the fraction
02243                         weight[arrayCtr] = 1./ TMath::Sqrt(plane->GetPlaneCharge());
02244                         if(plane->IsGolden()) weight[arrayCtr] /= TMath::Sqrt(10.);
02245                         
02246                         y[arrayCtr] = plane->GetCoG();
02247                         
02248                         //only do this if the shower plane is not golden - if it is, 
02249                         //you know where it should be reconstructed to
02250                         if( plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){
02251                                 
02252                                 if(hypotheses[planeCtr] == 1){
02253                                         weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()); 
02254                                 }
02255                                 //plane->GetCoG() returns the best hypothesis CoG so only change the value 
02256                                 //of y if you are wanting the 2nd or 3rd best hypotheses
02257                                 else if(hypotheses[planeCtr] == 2){ 
02258                                         y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
02259                                         weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio());
02260                                 }
02261                                 else if(hypotheses[planeCtr] == 3){ 
02262                                         y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
02263                                         weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());
02264                                 }
02265                         }
02266                         
02267                         MSG("AlgDeMuxBeam", Msg::kDebug) << "zpos = " << x[arrayCtr] << " tpos = " << y[arrayCtr]
02268                                 << " hypothesis = " << hypotheses[arrayCtr] << endl;
02269                         
02270                         ++arrayCtr;
02271                         
02272                 }//end if the plane isnt marked as sucking && it is in the window bounds
02273                 
02274                 ++planeCtr;
02275         } //end loop over planes to find variables for fit
02276         
02277         //cout << "filled arrays" << endl;
02278         
02279         //use the DmxUtilities function to find a linear fit
02280         util.FindLinearFit(x, y, weight, arrayCtr, a1, a2, chiSq);
02281         
02282         prevChiSq = (double)chiSq;
02283         //a1 = intercept;
02284         //a2 = slope;   
02285         MSG("AlgDeMuxBeam", Msg::kDebug) << "\tinitial fit " << a1 << "\t" << a2 << "\t" 
02286                                                                 << a3 << "\t" << a4 << "\t" << chiSq << endl; 
02287         
02288         planeItr.Reset();
02289         planeCtr = 0;
02290         arrayCtr = 0;
02291         //cout << "reset the iterator" << endl;
02292         
02293         //redo the fit as many times as there are planes, dropping each plane in 
02294         //turn to see if it produces a better chi^2 value - only do it if you have 
02295         //more than 3 planes - ie you can drop one and still do a reasonable fit
02296         //2 planes = great straight line every time.
02297         if(leverArm > 3){
02298                 
02299                 Double_t fitA1 = 0.;
02300                 Double_t fitA2 = 0.;
02301                 Double_t fitA3 = 0.;
02302                 Double_t fitA4 = 0.;
02303                 
02304                 for(Int_t i = 0; i < leverArm; i++){
02305                         arrayCtr = 0;
02306                         
02307                         //cout <<"in loop " << i << endl;
02308                         //loop over the plane iterator to fill the variables
02309                         while( (plane = planeItr()) && planeCtr < leverArm){
02310                                 
02311                                 //use those planes not marked as kTRUE - a mark of kTRUE means to drop the plane
02312                                 //from the set when finding the fit. ie if its true, the plane truely sucks
02313                                 if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){
02314                                         
02315                                         if(i != planeCtr || plane->IsGolden() ){
02316                                                 x[arrayCtr] = plane->GetZPosition();
02317                                                 y[arrayCtr] = plane->GetCoG();
02318                                                 
02319                                                 weight[arrayCtr] = 1. / TMath::Sqrt(plane->GetPlaneCharge());
02320                                                 if(plane->IsGolden()) weight[arrayCtr] /= TMath::Sqrt(10.);
02321                                                 
02322                                                 //only make changes if not a golden plane
02323                                                 if( plane->GetPlaneType() == DmxPlaneTypes::kShower  && !plane->IsGolden()){
02324                                                         //only change y if the hypothesis isnt the best one - the unset plane 
02325                                                         //returns the best cog from the GetCoG() method
02326                                                         
02327                                                         if(hypotheses[planeCtr] == 1){
02328                                                                 weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio());
02329                                                         }
02330                                                         else if(hypotheses[planeCtr] == 2){ 
02331                                                                 y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
02332                                                                 weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio());
02333                                                         }
02334                                                         else if(hypotheses[planeCtr] == 3){ 
02335                                                                 y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
02336                                                                 weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());
02337                                                         }
02338                                                 }
02339                                                 
02340                                                 ++arrayCtr;
02341                                                 //else MSG("AlgDeMuxBeam", Msg::kDebug) << "\tdrop plane " << plane->GetPlaneNumber() 
02342                                                 //<< endl; 
02343                                         } //end if plane is not the dropped one
02344                                 }//end if plane is in window bounds
02345                                 ++planeCtr;
02346                         } //end loop over planes to find variables for fit
02347             
02348                         //use the DmxUtilities function to find a linear fit
02349                         
02350                         util.FindLinearFit(x, y, weight, arrayCtr, fitA1, fitA2, chiSq);
02351                         MSG("AlgDeMuxBeam", Msg::kDebug) << "\tlinear fit drop plane " << chiSq << "\t" << fitA1 
02352                                 << "\t" << fitA2 << endl;
02353                         
02354                         planeItr.Reset();
02355                         
02356                         planeCtr = 0;
02357                         if(chiSq < prevChiSq){
02358                                 prevChiSq = chiSq;
02359                                 a1 = fitA1;
02360                                 a2 = fitA2;
02361                                 a3 = fitA3;
02362                                 a4 = fitA4;
02363                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "\trefined fit, drop plane " << i 
02364                                         << "\t" << prevChiSq << "\t" << a1 << "\t" << a2 << endl;
02365                         }
02366                         
02367                 }//end for loop to improve fit by dropping a plane
02368         }//end if enough planes to improve fit
02369         
02370         planeItr.ResetFirst();
02371         
02372         return;
02373 }

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 fMaxBadResidual, fMaxMuonRemovalFraction, fMaxResidual, Msg::kDebug, MSG, and RegionLinearFit().

Referenced by DeMuxFirstNPlanesTest(), UseMuonSlidingWindow(), and UseShowerSlidingWindow().

00900 {
00901         Float_t chiSqr;
00902         Float_t averageResidual = 0.;
00903         Float_t maxResidual = 0.;
00904         Float_t maxChiSqr = 0.;
00905         Float_t prevWeight = 0.;
00906         Int_t maxResidualLoc = 0;
00907         Int_t numBadResidual = 0;
00908         
00909         //now fit a straight line
00910         RegionLinearFit(slope, intercept, chiSqr, x, y, sigma, residual, numPoints);
00911         
00912         //find the fit for the first n planes and look for 
00913         //planes with large residuals.  repeat until planes with 
00914         //large residuals are removed
00915         Int_t nRemoved = 0;
00916         bool pointRemoved = true;
00917         while(nRemoved<=fMaxMuonRemovalFraction*numPoints && pointRemoved){
00918                 
00919                 pointRemoved = false;
00920                 RegionLinearFit(slope, intercept, chiSqr, x, y, sigma, residual, numPoints);
00921                 
00922                 averageResidual = 0.;
00923                 maxResidual = 0.;
00924                 numBadResidual = 0;
00925                 
00926                 for(Int_t i = 0; i<numPoints; ++i){
00927                         MSG("AlgDeMuxBeam", Msg::kDebug) << "zpos = " << x[i] << " tpos = " << y[i]
00928                                                                                         << " plane = " << planeNum[i]
00929                                                                                         << " residual = " << residual[i] << endl;
00930                         if(sigma[i]>0.){
00931                                 averageResidual += TMath::Abs(residual[i]);
00932                                 if(TMath::Abs(residual[i])>fMaxResidual) ++numBadResidual;
00933                         }
00934                         //find the point with the largest residual more than the allowed value
00935                         if(TMath::Abs(residual[i])>maxResidual){
00936                                 maxResidual = TMath::Abs(residual[i]);
00937                                 maxResidualLoc = i;
00938                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "biggest residual for plane " << planeNum[i] << endl;
00939                         }
00940                 }//end loop to find average residual and entry with the largest residual
00941                 
00942                 averageResidual /= (1.*(numPoints-nRemoved));
00943                 
00944                 MSG("AlgDeMuxBeam", Msg::kDebug) << "average residual for first n planes = " 
00945                         << averageResidual << "/" << fMaxResidual << endl;
00946                 
00947                 //is the average residual too big?  then remove the point with the 
00948                 //largest residual and try again.
00949                 if(averageResidual>fMaxResidual || numBadResidual>fMaxBadResidual){
00950                         ++nRemoved;
00951                         pointRemoved = true;
00952                         maxChiSqr = chiSqr;
00953                         //try removing one plane at a time from the fit and see which one causes the bad fit
00954                         for(Int_t j = 0; j<numPoints; ++j){
00955                                 prevWeight = sigma[j];
00956                                 sigma[j] = 0.;
00957                                 RegionLinearFit(slope, intercept, chiSqr, x, y, sigma, residual, numPoints);
00958                                 sigma[j] = prevWeight;
00959                                 if(chiSqr<maxChiSqr){
00960                                         maxChiSqr = chiSqr;
00961                                         maxResidualLoc = j;
00962                                 }
00963                         }//end loop over planes to get rid of the bad one
00964                         MSG("AlgDeMuxBeam", Msg::kDebug) << "!!!! removed plane " << planeNum[maxResidualLoc] 
00965                                 << " from fit !!!!" << endl;
00966                         sigma[maxResidualLoc] = 0.;
00967                 }
00968                 
00969         }//end loop to check fit
00970         
00971         return;
00972 }

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::kDebug, DmxPlaneTypes::kMuon, DmxPlaneTypes::kShower, MSG, and DmxPlane::SetGolden().

Referenced by FindFit().

01776 {
01777 
01778         MSG("AlgDeMuxBeam", Msg::kDebug) << "in FindPlanesToDropFromFit" << endl;
01779 
01780         planeItr.ResetFirst();
01781 
01782         Double_t diff1 = 0.;
01783         Double_t diff2 = 0.;
01784         Double_t diff3 = 0.;
01785         Double_t fitCoG = 0.;
01786         bool planesDropped = false;
01787 
01788         Int_t dropPlanes = 0;
01789         Int_t planeCtr = 0;
01790         //MSG("AlgDeMuxBeam", Msg::kDebug)<<"\tin FindPlanesToDropFromFit"; 
01791   
01792         DmxPlane *plane = 0;
01793         //find the plane farthest off from the fit, keep going until you have just 3 planes left
01794         while( (plane = planeItr()) && dropPlanes<(Int_t)(0.1*leverArm) && planeCtr<leverArm){
01795 
01796                 fitCoG = (a1 + (plane->GetZPosition())*a2 
01797                                   + TMath::Power(plane->GetZPosition(),2)*a3
01798                                   + TMath::Power(plane->GetZPosition(),3)*a4);
01799   
01800                 //MSG("AlgDeMuxBeam", Msg::kDebug)<<"\t" << plane->GetPlaneNumber() << endl;
01801       
01802                 if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){
01803                         if( plane->GetPlaneType() == DmxPlaneTypes::kShower ){
01804                                 diff1 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG());
01805                                 diff2 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG());
01806                                 diff3 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG());
01807                         }
01808                         else if( plane->GetPlaneType() == DmxPlaneTypes::kMuon ){
01809                                 diff1 = TMath::Abs(fitCoG - plane->GetCoG());
01810                                 diff2 = TMath::Abs(fitCoG - plane->GetCoG());
01811                                 diff3 = TMath::Abs(fitCoG - plane->GetCoG());
01812                         }
01813         
01814                         //if all the 3 best hypotheses are way off from the fit, that plane is most 
01815                         //likely messing up the fit.  so mark it as kTRUE - ie it, it truely sucks
01816                         if(diff1 > .504 && diff2 > .504 && diff3 > .504 && !plane->IsGolden() ){
01817                                 planeItr.GetSet()->GetMasks().SetMask(kTRUE, planeItr.Ptr());
01818                                 planesDropped = true;
01819                                 MSG("AlgDeMuxBeam", Msg::kDebug)<<"\tmarking plane " << plane->GetPlaneNumber() << endl;
01820                                 //decrement the number of planes now used for a fit.
01821                                 ++dropPlanes;
01822                         }
01823                         else if(diff1 > 0.504 && plane->IsGolden()){
01824           
01825                                 //if this was a golden plane but it just doesnt work, make a non-golden plane
01826                                 planeItr.GetSet()->GetMasks().SetMask(kTRUE, planeItr.Ptr());
01827                                 plane->SetGolden(false);
01828                                 planesDropped = true;
01829                                 MSG("AlgDeMuxBeam", Msg::kDebug)<<"\tmarking plane " << plane->GetPlaneNumber() << endl;
01830                                 //decrement the number of planes now used for a fit.
01831                                 ++dropPlanes;
01832                         }//end if dropping a golden plane
01833                 } //end if the mask hasnt already been set for this plane
01834       
01835                 ++planeCtr;
01836         }
01837         planeItr.ResetFirst();
01838   
01839         //cout << "finished drop planes from fit" << endl;
01840 
01841         return planesDropped;
01842 }

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

Definition at line 1541 of file AlgDeMuxBeam.cxx.

References DmxPlane::GetCoG(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), Msg::kDebug, DmxPlaneTypes::kMuon, DmxPlaneTypes::kShower, and MSG.

Referenced by RunAlg().

01542 {
01543         MSG("AlgDeMuxBeam", Msg::kDebug) << "In FindRegionMatedChargeFraction" << endl;
01544         
01545         //reset the validPlaneItr
01546         validPlaneItr.ResetFirst();
01547         DmxPlane *plane = 0;
01548         Float_t totalCharge = 0.;
01549         Float_t matedCharge = 0.;
01550         Float_t charge = 0.;
01551         
01552         while( (plane = validPlaneItr()) ){
01553                 
01554                 charge = plane->GetPlaneCharge();
01555                 totalCharge += charge;
01556                 
01557                 MSG("AlgDeMuxBeam", Msg::kDebug) << "plane = " << plane->GetPlaneNumber() 
01558                                                                                 << " charge = " << charge << endl;
01559 
01560                 if(plane->GetPlaneType() == DmxPlaneTypes::kShower)
01561                         matedCharge += dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio()*charge;
01562                 else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
01563                         
01564                         MSG("AlgDeMuxBeam", Msg::kDebug) << "muon plane cog = " << plane->GetCoG() 
01565                                                                                         << "/" << dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG()
01566                                                                                         << "/" << dynamic_cast<DmxMuonPlane *>(plane)->GetStripCoG()
01567                                                                                         << "/" << dynamic_cast<DmxMuonPlane *>(plane)->GetInitialStripCoG()
01568                                                                                         << endl;
01569                         
01570                         //since this is a valid muon plane, all the charge is mated, assuming that the initial
01571                         //center of gravity is the same as the set center of gravity
01572                         if(plane->GetCoG() == dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG())
01573                                 matedCharge += charge;
01574                         
01575                 }//end if muon plane
01576                         
01577                 MSG("AlgDeMuxBeam", Msg::kDebug) << " mated charge = " << matedCharge << endl;
01578                 
01579         }//end loop over valid planes
01580         
01581         
01582         //get the fraction of mated charge in this region.  check it against the 
01583         //passed in parameter
01584         Float_t matedChargeFrac = matedCharge;
01585         if(totalCharge > 0.) matedChargeFrac /= totalCharge;    
01586         else matedChargeFrac = -1.;
01587         
01588         MSG("AlgDeMuxBeam", Msg::kDebug) << "mated charge fraction in region = " 
01589                                                                         << matedChargeFrac << endl;
01590         
01591         return matedChargeFrac;
01592 }

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(), Msg::kDebug, DmxPlaneTypes::kShower, and MSG.

01849 {
01850         //   MSG("AlgDeMuxBeam", Msg::kDebug) << "In FindRMSOfPlanes" << endl;
01851         
01852         //initialize the variables for finding the RMS
01853         Double_t sigTPos2Sum = 0.;
01854         Double_t sigSum = 0.;
01855         Double_t sigTPosSum = 0.;
01856         Double_t signal = 0.;
01857         Double_t planeCoG = 0.;
01858         
01859         Int_t ctr = 0;
01860         
01861         planeItr.ResetFirst();
01862         //loop over the planes in the set and use only those that are 
01863         //between firstPlane and lastPlane
01864         
01865         DmxPlane *plane = 0;
01866         while( planeItr.IsValid() && ctr < leverArm){
01867                 plane = planeItr.Ptr();
01868                 signal = plane->GetPlaneCharge();
01869                 planeCoG = plane->GetCoG();
01870                 
01871                 if( plane->GetPlaneType() == DmxPlaneTypes::kShower ){
01872                         if( hypsToUse[ctr] == 1){
01873                                 signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio();
01874                                 planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
01875                         } 
01876                         else if( hypsToUse[ctr] == 2){
01877                                 signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
01878                                 planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
01879                         } 
01880                         else if( hypsToUse[ctr] == 3){
01881                                 signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
01882                                 planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
01883                         } 
01884                 }
01885                 
01886                 sigTPos2Sum += signal * planeCoG * planeCoG;
01887                 sigTPosSum += signal * planeCoG;
01888                 sigSum += signal;
01889                 
01890                 //increment the ctr
01891                 ++ctr;
01892                 
01893                 //move to the next plane in the list
01894                 planeItr.Next();
01895         }
01896         
01897         planeItr.Reset();
01898         
01899         //calculate the RMS using the formula from Bevington
01900         //[Sum(S_i * x_i^2) / Sum(S_i) - (av_x)^2] * (N / N -1) where
01901         //S_i is the signal at a distance x_i and av_x is the signal 
01902         //weighted average distance from the line for a digit and N is the
01903         //number of measures in the set
01904         
01905         Double_t weighter = (leverArm * 1.0) / (1.0 * (leverArm - 1));
01906         Double_t avTPos = sigTPosSum / sigSum;
01907         Double_t rmsSq = ((sigTPos2Sum / sigSum) - (avTPos * avTPos)) * weighter;
01908         
01909         cog = avTPos;
01910         
01911         //stupid round off error - make sure it doesnt make the rms^2 negative
01912         if(TMath::Abs(rmsSq)<1.e-3) rmsSq = 0.;
01913         
01914         rms = TMath::Sqrt(rmsSq);
01915         
01916         MSG("AlgDeMuxBeam", Msg::kDebug) << "cog = " << cog << " rms = " << rms << endl;
01917         
01918         
01919         //see what happens if you drop a plane from the set
01920         ctr = 0;
01921         
01922         //redo the fit using only the previously set planes and only if you have 
01923         //more than 3 planes in the set
01924         if(leverArm > 3){
01925                 
01926                 
01927                 Double_t rmsRef = 0.;
01928                 Double_t cogRef = 0.;
01929                 
01930                 for(Int_t i = 0; i < leverArm; i++){
01931                         ctr = 0;
01932                         sigTPosSum = 0;
01933                         sigTPos2Sum = 0;
01934                         sigSum = 0;
01935                         
01936                         //loop over the plane iterator to fill the variables
01937                         while( planeItr.IsValid() && ctr < leverArm){
01938                                 plane = planeItr.Ptr();
01939                                 
01940                                 if(i != ctr){
01941                                         
01942                                         signal = plane->GetPlaneCharge();
01943                                         planeCoG = plane->GetCoG();
01944                                         
01945                                         if( plane->GetPlaneType() == DmxPlaneTypes::kShower ){
01946                                                 if( hypsToUse[ctr] == 1){
01947                                                         signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio();
01948                                                         planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
01949                                                 } 
01950                                                 else if( hypsToUse[ctr] == 2){
01951                                                         signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
01952                                                         planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
01953                                                 } 
01954                                                 else if( hypsToUse[ctr] == 3){
01955                                                         signal *= dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
01956                                                         planeCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
01957                                                 } 
01958                                         }
01959                                         
01960                                         sigTPos2Sum += signal * planeCoG * planeCoG;
01961                                         sigTPosSum += signal * planeCoG;
01962                                         sigSum += signal;
01963                                         
01964                                         //        MSG("AlgDeMuxBeam", Msg::kDebug) << "inloop " << sigTPos2Sum << " " << sigTPosSum 
01965                                         //                                << " " << sigSum << endl;
01966                                         
01967                                 }//end if not the plane to drop
01968                                 ++ctr;
01969                                 
01970                                 planeItr.Next();
01971                         } //end loop over planes to find variables for fit
01972                         
01973                         planeItr.Reset();
01974                         
01975                         Double_t weighterRef = (leverArm * 1.0) / (1.0 * (leverArm - 1));
01976                         Double_t avTPosRef = sigTPosSum / sigSum;
01977                         Double_t rmsSqRef = ((sigTPos2Sum / sigSum) - (avTPosRef * avTPosRef)) * weighterRef;
01978                         
01979                         //     MSG("AlgDeMuxBeam", Msg::kDebug) << avTPosRef << " "  << sigTPos2Sum/sigSum << " "
01980                         //                          << weighterRef << endl;
01981                         
01982                         //stupid round off error - make sure it doesnt make the rms^2 negative
01983                         if(TMath::Abs(rmsSqRef)<1.e-3) rmsSqRef = 0.;
01984                         
01985                         cogRef = avTPosRef;
01986                         rmsRef = TMath::Sqrt(rmsSqRef);
01987                         
01988                         if(rmsRef < rms){  
01989                                 rms = rmsRef;
01990                                 cog = cogRef;
01991                                 
01992                                 //       MSG("AlgDeMuxBeam", Msg::kDebug) << "new cog = " << cog << " rms = " << rms << endl;
01993                                 
01994                         }
01995                         
01996                         
01997                 }//end loop over planes to drop one at a time
01998         }//end if at least 3 planes in set
01999         
02000         return;
02001 }

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, fFinalShowerRegionIntercept, fFinalShowerRegionSlope, fInitialMuonRegionIntercept, fInitialMuonRegionSlope, fMuonStartPlaneNumber, fMuonStartPlaneZPos, 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::kDebug, MSG, PlexSEIdAltL::Next(), PlexSEIdAltL::SetCurrentWeight(), and PlexSEIdAltL::SetFirst().

Referenced by RunAlg().

01602 {
01603         //get the last slope and intercept for the shower and muon regions and 
01604         //compare the predicted tpos at the muon start plane number.
01605         //also find the extent of the shower at that point.  if the muon track 
01606         //is more than the allowed value away from the shower, then reconcile 
01607         //the shower to the track.
01608         
01609         Float_t showerPredictedTPos = fFinalShowerRegionSlope*fMuonStartPlaneZPos;
01610         showerPredictedTPos += fFinalShowerRegionIntercept;
01611         
01612         Float_t muonPredictedTPos = fInitialMuonRegionSlope*fMuonStartPlaneZPos;
01613         muonPredictedTPos += fInitialMuonRegionIntercept;
01614         
01615         MSG("AlgDeMuxBeam", Msg::kDebug) << "shower prediction = " << showerPredictedTPos
01616                                                                         << " muon prediction = " << muonPredictedTPos 
01617                                                                         << " muon start zpos = " << fMuonStartPlaneZPos 
01618                                                                         << endl;
01619         
01620         CandDeMuxDigitHandle *digit = 0;
01621         DmxPlane *plane = 0;    
01622         
01623         UgliGeomHandle  ugh(fVldContext);
01624         UgliStripHandle ush;
01625 
01626         Float_t maxPlaneTPos = -5.;
01627         Float_t minPlaneTPos = 5.;
01628         Float_t digitTPos = 0.;
01629         Float_t setDigitTPos = 0.;
01630         
01631         if(TMath::Abs(muonPredictedTPos-showerPredictedTPos) > fShowerMuonTransverseDifference){
01632                 
01633                 //it looks like the muon and shower don't match up.  loop over the 
01634                 //CandDeMuxDigitList and see what the extent of the shower region is
01635                 while( (plane = planeItr()) ){
01636                         
01637                         maxPlaneTPos = -5.;
01638                         minPlaneTPos = 5.;
01639                         
01640                         
01641                         cdhit.GetSet()->Slice(plane->GetPlaneNumber());
01642                         
01643                         //got to tell the iterator to sort the slice
01644                         cdhit.GetSet()->Update();
01645                         while( (digit = cdhit()) ){
01646                                 if(digit->GetDeMuxDigitFlagWord() == 0){
01647                                         ush = ugh.GetStripHandle(digit->GetPlexSEIdAltL().GetBestSEId());
01648                                         digitTPos = ush.GetTPos();
01649                                         if(digitTPos > maxPlaneTPos) maxPlaneTPos = digitTPos;
01650                                         else if(digitTPos < minPlaneTPos) minPlaneTPos = digitTPos;
01651                                 }
01652                         }//end loop over digits in the plane
01653 
01654                         //reset the digit iterator
01655                         cdhit.ResetFirst();
01656                         
01657                         //check the demuxed cog of each plane in the shower region against 
01658                         //the predicted position based on the initial muon linear fit. also
01659                         //check that the predicted position is outside the minimum and maximum
01660                         //transverse positions of the reconstructed strips in the the plane.
01661                         //if the difference is more than the allowed value, or the predicted 
01662                         //location is outside the reconstructed strip bounds, then redo the reconstruction
01663                         muonPredictedTPos = fInitialMuonRegionSlope*plane->GetZPosition();
01664                         muonPredictedTPos += fInitialMuonRegionIntercept;
01665 
01666                         if(plane->GetPlaneNumber() < fMuonStartPlaneNumber
01667                            && TMath::Abs(plane->GetCoG()-muonPredictedTPos) > fShowerMuonTransverseDifference
01668                            && (muonPredictedTPos<minPlaneTPos || muonPredictedTPos>maxPlaneTPos)){
01669                                 
01670                                 //plane->SetStrips(muonPredictedTPos);
01671         
01672                                 //loop over the digits again and look at the SEIdAltLists - if there is 
01673                                 //an alternative within the accepted bound of the predicted position,
01674                                 //take it.
01675                                 while( (digit = cdhit()) ){ 
01676                                         
01677                                         ush = ugh.GetStripHandle(digit->GetPlexSEIdAltL().GetBestSEId());
01678                                         setDigitTPos = ush.GetTPos();
01679 
01680                                         //set the alternatives list to the first one
01681                                         digit->GetPlexSEIdAltL().SetFirst();
01682                                         
01683                                         while( digit->GetPlexSEIdAltL().IsValid() ){
01684                                                 ush = ugh.GetStripHandle(digit->GetPlexSEIdAltL().GetCurrentSEId());
01685                                                 digitTPos = ush.GetTPos();
01686                                                         
01687                                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "current tpos = " << digitTPos
01688                                                                                                                 << "/" << muonPredictedTPos << endl;
01689                                                 
01690                                                 if(TMath::Abs(digitTPos-muonPredictedTPos) < fDigitResetCut
01691                                                    || (TMath::Abs(setDigitTPos-muonPredictedTPos) > fShowerMuonTransverseDifference
01692                                                            && TMath::Abs(digitTPos-muonPredictedTPos) < fShowerMuonTransverseDifference) ){
01693                                                         //clear all weights for the digit
01694                                                         digit->GetPlexSEIdAltLWritable().ClearWeights();
01695                                                         MSG("AlgDeMuxBeam", Msg::kDebug) << "setting" << endl;
01696                                                         digit->GetPlexSEIdAltLWritable().SetCurrentWeight(1.0);
01697                                                         break;
01698                                                 }
01699                                                         
01700                                                 digit->GetPlexSEIdAltL().Next();
01701                                         }//end loop over alt list
01702                                 }//end loop over digits to find the ones near the predicted tpos
01703                         }//end if this plane should have some digits reset
01704                         
01705                         //clear the digit slice
01706                         cdhit.GetSet()->ClearSlice(false);
01707                                 
01708                 }//end loop over planeItr
01709                 
01710         }//end if the shower and muon regions have different predicted transverse positions
01711         
01712         return;
01713 }

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::kDebug, Msg::kWarning, and MSG.

Referenced by FindFitOverNPlanes().

01480 {
01481 
01482         //define the necessary sums
01483         Double_t sumXSqrDivSigmaSqr = 0.;
01484         Double_t sumYDivSigmaSqr = 0.;
01485         Double_t sumXDivSigmaSqr = 0.;
01486         Double_t sumXYDivSigmaSqr = 0.;
01487         Double_t sum1DivSigmaSqr = 0.;
01488         Double_t sigmaSqr = 0.;
01489         chiSqr = 0.;
01490         
01491         for(Int_t i=0; i<numPoints; ++i){
01492 
01493                 if(TMath::Abs(sigma[i])>0.){
01494                         sigmaSqr = sigma[i]*sigma[i];
01495                         sum1DivSigmaSqr += 1./sigmaSqr;
01496                         sumXDivSigmaSqr += x[i]/sigmaSqr;
01497                         sumYDivSigmaSqr += y[i]/sigmaSqr;
01498                         sumXYDivSigmaSqr += (x[i]*y[i])/sigmaSqr;
01499                         sumXSqrDivSigmaSqr += (x[i]*x[i])/sigmaSqr;
01500                 }
01501                 MSG("AlgDeMuxBeam", Msg::kDebug) << "x = " << x[i] << " y = " << y[i] << " sigma = " 
01502                                                                     << sigma[i] << " sum1DivSigmaSqr = " << sum1DivSigmaSqr
01503                                                                     << " sumXDivSigmaSqr = " << sumXDivSigmaSqr
01504                                                                     << " sumYDivSigmaSqr = " << sumYDivSigmaSqr
01505                                                                     << " sumXYDivSigmaSqr = " << sumXYDivSigmaSqr
01506                                                                     << " sumXSqrDivSigmaSqr = " << sumXSqrDivSigmaSqr << endl;
01507         }
01508         
01509         Float_t delta = sum1DivSigmaSqr*sumXSqrDivSigmaSqr - sumXDivSigmaSqr*sumXDivSigmaSqr;
01510         if(delta==0.){
01511                 MSG("AlgDeMuxBeam", Msg::kWarning) << "delta = 0 in linear fit - kludge an answer and bail" << endl;
01512                 intercept = 0.;
01513                 slope = 0.;
01514                 chiSqr = 10000.;
01515                 return;
01516         }
01517         intercept = (sumXSqrDivSigmaSqr*sumYDivSigmaSqr - sumXDivSigmaSqr*sumXYDivSigmaSqr)/delta;
01518         slope = (sum1DivSigmaSqr*sumXYDivSigmaSqr - sumXDivSigmaSqr*sumYDivSigmaSqr)/delta;
01519         
01520         //now get the chi^2 and residual for the fit - zero chiSqr to be safe
01521         chiSqr = 0.;
01522         Double_t numerator = 0.;
01523         Double_t denominator = 0.;
01524         for(Int_t i=0; i<numPoints; ++i){
01525                 residuals[i] = (y[i]-intercept-slope*x[i]);
01526                 numerator = TMath::Power(y[i]-intercept-slope*x[i],2.);
01527                 denominator = sigma[i]*sigma[i];
01528                 if(denominator>0.) chiSqr += numerator/denominator;
01529         }
01530         
01531         //fit to a straight line has 2 degrees of freedom
01532         if(numPoints>2) chiSqr /= 1.*(numPoints-2);
01533         
01534         return;
01535 }

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

Implements AlgBase.

Definition at line 91 of file AlgDeMuxBeam.cxx.

References bfld::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(), 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(), Msg::kDebug, KeyFromPlane1(), KeyOnPlane1(), KeyOnUView1(), KeyOnValidU1(), KeyOnValidV1(), KeyOnVView1(), Detector::kFar, DmxPlaneTypes::kMuon, CandDeMuxDigitList::kNonPhysicalStripSolution, CandDeMuxDigitList::kNoVertex, MSG, ReconcileShowerAndMuonRegions(), DmxStatus::ResetStatus(), DmxStatus::SetAverageTimingOffset(), DmxStatus::SetEndPlaneNumber(), DmxStatus::SetEventNumber(), DmxStatus::SetFigureOfMerit(), DmxStatus::SetMultipleMuon(), DmxStatus::SetMuonStartPlaneNumber(), DmxStatus::SetNoVertexFailure(), DmxStatus::SetNumberOfPlanes(), 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().

00092 {
00093 
00094   assert( ch.InheritsFrom("CandDigitListHandle") );
00095   CandDeMuxDigitListHandle &cdlh = dynamic_cast<CandDeMuxDigitListHandle&>(ch);
00096   CandDeMuxDigitHandleItr cdhit = cdlh.GetDaughterIterator();
00097 
00098   fVldContext = *(ch.GetVldContext());
00099 
00100   if (fVldContext.GetDetector() != Detector::kFar) {
00101     static int msglimit = 20; // no more than 20 warnings about the detector
00102     if (msglimit) {
00103       MSG("AlgDeMuxBeam",Msg::kDebug) 
00104         << "AlgDeMuxBeam::RunAlg() can not DeMux "
00105         << Detector::AsString(fVldContext.GetDetector())
00106         << " detector.  Skip futher DeMux processing."
00107         << endl;
00108       if (--msglimit == 0) 
00109         MSG("AlgDeMuxBeam",Msg::kDebug) 
00110           << " ... last message of this type." << endl;
00111     }
00112     return;
00113   }
00114 
00115   //get the AlgConfigDeMux object
00116   fPlanesInSet = acd.GetInt("PlanesInSet");
00117   fMinMuonPlanesRequired = acd.GetInt("MinMuonPlanesRequired");
00118   fMinMatedChargeFrac = acd.GetDouble("MinMatedChargeFraction");
00119   fMuonPlanesInSet = acd.GetInt("MuonPlanesInSet");
00120   fMuonSlopeLimit = acd.GetDouble("MuonTrackSlopeLimit");
00121   fMaxResidual = acd.GetDouble("MaxResidual");
00122   fMaxBadResidual = acd.GetInt("MaxBadResidual");
00123   fMaxMuonRemovalFraction = acd.GetDouble("MaxMuonRemovalFraction");
00124   fMaxStripAdjustment = acd.GetInt("MaxStripAdjustment");
00125   fStrayCut = acd.GetInt("StrayDeltaStripLimit");
00126   fNumberOfHypotheses = acd.GetInt("NumberOfHypotheses");
00127   fShowerMuonTransverseDifference = acd.GetDouble("ShowerMuonTransverseDifference");
00128   fDigitResetCut = acd.GetDouble("ResetDigitCut");
00129   
00130   //MSG("AlgDeMuxBeam", Msg::kDebug) << "window size = " << fPlanesInSet << endl;
00131 
00132   //get the DmxStatus object
00133   // Find the DmxStatus object or create one for needed scratch space
00134   DmxStatus *status = dynamic_cast<DmxStatus *>(gROOT->GetRootFolder()->FindObject("Loon/DeMux/DmxStatus"));
00135   bool tempStatus = false;
00136   if (status == 0) {
00137     MSG("AlgDeMuxBeam", Msg::kDebug) << "//root/Loon/DeMux/DmxStatus not found."
00138                               << " Create a temporary DmxStatus." << endl;
00139     status = new DmxStatus;      // Make a temporary DmxStatus if needed
00140     tempStatus = true;
00141   }
00142   else status->ResetStatus(); 
00143 
00144   status->SetEventNumber(ch.GetCandRecord()->GetCandHeader()->GetSnarl());
00145  
00146   MSG("AlgDeMuxBeam", Msg::kDebug) << "starting beam demuxing for event " 
00147                                                                   << status->GetEventNumber()
00148                                                                   << endl;
00149 
00150   DmxUtilities util;
00151   
00152   util.FillPlaneArray(status, cdlh, acd);
00153   const TObjArray *planeArray = status->GetPlaneArray();
00154   //DmxPlane *plane = dynamic_cast<DmxPlane *>(planeArray->First());
00155   //MSG("AlgDeMuxBeam", Msg::kDebug) <<"Event = " <<  status->GetEventNumber() << endl;
00156   
00157   //create the DmxPlaneItr over planes and program it to sort the planes
00158   DmxPlaneItr planeItr(planeArray);
00159   //create a KeyFunc to sort planes by number
00160   DmxPlaneKeyFunc *pnKF = planeItr.CreateKeyFunc();
00161   pnKF->SetFun(KeyOnPlane1);
00162   planeItr.GetSet()->AdoptSortKeyFunc(pnKF);
00163   pnKF = 0;
00164   
00165   planeItr.ResetFirst();
00166   status->SetNumberOfPlanes(planeItr.SizeSelect());
00167   status->SetVertexPlaneNumber(util.FindVertexPlane(planeItr));
00168   
00169   //demux the event if there is an identified vertex
00170   if( status->GetVertexPlaneNumber() != -1){
00171     
00172     planeItr.ResetFirst();
00173     planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), 500);
00174     planeItr.GetSet()->Update();
00175 
00176     //set the end plane number since there was a vertex plane
00177     status->SetEndPlaneNumber(util.FindEndPlane(planeItr));
00178     planeItr.GetSet()->ClearSlice(false);
00179     planeItr.ResetFirst();
00180 
00181     //set the z position of the vertex plane
00182     planeItr.GetSet()->Slice(status->GetVertexPlaneNumber());
00183     planeItr.GetSet()->Update();
00184     planeItr.ResetFirst();
00185     status->SetVertexPlaneZPosition(dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetZPosition());
00186 
00187     planeItr.GetSet()->ClearSlice(false);
00188     planeItr.ResetFirst();
00189 
00190     //set the plane number for the muon start plane
00191         planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber());
00192     planeItr.GetSet()->Update();
00193     //set the muon start plane number
00194     status->SetMuonStartPlaneNumber(util.FindBeamMuonStartPlane(planeItr, 
00195                                                                                                                                 acd.GetDouble("AveragePEGainConversion")));
00196         fMuonStartPlaneNumber = status->GetMuonStartPlaneNumber();
00197     planeItr.GetSet()->ClearSlice(false);
00198     planeItr.ResetFirst();
00199     
00200     //create a KeyFunc to sort digits by plane
00201     CandDeMuxDigitHandleKeyFunc *planeKF = cdhit.CreateKeyFunc();
00202     
00203     //program the KeyFunc with the sort function
00204     planeKF->SetFun(KeyFromPlane1);
00205     
00206     //get the NavSet from the iterator and pass the KeyFunc to it
00207     cdhit.GetSet()->AdoptSortKeyFunc(planeKF);
00208     
00209     //clear the KF pointer because we no longer own the KeyFunc
00210     planeKF = 0;
00211 
00212     //make the iterators for the window, valid planes, and vertex
00213     DmxPlaneItr vertexPlaneItr(planeArray);
00214     DmxPlaneItr validPlaneItr(planeArray);
00215        
00216     //create a KeyFunc to sort planes by number
00217     DmxPlaneKeyFunc *vertexPlaneKF = vertexPlaneItr.CreateKeyFunc();
00218     DmxPlaneKeyFunc *validPlaneKF = validPlaneItr.CreateKeyFunc();
00219        
00220     //program the KeyFunc with the sort function
00221     vertexPlaneKF->SetFun(KeyOnPlane1);
00222     validPlaneKF->SetFun(KeyOnPlane1);
00223        
00224     //get the NavSet from the iterator and pass the KeyFunc to it
00225     vertexPlaneItr.GetSet()->AdoptSortKeyFunc(vertexPlaneKF);
00226     validPlaneItr.GetSet()->AdoptSortKeyFunc(validPlaneKF);
00227        
00228     //clear the KF pointer because we no longer own the KeyFunc
00229     vertexPlaneKF = 0;
00230     validPlaneKF = 0;
00231   
00232     planeItr.ResetFirst();
00233     validPlaneItr.ResetFirst();
00234         
00235     //now program a key function to select on plane orientation - U first
00236     DmxPlaneKeyFunc *orientValidUKF = validPlaneItr.CreateKeyFunc();
00237     DmxPlaneKeyFunc *orientUKF = planeItr.CreateKeyFunc();
00238     DmxPlaneKeyFunc *orientValidVKF = validPlaneItr.CreateKeyFunc();
00239     DmxPlaneKeyFunc *orientVKF = planeItr.CreateKeyFunc();
00240         
00241         //initialize the muon and shower region mated fraction variables in the DmxStatus object
00242         status->SetUShowerMatedFraction(-1.);   
00243         status->SetUMuonMatedFraction(-1.);   
00244         status->SetVShowerMatedFraction(-1.);   
00245         status->SetVMuonMatedFraction(-1.);   
00246         
00247         DmxPlane *plane = 0;
00248         bool endPlaneValid = false;
00249 
00250         planeItr.GetSet()->Update();
00251         planeItr.ResetFirst();
00252         //see if it worked
00253         while( (plane = planeItr()) ){
00254                 MSG("AlgDeMuxBeam", Msg::kDebug) << "number = " << plane->GetPlaneNumber()
00255                                                                                 << " orientation = " << (Int_t)plane->GetPlaneView()
00256                                                                                 << " type = " 
00257                                                                                 << DmxPlaneTypes::GetPlaneTypeString(plane->GetPlaneType())
00258                                                                                 << " IsValid() = " << (Int_t)plane->IsValid()  << endl;
00259                 
00260                 if(plane->GetPlaneNumber() == fMuonStartPlaneNumber) 
00261                         fMuonStartPlaneZPos = plane->GetZPosition();
00262                 if(plane->GetPlaneNumber() == status->GetEndPlaneNumber() && plane->IsValid())
00263                         endPlaneValid = true;
00264         }
00265         planeItr.ResetFirst();
00266         
00267     Int_t viewCtr = 0;
00268         Int_t numValidMuonPlanes = 0;
00269         Float_t muonMatedFraction = -1.;
00270         Float_t showerMatedFraction = -1.;
00271     while( viewCtr < 2 ){
00272                 
00273                 //set the muon start plane number at the start of the loop for each 
00274                 //view as one view may have reset it to -1 if there were not enough valid 
00275                 //muon planes in that view after the muon start plane number
00276                 fMuonStartPlaneNumber = status->GetMuonStartPlaneNumber();
00277                 muonMatedFraction = -1.;
00278                 showerMatedFraction = -1.;
00279                 
00280                 MSG("AlgDeMuxBeam", Msg::kDebug) << "vertex plane = " << status->GetVertexPlaneNumber()
00281                                                                                 << " muon start plane = " << fMuonStartPlaneNumber
00282                                                                                 << " end plane = " << status->GetEndPlaneNumber() << endl;
00283                 
00284      
00285                 if( viewCtr == 0){
00286                         //program this function with the select function
00287                         orientUKF->SetFun(KeyOnUView1);
00288                         orientValidUKF->SetFun(KeyOnValidU1);
00289                 
00290                         //adopt it as a selection function
00291                         planeItr.GetSet()->AdoptSelectKeyFunc(orientUKF);
00292                         orientUKF = 0;
00293                         planeItr.Reset();
00294                         validPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidUKF);
00295                         orientValidUKF = 0;
00296                         validPlaneItr.Reset();
00297                 }
00298                 else if(viewCtr == 1){
00299                         //program this function with the select function
00300                         orientVKF->SetFun(KeyOnVView1);
00301                         orientValidVKF->SetFun(KeyOnValidV1);
00302                 
00303                         //adopt it as a selection function
00304                         planeItr.GetSet()->AdoptSelectKeyFunc(orientVKF);
00305                         orientVKF = 0;
00306                         planeItr.Reset();
00307                         validPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidVKF);
00308                         orientValidVKF = 0;
00309                         validPlaneItr.Reset();
00310                 }
00311                   
00312                 planeItr.ResetFirst();
00313       
00314                 //there must be at least 2 valid planes for the demuxing to work.
00315                 //slice the validPlaneItr to all planes in the event
00316                 validPlaneItr.GetSet()->Slice(status->GetVertexPlaneNumber(), 
00317                                                                           status->GetEndPlaneNumber());
00318                 validPlaneItr.GetSet()->Update();
00319                 
00320                 if(validPlaneItr.SizeSelect() >= 2){
00321         
00322                         planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), 
00323                                                                          status->GetEndPlaneNumber()+1);
00324                         planeItr.GetSet()->Update();
00325 
00326                         //demux the first n planes in the detector.  the method adjusts the size of n
00327                         //such that if there are < fPlanesInSet planes hit, n becomes the number of planes
00328                         //hit.  Otherwise, n is fPlanesInSet 
00329                         if((Int_t)validPlaneItr.SizeSelect() < fPlanesInSet)
00330                                 DeMuxFirstNPlanesTest(planeItr, validPlaneItr);
00331         
00332                         else if( (Int_t)validPlaneItr.SizeSelect() >= fPlanesInSet){
00333 
00334                                 //still have to demux the first n planes of the event the same for 
00335                                 //muon or shower like events.  loop over the planes to find the last
00336                                 //plane in the first set of n
00337                                 Int_t ctr = 0;
00338                                 Int_t lastPlane = 0;
00339                                 Int_t firstPlane = 0;
00340                                 
00341                                 while( validPlaneItr.IsValid() && ctr < fPlanesInSet){
00342             
00343                                         plane = validPlaneItr.Ptr();
00344                                         //MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber()
00345                                         //             << endl;
00346 
00347                                         if( ctr == 0 ) firstPlane = plane->GetPlaneNumber();
00348                                         lastPlane = plane->GetPlaneNumber(); 
00349             
00350                                         validPlaneItr.Next();
00351                                         ++ctr;
00352                                 }
00353 
00354                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "\tfirst plane = " << firstPlane 
00355                                                                                                  << "\tlast plane = " << lastPlane << endl;
00356 
00357                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "demuxing first N planes" << endl;
00358 
00359                                 //clear the slice of the validPlaneItr
00360                                 validPlaneItr.GetSet()->ClearSlice(false);
00361                                 planeItr.GetSet()->ClearSlice(false);
00362 
00363                                 //set the initial window 
00364                                 validPlaneItr.GetSet()->Slice(status->GetVertexPlaneNumber(), lastPlane+1);
00365                                 validPlaneItr.GetSet()->Update();
00366                                 
00367                                 //check to see if the number of valid planes == fPlanesInSet.  if so, then 
00368                                 //it means that there wont be any sliding windows, so set all planes in 
00369                                 //the event
00370                                 if(status->GetEndPlaneNumber()-lastPlane <= 2 && !endPlaneValid) 
00371                                         lastPlane = status->GetEndPlaneNumber();
00372                                 
00373                                 //set the plane iterator to be from the vertex to the last plane in the window,
00374                                 //since this is the first window done
00375                                 planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), lastPlane+1);
00376                                 planeItr.GetSet()->Update();
00377 
00378                                 DeMuxFirstNPlanesTest(planeItr, validPlaneItr);
00379                                 
00380                                 validPlaneItr.GetSet()->ClearSlice(false);
00381                                 planeItr.GetSet()->ClearSlice(false);
00382 
00383                                 //check to make sure that there are valid planes after the 
00384                                 //muon start plane
00385                                 validPlaneItr.GetSet()->Slice(firstPlane, status->GetEndPlaneNumber());
00386                                 validPlaneItr.GetSet()->Update();
00387                                 numValidMuonPlanes = 0;
00388                                 while( (plane = validPlaneItr() ) ){
00389                                         if(plane->GetPlaneNumber()>fMuonStartPlaneNumber 
00390                                            && fMuonStartPlaneNumber > 0 
00391                                            && plane->GetPlaneType()==DmxPlaneTypes::kMuon)
00392                                                 ++numValidMuonPlanes;
00393                                 }
00394                                 validPlaneItr.GetSet()->ClearSlice(false);
00395                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "number valid muon planes in view " << viewCtr
00396                                                                                                  << " = " << numValidMuonPlanes << endl;
00397                                 
00398                                 //if the number of valid muon planes is less than 2, then set the 
00399                                 //fMuonStartPlaneNumber to -1 since you will just demux the muon stuff 
00400                                 //with the shower region
00401                                 if(numValidMuonPlanes<fMinMuonPlanesRequired) fMuonStartPlaneNumber = -1;
00402                                 
00403                                 //demux the shower region if it exists and it extends beyond the 
00404                                 //first set of planes demuxed
00405                                 if(status->GetVertexPlaneNumber() != fMuonStartPlaneNumber
00406                                    && fMuonStartPlaneNumber > 0 && TMath::Abs(fMuonStartPlaneNumber-lastPlane)>2){
00407                                         //slice the valid plane iterator to do the shower region
00408                                         MSG("AlgDeMuxBeam", Msg::kDebug) << "slice shower region with track from " 
00409                                                                                                         << firstPlane+1 << " to " 
00410                                                                                                         << fMuonStartPlaneNumber << endl;
00411                                         validPlaneItr.GetSet()->Slice(firstPlane+1, fMuonStartPlaneNumber);
00412                                         validPlaneItr.GetSet()->Update();
00413                                         UseShowerSlidingWindow(planeItr, validPlaneItr, status);
00414                                         showerMatedFraction = FindRegionMatedChargeFraction(validPlaneItr);
00415                                         if(viewCtr == 0)
00416                                                 status->SetUShowerMatedFraction(showerMatedFraction);   
00417                                         if(viewCtr == 1)
00418                                                 status->SetVShowerMatedFraction(showerMatedFraction);   
00419                                                                                                         
00420                                         validPlaneItr.GetSet()->ClearSlice(false);
00421                                                                                 
00422                                         MSG("AlgDeMuxBeam", Msg::kDebug) <<"Done Setting Shower Region for snarl " 
00423                                                                                                         << status->GetEventNumber() << endl;
00424                                 }
00425                                 else if(fMuonStartPlaneNumber < 0 
00426                                                 && TMath::Abs(status->GetEndPlaneNumber()-lastPlane)>2){
00427                                         //slice the valid plane iterator to do the shower region - there is 
00428                                         //no track to speak of, so just slice the planes out to the end of 
00429                                         //the detector
00430                                         validPlaneItr.GetSet()->Slice(firstPlane+1, 486);
00431                                         validPlaneItr.GetSet()->Update();
00432                                         UseShowerSlidingWindow(planeItr, validPlaneItr, status);
00433                                         showerMatedFraction = FindRegionMatedChargeFraction(validPlaneItr);
00434                                         if(viewCtr == 0)
00435                                                 status->SetUShowerMatedFraction(showerMatedFraction);   
00436                                         if(viewCtr == 1)
00437                                                 status->SetVShowerMatedFraction(showerMatedFraction);   
00438                                         
00439                                         validPlaneItr.GetSet()->ClearSlice(false);
00440             
00441                                         MSG("AlgDeMuxBeam", Msg::kDebug) <<"Done Setting Shower Region for no track or small track " 
00442                                                                                                         << status->GetEventNumber() << endl;
00443                                 }
00444                                 //now do the muon region if it exists           
00445                                 //if( fMuonStartPlaneNumber == status->GetVertexPlaneNumber() 
00446                                 //      && numValidMuonPlanes>1){
00447                                 //      validPlaneItr.GetSet()->Slice(firstPlane+1, 
00448                                 //                                                                status->GetEndPlaneNumber());
00449                                 //      validPlaneItr.GetSet()->Update();
00450                                 //      UseMuonSlidingWindow(planeItr,validPlaneItr, status);
00451                                 //      validPlaneItr.GetSet()->ClearSlice(false);
00452                                         
00453                                 //      MSG("AlgDeMuxBeam", Msg::kDebug) << "done setting the muon for snarl " 
00454                                 //              << status->GetEventNumber()  << endl;
00455                                 //}
00456                                 if(fMuonStartPlaneNumber > 0  && numValidMuonPlanes>1){
00457                                         validPlaneItr.GetSet()->Slice(fMuonStartPlaneNumber, 
00458                                                                                                   status->GetEndPlaneNumber());
00459                                         validPlaneItr.GetSet()->Update();
00460                                         UseMuonSlidingWindow(planeItr,validPlaneItr, status);
00461                                         muonMatedFraction = FindRegionMatedChargeFraction(validPlaneItr);
00462                                         if(viewCtr == 0)
00463                                                 status->SetUMuonMatedFraction(muonMatedFraction);   
00464                                         if(viewCtr == 1)
00465                                                 status->SetVMuonMatedFraction(muonMatedFraction);   
00466                                         
00467                                         validPlaneItr.GetSet()->ClearSlice(false);
00468                                         
00469                                         MSG("AlgDeMuxBeam", Msg::kDebug) << "done setting the muon for snarl " 
00470                                                                                                         << status->GetEventNumber()  << endl;
00471                                 }
00472                                 
00473                                 //check and see that both tracklike and showerlike regions,
00474                                 //that the fraction of the charge in the muon region that is mated
00475                                 //is above the threshold and that there are a reasonable number of 
00476                                 //planes in the muon region for the fit
00477                                 //if so, see that the track and shower solutions are consistent
00478                                 planeItr.GetSet()->Update();
00479                                 planeItr.ResetFirst();
00480                                 if(status->GetVertexPlaneNumber() != fMuonStartPlaneNumber
00481                                    && muonMatedFraction > fMinMatedChargeFrac
00482                                    && numValidMuonPlanes >= fPlanesInSet)
00483                                         ReconcileShowerAndMuonRegions(planeItr, cdhit);
00484                                 
00485                                 //MSG("AlgDeMuxBeam", Msg::kDebug)<<"Done with enough planes to demux"<< endl;
00486                         }//end if at least fPlanesInSet valid planes to reconstruct
00487                 } //end if >=2 valid planes to reconstruct the event
00488                 else{
00489                         MSG("AlgDeMuxBeam", Msg::kDebug)<< "Event " << status->GetEventNumber() 
00490                                                                            << " not demuxed; not enough valid planes" << endl;
00491                         status->SetValidPlanesFailure(true); 
00492                 }
00493 
00494                 //what if there is more to the track?  perhaps the muon went through the
00495                 //coil hole and disappeared for several planes?  treat it in a way similar 
00496                 //to multiple cosmics
00497 
00498                 planeItr.GetSet()->ClearSlice(false);
00499                 validPlaneItr.GetSet()->ClearSlice(false);
00500                 vertexPlaneItr.GetSet()->ClearSlice(false);
00501                 //MSG("AlgDeMuxBeam", Msg::kDebug) << "slice vertex itr" << endl;
00502       
00503                 vertexPlaneItr.GetSet()->Slice(status->GetEndPlaneNumber()+1, 500);
00504                 vertexPlaneItr.GetSet()->Update();
00505                 //MSG("AlgDeMuxBeam", Msg::kDebug)<< "event =  " << status->GetEventNumber() 
00506                 //                   << "\tend plane = " << status->GetEndPlaneNumber() << endl;
00507       
00508 
00509                 //MSG("AlgDeMuxBeam", Msg::kDebug) << "check for new vertex" << endl;
00510 
00511                 Int_t nextVertex = -1;
00512                 if(vertexPlaneItr.SizeSelect() > 3)
00513                         nextVertex = util.FindVertexPlane(vertexPlaneItr);
00514                         vertexPlaneItr.GetSet()->ClearSlice(false);
00515       
00516                         while( nextVertex != -1){
00517         
00518                                 status->SetMultipleMuon(true); 
00519                                 vertexPlaneItr.GetSet()->ClearSlice(false);
00520                                 vertexPlaneItr.GetSet()->Slice(nextVertex, 500);
00521                                 vertexPlaneItr.GetSet()->Update();
00522 
00523                                 //get the next end plane
00524                                 Int_t endPlane = util.FindEndPlane(vertexPlaneItr);
00525                                 
00526                                 status->SetEndPlaneNumber(endPlane);
00527 
00528                                 MSG("AlgDeMuxBeam", Msg::kDebug)<< "\tin gap loop, vertex at plane " << nextVertex 
00529                                                                                    << "\tend plane = " << endPlane << "\tplanes in slice = " 
00530                                                                                    << vertexPlaneItr.SizeSelect() << endl;
00531         
00532                                 //if we are in this loop, then it is most likely because a muon went through
00533                                 //the coil hole.  so the thing to do is just treat all the remaining planes as 
00534                                 //being in the muon track and only demux that way
00535 
00536                                 validPlaneItr.GetSet()->Slice(nextVertex, endPlane);
00537                                 validPlaneItr.GetSet()->Update();
00538                                 validPlaneItr.ResetFirst();
00539                                 
00540                                 //loop over the valid planes and count the number of muon planes
00541                                 numValidMuonPlanes = 0;
00542                                 while( (plane = validPlaneItr()) ){
00543                                         MSG("AlgDeMuxBeam", Msg::kDebug) << "number = " << plane->GetPlaneNumber()
00544                                         << " orientation = " << (Int_t)plane->GetPlaneView()
00545                                         << " type = " 
00546                                         << DmxPlaneTypes::GetPlaneTypeString(plane->GetPlaneType())
00547                                         << " IsValid() = " << (Int_t)plane->IsValid()  << endl;
00548                                         
00549                                         if(plane->GetPlaneType()==DmxPlaneTypes::kMuon) ++numValidMuonPlanes;
00550                                 }
00551                                 
00552                                 validPlaneItr.ResetFirst();
00553                                 
00554                                 if(numValidMuonPlanes<2){
00555                                         nextVertex = -1;
00556                                         continue;
00557                                 }
00558                                 UseMuonSlidingWindow(planeItr,validPlaneItr, status);
00559                                 validPlaneItr.GetSet()->ClearSlice(false);
00560                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "done setting the muon for snarl " 
00561                                                                                     << status->GetEventNumber()  << endl;
00562         
00563                                 //look for another vertex after the end plane
00564                                 planeItr.GetSet()->ClearSlice(false);
00565                                 validPlaneItr.GetSet()->ClearSlice(false);
00566                                 vertexPlaneItr.GetSet()->ClearSlice(false);
00567                                 vertexPlaneItr.GetSet()->Slice(endPlane+1, 500);
00568                                 vertexPlaneItr.GetSet()->Update();
00569                                 //MSG("AlgDeMuxBeam", Msg::kDebug)<< "\tlook for new vertex, slice from " << endPlane+1 
00570                                 //               << "\tto plane 500, planes in slice = " 
00571                                 //               << vertexPlaneItr.SizeSelect() << endl;
00572                                 
00573                                 vertexPlaneItr.Reset();
00574         
00575                                 if(vertexPlaneItr.SizeSelect() > 0) 
00576                                         nextVertex = util.FindVertexPlane(vertexPlaneItr);
00577                                 else nextVertex = -1; 
00578                         }//end loop to find next vertex for spread out multiples
00579 
00580                         fStrayPlanes = util.FindPlanesOffFit(planeItr, fStrayCut);
00581       
00582                         MSG("AlgDeMuxBeam", Msg::kDebug) << status->GetEventNumber() << " " 
00583                                                                                 << fStrayPlanes << endl; 
00584       
00585                         status->SetFigureOfMerit(validPlaneItr.SizeSelect(), fStrayPlanes);
00586                         if(viewCtr == 0 && status->GetEventDeMuxed() ){ 
00587                                 status->SetUStrayPlanes(fStrayPlanes);
00588                                 status->SetUValidPlanes(validPlaneItr.SizeSelect());
00589                                 cdlh.SetNumValidPlanesU(validPlaneItr.SizeSelect());
00590                                 cdlh.SetNumStrayPlanesU(fStrayPlanes);
00591                         }
00592                         if(viewCtr == 1 && status->GetEventDeMuxed() ){ 
00593                                 status->SetVStrayPlanes(fStrayPlanes);
00594                                 status->SetVValidPlanes(validPlaneItr.SizeSelect());
00595                                 cdlh.SetNumValidPlanesV(validPlaneItr.SizeSelect());
00596                                 cdlh.SetNumStrayPlanesV(fStrayPlanes);
00597                         }
00598 
00599                         //done using the first view, so clear the selection function
00600                         //MSG("AlgDeMuxBeam", Msg::kDebug) <<"Done with View " << viewCtr << endl;
00601                         planeItr.GetSet()->ClearSlice(false);
00602                         planeItr.GetSet()->AdoptSelectKeyFunc(0);
00603                         planeItr.ResetFirst();
00604                         validPlaneItr.GetSet()->ClearSlice(false);
00605                         validPlaneItr.GetSet()->AdoptSelectKeyFunc(0);
00606                         validPlaneItr.Reset();
00607                         vertexPlaneItr.GetSet()->ClearSlice(false);
00608                         vertexPlaneItr.Reset();
00609 
00610                         viewCtr++;
00611                 } //end loop over views
00612 
00613                 //if the event was demuxed, check to see how the demuxing and timing jive
00614                 //otherwise see if you need to set the nonphysical solution flag
00615                 if( status->GetEventDeMuxed() ){
00616                         status->SetAverageTimingOffset(util.CheckFitWithTiming(validPlaneItr));
00617                         cdlh.SetAvgTimeOffset(status->GetAverageTimingOffset());
00618                 }
00619                 else if(status->GetNonPhysicalFailure() )  
00620                         cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kNonPhysicalStripSolution);
00621     
00622         }//end if there was a vertex
00623         else{ 
00624                 status->SetNoVertexFailure(true);
00625                 cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kNoVertex);
00626                 MSG("AlgDeMuxBeam", Msg::kDebug) << "no vertex found for event " 
00627                                                                         << status->GetEventNumber() << endl;
00628         }
00629   
00630         //get rid of the temporary status object if you made it.
00631         if(tempStatus) delete status;
00632 
00633         return;
00634 }

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::kDebug, MSG, and DmxPlane::SetStrips().

Referenced by DeMuxFirstNPlanes(), and DeMuxFirstNPlanesTest().

00863 {
00864         MSG("AlgDeMuxBeam", Msg::kDebug) << "In SetFirstNPlanes " << slope 
00865                                                                         << " " << intercept <<  endl;
00866         
00867         Float_t cog = 0.;
00868         planeItr.ResetFirst();
00869         DmxPlane *currentPlane = 0;
00870         
00871         while( (currentPlane = planeItr()) ){
00872                 
00873                 MSG("AlgDeMuxBeam", Msg::kDebug) << "on plane " << currentPlane->GetPlaneNumber()
00874                                                                                 << endl;
00875                 
00876                 cog = slope*currentPlane->GetZPosition() + intercept;
00877                 //set those strips if the plane isnt a valid muon plane 
00878                 //if(currentPlane->GetPlaneType() == DmxPlaneTypes::kMuon 
00879                 //   || !currentPlane->IsValid() ){
00880                         currentPlane->SetStrips(cog);
00881                         MSG("AlgDeMuxBeam", Msg::kDebug) << "setting first n planes " 
00882                                                                                         << currentPlane->GetPlaneNumber() << " to "
00883                                                                                         << cog << "/" << currentPlane->GetCoG() << endl;
00884                 //}
00885         
00886         }//end loop over planes
00887 
00888         planeItr.ResetFirst();
00889         return;
00890 }

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::kDebug, DmxPlaneTypes::kShower, MSG, and DmxPlane::SetStrips().

Referenced by UseMuonSlidingWindow().

01437 {
01438         MSG("AlgDeMuxBeam", Msg::kDebug) << "In SetMuonRegionWindow" << endl;
01439         
01440         Float_t cog = 0.;
01441         planeItr.ResetFirst();
01442         DmxPlane *currentPlane = 0;
01443         
01444         while( planeItr.IsValid() ){
01445                 //get the current plane
01446                 currentPlane = planeItr.Ptr();
01447                 
01448                 MSG("AlgDeMuxBeam", Msg::kDebug) << "on plane " 
01449                                                                                 << currentPlane->GetPlaneNumber()
01450                                                                                 << endl;
01451                 
01452                 cog = slope*currentPlane->GetZPosition() + intercept;
01453                 //set those strips if the plane isnt a valid muon plane in the 
01454                 //muon track or it has been determined that the plane shouldnt be a 
01455                 //valid muon plane
01456                 if(currentPlane->GetPlaneNumber() >= fMuonStartPlaneNumber
01457                    && (!currentPlane->IsValid()
01458                            || currentPlane->GetPlaneType() == DmxPlaneTypes::kShower
01459                            || setAll)
01460                    ){
01461                         currentPlane->SetStrips(cog);
01462                         MSG("AlgDeMuxBeam", Msg::kDebug) << "setting muon region plane " 
01463                                                                                         << currentPlane->GetPlaneNumber() 
01464                                                                                         << " to cog = " << cog << endl;
01465                 }//end if this plane to be set
01466                 
01467                 planeItr.Next();
01468         }//end loop over planes
01469         
01470         planeItr.ResetFirst();
01471         return;
01472 }

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::kDebug, MSG, and DmxPlane::SetStrips().

Referenced by UseShowerSlidingWindow().

01208 {
01209         MSG("AlgDeMuxBeam", Msg::kDebug) << "In SetShowerRegionWindow" << endl;
01210         
01211         Float_t cog = 0.;
01212         planeItr.ResetFirst();
01213         DmxPlane *currentPlane = 0;
01214         
01215         while( (currentPlane = planeItr()) ){
01216                 //get the current plane
01217                 MSG("AlgDeMuxBeam", Msg::kDebug) << "on plane " << currentPlane->GetPlaneNumber()
01218                                                                                 << endl;
01219                 
01220                 cog = slope*currentPlane->GetZPosition() + intercept;
01221                 if(fMuonStartPlaneNumber>0&&currentPlane->GetPlaneNumber()<fMuonStartPlaneNumber)
01222                         currentPlane->SetStrips(cog); 
01223                 else if(fMuonStartPlaneNumber<0)
01224                         currentPlane->SetStrips(cog); 
01225                 
01226                 MSG("AlgDeMuxBeam", Msg::kDebug) << "setting shower region plane " 
01227                                                                                 << currentPlane->GetPlaneNumber() 
01228                                                                                 << " to " << cog << endl;
01229                 
01230         }//end loop over planes
01231         
01232         planeItr.ResetFirst();
01233         return;
01234 }

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

Reimplemented from AlgBase.

Definition at line 637 of file AlgDeMuxBeam.cxx.

00638 {
00639 }

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

Definition at line 1243 of file AlgDeMuxBeam.cxx.

References fFinalShowerRegionIntercept, fFinalShowerRegionSlope, FindFitOverNPlanes(), fInitialMuonRegionIntercept, fInitialMuonRegionSlope, fMuonPlanesInSet, fMuonSlopeLimit, DmxPlane::GetCoG(), DmxStatus::GetEndPlaneNumber(), DmxStatus::GetMultipleMuon(), DmxStatus::GetMuonStartPlaneNumber(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlaneTypes::GetPlaneTypeString(), DmxPlane::GetZPosition(), Msg::kDebug, DmxPlaneTypes::kMuon, MSG, and SetMuonRegionWindow().

Referenced by RunAlg().

01246 {
01247         
01248         MSG("AlgDeMuxBeam", Msg::kDebug) << "In UseMuonSlidingWindowTest" << endl;
01249         
01250         DmxPlane *plane = 0;
01251         
01252         validPlaneItr.ResetFirst();
01253 
01254         Int_t loopCtr = 0;
01255         Int_t ctr = 0;
01256         Float_t slope = 0.;
01257         Float_t intercept = 0.;
01258         Float_t prevSlope = 0.;
01259         Float_t prevIntercept = 0.;
01260         
01261         //arrays to hold the information for fitting the straight lines
01262         //weight each plane by 1/charge
01263         Float_t tpos[10];
01264         Float_t zpos[10];
01265         Float_t weight[10];
01266         Float_t residual[10];
01267         Int_t planeNum[10];
01268         Int_t muonPlanesInSet = fMuonPlanesInSet;
01269         
01270         //find the number of valid muon planes
01271         Int_t validPlaneCnt = 0;
01272         while( (plane = validPlaneItr()) ){
01273                 MSG("AlgDeMuxBeam", Msg::kDebug) << " on plane " << plane->GetPlaneNumber() << " type " 
01274                                                                    << DmxPlaneTypes::GetPlaneTypeString(plane->GetPlaneType()) 
01275                                                                    << endl;
01276                 if(plane->GetPlaneType()==DmxPlaneTypes::kMuon) ++validPlaneCnt;
01277         }
01278         
01279         MSG("AlgDeMuxBeam", Msg::kDebug) << "num valid muon planes = " << validPlaneCnt << endl;
01280         
01281         //make sure you dont have fewer valid planes than the default expected
01282         if(validPlaneCnt<muonPlanesInSet) muonPlanesInSet = validPlaneCnt;
01283         
01284         validPlaneItr.ResetFirst();
01285         
01286         //flag to reset valid muon planes if they are determined to be misreconstructions
01287         bool setAll = false;
01288 
01289         //loop over the valid planes and get the first n valid muon planes
01290         while( validPlaneItr.IsValid() && ctr<muonPlanesInSet){
01291                 plane = validPlaneItr.Ptr();
01292                 if(plane->GetPlaneType()==DmxPlaneTypes::kMuon){
01293                         planeNum[ctr] = plane->GetPlaneNumber();
01294                         zpos[ctr] = plane->GetZPosition();
01295                         tpos[ctr] = plane->GetCoG();
01296                         weight[ctr] = 1./plane->GetPlaneCharge();
01297                         ++ctr;
01298                 }
01299                 validPlaneItr.Next();
01300         }//end loop to fill array for first n planes
01301         
01302         //find the fit for the first n planes and look for 
01303         //planes with large residuals.  repeat until planes with 
01304         //large residuals are removed
01305         FindFitOverNPlanes(zpos,tpos,weight,residual,slope,intercept,ctr,planeNum);
01306         
01307         //loop over the weights and if any are 0, then set the setAll flag
01308         for(Int_t i = 0; i<ctr; ++i){
01309                 if(weight[i]==0.) setAll = true;
01310         }
01311         
01312         if(validPlaneCnt == muonPlanesInSet) planeNum[ctr-1] = status->GetEndPlaneNumber();
01313         if(status->GetMuonStartPlaneNumber() < planeNum[0] && !status->GetMultipleMuon() ) 
01314                 planeNum[0] = status->GetMuonStartPlaneNumber();
01315         
01316         //check that the slope seems reasonable - if not use the last good 
01317         //slope and intercept
01318         if(TMath::Abs(slope) > fMuonSlopeLimit){
01319                 MSG("AlgDeMuxBeam", Msg::kDebug) << "slope = " << slope << "/" << fMuonSlopeLimit 
01320                                                                                 << "reset slope to " << fFinalShowerRegionSlope << endl;
01321                 slope = fFinalShowerRegionSlope;
01322                 intercept = fFinalShowerRegionIntercept;
01323         }
01324         
01325         //set the initial muon region slope and intercept variables
01326         fInitialMuonRegionSlope = slope;
01327         fInitialMuonRegionIntercept = intercept;
01328         prevSlope = slope;
01329         prevIntercept = intercept;
01330         
01331         //set the first planes to the slope and intercept found
01332         planeItr.GetSet()->Slice(planeNum[0], planeNum[ctr-1]+1);
01333         planeItr.GetSet()->Update();
01334         SetMuonRegionWindow(planeItr, slope, intercept, setAll);
01335         planeItr.GetSet()->ClearSlice(false);
01336         
01337         //back the iterator up fMuonPlanesInSet-1 places, unless the are only
01338         //fMuonPlanesInSet.  in that case, just return because the work here is 
01339         //done
01340         if(validPlaneCnt > muonPlanesInSet){
01341                 validPlaneItr.Prev();
01342                 MSG("AlgDeMuxBeam", Msg::kDebug) << validPlaneItr.Ptr()->GetPlaneNumber() << "/" << planeNum[1] << endl; 
01343                 while(validPlaneItr.Ptr()->GetPlaneNumber() > planeNum[1]){
01344                         validPlaneItr.Prev();
01345                         MSG("AlgDeMuxBeam", Msg::kDebug) << validPlaneItr.Ptr()->GetPlaneNumber() << "/" << planeNum[1] << endl; 
01346                 }
01347         }
01348         else return;
01349 
01350         //now demux the event 
01351         while(loopCtr < (validPlaneCnt - muonPlanesInSet) ){
01352                 
01353                 setAll = false;
01354                 
01355                 ctr = 0;
01356                         
01357                 MSG("AlgDeMuxBeam", Msg::kDebug) << "find new window" << endl;
01358                 
01359                 while( validPlaneItr.IsValid() && ctr < muonPlanesInSet){
01360                         
01361                         plane = validPlaneItr.Ptr();
01362                         
01363                         MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber()
01364                                                                                          << endl;
01365                         
01366                         if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){ 
01367         
01368                                 planeNum[ctr] = plane->GetPlaneNumber();
01369                                 tpos[ctr] = plane->GetCoG();
01370                                 zpos[ctr] = plane->GetZPosition();
01371                                 weight[ctr] = 1./plane->GetPlaneCharge();
01372                                 ++ctr;
01373                                 MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber()
01374                                                                                                  << " z = "  << zpos[ctr-1] << " cog = " << tpos[ctr-1]
01375                                                                                          << endl;
01376                         }
01377                         
01378                         validPlaneItr.Next();
01379                 }//end loop over valid planes
01380                 
01381                 //find the straight line to valid muon planes
01382                 FindFitOverNPlanes(zpos, tpos, weight, residual, slope, intercept, ctr, planeNum);
01383                 MSG("AlgDeMuxBeam", Msg::kDebug) << "slope = " << slope 
01384                                                                                  << " intercept = " << intercept << endl;
01385                 
01386                 //loop over the weights and if any are 0, then set the setAll flag
01387                 for(Int_t i = 0; i<ctr; ++i){
01388                         if(weight[i]==0.) setAll = true;
01389                 }
01390                                 
01391                 //check that the slope seems reasonable - if not use the last good 
01392                 //slope and intercept
01393                 if(TMath::Abs(slope) > fMuonSlopeLimit){
01394                         MSG("AlgDeMuxBeam", Msg::kDebug) << "slope = " << slope << "/" << fMuonSlopeLimit 
01395                                                                                         << "reset slope to " << prevSlope << endl;
01396                         slope = prevSlope;
01397                         intercept = prevIntercept;
01398                 }
01399                 
01400                 //check to see if this is the last loop.  if it is, make sure to set all unset planes
01401                 //just use the fit found for the last window of valid planes
01402                 if(loopCtr+1 == validPlaneCnt - muonPlanesInSet) 
01403                         planeNum[ctr-1] = status->GetEndPlaneNumber();
01404                 
01405                 MSG("AlgDeMuxBeam", Msg::kDebug) << "loopCtr = " << loopCtr 
01406                                                                                 << " valid planes = " << validPlaneCnt
01407                                                                                 << " last plane in set = " << planeNum[ctr-1] 
01408                                                                                 << " last plane in event " << status->GetEndPlaneNumber() << endl;
01409                                 
01410                 //set the planes from the last set plane to the current one
01411                 planeItr.GetSet()->Slice(planeNum[ctr-2]+1, planeNum[ctr-1]);
01412                 planeItr.GetSet()->Update();
01413                 SetMuonRegionWindow(planeItr, slope, intercept, setAll);
01414                 planeItr.GetSet()->ClearSlice(false);
01415                 
01416                 ++loopCtr;
01417                 
01418                 prevSlope = slope;
01419                 prevIntercept = intercept;
01420                 
01421                 //back the iterator up 
01422                 validPlaneItr.Prev();
01423                 while( validPlaneItr.Ptr()->GetPlaneNumber() > planeNum[1]) validPlaneItr.Prev();
01424                 
01425         }//end loop over windows
01426         
01427         planeItr.GetSet()->ClearSlice(false);
01428         return;
01429 }

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, fPlanesInSet, DmxPlane::GetCoG(), DmxStatus::GetEndPlaneNumber(), DmxStatus::GetMultipleMuon(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxStatus::GetVertexPlaneNumber(), DmxPlane::GetZPosition(), Msg::kDebug, DmxPlaneTypes::kShower, MSG, and SetShowerRegionWindow().

Referenced by RunAlg().

00978 {
00979         MSG("AlgDeMuxBeam", Msg::kDebug) << "In UseShowerSlidingWindow" << endl;
00980 
00981   //figure out how many windows your are going to need for the shower region.  
00982   //you need to find your lever arm
00983   //for the window too - it is either fPlanesInSet or the size of the validPlaneItr set.
00984   //do like in the cosmic algorithm, start from the vertex and go on.
00985   //only use the window until you get to the start of the muon track.
00986   validPlaneItr.ResetFirst();
00987 
00988   DmxPlane *plane = 0;
00989   
00990   validPlaneItr.ResetFirst();
00991   
00992   Int_t loopCtr = 0;
00993   Int_t ctr = 0;
00994   Float_t slope = 0.;
00995   Float_t intercept = 0.;
00996   
00997   //arrays to hold the information for fitting the straight lines
00998   //weight each plane by 1/charge
00999   Float_t tpos[10];
01000   Float_t zpos[10];
01001   Float_t weight[10];
01002   Float_t residual[10];
01003   Int_t planeNum[10];
01004   Int_t planesInSet = fPlanesInSet;
01005   Float_t minResidual = 1.e20;
01006   Float_t minResidualTPos = 0.;
01007   
01008   //find the number of valid planes - in the shower region both muon and shower type planes 
01009   //are ok to use
01010   Int_t validPlaneCnt = validPlaneItr.SizeSelect();
01011   MSG("AlgDeMuxBeam", Msg::kDebug) << "num valid planes = " << validPlaneCnt << endl;
01012   
01013   if(validPlaneCnt==0){
01014           MSG("AlgDeMuxBeam", Msg::kDebug) << "no valid planes in shower window" << endl;
01015           return;
01016   }
01017   
01018   //make sure you dont have fewer valid planes than the default expected
01019   if(validPlaneCnt<planesInSet) planesInSet = validPlaneCnt;
01020   
01021   validPlaneItr.ResetFirst();
01022   
01023   //flag to reset valid muon planes if they are determined to be misreconstructions
01024   bool setAll = false;
01025   
01026   //loop over the valid planes and get the first n valid muon planes
01027   while( validPlaneItr.IsValid() && ctr<planesInSet){
01028           plane = validPlaneItr.Ptr();
01029           planeNum[ctr] = plane->GetPlaneNumber();
01030           zpos[ctr] = plane->GetZPosition();
01031           tpos[ctr] = plane->GetCoG();
01032           weight[ctr] = 1./plane->GetPlaneCharge();
01033           ++ctr;
01034           validPlaneItr.Next();
01035   }//end loop to fill array for first n planes
01036   
01037   //find the fit for the first window of n planes.  check to see if the 
01038   //average residual is too big, if so then remove the last plane in the 
01039   //window and redo the fit.
01040   FindFitOverNPlanes(zpos, tpos, weight, residual, slope, intercept, ctr, planeNum);
01041           
01042   //check that the slope seems reasonable - if not use the last good 
01043   //slope and intercept
01044   if(TMath::Abs(slope) > fMuonSlopeLimit){
01045           MSG("AlgDeMuxBeam", Msg::kDebug) << "slope = " << slope << "/" << fMuonSlopeLimit 
01046                                                                           << "reset slope to " << fFinalShowerRegionSlope << endl;
01047           slope = fFinalShowerRegionSlope;
01048           intercept = fFinalShowerRegionIntercept;
01049   }
01050   
01051   //if there are no more planes to do the sliding window, check and see if there 
01052   //looks to be a track.  if no track, set the last plane to be set at the end of the 
01053   //event
01054   if(validPlaneCnt == planesInSet && fMuonStartPlaneNumber<0) 
01055         planeNum[ctr-1] = status->GetEndPlaneNumber();
01056   else if(validPlaneCnt == planesInSet && fMuonStartPlaneNumber>0)
01057           planeNum[ctr-1] = fMuonStartPlaneNumber-1;
01058   if(status->GetVertexPlaneNumber() < planeNum[0] && !status->GetMultipleMuon() ) 
01059         planeNum[0] = status->GetVertexPlaneNumber();
01060   
01061   //set the first planes to the slope and intercept found
01062   planeItr.GetSet()->Slice(planeNum[ctr-2], planeNum[ctr-1]+1);
01063   planeItr.GetSet()->Update();
01064   SetShowerRegionWindow(planeItr, slope, intercept);
01065   planeItr.GetSet()->ClearSlice(false);
01066   
01067   //back the iterator up fMuonPlanesInSet-1 places, unless the are only
01068   //fMuonPlanesInSet.  in that case, just return because the work here is 
01069   //done
01070   if(validPlaneCnt > planesInSet){
01071           validPlaneItr.Prev();
01072           MSG("AlgDeMuxBeam", Msg::kDebug) << validPlaneItr.Ptr()->GetPlaneNumber() << "/" << planeNum[1] << endl; 
01073           while(validPlaneItr.Ptr()->GetPlaneNumber() > planeNum[1]){
01074                   validPlaneItr.Prev();
01075                   MSG("AlgDeMuxBeam", Msg::kDebug) << validPlaneItr.Ptr()->GetPlaneNumber() << "/" << planeNum[1] << endl; 
01076           }
01077   }
01078   else return;
01079   
01080   //now demux the event 
01081   while(loopCtr < (validPlaneCnt - planesInSet) ){
01082           
01083           setAll = false;
01084           
01085           ctr = 0;
01086           
01087           MSG("AlgDeMuxBeam", Msg::kDebug) << "find new window" << endl;
01088           
01089           //get the tpos for each of the previously set planes
01090           while( validPlaneItr.IsValid() && ctr < planesInSet-1){
01091                   
01092                   plane = validPlaneItr.Ptr();
01093                   
01094                   MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber()
01095                           << endl;
01096                   
01097                   tpos[ctr] = plane->GetCoG();
01098                   zpos[ctr] = plane->GetZPosition();
01099                   weight[ctr] = 1./plane->GetPlaneCharge();
01100                   planeNum[ctr] = plane->GetPlaneNumber();
01101                   ++ctr;
01102                   MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber()
01103                                                                          << " z = "  << zpos[ctr-1] << " cog = " << tpos[ctr-1]
01104                                                                          << endl;
01105                   validPlaneItr.Next();
01106           }//end loop over valid planes
01107 
01108           //fit the previously set planes
01109           FindFitOverNPlanes(zpos, tpos, weight, residual, slope, intercept, ctr, planeNum);
01110 
01111           //now check out the last plane - if it is a shower plane, then loop over the 3 best hypotheses
01112           //to see which is closest to the line fit to the previously set planes
01113           plane = validPlaneItr.Ptr();
01114           
01115           MSG("AlgDeMuxBeam", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber()
01116                   << endl;
01117           
01118           //get the plane info that is independant of plane type
01119           zpos[ctr] = plane->GetZPosition();
01120           weight[ctr] = 1./plane->GetPlaneCharge();
01121           planeNum[ctr] = plane->GetPlaneNumber();
01122           
01123           minResidual = 1.e20;
01124           //loop over the best 3 hyps to see if any of them fit well with the track
01125           for(Int_t i = 0; i<3; ++i){
01126                   
01127                   tpos[ctr] = plane->GetCoG();
01128                   minResidualTPos = tpos[ctr];
01129                   
01130                   if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
01131                           if(i<1) tpos[ctr] = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
01132                           else if(i<2) tpos[ctr] = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
01133                           else if(i<3) tpos[ctr] = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
01134                   }
01135                   //muon planes only have one possible cog
01136                   else break;
01137                   
01138                   residual[ctr] = TMath::Abs(tpos[ctr] - slope*zpos[ctr]+intercept);
01139 
01140                   //check the residual for the last plane
01141                   if(residual[ctr]<minResidual){
01142                           minResidual = residual[ctr];
01143                           minResidualTPos = tpos[ctr];
01144                   }
01145           }//end loop over best 3 hypotheses for this plane
01146 
01147           //use the transverse position with the minimum residual
01148           tpos[ctr] = minResidualTPos;
01149           
01150           //increment the ctr for later
01151           ++ctr;
01152           
01153           //found the best fit of the hypotheses, so now fit the entire set
01154           FindFitOverNPlanes(zpos, tpos, weight, residual, slope, intercept, ctr, planeNum);
01155           MSG("AlgDeMuxBeam", Msg::kDebug) << "fit over all planes in window slope = " << slope 
01156                                                                           << " intercept = " << intercept << endl;
01157         
01158           //check that the slope seems reasonable - if not use the last good 
01159           //slope and intercept
01160           if(TMath::Abs(slope) > fMuonSlopeLimit){
01161                   MSG("AlgDeMuxBeam", Msg::kDebug) << "slope = " << slope << "/" << fMuonSlopeLimit 
01162                                                                                   << "reset slope to " << fFinalShowerRegionSlope << endl;
01163                   slope = fFinalShowerRegionSlope;
01164                   intercept = fFinalShowerRegionIntercept;
01165           }
01166 
01167           //check to see if this is the last loop.  if it is, make sure to set all unset planes
01168           //just use the fit found for the last window of valid planes - if no track region 
01169           //is present, just set all the remaining planes in the detector - what can it hurt?
01170           if(loopCtr+1 == validPlaneCnt - planesInSet && fMuonStartPlaneNumber<0.) 
01171                   planeNum[ctr-1] = 486;
01172           else if(loopCtr+1 == validPlaneCnt - planesInSet && fMuonStartPlaneNumber>0.) 
01173                   planeNum[ctr-1] = fMuonStartPlaneNumber-1;
01174           
01175           MSG("AlgDeMuxBeam", Msg::kDebug) << "loopCtr = " << loopCtr 
01176                                                                           << " valid planes = " << validPlaneCnt
01177                                                                           << " last plane in set = " << planeNum[ctr-1] 
01178                                                                           << " last plane in event " << status->GetEndPlaneNumber() << endl;
01179           
01180           planeItr.GetSet()->Slice(planeNum[ctr-2]+1, planeNum[ctr-1]+1);
01181           planeItr.GetSet()->Update();
01182           SetShowerRegionWindow(planeItr, slope, intercept);
01183           planeItr.GetSet()->ClearSlice(false);
01184           
01185           ++loopCtr;
01186           
01187           //set the final shower region slope and intercept to the current values
01188           fFinalShowerRegionSlope = slope;
01189           fFinalShowerRegionIntercept = intercept;
01190           
01191           //back the iterator up 
01192           validPlaneItr.Prev();
01193           while( validPlaneItr.Ptr()->GetPlaneNumber() > planeNum[1]) validPlaneItr.Prev();
01194           
01195   }//end loop over windows
01196   
01197   planeItr.GetSet()->ClearSlice(false);
01198     
01199   return;
01200 }


Member Data Documentation

Float_t AlgDeMuxBeam::fDigitResetCut [private]

Definition at line 44 of file AlgDeMuxBeam.h.

Referenced by ReconcileShowerAndMuonRegions(), and RunAlg().

Definition at line 43 of file AlgDeMuxBeam.h.

Referenced by ReconcileShowerAndMuonRegions(), and UseMuonSlidingWindow().

Definition at line 42 of file AlgDeMuxBeam.h.

Referenced by ReconcileShowerAndMuonRegions(), and UseMuonSlidingWindow().

Definition at line 45 of file AlgDeMuxBeam.h.

Referenced by FindFitOverNPlanes(), and RunAlg().

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().

Definition at line 48 of file AlgDeMuxBeam.h.

Referenced by DeMuxFirstNPlanesTest(), and RunAlg().

Definition at line 50 of file AlgDeMuxBeam.h.

Referenced by RunAlg().

Definition at line 49 of file AlgDeMuxBeam.h.

Referenced by RunAlg().

Definition at line 51 of file AlgDeMuxBeam.h.

Referenced by RunAlg(), and UseMuonSlidingWindow().

Float_t AlgDeMuxBeam::fMuonSlopeLimit [private]

Definition at line 52 of file AlgDeMuxBeam.h.

Referenced by RunAlg(), UseMuonSlidingWindow(), and UseShowerSlidingWindow().

Definition at line 54 of file AlgDeMuxBeam.h.

Referenced by ReconcileShowerAndMuonRegions(), and RunAlg().

Definition at line 58 of file AlgDeMuxBeam.h.

Referenced by DeMuxFirstNPlanesTest(), and RunAlg().

Int_t AlgDeMuxBeam::fPlanesInSet [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().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1