diff --git a/src/LauSimFitCoordinator.cc b/src/LauSimFitCoordinator.cc index b943917..0bcb97b 100644 --- a/src/LauSimFitCoordinator.cc +++ b/src/LauSimFitCoordinator.cc @@ -1,1050 +1,1050 @@ /* Copyright 2013 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 LauSimFitCoordinator.cc \brief File containing implementation of LauSimFitCoordinator class. */ #include #include #include #include "TMath.h" #include "TMatrixD.h" #include "TMessage.h" #include "TMonitor.h" #include "TObjArray.h" #include "TObjString.h" #include "TServerSocket.h" #include "TSocket.h" #include "TSystem.h" #include "LauAbsFitter.hh" #include "LauFitNtuple.hh" #include "LauFitter.hh" #include "LauFormulaPar.hh" #include "LauParameter.hh" #include "LauParamFixed.hh" #include "LauSimFitCoordinator.hh" ClassImp(LauSimFitCoordinator) LauSimFitCoordinator::LauSimFitCoordinator( UInt_t numTasks, UInt_t port ) : nTasks_(numTasks), reqPort_(port), socketMonitor_(0), messageFromTask_(0), fitNtuple_(0) { messagesToTasks_.resize( nTasks_ ); for ( UInt_t iTask(0); iTask < nTasks_; ++iTask ) { messagesToTasks_[iTask] = new TMessage(); } } LauSimFitCoordinator::~LauSimFitCoordinator() { delete socketMonitor_; socketMonitor_ = 0; // Tell all tasks that they are finished and delete corresponding socket TString msgStr("Finish"); TMessage message( kMESS_STRING ); message.WriteTString(msgStr); for ( std::vector::iterator iter = socketTasks_.begin(); iter != socketTasks_.end(); ++iter ) { (*iter)->Send(message); (*iter)->Close(); delete (*iter); } socketTasks_.clear(); // Remove the components created to apply constraints to fit parameters for (std::vector::iterator iter = conVars_.begin(); iter != conVars_.end(); ++iter){ if ( !(*iter)->isLValue() ){ delete (*iter); (*iter) = 0; } } conVars_.clear(); // Remove all fit parameters for ( std::vector::iterator iter = params_.begin(); iter != params_.end(); ++iter ) { delete *iter; } params_.clear(); for ( std::vector::iterator iter = vectorPar_.begin(); iter != vectorPar_.end(); ++iter ) { delete[] (*iter); } vectorPar_.clear(); delete messageFromTask_; messageFromTask_ = 0; for ( std::vector::iterator iter = messagesToTasks_.begin(); iter != messagesToTasks_.end(); ++iter ) { delete (*iter); } messagesToTasks_.clear(); delete fitNtuple_; } void LauSimFitCoordinator::initSockets() { if ( socketMonitor_ != 0 ) { std::cerr << "ERROR in LauSimFitCoordinator::initSockets : Sockets already initialised." << std::endl; return; } //initialise socket connection, then accept a connection and return a full-duplex communication socket. socketMonitor_ = new TMonitor(); TServerSocket *ss = new TServerSocket( reqPort_, kFALSE ); UInt_t actual_port = ss->GetLocalPort(); - std::cout << "INFO in LauSimFitCoordinator::initSockets : Waiting for connection with " << nTasks_ << " workers on port " << actual_port << std::endl; + std::cout << "INFO in LauSimFitCoordinator::initSockets : Waiting for connection with " << nTasks_ << " tasks on port " << actual_port << std::endl; socketTasks_.resize(nTasks_); for ( UInt_t iTask(0); iTaskAccept(); std::cout << " : Added task " << iTask << std::endl; } // tell the clients to start std::cout << "INFO in LauSimFitCoordinator::initSockets : Initialising tasks" << std::endl; for ( UInt_t iTask(0); iTaskuseAsymmFitErrors()); socketTasks_[iTask]->Send(message); socketMonitor_->Add(socketTasks_[iTask]); } std::cout << " : Now start fit\n" << std::endl; ss->Close(); delete ss; } /* * OLD VERSION THAT JUST GETS THE NAMES - COULD HAVE A SERIES OF EXCHANGES TO GET THE NAMES, INIT VALUES, RANGES, ETC. INSTEAD OF PASSING PARAMETERS * THIS INCREASES THE GENERALITY OF THE CODE, I.E. THERE IS NO NEED FOR THE TASKS TO KNOW ANY LAURA++ CLASS BUT THIS ONE, BUT MAKES IT RATHER MORE DENSE * FOR THE MOMENT I WILL STICK WITH THE METHOD OF PASSING LAUPARAMETER OBJECTS AROUND AND CONSIDER GOING BACK TO THIS GENERAL METHOD ONCE EVERYTHING IS WORKING * void LauSimFitCoordinator::getParametersFromTasksFirstTime() { taskIndices_.resize( nTasks_ ); TSocket* sActive(0); for ( UInt_t iTask(0); iTaskSend(message); // Wait to receive the response and check that it has come from the task we just requested from sActive = socketMonitor_->Select(); if ( sActive != socketTasks_[iTask] ) { std::cerr << "ERROR in LauSimFitCoordinator::getParametersFromTasksFirstTime : Received message from a different task than expected!" << std::endl; gSystem->Exit(1); } // Read the object and extract the parameter names socketTasks_[iTask]->Recv( messageFromTask_ ); TObjArray * objarray = dynamic_cast( messageFromTask_->ReadObject( messageFromTask_->GetClass() ) ); if ( ! objarray ) { std::cerr << "ERROR in LauSimFitCoordinator::getParametersFromTasksFirstTime : Error reading parameter names from task" << std::endl; gSystem->Exit(1); } Int_t nPars = objarray->GetEntries(); for ( Int_t iPar(0); iPar < nPars; ++iPar ) { TObjString* objstring = dynamic_cast( (*objarray)[iPar] ); if ( ! objstring ) { std::cerr << "ERROR in LauSimFitCoordinator::getParametersFromTasksFirstTime : Error reading parameter names from task" << std::endl; gSystem->Exit(1); } TString parname = objstring->GetString(); std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname ); if ( iter != parIndices_.end() ) { UInt_t index = iter->second; taskIndices_[iTask].push_back( index ); } else { UInt_t index = parIndices_.size(); parIndices_.insert( std::make_pair( parname, index ) ); parNames_.insert( std::make_pair( index, parname ) ); taskIndices_[iTask].push_back( index ); } } delete objarray; objarray = 0; delete messageFromTask_; messageFromTask_ = 0; } UInt_t nPars = parNames_.size(); parValues_.resize( nPars ); } */ void LauSimFitCoordinator::getParametersFromTasks() { if ( socketMonitor_ == 0 ) { std::cerr << "ERROR in LauSimFitCoordinator::getParametersFromTasks : Sockets not initialised." << std::endl; return; } if ( params_.empty() ) { this->getParametersFromTasksFirstTime(); // Add variables to Gaussian constrain to a list this->addConParameters(); } else { this->updateParametersFromTasks(); } } void LauSimFitCoordinator::updateParametersFromTasks() { TSocket* sActive(0); // Construct a message, requesting the list of parameter names TString msgStr = "Send Parameters"; TMessage message( kMESS_STRING ); message.WriteTString( msgStr ); for ( UInt_t iTask(0); iTaskSend(message); // Wait to receive the response and check that it has come from the task we just requested from sActive = socketMonitor_->Select(); if ( sActive != socketTasks_[iTask] ) { std::cerr << "ERROR in LauSimFitCoordinator::updateParametersFromTasks : Received message from a different task than expected!" << std::endl; gSystem->Exit(1); } // Read the object and extract the parameter names socketTasks_[iTask]->Recv( messageFromTask_ ); TObjArray * objarray = dynamic_cast( messageFromTask_->ReadObject( messageFromTask_->GetClass() ) ); if ( ! objarray ) { std::cerr << "ERROR in LauSimFitCoordinator::updateParametersFromTasks : Error reading parameter names from task" << std::endl; gSystem->Exit(1); } // We want to auto-delete the supplied parameters since we only copy their values in this case objarray->SetOwner(kTRUE); const UInt_t nPars = objarray->GetEntries(); if ( nPars != taskIndices_[iTask].size() ) { std::cerr << "ERROR in LauSimFitCoordinator::updateParametersFromTasks : Unexpected number of parameters received from task" << std::endl; gSystem->Exit(1); } for ( UInt_t iPar(0); iPar < nPars; ++iPar ) { LauParameter* parameter = dynamic_cast( (*objarray)[iPar] ); if ( ! parameter ) { std::cerr << "ERROR in LauSimFitCoordinator::updateParametersFromTasks : Error reading parameter from task" << std::endl; gSystem->Exit(1); } TString parname = parameter->name(); Double_t parvalue = parameter->initValue(); std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname ); if ( iter == parIndices_.end() ) { std::cerr << "ERROR in LauSimFitCoordinator::updateParametersFromTasks : Unexpected parameter name received from task" << std::endl; gSystem->Exit(1); } const UInt_t index = iter->second; if ( taskIndices_[iTask][iPar] != index ) { std::cerr << "ERROR in LauSimFitCoordinator::updateParametersFromTasks : Unexpected parameter received from task" << std::endl; gSystem->Exit(1); } params_[index]->initValue( parvalue ); parValues_[index] = parvalue; vectorPar_[iTask][iPar] = parvalue; this->checkParameter( parameter, index ); } delete objarray; objarray = 0; delete messageFromTask_; messageFromTask_ = 0; } } void LauSimFitCoordinator::getParametersFromTasksFirstTime() { taskIndices_.resize( nTasks_ ); taskFreeIndices_.resize( nTasks_ ); vectorPar_.resize( nTasks_ ); vectorRes_.resize( nTasks_ ); TSocket* sActive(0); // Construct a message, requesting the list of parameter names TString msgStr = "Send Parameters"; TMessage message( kMESS_STRING ); message.WriteTString( msgStr ); for ( UInt_t iTask(0); iTaskSend(message); // Wait to receive the response and check that it has come from the task we just requested from sActive = socketMonitor_->Select(); if ( sActive != socketTasks_[iTask] ) { std::cerr << "ERROR in LauSimFitCoordinator::getParametersFromTasksFirstTime : Received message from a different task than expected!" << std::endl; gSystem->Exit(1); } // Read the object and extract the parameter names socketTasks_[iTask]->Recv( messageFromTask_ ); TObjArray * objarray = dynamic_cast( messageFromTask_->ReadObject( messageFromTask_->GetClass() ) ); if ( ! objarray ) { std::cerr << "ERROR in LauSimFitCoordinator::getParametersFromTasksFirstTime : Error reading parameters from task" << std::endl; gSystem->Exit(1); } const UInt_t nPars = objarray->GetEntries(); vectorPar_[iTask] = new Double_t[nPars]; for ( UInt_t iPar(0); iPar < nPars; ++iPar ) { LauParameter* parameter = dynamic_cast( (*objarray)[iPar] ); if ( ! parameter ) { std::cerr << "ERROR in LauSimFitCoordinator::getParametersFromTasksFirstTime : Error reading parameter from task" << std::endl; gSystem->Exit(1); } TString parname = parameter->name(); Double_t parvalue = parameter->initValue(); Bool_t parfixed = parameter->fixed(); std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname ); if ( iter != parIndices_.end() ) { UInt_t index = iter->second; taskIndices_[iTask].push_back( index ); if ( ! parfixed ) { taskFreeIndices_[iTask].push_back( index ); } this->checkParameter( parameter, index ); } else { UInt_t index = parIndices_.size(); parIndices_.insert( std::make_pair( parname, index ) ); parNames_.insert( std::make_pair( index, parname ) ); taskIndices_[iTask].push_back( index ); if ( ! parfixed ) { taskFreeIndices_[iTask].push_back( index ); } params_.push_back( parameter ); parValues_.push_back( parvalue ); } vectorPar_[iTask][iPar] = parvalue; } delete objarray; objarray = 0; delete messageFromTask_; messageFromTask_ = 0; } } void LauSimFitCoordinator::printParInfo() const { for ( UInt_t iTask(0); iTask& indices = taskIndices_[iTask]; std::cout << "INFO in LauSimFitCoordinator::printParInfo : Task " << iTask << " has the following parameters:\n"; for ( std::vector::const_iterator iter = indices.begin(); iter != indices.end(); ++iter ) { const TString& parName = parNames_.find(*iter)->second; Double_t parValue = parValues_[*iter]; const LauParameter* par = params_[*iter]; if ( par->name() != parName || par->initValue() != parValue ) { std::cerr << "ERROR in LauSimFitCoordinator::printParInfo : Discrepancy in parameter name and value records, this is very strange!!" << std::endl; } std::cout << " : " << parName << " = " << parValue << " and has index " << *iter << "\n"; } std::cout << std::endl; } std::cout << "INFO in LauSimFitCoordinator::printParInfo : " << "There are " << params_.size() << " parameters in total" << std::endl; } void LauSimFitCoordinator::checkParameter( const LauParameter* param, UInt_t index ) const { const LauParameter* storedPar = params_[index]; TString parName = storedPar->name(); if ( param->name() != parName ) { std::cerr << "ERROR in LauSimFitCoordinator::checkParameter : Parameter name is different!! This shouldn't happen!!" << std::endl; } if ( param->initValue() != storedPar->initValue() ) { std::cerr << "WARNING in LauSimFitCoordinator::checkParameter : Initial value for parameter " << parName << " is different, will use the value first set: " << storedPar->initValue() << std::endl; } if ( param->minValue() != storedPar->minValue() ) { std::cerr << "WARNING in LauSimFitCoordinator::checkParameter : Minimum allowed value for parameter " << parName << " is different, will use the value first set: " << storedPar->minValue() << std::endl; } if ( param->maxValue() != storedPar->maxValue() ) { std::cerr << "WARNING in LauSimFitCoordinator::checkParameter : Maximum allowed value for parameter " << parName << " is different, will use the value first set: " << storedPar->maxValue() << std::endl; } if ( param->fixed() != storedPar->fixed() ) { std::cerr << "WARNING in LauSimFitCoordinator::checkParameter : Fixed/floating property of parameter " << parName << " is different, will use the value first set: " << (storedPar->fixed() ? "fixed" : "floating") << std::endl; } if ( param->secondStage() != storedPar->secondStage() ) { std::cerr << "WARNING in LauSimFitCoordinator::checkParameter : Second stage property of parameter " << parName << " is different, will use the value first set: " << (storedPar->secondStage() ? "true" : "false") << std::endl; } } void LauSimFitCoordinator::initialise() { this->initSockets(); } void LauSimFitCoordinator::runSimFit( const TString& fitNtupleFileName, const UInt_t nExp, const UInt_t firstExp, const Bool_t useAsymmErrors, const Bool_t doTwoStageFit ) { // Routine to perform the total fit. // First, initialise this->useAsymmFitErrors(useAsymmErrors); this->twoStageFit(doTwoStageFit); this->initialise(); std::cout << "INFO in LauSimFitCoordinator::runSimFit : First experiment = " << firstExp << std::endl; std::cout << "INFO in LauSimFitCoordinator::runSimFit : Number of experiments = " << nExp << std::endl; // Start the cumulative timer cumulTimer_.Start(); this->resetFitCounters(); // Create and setup the fit results ntuple std::cout << "INFO in LauSimFitCoordinator::runSimFit : Creating fit ntuple." << std::endl; if (fitNtuple_ != 0) {delete fitNtuple_; fitNtuple_ = 0;} fitNtuple_ = new LauFitNtuple(fitNtupleFileName, useAsymmErrors); // Loop over the number of experiments for (UInt_t iExp = firstExp; iExp < (firstExp+nExp); ++iExp) { // Start the timer to see how long each fit takes timer_.Start(); this->setCurrentExperiment( iExp ); // Instruct the tasks to read the data for this experiment Bool_t readOK = this->readData(); if ( ! readOK ) { std::cerr << "ERROR in LauSimFitCoordinator::runSimFit : One or more tasks reported problems with reading data for experiment " << iExp << ", skipping..." << std::endl; timer_.Stop(); continue; } // Instruct the tasks to perform the caching this->cacheInputData(); // Do the fit this->fitExpt(); // Stop the timer and see how long the program took so far timer_.Stop(); timer_.Print(); // Instruct the tasks to finalise the results this->finalise(); } // Print out total timing info. std::cout << "INFO in LauSimFitCoordinator::runSimFit : Cumulative timing:" << std::endl; cumulTimer_.Stop(); cumulTimer_.Print(); // Print out stats on OK fits. const UInt_t nOKFits = this->numberOKFits(); const UInt_t nBadFits = this->numberBadFits(); std::cout << "INFO in LauSimFitCoordinator::runSimFit : Number of OK Fits = " << nOKFits << std::endl; std::cout << "INFO in LauSimFitCoordinator::runSimFit : Number of Failed Fits = " << nBadFits << std::endl; Double_t fitEff(0.0); if (nExp != 0) {fitEff = nOKFits/(1.0*nExp);} std::cout << "INFO in LauSimFitCoordinator::runSimFit : Fit efficiency = " << fitEff*100.0 << "%." << std::endl; // Instruct the tasks to write out any fit results (ntuples etc...). this->writeOutResults(); } void LauSimFitCoordinator::withinAsymErrorCalc(const Bool_t inAsymErrCalc) { this->LauFitObject::withinAsymErrorCalc(inAsymErrCalc); if ( socketMonitor_ == 0 ) { std::cerr << "ERROR in LauSimFitCoordinator::withinAsymErrorCalc : Sockets not initialised." << std::endl; return; } // Construct a message, informing the tasks whether or not we are now within the asymmetric error calculation TString msgStr("Asym Error Calc"); const Bool_t asymErrorCalc( this->withinAsymErrorCalc() ); TMessage message( kMESS_STRING ); message.WriteTString( msgStr ); message.WriteBool( asymErrorCalc ); // Send the message to the tasks for ( UInt_t iTask(0); iTaskSend(message); } TSocket* sActive(0); UInt_t responsesReceived(0); while ( responsesReceived != nTasks_ ) { // Get the next queued response sActive = socketMonitor_->Select(); // Extract from the message the ID of the task and the number of events read Bool_t response(kTRUE); UInt_t iTask(0); sActive->Recv( messageFromTask_ ); messageFromTask_->ReadUInt( iTask ); messageFromTask_->ReadBool( response ); if ( response != asymErrorCalc ) { std::cerr << "WARNING in LauSimFitCoordinator::withinAsymErrorCalc : Problem informing task " << iTask << std::endl; } ++responsesReceived; } } Bool_t LauSimFitCoordinator::readData() { if ( socketMonitor_ == 0 ) { std::cerr << "ERROR in LauSimFitCoordinator::readData : Sockets not initialised." << std::endl; return kFALSE; } // Construct a message, requesting to read the data for the given experiment TString msgStr("Read Expt"); const UInt_t iExp( this->iExpt() ); TMessage message( kMESS_STRING ); message.WriteTString( msgStr ); message.WriteUInt( iExp ); // Send the message to the tasks for ( UInt_t iTask(0); iTaskSend(message); } TSocket* sActive(0); UInt_t responsesReceived(0); Bool_t ok(kTRUE); while ( responsesReceived != nTasks_ ) { // Get the next queued response sActive = socketMonitor_->Select(); // Extract from the message the ID of the task and the number of events read sActive->Recv( messageFromTask_ ); UInt_t iTask(0); UInt_t nEvents(0); messageFromTask_->ReadUInt( iTask ); messageFromTask_->ReadUInt( nEvents ); if ( nEvents <= 0 ) { std::cerr << "ERROR in LauSimFitCoordinator::readData : Task " << iTask << " reports no events found for experiment " << iExp << std::endl; ok = kFALSE; } else { std::cerr << "INFO in LauSimFitCoordinator::readData : Task " << iTask << " reports " << nEvents << " events found for experiment " << iExp << std::endl; } ++responsesReceived; } return ok; } Bool_t LauSimFitCoordinator::cacheInputData() { if ( socketMonitor_ == 0 ) { std::cerr << "ERROR in LauSimFitCoordinator::cacheInputData : Sockets not initialised." << std::endl; return kFALSE; } // Construct a message, requesting it to read the data for the given experiment TString msgStr("Cache"); TMessage message( kMESS_STRING ); message.WriteTString( msgStr ); for ( UInt_t iTask(0); iTaskSend(message); } TSocket* sActive(0); UInt_t responsesReceived(0); Bool_t allOK(kTRUE); while ( responsesReceived != nTasks_ ) { // Get the next queued response sActive = socketMonitor_->Select(); // Extract from the message the ID of the task and the success/failure flag sActive->Recv( messageFromTask_ ); UInt_t iTask(0); Bool_t ok(kTRUE); messageFromTask_->ReadUInt( iTask ); messageFromTask_->ReadBool( ok ); if ( ! ok ) { std::cerr << "ERROR in LauSimFitCoordinator::cacheInputData : Task " << iTask << " reports an error performing caching" << std::endl; allOK = kFALSE; } ++responsesReceived; } return allOK; } void LauSimFitCoordinator::checkInitFitParams() { this->getParametersFromTasks(); this->printParInfo(); } void LauSimFitCoordinator::fitExpt() { // Routine to perform the actual fit for the given experiment // Instruct the tasks to update initial fit parameters if required (e.g. if using random numbers). this->checkInitFitParams(); // Initialise the fitter LauFitter::fitter()->useAsymmFitErrors( this->useAsymmFitErrors() ); LauFitter::fitter()->twoStageFit( this->twoStageFit() ); LauFitter::fitter()->initialise( this, params_ ); this->startNewFit( LauFitter::fitter()->nParameters(), LauFitter::fitter()->nFreeParameters() ); // Now ready for minimisation step std::cout << "\nINFO in LauSimFitCoordinator::fitExpt : Start minimisation...\n"; LauAbsFitter::FitStatus fitResult = LauFitter::fitter()->minimise(); // If we're doing a two stage fit we can now release (i.e. float) // the 2nd stage parameters and re-fit if (this->twoStageFit()) { if ( fitResult.status != 3 ) { std::cerr << "ERROR in LauSimFitCoordinator:fitExpt : Not running second stage fit since first stage failed." << std::endl; LauFitter::fitter()->releaseSecondStageParameters(); } else { LauFitter::fitter()->releaseSecondStageParameters(); this->startNewFit( LauFitter::fitter()->nParameters(), LauFitter::fitter()->nFreeParameters() ); fitResult = LauFitter::fitter()->minimise(); } } const TMatrixD& covMat = LauFitter::fitter()->covarianceMatrix(); this->storeFitStatus( fitResult, covMat ); // Store the final fit results and errors into protected internal vectors that // all sub-classes can use within their own finalFitResults implementation // used below (e.g. putting them into an ntuple in a root file) LauFitter::fitter()->updateParameters(); } void LauSimFitCoordinator::setParsFromMinuit(Double_t* par, Int_t npar) { // This function sets the internal parameters based on the values // that Minuit is using when trying to minimise the total likelihood function. // MINOS reports different numbers of free parameters depending on the // situation, so disable this check if ( ! this->withinAsymErrorCalc() ) { const UInt_t nFreePars = this->nFreeParams(); if (static_cast(npar) != nFreePars) { std::cerr << "ERROR in LauSimFitCoordinator::setParsFromMinuit : Unexpected number of free parameters: " << npar << ".\n"; std::cerr << " Expected: " << nFreePars << ".\n" << std::endl; gSystem->Exit(EXIT_FAILURE); } } // Despite npar being the number of free parameters // the par array actually contains all the parameters, // free and floating... // Update all the parameters with their new values. // Change the value in the array to be sent out to the tasks and the // parameters themselves (so that constraints are correctly calculated) for (UInt_t i(0); inTotParams(); ++i) { if (!params_[i]->fixed()) { parValues_[i] = par[i]; params_[i]->value(par[i]); } } } Double_t LauSimFitCoordinator::getTotNegLogLikelihood() { if ( socketMonitor_ == 0 ) { std::cerr << "ERROR in LauSimFitCoordinator::getTotNegLogLikelihood : Sockets not initialised." << std::endl; return 0.0; } // Send current values of the parameters to the tasks. for ( UInt_t iTask(0); iTask& indices = taskIndices_[iTask]; std::vector& freeIndices = taskFreeIndices_[iTask]; UInt_t nPars = indices.size(); UInt_t nFreePars = freeIndices.size(); for ( UInt_t iPar(0); iPar < nPars; ++iPar ) { vectorPar_[iTask][iPar] = parValues_[ indices[iPar] ]; } TMessage* message = messagesToTasks_[iTask]; message->Reset( kMESS_ANY ); message->WriteUInt( nPars ); message->WriteUInt( nFreePars ); message->WriteFastArray( vectorPar_[iTask], nPars ); socketTasks_[iTask]->Send(*message); } Double_t negLogLike(0.0); TSocket *sActive(0); UInt_t responsesReceived(0); Bool_t allOK(kTRUE); while ( responsesReceived != nTasks_ ) { sActive = socketMonitor_->Select(); sActive->Recv(messageFromTask_); messageFromTask_->ReadDouble( vectorRes_[responsesReceived] ); Double_t& nLL = vectorRes_[responsesReceived]; if ( nLL == 0.0 || TMath::IsNaN(nLL) || !TMath::Finite(nLL) ) { allOK = kFALSE; } negLogLike += vectorRes_[responsesReceived]; ++responsesReceived; } // Calculate any penalty terms from Gaussian constrained variables if ( ! conVars_.empty() ){ negLogLike += this->getLogLikelihoodPenalty(); } const Double_t worstNegLogLike = -1.0*this->worstLogLike(); if ( ! allOK ) { std::cerr << "WARNING in LauSimFitCoordinator::getTotNegLogLikelihood : Strange NLL value returned by one or more tasks\n"; std::cerr << " : Returning worst NLL found so far to force MINUIT out of this region." << std::endl; negLogLike = worstNegLogLike; } else if ( negLogLike > worstNegLogLike ) { this->worstLogLike( -negLogLike ); } return negLogLike; } Double_t LauSimFitCoordinator::getLogLikelihoodPenalty() { Double_t penalty(0.0); for ( std::vector::const_iterator iter = conVars_.begin(); iter != conVars_.end(); ++iter ) { Double_t val = (*iter)->unblindValue(); Double_t mean = (*iter)->constraintMean(); Double_t width = (*iter)->constraintWidth(); Double_t term = ( val - mean )*( val - mean ); penalty += term/( 2*width*width ); } return penalty; } void LauSimFitCoordinator::addConParameters() { // Add penalties from the constraints to fit parameters // First, constraints on the fit parameters themselves for ( std::vector::const_iterator iter = params_.begin(); iter != params_.end(); ++iter ) { if ( (*iter)->gaussConstraint() ) { conVars_.push_back( *iter ); std::cout << "INFO in LauSimFitCoordinator::addConParameters : Added Gaussian constraint to parameter "<< (*iter)->name() << std::endl; } } // Second, constraints on arbitrary combinations const std::vector& storeCon = this->constraintsStore(); for ( std::vector::const_iterator iter = storeCon.begin(); iter != storeCon.end(); ++iter ) { const std::vector& names = (*iter).conPars_; std::vector params; for ( std::vector::const_iterator iternames = names.begin(); iternames != names.end(); ++iternames ) { for ( std::vector::const_iterator iterfit = params_.begin(); iterfit != params_.end(); ++iterfit ) { if ( (*iternames) == (*iterfit)->name() ){ params.push_back(*iterfit); } } } // If the parameters are not found, skip it if ( params.size() != (*iter).conPars_.size() ) { std::cerr << "WARNING in LauSimFitCoordinator::addConParameters: Could not find parameters to constrain in the formula... skipping" << std::endl; continue; } LauFormulaPar* formPar = new LauFormulaPar( (*iter).formula_, (*iter).formula_, params ); formPar->addGaussianConstraint( (*iter).mean_, (*iter).width_ ); conVars_.push_back(formPar); std::cout << "INFO in LauSimFitCoordinator::addConParameters : Added Gaussian constraint to formula\n"; std::cout << " : Formula: " << (*iter).formula_ << std::endl; for ( std::vector::iterator iterparam = params.begin(); iterparam != params.end(); ++iterparam ) { std::cout << " : Parameter: " << (*iterparam)->name() << std::endl; } } } Bool_t LauSimFitCoordinator::finalise() { if ( socketMonitor_ == 0 ) { std::cerr << "ERROR in LauSimFitCoordinator::finalise : Sockets not initialised." << std::endl; return kFALSE; } // Prepare the covariance matrices const TMatrixD& covMatrix = this->covarianceMatrix(); covMatrices_.resize( nTasks_ ); LauParamFixed pred; std::map freeParIndices; UInt_t counter(0); for ( UInt_t iPar(0); iPar < this->nTotParams(); ++iPar ) { const LauParameter* par = params_[iPar]; if ( ! pred(par) ) { freeParIndices.insert( std::make_pair(iPar,counter) ); ++counter; } } for ( UInt_t iTask(0); iTask freeIndices; freeIndices.reserve( nPar ); for ( UInt_t iPar(0); iPar < nPar; ++iPar ) { UInt_t index = taskIndices_[iTask][iPar]; std::map::iterator freeIter = freeParIndices.find(index); if ( freeIter == freeParIndices.end() ) { continue; } UInt_t freeIndex = freeIter->second; freeIndices.push_back( freeIndex ); } const UInt_t nFreePars = freeIndices.size(); TMatrixD& covMat = covMatrices_[iTask]; covMat.ResizeTo( nFreePars, nFreePars ); for ( UInt_t iPar(0); iPar < nFreePars; ++iPar ) { for ( UInt_t jPar(0); jPar < nFreePars; ++jPar ) { UInt_t i = freeIndices[iPar]; UInt_t j = freeIndices[jPar]; covMat( iPar, jPar ) = covMatrix( i, j ); } } } // The array to hold the parameters TObjArray array; // Send messages to all tasks containing the final parameters and fit status, NLL for ( UInt_t iTask(0); iTask& indices = taskIndices_[iTask]; UInt_t nPars = indices.size(); for ( UInt_t iPar(0); iPar < nPars; ++iPar ) { array.Add( params_[ indices[iPar] ] ); } const Int_t status = this->statusCode(); const Double_t NLL = this->nll(); const Double_t EDM = this->edm(); TMatrixD& covMat = covMatrices_[iTask]; TMessage* message = messagesToTasks_[iTask]; message->Reset( kMESS_OBJECT ); message->WriteInt( status ); message->WriteDouble( NLL ); message->WriteDouble( EDM ); message->WriteObject( &array ); message->WriteObject( &covMat ); socketTasks_[iTask]->Send(*message); } TSocket *sActive(0); UInt_t responsesReceived(0); Bool_t allOK(kTRUE); while ( responsesReceived != nTasks_ ) { // Get the next queued response sActive = socketMonitor_->Select(); // Extract from the message the ID of the task and the number of events read sActive->Recv( messageFromTask_ ); UInt_t iTask(0); Bool_t ok(kTRUE); messageFromTask_->ReadUInt( iTask ); messageFromTask_->ReadBool( ok ); if ( ok ) { TObjArray * objarray = dynamic_cast( messageFromTask_->ReadObject( messageFromTask_->GetClass() ) ); if ( ! objarray ) { std::cerr << "ERROR in LauSimFitCoordinator::finalise : Error reading finalised parameters from task" << std::endl; allOK = kFALSE; } else { // We want to auto-delete the supplied parameters since we only copy their values in this case objarray->SetOwner(kTRUE); const UInt_t nPars = objarray->GetEntries(); if ( nPars != taskIndices_[iTask].size() ) { std::cerr << "ERROR in LauSimFitCoordinator::finalise : Unexpected number of finalised parameters received from task" << std::endl; allOK = kFALSE; } else { for ( UInt_t iPar(0); iPar < nPars; ++iPar ) { LauParameter* parameter = dynamic_cast( (*objarray)[iPar] ); if ( ! parameter ) { std::cerr << "ERROR in LauSimFitCoordinator::finalise : Error reading parameter from task" << std::endl; allOK = kFALSE; continue; } TString parname = parameter->name(); std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname ); if ( iter == parIndices_.end() ) { std::cerr << "ERROR in LauSimFitCoordinator::finalise : Unexpected parameter name received from task" << std::endl; allOK = kFALSE; continue; } const UInt_t index = iter->second; if ( taskIndices_[iTask][iPar] != index ) { std::cerr << "ERROR in LauSimFitCoordinator::finalise : Unexpected parameter received from task" << std::endl; allOK = kFALSE; continue; } Double_t parvalue = parameter->value(); params_[index]->value( parvalue ); parValues_[index] = parvalue; vectorPar_[iTask][iPar] = parvalue; } } delete objarray; } } else { std::cerr << "ERROR in LauSimFitCoordinator::finalise : Task " << iTask << " reports an error performing finalisation" << std::endl; allOK = kFALSE; } ++responsesReceived; } // Fill our ntuple as well if ( fitNtuple_ != 0 ) { for ( std::vector::iterator iter = params_.begin(); iter != params_.end(); ++iter ) { if (!(*iter)->fixed()) { (*iter)->updatePull(); } } std::vector extraVars; fitNtuple_->storeParsAndErrors( params_, extraVars ); fitNtuple_->storeCorrMatrix(this->iExpt(), this->fitStatus(), this->covarianceMatrix()); fitNtuple_->updateFitNtuple(); } return allOK; } Bool_t LauSimFitCoordinator::writeOutResults() { if ( socketMonitor_ == 0 ) { std::cerr << "ERROR in LauSimFitCoordinator::writeOutResults : Sockets not initialised." << std::endl; return kFALSE; } // Construct a message, requesting to write out the fit results TString msgStr("Write Results"); TMessage message( kMESS_STRING ); message.WriteTString( msgStr ); // Send the message to the tasks for ( UInt_t iTask(0); iTaskSend(message); } TSocket *sActive(0); UInt_t responsesReceived(0); Bool_t allOK(kTRUE); while ( responsesReceived != nTasks_ ) { // Get the next queued response sActive = socketMonitor_->Select(); // Extract from the message the ID of the task and the number of events read sActive->Recv( messageFromTask_ ); UInt_t iTask(0); Bool_t ok(kTRUE); messageFromTask_->ReadUInt( iTask ); messageFromTask_->ReadBool( ok ); if ( ! ok ) { std::cerr << "ERROR in LauSimFitCoordinator::writeOutResults : Task " << iTask << " reports an error performing finalisation" << std::endl; allOK = kFALSE; } ++responsesReceived; } // Write out our ntuple as well if (fitNtuple_ != 0) { fitNtuple_->writeOutFitResults(); } return allOK; }