diff --git a/AmpTools/IUAmpTools/AmpToolsInterface.cc b/AmpTools/IUAmpTools/AmpToolsInterface.cc index 067f439d..f8bab08e 100644 --- a/AmpTools/IUAmpTools/AmpToolsInterface.cc +++ b/AmpTools/IUAmpTools/AmpToolsInterface.cc @@ -317,6 +317,46 @@ AmpToolsInterface::reinitializePars(){ minuitMinimizationManager()->resetErrors(); } +void +AmpToolsInterface::randomizeProductionPars( vector ampNames, float maxFitFraction ){ + + // shouldn't be callin' unless you're fittin' + if( m_functionality != kFull ) return; + + vector< AmplitudeInfo* > amps = m_configurationInfo->amplitudeList(); + + for( vector< AmplitudeInfo* >::iterator ampItr = amps.begin();ampItr != amps.end();++ampItr ){ + if( (**ampItr).fixed() ) continue; + bool skip=true; + for(unsigned int i=0;inumSignalEvents(); + // for the NI interface to use the cache to avoid a mess with MPI jobs and + // retrigging of NI calculation -- we don't require a precise number + // for this anyway! + double normInt = normIntInterface(reac)->normInt(ampName,ampName,true).real(); + double ampScale = intensityManager(reac)->getScale(ampName); + // fit fraction = ( scale^2 * prodPar^2 * normInt ) / numSignalEvents + double fitFraction = random( maxFitFraction ); + double prodParMag = sqrt( ( numSignalEvents * fitFraction ) /( ampScale * ampScale * normInt ) ); + double prodParPhase = ( (**ampItr).real() ? 0 : random( 2*PI ) ); + complex< double > prodPar( prodParMag*cos( prodParPhase ),prodParMag*sin( prodParPhase ) ); + m_parameterManager->setProductionParameter( ampName, prodPar ); + } + + // reset parameter steps after randomizing parameters + minuitMinimizationManager()->resetErrors(); +} + void AmpToolsInterface::randomizeProductionPars( float maxFitFraction ){ diff --git a/AmpTools/IUAmpTools/AmpToolsInterface.h b/AmpTools/IUAmpTools/AmpToolsInterface.h index 0a470aa3..a90a7b5c 100644 --- a/AmpTools/IUAmpTools/AmpToolsInterface.h +++ b/AmpTools/IUAmpTools/AmpToolsInterface.h @@ -241,6 +241,16 @@ class AmpToolsInterface{ void reinitializePars( ); + /** This function will randomly set the production parameters containing a string in ampNames in + * the likelihood calculator. It is useful when searching for multiple + * ambiguous solutions. Production parameters will be set such that + * starting fit fraction varies between zero and fraction specified + * by the optional argument. The phase will be set randomly between + * zero and 2 pi. + */ + + void randomizeProductionPars( vector ampNames, float maxFitFraction = 1 ); + /** This function will randomly set the production parameters in * the likelihood calculator. It is useful when searching for multiple * ambiguous solutions. Production parameters will be set such that diff --git a/AmpTools/IUAmpTools/FitResults.cc b/AmpTools/IUAmpTools/FitResults.cc index f188ae78..697d6fbb 100644 --- a/AmpTools/IUAmpTools/FitResults.cc +++ b/AmpTools/IUAmpTools/FitResults.cc @@ -502,6 +502,12 @@ FitResults::productionParameter( const string& ampName ) const { return complex< double >( m_parValues[iRe], m_parValues[iIm] ); } +void FitResults::setProductionParameter( const string& ampName, complex value ) { + int iRe = m_parIndex.find( realProdParName( ampName ) )->second; + int iIm = m_parIndex.find( imagProdParName( ampName ) )->second; + m_parValues.at(iRe) = value.real(); + m_parValues.at(iIm) = value.imag(); +} complex< double > FitResults::scaledProductionParameter( const string& ampName ) const { @@ -613,6 +619,43 @@ FitResults::ampList( const string& reaction ) const { return m_ampNames[reacIndex->second]; } +void +FitResults::writeFitResults( const string& outFile ){ + + ofstream output( outFile.c_str() ); + + output << m_numAmps[0]/4. << "\t0" << endl; + + for( unsigned int i = 0; i < m_parNames.size()-1; i++ ){ + string reactName = m_parNames[i].substr( 8, 4 ); + if( strcmp( reactName.c_str(), "amp1" ) != 0 ) continue; + string lastChars = m_parNames[i].substr( m_parNames[i].size() - 2, 2 ); + if( strcmp( lastChars.c_str(), "im" ) == 0 ){ + output << m_parNames[i].substr( 0, m_parNames[i].size() - 3 ); + if( m_parValues[i] == 0 && m_covMatrix[i][i] == 0 ) output << "+"; + output << "\t(" << m_parValues[i-1] << "," << m_parValues[i] << ")" << endl; + } + } + + for( unsigned int i = 0; i < m_parNames.size()-1; i++ ){ + string reactNamei = m_parNames[i].substr( 8, 4 ); + if( strcmp( reactNamei.c_str(), "amp1" ) != 0 ) continue; + string lastCharsi = m_parNames[i].substr( m_parNames[i].size() - 2, 2 ); + if( strcmp( lastCharsi.c_str(), "im" ) == 0 && m_parValues[i] == 0 && m_covMatrix[i][i] == 0 ) continue; + for( unsigned int j = 0; j < m_parNames.size()-1; j++ ){ + string reactNamej = m_parNames[j].substr( 8, 4 ); + if( strcmp( reactNamej.c_str(), "amp1" ) == 0 ){ + string lastCharsj = m_parNames[j].substr( m_parNames[j].size() - 2, 2 ); + if( strcmp( lastCharsj.c_str(), "im" ) == 0 && m_parValues[j] == 0 && m_covMatrix[j][j] == 0 ) continue; + output << m_covMatrix[i][j] << "\t"; + } + } + output << endl; + } + + output.close(); +} + void FitResults::writeResults( const string& outFile ) const { diff --git a/AmpTools/IUAmpTools/FitResults.h b/AmpTools/IUAmpTools/FitResults.h index d13fc4e0..f8af2f15 100644 --- a/AmpTools/IUAmpTools/FitResults.h +++ b/AmpTools/IUAmpTools/FitResults.h @@ -455,6 +455,9 @@ class FitResults */ void rotateResults(); + void writeFitResults( const string& fileName ); + void setProductionParameter( const string& ampName, complex value ); + protected: /**