AlgDeMuxCosmics Class Reference

#include <AlgDeMuxCosmics.h>

Inheritance diagram for AlgDeMuxCosmics:
AlgBase

List of all members.

Public Member Functions

 AlgDeMuxCosmics ()
virtual ~AlgDeMuxCosmics ()
virtual void Trace (const char *c) const
virtual void RunAlg (AlgConfig &acd, CandHandle &ch, CandContext &cx)

Private Member Functions

void FindFit (DmxPlaneItr &planeItr, Double_t &prevChiSq, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &a4, Int_t *hypotheses, Int_t direction, Int_t leverArm, Float_t offset)
Bool_t FindPlanesToDropFromFit (DmxPlaneItr &planeItr, Double_t a1, Double_t a2, Double_t a3, Double_t a4, Int_t direction, Int_t leverArm, Float_t offset)
void SetPlanesToDeterminedFit (DmxPlaneItr &planeItr, Double_t a1, Double_t a2, Double_t a3, Double_t a4, Float_t offset)
void SetPlanesToDeterminedFit (DmxPlaneItr &planeItr, Double_t a, Double_t b)
void FindCosmicSolution (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, DmxStatus *status, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &a4, Double_t &chiSq, Float_t offset)
void FindWindowCosmicSolution (DmxPlaneItr &validPlaneItr, DmxStatus *status, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &a4, Double_t &chiSq, Float_t offset)
void UseSingleFit (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, DmxStatus *status)
void UseSlidingWindow (DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, DmxPlaneItr &windowPlaneItr, DmxStatus *status, Int_t vertex, Int_t endPlane)

Private Attributes

Double_t fRequiredMatedSignalRatio
Double_t fSlopeRMS
Int_t fStrayPlanes
UInt_t fPlanesInSet
Int_t fStrayCut

Detailed Description

Definition at line 20 of file AlgDeMuxCosmics.h.


Constructor & Destructor Documentation

AlgDeMuxCosmics::AlgDeMuxCosmics (  ) 

Definition at line 74 of file AlgDeMuxCosmics.cxx.

00074                                  :
00075   fRequiredMatedSignalRatio(0.33),
00076   fSlopeRMS(0.),
00077   fStrayPlanes(0),
00078   fPlanesInSet(6),
00079   fStrayCut(6)
00080 {
00081   //default constructor
00082 }

AlgDeMuxCosmics::~AlgDeMuxCosmics (  )  [virtual]

Definition at line 86 of file AlgDeMuxCosmics.cxx.

00087 {
00088 
00089 }


Member Function Documentation

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

Definition at line 559 of file AlgDeMuxCosmics.cxx.

References FindFit(), FindPlanesToDropFromFit(), fPlanesInSet, DmxStatus::GetEventDirection(), and DmxUtilities::IsValidFit().

Referenced by UseSingleFit(), and UseSlidingWindow().

00563 {
00564   //variables to keep track of the best reconstruction
00565   Double_t fitA1 = 0.;
00566   Double_t fitA2 = 0.;
00567   Double_t fitA3 = 0.;
00568   Double_t fitA4 = 0.;
00569   Double_t useA1 = -10000.;
00570   Double_t useA2 = -10000.;
00571   Double_t useA3 = -10000.;
00572   Double_t useA4 = -10000.;
00573   Double_t bestChi2 = 1.e20;
00574   Double_t chi2 = 1.000001e20;
00575   Int_t direction = status->GetEventDirection();
00576 
00577   DmxUtilities util;
00578 
00579   //make the lever arm for the demuxer the # of valid planes or 6, whichever
00580   //is smaller
00581   UInt_t leverArm = validPlaneItr.SizeSelect();
00582     
00583   if(leverArm > fPlanesInSet){ leverArm = fPlanesInSet;}
00584   
00585   //MSG("DMXX", Msg::kDebug) << "in FindCosmicSolution" << endl;
00586 
00587   //MSG("DMX", Msg::kDebug) << "\tlever arm = " << leverArm << endl;
00588 
00589   //declare an array to tell the fitter which hypotheses to use
00590   Int_t hypsToUse[] = {0,0,0,0,0,0};
00591   
00592   for(Int_t plane0Hyp = 1; plane0Hyp < 4; plane0Hyp++){
00593     hypsToUse[0] = plane0Hyp;
00594     for(Int_t plane1Hyp = 1; plane1Hyp < 4; plane1Hyp++){
00595       hypsToUse[1] = plane1Hyp;
00596       if(leverArm == 2){
00597         
00598         FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
00599         
00600         if(chi2 < bestChi2){
00601           bestChi2 = chi2;
00602           useA1 = fitA1;
00603           useA2 = fitA2;
00604           useA3 = fitA3;
00605           useA4 = fitA4;
00606         }//end if its a better straight line fit
00607       }
00608       else if(leverArm >= 3){
00609         for(Int_t plane2Hyp = 1; plane2Hyp < 4; plane2Hyp++){
00610           hypsToUse[2] = plane2Hyp;
00611           if(leverArm == 3){
00612             
00613             FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
00614             
00615             if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
00616               bestChi2 = chi2;
00617               useA1 = fitA1;
00618               useA2 = fitA2;
00619               useA3 = fitA3;
00620               useA4 = fitA4;
00621 
00622               //see if you can drop some planes and improve the fit
00623               if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){
00624                 FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
00625                 
00626                 if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
00627                   bestChi2 = chi2;
00628                   useA1 = fitA1;
00629                   useA2 = fitA2;
00630                   useA3 = fitA3;
00631                   useA4 = fitA4;
00632                 }
00633                 validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
00634               }//end if planes should be dropped because they are way off
00635             }//end if its a better straight line fit
00636           }
00637           else if(leverArm >= 4){
00638             for(Int_t plane3Hyp = 1; plane3Hyp < 4; plane3Hyp++){
00639               hypsToUse[3] = plane3Hyp;
00640               if(leverArm == 4){
00641                 
00642                 FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
00643                 
00644                 if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
00645                   bestChi2 = chi2;
00646                   useA1 = fitA1;
00647                   useA2 = fitA2;
00648                   useA3 = fitA3;
00649                   useA4 = fitA4;
00650 
00651                   //see if you can drop some planes and improve the fit
00652                   if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){
00653                     FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
00654                     
00655                     if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
00656                       bestChi2 = chi2;
00657                       useA1 = fitA1;
00658                       useA2 = fitA2;
00659                       useA3 = fitA3;
00660                       useA4 = fitA4;
00661                     }
00662                     validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
00663                   }//end if planes should be dropped because they are way off
00664                 }//end if its a better straight line fit
00665               }
00666               else if(leverArm >= 5){
00667                 for(Int_t plane4Hyp = 1; plane4Hyp < 4; plane4Hyp++){
00668                   hypsToUse[4] = plane4Hyp;
00669                   
00670                   if(leverArm == 5){
00671                     
00672                     FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
00673                     
00674                     if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
00675                       bestChi2 = chi2;
00676                       useA1 = fitA1;
00677                       useA2 = fitA2;
00678                       useA3 = fitA3;
00679                       useA4 = fitA4;
00680 
00681                       //see if you can drop some planes and improve the fit
00682                       if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){
00683                         FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
00684                         
00685                         if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
00686                           bestChi2 = chi2;
00687                           useA1 = fitA1;
00688                           useA2 = fitA2;
00689                           useA3 = fitA3;
00690                           useA4 = fitA4;
00691                         }
00692                         validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
00693                       }//end if planes should be dropped because they are way off
00694                     }//end if its a better straight line fit
00695                   }
00696                   else if(leverArm == 6){
00697                     for(Int_t plane5Hyp = 1; plane5Hyp < 4; plane5Hyp++){
00698                       
00699                       hypsToUse[5] = plane5Hyp;
00700                       FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
00701                       if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2, fitA3, fitA4, offset) ){
00702                         bestChi2 = chi2;
00703                         useA1 = fitA1;
00704                         useA2 = fitA2;  
00705                         useA3 = fitA3;
00706                         useA4 = fitA4;
00707                      
00708                         //see if you can drop some planes and improve the fit
00709                         if( FindPlanesToDropFromFit(validPlaneItr, fitA1, fitA2, fitA3, fitA4, direction, leverArm, offset) ){
00710                           FindFit(validPlaneItr, chi2, fitA1, fitA2, fitA3, fitA4, hypsToUse, direction, leverArm, offset);
00711                           if(chi2 < bestChi2 && util.IsValidFit(planeItr, fitA1, fitA2,  fitA3, fitA4, offset) ){
00712                             bestChi2 = chi2;
00713                             useA1 = fitA1;
00714                             useA2 = fitA2;
00715                             useA3 = fitA3;
00716                             useA4 = fitA4;
00717                       
00718                           }
00719                           validPlaneItr.GetSet()->GetMasks().SetAllMasks(kFALSE);
00720                         } //end if planes should be dropped because they are way off
00721                       }//end if its a better straight line fit
00722                     }//end for loop over 6th plane's hyps
00723                   }//end if using 6 planes
00724                 }//end for loop over 5th plane's hyps
00725               }//end if at least 5 planes
00726             } //loop over 4th plane
00727           } //end if at least 4 planes
00728         } //loop over 3rd plane
00729       } //end if at least 3 planes
00730     } //loop over 2nd plane
00731   } //loop over 1st plane
00732   
00733   //force the fit to the one found
00734   if(useA1 != -10000. && useA2 != -10000. && useA3 != 10000. && useA4 != 10000.){
00735     a1 = useA1;
00736     a2 = useA2;
00737     a3 = useA3;
00738     a4 = useA4;
00739     chiSq = bestChi2;
00740   }
00741   
00742   //MSG("DMX", Msg::kDebug) << "finished in FindCosmicsSolution" << endl;
00743   return;
00744 }

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

Definition at line 1477 of file AlgDeMuxCosmics.cxx.

References DmxUtilities::FindLinearFit(), DmxPlane::GetCoG(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneType(), DmxPlane::GetZPosition(), DmxPlane::IsGolden(), and DmxPlaneTypes::kShower.

Referenced by FindCosmicSolution().

01481 {
01482 
01483   //MSG("DMXX", Msg::kDebug) << "in FindFit" << endl;
01484 
01485   //cout << "in find fit for first n planes" << endl;
01486   Double_t a1L = 0.;
01487   Double_t a2L = 0.;
01488   Double_t a1Q = 0.;
01489   Double_t a2Q = 0.;
01490   Double_t a3Q = 0.;
01491   
01492   Double_t chiSqL = 1.e20;
01493   Double_t chiSqQ = 1.e20;
01494 
01495   //declare the variables to do a least squares fit to a straight line.  
01496   //weight is the inverse fraction of the total charge on the plane that 
01497   //is on opposite ends of common strips. multiply the inverse by a 
01498   //sensible number to take account for digits with 1000+ adc's
01499 
01500   //the arrays are the number of planes in the set, ie the lever arm
01501   Double_t x[] = {0.,0.,0.,0.,0.,0.};
01502   Double_t y[] = {0.,0.,0.,0.,0.,0.};
01503   Double_t weight[] = {1.,1.,1.,1.,1.,1.};
01504 
01505   DmxUtilities util;
01506 
01507   //cout << "got util object" << endl;
01508   if(direction == -1){ planeItr.ResetLast(); }
01509   else if( direction == 1){ planeItr.ResetFirst(); }
01510 
01511   //planeItr.ResetFirst();
01512 
01513   Int_t planeCtr = 0;
01514   Int_t arrayCtr = 0;
01515   DmxPlane *plane = 0;
01516 
01517   
01518   //MSG("DMX", Msg::kDebug) << "\t" << hypotheses[0] << "\t" << hypotheses[1] 
01519   //             << "\t" << hypotheses[2] << "\t" << hypotheses[3] 
01520   //             << "\t" << hypotheses[4] << "\t" << hypotheses[5] << endl;
01521   
01522   //loop over the plane iterator to fill the variables
01523 
01524   while( (plane = planeItr()) && planeCtr < leverArm){
01525 
01526     if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){
01527       x[arrayCtr] = plane->GetZPosition() - offset;
01528      
01529       //get the maximum possible weight first, then for each hypothesis multiply by the square of the fraction
01530       weight[arrayCtr] = 1./ TMath::Sqrt(plane->GetPlaneCharge());
01531       if(plane->IsGolden()) weight[arrayCtr] /= TMath::Sqrt(10.);
01532       
01533       y[arrayCtr] = plane->GetCoG();
01534       
01535       //only do this if the shower plane is not golden - if it is, you know where it should be
01536       //reconstructed to
01537       if( plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden()){
01538 
01539         if(hypotheses[planeCtr] == 1){
01540           //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetCoG();
01541           weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()); 
01542         }
01543         //plane->GetCoG() returns the best hypothesis CoG so only change the value of y if you are wanting
01544         //the 2nd or 3rd best hypotheses
01545         else if(hypotheses[planeCtr] == 2){ 
01546           //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetCoG();
01547           y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
01548           weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio());
01549         }
01550         else if(hypotheses[planeCtr] == 3){ 
01551           //MSG("DMX", Msg::kDebug) << "\t" << dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetCoG();
01552           y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
01553           weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());
01554         }
01555       }
01556       ++arrayCtr;
01557       
01558     }//end if the plane isnt marked as sucking && it is in the window bounds
01559     
01560     ++planeCtr;
01561   } //end loop over planes to find variables for fit
01562 
01563   //cout << "filled arrays" << endl;
01564 
01565   Double_t chiSq = 0.;
01566   
01567   //use the DmxUtilities function to find a linear fit
01568   util.FindLinearFit(x, y, weight, arrayCtr, a1L, a2L, chiSqL);
01569   //util.FindQuadraticFit(x, y, weight, arrayCtr, a1Q, a2Q, a3Q, chiSqQ);
01570   //util.FindCubicFit(x, y, weight, arrayCtr, a1, a2, a3, a4, chiSq);
01571 
01572   //MSG("DMXX", Msg::kDebug) << "\tlinear fit " << chiSqL << "\t" << a1L << "\t" << a2L << endl;
01573   //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q 
01574   //                      << "\t" << a2Q << "\t" << a3Q << endl;
01575 
01576   //pick the better fit
01577   if(chiSqL<=chiSqQ){
01578     prevChiSq = chiSqL;
01579     a1 = a1L;
01580     a2 = a2L;
01581     a3 = 0.;
01582     a4 = 0.;
01583   }
01584   else{
01585     prevChiSq = chiSqQ;
01586     a1 = a1Q;
01587     a2 = a2Q;
01588     a3 = a3Q;
01589     a4 = 0.;
01590   }
01591   //MSG("DMX", Msg::kDebug) << "\tinitial fit " << a1 << "\t" << a2 << "\t" << a3 << "\t" << a4 << chiSq << endl; 
01592  
01593   planeItr.Reset();
01594   planeCtr = 0;
01595   arrayCtr = 0;
01596   //cout << "reset the iterator" << endl;
01597   
01598   //redo the fit as many times as there are planes, dropping each plane in 
01599   //turn to see if it produces a better chi^2 value - only do it if you have 
01600   //more than 3 planes - ie you can drop one and still do a reasonable fit
01601   //2 planes = great straight line every time.
01602   if(leverArm > 3){
01603     a1L = 0.;
01604     a2L = 0.;
01605     a1Q = 0.;
01606     a2Q = 0.;
01607     a3Q = 0.;
01608 
01609     Double_t fitA1 = 0.;
01610     Double_t fitA2 = 0.;
01611     Double_t fitA3 = 0.;
01612     Double_t fitA4 = 0.;
01613 
01614     for(Int_t i = 0; i < leverArm; i++){
01615       arrayCtr = 0;
01616 
01617       //cout <<"in loop " << i << endl;
01618       //loop over the plane iterator to fill the variables
01619       while( (plane = planeItr()) && planeCtr < leverArm){
01620         
01621         //use those planes not marked as kTRUE - a mark of kTRUE means to drop the plane
01622         //from the set when finding the fit. ie if its true, the plane truely sucks
01623         if( !planeItr.GetSet()->GetMasks().GetMask(plane) ){
01624 
01625           if(i != planeCtr){
01626             x[arrayCtr] = plane->GetZPosition() - offset;
01627             y[arrayCtr] = plane->GetCoG();
01628 
01629             //cout << "filling arrays " << arrayCtr << " " << x[arrayCtr] << " " << y[arrayCtr] 
01630             // << " " << weightSq[arrayCtr] << endl;
01631 
01632             weight[arrayCtr] = 1. / TMath::Sqrt(plane->GetPlaneCharge());
01633             if(plane->IsGolden()) weight[arrayCtr] /= TMath::Sqrt(10.);
01634       
01635             //only make changes if not a golden plane
01636             if( plane->GetPlaneType() == DmxPlaneTypes::kShower  && !plane->IsGolden()){
01637               //only change y if the hypothesis isnt the best one - the unset plane returns the best
01638               //cog from the GetCoG() method
01639               
01640               if(hypotheses[planeCtr] == 1){
01641                 weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio());
01642               }
01643               else if(hypotheses[planeCtr] == 2){ 
01644                 y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
01645                 weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio());
01646               }
01647               else if(hypotheses[planeCtr] == 3){ 
01648                 y[arrayCtr] = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
01649                 weight[arrayCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());
01650               }
01651             }
01652 
01653             ++arrayCtr;
01654             //else{ MSG("DMX", Msg::kDebug) << "\tdrop plane " << plane->GetPlaneNumber() << endl; }
01655           } //end if plane is not the dropped one
01656         }//end if plane is in window bounds
01657         ++planeCtr;
01658       } //end loop over planes to find variables for fit
01659             
01660       //use the DmxUtilities function to find a linear fit
01661       
01662       util.FindLinearFit(x, y, weight, arrayCtr, a1L, a2L, chiSqL);
01663       //util.FindQuadraticFit(x, y, weight, arrayCtr, a1Q, a2Q, a3Q, chiSqQ);
01664       //util.FindCubicFit(x, y, weight, arrayCtr, fitA1, fitA2, fitA3, fitA4, chiSq);
01665             
01666       //MSG("DMXX", Msg::kDebug) << "\tlinear fit " << chiSqL << "\t" << a1L << "\t" << a2L << endl;
01667       //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q 
01668       //                              << "\t" << a2Q << "\t" << a3Q << endl;
01669       //pick the better fit
01670       if(chiSqL<=chiSqQ){
01671         chiSq = chiSqL;
01672         fitA1 = a1L;
01673         fitA2 = a2L;
01674       }
01675       else{
01676         chiSq = chiSqQ;
01677         fitA1 = a1Q;
01678         fitA2 = a2Q;
01679         fitA3 = a3Q;
01680       }
01681     
01682       planeItr.Reset();
01683 
01684       planeCtr = 0;
01685       if(chiSq < prevChiSq){
01686         prevChiSq = chiSq;
01687         a1 = fitA1;
01688         a2 = fitA2;
01689         a3 = fitA3;
01690         a4 = fitA4;
01691         //MSG("DMX", Msg::kDebug) << "\trefined fit, drop plane " << i 
01692         //             << "\t" << prevChiSq << "\t" << a << "\t" << b << endl;
01693       }
01694       
01695     }//end for loop to improve fit by dropping a plane
01696   }//end if enough planes to improve fit
01697 
01698   planeItr.ResetFirst();
01699   //MSG("DMX", Msg::kDebug) << "\tfinal fit" << "\t" << a << "\t" << b << "\t" << prevChiSq << endl; 
01700 
01701   //cout << "finished find fit for first n planes" << endl;
01702   
01703   return;
01704 }

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

Definition at line 1790 of file AlgDeMuxCosmics.cxx.

References DmxPlane::GetCoG(), DmxPlane::GetPlaneType(), DmxPlane::GetZPosition(), DmxPlane::IsGolden(), DmxPlaneTypes::kMuon, DmxPlaneTypes::kShower, and DmxPlane::SetGolden().

Referenced by FindCosmicSolution().

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

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

Definition at line 747 of file AlgDeMuxCosmics.cxx.

References DmxUtilities::FindLinearFit(), fPlanesInSet, DmxPlane::GetCoG(), DmxStatus::GetEventDirection(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlane::GetZPosition(), DmxPlane::IsGolden(), Msg::kDebug, DmxPlaneTypes::kShower, MSG, and DmxPlane::SetGolden().

Referenced by UseSlidingWindow().

00753 {
00754   //MSG("DMXX", Msg::kDebug) << "in FindWindowCosmicSolution" << endl;
00755 
00756   //variables to keep track of the best reconstruction
00757   Double_t fitA1L = 0.;
00758   Double_t fitA2L = 0.;
00759   Double_t chi2L = 1.0e20;
00760 //   Double_t fitA1Q = 0.;
00761 //   Double_t fitA2Q = 0.;
00762 //   Double_t fitA3Q = 0.;
00763 //   Double_t chi2Q = 1.0e20;
00764   Double_t fitA1 = 0.;
00765   Double_t fitA2 = 0.;
00766   Double_t fitA3 = 0.;
00767   Double_t fitA4 = 0.;
00768   Double_t chi2 = 1.0e20;
00769   Double_t chi2Fit = 1.0e20;
00770 
00771   Int_t direction = status->GetEventDirection();
00772 
00773   //make the lever arm for the demuxer the # of valid planes or 6, whichever
00774   //is smaller
00775   Int_t leverArm = fPlanesInSet;
00776 
00777   //offset is the lowest z pos in the window. subtract it from each x array value
00778   //- this has the effect of changing your y axis location 
00779   //relative to the z = 0 location so that when doing higher order fits, you 
00780   //dont get really stupid numbers
00781  
00782   //declare the variables to do a least squares fit to a straight line.  
00783   //use three chiSq values to keep track of the 3 different hyps in the
00784   //unset plane - also have 3 different y values for those hyps and 3
00785   //different weights
00786   Double_t x[] = {0.,0.,0.,0.,0.,0.};
00787   Double_t y[] = {0.,0.,0.,0.,0.,0.};
00788   Double_t weight[] = {1.,1.,1.,1.,1.,1.};
00789   Double_t weight2 = 1.;
00790   Double_t weight3 = 1.;
00791   Double_t y2 = 0.;
00792   Double_t y3 = 0.;
00793 
00794   //variable to let you know if the unset plane is a golden plane or not
00795   bool lastGolden = false;
00796 
00797   Double_t a1L = 0.;
00798   Double_t a2L = 0.;
00799   Double_t a1LB = 0.;
00800   Double_t a2LB = 0.;
00801   Double_t a1LS = 0.;
00802   Double_t a2LS = 0.;
00803   Double_t a1LT = 0.;
00804   Double_t a2LT = 0.;
00805  //  Double_t a1QB = 0.;
00806 //   Double_t a2QB = 0.;
00807 //   Double_t a3QB = 0.;
00808 //   Double_t a1QS = 0.;
00809 //   Double_t a2QS = 0.;
00810 //   Double_t a3QS = 0.;
00811 //   Double_t a1QT = 0.;
00812 //   Double_t a2QT = 0.;
00813 //   Double_t a3QT = 0.;
00814 //   Double_t a1Q = 0.;
00815 //   Double_t a2Q = 0.;
00816 //   Double_t a3Q = 0.;
00817   
00818   Double_t chiSqL = 1.e20;
00819   Double_t chiSqLB = 1.e20;
00820   Double_t chiSqLS = 1.e20;
00821   Double_t chiSqLT = 1.e20;
00822 //   Double_t chiSqQB = 1.e20;
00823 //   Double_t chiSqQS = 1.e20;
00824 //   Double_t chiSqQT = 1.e20;
00825 //   Double_t chiSqQ = 1.e20;
00826 
00827   Double_t a1Fit = 0.;
00828   Double_t a2Fit = 0.;
00829   Double_t a3Fit = 0.;
00830   Double_t a4Fit = 0.;
00831 
00832   
00833   DmxUtilities util;
00834 
00835   if(direction == -1) validPlaneItr.ResetLast();
00836   else if( direction == 1) validPlaneItr.ResetFirst();
00837   
00838   Int_t planeCtr = 0;
00839   DmxPlane *plane = 0;
00840       
00841   while( (plane = validPlaneItr()) && planeCtr < leverArm){
00842 
00843     //cout << "plane = " << plane->GetPlaneNumber() << endl;
00844 
00845     x[planeCtr] = plane->GetZPosition() - offset;
00846     y[planeCtr] = plane->GetCoG();
00847     weight[planeCtr] = 1. / TMath::Sqrt(plane->GetPlaneCharge());
00848 
00849     if(plane->IsGolden()){
00850 
00851       //golden plane, so make the weight smaller to give it more
00852       //clout in the fit
00853       weight[planeCtr] /= TMath::Sqrt(10.);
00854       if(planeCtr == leverArm-1){
00855         lastGolden = true;
00856         y2 = y[planeCtr];
00857         y3 = y[planeCtr];
00858         weight2 = weight[planeCtr];
00859         weight3 = weight[planeCtr];
00860       }
00861     }
00862 
00863     //only make changes if not a golden plane
00864     if(plane->GetPlaneType() == DmxPlaneTypes::kShower && !plane->IsGolden() ){
00865 
00866       //make sure you dont divide by zero when getting the mated signal ratio of
00867       //the set hypothesis - could be that you set the plane to a hypothesis
00868       //with no mated signal, even though it is a valid plane
00869       if(planeCtr < leverArm-1 
00870          && dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio() > 0.){
00871         weight[planeCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio()); 
00872 //      MSG("Dmx", Msg::kDebug) << dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio() << endl;
00873       }
00874       else if( planeCtr == leverArm-1 ){          
00875         
00876         y2 = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
00877         y3 = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
00878 
00879         //dont have to worry about dividing by zero here because you havent
00880         //set the plane yet and you know that the best hypotheses in the planes
00881         //must have at least 0.3 mated signal
00882         weight2 = weight[planeCtr]/TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio()); 
00883         weight3 = weight[planeCtr]/TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio());         
00884 
00885         weight[planeCtr] /= TMath::Sqrt(dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio()); 
00886       }
00887     }
00888 
00889     ++planeCtr;
00890     
00891   } //end loop over planes to find variables for fit
00892       
00893   //do fit using function in DmxUtilities
00894 
00895   //the best hyp fit first
00896   util.FindLinearFit(x, y, weight, leverArm, a1LB, a2LB, chiSqLB);
00897   //util.FindQuadraticFit(x, y, weight, leverArm, a1QB, a2QB, a3QB, chiSqQB);
00898 
00899   //find fit for using the last plane's second best hyp
00900   y[leverArm-1] = y2;
00901   weight[leverArm-1] = weight2;
00902   util.FindLinearFit(x, y, weight, leverArm, a1LS, a2LS, chiSqLS);
00903   //util.FindQuadraticFit(x, y, weight, leverArm, a1QS, a2QS, a3QS, chiSqQS);
00904 
00905   //find fit using the last plane's third best hyp
00906   y[leverArm-1] = y3;
00907   weight[leverArm-1] = weight3;
00908   util.FindLinearFit(x, y, weight, leverArm, a1LT, a2LT, chiSqLT);
00909   //util.FindQuadraticFit(x, y, weight, leverArm, a1QT, a2QT, a3QT, chiSqQT);
00910 
00911   //util.FindCubicFit(x, y, weight, leverArm, a1, a2, a3, a4, chiSq);
00912 
00913  //  MSG("DMXX", Msg::kDebug) << "\tlinear fit1 " << chiSqLB << "\t" << a1LB << "\t" << a2LB << endl;
00914 //   MSG("DMXX", Msg::kDebug) << "\tlinear fit2 " << chiSqLS << "\t" << a1LS << "\t" << a2LS << endl;
00915 //   MSG("DMXX", Msg::kDebug) << "\tlinear fit3 " << chiSqLT << "\t" << a1LT << "\t" << a2LT << endl;
00916   //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q 
00917   //              << "\t" << a2Q << "\t" << a3Q << endl;
00918    
00919   //pick the better linear fit
00920   if(chiSqLB<=chiSqLS && chiSqLB<=chiSqLT){
00921     chiSqL = chiSqLB;
00922     a1L = a1LB;
00923     a2L = a2LB;
00924   }
00925   else if(chiSqLS<chiSqLB && chiSqLS<=chiSqLT){
00926     chiSqL = chiSqLS;
00927     a1L = a1LS;
00928     a2L = a2LS;
00929   }
00930   else if(chiSqLT<chiSqLB && chiSqLT<chiSqLS){
00931     chiSqL = chiSqLT;
00932     a1L = a1LT;
00933     a2L = a2LT;
00934   }
00935 
00936   //pick the better quadratic fit
00937   // if(chiSqQB<=chiSqQS && chiSqQB<=chiSqQT){
00938 //     chiSqQ = chiSqQB;
00939 //     a1Q = a1QB;
00940 //     a2Q = a2QB;
00941 //     a3Q = a3QB;
00942 //   }
00943 //   else if(chiSqQS<chiSqQB && chiSqQS<=chiSqQT){
00944 //     chiSqQ = chiSqQS;
00945 //     a1Q = a1QS;
00946 //     a2Q = a2QS;
00947 //     a3Q = a3QS;
00948 //   }
00949 //   else if(chiSqQT<chiSqQB && chiSqQT<chiSqQS){
00950 //     chiSqQ = chiSqQT;
00951 //     a1Q = a1QT;
00952 //     a2Q = a2QT;
00953 //     a3Q = a3QT;
00954 //   }
00955 
00956   //pick the better fit - linear or quadratic
00957   //if(chiSqL<=chiSqQ){
00958     chi2Fit = chiSqL;
00959     a1Fit = a1L;
00960     a2Fit = a2L;
00961     a3Fit = 0.;
00962     a4Fit = 0.;
00963     //}
00964  //  else{
00965 //     chi2Fit = chiSqQ;
00966 //     a1Fit = a1Q;
00967 //     a2Fit = a2Q;
00968 //     a3Fit = a3Q;
00969 //     a4Fit = 0.;
00970 //   }
00971 
00972     //MSG("DMXX", Msg::kDebug) << "\tinitial fit " << chi2Fit << "\t" << a1Fit << "\t" << a2Fit 
00973     //            << "\t" << a3Fit << "\t" << a4Fit << endl;
00974   
00975   //redo the fit using only the previously set planes and only if you have more than 3 planes in the set
00976   if(leverArm > 3){
00977     
00978     //do fit using function in DmxUtilities
00979     util.FindLinearFit(x, y, weight, leverArm-1, fitA1L, fitA2L, chi2L);
00980     //util.FindQuadraticFit(x, y, weight, planeCtr, fitA1Q, fitA2Q, fitA3Q, chi2Q);
00981     //util.FindCubicFit(x, y, weight, planeCtr, fitA1, fitA2, fitA3, fitA4, chiSq);
00982     
00983     //MSG("DMXX", Msg::kDebug) << "\trefined linear fit " << chi2L << "\t" << fitA1L << "\t" << fitA2L << endl;
00984     //MSG("DMXX", Msg::kDebug) << "\tquadratic fit " << chiSqQ << "\t" << a1Q 
00985     //              << "\t" << a2Q << "\t" << a3Q << endl;
00986   
00987     //pick the better fit - linear or quadratic
00988     //if(chi2L<=chi2Q){
00989     chi2 = chi2L;
00990     fitA1 = fitA1L;
00991     fitA2 = fitA2L;
00992     fitA3 = 0.;
00993     fitA4 = 0.;
00994     //}
00995     //  else{
00996     //       chi2 = chi2Q;
00997     //       fitA1 = fitA1Q;
00998     //       fitA2 = fitA2Q;
00999     //       fitA3 = fitA3Q;
01000     //       fitA4 = 0.;
01001     //     }
01002 
01003     
01004     if( lastGolden ){
01005 
01006       if(direction == 1) validPlaneItr.ResetLast();
01007       else if( direction == -1) validPlaneItr.ResetFirst();
01008       
01009       plane = validPlaneItr();
01010       
01011       MSG("AlgDmx", Msg::kDebug) << "last plane number = " << plane->GetPlaneNumber() 
01012                                  << " this is a golden plane " 
01013                                  << plane->IsGolden() << endl;
01014      
01015       
01016       //unset the IsGolden flag for the last plane if the difference between the
01017       //fit center of gravity and the plane's center of gravity is > 0.5m and
01018       //any of the following
01019       //1. the golden CoG is different from the previous set valid plane's
01020       //   CoG by more than 0.25m/plane (plane to plane spacing is 0.1188m)
01021       //2. chi^2 of fit without the golden plane is an order of magnitude
01022       //   smaller than the chi^2 of the fit with it
01023       //3. the average difference between the set CoG's and the golden CoG and
01024       //   the fit CoG's for all planes in the window is greater than 0.5m
01025 
01026       //cout << x[leverArm-1] << " " << x[leverArm-2] << " " << a1Fit << " " << a2Fit << endl;
01027 
01028       Double_t fitCoG = (a1Fit + a2Fit*(x[leverArm-1])
01029                          + a3Fit*TMath::Power(x[leverArm-1],2) 
01030                          + a4Fit*TMath::Power(x[leverArm-1],3));
01031 
01032       Double_t fitDiff = 0.;
01033       for(Int_t ii = 0; ii<leverArm; ii++){
01034         fitDiff += TMath::Abs((a1Fit + a2Fit*(x[ii]) + a3Fit*TMath::Power(x[ii],2) 
01035                                + a4Fit*TMath::Power(x[ii],3) - y[ii]));
01036       }
01037       Double_t planeDiff = TMath::Abs(x[leverArm-1]-x[leverArm-2])/0.1188;
01038       
01039       MSG("DmxAlg", Msg::kDebug) << "fit is " << fitCoG << " golden CoG = " 
01040                                 << plane->GetCoG() << endl
01041                                 <<" y diff = " << TMath::Abs(y[leverArm-2]-plane->GetCoG()) 
01042                                 << "/" << planeDiff*.25 << endl
01043                                 <<" fitDiff = " << fitDiff <<"/" << 0.5*leverArm << endl
01044                                 << "chi2Fit = " << chi2Fit << "/" << 10.*chi2 << endl;
01045       
01046       if(TMath::Abs(plane->GetCoG()-fitCoG)>0.5 && (fitDiff>0.5*leverArm
01047          || TMath::Abs(y[leverArm-2]-plane->GetCoG())>(0.25*planeDiff)
01048          || chi2Fit>10.*chi2)){
01049         
01050         plane->SetGolden(false);
01051         MSG("DmxAlg", Msg::kDebug) << plane->GetPlaneNumber() << " was golden " 
01052                                    << plane->IsGolden() << endl;
01053         lastGolden = false;
01054       }
01055     }//end check of last plane is golden
01056 
01057     //if the fit found by dropping the last plane has a better chi^2 than the 
01058     //fit using all planes, and that last plane was not golden, then pass the
01059     //fit from the first planes in the window out.  if the last plane was initially
01060     //golden, this only happens for the last time through.
01061     //
01062     if(chi2 < chi2Fit && !lastGolden){
01063     
01064       chi2Fit = chi2;
01065       a1Fit = fitA1;
01066       a2Fit = fitA2;
01067       a3Fit = fitA3;
01068       a4Fit = fitA4;
01069     }
01070    
01071   }//end for loop to improve fit by dropping a plane
01072 
01073   chiSq = chi2Fit;  
01074   a1 = a1Fit;
01075   a2 = a2Fit;
01076   a3 = a3Fit;
01077   a4 = a4Fit;
01078   
01079   //MSG("DMXX", Msg::kDebug) << "\tchosen linear fit " << chiSq << "\t" << a1 << "\t" << a2 << endl;
01080 
01081   validPlaneItr.Reset();
01082 
01083   return;
01084 }

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

Implements AlgBase.

Definition at line 94 of file AlgDeMuxCosmics.cxx.

References bfld::AsString(), DmxUtilities::CheckFitWithTiming(), DmxUtilities::FillPlaneArray(), DmxUtilities::FindEndPlane(), DmxUtilities::FindPlanesOffFit(), DmxUtilities::FindVertexPlane(), fPlanesInSet, fRequiredMatedSignalRatio, fSlopeRMS, fStrayCut, fStrayPlanes, DmxStatus::GetAverageTimingOffset(), CandRecord::GetCandHeader(), CandHandle::GetCandRecord(), VldContext::GetDetector(), Registry::GetDouble(), DmxStatus::GetEndPlaneNumber(), DmxStatus::GetEventDeMuxed(), DmxStatus::GetEventNumber(), DmxStatus::GetFigureOfMeritFailure(), Registry::GetInt(), DmxStatus::GetNonPhysicalFailure(), DmxStatus::GetPlaneArray(), CandHeader::GetSnarl(), DmxStatus::GetVertexPlaneNumber(), RecMinosHdr::GetVldContext(), DmxUtilities::IsOverlappingMultiple(), Msg::kDebug, CandDeMuxDigitList::kEventFailedFilter, KeyOnPlane(), KeyOnUView(), KeyOnValidU(), KeyOnValidV(), KeyOnVView(), Detector::kFar, CandDeMuxDigitList::kMultipleMuonEvent, CandDeMuxDigitList::kNonPhysicalStripSolution, CandDeMuxDigitList::kNoVertex, DmxPlaneTypes::kShower, CandDeMuxDigitList::kTooFewValidPlanes, MSG, DmxStatus::ResetStatus(), DmxStatus::SetAverageTimingOffset(), DmxStatus::SetEndPlaneNumber(), DmxStatus::SetEventNumber(), DmxStatus::SetFigureOfMerit(), DmxStatus::SetMultipleMuon(), DmxStatus::SetNoVertexFailure(), DmxStatus::SetNumberOfPlanes(), SetPlanesToDeterminedFit(), DmxStatus::SetUOverlappingMultiple(), DmxStatus::SetUStrayPlanes(), DmxStatus::SetUValidPlanes(), DmxStatus::SetValidPlanesFailure(), DmxStatus::SetVertexPlaneNumber(), DmxStatus::SetVertexPlaneZPosition(), DmxStatus::SetVOverlappingMultiple(), DmxStatus::SetVStrayPlanes(), DmxStatus::SetVValidPlanes(), UseSingleFit(), and UseSlidingWindow().

00095 {
00096 
00097   assert( ch.InheritsFrom("CandDigitListHandle") );
00098   CandDeMuxDigitListHandle &cdlh = dynamic_cast<CandDeMuxDigitListHandle&>(ch);
00099   
00100   MSG("AlgDeMuxCosmics", Msg::kDebug) << cdlh.GetNDaughters() << " digits in AlgDeMuxCosmics" << endl;
00101   
00102   VldContext vldc = ch.GetCandRecord()->GetCandHeader()->GetVldContext();
00103   if (vldc.GetDetector() != Detector::kFar) {
00104     static int msglimit = 20; // no more than 20 warnings about the detector
00105     if (msglimit) {
00106       MSG("DMX",Msg::kDebug) 
00107         << "AlgDeMuxCosmics::RunAlg() can not DeMux "
00108         << Detector::AsString(vldc.GetDetector())
00109         << " detector.  Skip futher DeMux processing."
00110         << endl;
00111       if (--msglimit == 0) 
00112         MSG("DMX",Msg::kDebug) 
00113           << " ... last message of this type." << endl;
00114     }
00115     return;
00116   }
00117 
00118   //get the AlgConfigDeMux object
00119   fPlanesInSet = acd.GetInt("PlanesInSet");
00120   fRequiredMatedSignalRatio = acd.GetDouble("RatioMatedSignalForValid");
00121   fStrayCut = acd.GetInt("StrayDeltaStripLimit");
00122   //MSG("DMX", Msg::kDebug) << "starting cosmics algorithm" << endl;
00123 
00124   //get the DmxStatus object
00125   // Find the DmxStatus object or create one for needed scratch space
00126   DmxStatus *status = dynamic_cast<DmxStatus *>(gROOT->GetRootFolder()->FindObject("Loon/DeMux/DmxStatus"));
00127   bool tempStatus = false;
00128   if (status == 0) {
00129     MSG("DMXX", Msg::kDebug) << "//root/Loon/DeMux/DmxStatus not found."
00130                               << " Create a temporary DmxStatus." << endl;
00131     status = new DmxStatus;      // Make a temporary DmxStatus if needed
00132     tempStatus = true;
00133   }
00134   else status->ResetStatus(); 
00135 
00136   status->SetEventNumber(ch.GetCandRecord()->GetCandHeader()->GetSnarl());
00137   
00138   //instantiate the DmxInitialize object and then fill DmxStatus with the event information     
00139   DmxUtilities util;
00140   
00141   util.FillPlaneArray(status, cdlh, acd);
00142   const TObjArray *planeArray = status->GetPlaneArray();
00143   
00144   //create the TObjectItr over planes and program it to sort the planes
00145   DmxPlaneItr planeItr(planeArray);
00146   DmxPlaneKeyFunc *pnKF = planeItr.CreateKeyFunc();
00147   pnKF->SetFun(KeyOnPlane);
00148   planeItr.GetSet()->AdoptSortKeyFunc(pnKF);
00149   pnKF = 0;
00150   
00151   planeItr.ResetFirst();
00152 
00153   // while(planeItr.IsValid()){
00154 //     cout << dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber() << endl;
00155 //     planeItr.Next();
00156 //   }
00157 //   planeItr.ResetFirst();
00158 
00159   status->SetNumberOfPlanes(planeItr.SizeSelect());
00160   status->SetVertexPlaneNumber(util.FindVertexPlane(planeItr));
00161   
00162   //demux the event if there is an identified vertex
00163   if( status->GetVertexPlaneNumber() != -1){
00164     
00165     planeItr.ResetFirst();
00166     planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), 500);
00167     planeItr.GetSet()->Update();
00168     //set the end plane number since there was a vertex plane
00169     status->SetEndPlaneNumber(util.FindEndPlane(planeItr));
00170     planeItr.GetSet()->ClearSlice(false);
00171     planeItr.ResetFirst();
00172 
00173     //set the z position of the vertex plane
00174     planeItr.GetSet()->Slice(status->GetVertexPlaneNumber());
00175     planeItr.GetSet()->Update();
00176     planeItr.ResetFirst();
00177     status->SetVertexPlaneZPosition(dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetZPosition());
00178 
00179     planeItr.GetSet()->ClearSlice(false);
00180     planeItr.ResetFirst();
00181 
00182     //make the iterators for the window, valid planes, and vertex
00183     DmxPlaneItr vertexPlaneItr(planeArray);
00184     DmxPlaneItr validPlaneItr(planeArray);
00185     DmxPlaneItr windowPlaneItr(planeArray);
00186        
00187     //create a KeyFunc to sort planes by number
00188     DmxPlaneKeyFunc *vertexPlaneKF = vertexPlaneItr.CreateKeyFunc();
00189     DmxPlaneKeyFunc *validPlaneKF = validPlaneItr.CreateKeyFunc();
00190     DmxPlaneKeyFunc *validWindowPlaneKF = windowPlaneItr.CreateKeyFunc();
00191        
00192     //program the KeyFunc with the sort function
00193     vertexPlaneKF->SetFun(KeyOnPlane);
00194     validPlaneKF->SetFun(KeyOnPlane);
00195     validWindowPlaneKF->SetFun(KeyOnPlane);
00196        
00197     //get the NavSet from the iterator and pass the KeyFunc to it
00198     vertexPlaneItr.GetSet()->AdoptSortKeyFunc(vertexPlaneKF);
00199     validPlaneItr.GetSet()->AdoptSortKeyFunc(validPlaneKF);
00200     windowPlaneItr.GetSet()->AdoptSortKeyFunc(validWindowPlaneKF);
00201        
00202     //clear the KF pointer because we no longer own the KeyFunc
00203     vertexPlaneKF = 0;
00204     validPlaneKF = 0;
00205     validWindowPlaneKF = 0;
00206     planeItr.ResetFirst();
00207     validPlaneItr.ResetFirst();
00208     windowPlaneItr.ResetFirst();
00209         
00210     //now program a key function to select on plane orientation - U first
00211     DmxPlaneKeyFunc *orientValidUKF = validPlaneItr.CreateKeyFunc();
00212     DmxPlaneKeyFunc *orientValidWindowUKF = windowPlaneItr.CreateKeyFunc();
00213     DmxPlaneKeyFunc *orientUKF = planeItr.CreateKeyFunc();
00214     DmxPlaneKeyFunc *orientValidVKF = validPlaneItr.CreateKeyFunc();
00215     DmxPlaneKeyFunc *orientValidWindowVKF = windowPlaneItr.CreateKeyFunc();
00216     DmxPlaneKeyFunc *orientVKF = planeItr.CreateKeyFunc();
00217     
00218     Int_t viewCtr = 0;
00219     while( viewCtr < 2 ){
00220       fSlopeRMS = 0.;
00221       fStrayPlanes = 0;
00222       
00223       if( viewCtr == 0){
00224         //program this function with the select function
00225         orientUKF->SetFun(KeyOnUView);
00226         orientValidUKF->SetFun(KeyOnValidU);
00227         orientValidWindowUKF->SetFun(KeyOnValidU);
00228                 
00229         //adopt it as a selection function
00230         planeItr.GetSet()->AdoptSelectKeyFunc(orientUKF);
00231         orientUKF = 0;
00232         planeItr.Reset();
00233         validPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidUKF);
00234         orientValidUKF = 0;
00235         validPlaneItr.Reset();
00236         windowPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidWindowUKF);
00237         orientValidWindowUKF = 0;
00238         windowPlaneItr.Reset();
00239       }
00240       else if(viewCtr == 1){
00241         //program this function with the select function
00242         orientVKF->SetFun(KeyOnVView);
00243         orientValidVKF->SetFun(KeyOnValidV);
00244         orientValidWindowVKF->SetFun(KeyOnValidV);
00245                 
00246         //adopt it as a selection function
00247         planeItr.GetSet()->AdoptSelectKeyFunc(orientVKF);
00248         orientVKF = 0;
00249         planeItr.Reset();
00250         validPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidVKF);
00251         orientValidVKF = 0;
00252         validPlaneItr.Reset();
00253         windowPlaneItr.GetSet()->AdoptSelectKeyFunc(orientValidWindowVKF);
00254         orientValidWindowVKF = 0;
00255         windowPlaneItr.Reset();
00256       }
00257 
00258       //  validPlaneItr.ResetLast();
00259 //      MSG("AlgDeMuxCosmics", Msg::kDebug)<<"event " << status->GetEventNumber() << endl;
00260 //        while(validPlaneItr.IsValid()){
00261 //      MSG("DMX1", Msg::kDebug)<<"\tvalid plane " 
00262 //                             << dynamic_cast<DmxPlane* >(validPlaneItr.Ptr())->GetPlaneNumber() << endl; 
00263 //      validPlaneItr.Prev();
00264 //        }
00265       
00266 //        validPlaneItr.Reset();
00267       
00268       //  DmxPlane *plane = 0;
00269       
00270 //        while( (plane = dynamic_cast<DmxPlane *>(validPlaneItr())) ){
00271 //      MSG("DMX1", Msg::kDebug) << "\t" << plane->GetPlaneNumber() <<endl;
00272 //        }
00273 //        validPlaneItr.Reset();
00274 
00275       //slice to those planes from the vertex on
00276       planeItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber());
00277       planeItr.GetSet()->Update();
00278 
00279       validPlaneItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber());
00280       validPlaneItr.GetSet()->Update();
00281 
00282       windowPlaneItr.GetSet()->Slice(status->GetVertexPlaneNumber(), status->GetEndPlaneNumber());
00283       windowPlaneItr.GetSet()->Update();
00284       
00285       //if only 2 planes set them to their best hyps
00286       if(validPlaneItr.SizeSelect() == 2){
00287         validPlaneItr.ResetFirst();
00288 
00289         //variables for straight line fit
00290         Float_t cogBeg = -1., zBeg = -1., cogEnd = -1., zEnd = -1.;
00291 
00292         while( validPlaneItr.IsValid() ){
00293           
00294           //only set it if it is a shower plane, muon planes are set in the constructor
00295           if(validPlaneItr.Ptr()->GetPlaneType() == DmxPlaneTypes::kShower
00296              && !validPlaneItr.Ptr()->IsGolden())
00297             dynamic_cast<DmxShowerPlane *>(validPlaneItr.Ptr())->SetStrips("best");
00298 
00299           //cout << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetCoG() << "\t" 
00300           //   << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetPlaneNumber() 
00301           //   << "\t" << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetZPosition() << endl;
00302           
00303           //set the variables for the straighe line fit between the two planes
00304           if(cogBeg == -1.){
00305             cogBeg = validPlaneItr.Ptr()->GetCoG();
00306             zBeg = validPlaneItr.Ptr()->GetZPosition();
00307           }
00308           else{
00309             cogEnd = validPlaneItr.Ptr()->GetCoG();
00310             zEnd = validPlaneItr.Ptr()->GetZPosition();
00311           }
00312           
00313           validPlaneItr.Next();
00314         }//end loop over planes
00315         
00316         //set the nonvalid planes if they exist
00317         if(planeItr.SizeSelect()>2 && cogEnd-cogBeg != 0. && cogBeg != -1. && cogEnd != -1.){
00318           Double_t slope = (cogEnd-cogBeg)/(zEnd - zBeg);
00319           Double_t intercept = cogBeg - slope*zBeg;
00320           //cout << slope << "\t" << intercept << endl;
00321           SetPlanesToDeterminedFit(planeItr, intercept, slope);
00322           fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);
00323         }
00324 
00325         //fStrayPlanes -= util.CheckFit(validPlaneItr);
00326         validPlaneItr.ResetFirst();
00327 
00328       }//end if only 2 valid planes
00329       else if(validPlaneItr.SizeSelect() >2 ){
00330 
00331         //if there are only enough planes for 1 instances of the sliding window + initial fit,
00332         //change fPlanesInSet to give you 3 + the initial fit, but at least 4 always
00333         if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) && fPlanesInSet==6) fPlanesInSet -= 2;
00334         else if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) 
00335                 && fPlanesInSet==5) fPlanesInSet -= 1;
00336 
00337         if(validPlaneItr.SizeSelect() > fPlanesInSet ){
00338           UseSlidingWindow(planeItr, validPlaneItr, windowPlaneItr, status, 
00339                            status->GetVertexPlaneNumber(),
00340                            status->GetEndPlaneNumber());        
00341         }
00342         else if(validPlaneItr.SizeSelect() <= fPlanesInSet){
00343           UseSingleFit(planeItr, validPlaneItr, status);
00344         }
00345         
00346         //check the end points of the event - send the validPlaneItr to the function
00347         //fStrayPlanes -= util.CheckFit(validPlaneItr);
00348         
00349         //what about multiple muons?  look for another vertex after the end plane
00350         planeItr.GetSet()->ClearSlice(false);
00351         validPlaneItr.GetSet()->ClearSlice(false);
00352         windowPlaneItr.GetSet()->ClearSlice(false);
00353         vertexPlaneItr.GetSet()->Slice(status->GetEndPlaneNumber()+1, 500);
00354         vertexPlaneItr.GetSet()->Update();
00355 
00356         //MSG("DMX1", Msg::kDebug)<< "event =  " << status->GetEventNumber() 
00357         //                   << "\tend plane = " << status->GetEndPlaneNumber() << endl;
00358         
00359         Int_t nextVertex = -1;
00360         if(vertexPlaneItr.SizeSelect() > 3){ nextVertex = util.FindVertexPlane(vertexPlaneItr);}
00361         
00362         while( nextVertex != -1){
00363           
00364           status->SetMultipleMuon(true); 
00365           cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kMultipleMuonEvent);
00366           vertexPlaneItr.GetSet()->ClearSlice(false);
00367           vertexPlaneItr.GetSet()->Slice(nextVertex, 500);
00368           vertexPlaneItr.GetSet()->Update();
00369           
00370           //get the next end plane
00371           Int_t endPlane = util.FindEndPlane(vertexPlaneItr);
00372 
00373           //MSG("DMX1", Msg::kDebug)<< "\tin double muon loop, vertex at plane " << nextVertex 
00374           //             << "\tend plane = " << endPlane << "\tplanes in slice = " 
00375           //             << vertexPlaneItr.SizeSelect() << endl;
00376           
00377           planeItr.GetSet()->Slice(nextVertex, endPlane);
00378           planeItr.GetSet()->Update();
00379           
00380           validPlaneItr.GetSet()->Slice(nextVertex, endPlane);
00381           validPlaneItr.GetSet()->Update();
00382           
00383           windowPlaneItr.GetSet()->Slice(nextVertex, endPlane);
00384           windowPlaneItr.GetSet()->Update();
00385                     
00386           //if only 2 planes set them to their best hyps
00387           if(validPlaneItr.SizeSelect() == 2){
00388             validPlaneItr.ResetFirst();
00389             
00390             //variables for straight line fit
00391             Float_t cogBeg = -1., zBeg = -1., cogEnd = -1., zEnd = -1.;
00392             
00393             while( validPlaneItr.IsValid() ){
00394               
00395               //only set it if it is a shower plane, muon planes are set in the constructor
00396               if(validPlaneItr.Ptr()->GetPlaneType() == DmxPlaneTypes::kShower
00397                  && !validPlaneItr.Ptr()->IsGolden())
00398                 dynamic_cast<DmxShowerPlane *>(validPlaneItr.Ptr())->SetStrips("best");
00399               
00400               //set the variables for the straighe line fit between the two planes
00401               if(cogBeg == -1.){
00402                 cogBeg = validPlaneItr.Ptr()->GetCoG();
00403                 zBeg = validPlaneItr.Ptr()->GetZPosition();
00404               }
00405               else{
00406                 cogEnd = validPlaneItr.Ptr()->GetCoG();
00407                 zEnd = validPlaneItr.Ptr()->GetZPosition();
00408               }
00409               
00410               validPlaneItr.Next();
00411             }//end loop over planes
00412             
00413             if(planeItr.SizeSelect()>2 && cogEnd-cogBeg != 0. && cogBeg != -1. && cogEnd != -1.){
00414               Double_t slope = (cogEnd-cogBeg)/(zEnd - zBeg);
00415               Double_t intercept = cogBeg - slope*zBeg;
00416               SetPlanesToDeterminedFit(planeItr, intercept, slope);
00417               fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);
00418             }
00419             
00420             //fStrayPlanes -= util.CheckFit(validPlaneItr);
00421             validPlaneItr.ResetFirst();
00422             
00423           }//end if only 2 valid planes
00424           else if(validPlaneItr.SizeSelect() > 2){
00425             //if there are only enough planes for 1 instances of the sliding window + initial fit,
00426             //change fPlanesInSet to give you 3 + the initial fit
00427             if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) && fPlanesInSet==6) fPlanesInSet -= 2;
00428             else if(validPlaneItr.SizeSelect() == (fPlanesInSet+1) 
00429                     && fPlanesInSet==5) fPlanesInSet -= 1;
00430                     
00431             if(validPlaneItr.SizeSelect() > fPlanesInSet ){
00432               UseSlidingWindow(planeItr, validPlaneItr, windowPlaneItr, status, 
00433                                nextVertex, endPlane);   
00434             }
00435             else if(validPlaneItr.SizeSelect() <= fPlanesInSet ){
00436               UseSingleFit(planeItr, validPlaneItr, status);
00437             }
00438           }//end if more than 2 valid planes in event
00439 
00440           //check the end points of the event - send the validPlaneItr to the function
00441           //fStrayPlanes -= util.CheckFit(validPlaneItr);
00442 
00443           //look for another vertex after the end plane
00444           planeItr.GetSet()->ClearSlice(false);
00445           validPlaneItr.GetSet()->ClearSlice(false);
00446           windowPlaneItr.GetSet()->ClearSlice(false);
00447           vertexPlaneItr.GetSet()->ClearSlice(false);
00448 
00449           vertexPlaneItr.GetSet()->Slice(endPlane+1, 500);
00450           vertexPlaneItr.GetSet()->Update();
00451 
00452           //MSG("DMX1", Msg::kDebug)<< "\tlook for new vertex, slice from " << endPlane+1 
00453           //             << "\tto plane 500, planes in slice = " 
00454           //             << vertexPlaneItr.SizeSelect() << endl;
00455           
00456           vertexPlaneItr.Reset();
00457           
00458           if(vertexPlaneItr.SizeSelect() > 0){
00459             nextVertex = util.FindVertexPlane(vertexPlaneItr);
00460           }
00461           else{ nextVertex = -1; }
00462         }//end loop to find next vertex for spread out multiples
00463       }//end if more than 2 valid planes in set
00464       else if(validPlaneItr.SizeSelect() < 2){ 
00465         MSG("DMX", Msg::kDebug)<< "Event " << status->GetEventNumber() 
00466                                << " not demuxed; only " 
00467                                << validPlaneItr.SizeSelect() 
00468                                << " valid planes in view" << endl;
00469 
00470         status->SetValidPlanesFailure(true); 
00471         cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kTooFewValidPlanes);
00472       } //not enough planes to demux
00473 
00474       //done using the first view, so clear the selection function
00475       //MSG("Dmx", Msg::kDebug) <<"Done with View " << viewCtr << endl;
00476       
00477       status->SetFigureOfMerit(validPlaneItr.SizeSelect(), fStrayPlanes);
00478       if(viewCtr == 0 && status->GetEventDeMuxed() ){ 
00479         status->SetUStrayPlanes(fStrayPlanes);
00480         status->SetUValidPlanes(validPlaneItr.SizeSelect());
00481         cdlh.SetNumValidPlanesU(validPlaneItr.SizeSelect());
00482         cdlh.SetNumStrayPlanesU(fStrayPlanes);
00483         if( status->GetFigureOfMeritFailure() ){
00484           //see if the event failed the figure of Merit test, if so, see if it is an 
00485           //overlapping muon
00486           cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kEventFailedFilter);
00487           if(util.IsOverlappingMultiple(validPlaneItr, 1.*status->GetVertexPlaneNumber(), 
00488                                         acd.GetDouble("HoughInterceptRMSLimit"), 
00489                                         acd.GetDouble("HoughSlopeRMSLimit"),
00490                                         acd.GetDouble("HoughPeakLimit"))){
00491             status->SetUOverlappingMultiple(true);
00492             cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kMultipleMuonEvent);
00493           }
00494         }
00495       }
00496       if(viewCtr == 1 && status->GetEventDeMuxed() ){ 
00497         status->SetVStrayPlanes(fStrayPlanes);
00498         status->SetVValidPlanes(validPlaneItr.SizeSelect());
00499         cdlh.SetNumValidPlanesV(validPlaneItr.SizeSelect());
00500         cdlh.SetNumStrayPlanesV(fStrayPlanes);
00501         //see if the event failed the figure of Merit test, if so, see if it is an 
00502         //overlapping muon
00503         if( status->GetFigureOfMeritFailure() ){
00504           cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kEventFailedFilter);
00505           if(util.IsOverlappingMultiple(validPlaneItr, 1.*status->GetVertexPlaneNumber(), 
00506                                         acd.GetDouble("HoughInterceptRMSLimit"), 
00507                                         acd.GetDouble("HoughSlopeRMSLimit"),
00508                                         acd.GetDouble("HoughPeakLimit"))){
00509             status->SetVOverlappingMultiple(true);
00510             cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kMultipleMuonEvent);
00511           }
00512         }
00513       }
00514       vertexPlaneItr.GetSet()->ClearSlice(false);
00515       planeItr.GetSet()->ClearSlice(false);
00516       planeItr.GetSet()->AdoptSelectKeyFunc(0);
00517       planeItr.Reset();
00518       validPlaneItr.GetSet()->ClearSlice(false);
00519       validPlaneItr.GetSet()->AdoptSelectKeyFunc(0);
00520       validPlaneItr.Reset();
00521       windowPlaneItr.GetSet()->ClearSlice(false);
00522       windowPlaneItr.GetSet()->AdoptSelectKeyFunc(0);
00523       windowPlaneItr.Reset();
00524 
00525       ++viewCtr;
00526     } //end loop over views
00527 
00528     //if the event was demuxed, check to see how the demuxing and timing jive
00529     //otherwise see if you need to set the nonphysical solution flag
00530     if( status->GetEventDeMuxed() ){
00531       status->SetAverageTimingOffset(util.CheckFitWithTiming(validPlaneItr));
00532       cdlh.SetAvgTimeOffset(status->GetAverageTimingOffset());
00533     }
00534     else if(status->GetNonPhysicalFailure() )  
00535       cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kNonPhysicalStripSolution);
00536          
00537     
00538   }//end if there was a vertex
00539   else{ 
00540     status->SetNoVertexFailure(true);
00541     cdlh.SetDeMuxDigitListFlagBit(CandDeMuxDigitList::kNoVertex);
00542     MSG("DMX", Msg::kDebug) << "no vertex found for event " << status->GetEventNumber() << endl;
00543   }
00544 
00545   //get rid of the temporary status object if you made it.
00546   if(tempStatus) delete status;
00547 
00548   return;
00549 }

void AlgDeMuxCosmics::SetPlanesToDeterminedFit ( DmxPlaneItr &  planeItr,
Double_t  a,
Double_t  b 
) [private]

Definition at line 1752 of file AlgDeMuxCosmics.cxx.

References DmxPlane::GetZPosition(), DmxPlane::IsGolden(), and DmxPlane::SetStrips().

01753 {
01754 
01755   //MSG("DMXX", Msg::kDebug) << "in SetPlanesToDeterminedFit" << endl;
01756 
01757   planeItr.ResetFirst();
01758 
01759   DmxPlane *plane = 0;
01760 
01761   while( planeItr.IsValid() ){
01762 
01763     plane = planeItr.Ptr();
01764       
01765     //if(plane->GetPlaneNumber()>=firstPlane && plane->GetPlaneNumber()<=lastPlane){
01766     //MSG("DMXX", Msg::kDebug)<<"\tin SetPlanesToDeterminedFit " << plane->GetPlaneNumber() 
01767     //             <<" has golden flag " << (Int_t)plane->IsGolden() << endl;
01768     Float_t fitCoG = a + (plane->GetZPosition() * b);
01769     //MSG("DMXX", Msg::kDebug)<< a << "\t" << b << "\t" << fitCoG << endl;
01770     
01771     //golden planes have already been set in their constructors
01772     if( !plane->IsGolden() ){
01773       //MSG("DMXX", Msg::kDebug)<<"\tsetting plane " << plane->GetPlaneNumber() << endl;
01774       plane->SetStrips(fitCoG); 
01775     }
01776     //}//end if plane is in the window bounds
01777     
01778     planeItr.Next();
01779   }
01780 
01781   planeItr.Reset();
01782   
01783   return;
01784 }

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

Definition at line 1709 of file AlgDeMuxCosmics.cxx.

References DmxPlane::GetZPosition(), DmxPlane::IsGolden(), and DmxPlane::SetStrips().

Referenced by RunAlg(), UseSingleFit(), and UseSlidingWindow().

01711 {
01712   //MSG("DMXX", Msg::kDebug) << "in SetPlanesToDeterminedFit" << endl;
01713 
01714   planeItr.ResetFirst();
01715 
01716   DmxPlane *plane = 0;
01717 
01718   while( planeItr.IsValid() ){
01719 
01720     plane = planeItr.Ptr();
01721     
01722     //MSG("DMXX", Msg::kDebug)<<"\tin SetPlanesToDeterminedFit " << plane->GetPlaneNumber() 
01723     //             <<" has golden flag " << (Int_t)plane->IsGolden() << endl;
01724     
01725     Float_t fitCoG = (a1 + (plane->GetZPosition()-offset)*a2 
01726                       + TMath::Power(plane->GetZPosition()-offset,2)*a3
01727                       + TMath::Power(plane->GetZPosition()-offset,3)*a4);
01728 
01729     //MSG("DMXX", Msg::kDebug)<< a1 << "\t" << a2 << "\t" << a3 
01730     //             << "\t" << a4 << "\t" << fitCoG 
01731     //             << "\t" << plane->GetZPosition() << endl;
01732     
01733     //golden planes have already been set in their constructors
01734     if( !plane->IsGolden() ){
01735       //MSG("DMXX", Msg::kDebug)<<"\tsetting plane " << plane->GetPlaneNumber() 
01736       //                     << " to " << fitCoG << endl;
01737       plane->SetStrips(fitCoG); 
01738     }
01739         
01740     planeItr.Next();
01741   }
01742 
01743   planeItr.Reset();
01744   
01745   //cout << "finished set planes for fit" << endl;
01746   return;
01747 }

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

Reimplemented from AlgBase.

Definition at line 552 of file AlgDeMuxCosmics.cxx.

00553 {
00554 }

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

Definition at line 1088 of file AlgDeMuxCosmics.cxx.

References FindCosmicSolution(), DmxUtilities::FindPlanesOffFit(), fStrayCut, fStrayPlanes, DmxPlane::GetCoG(), DmxStatus::GetEventNumber(), DmxPlane::GetZPosition(), DmxPlane::IsGolden(), Msg::kDebug, MSG, DmxStatus::SetEventDirection(), DmxPlane::SetGolden(), DmxStatus::SetNonPhysicalFailure(), and SetPlanesToDeterminedFit().

Referenced by RunAlg().

01090 {
01091   //MSG("DMXX", Msg::kDebug) << "in UseSingleFit" << endl;
01092 
01093   DmxUtilities util = DmxUtilities();
01094   
01095   Double_t a1F = -10000.;
01096   Double_t a2F = -10000.;
01097   Double_t a3F = -10000.;
01098   Double_t a4F = -10000.;
01099   Double_t chiSqF = 1.1e20;
01100   Double_t a1B = -10000.;
01101   Double_t a2B = -10000.;
01102   Double_t a3B = -10000.;
01103   Double_t a4B = -10000.;
01104   Double_t chiSqB = 1.1e20;
01105  
01106   validPlaneItr.ResetFirst();
01107   Float_t offset = validPlaneItr.Ptr()->GetZPosition();
01108   validPlaneItr.ResetLast();
01109   //set the event direction forwards
01110   status->SetEventDirection(1);
01111   
01112   //cout << "calling FindCosmicSolution" << endl;
01113   FindCosmicSolution(planeItr, validPlaneItr, status, a1F, a2F, a3F, a4F, chiSqF, offset);
01114   
01115   //set the event direction to backwards
01116   status->SetEventDirection(-1);
01117   
01118   FindCosmicSolution(planeItr, validPlaneItr, status, a1B, a2B, a3B, a4B, chiSqB, offset);
01119   
01120   DmxPlane *plane = 0;
01121   if( chiSqF <= chiSqB && a1F != -10000. && a2F != -10000. 
01122       && a3F != -10000. && a4F != 10000.){
01123     
01124     //loop over the planes and if a golden plane cog is more than 0.25m off from the fit
01125     //set it as non-golden.
01126     validPlaneItr.Reset();
01127 
01128     while( (plane = planeItr()) ){
01129       Double_t fitCoG = (a1F + a2F*plane->GetZPosition() 
01130                          + a3F*TMath::Power(plane->GetZPosition(),2) 
01131                          + a4F*TMath::Power(plane->GetZPosition(),3));
01132       if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false);
01133     }
01134 
01135     SetPlanesToDeterminedFit(planeItr, a1F, a2F, a3F, a4F, offset); 
01136    
01137     fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);
01138   }
01139   else if(chiSqB < chiSqF && a1B != -10000. && a2B != -10000. 
01140           && a3B != 10000. && a4B != 10000.){
01141    
01142     //loop over the planes and if a golden plane cog is more than 0.25m off from the fit
01143     //set it as non-golden.
01144     validPlaneItr.Reset();
01145   
01146     while( (plane = planeItr()) ){
01147       Double_t fitCoG = (a1B + a2B*plane->GetZPosition() 
01148                          + a3B*TMath::Power(plane->GetZPosition(),2) 
01149                          + a4B*TMath::Power(plane->GetZPosition(),3));
01150       if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false);
01151     }
01152 
01153     SetPlanesToDeterminedFit(planeItr, a1B, a2B, a3B, a4B, offset); 
01154     
01155     fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);
01156 }
01157   else if( a1B == -10000. && a2B == -10000. && a3B != 10000. && a4B != 10000.
01158            && a1F == -10000. && a2F == -10000. && a3F != -10000. && a4F != 10000.){
01159     
01160     MSG("DMX", Msg::kDebug)<< "Event " << status->GetEventNumber() 
01161                            << " not demuxed; can't find a physical solution" 
01162                            << endl;
01163     status->SetNonPhysicalFailure(true); 
01164   }
01165   
01166   
01167   return;
01168 }

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

Definition at line 1172 of file AlgDeMuxCosmics.cxx.

References FindCosmicSolution(), DmxUtilities::FindPlanesOffFit(), FindWindowCosmicSolution(), fPlanesInSet, fStrayCut, fStrayPlanes, DmxPlane::GetCoG(), DmxStatus::GetEventNumber(), DmxPlane::GetPlaneNumber(), DmxStatus::GetVertexZPosition(), DmxPlane::GetZPosition(), DmxPlane::IsGolden(), Msg::kDebug, MSG, DmxStatus::SetEventDirection(), DmxPlane::SetGolden(), DmxStatus::SetNonPhysicalFailure(), and SetPlanesToDeterminedFit().

Referenced by RunAlg().

01175 {
01176   //MSG("DMXX", Msg::kDebug) << "in UseSlidingWindow" << endl;
01177 
01178   //start the sliding window - keep the window to 6 valid planes 
01179   //and use the fit parameters you found for the previous set.  if the 
01180   //new plane is a shower plane, set to the hypothesis nearest the fit.  
01181   //if none of the 3 best hypotheses is within 12 strips of the fit, then
01182   //just use the same fit parameters.  redo the fit, and then set all 
01183   //non-set planes in the window. move the the window down another plane 
01184   //and do it again.
01185 
01186   DmxUtilities util;
01187 
01188   UInt_t validPlaneCnt = validPlaneItr.SizeSelect();
01189   UInt_t loopCtr = 0;
01190   UInt_t ctr = 0;
01191   Int_t firstPlane = 0;
01192   Int_t almostLastPlane = 0;
01193   Int_t lastPlane = 0;
01194   Double_t a1 = -10000.;
01195   Double_t a2 = -10000.;
01196   Double_t a3 = -10000.;
01197   Double_t a4 = -10000.;
01198   Double_t chiSq = 1.1e20;
01199   Float_t offset = 0.;
01200   DmxPlane *plane = 0;
01201   Int_t direction = 1;  //status->GetEventDirection();
01202 
01203   // MSG("dmx", Msg::kDebug) << "Event = " << status->GetEventNumber() << "\tvertex = " << vertex 
01204 //               << "\tend = " << endPlane << endl;
01205 //   while(validPlaneItr.IsValid() ){
01206 //     MSG("dmx", Msg::kDebug) << "\t" << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetPlaneNumber()
01207 //                         << "\t" << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->IsGolden()
01208 //                         <<endl;
01209 //     validPlaneItr.Next();
01210 //   }
01211 
01212   validPlaneItr.ResetFirst();
01213   windowPlaneItr.ResetFirst();
01214 
01215   //now demux the event the <= makes sure that you find an initial fit for the first
01216   //set of planes, then go from there
01217   while(loopCtr <= (validPlaneCnt - fPlanesInSet) ){
01218 
01219       ctr = 0;
01220       
01221       //MSG("dmx", Msg::kDebug) << "find new window" << endl;
01222       //MSG("dmx", Msg::kDebug) << "\tdirection = " << direction;
01223 
01224       while( validPlaneItr.IsValid() && ctr < fPlanesInSet){
01225         
01226         plane = validPlaneItr.Ptr();
01227         //MSG("dmx", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber()
01228         //             << endl;
01229         
01230         if(direction == 1){ 
01231           if( ctr == 0){ 
01232             firstPlane = plane->GetPlaneNumber(); 
01233             offset = plane->GetZPosition();
01234           }
01235           if( ctr == fPlanesInSet-2){ almostLastPlane = plane->GetPlaneNumber(); }
01236           if( ctr == fPlanesInSet-1){ lastPlane = plane->GetPlaneNumber(); }
01237           validPlaneItr.Next();
01238         }
01239         else if(direction == -1){ 
01240           if( ctr == 0){ lastPlane = plane->GetPlaneNumber(); }
01241           if( ctr == fPlanesInSet-2){ almostLastPlane = plane->GetPlaneNumber(); }
01242           if( ctr == fPlanesInSet-1){ 
01243             firstPlane = plane->GetPlaneNumber(); 
01244             offset = plane->GetZPosition();
01245           }
01246           validPlaneItr.Prev();
01247         }
01248         ++ctr;
01249         
01250       }
01251 
01252       //MSG("dmx", Msg::kDebug) << "end find new window" << endl;
01253       ctr = 0;
01254             
01255       //MSG("dmx", Msg::kDebug) << "\tfirst plane = " << firstPlane << "\talmost last plane = " 
01256       //                     << almostLastPlane << "\tlast plane = " << lastPlane << endl;
01257       
01258       //move the window
01259       windowPlaneItr.GetSet()->ClearSlice(false);
01260       planeItr.GetSet()->ClearSlice(false);
01261 
01262       windowPlaneItr.GetSet()->Slice(firstPlane, lastPlane);
01263       windowPlaneItr.GetSet()->Update();
01264       planeItr.GetSet()->Slice(firstPlane, lastPlane);
01265       planeItr.GetSet()->Update();
01266 
01267       if( loopCtr == 0){
01268         
01269         //set the initial window - try it both ways and take the best one as
01270         //determined by the chiSq
01271         Double_t a1B = -10000.;
01272         Double_t a2B = -10000.;
01273         Double_t a3B = -10000.;
01274         Double_t a4B = -10000.;
01275         Double_t chiSqB = 1.1e20;
01276         
01277         //set the event direction forwards
01278         status->SetEventDirection(1);
01279 
01280         //you are looking at the first planes of the event, so your offset should be the
01281         //z position of the vertex
01282         Float_t offsetF = status->GetVertexZPosition();
01283 
01284         FindCosmicSolution(planeItr, windowPlaneItr, status, a1, a2, a3, a4, chiSq, offsetF);
01285 
01286         //MSG("dmx", Msg::kDebug) << "\tforward direction" << a1 << "\t" << a2 
01287         //             << "\t" << chiSq << "\t" << offset << endl;
01288         
01289         //set the event direction to backwards
01290         Int_t forwardLastPlane = lastPlane;
01291 
01292         status->SetEventDirection(-1);
01293         direction = -1;
01294         validPlaneItr.ResetLast();
01295         windowPlaneItr.GetSet()->ClearSlice(false);
01296         planeItr.GetSet()->ClearSlice(false);
01297         windowPlaneItr.ResetLast();
01298         ctr = 0;
01299         
01300         //MSG("dmx", Msg::kDebug) << "check reverse direction" << endl;
01301         while( validPlaneItr.IsValid() && ctr < fPlanesInSet){
01302 
01303           plane = validPlaneItr.Ptr();
01304           // MSG("dmx", Msg::kDebug) << "plane number = " << plane->GetPlaneNumber()
01305           //              << endl;
01306           
01307           if( ctr == 0){ lastPlane = plane->GetPlaneNumber(); }
01308           if( ctr == fPlanesInSet-2){ almostLastPlane = plane->GetPlaneNumber(); }
01309           if( ctr == fPlanesInSet-1){ 
01310             firstPlane = plane->GetPlaneNumber(); 
01311             offset = plane->GetZPosition();
01312           }
01313           ++ctr;
01314           validPlaneItr.Prev();
01315         }
01316         //MSG("dmx", Msg::kDebug) << "end check reverse direction" << endl;
01317         ctr = 0;
01318         
01319         //MSG("dmx", Msg::kDebug) << "first plane = " << firstPlane << "\talmost last plane = " 
01320         //              << almostLastPlane << "\tlast plane = " << lastPlane << endl;
01321                 
01322         windowPlaneItr.GetSet()->Slice(firstPlane, lastPlane);
01323         windowPlaneItr.GetSet()->Update();
01324         windowPlaneItr.ResetLast();
01325         planeItr.GetSet()->Slice(firstPlane, lastPlane);
01326         planeItr.GetSet()->Update();
01327 
01328         FindCosmicSolution(planeItr, windowPlaneItr, status, a1B, a2B, a3B, a4B, chiSqB, offset);
01329         
01330         //MSG("dmx", Msg::kDebug) << "\tbackward direction" << a1B << "\t" << a2B 
01331         //             << "\t" << chiSqB << "\t" << offset << endl;
01332         
01333         if( chiSq <= chiSqB && a1 != -10000. && a2 != -10000. && a3 != 10000. && a4 != 10000.){
01334           planeItr.GetSet()->ClearSlice(false);
01335 
01336           //slice from the vertex to the last plane in the forward direction to get them
01337           //all
01338           planeItr.GetSet()->Slice(vertex, forwardLastPlane);
01339           planeItr.GetSet()->Update();
01340 
01341           //loop over the planes and if a golden plane cog is more than 0.25m off from the fit
01342           //set it as non-golden.
01343           windowPlaneItr.GetSet()->ClearSlice(false);
01344           windowPlaneItr.GetSet()->Slice(vertex, forwardLastPlane);
01345           windowPlaneItr.GetSet()->Update();
01346 
01347           while( (plane = windowPlaneItr()) ){
01348             Double_t fitCoG = (a1 + a2*plane->GetZPosition() 
01349                                + a3*TMath::Power(plane->GetZPosition(),2) 
01350                                + a4*TMath::Power(plane->GetZPosition(),3));
01351             if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false);
01352           }
01353 
01354           //MSG("dmx", Msg::kDebug) << "\tloop = " << loopCtr << "\tsetting planes " << vertex 
01355           //     << "\tto " << forwardLastPlane << endl;
01356 
01357           SetPlanesToDeterminedFit(planeItr, a1, a2, a3, a4, offsetF); 
01358 
01359           fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);
01360             
01361           status->SetEventDirection(1);
01362           direction = 1;
01363           validPlaneItr.ResetFirst();
01364           ctr = 0;
01365           
01366           //get the validPlaneItr back to where it should be for this direction
01367           while(ctr<fPlanesInSet){
01368             validPlaneItr.Next();
01369             ++ctr;
01370           }
01371           
01372           ctr = 0;
01373           
01374         }           
01375         else if(chiSqB < chiSq && a1B != -10000. && a2B != -10000. && a3B != 10000. && a4B != 10000.){
01376           
01377           //set a, b, and chiSq to the backwards values
01378           a1 = a1B;
01379           a2 = a2B;
01380           a3 = a3B;
01381           a4 = a4B;
01382           chiSq = chiSqB;
01383 
01384           //slice from the end plane to the first plane
01385           planeItr.GetSet()->ClearSlice(false);
01386           planeItr.GetSet()->Slice(firstPlane, endPlane);
01387           planeItr.GetSet()->Update();
01388 
01389           //loop over the planes and if a golden plane cog is more than 0.25m off from the fit
01390           //set it as non-golden.
01391           windowPlaneItr.ResetFirst();
01392           while( (plane = windowPlaneItr()) ){
01393             Double_t fitCoG = (a1 + a2*plane->GetZPosition() 
01394                                + a3*TMath::Power(plane->GetZPosition(),2) 
01395                                + a4*TMath::Power(plane->GetZPosition(),3));
01396             if(plane->IsGolden() && TMath::Abs(plane->GetCoG()-fitCoG)>0.25) plane->SetGolden(false);
01397           }
01398           //MSG("dmx", Msg::kDebug) << "\tloop = " << loopCtr << "\tsetting planes " << firstPlane 
01399           //             << "\tto " << endPlane << endl;
01400 
01401           SetPlanesToDeterminedFit(planeItr, a1, a2, a3, a4, offset);
01402 
01403           fStrayPlanes += util.FindPlanesOffFit(planeItr, fStrayCut);
01404 
01405           //event direction is already backwards and the validPlaneItr is where it should be
01406         }
01407         else if( a1B == -10000. && a2B == -10000. && a1 == -10000. 
01408                  && a2 == -10000. && a3 == -10000. && a4 == -10000. 
01409                  && a3B == -10000. && a4B == -10000.){
01410            
01411           MSG("dmx", Msg::kDebug)<< "Event " << status->GetEventNumber() 
01412                                  << " not demuxed;" 
01413                                  << " can't find a physical solution" << endl;
01414           status->SetNonPhysicalFailure(true);
01415         }
01416         //MSG("dmx", Msg::kDebug) << "\tloop 0 fit - " << a1 << "\t" << a2 << "\t" << a3 << "\tchi^2 = " 
01417         //             << chiSq << "\tstray planes = " << fStrayPlanes << endl;
01418       }//end if initial set
01419       else{
01420         //set the window
01421         
01422         FindWindowCosmicSolution(windowPlaneItr, status, a1, a2, a3, a4, chiSq, offset);
01423 
01424         //MSG("dmx", Msg::kDebug) << "\tloop " << loopCtr << " " << firstPlane << " " << lastPlane 
01425         //             << "\t" << a1 << "\t" << a2 << "\t" << a3 << "\tchi^2 = " 
01426         //             << chiSq << "\tstray planes = " << fStrayPlanes << endl;
01427         
01428         planeItr.GetSet()->ClearSlice(false);
01429         if(direction == 1 ){
01430           firstPlane = almostLastPlane + 1;
01431           //if its the last loop, select the rest of the unset planes to set, otherwise
01432           //stay in the loop bounds
01433           if(loopCtr == (validPlaneItr.SizeSelect()-fPlanesInSet)) lastPlane = endPlane;
01434         }
01435         else if(direction == -1 ){
01436           lastPlane = almostLastPlane-1;
01437           if(loopCtr == (validPlaneItr.SizeSelect()-fPlanesInSet)) firstPlane = vertex;
01438         }
01439         planeItr.GetSet()->Slice(firstPlane, lastPlane);
01440         planeItr.GetSet()->Update();
01441 
01442         //MSG("dmx", Msg::kDebug) << "\tloop = " << loopCtr << "\tsetting planes " << firstPlane 
01443         //             << "\tto " << lastPlane << endl;
01444         
01445         SetPlanesToDeterminedFit(planeItr, a1, a2, a3, a4, offset);
01446         fStrayPlanes += util.FindPlanesOffFit(planeItr,  fStrayCut);
01447         //status->SetHistogramEntries(b, a, plane->GetPlaneView());
01448       }
01449 
01450 
01451       ctr = 0;
01452 
01453       //set the validPlaneItr to the opposite direction and back up 5 planes
01454       //MSG("dmx", Msg::kDebug) << "move iterator to first plane in new window" << endl;
01455       while(ctr<fPlanesInSet-1 && validPlaneItr.IsValid()){
01456         
01457         //MSG("dmx", Msg::kDebug) << "plane number = " 
01458         //             << dynamic_cast<DmxPlane *>(validPlaneItr.Ptr())->GetPlaneNumber()
01459         //             << endl;
01460         ++ctr;
01461         if(direction == 1){ validPlaneItr.Prev(); }
01462         else if(direction == -1){ validPlaneItr.Next(); }
01463       }
01464       //MSG("dmx", Msg::kDebug) << "end move iterator to first plane in new window" << endl;
01465       ctr = 0;
01466  
01467       //should go through the loop and get the new window
01468       ++loopCtr;
01469   }
01470 
01471   return;
01472 }


Member Data Documentation

Definition at line 37 of file AlgDeMuxCosmics.h.

Referenced by RunAlg().

Double_t AlgDeMuxCosmics::fSlopeRMS [private]

Definition at line 38 of file AlgDeMuxCosmics.h.

Referenced by RunAlg().

Int_t AlgDeMuxCosmics::fStrayCut [private]

Definition at line 41 of file AlgDeMuxCosmics.h.

Referenced by RunAlg(), UseSingleFit(), and UseSlidingWindow().

Definition at line 39 of file AlgDeMuxCosmics.h.

Referenced by RunAlg(), UseSingleFit(), and UseSlidingWindow().


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

Generated on 2 Nov 2017 for loon by  doxygen 1.6.1