laura is hosted by Hepforge, IPPP Durham
Laura++  v2r0
A maximum likelihood fitting package for performing Dalitz-plot analysis.
LauSimFitMaster.cc
Go to the documentation of this file.
1 
2 // Copyright University of Warwick 2013 - 2013.
3 // Distributed under the Boost Software License, Version 1.0.
4 // (See accompanying file LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5 
6 // Authors:
7 // Thomas Latham
8 // John Back
9 // Paul Harrison
10 
15 #include <cstdlib>
16 #include <iostream>
17 
18 #include "TMessage.h"
19 #include "TMonitor.h"
20 #include "TObjArray.h"
21 #include "TObjString.h"
22 #include "TServerSocket.h"
23 #include "TSocket.h"
24 #include "TSystem.h"
25 #include "TVectorD.h"
26 
27 #include "LauAbsFitter.hh"
28 #include "LauFitNtuple.hh"
29 #include "LauFitter.hh"
30 #include "LauParameter.hh"
31 #include "LauParamFixed.hh"
32 #include "LauSimFitMaster.hh"
33 
34 
35 ClassImp(LauSimFitMaster)
36 
37 
38 LauSimFitMaster::LauSimFitMaster( UInt_t numSlaves, UInt_t port ) :
39  nSlaves_(numSlaves),
40  reqPort_(port),
41  nParams_(0),
42  nFreeParams_(0),
43  withinAsymErrorCalc_(kFALSE),
44  numberOKFits_(0),
45  numberBadFits_(0),
46  fitStatus_(0),
47  iExpt_(0),
48  NLL_(0.0),
49  socketMonitor_(0),
50  messageFromSlave_(0),
51  fitNtuple_(0)
52 {
53  messagesToSlaves_.resize( nSlaves_ );
54  for ( UInt_t iSlave(0); iSlave < nSlaves_; ++iSlave ) {
55  messagesToSlaves_[iSlave] = new TMessage();
56  }
57 }
58 
60 {
61  delete socketMonitor_; socketMonitor_ = 0;
62  TString msgStr("Finish");
63  TMessage message( kMESS_STRING );
64  message.WriteTString(msgStr);
65  for ( std::vector<TSocket*>::iterator iter = sSlaves_.begin(); iter != sSlaves_.end(); ++iter ) {
66  (*iter)->Send(message);
67  (*iter)->Close();
68  delete (*iter);
69  }
70  sSlaves_.clear();
71  for ( std::vector<LauParameter*>::iterator iter = params_.begin(); iter != params_.end(); ++iter ) {
72  delete *iter;
73  }
74  params_.clear();
75  for ( std::vector<Double_t*>::iterator iter = vectorPar_.begin(); iter != vectorPar_.end(); ++iter ) {
76  delete[] (*iter);
77  }
78  vectorPar_.clear();
80  for ( std::vector<TMessage*>::iterator iter = messagesToSlaves_.begin(); iter != messagesToSlaves_.end(); ++iter ) {
81  delete (*iter);
82  }
83  messagesToSlaves_.clear();
84  delete fitNtuple_;
85 }
86 
88 {
89  if ( socketMonitor_ != 0 ) {
90  std::cerr << "ERROR in LauSimFitMaster::initSockets : Sockets already initialised." << std::endl;
91  return;
92  }
93 
94  //initialise socket connection, then accept a connection and return a full-duplex communication socket.
95  socketMonitor_ = new TMonitor();
96 
97  TServerSocket *ss = new TServerSocket( reqPort_, kFALSE );
98  UInt_t actual_port = ss->GetLocalPort();
99 
100  std::cout << "INFO in LauSimFitMaster::initSockets : Waiting for connection with " << nSlaves_ << " workers on port " << actual_port << std::endl;
101 
102  sSlaves_.resize(nSlaves_);
103  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
104  sSlaves_[iSlave] = ss->Accept();
105  std::cout << " : Added slave " << iSlave << std::endl;
106  }
107 
108  // tell the clients to start
109  std::cout << "INFO in LauSimFitMaster::initSockets : Initialising slaves" << std::endl;
110  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
111 
112  TMessage message( kMESS_ANY );
113  message.WriteUInt(iSlave);
114  message.WriteUInt(nSlaves_);
115 
116  sSlaves_[iSlave]->Send(message);
117 
118  socketMonitor_->Add(sSlaves_[iSlave]);
119  }
120  std::cout << " : Now start fit\n" << std::endl;
121 
122  ss->Close();
123  delete ss;
124 }
125 
126 /*
127  * 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
128  * THIS INCREASES THE GENERALITY OF THE CODE, I.E. THERE IS NO NEED FOR THE SLAVES TO KNOW ANY LAURA++ CLASS BUT THIS ONE, BUT MAKES IT RATHER MORE DENSE
129  * 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
130  *
131 void LauSimFitMaster::getParametersFromSlavesFirstTime()
132 {
133  slaveIndices_.resize( nSlaves_ );
134 
135  TSocket* sActive(0);
136  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
137  // Send a message to the slave, requesting the list of parameter names
138  TString msgStr = "Parameter Names";
139  TMessage message( kMESS_STRING );
140  message.WriteTString( msgStr );
141  sSlaves_[iSlave]->Send(message);
142 
143  // Wait to receive the response and check that it has come from the slave we just requested from
144  sActive = socketMonitor_->Select();
145  if ( sActive != sSlaves_[iSlave] ) {
146  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Received message from a different slave than expected!" << std::endl;
147  gSystem->Exit(1);
148  }
149 
150  // Read the object and extract the parameter names
151  sSlaves_[iSlave]->Recv( messageFromSlave_ );
152  TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
153  if ( ! objarray ) {
154  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameter names from slave" << std::endl;
155  gSystem->Exit(1);
156  }
157 
158  Int_t nPars = objarray->GetEntries();
159  for ( Int_t iPar(0); iPar < nPars; ++iPar ) {
160  TObjString* objstring = dynamic_cast<TObjString*>( (*objarray)[iPar] );
161  if ( ! objstring ) {
162  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameter names from slave" << std::endl;
163  gSystem->Exit(1);
164  }
165  TString parname = objstring->GetString();
166 
167  std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
168  if ( iter != parIndices_.end() ) {
169  UInt_t index = iter->second;
170  slaveIndices_[iSlave].push_back( index );
171  } else {
172  UInt_t index = parIndices_.size();
173  parIndices_.insert( std::make_pair( parname, index ) );
174  parNames_.insert( std::make_pair( index, parname ) );
175  slaveIndices_[iSlave].push_back( index );
176  }
177  }
178 
179  delete objarray; objarray = 0;
180  delete messageFromSlave_; messageFromSlave_ = 0;
181  }
182 
183  UInt_t nPars = parNames_.size();
184  parValues_.resize( nPars );
185 }
186 */
187 
189 {
190  if ( socketMonitor_ == 0 ) {
191  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlaves : Sockets not initialised." << std::endl;
192  return;
193  }
194 
195  if ( params_.empty() ) {
197  } else {
199  }
200 }
201 
203 {
204  TSocket* sActive(0);
205 
206  // Construct a message, requesting the list of parameter names
207  TString msgStr = "Send Parameters";
208  TMessage message( kMESS_STRING );
209  message.WriteTString( msgStr );
210 
211  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
212  // Send the message to the slave
213  sSlaves_[iSlave]->Send(message);
214 
215  // Wait to receive the response and check that it has come from the slave we just requested from
216  sActive = socketMonitor_->Select();
217  if ( sActive != sSlaves_[iSlave] ) {
218  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Received message from a different slave than expected!" << std::endl;
219  gSystem->Exit(1);
220  }
221 
222  // Read the object and extract the parameter names
223  sSlaves_[iSlave]->Recv( messageFromSlave_ );
224  TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
225  if ( ! objarray ) {
226  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Error reading parameter names from slave" << std::endl;
227  gSystem->Exit(1);
228  }
229 
230  // We want to auto-delete the supplied parameters since we only copy their values in this case
231  objarray->SetOwner(kTRUE);
232 
233  const UInt_t nPars = objarray->GetEntries();
234  if ( nPars != slaveIndices_[iSlave].size() ) {
235  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Unexpected number of parameters received from slave" << std::endl;
236  gSystem->Exit(1);
237  }
238 
239  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
240  LauParameter* parameter = dynamic_cast<LauParameter*>( (*objarray)[iPar] );
241  if ( ! parameter ) {
242  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Error reading parameter from slave" << std::endl;
243  gSystem->Exit(1);
244  }
245 
246  TString parname = parameter->name();
247  Double_t parvalue = parameter->initValue();
248 
249  std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
250  if ( iter == parIndices_.end() ) {
251  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Unexpected parameter name received from slave" << std::endl;
252  gSystem->Exit(1);
253  }
254 
255  const UInt_t index = iter->second;
256  if ( slaveIndices_[iSlave][iPar] != index ) {
257  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Unexpected parameter received from slave" << std::endl;
258  gSystem->Exit(1);
259  }
260 
261  params_[index]->initValue( parvalue );
262  parValues_[index] = parvalue;
263  vectorPar_[iSlave][iPar] = parvalue;
264  this->checkParameter( parameter, index );
265  }
266 
267  delete objarray; objarray = 0;
269  }
270 }
271 
273 {
274  slaveIndices_.resize( nSlaves_ );
275  vectorPar_.resize( nSlaves_ );
276  vectorRes_.resize( nSlaves_ );
277 
278  TSocket* sActive(0);
279 
280  // Construct a message, requesting the list of parameter names
281  TString msgStr = "Send Parameters";
282  TMessage message( kMESS_STRING );
283  message.WriteTString( msgStr );
284 
285  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
286  // Send the message to the slave
287  sSlaves_[iSlave]->Send(message);
288 
289  // Wait to receive the response and check that it has come from the slave we just requested from
290  sActive = socketMonitor_->Select();
291  if ( sActive != sSlaves_[iSlave] ) {
292  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Received message from a different slave than expected!" << std::endl;
293  gSystem->Exit(1);
294  }
295 
296  // Read the object and extract the parameter names
297  sSlaves_[iSlave]->Recv( messageFromSlave_ );
298  TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
299  if ( ! objarray ) {
300  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameters from slave" << std::endl;
301  gSystem->Exit(1);
302  }
303 
304  const UInt_t nPars = objarray->GetEntries();
305 
306  vectorPar_[iSlave] = new Double_t[nPars];
307 
308  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
309  LauParameter* parameter = dynamic_cast<LauParameter*>( (*objarray)[iPar] );
310  if ( ! parameter ) {
311  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameter from slave" << std::endl;
312  gSystem->Exit(1);
313  }
314 
315  TString parname = parameter->name();
316  Double_t parvalue = parameter->initValue();
317 
318  std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
319  if ( iter != parIndices_.end() ) {
320  UInt_t index = iter->second;
321  slaveIndices_[iSlave].push_back( index );
322  this->checkParameter( parameter, index );
323  } else {
324  UInt_t index = parIndices_.size();
325  parIndices_.insert( std::make_pair( parname, index ) );
326  parNames_.insert( std::make_pair( index, parname ) );
327  slaveIndices_[iSlave].push_back( index );
328  params_.push_back( parameter );
329  parValues_.push_back( parvalue );
330  }
331  vectorPar_[iSlave][iPar] = parvalue;
332  }
333 
334  delete objarray; objarray = 0;
336  }
337 
338  nParams_ = params_.size();
339 }
340 
342 {
343  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
344  const std::vector<UInt_t>& indices = slaveIndices_[iSlave];
345 
346  std::cout << "INFO in LauSimFitMaster::printParInfo : Slave " << iSlave << " has the following parameters:\n";
347  for ( std::vector<UInt_t>::const_iterator iter = indices.begin(); iter != indices.end(); ++iter ) {
348  const TString& parName = parNames_.find(*iter)->second;
349  Double_t parValue = parValues_[*iter];
350  const LauParameter* par = params_[*iter];
351  if ( par->name() != parName || par->initValue() != parValue ) {
352  std::cerr << "ERROR in LauSimFitMaster::printParInfo : Discrepancy in parameter name and value records, this is very strange!!" << std::endl;
353  }
354 
355  std::cout << " : " << parName << " = " << parValue << " and has index " << *iter << "\n";
356  }
357 
358  std::cout << std::endl;
359  }
360 
361  std::cout << "INFO in LauSimFitMaster::printParInfo : " << "There are " << nParams_ << " parameters in total" << std::endl;
362 }
363 
364 void LauSimFitMaster::checkParameter( const LauParameter* param, UInt_t index ) const
365 {
366  const LauParameter* storedPar = params_[index];
367 
368  TString parName = storedPar->name();
369  if ( param->name() != parName ) {
370  std::cerr << "ERROR in LauSimFitMaster::checkParameter : Parameter name is different!! This shouldn't happen!!" << std::endl;
371  }
372  if ( param->initValue() != storedPar->initValue() ) {
373  std::cerr << "WARNING in LauSimFitMaster::checkParameter : Initial value for parameter " << parName << " is different, will use the value first set: " << storedPar->initValue() << std::endl;
374  }
375  if ( param->minValue() != storedPar->minValue() ) {
376  std::cerr << "WARNING in LauSimFitMaster::checkParameter : Minimum allowed value for parameter " << parName << " is different, will use the value first set: " << storedPar->minValue() << std::endl;
377  }
378  if ( param->maxValue() != storedPar->maxValue() ) {
379  std::cerr << "WARNING in LauSimFitMaster::checkParameter : Maximum allowed value for parameter " << parName << " is different, will use the value first set: " << storedPar->maxValue() << std::endl;
380  }
381  if ( param->fixed() != storedPar->fixed() ) {
382  std::cerr << "WARNING in LauSimFitMaster::checkParameter : Fixed/floating property of parameter " << parName << " is different, will use the value first set: " << (storedPar->fixed() ? "fixed" : "floating") << std::endl;
383  }
384  if ( param->firstStage() != storedPar->firstStage() ) {
385  std::cerr << "WARNING in LauSimFitMaster::checkParameter : First stage property of parameter " << parName << " is different, will use the value first set: " << (storedPar->firstStage() ? "true" : "false") << std::endl;
386  }
387  if ( param->secondStage() != storedPar->secondStage() ) {
388  std::cerr << "WARNING in LauSimFitMaster::checkParameter : Second stage property of parameter " << parName << " is different, will use the value first set: " << (storedPar->secondStage() ? "true" : "false") << std::endl;
389  }
390 }
391 
393 {
394  this->initSockets();
395 }
396 
397 void LauSimFitMaster::runSimFit( const TString& fitNtupleFileName, UInt_t nExpt, UInt_t firstExpt, Bool_t useAsymmErrors, Bool_t twoStageFit )
398 {
399  // Routine to perform the total fit.
400 
401  // First, initialise
402  this->initialise();
403 
404  std::cout << "INFO in LauSimFitMaster::runSimFit : First experiment = " << firstExpt << std::endl;
405  std::cout << "INFO in LauSimFitMaster::runSimFit : Number of experiments = " << nExpt << std::endl;
406 
407  // Start the cumulative timer
408  cumulTimer_.Start();
409 
410  numberOKFits_ = 0, numberBadFits_ = 0;
411  fitStatus_ = -1;
412 
413  // Create and setup the fit results ntuple
414  std::cout << "INFO in LauSimFitMaster::runSimFit : Creating fit ntuple." << std::endl;
415  if (fitNtuple_ != 0) {delete fitNtuple_; fitNtuple_ = 0;}
416  fitNtuple_ = new LauFitNtuple(fitNtupleFileName, useAsymmErrors);
417 
418  // Loop over the number of experiments
419  for (iExpt_ = firstExpt; iExpt_ < (firstExpt+nExpt); ++iExpt_) {
420 
421  // Start the timer to see how long each fit takes
422  timer_.Start();
423 
424  // Instruct the slaves to read the data for this experiment
425  Bool_t readOK = this->readData();
426  if ( ! readOK ) {
427  std::cerr << "ERROR in LauSimFitMaster::runSimFit : One or more slaves reported problems with reading data for experiment " << iExpt_ << ", skipping..." << std::endl;
428  timer_.Stop();
429  continue;
430  }
431 
432  // Instruct the slaves to perform the caching
433  this->cacheInputData();
434 
435  // Do the fit
436  this->fitExpt( useAsymmErrors, twoStageFit );
437 
438  // Stop the timer and see how long the program took so far
439  timer_.Stop();
440  timer_.Print();
441 
442  // Instruct the slaves to finalise the results
443  this->finalise();
444 
445  // Keep track of how many fits succeeded or failed
446  if (fitStatus_ == 3) {
447  ++numberOKFits_;
448  } else {
449  ++numberBadFits_;
450  }
451  }
452 
453  // Print out total timing info.
454  std::cout << "INFO in LauSimFitMaster::runSimFit : Cumulative timing:" << std::endl;
455  cumulTimer_.Stop();
456  cumulTimer_.Print();
457 
458  // Print out stats on OK fits.
459  std::cout << "INFO in LauSimFitMaster::runSimFit : Number of OK Fits = " << numberOKFits_ << std::endl;
460  std::cout << "INFO in LauSimFitMaster::runSimFit : Number of Failed Fits = " << numberBadFits_ << std::endl;
461  Double_t fitEff(0.0);
462  if (nExpt != 0) {fitEff = numberOKFits_/(1.0*nExpt);}
463  std::cout << "INFO in LauSimFitMaster::runSimFit : Fit efficiency = " << fitEff*100.0 << "%." << std::endl;
464 
465  // Instruct the slaves to write out any fit results (ntuples etc...).
466  this->writeOutResults();
467 }
468 
470 {
471  if ( socketMonitor_ == 0 ) {
472  std::cerr << "ERROR in LauSimFitMaster::readData : Sockets not initialised." << std::endl;
473  return kFALSE;
474  }
475 
476  // Construct a message, requesting to read the data for the given experiment
477  TString msgStr("Read Expt");
478  TMessage message( kMESS_STRING );
479  message.WriteTString( msgStr );
480  message.WriteUInt( iExpt_ );
481 
482  // Send the message to the slaves
483  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
484  sSlaves_[iSlave]->Send(message);
485  }
486 
487  TSocket* sActive(0);
488  UInt_t responsesReceived(0);
489  Bool_t ok(kTRUE);
490  while ( responsesReceived != nSlaves_ ) {
491 
492  // Get the next queued response
493  sActive = socketMonitor_->Select();
494 
495  // Extract from the message the ID of the slave and the number of events read
496  sActive->Recv( messageFromSlave_ );
497  UInt_t iSlave(0);
498  UInt_t nEvents(0);
499  messageFromSlave_->ReadUInt( iSlave );
500  messageFromSlave_->ReadUInt( nEvents );
501 
502  if ( nEvents <= 0 ) {
503  std::cerr << "ERROR in LauSimFitMaster::readData : Slave " << iSlave << " reports no events found for experiment " << iExpt_ << std::endl;
504  ok = kFALSE;
505  } else {
506  std::cerr << "INFO in LauSimFitMaster::readData : Slave " << iSlave << " reports " << nEvents << " events found for experiment " << iExpt_ << std::endl;
507  }
508 
509  ++responsesReceived;
510  }
511 
512  return ok;
513 }
514 
516 {
517  if ( socketMonitor_ == 0 ) {
518  std::cerr << "ERROR in LauSimFitMaster::cacheInputData : Sockets not initialised." << std::endl;
519  return kFALSE;
520  }
521 
522  // Construct a message, requesting it to read the data for the given experiment
523  TString msgStr("Cache");
524  TMessage message( kMESS_STRING );
525  message.WriteTString( msgStr );
526 
527  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
528  // Send the message to the slave
529  sSlaves_[iSlave]->Send(message);
530  }
531 
532  TSocket* sActive(0);
533  UInt_t responsesReceived(0);
534  Bool_t allOK(kTRUE);
535  while ( responsesReceived != nSlaves_ ) {
536 
537  // Get the next queued response
538  sActive = socketMonitor_->Select();
539 
540  // Extract from the message the ID of the slave and the success/failure flag
541  sActive->Recv( messageFromSlave_ );
542  UInt_t iSlave(0);
543  Bool_t ok(kTRUE);
544  messageFromSlave_->ReadUInt( iSlave );
545  messageFromSlave_->ReadBool( ok );
546 
547  if ( ! ok ) {
548  std::cerr << "ERROR in LauSimFitMaster::cacheInputData : Slave " << iSlave << " reports an error performing caching" << std::endl;
549  allOK = kFALSE;
550  }
551 
552  ++responsesReceived;
553  }
554 
555  return allOK;
556 }
557 
559 {
560  this->getParametersFromSlaves();
561  this->printParInfo();
562 }
563 
564 void LauSimFitMaster::fitExpt( Bool_t useAsymmErrors, Bool_t twoStageFit )
565 {
566  // Routine to perform the actual fit for the given experiment
567 
568  // Instruct the salves to update initial fit parameters if required (e.g. if using random numbers).
569  this->checkInitFitParams();
570 
571  // Initialise the fitter
572  LauFitter::fitter()->useAsymmFitErrors( useAsymmErrors );
573  LauFitter::fitter()->twoStageFit( twoStageFit );
575 
578 
579  // Now ready for minimisation step
580  std::cout << "\nINFO in LauSimFitMaster::fitExpt : Start minimisation...\n";
581  std::pair<Int_t,Double_t> fitResult = LauFitter::fitter()->minimise();
582 
583  fitStatus_ = fitResult.first;
584  NLL_ = fitResult.second;
585 
586  // If we're doing a two stage fit we can now release (i.e. float)
587  // the 2nd stage parameters and re-fit
588  if ( twoStageFit ) {
589 
590  if ( fitStatus_ != 3 ) {
591  std::cerr << "ERROR in LauSimFitMaster:fitExpt : Not running second stage fit since first stage failed." << std::endl;
593  } else {
598  fitResult = LauFitter::fitter()->minimise();
599  }
600  }
601 
602  fitStatus_ = fitResult.first;
603  NLL_ = fitResult.second;
604  const TMatrixD& covMat = LauFitter::fitter()->covarianceMatrix();
605  covMatrix_.Clear();
606  covMatrix_.ResizeTo( covMat.GetNrows(), covMat.GetNcols() );
607  covMatrix_.SetMatrixArray( covMat.GetMatrixArray() );
608 
609  // Store the final fit results and errors into protected internal vectors that
610  // all sub-classes can use within their own finalFitResults implementation
611  // used below (e.g. putting them into an ntuple in a root file)
614 }
615 
616 void LauSimFitMaster::setParsFromMinuit(Double_t* par, Int_t npar)
617 {
618  // This function sets the internal parameters based on the values
619  // that Minuit is using when trying to minimise the total likelihood function.
620 
621  // MINOS reports different numbers of free parameters depending on the
622  // situation, so disable this check
623  if ( ! withinAsymErrorCalc_ ) {
624  if (static_cast<UInt_t>(npar) != nFreeParams_) {
625  std::cerr << "ERROR in LauSimFitMaster::setParsFromMinuit : Unexpected number of free parameters: " << npar << ".\n";
626  std::cerr << " Expected: " << nFreeParams_ << ".\n" << std::endl;
627  gSystem->Exit(EXIT_FAILURE);
628  }
629  }
630 
631  // Despite npar being the number of free parameters
632  // the par array actually contains all the parameters,
633  // free and floating...
634  // Update all the parameters with their new values.
635  for (UInt_t i(0); i<nParams_; ++i) {
636  parValues_[i] = par[i];
637  }
638 }
639 
641 {
642  if ( socketMonitor_ == 0 ) {
643  std::cerr << "ERROR in LauSimFitMaster::getTotNegLogLikelihood : Sockets not initialised." << std::endl;
644  return 0.0;
645  }
646 
647  // Send current values of the parameters to the slaves.
648  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
649 
650  std::vector<UInt_t>& indices = slaveIndices_[iSlave];
651  UInt_t nPars = indices.size();
652  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
653  vectorPar_[iSlave][iPar] = parValues_[ indices[iPar] ];
654  }
655 
656  TMessage* message = messagesToSlaves_[iSlave];
657  message->Reset( kMESS_ANY );
658  message->WriteUInt( nPars );
659  message->WriteFastArray( vectorPar_[iSlave], nPars );
660 
661  sSlaves_[iSlave]->Send(*message);
662  }
663 
664  Double_t negLogLike(0.0);
665  TSocket *sActive(0);
666  UInt_t responsesReceived(0);
667  while ( responsesReceived != nSlaves_ ) {
668 
669  sActive = socketMonitor_->Select();
670  sActive->Recv(messageFromSlave_);
671 
672  messageFromSlave_->ReadDouble( vectorRes_[responsesReceived] );
673 
674  negLogLike += vectorRes_[responsesReceived];
675 
676  ++responsesReceived;
677  }
678 
679  return negLogLike;
680 }
681 
683 {
684  if ( socketMonitor_ == 0 ) {
685  std::cerr << "ERROR in LauSimFitMaster::finalise : Sockets not initialised." << std::endl;
686  return kFALSE;
687  }
688 
689  // Prepare the covariance matrices
690  covMatrices_.resize( nSlaves_ );
691 
692  LauParamFixed pred;
693 
694  std::map<UInt_t,UInt_t> freeParIndices;
695 
696  UInt_t counter(0);
697  for ( UInt_t iPar(0); iPar < nParams_; ++iPar ) {
698  const LauParameter* par = params_[iPar];
699  if ( ! pred(par) ) {
700  freeParIndices.insert( std::make_pair(iPar,counter) );
701  ++counter;
702  }
703  }
704 
705  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
706  const UInt_t nPar = slaveIndices_[iSlave].size();
707 
708  std::vector<UInt_t> freeIndices;
709  freeIndices.reserve( nPar );
710 
711  for ( UInt_t iPar(0); iPar < nPar; ++iPar ) {
712  UInt_t index = slaveIndices_[iSlave][iPar];
713  std::map<UInt_t,UInt_t>::iterator freeIter = freeParIndices.find(index);
714  if ( freeIter == freeParIndices.end() ) {
715  continue;
716  }
717  UInt_t freeIndex = freeIter->second;
718  freeIndices.push_back( freeIndex );
719  }
720 
721  const UInt_t nFreePars = freeIndices.size();
722  TMatrixD& covMat = covMatrices_[iSlave];
723  covMat.ResizeTo( nFreePars, nFreePars );
724 
725  for ( UInt_t iPar(0); iPar < nFreePars; ++iPar ) {
726  for ( UInt_t jPar(0); jPar < nFreePars; ++jPar ) {
727  UInt_t i = freeIndices[iPar];
728  UInt_t j = freeIndices[jPar];
729  covMat( iPar, jPar ) = covMatrix_( i, j );
730  }
731  }
732  }
733 
734  // The array to hold the parameters
735  TObjArray array;
736 
737  // Send messages to all slaves containing the final parameters and fit status, NLL
738  // TODO - at present we lose the information on the correlations between the parameters that are unique to each slave
739  // - so should we store the full correlation matrix in an ntuple? along with all the parameters?
740  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
741 
742  array.Clear();
743 
744  std::vector<UInt_t>& indices = slaveIndices_[iSlave];
745  UInt_t nPars = indices.size();
746  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
747  array.Add( params_[ indices[iPar] ] );
748  }
749 
750  TMatrixD& covMat = covMatrices_[iSlave];
751 
752  TMessage* message = messagesToSlaves_[iSlave];
753  message->Reset( kMESS_OBJECT );
754  message->WriteInt( fitStatus_ );
755  message->WriteDouble( NLL_ );
756  message->WriteObject( &array );
757  message->WriteObject( &covMat );
758 
759  sSlaves_[iSlave]->Send(*message);
760  }
761 
762  TSocket *sActive(0);
763  UInt_t responsesReceived(0);
764  Bool_t allOK(kTRUE);
765  while ( responsesReceived != nSlaves_ ) {
766 
767  // Get the next queued response
768  sActive = socketMonitor_->Select();
769 
770  // Extract from the message the ID of the slave and the number of events read
771  sActive->Recv( messageFromSlave_ );
772  UInt_t iSlave(0);
773  Bool_t ok(kTRUE);
774  messageFromSlave_->ReadUInt( iSlave );
775  messageFromSlave_->ReadBool( ok );
776 
777  if ( ok ) {
778  TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
779  if ( ! objarray ) {
780  std::cerr << "ERROR in LauSimFitMaster::finalise : Error reading finalised parameters from slave" << std::endl;
781  allOK = kFALSE;
782  } else {
783  // We want to auto-delete the supplied parameters since we only copy their values in this case
784  objarray->SetOwner(kTRUE);
785 
786  const UInt_t nPars = objarray->GetEntries();
787  if ( nPars != slaveIndices_[iSlave].size() ) {
788  std::cerr << "ERROR in LauSimFitMaster::finalise : Unexpected number of finalised parameters received from slave" << std::endl;
789  allOK = kFALSE;
790  } else {
791  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
792  LauParameter* parameter = dynamic_cast<LauParameter*>( (*objarray)[iPar] );
793  if ( ! parameter ) {
794  std::cerr << "ERROR in LauSimFitMaster::finalise : Error reading parameter from slave" << std::endl;
795  allOK = kFALSE;
796  continue;
797  }
798 
799  TString parname = parameter->name();
800 
801  std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
802  if ( iter == parIndices_.end() ) {
803  std::cerr << "ERROR in LauSimFitMaster::finalise : Unexpected parameter name received from slave" << std::endl;
804  allOK = kFALSE;
805  continue;
806  }
807 
808  const UInt_t index = iter->second;
809  if ( slaveIndices_[iSlave][iPar] != index ) {
810  std::cerr << "ERROR in LauSimFitMaster::finalise : Unexpected parameter received from slave" << std::endl;
811  allOK = kFALSE;
812  continue;
813  }
814 
815  Double_t parvalue = parameter->value();
816  params_[index]->value( parvalue );
817  parValues_[index] = parvalue;
818  vectorPar_[iSlave][iPar] = parvalue;
819  }
820  }
821  delete objarray;
822  }
823  } else {
824  std::cerr << "ERROR in LauSimFitMaster::finalise : Slave " << iSlave << " reports an error performing finalisation" << std::endl;
825  allOK = kFALSE;
826  }
827 
828  ++responsesReceived;
829  }
830 
831  // Fill our ntuple as well
832  if ( fitNtuple_ != 0 ) {
833  for ( std::vector<LauParameter*>::iterator iter = params_.begin(); iter != params_.end(); ++iter ) {
834  if (!(*iter)->fixed()) {
835  (*iter)->updatePull();
836  }
837  }
838  std::vector<LauParameter> extraVars;
839  fitNtuple_->storeParsAndErrors( params_, extraVars );
842  }
843 
844  return allOK;
845 }
846 
848 {
849  if ( socketMonitor_ == 0 ) {
850  std::cerr << "ERROR in LauSimFitMaster::writeOutResults : Sockets not initialised." << std::endl;
851  return kFALSE;
852  }
853 
854  // Construct a message, requesting to write out the fit results
855  TString msgStr("Write Results");
856  TMessage message( kMESS_STRING );
857  message.WriteTString( msgStr );
858 
859  // Send the message to the slaves
860  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
861  sSlaves_[iSlave]->Send(message);
862  }
863 
864  TSocket *sActive(0);
865  UInt_t responsesReceived(0);
866  Bool_t allOK(kTRUE);
867  while ( responsesReceived != nSlaves_ ) {
868 
869  // Get the next queued response
870  sActive = socketMonitor_->Select();
871 
872  // Extract from the message the ID of the slave and the number of events read
873  sActive->Recv( messageFromSlave_ );
874  UInt_t iSlave(0);
875  Bool_t ok(kTRUE);
876  messageFromSlave_->ReadUInt( iSlave );
877  messageFromSlave_->ReadBool( ok );
878 
879  if ( ! ok ) {
880  std::cerr << "ERROR in LauSimFitMaster::writeOutResults : Slave " << iSlave << " reports an error performing finalisation" << std::endl;
881  allOK = kFALSE;
882  }
883 
884  ++responsesReceived;
885  }
886 
887  // Write out our ntuple as well
888  if (fitNtuple_ != 0) {
890  }
891 
892  return allOK;
893 }
894 
TMonitor * socketMonitor_
Parallel setup monitor.
std::vector< std::vector< UInt_t > > slaveIndices_
Lists of indices for each slave.
Bool_t readData()
Instruct the slaves to read the input data for the given experiment.
void getParametersFromSlaves()
Determine/update the parameter initial values from all slaves.
Bool_t fixed() const
Check whether the parameter is fixed or floated.
virtual void updateParameters()=0
Update the values and errors of the parameters based on the fit minimum.
TStopwatch timer_
The fit timer.
std::map< UInt_t, TString > parNames_
Reverse map of index in the values vector to parameter names.
File containing declaration of LauParamFixed class.
LauFitNtuple * fitNtuple_
The fit results ntuple.
Double_t maxValue() const
The maximum value allowed for the parameter.
Bool_t writeOutResults()
Instruct the slaves to write out the fit results.
std::vector< LauParameter * > params_
Parameters.
const TString & name() const
The parameter name.
virtual Bool_t twoStageFit() const =0
Determine whether the two-stage fit is enabled.
File containing declaration of LauAbsFitter class.
Double_t minValue() const
The minimum value allowed for the parameter.
void getParametersFromSlavesFirstTime()
Determine the parameter names and initial values from all slaves.
UInt_t nParams_
The number of fit parameters.
std::vector< TSocket * > sSlaves_
Sockets for each of the slaves.
const UInt_t reqPort_
The requested port.
void updateFitNtuple()
Update the fit ntuple.
std::vector< Double_t > vectorRes_
Likelihood values returned from the slaves.
std::map< TString, UInt_t > parIndices_
Map of parameter names to index in the values vector.
virtual void initialise(LauFitObject *fitObj, const std::vector< LauParameter * > &parameters)=0
Initialise the fitter, setting the information on the parameters.
virtual std::pair< Int_t, Double_t > minimise()=0
Perform the minimisation of the fit function.
Bool_t withinAsymErrorCalc_
Flag to indicate if the asymmetric error calculation (e.g. MINOS) is currently running.
File containing LauFitter namespace.
void initSockets()
Initialise socket connections for the slaves.
UInt_t numberBadFits_
The number of fit failures.
Bool_t firstStage() const
Check whether the parameter should be floated only in the first stage of a two stage fit...
void writeOutFitResults()
Write out fit results.
TMatrixD covMatrix_
The fit covariance matrix.
void checkInitFitParams()
Instruct the slaves to update the initial fit parameter values, if required.
File containing declaration of LauParameter class.
Predicate to allow counting of the number of fixed parameters.
std::vector< Double_t * > vectorPar_
Parameter values to send to the slaves.
Bool_t secondStage() const
Check whether the parameter should be floated only in the second stage of a two stage fit...
void updateParametersFromSlaves()
Update and verify the parameter initial values from all slaves.
virtual void fixFirstStageParameters()=0
Fix parameters marked as &quot;first stage&quot;.
std::vector< TMatrixD > covMatrices_
The covariance sub-matrices for each slave.
virtual UInt_t nParameters() const =0
Get the total number of fit parameters.
std::vector< TMessage * > messagesToSlaves_
Messages to slaves.
virtual void releaseSecondStageParameters()=0
Release parameters marked as &quot;second stage&quot;.
void storeCorrMatrix(UInt_t iExpt, Double_t NLL, Int_t fitStatus, const TMatrixD &covMatrix)
Store the correlation matrix and other fit information.
Definition: LauFitNtuple.cc:63
void checkParameter(const LauParameter *param, UInt_t index) const
Check for compatibility between two same-named parameters, which should therefore be identical...
Class for defining the fit parameter objects.
Definition: LauParameter.hh:32
void fitExpt(Bool_t useAsymmErrors, Bool_t twoStageFit)
Perform the fit for the current experiment.
Int_t fitStatus_
The status of the current fit.
The master process for simultaneous/combined fits.
static LauAbsFitter * fitter()
Method that provides access to the singleton fitter.
Definition: LauFitter.cc:34
Bool_t finalise()
Return the final parameters to the slaves and instruct them to perform their finalisation.
TStopwatch cumulTimer_
The total fit timer.
virtual UInt_t nFreeParameters() const =0
Get the number of floating fit parameters.
virtual void setParsFromMinuit(Double_t *par, Int_t npar)
This function sets the parameter values from Minuit.
std::vector< Double_t > parValues_
Parameter values.
const UInt_t nSlaves_
The number of slaves.
File containing declaration of LauSimFitMaster class.
void storeParsAndErrors(const std::vector< LauParameter * > &fitVars, const std::vector< LauParameter > &extraVars)
Store parameters and their errors.
void runSimFit(const TString &fitNtupleFileName, UInt_t nExpt, UInt_t firstExpt=0, Bool_t useAsymmErrors=kFALSE, Bool_t twoStageFit=kFALSE)
Run the fit.
Double_t initValue() const
The initial value of the parameter.
UInt_t nFreeParams_
The number of free fit parameters.
TMessage * messageFromSlave_
Message from slaves to the master.
virtual Bool_t useAsymmFitErrors() const =0
Determine whether calculation of asymmetric errors is enabled.
Bool_t cacheInputData()
Instruct the slaves to perform the caching.
Double_t NLL_
The negative log-likelihood.
Double_t value() const
The value of the parameter.
void initialise()
Initialise.
friend class LauFitNtuple
LauFitNtuple is a friend class.
virtual Double_t getTotNegLogLikelihood()
Calculate the new value of the negative log likelihood.
virtual void releaseFirstStageParameters()=0
Release parameters marked as &quot;first stage&quot;.
UInt_t iExpt_
The experiment number of the current fit.
virtual ~LauSimFitMaster()
Destructor.
UInt_t numberOKFits_
The number of successful fits.
File containing declaration of LauFitNtuple class.
virtual const TMatrixD & covarianceMatrix() const =0
Retrieve the fit covariance matrix.
void printParInfo() const
Print information on the parameters.