#include #include #include #include #include #include #include #include using namespace RAT; using namespace RAT::Methods; using namespace RAT::DS; #include #include using namespace std; #include void PositionTimeDirectionLikelihood::Initialise( const std::string& ) { fPolarDirection = false; } void PositionTimeDirectionLikelihood::SetI( const std::string& param, const int value ) { if( param=="polarDirection" ) { if( value>0 ) fPolarDirection = true; else fPolarDirection = false; } else throw Processor::ParamUnknown( param ); } void PositionTimeDirectionLikelihood::DefaultSeed() { FitVertex vertex; // Initialise the SeedResult at the centre and with arbitrary // 3000 mm errors in all directions (note units are mm and ns) const TVector3 seedPosition( 0.0, 0.0, 0.0 ); const ROOT::Math::XYZVectorF seedPositionError( 3000.0, 3000.0, 3000.0 ); // 230 ns is roughly the peak from ET1D and GV1D PDFs const double seedTime = 230.0; const double seedTimeError = 100.0; // Take the direction seed as +ve x direction unit vector and arbitrary error. const TVector3 seedDirection( 1.0, 0.0, 0.0 ); const TVector3 seedDirectionError( 1.0, 1.0, 1.0 ); vertex.SetPosition( seedPosition, false, true ); vertex.SetPositionErrors( seedPositionError, false, true ); vertex.SetTime( seedTime, false, true ); vertex.SetTimeErrors( seedTimeError, false, true ); vertex.SetDirection( seedDirection, false, true ); vertex.SetDirectionErrors( seedDirectionError, false, true ); fSeedResult.SetVertex( 0, vertex ); } FitResult PositionTimeDirectionLikelihood::GetBestFit() { fFitResult.Reset(); if( fPMTData.empty() ) return fFitResult; CopySeedToResult(); SelectPMTData( fFitResult.GetVertex( 0 ) ); fFitResult.SetFOM( "SelectedNHit", static_cast( fSelectedPMTData.size() ) ); if( fSelectedPMTData.empty() ) throw MethodFailError( "PositionTimeDirectionLikelihood: No hits to fit" ); fOptimiser->SetComponent( this ); // Call the designated optimiser to fit and save the fom const double logL = fOptimiser->Maximise(); fFitResult.SetFOM( "LogL", logL ); // Now save the optimised values SetParams( fOptimiser->GetParams() ); SetPositiveErrors( fOptimiser->GetPositiveErrors() ); SetNegativeErrors( fOptimiser->GetNegativeErrors() ); // The direction should be a unit vector and the errors // should be scaled accordingly. // This is not performed in SetParams as the optimiser adjusts and // uses via that function; the method should not force the optimisers decisions. FitVertex vertex = fFitResult.GetVertex( 0 ); double scaling = vertex.GetDirection().Mag(); bool validity = fOptimiser->GetValid(); if( scaling == 0.0 ) { // Cannot have a direction with no magnitude // Leave as 0, 0, 0 but set as invalid validity = false; // Ensure we don't divide by 0 scaling = 1.0; } vertex.SetDirection( vertex.GetDirection().Unit(), validity ); vertex.SetPositiveDirectionError( TVector3( vertex.GetPositiveDirectionError().X() / scaling, vertex.GetPositiveDirectionError().Y() / scaling, vertex.GetPositiveDirectionError().Z() / scaling ), validity ); vertex.SetNegativeDirectionError( TVector3( vertex.GetNegativeDirectionError().X() / scaling, vertex.GetNegativeDirectionError().Y() / scaling, vertex.GetNegativeDirectionError().Z() / scaling ), validity ); fFitResult.SetVertex( 0, vertex ); return fFitResult; } double PositionTimeDirectionLikelihood::operator()( const std::vector& params ) { SetParams(params); double logLike = 0.0; // loglike sum from pdf for( vector::const_iterator iPMT = fSelectedPMTData.begin(); iPMT != fSelectedPMTData.end(); ++iPMT ) logLike += log( fPDF->GetProbability( *iPMT, fFitResult.GetVertex( 0 ) ) ); return logLike; } vector PositionTimeDirectionLikelihood::GetParams() const { vector params; FitVertex vertex = fFitResult.GetVertex(0); params.push_back( vertex.GetPosition().x() ); params.push_back( vertex.GetPosition().y() ); params.push_back( vertex.GetPosition().z() ); params.push_back( vertex.GetTime() ); if( fPolarDirection == true ) { TVector3 direction = vertex.GetDirection().Unit(); // should always be a unit vector params.push_back( direction.Theta() ); params.push_back( direction.Phi() ); } else { params.push_back( vertex.GetDirection().x() ); params.push_back( vertex.GetDirection().y() ); params.push_back( vertex.GetDirection().z() ); } return params; } vector PositionTimeDirectionLikelihood::GetPositiveErrors() const { vector errors; FitVertex vertex = fFitResult.GetVertex(0); errors.push_back( vertex.GetPositivePositionError().x() ); errors.push_back( vertex.GetPositivePositionError().y() ); errors.push_back( vertex.GetPositivePositionError().z() ); errors.push_back( vertex.GetPositiveTimeError() ); if( fPolarDirection == true ) { errors.push_back( vertex.GetPositiveDirectionError().Theta() ); errors.push_back( vertex.GetPositiveDirectionError().Phi() ); } else { errors.push_back( vertex.GetPositiveDirectionError().x() ); errors.push_back( vertex.GetPositiveDirectionError().y() ); errors.push_back( vertex.GetPositiveDirectionError().z() ); } return errors; } vector PositionTimeDirectionLikelihood::GetNegativeErrors() const { vector errors; FitVertex vertex = fFitResult.GetVertex(0); errors.push_back( vertex.GetNegativePositionError().x() ); errors.push_back( vertex.GetNegativePositionError().y() ); errors.push_back( vertex.GetNegativePositionError().z() ); errors.push_back( vertex.GetNegativeTimeError() ); if( fPolarDirection == true ) { errors.push_back( vertex.GetNegativeDirectionError().Theta() ); errors.push_back( vertex.GetNegativeDirectionError().Phi() ); } else { errors.push_back( vertex.GetNegativeDirectionError().x() ); errors.push_back( vertex.GetNegativeDirectionError().y() ); errors.push_back( vertex.GetNegativeDirectionError().z() ); } return errors; } void PositionTimeDirectionLikelihood::SetParams( const std::vector& params ) { FitVertex vertex = fFitResult.GetVertex(0); vertex.SetPosition( TVector3( params[0], params[1], params[2] ), fOptimiser->GetValid() ); vertex.SetTime( params[3], fOptimiser->GetValid() ); if( fPolarDirection == true ) { TVector3 dir; // unit vector for the direction dir.SetMagThetaPhi(1.0, params[4], params[5]); vertex.SetDirection( dir, fOptimiser->GetValid() ); } else vertex.SetDirection( TVector3( params[4], params[5], params[6] ), fOptimiser->GetValid() ); fFitResult.SetVertex( 0, vertex ); } void PositionTimeDirectionLikelihood::SetPositiveErrors( const std::vector& errors ) { FitVertex vertex = fFitResult.GetVertex(0); vertex.SetPositivePositionError( ROOT::Math::XYZVectorF(errors[0], errors[1], errors[2]), fOptimiser->GetValid() ); vertex.SetPositiveTimeError( errors[3], fOptimiser->GetValid() ); if( fPolarDirection == true ) { TVector3 dir; // unit vector for the direction dir.SetMagThetaPhi(1.0, errors[4], errors[5]); vertex.SetPositiveDirectionError( dir, fOptimiser->GetValid() ); } else vertex.SetPositiveDirectionError( TVector3(errors[4], errors[5], errors[6]) ); fFitResult.SetVertex( 0, vertex ); } void PositionTimeDirectionLikelihood::SetNegativeErrors( const std::vector& errors ) { FitVertex vertex = fFitResult.GetVertex(0); vertex.SetNegativePositionError( ROOT::Math::XYZVectorF(errors[0], errors[1], errors[2]), fOptimiser->GetValid() ); vertex.SetNegativeTimeError( errors[3], fOptimiser->GetValid() ); if( fPolarDirection == true ) { TVector3 dir; dir.SetMagThetaPhi(1.0, errors[4], errors[5]); vertex.SetNegativeDirectionError( dir, fOptimiser->GetValid() ); } else vertex.SetNegativeDirectionError( TVector3(errors[4], errors[5], errors[6]) ); fFitResult.SetVertex( 0, vertex ); }