NtpTimingFit Class Reference

#include <NtpTimingFit.h>

List of all members.

Public Member Functions

void SetDetector (Detector::Detector_t myDetector)
void UseTimeWalk (Bool_t yesno=1)
void UseTimeConstants (Bool_t yesno=1)
void UseSpectrometer (Bool_t yesno=1)
void UseCalculator (Bool_t yesno=1)
void UseBeamMode ()
void UseCosmicMode ()
void AddEvent (NtpStRecord *ntpRecord, Int_t ievent=0)
void AddData (Int_t plane, Int_t strip, Int_t eastwest, Int_t view, Int_t digits, Double_t u_metres, Double_t v_metres, Double_t z_metres, Double_t t_nanosec, Double_t q_pes)
void AddData (Int_t plane, Int_t strip, Int_t eastwest, Int_t view, Int_t digits, Double_t ds_metres, Double_t dt_nanosecs, Double_t q_pes)
void Reset ()
void ResetTimeStrips ()
void ResetTimePoints ()
void PrintTimeStrips ()
void PrintTimePoints ()
Int_t GetPlane (UInt_t ihit)
Int_t GetStrip (UInt_t ihit)
Int_t GetEnd (UInt_t ihit)
Int_t GetDigits (UInt_t ihit)
Double_t GetS (UInt_t ihit)
Double_t GetT (UInt_t ihit)
Double_t GetQ (UInt_t ihit)
Double_t GetW (UInt_t ihit)
UInt_t GetNumPoints ()
Int_t GetMinPlane ()
Int_t GetMaxPlane ()
Int_t GetTrackPlanes ()
Double_t GetMinZ ()
Double_t GetMinTime ()
Double_t GetTrigTime ()
Double_t GetGateTime ()
Double_t GetLength ()
Bool_t CheckData ()
Bool_t RunLinearFit (Double_t &invbeta, Double_t &vtime)
Bool_t RunLinearFit (Double_t &invbeta, Double_t &vtime, Double_t &err, Double_t &rms, Double_t &chi2)
Bool_t RunConstrainedFit (Double_t &invbeta, Double_t &vtime)
Bool_t RunConstrainedFit (Double_t &invbeta, Double_t &vtime, Double_t &err, Double_t &rms, Double_t &chi2)
Bool_t RunConstrainedFit (Int_t dir, Double_t &invbeta, Double_t &vtime, Double_t &err, Double_t &rms, Double_t &chi2)
Bool_t CalcSeedTime (Int_t dir, Double_t &vtime, Double_t &window)
void SetDebug (Bool_t onoff=1)

Static Public Member Functions

static NtpTimingFitInstance ()
static void SelectBeamMode ()
static void SelectCosmicMode ()
static Bool_t RunFit (NtpStRecord *ntpRecord, Int_t ievent, Double_t &vtime, Double_t &err, Double_t &chi2)
static Bool_t RunBeamFit (NtpStRecord *ntpRecord, Int_t ievent, Double_t &vtime, Double_t &err, Double_t &chi2)
static Bool_t RunCosmicFit (NtpStRecord *ntpRecord, Int_t ievent, Double_t &vtime, Double_t &err, Double_t &chi2)

Private Member Functions

 NtpTimingFit ()
 ~NtpTimingFit ()
void CalcTimePoints ()
Double_t GetTimeConstant (Int_t plane, Int_t strip, Int_t end)
Double_t GetTimeWalk (Int_t ndigits, Double_t q)
Double_t GetTimeWeight (Double_t q)

Private Attributes

Double_t fTimeResolution
Int_t fMinPlaneWindow
Int_t fMaxPlaneWindow
TimeStrip fTimeStrip
TimePoint fTimePoint
std::vector< TimeStripvecTimeStrips
std::vector< TimePointvecTimePoints
std::vector< Double_t > vecResiduals
std::vector< Double_t > vecResidualsForFit
Int_t fDetector
Int_t fMinPlane
Int_t fMaxPlane
Int_t fTrackPlanes
Double_t fMinZ
Double_t fMinTime
Double_t fTrigTime
Double_t fGateTime
Double_t fCrateT0
Double_t fLength
Double_t fDelta
Bool_t fDebug
Bool_t fFoundMode
Bool_t fUseBeamMode
Bool_t fUseTimeWalk
Bool_t fUseTimeConstants
Bool_t fUseSpectrometer
Bool_t fUseCalculator

Detailed Description

Definition at line 36 of file NtpTimingFit.h.


Constructor & Destructor Documentation

NtpTimingFit::NtpTimingFit (  )  [private]

Definition at line 31 of file NtpTimingFit.cxx.

References fCrateT0, fDebug, fDetector, fFoundMode, fLength, fMaxPlane, fMaxPlaneWindow, fMinPlane, fMinPlaneWindow, fMinTime, fMinZ, fTimeResolution, fTrackPlanes, fTrigTime, fUseBeamMode, fUseCalculator, fUseSpectrometer, fUseTimeConstants, fUseTimeWalk, and Detector::kUnknown.

00032 {
00033   std::cout << " *** Loading Timing Fit *** " << std::endl;
00034 
00035   fDetector = Detector::kUnknown;
00036 
00037   fDebug = 0;
00038   
00039   fMinPlane = -1;
00040   fMaxPlane = -1;
00041   fTrackPlanes = 0;
00042   
00043   fMinZ = -99999.9;
00044   fLength = 0.0;
00045 
00046   fTimeResolution = 50.0;
00047  
00048   fMinPlaneWindow = 0;
00049   fMaxPlaneWindow = 500;
00050 
00051   fMinTime = -99999.9;
00052   fTrigTime = 0.0;
00053   fCrateT0 = 0.0;
00054 
00055   fFoundMode = 0;
00056 
00057   fUseBeamMode = 1;       // forward-going
00058   fUseTimeWalk = 0;       // apply time walk
00059   fUseTimeConstants = 0;  // apply time constants
00060   fUseSpectrometer = 1;   // use spectrometer
00061 
00062   fUseCalculator = 0;     // re-calculate dS
00063 }

NtpTimingFit::~NtpTimingFit (  )  [private]

Definition at line 65 of file NtpTimingFit.cxx.

00066 {
00067 
00068 }


Member Function Documentation

void NtpTimingFit::AddData ( Int_t  plane,
Int_t  strip,
Int_t  eastwest,
Int_t  view,
Int_t  digits,
Double_t  ds_metres,
Double_t  dt_nanosecs,
Double_t  q_pes 
)

Definition at line 447 of file NtpTimingFit.cxx.

References TimePoint::digits, TimePoint::ds, TimePoint::dt, TimePoint::end, fLength, fTimePoint, TimePoint::plane, TimePoint::q, TimePoint::strip, vecTimePoints, and TimePoint::view.

00448 {
00449   if( plane>=0 && plane<500
00450    && q_pes>2.0 ){  // <-- 2PE cut
00451 
00452     fTimePoint.plane  = plane;
00453     fTimePoint.strip  = strip;
00454     fTimePoint.end    = eastwest;
00455     fTimePoint.view   = view;
00456     fTimePoint.digits = digits;
00457     fTimePoint.ds     = ds_metres;
00458     fTimePoint.dt     = dt_nanosecs;
00459     fTimePoint.q      = q_pes;
00460 
00461     if( ds_metres>fLength ){
00462       fLength = ds_metres;
00463     }
00464 
00465     vecTimePoints.push_back(fTimePoint);
00466   }
00467 }

void NtpTimingFit::AddData ( Int_t  plane,
Int_t  strip,
Int_t  eastwest,
Int_t  view,
Int_t  digits,
Double_t  u_metres,
Double_t  v_metres,
Double_t  z_metres,
Double_t  t_nanosec,
Double_t  q_pes 
)

Definition at line 421 of file NtpTimingFit.cxx.

References TimeStrip::digits, TimeStrip::end, fMaxPlane, fMinPlane, fMinTime, fMinZ, fTimeStrip, TimeStrip::plane, TimeStrip::q, TimeStrip::strip, TimeStrip::t, TimeStrip::u, TimeStrip::v, vecTimeStrips, TimeStrip::view, and TimeStrip::z.

Referenced by AddEvent(), and CalcTimePoints().

00422 {
00423   if( plane>=0 && plane<500
00424    && q_pes>2.0 ){  // <-- 2PE cut
00425 
00426     fTimeStrip.plane  = plane;
00427     fTimeStrip.strip  = strip;
00428     fTimeStrip.end    = eastwest;
00429     fTimeStrip.view   = view;
00430     fTimeStrip.digits = digits;
00431     fTimeStrip.u      = u_metres; 
00432     fTimeStrip.v      = v_metres; 
00433     fTimeStrip.z      = z_metres; 
00434     fTimeStrip.t      = t_nanosec; 
00435     fTimeStrip.q      = q_pes;
00436 
00437     if( fMinZ<0    || z_metres<fMinZ ) fMinZ = z_metres;
00438     if( fMinTime<0 || t_nanosec<fMinTime ) fMinTime = t_nanosec;
00439 
00440     if( fMinPlane<0 || plane<fMinPlane ) fMinPlane = plane;
00441     if( fMaxPlane<0 || plane>fMaxPlane ) fMaxPlane = plane;
00442 
00443     vecTimeStrips.push_back(fTimeStrip);
00444   }
00445 }

void NtpTimingFit::AddEvent ( NtpStRecord ntpRecord,
Int_t  ievent = 0 
)

Definition at line 205 of file NtpTimingFit.cxx.

References AddData(), CalcTimePoints(), NtpSRTimeStatus::crate_t0_ns, digits(), NtpSRTrack::ds, NtpSRTrack::end, NtpStRecord::evt, NtpStRecord::evthdr, fCrateT0, fDebug, fGateTime, fMaxPlaneWindow, fMinPlaneWindow, fMinTime, fMinZ, fTrackPlanes, fTrigTime, fUseCalculator, fUseTimeConstants, fUseTimeWalk, VldContext::GetDetector(), RecRecordImp< T >::GetHeader(), RecDataHeader::GetRun(), RecPhysicsHeader::GetSnarl(), GetTimeConstant(), GetTimeWalk(), RecHeader::GetVldContext(), Detector::kNear, n, NtpSRStrip::ndigit, NtpSREventSummary::nevent, Munits::ns, NtpSRTrack::nstrip, NtpSRPulseHeight::pe, NtpSRStrip::ph0, NtpSRStrip::ph1, NtpSRStrip::plane, NtpSRVertex::plane, NtpSRStrip::planeview, NtpSelection::PrimaryTrack(), PrintTimePoints(), PrintTimeStrips(), Reset(), SetDetector(), NtpSRTimeStatus::sgate_53mhz, NtpStRecord::stp, NtpSRTrack::stp, NtpSRTrack::stpds, NtpSRTrack::stpt0, NtpSRTrack::stpt1, NtpSRTrack::stpu, NtpSRTrack::stpv, NtpSRTrack::stpz, NtpSRStrip::strip, NtpStRecord::timestatus, NtpSREventSummary::trigtime, NtpStRecord::trk, and NtpSRTrack::vtx.

Referenced by NtpTimingNtuple::Fill(), NtpEventDisplayTiming::MakePicture(), and RunFit().

00206 {  
00207   // reset timing fit
00208   // ================
00209   this->Reset();  
00210 
00211   // sanity check
00212   // ============
00213   if( ntpRecord==0 || ievent<0 ) return;
00214 
00215   // run/snarl/event
00216   // ===============
00217   if( fDebug ) std::cout << " --- NtpTimingFit::AddEvent(...) --- " << std::endl;
00218   if( fDebug ) std::cout << "   Event: " << ntpRecord->GetHeader().GetRun() 
00219                                   << "/" << ntpRecord->GetHeader().GetSnarl() 
00220                                   << "/" << ievent << std::endl;
00221 
00222   // detector type
00223   Detector::Detector_t myDetector = ntpRecord->GetHeader().GetVldContext().GetDetector();
00224   SetDetector(myDetector);
00225 
00226   // crate T0
00227   fCrateT0 = ntpRecord->timestatus.crate_t0_ns; // nsec
00228 
00229   // trigger time
00230   fTrigTime = 1.0e9*ntpRecord->evthdr.trigtime  // sec -> nsec
00231             - fCrateT0;
00232 
00233   // ND only: Latch time of sgate on 53.1 MHz Clock
00234   if( myDetector==Detector::kNear ){
00235     fGateTime = (1000./53.1)*ntpRecord->timestatus.sgate_53mhz; // nsec
00236   }
00237   else{
00238     fGateTime = 0.0;
00239   }
00240 
00241   // primary event
00242   // =============
00243   Int_t nevent = ntpRecord->evthdr.nevent;  
00244   if( nevent<=0 ){
00245     if( fDebug ) std::cout << " <warning> This snarl has no events! [return] " << std::endl;
00246     return;
00247   }
00248 
00249   TClonesArray* eventArray  = (TClonesArray*)(ntpRecord->evt);
00250 
00251   if( ievent<0 || ievent>eventArray->GetLast() ){
00252     if( fDebug ) std::cout << " <warning> Have overshot this snarl! [return] " << std::endl;
00253     return;
00254   }
00255 
00256   //
00257   // NtpSREvent* event = (NtpSREvent*)(eventArray->At(ievent));
00258   //
00259 
00260   // primary track
00261   // =============
00262   TClonesArray* trackArray  = (TClonesArray*)(ntpRecord->trk);
00263 
00264   Int_t primtrkid = NtpSelection::PrimaryTrack( ntpRecord, ievent );
00265 
00266   if( primtrkid<0 ){
00267     if( fDebug ) std::cout << " <warning> Event has no Primary Track [return] " << std::endl;
00268     return;
00269   }
00270 
00271   NtpSRTrack* track = (NtpSRTrack*)(trackArray->At(primtrkid));
00272 
00273   // loop over track strips
00274   // ======================
00275   Int_t track_nstrip = track->nstrip;
00276   Int_t vtxPlane = track->vtx.plane;
00277   Int_t endPlane = track->end.plane;
00278   Int_t trkSpanPlanes = 1+abs(endPlane-vtxPlane);
00279   Float_t s0 = track->ds; 
00280 
00281   if( fDebug ) std::cout << "   Fitting Track: [" << ievent << "/" << primtrkid << "] VtxPlane=" << vtxPlane << " EndPlane=" << endPlane << std::endl;
00282 
00283   if( trkSpanPlanes<5 ){
00284     std::cout << " <warning> Primary Track Too Short (" << trkSpanPlanes << " planes) [fail] " << std::endl;
00285     return;
00286   }
00287 
00288   if( !(vtxPlane>fMinPlaneWindow && vtxPlane<fMaxPlaneWindow)
00289    && !(endPlane>fMinPlaneWindow && endPlane<fMaxPlaneWindow) ){
00290     std::cout << " <warning> Not in Region of Interest (Vtx=" << vtxPlane<< ",End=" << endPlane << ") [fail] " << std::endl;
00291     return;
00292   }
00293 
00294   // strip information
00295   TClonesArray* stripArray  = (TClonesArray*)(ntpRecord->stp);
00296   Int_t* track_stp = track->stp;
00297   Float_t* track_stpu = track->stpu;
00298   Float_t* track_stpv = track->stpv;
00299   Float_t* track_stpz = track->stpz;
00300   Float_t* track_stpds = track->stpds;
00301   Double_t* track_stpt0 = track->stpt0; 
00302   Double_t* track_stpt1 = track->stpt1; 
00303 
00304   // some other arrays
00305   // Float_t* track_stpx = track->stpx;
00306   // Float_t* track_stpy = track->stpy;
00307   // Double_t* track_stptcal0t0 = track->stptcal0t0; 
00308   // Double_t* track_stptcal1t0 = track->stptcal1t0; 
00309 
00310   // loop over strips
00311   Int_t trkHitPlanes = 0;
00312 
00313   Int_t N = 1+fMaxPlaneWindow-fMinPlaneWindow;
00314   Double_t* ph = new Double_t[N];
00315 
00316   for( Int_t n=0; n<N; n++ ){
00317     ph[n] = 0.0;
00318   }
00319 
00320   for( Int_t ns=0; ns<track_nstrip; ns++ ){
00321     Int_t ptr = track_stp[ns];
00322     if( ptr>=0 ){
00323       NtpSRStrip* strip = (NtpSRStrip*)(stripArray->At(ptr));
00324       Int_t planeNumber = strip->plane;
00325       Double_t q = strip->ph0.pe + strip->ph1.pe;
00326       
00327       if( planeNumber>=fMinPlaneWindow 
00328        && planeNumber<fMaxPlaneWindow ){
00329         ph[planeNumber-fMinPlaneWindow] += q;
00330       }
00331     }
00332   }
00333 
00334   for( Int_t n=0; n<N; n++ ){
00335     if( ph[n]>0.0 ) trkHitPlanes++;
00336   }
00337 
00338   delete [] ph;
00339 
00340   // set track properties
00341   fTrackPlanes = trkSpanPlanes;
00342 
00343   // loop over timing information
00344   for( Int_t ns=0; ns<track_nstrip; ns++ ){
00345     Int_t ptr = track_stp[ns];
00346     if( ptr>=0 ){
00347 
00348       // Information from Strip
00349       NtpSRStrip* strip = (NtpSRStrip*)(stripArray->At(ptr));
00350       Int_t planeNumber = strip->plane;
00351       Int_t stripNumber = strip->strip;
00352       Int_t viewNumber  = strip->planeview;
00353       Int_t digits = strip->ndigit;
00354       Double_t q0 = strip->ph0.pe;
00355       Double_t q1 = strip->ph1.pe;
00356 
00357       // Information from Track [Positions]
00358       Float_t u = track_stpu[ns];
00359       Float_t v = track_stpv[ns];
00360       Float_t z = track_stpz[ns];
00361       
00362       // Information from Track [Corrected Times]
00363       Double_t calT0 = 1.0e9*track_stpt0[ns];  // sec -> nsec
00364       Double_t calT1 = 1.0e9*track_stpt1[ns];  // sec -> nsec
00365 
00366       // Apply Time Walk Correction
00367       if( fUseTimeWalk ){
00368         if( q0>0.0 ) calT0 -= GetTimeWalk(digits,q0);
00369         if( q1>0.0 ) calT1 -= GetTimeWalk(digits,q1);
00370       }
00371 
00372       if( fUseTimeConstants ){
00373         if( q0>0.0 ) calT0 -= GetTimeConstant(planeNumber,stripNumber,0);
00374         if( q1>0.0 ) calT1 -= GetTimeConstant(planeNumber,stripNumber,1);
00375       }
00376 
00377       // Information from Strips [Distance Along Track]
00378       Float_t s = s0 - track_stpds[ns];        // from end -> from vertex
00379 
00380       // East 
00381       if( q0>2.0 ){ // <-- 2PE cut
00382         this->AddData( planeNumber, stripNumber, 0, viewNumber, digits,
00383                        u, v, z,
00384                        calT0, q0);
00385         if( s>=0.0 ){
00386           this->AddData( planeNumber, stripNumber, 0, viewNumber, digits,
00387                          s, calT0, q0 );
00388         }
00389       }
00390 
00391       // West
00392       if( q1>2.0 ){ // <-- 2PE cut
00393         this->AddData( planeNumber, stripNumber, 1, viewNumber, digits,
00394                        u, v, z,
00395                        calT1, q1);
00396         if( s>=0.0 ){
00397           this->AddData( planeNumber, stripNumber, 1, viewNumber, digits,
00398                          s, calT1, q1 );
00399         }
00400       }
00401     }
00402   }
00403 
00404   // recalculate time points
00405   // =======================
00406   if( fUseCalculator ){
00407     this->CalcTimePoints();
00408   }
00409 
00410   // print
00411   if( fDebug ){
00412     this->PrintTimeStrips();
00413     this->PrintTimePoints();
00414   }
00415 
00416   if( fDebug ) std::cout << "   Results: TrigTime=" << fTrigTime << " MinTime=" << fMinTime << " MinZ=" << fMinZ << std::endl;
00417 
00418   return;
00419 }

Bool_t NtpTimingFit::CalcSeedTime ( Int_t  dir,
Double_t &  vtime,
Double_t &  window 
)

Definition at line 787 of file NtpTimingFit.cxx.

References CheckData(), CompareTimes(), fDebug, fDelta, fMaxPlaneWindow, fMinPlaneWindow, fTimeResolution, GetNumPoints(), GetPlane(), GetS(), GetT(), n, vecResiduals, and vecResidualsForFit.

Referenced by NtpTimingNtuple::Fill(), NtpEventDisplayTiming::MakePicture(), and RunConstrainedFit().

00788 {
00789   // default values 
00790   vtime = 0.0;  
00791   window = 3.0*fTimeResolution;
00792 
00793   // sanity check
00794   if( dir==0 ) return 0;
00795   
00796   // more checks
00797   if( CheckData()==0 ) return 0;  
00798 
00799   if( fDebug ) std::cout << " *** Calculate Seed Time *** " << std::endl;  
00800 
00801 
00802   // calculate residuals
00803   // ===================
00804   vecResiduals.clear();
00805 
00806   Int_t plane = 0;
00807 
00808   Double_t t = 0.0;
00809   Double_t s = 0.0;
00810 
00811   Double_t ibeta = 0.0;
00812   if( dir>0 ) ibeta = +1.0;
00813   if( dir<0 ) ibeta = -1.0; 
00814 
00815   for( UInt_t n=0; n<GetNumPoints(); n++ ){
00816     plane = GetPlane(n);
00817 
00818     if( plane>fMinPlaneWindow 
00819      && plane<fMaxPlaneWindow ){
00820       s = GetS(n);
00821       t = GetT(n);
00822         
00823       fDelta = t - (ibeta/0.3)*s;
00824       vecResiduals.push_back(fDelta);
00825     }
00826   }
00827 
00828   // time ordering
00829   std::sort(vecResiduals.begin(), vecResiduals.end(), CompareTimes);
00830 
00831   // print list
00832   if( fDebug ){
00833      std::cout << " --- all residuals --- " << std::endl;
00834      for( UInt_t n=0; n<vecResiduals.size(); n++ ){
00835        std::cout << "  [" << n << "] " << vecResiduals.at(n) - vecResiduals.at(0) << std::endl;
00836      }
00837   }
00838 
00839   // sanity check
00840   if( vecResiduals.size()<4 ) return 0;
00841 
00842 
00843   // select good residuals
00844   // =====================
00845   vecResidualsForFit.clear();
00846 
00847   for( UInt_t n=0; n<vecResiduals.size(); n++ ){
00848     fDelta = vecResiduals.at(n);
00849 
00850     if( ( n>0
00851        && fabs(fDelta-vecResiduals.at(n-1))<3.0*fTimeResolution )
00852      || ( n<vecResiduals.size()-1 
00853        && fabs(fDelta-vecResiduals.at(n+1))<3.0*fTimeResolution ) ){
00854       vecResidualsForFit.push_back(fDelta);
00855     }
00856     
00857   }
00858 
00859   // print list
00860   if( fDebug ){
00861      std::cout << " --- good residuals --- " << std::endl;
00862      for( UInt_t n=0; n<vecResidualsForFit.size(); n++ ){
00863        std::cout << "  [" << n << "] " << vecResidualsForFit.at(n) - vecResiduals.at(0) << std::endl;
00864      }
00865    }
00866 
00867   // sanity check
00868   if( vecResidualsForFit.size()<4 ) return 0;
00869 
00870 
00871   // find midpoint of good residuals
00872   // ===============================
00873   UInt_t ctr = (UInt_t)(0.50*vecResidualsForFit.size());
00874 
00875   if( vecResidualsForFit.size()%2==0 ){
00876     vtime = 0.5*( vecResidualsForFit.at(ctr-1)+vecResidualsForFit.at(ctr) );
00877   }
00878   else{
00879     vtime = vecResidualsForFit.at(ctr);
00880   }
00881 
00882   // calculate window
00883   // ================
00884   UInt_t ctr1 = (UInt_t)( 0.16*vecResidualsForFit.size() );
00885   UInt_t ctr2 = (UInt_t)( 0.84*vecResidualsForFit.size()+1.0 );
00886 
00887   if( ctr2>vecResidualsForFit.size()-1 ) ctr2 = vecResidualsForFit.size()-1;
00888 
00889   window = fTimeResolution 
00890          + 1.5*fabs( vecResidualsForFit.at(ctr2) - vecResidualsForFit.at(ctr1) );
00891 
00892   if( fDebug ) std::cout << "   Results: SeedTime=" << vtime << " Window=" << window << std::endl;
00893 
00894   return 1;
00895 }

void NtpTimingFit::CalcTimePoints (  )  [private]

Definition at line 556 of file NtpTimingFit.cxx.

References AddData(), digits(), TimeStrip::digits, TimeStrip::end, fDebug, fTimeStrip, n, TimeStrip::plane, TimeStrip::q, ResetTimePoints(), TimeStrip::strip, TimeStrip::t, TimeStrip::u, TimeStrip::v, vecTimeStrips, TimeStrip::view, MinosMaterial::Z(), and TimeStrip::z.

Referenced by AddEvent().

00557 {
00558   // reset time points
00559   this->ResetTimePoints();
00560 
00561   // sanity check, do we have any data?
00562   if( vecTimeStrips.size()==0 ){
00563     if( fDebug ) std::cout << " <warning> No Time Strips [return]" << std::endl;
00564     return;
00565   }
00566 
00567   // calculate 3D positions on track
00568   Double_t* U = new Double_t[500];
00569   Double_t* V = new Double_t[500];
00570   Double_t* Z = new Double_t[500];
00571   Double_t* Q = new Double_t[500];
00572   Double_t* S = new Double_t[500];
00573   
00574   for( Int_t i=0; i<500; i++ ){
00575     U[i] = 0.0;  
00576     V[i] = 0.0;
00577     Z[i] = 0.0;
00578     Q[i] = 0.0;
00579     S[i] = 0.0;
00580   }
00581 
00582   Int_t plane = 0;
00583   Int_t strip = 0;
00584   Int_t end = 0;
00585   Int_t view = 0;
00586   Int_t digits = 0;
00587 
00588   Double_t u = 0.0;
00589   Double_t v = 0.0;
00590   Double_t z = 0.0;
00591   Double_t t = 0.0;
00592   Double_t q = 0.0;
00593   Double_t s = 0.0;
00594 
00595   for( UInt_t n=0; n<vecTimeStrips.size(); n++ ){
00596     fTimeStrip = (TimeStrip)(vecTimeStrips.at(n)); 
00597 
00598     plane = fTimeStrip.plane;
00599     u     = fTimeStrip.u;
00600     v     = fTimeStrip.v;
00601     z     = fTimeStrip.z;
00602     q     = fTimeStrip.q;
00603     
00604     if( plane>=0 && plane<500 
00605      && q>2.0 ){  // <-- 2PE cut
00606       U[plane] += q*u;
00607       V[plane] += q*v;
00608       Z[plane] += q*z;
00609       Q[plane] += q;
00610     }
00611   }
00612 
00613   for( Int_t i=0; i<500; i++ ){
00614     if( Q[i]>0.0 ){
00615       U[i] /= Q[i];  
00616       V[i] /= Q[i];
00617       Z[i] /= Q[i];
00618     }
00619   }  
00620 
00621   // could try to smooth out positions here
00622 
00623 
00624 
00625   // calculate 3D displacements along track
00626   Bool_t firstHit = 1;
00627   Int_t j  = +1; 
00628   Int_t jm = -1;
00629 
00630   Double_t du = 0.0;
00631   Double_t dv = 0.0;
00632   Double_t dz = 0.0;
00633   Double_t ds = 0.0;
00634 
00635   for( Int_t i=0; i<500; i++){
00636     if( Q[i]>0.0 ){
00637       if( firstHit ){
00638         firstHit = 0;
00639       }
00640       else{
00641         j=1; jm=-1;
00642         while( i-j>=0 && jm<0 ){
00643           if( Q[i-j]>0.0 ) jm=j;
00644           j++;
00645         }          
00646         if( jm>0 ){
00647           du = U[i] - U[i-jm]; 
00648           dv = V[i] - V[i-jm];
00649           dz = Z[i] - Z[i-jm];
00650           ds = sqrt(du*du+dv*dv+dz*dz);
00651           S[i] = S[i-jm] + ds;
00652         }
00653       }
00654     }
00655   } 
00656 
00657   // now create time points
00658   for( UInt_t n=0; n<vecTimeStrips.size(); n++ ){
00659     fTimeStrip = (TimeStrip)(vecTimeStrips.at(n)); 
00660 
00661     plane  = fTimeStrip.plane;
00662     strip  = fTimeStrip.strip;
00663     end    = fTimeStrip.end;
00664     view   = fTimeStrip.view;
00665     digits = fTimeStrip.digits;
00666     q      = fTimeStrip.q;
00667     t      = fTimeStrip.t;
00668 
00669     if( plane>=0 && plane<500 
00670      && q>2.0 ){ // <-- 2PE cut
00671       s = S[plane];
00672       this->AddData( plane, strip, end, view, digits,
00673                      s, t, q );
00674     }
00675   }
00676 
00677   // delete arrays
00678   delete [] U;
00679   delete [] V;
00680   delete [] Z;
00681   delete [] Q;
00682   delete [] S;
00683 
00684   return;
00685 }

Bool_t NtpTimingFit::CheckData (  ) 

Definition at line 741 of file NtpTimingFit.cxx.

References fDetector, fFoundMode, GetLength(), GetNumPoints(), and Detector::kUnknown.

Referenced by CalcSeedTime(), RunConstrainedFit(), and RunLinearFit().

00742 {
00743   // check fitter
00744   if( fFoundMode==0 ){
00745     std::cout << " <warning> Need to Select Beam or Cosmic Mode [abort] " << std::endl;
00746     std::cout << "   NtpTimingFit::SelectBeamMode() OR NtpTimingFit::SelectCosmicMode() " << std::endl;
00747     assert( fFoundMode );
00748     return 0;
00749   }
00750 
00751   // unknown detector
00752   if( fDetector == Detector::kUnknown ){
00753     std::cout << " <warning> Unknown Detector (" << fDetector << ") [fail]" << std::endl;
00754     return 0;
00755   }
00756 
00757   // sanity check
00758   if( GetNumPoints()==0 ){
00759     return 0;
00760   }
00761 
00762   // sanity check
00763   if( GetNumPoints()<5 ){
00764     std::cout << " <warning> Insufficient Points for Fit (" <<  GetNumPoints() << " points) [fail] " << std::endl;
00765     return 0;
00766   }
00767 
00768   // sanity check  
00769   if( GetLength()<0.25 ){
00770     std::cout << " <warning> Track Too Short (" << GetLength() << " metres) [fail] " << std::endl;
00771     return 0;
00772   }
00773 
00774   if( GetLength()>75.0 ){
00775     std::cout << " <warning> Track Too Long (" << GetLength() << " metres) [fail] " << std::endl;
00776     return 0;
00777   }
00778 
00779   return 1;
00780 }

Int_t NtpTimingFit::GetDigits ( UInt_t  ihit  ) 

Definition at line 499 of file NtpTimingFit.cxx.

References TimePoint::digits, fTimePoint, GetNumPoints(), and vecTimePoints.

Referenced by NtpTimingNtuple::Fill().

00500 {
00501   if( ihit<GetNumPoints() ){
00502     fTimePoint = (TimePoint)(vecTimePoints.at(ihit));
00503     return fTimePoint.digits;
00504   }
00505 
00506   return -1;
00507 }

Int_t NtpTimingFit::GetEnd ( UInt_t  ihit  ) 

Definition at line 489 of file NtpTimingFit.cxx.

References TimePoint::end, fTimePoint, GetNumPoints(), and vecTimePoints.

Referenced by NtpTimingNtuple::Fill().

00490 {
00491   if( ihit<GetNumPoints() ){
00492     fTimePoint = (TimePoint)(vecTimePoints.at(ihit));
00493     return fTimePoint.end;
00494   }
00495 
00496   return -1;
00497 }

Double_t NtpTimingFit::GetGateTime (  )  [inline]

Definition at line 102 of file NtpTimingFit.h.

References fGateTime.

Referenced by NtpTimingNtuple::Fill().

00102 { return fGateTime; }

Double_t NtpTimingFit::GetLength (  )  [inline]

Definition at line 104 of file NtpTimingFit.h.

References fLength.

Referenced by CheckData(), NtpTimingNtuple::Fill(), and NtpEventDisplayTiming::MakePicture().

00104 { return fLength; }

Int_t NtpTimingFit::GetMaxPlane (  )  [inline]

Definition at line 96 of file NtpTimingFit.h.

References fMaxPlane.

Referenced by NtpTimingNtuple::Fill().

00096 { return fMaxPlane; }

Int_t NtpTimingFit::GetMinPlane (  )  [inline]

Definition at line 95 of file NtpTimingFit.h.

References fMinPlane.

Referenced by NtpTimingNtuple::Fill().

00095 { return fMinPlane; }

Double_t NtpTimingFit::GetMinTime (  )  [inline]

Definition at line 100 of file NtpTimingFit.h.

References fMinTime.

Referenced by NtpTimingNtuple::Fill().

00100 { return fMinTime; }

Double_t NtpTimingFit::GetMinZ (  )  [inline]

Definition at line 99 of file NtpTimingFit.h.

References fMinZ.

Referenced by NtpTimingNtuple::Fill().

00099 { return fMinZ; }

UInt_t NtpTimingFit::GetNumPoints (  ) 
Int_t NtpTimingFit::GetPlane ( UInt_t  ihit  ) 

Definition at line 469 of file NtpTimingFit.cxx.

References fTimePoint, GetNumPoints(), TimePoint::plane, and vecTimePoints.

Referenced by CalcSeedTime(), NtpTimingNtuple::Fill(), RunConstrainedFit(), and RunLinearFit().

00470 {
00471   if( ihit<GetNumPoints() ){
00472     fTimePoint = (TimePoint)(vecTimePoints.at(ihit));
00473     return fTimePoint.plane;
00474   }
00475 
00476   return -1;
00477 }

Double_t NtpTimingFit::GetQ ( UInt_t  ihit  ) 

Definition at line 529 of file NtpTimingFit.cxx.

References fTimePoint, GetNumPoints(), TimePoint::q, and vecTimePoints.

Referenced by NtpTimingNtuple::Fill(), GetW(), RunConstrainedFit(), and RunLinearFit().

00530 {
00531   if( ihit<GetNumPoints() ){
00532     fTimePoint = (TimePoint)(vecTimePoints.at(ihit));
00533     return fTimePoint.q;
00534   }
00535 
00536   return 0.0;
00537 }

Double_t NtpTimingFit::GetS ( UInt_t  ihit  ) 

Definition at line 509 of file NtpTimingFit.cxx.

References TimePoint::ds, fTimePoint, GetNumPoints(), and vecTimePoints.

Referenced by CalcSeedTime(), NtpTimingNtuple::Fill(), NtpEventDisplayTiming::MakePicture(), RunConstrainedFit(), and RunLinearFit().

00510 {
00511   if( ihit<GetNumPoints() ){
00512     fTimePoint = (TimePoint)(vecTimePoints.at(ihit));
00513     return fTimePoint.ds;
00514   }
00515 
00516   return 0.0;
00517 }

Int_t NtpTimingFit::GetStrip ( UInt_t  ihit  ) 

Definition at line 479 of file NtpTimingFit.cxx.

References fTimePoint, GetNumPoints(), TimePoint::strip, and vecTimePoints.

Referenced by NtpTimingNtuple::Fill().

00480 {
00481   if( ihit<GetNumPoints() ){
00482     fTimePoint = (TimePoint)(vecTimePoints.at(ihit));
00483     return fTimePoint.strip;
00484   }
00485 
00486   return -1;
00487 }

Double_t NtpTimingFit::GetT ( UInt_t  ihit  ) 

Definition at line 519 of file NtpTimingFit.cxx.

References TimePoint::dt, fTimePoint, GetNumPoints(), and vecTimePoints.

Referenced by CalcSeedTime(), NtpTimingNtuple::Fill(), NtpEventDisplayTiming::MakePicture(), RunConstrainedFit(), and RunLinearFit().

00520 {
00521   if( ihit<GetNumPoints() ){
00522     fTimePoint = (TimePoint)(vecTimePoints.at(ihit));
00523     return fTimePoint.dt;
00524   }
00525 
00526   return 0.0;
00527 }

Double_t NtpTimingFit::GetTimeConstant ( Int_t  plane,
Int_t  strip,
Int_t  end 
) [private]

Definition at line 1176 of file NtpTimingFit.cxx.

References fDetector, fUseTimeConstants, Detector::kFar, and Detector::kNear.

Referenced by AddEvent().

01177 {
01178   // Sanity Check
01179   // ============
01180   if( fUseTimeConstants==0 ) return 0.0;
01181 
01182   // Far Detector [already done]
01183   // ===========================
01184   if( fDetector==Detector::kFar ){
01185     return 0.0;
01186   }  
01187 
01188   // Near Detector
01189   // =============
01190   if( fDetector==Detector::kNear ){
01191     return NtpTimeCalibrator::GetTimeConstant(plane,strip,end);
01192   }
01193 
01194   return 0.0;
01195 }

Double_t NtpTimingFit::GetTimeWalk ( Int_t  ndigits,
Double_t  q 
) [private]

Definition at line 1197 of file NtpTimingFit.cxx.

References fDetector, fUseTimeWalk, Detector::kFar, and Detector::kNear.

Referenced by AddEvent().

01198 {  
01199   // Sanity Check
01200   // ============
01201   if( fUseTimeWalk==0 ) return 0.0;
01202 
01203   // Pulse Height
01204   // ============
01205   Double_t qpes = q;
01206   if( qpes<=1.0 ) qpes = 1.0;
01207   if( qpes>=80.0 ) qpes = 80.0;
01208 
01209   // Far Detector [already done]
01210   // ===========================
01211   if( fDetector==Detector::kFar ){
01212     return 0.0;
01213   }
01214 
01215   // Near Detector
01216   // =============
01217   if( fDetector==Detector::kNear ){
01218     return NtpTimeCalibrator::GetTimeWalk(ndigits,qpes);
01219   }
01220 
01221   return 0.0;
01222 }

Double_t NtpTimingFit::GetTimeWeight ( Double_t  q  )  [private]

Definition at line 1224 of file NtpTimingFit.cxx.

References fDetector, Detector::kFar, and Detector::kNear.

Referenced by GetW(), RunConstrainedFit(), and RunLinearFit().

01225 {   
01226   // Pulse Height
01227   // ============
01228   Double_t qpes = q;
01229   if( qpes<=0 ) return 0.0;
01230   if( qpes>25.0 ) qpes = 25.0;
01231 
01232   // Parameterisation
01233   // ================
01234   Double_t A = 0.0;
01235   Double_t B = 0.0;
01236   Double_t C = 0.0;
01237   Double_t S = 0.0;
01238 
01239   // Far Detector
01240   // ============
01241   if( fDetector==Detector::kFar ){
01242     A = +0.58;
01243     B = +0.69;
01244     C = -0.33;    
01245     S = 1.00;
01246   }
01247 
01248   // Near Detector
01249   // =============
01250   if( fDetector==Detector::kNear ){
01251     A = +1.13; 
01252     B = +0.63; 
01253     C = -0.21;
01254     S = 1.33;  // But, S = 1.10 with extra time-walk correction!
01255   }
01256 
01257   // Calculate Timing Resolution
01258   // ===========================
01259   Double_t sigma = S*(A+exp(B+C*q))/0.3;
01260 
01261   if( sigma>0.0 ) return 1.0/(sigma*sigma);
01262   else return 0.0;
01263 }

Int_t NtpTimingFit::GetTrackPlanes (  )  [inline]

Definition at line 97 of file NtpTimingFit.h.

References fTrackPlanes.

Referenced by NtpTimingNtuple::Fill().

00097 { return fTrackPlanes; }

Double_t NtpTimingFit::GetTrigTime (  )  [inline]

Definition at line 101 of file NtpTimingFit.h.

References fTrigTime.

Referenced by NtpTimingNtuple::Fill().

00101 { return fTrigTime; }

Double_t NtpTimingFit::GetW ( UInt_t  ihit  ) 

Definition at line 539 of file NtpTimingFit.cxx.

References GetQ(), and GetTimeWeight().

Referenced by NtpTimingNtuple::Fill().

00540 {
00541   Double_t q = GetQ(ihit);
00542 
00543   if( q>0.0 ){
00544     return this->GetTimeWeight(q);
00545   }
00546   else{
00547     return 0.0;
00548   }
00549 }

NtpTimingFit * NtpTimingFit::Instance ( void   )  [static]

Definition at line 22 of file NtpTimingFit.cxx.

Referenced by NtpTimingNtuple::Fill(), NtpEventDisplayTiming::MakePicture(), RunFit(), SelectBeamMode(), and SelectCosmicMode().

00023 {
00024   if( !fgTimingFit ){
00025     fgTimingFit = new NtpTimingFit();
00026   }
00027 
00028   return fgTimingFit;
00029 }

void NtpTimingFit::PrintTimePoints (  ) 

Definition at line 195 of file NtpTimingFit.cxx.

References TimePoint::digits, TimePoint::ds, TimePoint::dt, TimePoint::end, fMinTime, fTimePoint, n, TimePoint::plane, TimePoint::q, TimePoint::strip, and vecTimePoints.

Referenced by AddEvent().

00196 {
00197   std::cout << " *** Printing Time Points *** " << std::endl;
00198   
00199   for( UInt_t n=0; n<vecTimePoints.size(); n++ ){
00200     fTimePoint = (TimePoint)(vecTimePoints.at(n));
00201     std::cout << " [" << n << "] (" << fTimePoint.plane << "|" << fTimePoint.strip << "|" << fTimePoint.end << ") (" << fTimePoint.digits << ") (" << fTimePoint.ds << "," << fTimePoint.dt-fMinTime << "," << fTimePoint.q << ")" << std::endl;
00202   }
00203 }

void NtpTimingFit::PrintTimeStrips (  ) 

Definition at line 183 of file NtpTimingFit.cxx.

References TimeStrip::digits, TimeStrip::end, fMinTime, fTimeStrip, n, TimeStrip::plane, TimeStrip::q, TimeStrip::strip, TimeStrip::t, TimeStrip::u, TimeStrip::v, vecTimeStrips, TimeStrip::view, and TimeStrip::z.

Referenced by AddEvent().

00184 {
00185   std::cout << " *** Printing Time Strips *** " << std::endl;
00186   
00187   for( UInt_t n=0; n<vecTimeStrips.size(); n++ ){
00188     fTimeStrip = (TimeStrip)(vecTimeStrips.at(n));
00189     std::cout << " [" << n << "] (" << fTimeStrip.view << ") (" << fTimeStrip.plane << "|" << fTimeStrip.strip << "|" << fTimeStrip.end << ") (" << fTimeStrip.digits << ") (" << fTimeStrip.u << "," << fTimeStrip.v << "," << fTimeStrip.z << ") (" << fTimeStrip.t-fMinTime << "," << fTimeStrip.q << ") " << std::endl;
00190   }
00191 
00192   return;
00193 }

void NtpTimingFit::Reset (  ) 

Definition at line 148 of file NtpTimingFit.cxx.

References ResetTimeStrips().

Referenced by AddEvent(), NtpTimingNtuple::Fill(), NtpEventDisplayTiming::MakePicture(), and RunFit().

00149 {
00150   this->ResetTimeStrips();
00151 
00152   return;
00153 }

void NtpTimingFit::ResetTimePoints (  ) 

Definition at line 174 of file NtpTimingFit.cxx.

References fLength, and vecTimePoints.

Referenced by CalcTimePoints(), and ResetTimeStrips().

00175 {
00176   fLength = 0.0;
00177 
00178   vecTimePoints.clear();
00179 
00180   return;
00181 }

void NtpTimingFit::ResetTimeStrips (  ) 

Definition at line 155 of file NtpTimingFit.cxx.

References fCrateT0, fGateTime, fMaxPlane, fMinPlane, fMinTime, fMinZ, fTrackPlanes, fTrigTime, ResetTimePoints(), and vecTimeStrips.

Referenced by Reset().

00156 {  
00157   fMinPlane = -1;
00158   fMaxPlane = -1;
00159   fTrackPlanes = 0;
00160 
00161   fMinZ = -99999.9;
00162   fMinTime = -99999.9;
00163   fTrigTime = 0.0;
00164   fGateTime = 0.0;
00165   fCrateT0 = 0.0;
00166 
00167   vecTimeStrips.clear();
00168 
00169   this->ResetTimePoints();
00170 
00171   return;
00172 }

Bool_t NtpTimingFit::RunBeamFit ( NtpStRecord ntpRecord,
Int_t  ievent,
Double_t &  vtime,
Double_t &  err,
Double_t &  chi2 
) [static]

Definition at line 80 of file NtpTimingFit.cxx.

References RunFit(), and SelectBeamMode().

Referenced by NuExtraction::ExtractVertexTimeFit().

00081 {
00082   NtpTimingFit::SelectBeamMode();
00083 
00084   return NtpTimingFit::RunFit( ntpRecord, ievent,
00085                                vtime, err, chi2 );
00086 }

Bool_t NtpTimingFit::RunConstrainedFit ( Int_t  dir,
Double_t &  invbeta,
Double_t &  vtime,
Double_t &  err,
Double_t &  rms,
Double_t &  chi2 
)

Definition at line 897 of file NtpTimingFit.cxx.

References CalcSeedTime(), CheckData(), fDebug, fMaxPlaneWindow, fMinPlaneWindow, fTimeResolution, GetNumPoints(), GetPlane(), GetQ(), GetS(), GetT(), GetTimeWeight(), and n.

00898 { 
00899   // default values 
00900   invbeta = 0.0;
00901   vtime   = 0.0;
00902   err     = 0.0;
00903   rms     = 0.0;
00904   chi2    = 0.0;
00905 
00906   // sanity check
00907   if( dir==0 ) return 0;
00908   
00909   // more checks
00910   if( CheckData()==0 ) return 0;
00911 
00912   // run fit
00913   Double_t input_vtime  = 0.0;
00914   Double_t input_window = 10.0*fTimeResolution;
00915   
00916   Bool_t input_pass = this->CalcSeedTime( dir, input_vtime, input_window );
00917   if( input_pass==0 ) return 0;
00918 
00919   // run constrained fit
00920   if( dir>0 ){
00921     if( fDebug ) std::cout << " *** Run Constrained Fit (Forwards) *** " << std::endl;
00922   }
00923   if( dir<0 ){
00924     if( fDebug ) std::cout << " *** Run Constrained Fit (Backwards) *** " << std::endl;
00925   }
00926 
00927   Int_t plane = 0;
00928   Int_t ctr = 0;
00929 
00930   Double_t delta = 0.0;
00931 
00932   Double_t t = 0.0;
00933   Double_t s = 0.0;
00934   Double_t q = 0.0;
00935   Double_t w = 0.0;
00936 
00937   Double_t Swtt = 0.0;
00938   Double_t Swt = 0.0;
00939   Double_t Sw = 0.0;
00940   Double_t S = 0.0;
00941   Double_t Stot = 0.0;
00942 
00943   Double_t t0    = input_vtime;
00944   Double_t tfit  = 0.0;
00945   Double_t terr  = 0.0;
00946   Double_t sigma = -99999.9;
00947   Double_t chisq = -99999.9;
00948   Double_t var   = 0.0; 
00949   Double_t frac  = 0.0;
00950   
00951   Double_t window = input_window;
00952 
00953   Double_t ibeta = 0.0;
00954   if( dir>0 ) ibeta = +1.0;
00955   if( dir<0 ) ibeta = -1.0;
00956 
00957   Double_t dof = (double)(GetNumPoints()-1.0);
00958 
00959   for( UInt_t n=0; n<GetNumPoints(); n++ ){
00960     plane = GetPlane(n);
00961     
00962     if( plane>fMinPlaneWindow 
00963      && plane<fMaxPlaneWindow ){
00964       s = GetS(n);
00965       t = GetT(n) - t0;
00966       q = GetQ(n);
00967       w = GetTimeWeight(q);
00968 
00969       delta = t - (ibeta/0.3)*s;
00970 
00971       if( fDebug ){
00972         std::cout << " [" << n << "]\t s=" << s << "\t t=" << t << "\t delta=" << delta << std::endl;
00973       }
00974 
00975       Stot += 1.0;
00976       if( fabs(delta)<window ){
00977         Swtt += w*delta*delta; 
00978         Swt  += w*delta; 
00979         Sw   += w; 
00980         S    += 1.0; 
00981         ctr++;
00982       }
00983     }
00984   }
00985 
00986   if( S>1.0 ){
00987     if( Sw>0.0 ){
00988       tfit  = Swt/Sw; 
00989       terr  = 1.0/sqrt(Sw);
00990       var   = (Swtt/Sw)-(Swt/Sw)*(Swt/Sw);
00991       dof   = S - 1.0;
00992       
00993       if( var>0.0 ){
00994         frac  = S/Stot;
00995         sigma = sqrt(var);
00996         chisq = Sw*var/dof;
00997       }
00998     } 
00999   }
01000        
01001   // return parameters  
01002   invbeta = ibeta;
01003   vtime   = t0 + tfit;
01004   err     = terr;
01005   rms     = sigma;
01006   chi2    = chisq;
01007 
01008   if( fDebug ) std::cout << "   Results: 1/beta=" << invbeta << " vtime=" << vtime << " err=" << err << " rms=" << rms << " chi2=" << chi2 << " pts=" << ctr << "/" << GetNumPoints() << std::endl;
01009 
01010   if( frac>0.0 ) return 1;  else return 0;
01011 }

Bool_t NtpTimingFit::RunConstrainedFit ( Double_t &  invbeta,
Double_t &  vtime,
Double_t &  err,
Double_t &  rms,
Double_t &  chi2 
)

Definition at line 697 of file NtpTimingFit.cxx.

References fUseBeamMode, and RunConstrainedFit().

00698 {
00699   Double_t invbeta_fwd = 0.0;
00700   Double_t vtime_fwd   = 0.0;
00701   Double_t err_fwd     = 0.0;
00702   Double_t rms_fwd     = 0.0;
00703   Double_t chi2_fwd    = 0.0;
00704   Bool_t pass_fwd = 0;
00705 
00706   Double_t invbeta_bwd = 0.0;
00707   Double_t vtime_bwd   = 0.0;
00708   Double_t err_bwd     = 0.0;
00709   Double_t rms_bwd     = 0.0;
00710   Double_t chi2_bwd    = 0.0;
00711   Bool_t pass_bwd = 0;
00712 
00713   pass_fwd = this->RunConstrainedFit( +1, invbeta_fwd, vtime_fwd, err_fwd, rms_fwd, chi2_fwd );
00714 
00715   if( fUseBeamMode ){
00716     if( pass_fwd ){
00717       invbeta = invbeta_fwd; vtime = vtime_fwd; err = err_fwd; rms = rms_fwd; chi2 = chi2_fwd;
00718       return 1;
00719     }
00720     else{
00721       return 0;
00722     }
00723   }
00724 
00725   pass_bwd = this->RunConstrainedFit( -1, invbeta_bwd, vtime_bwd, err_bwd, rms_bwd, chi2_bwd );
00726       
00727   if( pass_fwd && pass_bwd ){
00728     if( chi2_fwd<=chi2_bwd ){
00729       invbeta = invbeta_fwd; vtime = vtime_fwd; err = err_fwd; rms = rms_fwd; chi2 = chi2_fwd;
00730       return 1;
00731     }
00732     else {
00733       invbeta = invbeta_bwd; vtime = vtime_bwd; err = err_bwd; rms = rms_bwd; chi2 = chi2_bwd;
00734       return 1;
00735     }
00736   }
00737 
00738   return 0;
00739 }

Bool_t NtpTimingFit::RunConstrainedFit ( Double_t &  invbeta,
Double_t &  vtime 
)

Definition at line 687 of file NtpTimingFit.cxx.

References err().

Referenced by NtpTimingNtuple::Fill(), NtpEventDisplayTiming::MakePicture(), RunConstrainedFit(), RunFit(), and RunLinearFit().

00688 {
00689   Double_t err  = 0.0;
00690   Double_t rms  = 0.0;
00691   Double_t chi2 = 0.0;
00692 
00693   return this->RunConstrainedFit(invbeta, vtime,
00694                                  err, rms, chi2);
00695 }

Bool_t NtpTimingFit::RunCosmicFit ( NtpStRecord ntpRecord,
Int_t  ievent,
Double_t &  vtime,
Double_t &  err,
Double_t &  chi2 
) [static]

Definition at line 88 of file NtpTimingFit.cxx.

References RunFit(), and SelectCosmicMode().

00089 {
00090   NtpTimingFit::SelectCosmicMode();
00091 
00092   return NtpTimingFit::RunFit( ntpRecord, ievent,
00093                                vtime, err, chi2 );
00094 }

Bool_t NtpTimingFit::RunFit ( NtpStRecord ntpRecord,
Int_t  ievent,
Double_t &  vtime,
Double_t &  err,
Double_t &  chi2 
) [static]

Definition at line 96 of file NtpTimingFit.cxx.

References AddEvent(), Instance(), passfail(), Reset(), and RunConstrainedFit().

Referenced by RunBeamFit(), and RunCosmicFit().

00097 {
00098   NtpTimingFit::Instance()->Reset();
00099 
00100   NtpTimingFit::Instance()->AddEvent(ntpRecord,ievent);
00101 
00102   Double_t output_invbeta = 0.0;
00103   Double_t output_vtime   = 0.0; 
00104   Double_t output_err     = 0.0;
00105   Double_t output_rms     = 0.0;
00106   Double_t output_chi2    = 0.0;
00107     
00108   Bool_t passfail = 0;
00109 
00110   passfail = NtpTimingFit::Instance()->RunConstrainedFit( output_invbeta, output_vtime, output_err,
00111                                                           output_rms, output_chi2 );
00112   if( passfail ){
00113     vtime = 1.0e-9*output_vtime;  // nsec -> sec
00114     err   = 1.0e-9*output_err;    // nsec -> sec
00115     chi2  =        output_chi2;
00116   }
00117   
00118   return passfail;
00119 }

Bool_t NtpTimingFit::RunLinearFit ( Double_t &  invbeta,
Double_t &  vtime,
Double_t &  err,
Double_t &  rms,
Double_t &  chi2 
)

Definition at line 1023 of file NtpTimingFit.cxx.

References CheckData(), fDebug, fMaxPlaneWindow, fMinPlaneWindow, GetNumPoints(), GetPlane(), GetQ(), GetS(), GetT(), GetTimeWeight(), n, and RunConstrainedFit().

01024 {    
01025   // default values 
01026   invbeta = 0.0;
01027   vtime   = 0.0;
01028   err     = 0.0;
01029   rms     = 0.0;
01030   chi2    = 0.0;
01031 
01032   // more checks
01033   if( CheckData()==0 ) return 0;
01034 
01035   // run fit
01036   Double_t input_invbeta = 0.0;
01037   Double_t input_offset  = 0.0;
01038   Double_t input_err     = 0.0;
01039   Double_t input_rms     = 0.0;
01040   Double_t input_chi2    = 0.0;
01041 
01042   Bool_t input_pass = this->RunConstrainedFit( input_invbeta, input_offset, input_err, input_rms, input_chi2 );
01043   if( input_pass==0 ) return 0;
01044 
01045   Double_t slope  = input_invbeta/0.3;
01046   Double_t offset = input_offset;
01047 
01048   if( fDebug ) std::cout << " *** Run Linear Fit *** " << std::endl;  
01049 
01050   Int_t plane = 0;
01051   Int_t ctr = 0;
01052 
01053   Double_t delta0 = 0.0;
01054   Double_t delta = 0.0;
01055 
01056   Double_t t = 0.0;
01057   Double_t s = 0.0;
01058   Double_t q = 0.0;
01059   Double_t w = 0.0;
01060 
01061   Double_t Sws  = 0.0; 
01062   Double_t Swt  = 0.0; 
01063   Double_t Swss = 0.0; 
01064   Double_t Swst = 0.0; 
01065   Double_t Swtt = 0.0; 
01066   Double_t Sw   = 0.0;
01067   Double_t S    = 0.0; 
01068 
01069   Double_t Qwxx = 0.0;
01070   Double_t Qwx  = 0.0;
01071   Double_t Qw   = 0.0;
01072   Double_t Q    = 0.0;
01073   Double_t Qtot = 0.0;
01074 
01075   Double_t t0 = input_offset;
01076 
01077   Double_t error  = 0.0;
01078   Double_t sigma = -99999.9;
01079   Double_t chisq = -99999.9;
01080   Double_t var   = 0.0; 
01081   Double_t frac  = 0.0;
01082 
01083   Double_t window = 2.0*input_rms;
01084 
01085   Double_t dof = (double)(GetNumPoints()-2.0);
01086 
01087   // perform linear regression
01088   for( UInt_t n=0; n<GetNumPoints(); n++ ){     
01089     plane = GetPlane(n);
01090     
01091     if( plane>fMinPlaneWindow 
01092      && plane<fMaxPlaneWindow ){
01093       s = GetS(n);
01094       t = GetT(n) - t0;
01095       q = GetQ(n);
01096       w = GetTimeWeight(q);
01097 
01098       delta0 = t - slope*s;
01099 
01100       if( fDebug ){
01101         std::cout << " [" << n << "]\t s=" << s << "\t t=" << t << "\t delta=" << delta0 << std::endl;
01102       }
01103 
01104       if( fabs(delta0)<window ){
01105         Sws  += w*s; 
01106         Swt  += w*t; 
01107         Swss += w*s*s; 
01108         Swst += w*s*t; 
01109         Swtt += w*t*t; 
01110         Sw   += w; 
01111         S    += 1.0;
01112       }
01113     }
01114   }
01115 
01116   if( S>2.0 ){
01117     if( Sw>0.0 ){
01118       if( Sw*Swss-Sws*Sws !=0.0 ){
01119         slope  = (Sw*Swst-Sws*Swt)/(Sw*Swss-Sws*Sws);
01120         offset = (Swt*Swss-Sws*Swst)/(Sw*Swss-Sws*Sws);
01121       }
01122     }
01123   }
01124 
01125   // calculate residuals
01126   for( UInt_t n=0; n<GetNumPoints(); n++ ){     
01127     plane = GetPlane(n);
01128     
01129     if( plane>fMinPlaneWindow 
01130      && plane<fMaxPlaneWindow ){
01131       s = GetS(n);
01132       t = GetT(n) - t0;
01133       q = GetQ(n);
01134       w = GetTimeWeight(q);
01135 
01136       delta0 = t - slope*s;
01137       delta  = (t-offset) - slope*s;
01138 
01139       Qtot += 1.0;
01140       if( fabs(delta0)<window ){
01141         Qwxx += w*delta*delta;
01142         Qwx  += w*delta;
01143         Qw   += w;
01144         Q    += 1.0;
01145         ctr++;
01146       }
01147     }
01148   }
01149 
01150   if( Q>2.0 ){
01151     if( Qw>0.0 ){
01152       var = (Qwxx/Qw)-(Qwx/Qw)*(Qwx/Qw);
01153       dof = Q - 2.0;
01154 
01155       if( var>0.0 ){
01156         frac  = Q/Qtot;
01157         sigma = sqrt(var);
01158         error = sqrt(var/dof); // note: this is an approximation
01159         chisq = Qw*var/dof;
01160       }
01161     }
01162   }
01163 
01164   // return parameters
01165   invbeta = 0.3*slope;
01166   vtime   = t0 + offset;
01167   err     = error;
01168   rms     = sigma;
01169   chi2    = chisq;
01170 
01171   if( fDebug ) std::cout << "   Results: 1/beta=" << invbeta << " vtime=" << vtime << " err=" << err << " rms=" << rms << " chi2=" << chi2 << " pts=" << ctr << "/" << GetNumPoints() << std::endl;
01172 
01173   if( frac>0.0 ) return 1;  else return 0;
01174 }

Bool_t NtpTimingFit::RunLinearFit ( Double_t &  invbeta,
Double_t &  vtime 
)

Definition at line 1013 of file NtpTimingFit.cxx.

References err().

Referenced by NtpTimingNtuple::Fill(), and NtpEventDisplayTiming::MakePicture().

01014 {
01015   Double_t err  = 0.0;
01016   Double_t rms  = 0.0;
01017   Double_t chi2 = 0.0; 
01018 
01019   return this->RunLinearFit(invbeta, vtime,
01020                            err, rms, chi2);
01021 }

void NtpTimingFit::SelectBeamMode (  )  [static]

Definition at line 70 of file NtpTimingFit.cxx.

References Instance(), and UseBeamMode().

Referenced by RunBeamFit().

00071 {
00072   NtpTimingFit::Instance()->UseBeamMode();
00073 }

void NtpTimingFit::SelectCosmicMode (  )  [static]

Definition at line 75 of file NtpTimingFit.cxx.

References Instance(), and UseCosmicMode().

Referenced by RunCosmicFit().

00076 {
00077   NtpTimingFit::Instance()->UseCosmicMode();
00078 } 

void NtpTimingFit::SetDebug ( Bool_t  onoff = 1  )  [inline]

Definition at line 117 of file NtpTimingFit.h.

References fDebug.

00117 { fDebug = onoff; }

void NtpTimingFit::SetDetector ( Detector::Detector_t  myDetector  ) 

Definition at line 121 of file NtpTimingFit.cxx.

References fDetector, fMaxPlaneWindow, fMinPlaneWindow, fTimeResolution, fUseSpectrometer, Detector::kFar, and Detector::kNear.

Referenced by AddEvent().

00122 {
00123   if( fDetector != myDetector ){
00124     fDetector = myDetector;
00125     std::cout << "   Setting Detector: " << fDetector << std::endl;
00126 
00127     fTimeResolution = 50.0;
00128 
00129     if( fDetector==Detector::kNear ){
00130       fTimeResolution = 15.0;
00131       fMinPlaneWindow = 0;
00132       if( fUseSpectrometer ){
00133         fMaxPlaneWindow = 300;
00134       }
00135       else{
00136         fMaxPlaneWindow = 125;
00137       }
00138     }
00139 
00140     if( fDetector==Detector::kFar ){
00141       fTimeResolution = 2.5; 
00142       fMinPlaneWindow = 0;
00143       fMaxPlaneWindow = 500;
00144     }
00145   }
00146 }

void NtpTimingFit::UseBeamMode (  )  [inline]

Definition at line 66 of file NtpTimingFit.h.

References fFoundMode, and fUseBeamMode.

Referenced by SelectBeamMode().

00066 { fUseBeamMode = 1;  fFoundMode = 1; }

void NtpTimingFit::UseCalculator ( Bool_t  yesno = 1  )  [inline]

Definition at line 64 of file NtpTimingFit.h.

References fUseCalculator.

00064 { fUseCalculator    = yesno; }

void NtpTimingFit::UseCosmicMode (  )  [inline]

Definition at line 67 of file NtpTimingFit.h.

References fFoundMode, and fUseBeamMode.

Referenced by SelectCosmicMode().

00067 { fUseBeamMode = 0;  fFoundMode = 1; }

void NtpTimingFit::UseSpectrometer ( Bool_t  yesno = 1  )  [inline]

Definition at line 63 of file NtpTimingFit.h.

References fUseSpectrometer.

00063 { fUseSpectrometer  = yesno; }

void NtpTimingFit::UseTimeConstants ( Bool_t  yesno = 1  )  [inline]

Definition at line 61 of file NtpTimingFit.h.

References fUseTimeConstants.

00061 { fUseTimeConstants = yesno; }

void NtpTimingFit::UseTimeWalk ( Bool_t  yesno = 1  )  [inline]

Definition at line 60 of file NtpTimingFit.h.

References fUseTimeWalk.

00060 { fUseTimeWalk      = yesno; }


Member Data Documentation

Double_t NtpTimingFit::fCrateT0 [private]

Definition at line 151 of file NtpTimingFit.h.

Referenced by AddEvent(), NtpTimingFit(), and ResetTimeStrips().

Bool_t NtpTimingFit::fDebug [private]
Double_t NtpTimingFit::fDelta [private]

Definition at line 153 of file NtpTimingFit.h.

Referenced by CalcSeedTime().

Int_t NtpTimingFit::fDetector [private]
Bool_t NtpTimingFit::fFoundMode [private]

Definition at line 156 of file NtpTimingFit.h.

Referenced by CheckData(), NtpTimingFit(), UseBeamMode(), and UseCosmicMode().

Double_t NtpTimingFit::fGateTime [private]

Definition at line 150 of file NtpTimingFit.h.

Referenced by AddEvent(), GetGateTime(), and ResetTimeStrips().

Double_t NtpTimingFit::fLength [private]

Definition at line 152 of file NtpTimingFit.h.

Referenced by AddData(), GetLength(), NtpTimingFit(), and ResetTimePoints().

Int_t NtpTimingFit::fMaxPlane [private]

Definition at line 144 of file NtpTimingFit.h.

Referenced by AddData(), GetMaxPlane(), NtpTimingFit(), and ResetTimeStrips().

Int_t NtpTimingFit::fMinPlane [private]

Definition at line 143 of file NtpTimingFit.h.

Referenced by AddData(), GetMinPlane(), NtpTimingFit(), and ResetTimeStrips().

Double_t NtpTimingFit::fMinTime [private]
Double_t NtpTimingFit::fMinZ [private]

Definition at line 147 of file NtpTimingFit.h.

Referenced by AddData(), AddEvent(), GetMinZ(), NtpTimingFit(), and ResetTimeStrips().

Definition at line 133 of file NtpTimingFit.h.

Referenced by AddData(), GetDigits(), GetEnd(), GetPlane(), GetQ(), GetS(), GetStrip(), GetT(), and PrintTimePoints().

Double_t NtpTimingFit::fTimeResolution [private]

Definition at line 127 of file NtpTimingFit.h.

Referenced by CalcSeedTime(), NtpTimingFit(), RunConstrainedFit(), and SetDetector().

Definition at line 132 of file NtpTimingFit.h.

Referenced by AddData(), CalcTimePoints(), and PrintTimeStrips().

Int_t NtpTimingFit::fTrackPlanes [private]

Definition at line 145 of file NtpTimingFit.h.

Referenced by AddEvent(), GetTrackPlanes(), NtpTimingFit(), and ResetTimeStrips().

Double_t NtpTimingFit::fTrigTime [private]

Definition at line 149 of file NtpTimingFit.h.

Referenced by AddEvent(), GetTrigTime(), NtpTimingFit(), and ResetTimeStrips().

Bool_t NtpTimingFit::fUseBeamMode [private]

Definition at line 158 of file NtpTimingFit.h.

Referenced by NtpTimingFit(), RunConstrainedFit(), UseBeamMode(), and UseCosmicMode().

Bool_t NtpTimingFit::fUseCalculator [private]

Definition at line 162 of file NtpTimingFit.h.

Referenced by AddEvent(), NtpTimingFit(), and UseCalculator().

Definition at line 161 of file NtpTimingFit.h.

Referenced by NtpTimingFit(), SetDetector(), and UseSpectrometer().

Definition at line 160 of file NtpTimingFit.h.

Referenced by AddEvent(), GetTimeConstant(), NtpTimingFit(), and UseTimeConstants().

Bool_t NtpTimingFit::fUseTimeWalk [private]

Definition at line 159 of file NtpTimingFit.h.

Referenced by AddEvent(), GetTimeWalk(), NtpTimingFit(), and UseTimeWalk().

std::vector<Double_t> NtpTimingFit::vecResiduals [private]

Definition at line 138 of file NtpTimingFit.h.

Referenced by CalcSeedTime().

std::vector<Double_t> NtpTimingFit::vecResidualsForFit [private]

Definition at line 139 of file NtpTimingFit.h.

Referenced by CalcSeedTime().

std::vector<TimePoint> NtpTimingFit::vecTimePoints [private]
std::vector<TimeStrip> NtpTimingFit::vecTimeStrips [private]

Definition at line 135 of file NtpTimingFit.h.

Referenced by AddData(), CalcTimePoints(), PrintTimeStrips(), and ResetTimeStrips().


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

Generated on 22 Nov 2017 for loon by  doxygen 1.6.1