DmxUtilities Class Reference

#include <DmxUtilities.h>

List of all members.

Public Member Functions

 DmxUtilities ()
virtual ~DmxUtilities ()
void FillPlaneArray (DmxStatus *status, CandDeMuxDigitListHandle &cdlh, AlgConfig &acd)
Float_t CheckForXTalk (Int_t pixel, Float_t *planeSet)
void FillHitPixels (CandDeMuxDigitHandleItr crdhi, Int_t chip, Float_t *eastSet, Float_t *westSet)
Int_t GetDigitPixel (Int_t channel)
Int_t FindVertexPlane (DmxPlaneItr &planeItr)
Int_t FindEndPlane (DmxPlaneItr &planeItr)
Int_t FindMuonStartPlane (DmxPlaneItr &planeItr, Float_t peGainConversion)
Int_t FindBeamMuonStartPlane (DmxPlaneItr &planeItr, Float_t peGainConversion)
Bool_t IsValidFit (DmxPlaneItr &planeItr, Double_t a1, Double_t a2, Double_t a3, Double_t a4, Float_t offset)
Bool_t IsValidFit (DmxPlaneItr &planeItr, Double_t a1, Double_t a2, Double_t a3, Double_t a4)
Int_t FindPlanesOffFit (DmxPlaneItr &planeItr, Int_t strayCut)
Bool_t IsOverlappingMultiple (DmxPlaneItr &planeItr, Float_t vertexZ, Float_t slopeCutOff, Float_t interceptCutOff, Float_t peakCutOff)
Int_t CheckFit (DmxPlaneItr &planeItr)
Float_t CheckFitWithTiming (DmxPlaneItr &planeItr)
void FindCubicFit (Double_t *x, Double_t *y, Double_t *weight, Int_t nPoints, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &a4, Double_t &chiSq)
void FindQuadraticFit (Double_t *x, Double_t *y, Double_t *weight, Int_t nPoints, Double_t &a1, Double_t &a2, Double_t &a3, Double_t &chiSq)
void FindLinearFit (Double_t *x, Double_t *y, Double_t *weight, Int_t nPoints, Double_t &a1, Double_t &a2, Double_t &chiSq)


Detailed Description

Definition at line 25 of file DmxUtilities.h.


Constructor & Destructor Documentation

DmxUtilities::DmxUtilities (  ) 

Definition at line 42 of file DmxUtilities.cxx.

00043 {
00044 }

DmxUtilities::~DmxUtilities (  )  [virtual]

Definition at line 235 of file DmxUtilities.cxx.

00236 {
00237   //MSG("Dmx", Msg::kDebug) << "deleting DmxUtilities object " << endl;
00238 }


Member Function Documentation

Int_t DmxUtilities::CheckFit ( DmxPlaneItr &  planeItr  ) 

Definition at line 764 of file DmxUtilities.cxx.

References DmxPlane::GetCoG(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlane::IsStray(), DmxPlaneTypes::kMuon, DmxPlaneTypes::kShower, and DmxPlane::SetStrips().

Referenced by AlgDeMuxGolden::RunAlg().

00765 {
00766   DmxPlane *plane = 0;
00767 
00768   //compare the best hyps in valid shower planes to the fit hyp for the amount of signal
00769   //mated in the plane.  pick the one with the most mated signal, as long as the difference
00770   //is significant - like 0.2 or more.  also only use the best hyps if they have more than
00771   //85% of the signal mated in them
00772   Float_t firstBestMated = 0.;
00773   Float_t firstSMated = 0.;
00774   Float_t firstTMated = 0.;
00775   Float_t firstMated = 0.;
00776   Float_t lastBestMated = 0.;
00777   Float_t lastSMated = 0.;
00778   Float_t lastTMated = 0.;
00779   Float_t lastMated = 0.;
00780   
00781   Int_t ctr = 0;
00782   Int_t totPlanes = planeItr.SizeSelect();
00783   Int_t fixedPlanes = 0;
00784   planeItr.ResetFirst();
00785 
00786   //make an array to hold the CoG positions of the planes 
00787   Float_t firstPlaneCoG = -1.;
00788   Float_t firstPlaneBCoG = -1.;
00789   Float_t firstPlaneSCoG = -1.;
00790   Float_t firstPlaneTCoG = -1.;
00791   Float_t lastPlaneCoG = -1.;
00792   Float_t lastPlaneBCoG = -1.;
00793   Float_t lastPlaneSCoG = -1.;
00794   Float_t lastPlaneTCoG = -1.;
00795   Float_t firstPlaneSpace = 0.;
00796   Float_t lastPlaneSpace = 0.;
00797 
00798   while(planeItr.IsValid()){
00799     
00800     plane = planeItr.Ptr();
00801     if(ctr == 0){
00802       firstPlaneCoG  = plane->GetCoG();
00803       firstPlaneSpace = plane->GetPlaneNumber();
00804       if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
00805    
00806         firstPlaneBCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
00807         firstBestMated = dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio();
00808       
00809         firstPlaneSCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
00810         firstSMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
00811         
00812         firstPlaneTCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
00813         firstTMated = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
00814       
00815         firstMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio();
00816       }
00817       else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
00818         firstPlaneBCoG = firstPlaneCoG;
00819         firstPlaneSCoG = firstPlaneCoG;
00820         firstPlaneTCoG = firstPlaneCoG;
00821         firstMated = 1.;
00822         firstBestMated = 1.;
00823         firstSMated = 1.;
00824         firstTMated = 1.;
00825       }
00826     }
00827     else if(ctr == 1) firstPlaneSpace -= plane->GetPlaneNumber();
00828     else if(ctr == totPlanes-2) lastPlaneSpace = plane->GetPlaneNumber();
00829     else if(ctr == totPlanes-1){
00830       lastPlaneCoG  = plane->GetCoG();
00831       lastPlaneSpace -= plane->GetPlaneNumber();
00832       if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
00833    
00834         lastPlaneBCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
00835         lastBestMated = dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis()->GetMatedSignalRatio();
00836       
00837         lastPlaneSCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
00838         lastSMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis()->GetMatedSignalRatio();
00839         
00840         lastPlaneTCoG = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
00841         lastTMated = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis()->GetMatedSignalRatio();
00842       
00843         lastMated = dynamic_cast<DmxShowerPlane *>(plane)->GetSetHypothesis()->GetMatedSignalRatio();
00844       }
00845       else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
00846         lastPlaneBCoG = lastPlaneCoG;
00847         lastPlaneSCoG = lastPlaneCoG;
00848         lastPlaneTCoG = lastPlaneCoG;
00849         lastMated = 1.;
00850         lastBestMated = 1.;
00851         lastSMated = 1.;
00852         lastTMated = 1.;
00853       }
00854     }
00855     ++ctr;
00856     planeItr.Next();
00857   }
00858   planeItr.ResetFirst();
00859 
00860   //which is the better hypothesis among the 3 best?
00861   bool useFirstB = false;
00862   bool useFirstS = false;
00863   bool useFirstT = false;
00864   bool useLastB = false;
00865   bool useLastS = false;
00866   bool useLastT = false;
00867 
00868   if(firstBestMated>=firstSMated && firstBestMated>=firstTMated 
00869      && TMath::Abs(firstBestMated-firstMated)>0.2) useFirstB = true;
00870   else if(firstSMated>firstBestMated && firstSMated>=firstTMated 
00871      && TMath::Abs(firstSMated-firstMated)>0.2) useFirstS = true;
00872   else if(firstTMated>firstSMated && firstTMated>firstBestMated 
00873      && TMath::Abs(firstTMated-firstMated)>0.2) useFirstT = true;
00874  
00875   if(lastBestMated>=lastSMated && lastBestMated>=lastTMated 
00876      && TMath::Abs(lastBestMated-lastMated)>0.2) useLastB = true;
00877   else if(lastSMated>lastBestMated && lastSMated>=lastTMated 
00878      && TMath::Abs(lastSMated-lastMated)>0.2) useLastS = true;
00879   else if(lastTMated>lastSMated && lastTMated>lastBestMated 
00880      && TMath::Abs(lastTMated-lastMated)>0.2) useLastT = true;
00881 
00882 
00883 //    cout << "first plane: " << firstPlaneCoG << " " 
00884 //      << firstPlaneBCoG << " " << firstPlaneSCoG 
00885 //         << " " << firstPlaneTCoG << " " << firstMated 
00886 //      << " " << firstBestMated << " " 
00887 //         << firstSMated << " " << firstTMated << " " 
00888 //      << useFirstB << " " << useFirstS << " " 
00889 //         << useFirstT << " " << firstPlaneSpace << endl;
00890 
00891 //    cout << "last plane: " << lastPlaneCoG << " " << lastPlaneBCoG 
00892 //      << " " << lastPlaneSCoG << " " << lastPlaneTCoG << " " << lastMated 
00893 //      << " " << lastBestMated << " " << lastSMated << " " << lastTMated 
00894 //      << " " << useLastB << " " << useLastS << " " << useLastT 
00895 //      << " " << lastPlaneSpace << endl;
00896 
00897   //the first and last planes in the view should be reset to one of the
00898   //best hypothesis
00899   //
00900   //1.  the difference between the best hypothesis and the nearest neighboring
00901   //    valid plane hypothesis is less than 0.25m/plane
00902   //
00903   //    and
00904   //
00905   //2.  the best hypothesis has a lot more signal mated
00906 
00907   //plane iterator was previously set to the first plane in the set
00908     
00909   if(useFirstB 
00910      && TMath::Abs(firstPlaneSpace)*0.25>TMath::Abs(firstPlaneCoG-firstPlaneBCoG)){
00911     plane = planeItr.Ptr();
00912     if(plane->IsStray()){
00913       if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
00914         dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("best");
00915       else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon) plane->SetStrips();
00916       plane->SetStray(false);
00917       ++fixedPlanes;
00918     }
00919   }
00920   else if(useFirstS 
00921           && TMath::Abs(firstPlaneSpace)*0.25>TMath::Abs(firstPlaneCoG-firstPlaneSCoG)){
00922     plane = planeItr.Ptr();
00923     if(plane->IsStray()){
00924       if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
00925         dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("second");
00926       plane->SetStray(false);
00927       ++fixedPlanes;
00928     }
00929   }
00930   else if(useFirstT 
00931           && TMath::Abs(firstPlaneSpace)*0.25>TMath::Abs(firstPlaneCoG-firstPlaneTCoG)){
00932     plane = planeItr.Ptr();
00933     if(plane->IsStray()){
00934       if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
00935         dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("third");
00936       plane->SetStray(false);
00937       ++fixedPlanes;
00938     }
00939   }
00940 
00941   //set the iterator to the last plane
00942   planeItr.ResetLast();
00943 
00944   if(useLastB && TMath::Abs(lastPlaneSpace)*0.25>TMath::Abs(lastPlaneCoG-lastPlaneBCoG)){
00945     plane = planeItr.Ptr();
00946     if(plane->IsStray()){
00947       if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
00948         dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("best");
00949       else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon) plane->SetStrips();
00950       plane->SetStray(false);
00951       ++fixedPlanes;
00952     }
00953   }
00954   else if(useLastS 
00955           && TMath::Abs(lastPlaneSpace)*0.25>TMath::Abs(lastPlaneCoG-lastPlaneSCoG)){
00956     plane = planeItr.Ptr();
00957     if(plane->IsStray()){
00958       if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
00959         dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("second");
00960       plane->SetStray(false);
00961       ++fixedPlanes;
00962     }
00963   }
00964   else if(useLastT 
00965           && TMath::Abs(lastPlaneSpace)*0.25>TMath::Abs(lastPlaneCoG-lastPlaneTCoG)){
00966     plane = planeItr.Ptr();
00967     if(plane->IsStray()){
00968       if(plane->GetPlaneType() == DmxPlaneTypes::kShower) 
00969         dynamic_cast<DmxShowerPlane *>(plane)->SetStrips("third");
00970       plane->SetStray(false);
00971       ++fixedPlanes;
00972     }
00973   }
00974   
00975   return fixedPlanes;
00976 }

Float_t DmxUtilities::CheckFitWithTiming ( DmxPlaneItr &  planeItr  ) 

Definition at line 1160 of file DmxUtilities.cxx.

References DmxPlane::GetCoG(), DmxPlane::GetPlaneView(), and DmxPlane::GetTimingOffset().

Referenced by AlgDeMuxCosmics::RunAlg(), and AlgDeMuxBeam::RunAlg().

01161 {
01162   
01163   //the idea here is to get the offset from the center line from timing for a plane
01164   //in one view and compare that to the center of gravity from the demuxer of the neighboring
01165   //plane in the other view.
01166 
01167   //take the average of the centers of gravity for the planes on either side of the plane you 
01168   //are interested in 
01169 
01170   planeItr.ResetFirst();
01171 
01172   Int_t planeCtr = 0;
01173   Int_t numPlanes = planeItr.SizeSelect();
01174   Int_t forwardCtr = 0;
01175   Int_t backwardCtr = 0;
01176   bool opposite = false;
01177   Float_t offset = 0.;
01178   Float_t avOffsetDiff = 0.;
01179   Float_t avCoG = 0.;
01180   Float_t avCoGCnt = 0.;
01181 
01182   DmxPlane *plane = 0;
01183   DmxPlane *nextPlane = 0;
01184   DmxPlane *prevPlane = 0;
01185 
01186   while( planeItr.IsValid() ){
01187 
01188     plane = planeItr.Ptr();
01189     //MSG("DmxUtil", Msg::kDebug) << "plane of offset = " << plane->GetPlaneNumber() << endl; 
01190 
01191     offset = plane->GetTimingOffset();
01192     //MSG("DmxUtil", Msg::kDebug) << " offset = " << offset << " is valid = " << (Int_t)plane->IsValid() << endl;   
01193     forwardCtr = 0;
01194     backwardCtr = 0;
01195     avCoG = 0.;
01196     avCoGCnt = 0.;
01197     //this is the first plane, so you can only look at the one after it
01198     if(planeCtr == 0 && offset != -10.){
01199 
01200       //advance to the next plane in the opposite view
01201       while(planeItr.IsValid() && opposite == false){
01202         planeItr.Next();
01203         ++forwardCtr;
01204         nextPlane = planeItr.Ptr();
01205         if(nextPlane->GetPlaneView() != plane->GetPlaneView()){
01206           //MSG("DmxUtil", Msg::kDebug) << "next opposite plane to vertex = " << nextPlane->GetPlaneNumber() 
01207           //                 << " cog = " << nextPlane->GetCoG() << endl;
01208           opposite = true;
01209           avOffsetDiff += TMath::Abs(TMath::Abs(offset)-TMath::Abs(nextPlane->GetCoG()));
01210         }
01211       }//end loop to find next plane in opposite view
01212 
01213       while(forwardCtr > 0){
01214         planeItr.Prev();
01215         --forwardCtr;
01216       }//end loop to return planeItr to original plane
01217     
01218     }//end if first plane in set
01219     else if(planeCtr > 0 && planeCtr<numPlanes-1 && offset != -10.){
01220 
01221       opposite = false;
01222       //advance to the next plane in the opposite view
01223       while(planeItr.IsValid() && !opposite){
01224         planeItr.Next();
01225         ++forwardCtr;
01226         if(planeItr.IsValid()){
01227           nextPlane = planeItr.Ptr();
01228           if(nextPlane->GetPlaneView() != plane->GetPlaneView()){
01229             //MSG("DmxUtil", Msg::kDebug) << "next opposite plane = " << nextPlane->GetPlaneNumber() 
01230             //                 << " cog = " << nextPlane->GetCoG() << endl;
01231             opposite = true;
01232             avCoG += nextPlane->GetCoG();
01233             avCoGCnt += 1.;
01234           }
01235         }
01236       }//end loop to find next plane in opposite view
01237 
01238       while(forwardCtr > 0){
01239         planeItr.Prev();
01240         --forwardCtr;
01241       }//end loop to return planeItr to original plane
01242     
01243       opposite = false;
01244       
01245       //if(planeItr.IsValid() )MSG("DmxUtil", Msg::kDebug) << "start backward look from plane " 
01246       //                                                << dynamic_cast<DmxPlane *>(planeItr.Ptr())->GetPlaneNumber()
01247       //                                                << endl;
01248       while(planeItr.IsValid() && !opposite){
01249         planeItr.Prev();
01250         ++backwardCtr;
01251         if(planeItr.IsValid()){
01252           prevPlane = planeItr.Ptr();
01253           if(prevPlane->GetPlaneView() != plane->GetPlaneView()){
01254             //MSG("DmxUtil", Msg::kDebug) << "prev opposite plane = " << prevPlane->GetPlaneNumber() 
01255             //                 << " cog = " << prevPlane->GetCoG() << endl;
01256             opposite = true;
01257             avCoG += prevPlane->GetCoG();
01258             avCoGCnt += 1.;
01259           }
01260         }
01261       }//end loop to find prev plane in opposite view
01262 
01263       while(backwardCtr > 0){
01264         planeItr.Next();
01265         --backwardCtr;
01266       }//end loop to return planeItr to original plane
01267     
01268       //get the average CoG for the planes on either side and subtract the timing offset
01269       if(avCoGCnt > 0.) avOffsetDiff += TMath::Abs(TMath::Abs(offset)-TMath::Abs(avCoG/avCoGCnt));
01270 
01271     }//end if plane between end planes
01272     else if(planeCtr == numPlanes-1 && offset != -10.){
01273       opposite = false;
01274       while(planeItr.IsValid() && !opposite){
01275         planeItr.Prev();
01276         ++backwardCtr;
01277         if(planeItr.IsValid()){
01278           prevPlane = planeItr.Ptr();
01279           if(prevPlane->GetPlaneView() != plane->GetPlaneView()){
01280             //MSG("DmxUtil", Msg::kDebug) << "prev opposite plane to end = " <<prevPlane->GetPlaneNumber() 
01281             //                 << " cog = " << prevPlane->GetCoG() << endl;
01282             opposite = true;
01283             avOffsetDiff += TMath::Abs(TMath::Abs(offset)-TMath::Abs(prevPlane->GetCoG()));
01284           }
01285         }
01286       }//end loop to find prev plane in opposite view
01287 
01288       while(backwardCtr > 0){
01289         planeItr.Next();
01290         --backwardCtr;
01291       }//end loop to return planeItr to original plane
01292 
01293     }//end if plane is end plane
01294 
01295     ++planeCtr;
01296     planeItr.Next();
01297   }//end loop over planes
01298 
01299   if(numPlanes>0) avOffsetDiff /= (1.0*numPlanes);
01300 
01301   planeItr.ResetFirst();
01302 
01303   return avOffsetDiff;
01304 }

Float_t DmxUtilities::CheckForXTalk ( Int_t  pixel,
Float_t *  planeSet 
)

Definition at line 241 of file DmxUtilities.cxx.

Referenced by DmxDeMuxModule::Ana(), DmxDeMuxCosmicsModule::Ana(), and FillPlaneArray().

00242 {
00243   
00244   //go through and look to see if the digit is from cross talk
00245   Float_t crossTalk = -1.;
00246 
00247   if(pixel == 0){
00248     crossTalk = planeSet[1] + planeSet[4];
00249   }
00250   if(pixel == 1) {
00251     crossTalk = planeSet[0] + planeSet[2] + planeSet[5];
00252   }
00253   if(pixel == 2){
00254     crossTalk = planeSet[1] + planeSet[3] + planeSet[6];
00255   }
00256   if(pixel == 3){
00257     crossTalk = planeSet[2] + planeSet[7];
00258   }
00259   if(pixel == 4){
00260     crossTalk = planeSet[0] + planeSet[5] + planeSet[8];
00261   }
00262   if(pixel == 5){
00263     crossTalk = planeSet[1] + planeSet[4] + planeSet[6] + planeSet[9];
00264   }
00265   if(pixel == 6){
00266     crossTalk = planeSet[2] + planeSet[5] + planeSet[7] + planeSet[10];
00267   }
00268   if(pixel == 7){
00269     crossTalk = planeSet[3] + planeSet[6] + planeSet[11];
00270   }
00271   if(pixel == 8){
00272     crossTalk = planeSet[4] + planeSet[9] + planeSet[12];
00273   }
00274   if(pixel == 9){
00275     crossTalk = planeSet[5] + planeSet[8] + planeSet[10] + planeSet[13];
00276   }
00277   if(pixel == 10){
00278     crossTalk = planeSet[6] + planeSet[9] + planeSet[11] + planeSet[14];
00279   }
00280   if(pixel == 11){
00281     crossTalk = planeSet[7] + planeSet[10] + planeSet[15];
00282     //MSG("DmxUtil", Msg::kDebug)<< fUtilities.GetEventNumber() << "\t" 
00283     //               << pixel << "\t" << "\t" << planeSet[11] << "\t" 
00284     //               <<planeSet[7]/planeSet[11] << "\t" 
00285     //               << planeSet[10]/planeSet[11] << "\t" 
00286     //               << planeSet[15]/planeSet[11] << endl; 
00287   }
00288   if(pixel == 12){
00289     crossTalk = planeSet[8] + planeSet[13];
00290   }
00291   if(pixel == 13){
00292     crossTalk = planeSet[9] + planeSet[12] + planeSet[14];
00293   }
00294   if(pixel == 14){
00295     crossTalk = planeSet[10] + planeSet[13] + planeSet[15];
00296     //MSG("DmxUtil", Msg::kDebug)<< fUtilities.GetEventNumber() << "\t" 
00297     //               << pixel << "\t" << "\t" << planeSet[14] << "\t" 
00298     //               <<planeSet[13]/planeSet[14] << "\t" 
00299     //               << planeSet[10]/planeSet[14] << "\t" 
00300     //               << planeSet[15]/planeSet[14] << endl; 
00301   }
00302   if(pixel == 15){
00303     crossTalk = planeSet[11] + planeSet[14];
00304   }
00305   
00306   return crossTalk;
00307 }

void DmxUtilities::FillHitPixels ( CandDeMuxDigitHandleItr  crdhi,
Int_t  chip,
Float_t *  eastSet,
Float_t *  westSet 
)

Definition at line 333 of file DmxUtilities.cxx.

References digit(), GetDigitPixel(), StripEnd::kEast, and StripEnd::kWest.

Referenced by DmxDeMuxModule::Ana(), DmxDeMuxCosmicsModule::Ana(), and FillPlaneArray().

00335 {
00336   Int_t pixel = -1;
00337   //loop over the digits and see if they have been set
00338   while( crdhi.IsValid() ){
00339             
00340     //get the digit and it's strip, plane, and side
00341     CandDeMuxDigitHandle *digit = crdhi.Ptr();
00342     Int_t digitChip = digit->GetChannelId().GetVaChip();
00343     Int_t channel = digit->GetChannelId().GetVaChannel();
00344     pixel = GetDigitPixel(channel);
00345     Float_t charge = digit->GetCharge();
00346 
00347     if( digitChip == chip){
00348       if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast){
00349         eastSet[pixel] += charge;
00350       }
00351       if(digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest){
00352         westSet[pixel] += charge;
00353       }
00354     }
00355     //MSG("Dmx2", Msg::kDebug)<< pixel << "\t" << channel<< "\t" << chip << "\t" 
00356     //             << planeSet[pixel] << endl;
00357    
00358     crdhi.Next();
00359   }
00360   
00361   crdhi.ResetFirst();
00362   return;
00363 }

void DmxUtilities::FillPlaneArray ( DmxStatus status,
CandDeMuxDigitListHandle cdlh,
AlgConfig acd 
)

Definition at line 47 of file DmxUtilities.cxx.

References CheckForXTalk(), digit(), FillHitPixels(), CandHandle::GetDaughterIterator(), GetDigitPixel(), Registry::GetDouble(), Registry::GetInt(), CandHandle::GetNDaughters(), PlexSEIdAltL::GetPlane(), CandDigitHandle::GetPlexSEIdAltL(), CandDigitHandle::GetPlexSEIdAltLWritable(), PlexSEIdAltL::IsVetoShield(), Msg::kDebug, StripEnd::kEast, CalHelpers::KeyFromPlane(), CandDeMuxDigit::kLowSignal, StripEnd::kWest, CandDeMuxDigit::kXTalk, MSG, PlexSEIdAltL::SetDemuxVetoFlag(), and DmxStatus::SetPlaneArray().

Referenced by AlgDeMuxGolden::RunAlg(), AlgDeMuxCosmics::RunAlg(), and AlgDeMuxBeam::RunAlg().

00049 {
00050                 
00051         TObjArray *planeArray = new TObjArray();
00052         MSG("DmxUtil", Msg::kDebug) << cdlh.GetNDaughters() << " digits in DmxUtil" << endl;
00053 
00054         //make two iterators - one over the whole set and one to iterate over each plane
00055         //plane in turn
00056         CandDeMuxDigitHandleItr fCdhit = cdlh.GetDaughterIterator();
00057         
00058         //create the KeyFunc
00059         CandDeMuxDigitHandleKeyFunc *planeKF = fCdhit.CreateKeyFunc();
00060         
00061         //program the KeyFunc with the sort function
00062         planeKF->SetFun(KeyFromPlane);
00063         
00064         //get the NavSet from the iterator and pass the KeyFunc to it
00065         fCdhit.GetSet()->AdoptSortKeyFunc(planeKF);
00066         
00067         //clear the KF pointer because we no longer own the KeyFunc
00068         planeKF = 0;
00069 
00070         //initialize variable to keep track of the current plane number
00071         Int_t planeNum = 0;
00072 
00073         //MSG("DmxUtil", Msg::kDebug) << "current plane = " << planeNum << endl;
00074 
00075         CandDeMuxDigitHandle *currentDigit = 0;
00076         CandDeMuxDigitHandle *digit = 0;
00077 
00078         //go through the sorted Digits
00079         while( (currentDigit = fCdhit()) ){
00080 
00081                 Int_t plane = 0;
00082 
00083                 //find the plane number corresponding to that Digit
00084                 plane = currentDigit->GetPlexSEIdAltL().GetPlane();
00085     
00086                 //MSG("DmxUtil", Msg::kDebug) << "plane number of current digit = " << plane 
00087                 //                  << endl;
00088     
00089                 //if the plane number of the current digit is different from the
00090                 //current plane number, then add another DmxPlaneDebug object to 
00091                 //the TObjArray container for that plane number and increment the
00092                 //fNumberOfPlanes variable
00093     
00094                 if(plane != planeNum && !currentDigit->GetPlexSEIdAltL().IsVetoShield()){
00095       
00096                         planeNum = plane;
00097 
00098                         //make a new iterator for this plane's digit.  dont move this outside the
00099                         //loop because then you have one iterator for all planes versus 1 iterator
00100                         //per plane.  makes a big difference.
00101                         CandDeMuxDigitHandleItr planeCdhit(cdlh.GetDaughterIterator());
00102                         CandDeMuxDigitHandleKeyFunc *plane1KF = planeCdhit.CreateKeyFunc();
00103                         plane1KF->SetFun(KeyFromPlane);
00104                         planeCdhit.GetSet()->AdoptSortKeyFunc(plane1KF,true,false);
00105                         plane1KF = 0;
00106 
00107                         //get the slice for this plane
00108                         planeCdhit.GetSet()->Slice(planeNum);
00109                         planeCdhit.GetSet()->Update();
00110                         
00111                         Int_t numDigits = 0;
00112                         if( acd.GetInt("UseCandDigitMask")==1 ){
00113                                 //loop over the digits in this plane and figure out which ones might be 
00114                                 //from cross talk.  set the cutoff value for signal fraction to be in 
00115                                 //AlgConfigDeMux
00116                                 Float_t eastSet0[] = {0., 0., 0., 0., 0., 0., 0., 0., 
00117                                                                           0., 0., 0., 0., 0., 0., 0., 0.};
00118                                 Float_t westSet0[] = {0., 0., 0., 0., 0., 0., 0., 0., 
00119                                                   0., 0., 0., 0., 0., 0., 0., 0.};
00120                                 Float_t eastSet1[] = {0., 0., 0., 0., 0., 0., 0., 0., 
00121                                                                           0., 0., 0., 0., 0., 0., 0., 0.};
00122                                 Float_t westSet1[] = {0., 0., 0., 0., 0., 0., 0., 0., 
00123                                                                           0., 0., 0., 0., 0., 0., 0., 0.};
00124                                 Float_t eastSet2[] = {0., 0., 0., 0., 0., 0., 0., 0., 
00125                                                                           0., 0., 0., 0., 0., 0., 0., 0.};
00126                                 Float_t westSet2[] = {0., 0., 0., 0., 0., 0., 0., 0., 
00127                                                                           0., 0., 0., 0., 0., 0., 0., 0.};
00128         
00129                                 FillHitPixels(planeCdhit, 0, eastSet0, westSet0);
00130                                 FillHitPixels(planeCdhit, 1, eastSet1, westSet1);
00131                                 FillHitPixels(planeCdhit, 2, eastSet2, westSet2);
00132                                 Int_t pixel = -1;
00133                                 Float_t crossTalk = 0.; 
00134 
00135                                 numDigits = planeCdhit.SizeSelect();
00136 
00137                                 //loop over the digits in this plane to identify low signal and xtalk
00138                                 //digits.  also set the DeMuxVetoFlag for all digits in this plane - if the
00139                                 //digits have a strip end in the demuxed hypothesis, the DeMuxVetoFlag will
00140                                 //be unset when the digit is demuxed.
00141                                 while( (digit = planeCdhit()) ){
00142 
00143                                         pixel =  GetDigitPixel(digit->GetChannelId().GetVaChannel());
00144                                         crossTalk = 0.;
00145                 
00146                                         currentDigit->GetPlexSEIdAltLWritable().SetDemuxVetoFlag(1);
00147           
00148                                         if(digit->GetCharge()
00149                                            <acd.GetDouble("AveragePEGainConversion")*acd.GetDouble("IgnoreSignalLimit")){
00150                                                 digit->SetDeMuxDigitFlagBit(CandDeMuxDigit::kLowSignal);
00151                                                 //MSG("DmxUtil", Msg::kDebug) << "digit with charge " << digit->GetCharge() 
00152                                                 //                     << " has too little charge" << endl;
00153                                                 --numDigits;
00154                                         }
00155             
00156                                         if( digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kEast 
00157                                                 && digit->GetDeMuxDigitFlagWord() == 0){
00158 
00159                                                 if(digit->GetChannelId().GetVaChip() == 0)
00160                                                         crossTalk = CheckForXTalk(pixel, eastSet0);
00161                                                 else if(digit->GetChannelId().GetVaChip() == 1)
00162                                                         crossTalk = CheckForXTalk(pixel, eastSet1);
00163                                                 else if(digit->GetChannelId().GetVaChip() == 2)
00164                                                         crossTalk = CheckForXTalk(pixel, eastSet2);
00165 
00166                                                 //mark it if its cross talk and subtract its charge to the total for the side
00167                                                 if(digit->GetCharge()<(acd.GetDouble("XTalkFractionLimit")*crossTalk)){
00168                                                         digit->SetDeMuxDigitFlagBit(CandDeMuxDigit::kXTalk);
00169                                                         --numDigits;
00170                                                         //MSG("DmxUtil", Msg::kDebug) << "digit with charge = " << digit->GetCharge() 
00171                                                         //                << " is xtalk" << endl;
00172                                                 }
00173                                         }//end if digit is from the east side and not already flagged
00174                                         else if( digit->GetPlexSEIdAltL().GetEnd() == StripEnd::kWest &&
00175                                                          digit->GetDeMuxDigitFlagWord() == 0){
00176 
00177                                                 if(digit->GetChannelId().GetVaChip() == 0)
00178                                                         crossTalk = CheckForXTalk(pixel, westSet0);
00179                                                 else if(digit->GetChannelId().GetVaChip() == 1)
00180                                                         crossTalk = CheckForXTalk(pixel, westSet1);
00181                                                 else if(digit->GetChannelId().GetVaChip() == 2)
00182                                                         crossTalk = CheckForXTalk(pixel, westSet2);
00183           
00184                                                  //mark it if its not cross talk and add its charge to the total for the side
00185                                                 if(digit->GetCharge()<(acd.GetDouble("XTalkFractionLimit")*crossTalk)){
00186                                                         digit->SetDeMuxDigitFlagBit(CandDeMuxDigit::kXTalk);
00187                                                         --numDigits;
00188                                                         //MSG("DmxUtil", Msg::kDebug) << "digit with charge = " << digit->GetCharge() 
00189                                                         //                << " is xtalk" << endl;
00190                                                 }
00191                                         }//end if digit is from the west side and not already marked
00192           
00193                                 }//end loop over digits to mark xtalk digits
00194                         }
00195                         //if not marking xtalk digits, 
00196                         //then set numDigits to size of the whole set
00197                         else numDigits = planeCdhit.SizeSelect(); 
00198 
00199                         planeCdhit.ResetFirst();
00200       
00201                         //determine the type of plane object by the number of digits it
00202                         //has.  if > 2, its a shower plane, if  <= 2 its 
00203                         //a muon plane.  dont do a thing if the plane is a veto shield plane
00204      
00205                         if(numDigits > 2 && planeNum>-1){
00206                                 planeArray->AddLast(new DmxShowerPlane(acd, planeCdhit, planeNum));
00207                                 //MSG("DmxUtil", Msg::kDebug) << "Plane " << planeNum << " is a Shower Plane" 
00208                                 //      << " with " << numDigits << " digits" << endl;
00209                         }
00210                         else if(numDigits<=2 && planeNum>-1){ 
00211                                 planeArray->AddLast(new DmxMuonPlane(acd, planeCdhit, planeNum));
00212                                 //MSG("DmxUtil", Msg::kDebug) << "Plane " << planeNum << " is a Muon Plane" 
00213                                 //      << " with " << numDigits << " digits" << endl;
00214                         }
00215 
00216                         //reset the iterator over each plane's digits
00217                         planeCdhit.GetSet()->ClearSlice(false);
00218                         
00219                 }//end if on a new plane
00220 
00221         }//end loop over all digits
00222 
00223         //cout << planeArray->GetLast() << endl;
00224         status->SetPlaneArray(planeArray);
00225 
00226         //zero the pointer to the planeArray because the DmxStatus object now owns the 
00227         //array
00228         planeArray = 0;
00229 
00230         //MSG("DmxUtil", Msg::kDebug) << "done filling plane array" << endl;
00231         return;
00232 }

Int_t DmxUtilities::FindBeamMuonStartPlane ( DmxPlaneItr &  planeItr,
Float_t  peGainConversion 
)

Definition at line 575 of file DmxUtilities.cxx.

References DmxPlane::GetNumberOfStrips(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlaneTypes::GetPlaneTypeString(), DmxPlane::IsValid(), Msg::kDebug, DmxPlaneTypes::kShower, and MSG.

Referenced by AlgDeMuxBeam::RunAlg().

00576 {
00577         //now locate where the muon starts, if anywhere - define the start
00578         //plane to be the first muon plane in a set of 4/5 muon planes with 
00579         //3 of the 4 having at least 3 pe's each
00580         
00581         //this method assumes that the planeItr was sliced to be from the vertex to the 
00582         //end planes of the event
00583         
00584         Int_t fMuonStartPlaneNumber = -1;
00585         
00586         Int_t ctr = 0;
00587         Int_t muonCtr = 0;
00588         Int_t startPlane = 0;
00589         Int_t endPlane = 0;
00590         Int_t peCtr = 0;
00591         Int_t lastEventPlane = 0;
00592         
00593         DmxPlane *currentPlane = 0;
00594         
00595         while( (currentPlane = planeItr()) ){
00596                 
00597                 MSG("DmxUtilities", Msg::kDebug) << "plane " << currentPlane->GetPlaneNumber()
00598                 << " type " << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType())
00599                 << " charge " << currentPlane->GetPlaneCharge()
00600                 << " strips " << currentPlane->GetNumberOfStrips()
00601                 << " is valid " << (Int_t)currentPlane->IsValid() << endl;
00602                 
00603                 lastEventPlane = currentPlane->GetPlaneNumber();
00604         }
00605         
00606         planeItr.ResetFirst();
00607         
00608         while(planeItr.IsValid() && fMuonStartPlaneNumber == -1){
00609                 ctr = 0;
00610                 
00611                 muonCtr = 0;
00612                 startPlane = planeItr.Ptr()->GetPlaneNumber();
00613                 endPlane = startPlane+10;
00614                 peCtr = 0;
00615                 
00616                 MSG("DmxUtilities", Msg::kDebug) << "plane number " << startPlane 
00617                         << " charge " << planeItr.Ptr()->GetPlaneCharge()
00618                         << " type " 
00619                         << DmxPlaneTypes::GetPlaneTypeString(planeItr.Ptr()->GetPlaneType()) 
00620                         << endl;
00621                 planeItr.Next();
00622                 
00623                 //check 4 planes at a time
00624                 while(ctr < 4 && planeItr.IsValid() ){
00625                         currentPlane = planeItr.Ptr();
00626                         
00627                         MSG("DmxUtilities", Msg::kDebug) << "plane number " << currentPlane->GetPlaneNumber()
00628                                 << " charge " << currentPlane->GetPlaneCharge()
00629                                 << " type " 
00630                                 << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType()) 
00631                                 << " num strips " << currentPlane->GetNumberOfStrips()
00632                                 << endl;
00633                         endPlane = currentPlane->GetPlaneNumber();
00634                         
00635                         //check to see that its a muon plane or only has 2 strips hit
00636                         if( currentPlane->GetNumberOfStrips() < 3 )     ++muonCtr;
00637                         
00638                         if(currentPlane->GetPlaneCharge() > (peGainConversion * 3.)) ++peCtr;
00639                         
00640                         MSG("DmxUtilities", Msg::kDebug) << " muon planes = " << muonCtr 
00641                                 << " peCtr = " << peCtr << endl;
00642                         planeItr.Next();
00643                         ++ctr;
00644                 }
00645                 
00646                 if(muonCtr == 4 && peCtr >= 3 && endPlane-startPlane<=5) 
00647                         fMuonStartPlaneNumber = startPlane;
00648                 
00649                 //if there are still planes left to examine, the iterator will be valid
00650                 //otherwise, you ran out of planes, so break the cycle
00651                 if(planeItr.IsValid() ){
00652                         planeItr.Prev();
00653                         planeItr.Prev();
00654                         planeItr.Prev();
00655                         planeItr.Prev();
00656                 }
00657         }//end loop to find the muon track starting plane
00658         
00659         //found the most likely start plane for the muon.  but lets make sure 
00660         //we arent being fooled into thinking a shower event is really a muon event.  
00661         //loop over the planes and see the number of planes that are shower like and if 
00662         //there are lots more shower planes than muon planes
00663         planeItr.GetSet()->ClearSlice(false);
00664         
00665         planeItr.GetSet()->Slice(fMuonStartPlaneNumber, lastEventPlane);
00666         planeItr.GetSet()->Update();
00667         planeItr.ResetFirst();
00668         
00669         Int_t numShowerPlanes = 0;
00670         Int_t numMuonPlanes = 0;
00671         while( (currentPlane = planeItr()) ){
00672                 
00673                 if(currentPlane->GetPlaneType() == DmxPlaneTypes::kShower
00674                    && currentPlane->GetNumberOfStrips()>2) ++numShowerPlanes;
00675                 else ++numMuonPlanes;
00676                 
00677         }
00678         
00679         if(1.*numShowerPlanes>0.25*numMuonPlanes) fMuonStartPlaneNumber = -1;
00680 
00681         planeItr.GetSet()->ClearSlice(false);
00682         
00683         return fMuonStartPlaneNumber;
00684 }

void DmxUtilities::FindCubicFit ( Double_t *  x,
Double_t *  y,
Double_t *  weight,
Int_t  nPoints,
Double_t &  a1,
Double_t &  a2,
Double_t &  a3,
Double_t &  a4,
Double_t &  chiSq 
)

Definition at line 1309 of file DmxUtilities.cxx.

01312 {
01313 
01314   //use method of determinants to do the least squares fit
01315   //this is from bevington chapter 7
01316   //
01317   //f_1(x) = 1
01318   //f_2(x) = x
01319   //f_3(x) = x^2
01320   //f_4(x) = x^3
01321   
01322   //need to find the sums
01323   //SIGMA y_i*f_1(x_i)/weight_i
01324   //SIGMA y_i*f_2(x_i)/weight_i
01325   //SIGMA y_i*f_3(x_i)/weight_i
01326   //SIGMA y_i*f_4(x_i)/weight_i
01327   //SIGMA f_1(x_i)*f_1(x_i)/weight_i
01328   //SIGMA f_1(x_i)*f_2(x_i)/weight_i
01329   //SIGMA f_1(x_i)*f_3(x_i)/weight_i
01330   //SIGMA f_1(x_i)*f_4(x_i)/weight_i
01331   //SIGMA f_2(x_i)*f_1(x_i)/weight_i
01332   //SIGMA f_2(x_i)*f_2(x_i)/weight_i
01333   //SIGMA f_2(x_i)*f_3(x_i)/weight_i
01334   //SIGMA f_2(x_i)*f_4(x_i)/weight_i
01335   //SIGMA f_3(x_i)*f_1(x_i)/weight_i
01336   //SIGMA f_3(x_i)*f_2(x_i)/weight_i
01337   //SIGMA f_3(x_i)*f_3(x_i)/weight_i
01338   //SIGMA f_3(x_i)*f_4(x_i)/weight_i
01339   //SIGMA f_4(x_i)*f_1(x_i)/weight_i
01340   //SIGMA f_4(x_i)*f_2(x_i)/weight_i
01341   //SIGMA f_4(x_i)*f_3(x_i)/weight_i
01342   //SIGMA f_4(x_i)*f_4(x_i)/weight_i
01343   //
01344   //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need
01345   //10 instead of 16 variables
01346 
01347   Double_t yf1divWeight = 0.;
01348   Double_t yf2divWeight = 0.;
01349   Double_t yf3divWeight = 0.;
01350   Double_t yf4divWeight = 0.;
01351   Double_t f1f1divWeight = 0.;
01352   Double_t f1f2divWeight = 0.;
01353   Double_t f1f3divWeight = 0.;
01354   Double_t f1f4divWeight = 0.;
01355   Double_t f2f2divWeight = 0.;
01356   Double_t f2f3divWeight = 0.;
01357   Double_t f2f4divWeight = 0.;
01358   Double_t f3f3divWeight = 0.;
01359   Double_t f3f4divWeight = 0.;
01360   Double_t f4f4divWeight = 0.;
01361   
01362   for(Int_t i=0; i<nPoints; i++){
01363     yf1divWeight += y[i]/(weight[i]*weight[i]);
01364     yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]);
01365     yf3divWeight += (y[i]*x[i]*x[i])/(weight[i]*weight[i]);
01366     yf4divWeight += (y[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01367     f1f1divWeight += 1./(weight[i]*weight[i]);
01368     f1f2divWeight += x[i]/(weight[i]*weight[i]);
01369     f1f3divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
01370     f1f4divWeight += (x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01371     f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
01372     f2f3divWeight += (x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01373     f2f4divWeight += (x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01374     f3f3divWeight += (x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01375     f3f4divWeight += (x[i]*x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01376     f4f4divWeight += (x[i]*x[i]*x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01377   }
01378 
01379   //make a bunch of matricies to hold these values
01380   TMatrix deltaMatrix(4,4);
01381   TMatrix a1Matrix(4,4);
01382   TMatrix a2Matrix(4,4);
01383   TMatrix a3Matrix(4,4);
01384   TMatrix a4Matrix(4,4);
01385   
01386   //fill the matricies
01387   deltaMatrix(0,0) = f1f1divWeight;
01388   deltaMatrix(0,1) = f1f2divWeight;
01389   deltaMatrix(0,2) = f1f3divWeight;
01390   deltaMatrix(0,3) = f1f4divWeight;
01391   deltaMatrix(1,0) = f1f2divWeight;
01392   deltaMatrix(1,1) = f2f2divWeight;
01393   deltaMatrix(1,2) = f2f3divWeight;
01394   deltaMatrix(1,3) = f2f4divWeight;
01395   deltaMatrix(2,0) = f1f3divWeight;
01396   deltaMatrix(2,1) = f2f3divWeight;
01397   deltaMatrix(2,2) = f3f3divWeight;
01398   deltaMatrix(2,3) = f3f4divWeight;
01399   deltaMatrix(3,0) = f1f4divWeight;
01400   deltaMatrix(3,1) = f2f4divWeight;
01401   deltaMatrix(3,2) = f3f4divWeight;
01402   deltaMatrix(3,3) = f4f4divWeight;
01403   
01404   a1Matrix(0,0) = yf1divWeight;
01405   a1Matrix(0,1) = f1f2divWeight;
01406   a1Matrix(0,2) = f1f3divWeight;
01407   a1Matrix(0,3) = f1f4divWeight;
01408   a1Matrix(1,0) = yf2divWeight;
01409   a1Matrix(1,1) = f2f2divWeight;
01410   a1Matrix(1,2) = f2f3divWeight;
01411   a1Matrix(1,3) = f2f4divWeight;
01412   a1Matrix(2,0) = yf3divWeight;
01413   a1Matrix(2,1) = f2f3divWeight;
01414   a1Matrix(2,2) = f3f3divWeight;
01415   a1Matrix(2,3) = f3f4divWeight;
01416   a1Matrix(3,0) = yf4divWeight;
01417   a1Matrix(3,1) = f2f4divWeight;
01418   a1Matrix(3,2) = f3f4divWeight;
01419   a1Matrix(3,3) = f4f4divWeight;
01420   
01421   a2Matrix(0,0) = f1f1divWeight;
01422   a2Matrix(0,1) = yf1divWeight;
01423   a2Matrix(0,2) = f1f3divWeight;
01424   a2Matrix(0,3) = f1f4divWeight;
01425   a2Matrix(1,0) = f1f2divWeight;
01426   a2Matrix(1,1) = yf2divWeight;
01427   a2Matrix(1,2) = f2f3divWeight;
01428   a2Matrix(1,3) = f2f4divWeight;
01429   a2Matrix(2,0) = f1f3divWeight;
01430   a2Matrix(2,1) = yf3divWeight;
01431   a2Matrix(2,2) = f3f3divWeight;
01432   a2Matrix(2,3) = f3f4divWeight;
01433   a2Matrix(3,0) = f1f4divWeight;
01434   a2Matrix(3,1) = yf4divWeight;
01435   a2Matrix(3,2) = f3f4divWeight;
01436   a2Matrix(3,3) = f4f4divWeight;
01437   
01438   a3Matrix(0,0) = f1f1divWeight;
01439   a3Matrix(0,1) = f1f2divWeight;
01440   a3Matrix(0,2) = yf1divWeight;
01441   a3Matrix(0,3) = f1f4divWeight;
01442   a3Matrix(1,0) = f1f2divWeight;
01443   a3Matrix(1,1) = f2f2divWeight;
01444   a3Matrix(1,2) = yf2divWeight;
01445   a3Matrix(1,3) = f2f4divWeight;
01446   a3Matrix(2,0) = f1f3divWeight;
01447   a3Matrix(2,1) = f2f3divWeight;
01448   a3Matrix(2,2) = yf3divWeight;
01449   a3Matrix(2,3) = f3f4divWeight;
01450   a3Matrix(3,0) = f1f4divWeight;
01451   a3Matrix(3,1) = f2f4divWeight;
01452   a3Matrix(3,2) = yf4divWeight;
01453   a3Matrix(3,3) = f4f4divWeight;
01454   
01455   a4Matrix(0,0) = f1f1divWeight;
01456   a4Matrix(0,1) = f1f2divWeight;
01457   a4Matrix(0,2) = f1f3divWeight;
01458   a4Matrix(0,3) = yf1divWeight;
01459   a4Matrix(1,0) = f1f2divWeight;
01460   a4Matrix(1,1) = f2f2divWeight;
01461   a4Matrix(1,2) = f2f3divWeight;
01462   a4Matrix(1,3) = yf2divWeight;
01463   a4Matrix(2,0) = f1f3divWeight;
01464   a4Matrix(2,1) = f2f3divWeight;
01465   a4Matrix(2,2) = f3f3divWeight;
01466   a4Matrix(2,3) = yf3divWeight;
01467   a4Matrix(3,0) = f1f4divWeight;
01468   a4Matrix(3,1) = f2f4divWeight;
01469   a4Matrix(3,2) = f3f4divWeight;
01470   a4Matrix(3,3) = yf4divWeight;
01471   
01472   a1 = -10000.;
01473   a2 = -10000.;
01474   a3 = -10000.;
01475   a4 = -10000.;
01476   Double_t delta = deltaMatrix.Determinant();
01477   
01478   if(delta != 0.){
01479     a1 = a1Matrix.Determinant()/delta;
01480     a2 = a2Matrix.Determinant()/delta;
01481     a3 = a3Matrix.Determinant()/delta;
01482     a4 = a4Matrix.Determinant()/delta;
01483   }
01484 
01485   //find the chi^2 value for the fit
01486   chiSq = 0.;
01487   Double_t chi = 0.;
01488   for(Int_t j = 0; j<nPoints; j++){
01489     
01490     chi = (y[j] - (a1 + a2*x[j] + a3*(x[j]*x[j]) + a4*(x[j]*x[j]*x[j])))/weight[j];
01491     chiSq += chi*chi;
01492   
01493   }  
01494   
01495   //get the chiSq per degree of freedom
01496 
01497   chiSq /= (1.*nPoints);
01498 
01499   return;
01500 }

Int_t DmxUtilities::FindEndPlane ( DmxPlaneItr &  planeItr  ) 

Definition at line 463 of file DmxUtilities.cxx.

Referenced by AlgDeMuxGolden::RunAlg(), AlgDeMuxCosmics::RunAlg(), and AlgDeMuxBeam::RunAlg().

00464 {
00465   Int_t endPlaneNumber = 500;
00466 
00467   while( planeItr.IsValid() && endPlaneNumber == 500){
00468     Int_t planeNumber = planeItr.Ptr()->GetPlaneNumber();
00469     
00470     planeItr.Next();
00471     if( planeItr.IsValid() ){
00472       //if the difference between the current plane number and the previous is >= 10
00473       //call the previous plane the last one in the event
00474       if( TMath::Abs(planeItr.Ptr()->GetPlaneNumber() 
00475                      - planeNumber) >= 5 ){ endPlaneNumber = planeNumber; }
00476     }
00477     else{
00478       //didn't find a cutoff plane so call the last plane it
00479       endPlaneNumber = planeNumber;
00480     }
00481   }
00482 
00483   return endPlaneNumber;
00484 }

void DmxUtilities::FindLinearFit ( Double_t *  x,
Double_t *  y,
Double_t *  weight,
Int_t  nPoints,
Double_t &  a1,
Double_t &  a2,
Double_t &  chiSq 
)

Definition at line 1634 of file DmxUtilities.cxx.

References Msg::kDebug, and MSG.

Referenced by AlgDeMuxCosmics::FindFit(), AlgDeMuxBeam::FindFitForHypothesisSet(), and AlgDeMuxCosmics::FindWindowCosmicSolution().

01637 {
01638 
01639         //use method of determinants to do the least squares fit
01640         //this is from bevington chapter 6
01641         //
01642         //f_1(x) = 1
01643         //f_2(x) = x
01644   
01645         //need to find the sums
01646         //SIGMA y_i*f_1(x_i)/weight_i
01647         //SIGMA y_i*f_2(x_i)/weight_i
01648         //SIGMA f_1(x_i)*f_1(x_i)/weight_i
01649         //SIGMA f_1(x_i)*f_2(x_i)/weight_i
01650         //SIGMA f_2(x_i)*f_1(x_i)/weight_i
01651         //SIGMA f_2(x_i)*f_2(x_i)/weight_i
01652         //
01653         //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need
01654         //3 instead of 4 variables
01655 
01656         Double_t yf1divWeight = 0.;
01657         Double_t yf2divWeight = 0.;
01658         Double_t f1f1divWeight = 0.;
01659         Double_t f1f2divWeight = 0.;
01660         Double_t f2f2divWeight = 0.;
01661   
01662         for(Int_t i=0; i<nPoints; ++i){
01663                 MSG("DmxUtilities", Msg::kDebug) << " x = " << x[i]
01664                                                                                 << " y = " << y[i] 
01665                                                                                 << " weight = " << weight[i] << endl;
01666                 if(weight[i]>0.){
01667                         yf1divWeight += y[i]/(weight[i]*weight[i]);
01668                         yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]);
01669                         f1f1divWeight += 1./(weight[i]*weight[i]);
01670                         f1f2divWeight += x[i]/(weight[i]*weight[i]);
01671                         f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]);                     
01672                 }
01673         }
01674 
01675         //make a bunch of matricies to hold these values
01676         TMatrix deltaMatrix(2,2);
01677         TMatrix a1Matrix(2,2);
01678         TMatrix a2Matrix(2,2);
01679   
01680         //fill the matricies
01681         deltaMatrix(0,0) = f1f1divWeight;
01682         deltaMatrix(0,1) = f1f2divWeight;
01683         deltaMatrix(1,0) = f1f2divWeight;
01684         deltaMatrix(1,1) = f2f2divWeight;
01685   
01686         a1Matrix(0,0) = yf1divWeight;
01687         a1Matrix(0,1) = f1f2divWeight;
01688         a1Matrix(1,0) = yf2divWeight;
01689         a1Matrix(1,1) = f2f2divWeight;
01690   
01691         a2Matrix(0,0) = f1f1divWeight;
01692         a2Matrix(0,1) = yf1divWeight;
01693         a2Matrix(1,0) = f1f2divWeight;
01694         a2Matrix(1,1) = yf2divWeight;
01695   
01696   
01697         a1 = -10000.;
01698         a2 = -10000.;
01699         Double_t delta = deltaMatrix.Determinant();
01700   
01701         if(delta != 0.){
01702                 a1 = a1Matrix.Determinant()/delta;
01703                 a2 = a2Matrix.Determinant()/delta;
01704         }
01705 
01706         //find the chi^2 value for the fit
01707         chiSq = 0.;
01708         Double_t chi = 0.;
01709         for(Int_t j = 0; j<nPoints; ++j){
01710                 
01711                 chi = (y[j] - (a1 + a2*x[j]))/weight[j];
01712                 chiSq += chi*chi;
01713   
01714         }  
01715   
01716         //get the chiSq per degree of freedom
01717 
01718         chiSq /= (1.*nPoints);
01719 
01720         return;
01721 }

Int_t DmxUtilities::FindMuonStartPlane ( DmxPlaneItr &  planeItr,
Float_t  peGainConversion 
)

Definition at line 487 of file DmxUtilities.cxx.

References DmxPlane::GetNumberOfStrips(), DmxPlane::GetPlaneCharge(), DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlaneTypes::GetPlaneTypeString(), DmxPlane::IsValid(), Msg::kDebug, and MSG.

00488 {
00489         //now locate where the muon starts, if anywhere - define the start
00490         //plane to be the first muon plane in a set of 4/5 muon planes with 
00491         //3 of the 4 having at least 3 pe's each
00492   
00493         //this method assumes that the planeItr was sliced to be from the vertex to the 
00494         //end planes of the event
00495 
00496         Int_t fMuonStartPlaneNumber = -1;
00497 
00498         Int_t ctr = 0;
00499         Int_t muonCtr = 0;
00500         Int_t startPlane = 0;
00501         Int_t endPlane = 0;
00502         Int_t peCtr = 0;
00503         
00504         DmxPlane *currentPlane = 0;
00505         
00506         while( (currentPlane = planeItr()) ){
00507                 
00508                 MSG("DmxUtilities", Msg::kDebug) << "plane " << currentPlane->GetPlaneNumber()
00509                                                                                 << " type " << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType())
00510                                                                                 << " charge " << currentPlane->GetPlaneCharge()
00511                                                                                 << " strips " << currentPlane->GetNumberOfStrips()
00512                                                                                 << " is valid " << (Int_t)currentPlane->IsValid() << endl;
00513         }
00514 
00515         planeItr.ResetFirst();
00516         
00517         while(planeItr.IsValid() && fMuonStartPlaneNumber == -1){
00518                 ctr = 0;
00519 
00520                 muonCtr = 0;
00521                 startPlane = planeItr.Ptr()->GetPlaneNumber();
00522                 endPlane = startPlane+10;
00523                 peCtr = 0;
00524     
00525                 MSG("DmxUtilities", Msg::kDebug) << "plane number " << startPlane 
00526                                                                                 << " charge " << planeItr.Ptr()->GetPlaneCharge()
00527                                                                                 << " type " 
00528                                                                             << DmxPlaneTypes::GetPlaneTypeString(planeItr.Ptr()->GetPlaneType()) 
00529                                                                                 << endl;
00530                 planeItr.Next();
00531                         
00532                 //check 4 planes at a time
00533                 while(ctr < 4 && planeItr.IsValid() ){
00534                         currentPlane = planeItr.Ptr();
00535 
00536                         MSG("DmxUtilities", Msg::kDebug) << "plane number " << currentPlane->GetPlaneNumber()
00537                                                                                         << " charge " << currentPlane->GetPlaneCharge()
00538                                                                                         << " type " 
00539                                                                                         << DmxPlaneTypes::GetPlaneTypeString(currentPlane->GetPlaneType()) 
00540                                                                                         << " num strips " << currentPlane->GetNumberOfStrips()
00541                                                                                         << endl;
00542                         endPlane = currentPlane->GetPlaneNumber();
00543 
00544                         //check to see that its a muon plane or only has 2 strips hit
00545                         if( currentPlane->GetNumberOfStrips() < 3 )     ++muonCtr;
00546                         
00547                         if(currentPlane->GetPlaneCharge() > (peGainConversion * 3.)) ++peCtr;
00548                         
00549                         MSG("DmxUtilities", Msg::kDebug) << " muon planes = " << muonCtr 
00550                                                                                         << " peCtr = " << peCtr << endl;
00551                         planeItr.Next();
00552                         ++ctr;
00553                 }
00554     
00555                 if(muonCtr == 4 && peCtr >= 3 && endPlane-startPlane<=5) 
00556                         fMuonStartPlaneNumber = startPlane;
00557     
00558                 //if there are still planes left to examine, the iterator will be valid
00559                 //otherwise, you ran out of planes, so break the cycle
00560                 if(planeItr.IsValid() ){
00561                         planeItr.Prev();
00562                         planeItr.Prev();
00563                         planeItr.Prev();
00564                         planeItr.Prev();
00565                 }
00566         }//end loop to find the muon track starting plane
00567   
00568         planeItr.GetSet()->ClearSlice(false);
00569         planeItr.Reset();
00570   
00571         return fMuonStartPlaneNumber;
00572 }

Int_t DmxUtilities::FindPlanesOffFit ( DmxPlaneItr &  planeItr,
Int_t  strayCut 
)

Definition at line 1089 of file DmxUtilities.cxx.

References DmxPlane::GetCoG(), DmxPlane::GetPlaneType(), DmxPlane::IsValid(), and DmxPlaneTypes::kShower.

Referenced by AlgDeMuxCosmics::RunAlg(), AlgDeMuxBeam::RunAlg(), AlgDeMuxGolden::UseGoldenFit(), AlgDeMuxCosmics::UseSingleFit(), and AlgDeMuxCosmics::UseSlidingWindow().

01090 {
01091   
01092   planeItr.Reset();
01093 
01094   Int_t strayPlanes = 0;
01095 
01096   DmxPlane *plane = 0;
01097 
01098   //MSG("DmxUtil", Msg::kDebug) << "checking for stray planes in event " << endl;
01099   while( planeItr.IsValid() ){
01100     
01101     plane = planeItr.Ptr();
01102     
01103     if( plane->IsValid() ){
01104       //&& plane->GetPlaneNumber()>= firstPlane && plane->GetPlaneNumber()<= lastPlane){
01105       Double_t fitCoG = -1.;
01106       
01107       
01108       if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
01109         Double_t diff1 = 0.;
01110         Double_t diff2 = 0.;
01111         Double_t diff3 = 0.;
01112         fitCoG = plane->GetCoG();
01113 
01114         if( dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis() ){
01115           diff1 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG());
01116           //MSG("DmxUtil", Msg::kDebug)<<"best " << dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG()<<endl;
01117         }
01118         if( dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis() ){
01119           diff2 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG());
01120           //MSG("DmxUtil", Msg::kDebug)<<"2nd " <<dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG()<<endl;
01121         }
01122         if( dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis() ){
01123           diff3 = TMath::Abs(fitCoG-dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG());
01124           //MSG("DmxUtil", Msg::kDebug)<<"3rd " << dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG()<<endl;
01125         }
01126         
01127         //fStrayCut is in strip numbers, the 0.042 factor is the width of a strip
01128         if(diff1>strayCut*.042 && diff2>strayCut*.042 && diff3>strayCut*.042){ 
01129           //MSG("DmxUtil", Msg::kDebug) << "plane = " << plane->GetPlaneNumber() << " is stray "
01130           //              << diff1 << " " << diff2 << " " << diff3 << endl;
01131       
01132           ++strayPlanes; 
01133           plane->SetStray(true);
01134           // MSG("DmxUtil", Msg::kDebug) << "\tstray plane = " << plane->GetPlaneNumber()  
01135           //              << "\tfit cog = " << fitCoG << endl;
01136         }
01137       }//end if shower plane
01138      //  else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
01139 //      fitCoG = plane->GetCoG();
01140 //      if( TMath::Abs(fitCoG-dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG()) > strayCut*0.042){ 
01141 //        //MSG("DmxUtil", Msg::kDebug)<<"muon " << plane->GetCoG()<<endl;
01142 //        ++strayPlanes; 
01143 //        plane->SetStray(true);
01144 //        MSG("DmxUtil", Msg::kDebug) << "\tstray plane = " << plane->GetPlaneNumber() 
01145 //                                << "\tfit cog = " << fitCoG << " initial cog = " 
01146 //                                << dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG()
01147 //                                <<endl;
01148 //      }
01149 //       }//end if muon plane
01150     }//end if it is a valid plane
01151     planeItr.Next();
01152   }//end loop over planes
01153 
01154   planeItr.Reset();
01155 
01156   return strayPlanes;
01157 }

void DmxUtilities::FindQuadraticFit ( Double_t *  x,
Double_t *  y,
Double_t *  weight,
Int_t  nPoints,
Double_t &  a1,
Double_t &  a2,
Double_t &  a3,
Double_t &  chiSq 
)

Definition at line 1505 of file DmxUtilities.cxx.

01508 {
01509 
01510   //use method of determinants to do the least squares fit
01511   //this is from bevington chapter 7
01512   //
01513   //f_1(x) = 1
01514   //f_2(x) = x
01515   //f_3(x) = x^2
01516   
01517   //need to find the sums
01518   //SIGMA y_i*f_1(x_i)/weight_i
01519   //SIGMA y_i*f_2(x_i)/weight_i
01520   //SIGMA y_i*f_3(x_i)/weight_i
01521   //SIGMA f_1(x_i)*f_1(x_i)/weight_i
01522   //SIGMA f_1(x_i)*f_2(x_i)/weight_i
01523   //SIGMA f_1(x_i)*f_3(x_i)/weight_i
01524   //SIGMA f_2(x_i)*f_1(x_i)/weight_i
01525   //SIGMA f_2(x_i)*f_2(x_i)/weight_i
01526   //SIGMA f_2(x_i)*f_3(x_i)/weight_i
01527   //SIGMA f_3(x_i)*f_1(x_i)/weight_i
01528   //SIGMA f_3(x_i)*f_2(x_i)/weight_i
01529   //SIGMA f_3(x_i)*f_3(x_i)/weight_i
01530 
01531   //some of these are the same, ie f_2*f_1 = f_1*f_2 so only need
01532   //6 instead of 9 variables
01533 
01534   Double_t yf1divWeight = 0.;
01535   Double_t yf2divWeight = 0.;
01536   Double_t yf3divWeight = 0.;
01537   Double_t f1f1divWeight = 0.;
01538   Double_t f1f2divWeight = 0.;
01539   Double_t f1f3divWeight = 0.;
01540   Double_t f2f2divWeight = 0.;
01541   Double_t f2f3divWeight = 0.;
01542   Double_t f3f3divWeight = 0.;
01543   
01544   for(Int_t i=0; i<nPoints; i++){
01545     yf1divWeight += y[i]/(weight[i]*weight[i]);
01546     yf2divWeight += (y[i]*x[i])/(weight[i]*weight[i]);
01547     yf3divWeight += (y[i]*x[i]*x[i])/(weight[i]*weight[i]);
01548     f1f1divWeight += 1./(weight[i]*weight[i]);
01549     f1f2divWeight += x[i]/(weight[i]*weight[i]);
01550     f1f3divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
01551     f2f2divWeight += (x[i]*x[i])/(weight[i]*weight[i]);
01552     f2f3divWeight += (x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01553     f3f3divWeight += (x[i]*x[i]*x[i]*x[i])/(weight[i]*weight[i]);
01554   }
01555 
01556   //make a bunch of matricies to hold these values
01557   TMatrix deltaMatrix(3,3);
01558   TMatrix a1Matrix(3,3);
01559   TMatrix a2Matrix(3,3);
01560   TMatrix a3Matrix(3,3);
01561     
01562   //fill the matricies
01563   deltaMatrix(0,0) = f1f1divWeight;
01564   deltaMatrix(0,1) = f1f2divWeight;
01565   deltaMatrix(0,2) = f1f3divWeight;
01566   deltaMatrix(1,0) = f1f2divWeight;
01567   deltaMatrix(1,1) = f2f2divWeight;
01568   deltaMatrix(1,2) = f2f3divWeight;
01569   deltaMatrix(2,0) = f1f3divWeight;
01570   deltaMatrix(2,1) = f2f3divWeight;
01571   deltaMatrix(2,2) = f3f3divWeight;
01572   
01573   a1Matrix(0,0) = yf1divWeight;
01574   a1Matrix(0,1) = f1f2divWeight;
01575   a1Matrix(0,2) = f1f3divWeight;
01576   a1Matrix(1,0) = yf2divWeight;
01577   a1Matrix(1,1) = f2f2divWeight;
01578   a1Matrix(1,2) = f2f3divWeight;
01579   a1Matrix(2,0) = yf3divWeight;
01580   a1Matrix(2,1) = f2f3divWeight;
01581   a1Matrix(2,2) = f3f3divWeight;
01582   
01583   a2Matrix(0,0) = f1f1divWeight;
01584   a2Matrix(0,1) = yf1divWeight;
01585   a2Matrix(0,2) = f1f3divWeight;
01586   a2Matrix(1,0) = f1f2divWeight;
01587   a2Matrix(1,1) = yf2divWeight;
01588   a2Matrix(1,2) = f2f3divWeight;
01589   a2Matrix(2,0) = f1f3divWeight;
01590   a2Matrix(2,1) = yf3divWeight;
01591   a2Matrix(2,2) = f3f3divWeight;
01592     
01593   a3Matrix(0,0) = f1f1divWeight;
01594   a3Matrix(0,1) = f1f2divWeight;
01595   a3Matrix(0,2) = yf1divWeight;
01596   a3Matrix(1,0) = f1f2divWeight;
01597   a3Matrix(1,1) = f2f2divWeight;
01598   a3Matrix(1,2) = yf2divWeight;
01599   a3Matrix(2,0) = f1f3divWeight;
01600   a3Matrix(2,1) = f2f3divWeight;
01601   a3Matrix(2,2) = yf3divWeight;
01602     
01603   a1 = -10000.;
01604   a2 = -10000.;
01605   a3 = -10000.;
01606   Double_t delta = deltaMatrix.Determinant();
01607   
01608   if(delta != 0.){
01609     a1 = a1Matrix.Determinant()/delta;
01610     a2 = a2Matrix.Determinant()/delta;
01611     a3 = a3Matrix.Determinant()/delta;
01612   }
01613 
01614   //find the chi^2 value for the fit
01615   chiSq = 0.;
01616   Double_t chi = 0.;
01617   for(Int_t j = 0; j<nPoints; j++){
01618     
01619     chi = (y[j] - (a1 + a2*x[j] + a3*(x[j]*x[j])))/weight[j];
01620     chiSq += chi*chi;
01621   
01622   }  
01623   
01624   //get the chiSq per degree of freedom
01625 
01626   chiSq /= (1.*nPoints);
01627 
01628   return;
01629 }

Int_t DmxUtilities::FindVertexPlane ( DmxPlaneItr &  planeItr  ) 

Definition at line 367 of file DmxUtilities.cxx.

References MuELoss::a, Msg::kDebug, and MSG.

Referenced by AlgDeMuxGolden::RunAlg(), AlgDeMuxCosmics::RunAlg(), and AlgDeMuxBeam::RunAlg().

00368 {
00369   MSG("DmxUtil", Msg::kDebug) << "In FindVertexPlane" << endl;
00370 
00371   planeItr.ResetFirst();
00372 
00373   Int_t vertexPlaneNumber = -1;
00374   Int_t firstPlaneNumber = 0;
00375   if( planeItr.IsValid() )
00376     firstPlaneNumber = planeItr.Ptr()->GetPlaneNumber();
00377 
00378   planeItr.ResetFirst();
00379 
00380   //find the vertex.  make it the first plane in a set of 5 planes where at least 3
00381   //of the planes are valid.  use the number 3 since it is the least number of planes
00382   //you can have and still drop one to refine a fit.  use 5 because the trigger is 4/5
00383   Int_t beginPlaneNumber = 0;
00384   Int_t endPlaneNumber = 0;
00385   Int_t validCtr = 0;
00386   Int_t cnt = 0;
00387 
00388   while( planeItr.IsValid() ){ 
00389 
00390     beginPlaneNumber = planeItr.Ptr()->GetPlaneNumber();
00391     endPlaneNumber = 0;
00392     validCtr = 0;
00393     
00394     cnt = 0;
00395     while( planeItr.IsValid() && cnt < 5 && validCtr < 3){
00396       if( planeItr.Ptr()->IsValid() ){ ++validCtr;}
00397       ++cnt;
00398       endPlaneNumber = planeItr.Ptr()->GetPlaneNumber();
00399       MSG("DmxUtil", Msg::kDebug)<<"Start Plane = " << beginPlaneNumber 
00400                                 << "\tThis Plane = " 
00401                                 << endPlaneNumber << "\tValid = " 
00402                                 << (Int_t)planeItr.Ptr()->IsValid() 
00403                                 << "\tValidCtr = " << validCtr << endl;
00404       planeItr.Next();
00405     }
00406         
00407     if( validCtr >= 3 && (endPlaneNumber-beginPlaneNumber) <= 5 ){  
00408       vertexPlaneNumber = beginPlaneNumber;
00409       MSG("DmxUtil", Msg::kDebug)<<"\tVertex Plane = " << vertexPlaneNumber << endl;
00410       break;
00411     }
00412     else{
00413       for(Int_t a = 0; a < cnt-1; a++){
00414         planeItr.Prev();
00415       }
00416     }
00417   } // end loop to find vertex
00418   
00419   //see if you can find any contiguos planes before the vertex, if so, make them the vertex
00420   //slice from 0 to the vertex and then iterate backwards.  if there is a gap in signal
00421   //of more than 2 planes, stop
00422   
00423   planeItr.Reset();
00424   planeItr.GetSet()->ClearSlice(false);
00425 
00426   Int_t nextPlane = 0;
00427 
00428   if( vertexPlaneNumber > -1 && firstPlaneNumber < vertexPlaneNumber){
00429     planeItr.GetSet()->Slice(firstPlaneNumber, vertexPlaneNumber);
00430     planeItr.GetSet()->Update();
00431 
00432     planeItr.ResetLast();
00433 
00434     while( planeItr.IsValid() ){
00435     
00436       Int_t currentPlane = planeItr.Ptr()->GetPlaneNumber();
00437       MSG("DmxUtil", Msg::kDebug) << "\tcurrentPlaneNumber = " << currentPlane;
00438       planeItr.Prev();
00439       nextPlane = 0;
00440       if(planeItr.IsValid()){
00441         nextPlane = planeItr.Ptr()->GetPlaneNumber();
00442         MSG("DmxUtil", Msg::kDebug) << "\tnextPlaneNumber = " << nextPlane << endl;
00443       
00444         if( vertexPlaneNumber - nextPlane < 3)
00445           vertexPlaneNumber = nextPlane;
00446         else break;
00447       }
00448       
00449     }//end loop over planes
00450 
00451     planeItr.GetSet()->ClearSlice(false);
00452     planeItr.ResetFirst();
00453   }
00454   
00455   planeItr.GetSet()->ClearSlice(false);
00456   planeItr.ResetFirst();
00457   
00458   MSG("DmxUtil", Msg::kDebug) << "\tvertexPlaneNumber = " << vertexPlaneNumber << endl;
00459   return vertexPlaneNumber;
00460 }

Int_t DmxUtilities::GetDigitPixel ( Int_t  channel  ) 

Definition at line 310 of file DmxUtilities.cxx.

Referenced by DmxDeMuxModule::Ana(), DmxDeMuxCosmicsModule::Ana(), FillHitPixels(), and FillPlaneArray().

00311 {
00312         
00313   //find the pixel that the digit came from.  make stupid arrays
00314   //to hold the pixel to va channel info, as for some reason it 
00315   //isnt available from the framework
00316   Int_t vaChip0[] = {-1, -1, 14, 0, 15, 1, 10, 4, 11, 
00317                      5, 6, 8, 7, 9, 2, 13, 3, 12};
00318   //Int_t vaChip1[] = {-1, -1, 22, 8, 23, 9, 18, 12, 19, 
00319   //         13, 14, 16, 15, 17, 10, 21, 11, 20};
00320   //Int_t vaChip2[] = {-1, -1, 6, 16, 7, 17, 2, 20, 3, 21, 
00321   //         22, 0, 23, 1, 18, 5, 19, 4};
00322     
00323   Int_t pixel = vaChip0[channel];
00324   //if(chip == 0){ pixel = vaChip0[channel];}
00325   //if(chip == 1){ pixel = vaChip1[channel];}
00326   //if(chip == 2){ pixel = vaChip2[channel];}
00327  
00328   //MSG("Dmx2", Msg::kDebug)<< "pixel = " << pixel << endl;
00329   return pixel;
00330 }

Bool_t DmxUtilities::IsOverlappingMultiple ( DmxPlaneItr &  planeItr,
Float_t  vertexZ,
Float_t  slopeCutOff,
Float_t  interceptCutOff,
Float_t  peakCutOff 
)

Definition at line 980 of file DmxUtilities.cxx.

References DmxPlane::GetPlaneNumber(), DmxPlane::GetPlaneType(), DmxPlaneTypes::kMuon, DmxPlaneTypes::kShower, and Munits::m.

Referenced by AlgDeMuxGolden::RunAlg(), and AlgDeMuxCosmics::RunAlg().

00983 {
00984 
00985   planeItr.ResetFirst();
00986   DmxPlane *plane = 0;
00987 
00988   //make a histogram for each supermodule.  this is because the method is trickier in 
00989   //tpos-z space than strip-plane space, so just deal with the gap by pretending its not
00990   //there
00991   TH2F houghSpaceSM1 = TH2F("houghSpaceSM1", "Hough Space for Multiples", 
00992                             30, -30, 30, 500, -500, 500);
00993   TH2F houghSpaceSM2 = TH2F("houghSpaceSM2", "Hough Space for Multiples", 
00994                             30, -30, 30, 500, -500, 500);
00995   TH2F houghSpaceCutSM1 = TH2F("hS75SM1", "Hough Space for U with cut", 
00996                                30, -30, 30, 500, -500, 500);
00997   TH2F houghSpaceCutSM2 = TH2F("hS75SM2", "Hough Space for U with cut", 
00998                                30, -30, 30, 500, -500, 500);
00999       
01000   Float_t planeZ = 0.;
01001   Float_t cog = 0.;
01002   while( (plane = planeItr()) ){
01003     
01004     planeZ = 1.*plane->GetPlaneNumber();
01005     
01006     if(plane->GetPlaneType() == DmxPlaneTypes::kShower){
01007      
01008       for(Int_t j = 0; j <3; j++){
01009         if( j==0 && dynamic_cast<DmxShowerPlane *>(plane)->GetBestHypothesis() ){
01010           cog = dynamic_cast<DmxShowerPlane *>(plane)->GetBestCoG();
01011         }
01012         else if( j==1 && dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestHypothesis() ){
01013           cog = dynamic_cast<DmxShowerPlane *>(plane)->GetSecondBestCoG();
01014         }
01015         else if( j==2 && dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestHypothesis() ){
01016           cog = dynamic_cast<DmxShowerPlane *>(plane)->GetThirdBestCoG();
01017         }
01018         //loop over possible slopes
01019         for(Int_t m = -29; m<30; m+=2){
01020           
01021           Float_t c = cog - m*1.*(planeZ - vertexZ);
01022           
01023           if(planeZ<249) houghSpaceSM1.Fill(m,c);
01024           else if(planeZ>249) houghSpaceSM2.Fill(m,c);
01025         
01026         }//end loop over slopes
01027       }//end loop over hypotheses
01028       
01029     }//end if shower plane
01030     else if(plane->GetPlaneType() == DmxPlaneTypes::kMuon){
01031       cog = dynamic_cast<DmxMuonPlane *>(plane)->GetInitialCoG();
01032       
01033       //loop over possible slopes
01034       for(Int_t m = -29; m<30; m+=2){
01035         
01036         Float_t c = cog - m*1.*(planeZ - vertexZ);
01037 
01038         //fill the space with this entry 3 times so it is weighted the same as a shower plane
01039         if(planeZ<249){
01040           houghSpaceSM1.Fill(m,c);
01041           houghSpaceSM1.Fill(m,c);
01042           houghSpaceSM1.Fill(m,c);
01043         }
01044         else if(planeZ>249){
01045           houghSpaceSM2.Fill(m,c);
01046           houghSpaceSM2.Fill(m,c);
01047           houghSpaceSM2.Fill(m,c);
01048         }
01049       }//end loop over slopes
01050             
01051     }//end if muon plane
01052   }//end loop over planes
01053 
01054   planeItr.ResetFirst();
01055 
01056   //now fill a second histogram with the bins that are at least peakCutOff of the maximum
01057   //bin in the previous histogram
01058   Stat_t binContent = 0;
01059   Double_t maximumSM1 = houghSpaceSM1.GetMaximum();     
01060   Double_t maximumSM2 = houghSpaceSM2.GetMaximum();     
01061         
01062   for(Int_t x = 1; x < 30+1; x++){
01063     //cout << x << endl;
01064     for(Int_t y = 1; y < 500+1; y++){
01065       binContent = houghSpaceSM1.GetBinContent(x,y);
01066       if(binContent>= peakCutOff*maximumSM1) houghSpaceCutSM1.SetBinContent(x,y,binContent);
01067       binContent = houghSpaceSM2.GetBinContent(x,y);
01068       if(binContent>= peakCutOff*maximumSM2) houghSpaceCutSM2.SetBinContent(x,y,binContent);
01069           
01070     }
01071   }
01072 
01073   //if the rms in the slope is less than the value determined by tests
01074   //and the rms in the intercept is greater than the cutoff, mark the event
01075   //as a possible overlapping multiple
01076   bool multiple = false;
01077   if(houghSpaceCutSM1.GetRMS(1)<=slopeCutOff 
01078      && houghSpaceCutSM1.GetRMS(2)>=interceptCutOff)
01079     multiple = true;
01080   else if(houghSpaceCutSM2.GetRMS(1)<=slopeCutOff 
01081           && houghSpaceCutSM2.GetRMS(2)>=interceptCutOff)
01082     multiple = true;
01083 
01084   return multiple;
01085 }

Bool_t DmxUtilities::IsValidFit ( DmxPlaneItr &  planeItr,
Double_t  a1,
Double_t  a2,
Double_t  a3,
Double_t  a4 
)

Definition at line 726 of file DmxUtilities.cxx.

References DmxPlane::GetZPosition().

00728 {
00729         Int_t nonPhysicalCtr = 0;
00730         bool validFit = true;
00731         Int_t planesInEvent = planeItr.SizeSelect();
00732         //MSG("DmxUtil", Msg::kDebug)<<"\tin IsValidFit " << a << "\t" << b << endl;
00733         planeItr.Reset();
00734 
00735         DmxPlane *plane = 0;
00736         Double_t fitCoG = 0.;
00737 
00738         while( planeItr.IsValid() ){
00739    
00740                 plane = planeItr.Ptr();
00741     
00742                 fitCoG = (a1 + a2*(plane->GetZPosition()) 
00743                                   + a3*TMath::Power(plane->GetZPosition(),2)
00744                                   + a4*TMath::Power(plane->GetZPosition(),3));
00745                 //MSG("DmxUtil", Msg::kDebug)<<"\t" << plane->GetZPosition() << endl;
00746       
00747                 //if the fit CoG for the plane is outside of the detector, it is non Physical
00748                 if( TMath::Abs(fitCoG) >= 4.5) ++nonPhysicalCtr;
00749   
00750                 planeItr.Next(); 
00751         }
00752 
00753         planeItr.Reset();
00754   
00755         if(planesInEvent<8 && (1.*nonPhysicalCtr)/(1.*planesInEvent)>0.33) validFit = false;
00756         else if(planesInEvent>=8 && nonPhysicalCtr>=4) validFit = false;
00757   
00758         return validFit;
00759 }

Bool_t DmxUtilities::IsValidFit ( DmxPlaneItr &  planeItr,
Double_t  a1,
Double_t  a2,
Double_t  a3,
Double_t  a4,
Float_t  offset 
)

Definition at line 689 of file DmxUtilities.cxx.

References DmxPlane::GetZPosition().

Referenced by AlgDeMuxCosmics::FindCosmicSolution(), and AlgDeMuxBeam::FindFit().

00691 {
00692         Int_t nonPhysicalCtr = 0;
00693         bool validFit = true;
00694         Int_t planesInEvent = planeItr.SizeSelect();
00695         //MSG("DmxUtil", Msg::kDebug)<<"\tin IsValidFit " << a << "\t" << b << endl;
00696         planeItr.Reset();
00697 
00698         DmxPlane *plane = 0;
00699         Double_t fitCoG = 0.;
00700         
00701         while( planeItr.IsValid() ){
00702    
00703                 plane = planeItr.Ptr();
00704     
00705                 fitCoG = (a1 + a2*(plane->GetZPosition()-offset) 
00706                                   + a3*TMath::Power(plane->GetZPosition()-offset,2)
00707                                   + a4*TMath::Power(plane->GetZPosition()-offset,3));
00708       //MSG("DmxUtil", Msg::kDebug)<<"\t" << plane->GetZPosition() << endl;
00709       
00710       //if the fit CoG for the plane is outside of the detector, it is non Physical
00711       if( TMath::Abs(fitCoG) >= 4.5) ++nonPhysicalCtr;
00712   
00713     planeItr.Next(); 
00714   }
00715 
00716   planeItr.Reset();
00717   
00718   if(planesInEvent<8 && (1.*nonPhysicalCtr)/(1.*planesInEvent)>0.33) validFit = false;
00719   else if(planesInEvent>=8 && nonPhysicalCtr>=4) validFit = false;
00720   
00721   return validFit;
00722 }


The documentation for this class was generated from the following files:
Generated on Mon Sep 1 00:51:29 2014 for loon by  doxygen 1.4.7