#include #include #include #include #include #include #include #include #include #include #include using namespace RAT; using namespace RAT::PMTSelectors; #include #include #include using namespace std; void TimeResidualCut::Initialise( const std::string& param ) { // Default values fLowCut = -99999.0; fHighCut = 99999.0; fUseGroup = 0; fUseRefracted = false; } void TimeResidualCut::BeginOfRun(DS::Run& ) { fLightPath = DU::Utility::Get()->GetLightPathCalculator(); // Values for Isotropy classifier DB *db = DB::Get(); string index; if( fLowCut == -99999.0 || fHighCut == 99999.0 ){ try{ index = db->GetLink("GEO","inner_av")->GetS("material"); DBLinkPtr dbLink = db->GetLink("SELECTOR_ISOTROPY",index); fLowCut = dbLink->GetD("fLowCut"); fHighCut = dbLink->GetD("fHighCut"); DBLinkGroup grp = db->GetLinkGroup("SELECTOR_ISOTROPY"); DBLinkGroup::iterator it; for(it = grp.begin(); it != grp.end(); it++) { if(index == it->first) { dbLink = db->GetLink("SELECTOR_ISOTROPY",index); fLowCut = dbLink->GetD("fLowCut"); fHighCut = dbLink->GetD("fHighCut"); } } } catch(DBNotFoundError & ){ warn << "TimeResidualCut::BeginOfRun: No full definition for " + index + " using defaults instead.\n"; warn << "TimeResidualCut::BeginOfRun: Using fLowCut = " << fLowCut << ", fHighCut = " << fHighCut << "\n"; } } } void TimeResidualCut::SetI( const std::string& param, int value ) { if( param == string( "useGroup" ) ) { if( value == 0 ) fUseGroup = false; else fUseGroup = true; } else if( param == string( "useRefracted" ) ) { if( value == 0 ) fUseRefracted = false; else fUseRefracted = true; } else throw Processor::ParamUnknown( param ); } void TimeResidualCut::SetD( const std::string& param, double value ) { if( param == string( "lowLimit" ) ) fLowCut = value; else if( param == string( "highLimit" ) ) fHighCut = value; else throw Processor::ParamUnknown( param ); } double TimeResidualCut::CalcTransitTime( const TVector3& fitPosition, const TVector3& pmtPosition) { if(fUseRefracted) fLightPath.CalcByPosition( fitPosition, pmtPosition, 3.103125 * 1e-6, 10.0 ); else fLightPath.CalcByPosition( fitPosition, pmtPosition ); double distInInnerAV = fLightPath.GetDistInInnerAV(); double distInAV = fLightPath.GetDistInAV(); double distInWater = fLightPath.GetDistInWater(); double transitTime; if( fUseGroup ) transitTime = DU::Utility::Get()->GetGroupVelocity().CalcByDistance( distInInnerAV, distInAV, distInWater ); else transitTime = DU::Utility::Get()->GetEffectiveVelocity().CalcByDistance( distInInnerAV, distInAV, distInWater ); return transitTime; } vector TimeResidualCut::GetSelectedPMTs( const std::vector& data, const DS::FitVertex& vertex ) { // Get the fitted position // Will throw a DataNotFound error if no information is available const TVector3 fitPosition = vertex.GetPosition(); // TODO: should this work for invalid positions? const double fitTime = vertex.GetTime(); vector selectedPMTs; for( vector::const_iterator iPMT = data.begin(); iPMT != data.end(); ++iPMT ) { const double hitTime = iPMT->GetTime(); const TVector3 pmtPosition = DU::Utility::Get()->GetPMTInfo().GetPosition( iPMT->GetID() ); double transitTime = CalcTransitTime( fitPosition, pmtPosition); const double residual = hitTime - transitTime - fitTime; if( residual < fHighCut && residual > fLowCut ) selectedPMTs.push_back( *iPMT ); } return selectedPMTs; }