Main Page | Modules | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages
DeMux

AlgDeMuxBeam.cxx

Go to the documentation of this file.
00001 
00002 // $Id: AlgDeMuxBeam.cxx,v 1.46 2007/11/11 07:36:22 rhatcher Exp $
00003 //
00004 // AlgDeMuxBeam
00005 //
00006 // An Algorithm class that DeMultiplexes an event and adjusts the
00007 //weights in the PlexSEIdAltLists of the CandDigits from that event
00008 //
00009 // Author:  B. Rebel 11/2000
00011 
00012 #include <cassert>
00013 #include "CandData/CandHeader.h"
00014 #include "Candidate/CandContext.h"
00015 #include "Candidate/CandHandle.h"
00016 #include "CandDigit/CandDeMuxDigitHandle.h"
00017 #include "CandDigit/CandDeMuxDigitListHandle.h"
00018 #include "MessageService/MsgService.h"
00019 #include "Algorithm/AlgConfig.h"
00020 #include "DeMux/DmxStatus.h"
00021 #include "DeMux/DmxUtilities.h"
00022 #include "Algorithm/AlgFactory.h"
00023 #include "Algorithm/AlgHandle.h"
00024 #include "DeMux/AlgDeMuxBeam.h"
00025 #include "DeMux/DmxPlane.h"
00026 #include "DeMux/DmxPlaneTypes.h"
00027 #include "DeMux/DmxShowerPlane.h"
00028 #include "DeMux/DmxMuonPlane.h"
00029 #include "Navigation/NavKey.h"
00030 #include "Navigation/NavSet.h"
00031 #include "TObjArray.h"
00032 #include "UgliGeometry/UgliGeomHandle.h"
00033 #include "UgliGeometry/UgliStripHandle.h"
00034 #include "DeMux/DmxHypothesis.h"
00035 
00036 #include "TFolder.h"
00037 #include "TROOT.h"
00038 #include "TMath.h"
00039 
00040 //define the NavKey to sort the digits by plane
00041 static NavKey KeyFromPlane1( const CandDeMuxDigitHandle *cdh){
00042   return cdh->GetPlexSEIdAltL().GetPlane();}
00043 
00044 //define the NavKey to sort the planes
00045 static NavKey KeyOnPlane1( const DmxPlane *plane){
00046   return plane->GetPlaneNumber();
00047 }
00048 
00049 //define the NavKeys to select U and V planes
00050 NavKey KeyOnUView1( const DmxPlane *plane){
00051   return plane->GetPlaneView() == PlaneView::kU;
00052 }
00053 
00054 NavKey KeyOnVView1( const DmxPlane *plane){
00055   return plane->GetPlaneView() == PlaneView::kV;
00056 }
00057 
00058 //define the NavKeys to select only valid U and V planes
00059 NavKey KeyOnValidU1( const DmxPlane *plane){
00060   return plane->GetPlaneView() == PlaneView::kU && (Int_t)plane->IsValid() == 1;
00061 }
00062 
00063 NavKey KeyOnValidV1( const DmxPlane *plane){
00064   return plane->GetPlaneView() == PlaneView::kV && (Int_t)plane->IsValid() == 1;
00065 }
00066 
00067 ClassImp(AlgDeMuxBeam)
00068   
00069 CVSID("$Id: AlgDeMuxBeam.cxx,v 1.46 2007/11/11 07:36:22 rhatcher Exp $");
00070 
00071 //-----------------------------------------------------------------------
00072 AlgDeMuxBeam::AlgDeMuxBeam() :
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 }
00082 
00083 //-----------------------------------------------------------------------
00084 
00085 AlgDeMuxBeam::~AlgDeMuxBeam()
00086 {
00087 }
00088 
00089 //-----------------------------------------------------------------------
00090 //run the demuxing on an event
00091 void AlgDeMuxBeam::RunAlg(AlgConfig &acd, CandHandle &ch, CandContext & /* cx */)
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 }
00635 
00636 //______________________________________________________________________
00637 void AlgDeMuxBeam::Trace(const char * /* c */) const
00638 {
00639 }
00640 
00641 //-----------------------------------------------------------------------
00642 //this method takes the array of demuxPlanes and demuxCoGs and finds 
00643 //the best location for the event in the detector.  the other arguments
00644 //it takes are an iterator over the plane objects, the array to hold the 
00645 //best values, and the number of planes to do.
00646 void AlgDeMuxBeam::DeMuxFirstNPlanesTest(DmxPlaneItr &planeItr, 
00647                                                                                  DmxPlaneItr &validPlaneItr)
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 }
00856 
00857 //-----------------------------------------------------------------------
00858 //this method sets all the planes in a slice to a given hypothesis. its 
00859 //arguments are a DmxPlaneItr of DmxPlane objects, the bounding planes of the window,
00860 //and the slope and intercept of the fit
00861 void AlgDeMuxBeam::SetFirstNPlanes(DmxPlaneItr &planeItr, Float_t slope, 
00862                                                                    Float_t intercept)
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 }
00891 
00892 //---------------------------------------------------------------------
00893 //this method finds a linear fit for the x and y data and refines the fit
00894 //based on the residuals.  it can iteratively look for planes causing bad 
00895 //fits and remove them to refine the fit.
00896 void AlgDeMuxBeam::FindFitOverNPlanes(Float_t *x, Float_t *y, 
00897                                                                           Float_t *sigma, Float_t *residual, 
00898                                                                           Float_t &slope, Float_t &intercept,
00899                                                                           Int_t numPoints, Int_t *planeNum)
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 }
00973 
00974 //---------------------------------------------------------------------
00975 void AlgDeMuxBeam::UseShowerSlidingWindow(DmxPlaneItr &planeItr, 
00976                                                                                   DmxPlaneItr &validPlaneItr, 
00977                                                                                   DmxStatus *status)
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 }
01201 
01202 //-----------------------------------------------------------------------
01203 //this method sets all the planes in a slice to a given hypothesis. its 
01204 //arguments are a DmxPlaneItr of DmxPlane objects, the bounding planes of the window
01205 //and the center of gravity of the fit
01206 void AlgDeMuxBeam::SetShowerRegionWindow(DmxPlaneItr &planeItr, Float_t slope, 
01207                                                                                  Float_t intercept)
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 }
01235 
01236 //-----------------------------------------------------------------------------
01237 //method for setting planes in the muon region of the event
01238 //fit a straight line to a window of n valid muon planes.  the first 
01239 //n-1 valid muon planes have been previously set, the n-th plane should 
01240 //lie along the line defined by the first planes.  if it does not, it is 
01241 //not really a valid plane and should be set along the line defined by the 
01242 //previously set planes
01243 void AlgDeMuxBeam::UseMuonSlidingWindow(DmxPlaneItr &planeItr, 
01244                                                                                 DmxPlaneItr &validPlaneItr, 
01245                                                                                 DmxStatus *status)
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 }
01430 
01431 //-----------------------------------------------------------------------
01432 //this method sets all the planes in a slice to a given hypothesis. its 
01433 //arguments are a DmxPlaneItr of DmxPlane objects, the bounding planes of the window,
01434 //and the slope and intercept of the fit
01435 void AlgDeMuxBeam::SetMuonRegionWindow(DmxPlaneItr &planeItr, Float_t slope, 
01436                                                                            Float_t intercept, bool setAll)
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 }
01473 
01474 //----------------------------------------------------------------------
01475 //this is a standard least squares linear fit, see bevington eq. (6.12) for example
01476 void AlgDeMuxBeam::RegionLinearFit(Float_t &slope, Float_t &intercept,
01477                                                                            Float_t &chiSqr,
01478                                                                            Float_t *x, Float_t *y, Float_t *sigma,
01479                                                                            Float_t *residuals, Int_t numPoints)
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 }
01536 
01537 //---------------------------------------------------------------------
01538 //this method checks the DeMux solution in the region of the event defined 
01539 //by validPlaneItr by seeing how much of the total signal from valid planes 
01540 //in the region is mated - both muon and shower planes
01541 Float_t AlgDeMuxBeam::FindRegionMatedChargeFraction(DmxPlaneItr &validPlaneItr)
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 }
01593 
01594 //---------------------------------------------------------------------
01595 //this method compares the predicted transverse position of both the 
01596 //muon region and the shower region at the muon start plane.  if the
01597 //predictions dont match up, a loop over the digits in the shower 
01598 //region is done to see if the extent of the shower region the 
01599 //planeItr is pre-sliced to be over the planes in the shower region.
01600 void AlgDeMuxBeam::ReconcileShowerAndMuonRegions(DmxPlaneItr &planeItr,
01601                                                                                                  CandDeMuxDigitHandleItr &cdhit)
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 }
01714 
01715 
01716 //---------------------------------------------------------------------
01717 Float_t AlgDeMuxBeam::FindDigitCoG(DmxPlaneItr &validPlaneItr, 
01718                                                                    CandDeMuxDigitHandleItr &cdhit)
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 }
01768 
01769 //------------------------------------------------------------------------------
01770 //function to identify planes whose 3 best hypotheses are all way off from the fit line
01771 //make a cut for them to be at least 12 strips off. only pass this an iterator over valid 
01772 //planes
01773 Bool_t AlgDeMuxBeam::FindPlanesToDropFromFit(DmxPlaneItr &planeItr, Double_t a1, Double_t a2,
01774                                                                                          Double_t a3, Double_t a4,
01775                                                                                          Int_t leverArm)
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 }
01843 
01844 //-----------------------------------------------------------------------
01845 //This method finds the RMS of a possible reconstruction.  It is called in
01846 //the FindEventInDetector method
01847 void AlgDeMuxBeam::FindRMSOfPlanes(DmxPlaneItr &planeItr, Int_t leverArm, 
01848                                                                    Int_t *hypsToUse, Float_t &cog, Float_t &rms)
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 }
02002 
02003 //-----------------------------------------------------------------------
02004 //this method takes the array of demuxPlanes and demuxCoGs and finds 
02005 //the best location for the event in the detector.  the other arguments
02006 //it takes are an iterator over the plane objects, the array to hold the 
02007 //best values, and the number of planes to do.
02008 void AlgDeMuxBeam::DeMuxFirstNPlanes(DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, 
02009                                                                          CandDeMuxDigitHandleItr &cdhit)
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 }
02147 
02148 //----------------------------------------------------------------------
02149 void AlgDeMuxBeam::FindFit(DmxPlaneItr &planeItr, DmxPlaneItr &validPlaneItr, 
02150                                                    Double_t &a1, Double_t &a2, Double_t &a3, 
02151                                                    Double_t &a4, Double_t &bestChi2, Int_t leverArm, 
02152                                                    Int_t *hypotheses)
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 }
02194 
02195 //----------------------------------------------------------------------
02196 //method to find straight line fit over 5 planes
02197 void AlgDeMuxBeam::FindFitForHypothesisSet(DmxPlaneItr &planeItr, Double_t &prevChiSq, 
02198                                                                                    Double_t &a1, Double_t &a2, Double_t &a3, 
02199                                                                                    Double_t &a4, Int_t* hypotheses, 
02200                                                                                    Int_t leverArm)
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 }
02374 
02375
Generated on Mon Feb 15 11:06:17 2010 for loon by  doxygen 1.3.9.1