Changeset View
Changeset View
Standalone View
Standalone View
src/LauASqMaxFinder.cc
- This file was added.
/* | |||||
Copyright 2022 University of Warwick | |||||
Licensed under the Apache License, Version 2.0 (the "License"); | |||||
you may not use this file except in compliance with the License. | |||||
You may obtain a copy of the License at | |||||
http://www.apache.org/licenses/LICENSE-2.0 | |||||
Unless required by applicable law or agreed to in writing, software | |||||
distributed under the License is distributed on an "AS IS" BASIS, | |||||
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. | |||||
See the License for the specific language governing permissions and | |||||
limitations under the License. | |||||
*/ | |||||
/* | |||||
Laura++ package authors: | |||||
John Back | |||||
Paul Harrison | |||||
Thomas Latham | |||||
*/ | |||||
/*! \file LauASqMaxFinder.cc | |||||
\brief File containing implementation of LauASqMaxFinder class. | |||||
*/ | |||||
#include "LauASqMaxFinder.hh" | |||||
#include "LauAbsResonance.hh" | |||||
#include "LauFitter.hh" | |||||
#include "LauIsobarDynamics.hh" | |||||
#include "LauKinematics.hh" | |||||
#include "LauParameter.hh" | |||||
#include "TSystem.h" | |||||
#include <array> | |||||
#include <iostream> | |||||
#include <memory> | |||||
ClassImp( LauASqMaxFinder ) | |||||
LauASqMaxFinder::LauASqMaxFinder(LauIsobarDynamics& iso) : | |||||
iso_{iso} | |||||
{ | |||||
const LauKinematics* kinematics { iso_.getKinematics() }; | |||||
const Double_t m13SqMin{kinematics->getm13SqMin()}; | |||||
const Double_t m13SqMax{kinematics->getm13SqMax()}; | |||||
const Double_t m23SqMin{kinematics->getm23SqMin()}; | |||||
const Double_t m23SqMax{kinematics->getm23SqMax()}; | |||||
ownedParams_.resize(2); | |||||
ownedParams_[0] = std::make_unique<LauParameter>("m13Sq", 0.5*(m13SqMin+m13SqMax), m13SqMin, m13SqMax, kFALSE ); | |||||
ownedParams_[1] = std::make_unique<LauParameter>("m23Sq", 0.5*(m23SqMin+m23SqMax), m23SqMin, m23SqMax, kFALSE ); | |||||
params_.resize(2); | |||||
params_[0] = ownedParams_[0].get(); | |||||
params_[1] = ownedParams_[1].get(); | |||||
} | |||||
void LauASqMaxFinder::setParsFromMinuit(Double_t* par, Int_t npar) | |||||
{ | |||||
if ( static_cast<std::size_t>(npar) != params_.size() ) { | |||||
std::cerr << "ERROR in LauASqMaxFinder::setParsFromMinuit : Unexpected number of free parameters: " << npar << ".\n"; | |||||
std::cerr << " Expected: " << params_.size() << ".\n" << std::endl; | |||||
gSystem->Exit(EXIT_FAILURE); | |||||
} | |||||
for ( Int_t i{0}; i < npar; ++i ) { | |||||
params_[i]->value( par[i] ); | |||||
} | |||||
return; | |||||
} | |||||
Double_t LauASqMaxFinder::getTotNegLogLikelihood() | |||||
{ | |||||
const Double_t m13Sq { params_[0]->unblindValue() }; | |||||
const Double_t m23Sq { params_[1]->unblindValue() }; | |||||
if ( ! iso_.getKinematics()->withinDPLimits(m13Sq, m23Sq) ) { | |||||
if ( printMinimisationWarnings_ ) { | |||||
std::cerr << "WARNING in LauASqMaxFinder::getTotNegLogLikelihood : Fit has wandered outside the phase space boundary\n"; | |||||
std::cerr << " : Returning 0.0 value (further warnings of this type will be suppressed)" << std::endl; | |||||
printMinimisationWarnings_ = kFALSE; | |||||
} | |||||
return 0.0; | |||||
} | |||||
iso_.calcLikelihoodInfo( m13Sq, m23Sq ); | |||||
return -iso_.getEvtIntensity(); | |||||
} | |||||
Double_t LauASqMaxFinder::find() | |||||
{ | |||||
Double_t ASqMax{-1.0}; | |||||
const LauKinematics* kinematics { iso_.getKinematics() }; | |||||
LauKinematics localKin { | |||||
kinematics->getm1(), | |||||
kinematics->getm2(), | |||||
kinematics->getm3(), | |||||
kinematics->getmParent() | |||||
}; | |||||
const Double_t m13SqMin { localKin.getm13Min() }; | |||||
const Double_t m13SqMax { localKin.getm13Max() }; | |||||
const Double_t m23SqMin { localKin.getm23Min() }; | |||||
const Double_t m23SqMax { localKin.getm23Max() }; | |||||
LauFitter::setFitterMaxPars( 2 ); | |||||
LauFitter::setFitterVerbosity( LauOutputLevel::None ); | |||||
for ( UInt_t iRes{0}; iRes < iso_.getnTotAmp(); ++iRes ) { | |||||
// Find the center of the resonance and set that as the start of the search | |||||
const LauAbsResonance* res { const_cast<const LauIsobarDynamics&>(iso_).getResonance(iRes) }; | |||||
const Int_t bachelor { res->getPairInt() }; | |||||
const Double_t mass { res->getMass() }; | |||||
const Double_t width { res->getWidth() }; | |||||
const std::array<Double_t,3> massVals { mass - width, mass, mass + width }; | |||||
const std::array<Double_t,3> cosHelVals { -1.0, 0.0, 1.0 }; | |||||
for ( Double_t massVal : massVals ) { | |||||
for ( Double_t cosHel : cosHelVals ) { | |||||
switch ( bachelor ) { | |||||
case 0 : | |||||
localKin.updateKinematics( 0.5*(m13SqMax+m13SqMin), 0.5*(m23SqMax+m23SqMin) ); | |||||
break; | |||||
case 1 : | |||||
localKin.updateKinematicsFrom23( massVal, cosHel ); | |||||
break; | |||||
case 2 : | |||||
localKin.updateKinematicsFrom13( massVal, cosHel ); | |||||
break; | |||||
case 3 : | |||||
localKin.updateKinematicsFrom12( massVal, cosHel ); | |||||
break; | |||||
default : | |||||
break; | |||||
} | |||||
params_[0]->initValue( localKin.getm13Sq() ); | |||||
params_[1]->initValue( localKin.getm23Sq() ); | |||||
params_[0]->error(0.1); | |||||
params_[1]->error(0.1); | |||||
LauFitter::fitter().initialise( this, params_ ); | |||||
const LauAbsFitter::FitStatus fitResult { LauFitter::fitter().minimise() }; | |||||
if ( fitResult.status != 0 and -fitResult.NLL > ASqMax ) { | |||||
ASqMax = -fitResult.NLL; | |||||
} | |||||
} | |||||
} | |||||
} | |||||
LauFitter::destroyFitter(); | |||||
if ( ASqMax < 0.0 ) { | |||||
std::cerr << "WARNING in LauASqMaxFinder::find : Could not find ASqMax\n"; | |||||
std::cerr << " : Returning -1.0" << std::endl; | |||||
} | |||||
return ASqMax; | |||||
} |