#include #include #include #include #include #include #include using namespace RAT; using namespace RAT::Classifiers; using namespace RAT::DS; #include using namespace ROOT; #include #include using namespace std; void BiPoCumulTimeResid::Initialise(const std::string& param) { fIndex = param; } void BiPoCumulTimeResid::BeginOfRun(DS::Run&) { DB* db = DB::Get(); string materialIndex; string defaultIndex = "te_diol_dda_0p5_labppo_scintillator_bisMSB_Jan2018"; string usedIndex = defaultIndex; if( fIndex.empty() == true ) { try { materialIndex = db->GetLink("GEO", "inner_av")->GetS("material"); } catch (DBNotFoundError & ){ try{ // If we are in a partial-fill geometry, use the upper // material materialIndex = db->GetLink("GEO","inner_av")->GetS("material_top"); } catch(DBNotFoundError & ) { Log::Die("BiPoCumulTimeResid: inner_av material and material_top " "undefined. Try something like:\n" "/rat/db/set GEO[inner_av] material " "\"te_labppo_scintillator\" or " "/rat/db/set GEO[inner_av] material_top " "\"te_labppo_scintillator\""); } } } else { materialIndex = fIndex; } // Get entry for the given material index DBLinkPtr dbLink = db->GetLink("CLASSIFIER_BIPO_CUMULTIMERESID", defaultIndex); DBLinkGroup grp = db->GetLinkGroup("CLASSIFIER_BIPO_CUMULTIMERESID"); DBLinkGroup::iterator it; for(it=grp.begin();it!=grp.end();++it) { if(materialIndex == it->first) { dbLink = db->GetLink("CLASSIFIER_BIPO_CUMULTIMERESID", materialIndex); usedIndex = materialIndex; } } if( usedIndex != materialIndex ) warn << "BiPoCumulTimeResid::No PDF available for " << materialIndex << ", using the default entry (" << defaultIndex << ")" << newline; fMinTimeResidual = dbLink->GetD("min_time_residual"); fMaxTimeResidual = dbLink->GetD("max_time_residual"); fCumulativeFractions = dbLink->GetDArray("cumulative_fractions"); fLightPath = DU::Utility::Get()->GetLightPathCalculator(); } void BiPoCumulTimeResid::DefaultSeed() { fEventPosition = TVector3(); fEventTime = 0.0; } void BiPoCumulTimeResid::SetSeed(const DS::FitResult& seed) { FitVertex seedVertex; try {seedVertex = seed.GetVertex(0);} catch( FitResult::NoVertexError& error ) {return;} try {fEventPosition = seedVertex.GetPosition();} catch( FitVertex::NoValueError& error ) {return;} try {fEventTime = seedVertex.GetTime();} catch( FitVertex::NoValueError& error ) {return;} } DS::ClassifierResult BiPoCumulTimeResid::GetClassification() { fClassifierResult.Reset(); if( fPMTData.empty() ) return fClassifierResult; // Plot the Time Residual distribution over all the PMTs in the Event // Make sure that the number of bins in this distribution is the same as the number of entries in the time residuals above TH1D* timeResidPlot = new TH1D("timeResidPlot", "timeResidPlot", fCumulativeFractions.size(), fMinTimeResidual, fMaxTimeResidual); const int eventNhits = fPMTData.size(); for (int j = 0; j < eventNhits; j++) { const TVector3 pmtPosition = DU::Utility::Get()->GetPMTInfo().GetPosition(fPMTData[j].GetID()); const double pmtTime = fPMTData[j].GetTime(); fLightPath.CalcByPosition(fEventPosition, pmtPosition); double distInInnerAV = fLightPath.GetDistInInnerAV(); double distInAV = fLightPath.GetDistInAV(); double distInWater = fLightPath.GetDistInWater(); const double flightTime = DU::Utility::Get()->GetEffectiveVelocity().CalcByDistance(distInInnerAV, distInAV, distInWater); const double timeResid = pmtTime - flightTime - fEventTime; timeResidPlot->Fill(timeResid); } // Normalise the timeResidPlot timeResidPlot->Scale(1.0 / static_cast(eventNhits)); // Fill a vector of Cumulative Time Residuals Fractions analogous to the initialised vector above vector eventCumulFracts; for (int k = 1; k <= timeResidPlot->GetNbinsX(); k++) {eventCumulFracts.push_back(timeResidPlot->Integral(0, k));} delete timeResidPlot; // Calculate the GAMMA Value - the normalised sum of the squares of the per-bin difference between the initialised and event Cumulative Time Residual Fractions double gamma = 0.0; for (unsigned int l = 0; l < fCumulativeFractions.size(); l++) { double deltaValue = fCumulativeFractions[l] - eventCumulFracts[l]; gamma += (deltaValue * deltaValue); } gamma *= (1 / (double)eventNhits); fClassifierResult.SetClassification("gamma", gamma); fClassifierResult.SetValid(true); return fClassifierResult; }