laura is hosted by Hepforge, IPPP Durham
Laura++  v2r2
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 "LauFormulaPar.hh"
31 #include "LauParameter.hh"
32 #include "LauParamFixed.hh"
33 #include "LauSimFitMaster.hh"
34 
35 
37 
38 
39 LauSimFitMaster::LauSimFitMaster( UInt_t numSlaves, UInt_t port ) :
40  nSlaves_(numSlaves),
41  reqPort_(port),
42  nParams_(0),
43  nFreeParams_(0),
44  withinAsymErrorCalc_(kFALSE),
45  numberOKFits_(0),
46  numberBadFits_(0),
47  fitStatus_(0),
48  iExpt_(0),
49  NLL_(0.0),
50  socketMonitor_(0),
51  messageFromSlave_(0),
52  fitNtuple_(0)
53 {
54  messagesToSlaves_.resize( nSlaves_ );
55  for ( UInt_t iSlave(0); iSlave < nSlaves_; ++iSlave ) {
56  messagesToSlaves_[iSlave] = new TMessage();
57  }
58 }
59 
61 {
62  delete socketMonitor_; socketMonitor_ = 0;
63 
64  // Tell all slaves that they are finished and delete corresponding socket
65  TString msgStr("Finish");
66  TMessage message( kMESS_STRING );
67  message.WriteTString(msgStr);
68  for ( std::vector<TSocket*>::iterator iter = sSlaves_.begin(); iter != sSlaves_.end(); ++iter ) {
69  (*iter)->Send(message);
70  (*iter)->Close();
71  delete (*iter);
72  }
73  sSlaves_.clear();
74 
75  // Remove the components created to apply constraints to fit parameters
76  for (std::vector<LauAbsRValue*>::iterator iter = conVars_.begin(); iter != conVars_.end(); ++iter){
77  if ( !(*iter)->isLValue() ){
78  delete (*iter);
79  (*iter) = 0;
80  }
81  }
82  conVars_.clear();
83 
84  // Remove all fit parameters
85  for ( std::vector<LauParameter*>::iterator iter = params_.begin(); iter != params_.end(); ++iter ) {
86  delete *iter;
87  }
88  params_.clear();
89 
90  for ( std::vector<Double_t*>::iterator iter = vectorPar_.begin(); iter != vectorPar_.end(); ++iter ) {
91  delete[] (*iter);
92  }
93  vectorPar_.clear();
94 
96 
97  for ( std::vector<TMessage*>::iterator iter = messagesToSlaves_.begin(); iter != messagesToSlaves_.end(); ++iter ) {
98  delete (*iter);
99  }
100  messagesToSlaves_.clear();
101 
102  delete fitNtuple_;
103 }
104 
106 {
107  if ( socketMonitor_ != 0 ) {
108  std::cerr << "ERROR in LauSimFitMaster::initSockets : Sockets already initialised." << std::endl;
109  return;
110  }
111 
112  //initialise socket connection, then accept a connection and return a full-duplex communication socket.
113  socketMonitor_ = new TMonitor();
114 
115  TServerSocket *ss = new TServerSocket( reqPort_, kFALSE );
116  UInt_t actual_port = ss->GetLocalPort();
117 
118  std::cout << "INFO in LauSimFitMaster::initSockets : Waiting for connection with " << nSlaves_ << " workers on port " << actual_port << std::endl;
119 
120  sSlaves_.resize(nSlaves_);
121  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
122  sSlaves_[iSlave] = ss->Accept();
123  std::cout << " : Added slave " << iSlave << std::endl;
124  }
125 
126  // tell the clients to start
127  std::cout << "INFO in LauSimFitMaster::initSockets : Initialising slaves" << std::endl;
128  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
129 
130  TMessage message( kMESS_ANY );
131  message.WriteUInt(iSlave);
132  message.WriteUInt(nSlaves_);
133 
134  sSlaves_[iSlave]->Send(message);
135 
136  socketMonitor_->Add(sSlaves_[iSlave]);
137  }
138  std::cout << " : Now start fit\n" << std::endl;
139 
140  ss->Close();
141  delete ss;
142 }
143 
144 /*
145  * 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
146  * 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
147  * 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
148  *
149 void LauSimFitMaster::getParametersFromSlavesFirstTime()
150 {
151  slaveIndices_.resize( nSlaves_ );
152 
153  TSocket* sActive(0);
154  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
155  // Send a message to the slave, requesting the list of parameter names
156  TString msgStr = "Parameter Names";
157  TMessage message( kMESS_STRING );
158  message.WriteTString( msgStr );
159  sSlaves_[iSlave]->Send(message);
160 
161  // Wait to receive the response and check that it has come from the slave we just requested from
162  sActive = socketMonitor_->Select();
163  if ( sActive != sSlaves_[iSlave] ) {
164  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Received message from a different slave than expected!" << std::endl;
165  gSystem->Exit(1);
166  }
167 
168  // Read the object and extract the parameter names
169  sSlaves_[iSlave]->Recv( messageFromSlave_ );
170  TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
171  if ( ! objarray ) {
172  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameter names from slave" << std::endl;
173  gSystem->Exit(1);
174  }
175 
176  Int_t nPars = objarray->GetEntries();
177  for ( Int_t iPar(0); iPar < nPars; ++iPar ) {
178  TObjString* objstring = dynamic_cast<TObjString*>( (*objarray)[iPar] );
179  if ( ! objstring ) {
180  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameter names from slave" << std::endl;
181  gSystem->Exit(1);
182  }
183  TString parname = objstring->GetString();
184 
185  std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
186  if ( iter != parIndices_.end() ) {
187  UInt_t index = iter->second;
188  slaveIndices_[iSlave].push_back( index );
189  } else {
190  UInt_t index = parIndices_.size();
191  parIndices_.insert( std::make_pair( parname, index ) );
192  parNames_.insert( std::make_pair( index, parname ) );
193  slaveIndices_[iSlave].push_back( index );
194  }
195  }
196 
197  delete objarray; objarray = 0;
198  delete messageFromSlave_; messageFromSlave_ = 0;
199  }
200 
201  UInt_t nPars = parNames_.size();
202  parValues_.resize( nPars );
203 }
204 */
205 
207 {
208  if ( socketMonitor_ == 0 ) {
209  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlaves : Sockets not initialised." << std::endl;
210  return;
211  }
212 
213  if ( params_.empty() ) {
215 
216  // Add variables to Gaussian constrain to a list
217  this->addConParameters();
218  } else {
220  }
221 }
222 
224 {
225  TSocket* sActive(0);
226 
227  // Construct a message, requesting the list of parameter names
228  TString msgStr = "Send Parameters";
229  TMessage message( kMESS_STRING );
230  message.WriteTString( msgStr );
231 
232  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
233  // Send the message to the slave
234  sSlaves_[iSlave]->Send(message);
235 
236  // Wait to receive the response and check that it has come from the slave we just requested from
237  sActive = socketMonitor_->Select();
238  if ( sActive != sSlaves_[iSlave] ) {
239  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Received message from a different slave than expected!" << std::endl;
240  gSystem->Exit(1);
241  }
242 
243  // Read the object and extract the parameter names
244  sSlaves_[iSlave]->Recv( messageFromSlave_ );
245  TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
246  if ( ! objarray ) {
247  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Error reading parameter names from slave" << std::endl;
248  gSystem->Exit(1);
249  }
250 
251  // We want to auto-delete the supplied parameters since we only copy their values in this case
252  objarray->SetOwner(kTRUE);
253 
254  const UInt_t nPars = objarray->GetEntries();
255  if ( nPars != slaveIndices_[iSlave].size() ) {
256  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Unexpected number of parameters received from slave" << std::endl;
257  gSystem->Exit(1);
258  }
259 
260  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
261  LauParameter* parameter = dynamic_cast<LauParameter*>( (*objarray)[iPar] );
262  if ( ! parameter ) {
263  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Error reading parameter from slave" << std::endl;
264  gSystem->Exit(1);
265  }
266 
267  TString parname = parameter->name();
268  Double_t parvalue = parameter->initValue();
269 
270  std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
271  if ( iter == parIndices_.end() ) {
272  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Unexpected parameter name received from slave" << std::endl;
273  gSystem->Exit(1);
274  }
275 
276  const UInt_t index = iter->second;
277  if ( slaveIndices_[iSlave][iPar] != index ) {
278  std::cerr << "ERROR in LauSimFitMaster::updateParametersFromSlaves : Unexpected parameter received from slave" << std::endl;
279  gSystem->Exit(1);
280  }
281 
282  params_[index]->initValue( parvalue );
283  parValues_[index] = parvalue;
284  vectorPar_[iSlave][iPar] = parvalue;
285  this->checkParameter( parameter, index );
286  }
287 
288  delete objarray; objarray = 0;
290  }
291 }
292 
294 {
295  slaveIndices_.resize( nSlaves_ );
296  vectorPar_.resize( nSlaves_ );
297  vectorRes_.resize( nSlaves_ );
298 
299  TSocket* sActive(0);
300 
301  // Construct a message, requesting the list of parameter names
302  TString msgStr = "Send Parameters";
303  TMessage message( kMESS_STRING );
304  message.WriteTString( msgStr );
305 
306  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
307  // Send the message to the slave
308  sSlaves_[iSlave]->Send(message);
309 
310  // Wait to receive the response and check that it has come from the slave we just requested from
311  sActive = socketMonitor_->Select();
312  if ( sActive != sSlaves_[iSlave] ) {
313  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Received message from a different slave than expected!" << std::endl;
314  gSystem->Exit(1);
315  }
316 
317  // Read the object and extract the parameter names
318  sSlaves_[iSlave]->Recv( messageFromSlave_ );
319  TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
320  if ( ! objarray ) {
321  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameters from slave" << std::endl;
322  gSystem->Exit(1);
323  }
324 
325  const UInt_t nPars = objarray->GetEntries();
326 
327  vectorPar_[iSlave] = new Double_t[nPars];
328 
329  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
330  LauParameter* parameter = dynamic_cast<LauParameter*>( (*objarray)[iPar] );
331  if ( ! parameter ) {
332  std::cerr << "ERROR in LauSimFitMaster::getParametersFromSlavesFirstTime : Error reading parameter from slave" << std::endl;
333  gSystem->Exit(1);
334  }
335 
336  TString parname = parameter->name();
337  Double_t parvalue = parameter->initValue();
338 
339  std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
340  if ( iter != parIndices_.end() ) {
341  UInt_t index = iter->second;
342  slaveIndices_[iSlave].push_back( index );
343  this->checkParameter( parameter, index );
344  } else {
345  UInt_t index = parIndices_.size();
346  parIndices_.insert( std::make_pair( parname, index ) );
347  parNames_.insert( std::make_pair( index, parname ) );
348  slaveIndices_[iSlave].push_back( index );
349  params_.push_back( parameter );
350  parValues_.push_back( parvalue );
351  }
352  vectorPar_[iSlave][iPar] = parvalue;
353  }
354 
355  delete objarray; objarray = 0;
357  }
358 
359  nParams_ = params_.size();
360 }
361 
363 {
364  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
365  const std::vector<UInt_t>& indices = slaveIndices_[iSlave];
366 
367  std::cout << "INFO in LauSimFitMaster::printParInfo : Slave " << iSlave << " has the following parameters:\n";
368  for ( std::vector<UInt_t>::const_iterator iter = indices.begin(); iter != indices.end(); ++iter ) {
369  const TString& parName = parNames_.find(*iter)->second;
370  Double_t parValue = parValues_[*iter];
371  const LauParameter* par = params_[*iter];
372  if ( par->name() != parName || par->initValue() != parValue ) {
373  std::cerr << "ERROR in LauSimFitMaster::printParInfo : Discrepancy in parameter name and value records, this is very strange!!" << std::endl;
374  }
375 
376  std::cout << " : " << parName << " = " << parValue << " and has index " << *iter << "\n";
377  }
378 
379  std::cout << std::endl;
380  }
381 
382  std::cout << "INFO in LauSimFitMaster::printParInfo : " << "There are " << nParams_ << " parameters in total" << std::endl;
383 }
384 
385 void LauSimFitMaster::checkParameter( const LauParameter* param, UInt_t index ) const
386 {
387  const LauParameter* storedPar = params_[index];
388 
389  TString parName = storedPar->name();
390  if ( param->name() != parName ) {
391  std::cerr << "ERROR in LauSimFitMaster::checkParameter : Parameter name is different!! This shouldn't happen!!" << std::endl;
392  }
393  if ( param->initValue() != storedPar->initValue() ) {
394  std::cerr << "WARNING in LauSimFitMaster::checkParameter : Initial value for parameter " << parName << " is different, will use the value first set: " << storedPar->initValue() << std::endl;
395  }
396  if ( param->minValue() != storedPar->minValue() ) {
397  std::cerr << "WARNING in LauSimFitMaster::checkParameter : Minimum allowed value for parameter " << parName << " is different, will use the value first set: " << storedPar->minValue() << std::endl;
398  }
399  if ( param->maxValue() != storedPar->maxValue() ) {
400  std::cerr << "WARNING in LauSimFitMaster::checkParameter : Maximum allowed value for parameter " << parName << " is different, will use the value first set: " << storedPar->maxValue() << std::endl;
401  }
402  if ( param->fixed() != storedPar->fixed() ) {
403  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;
404  }
405  if ( param->firstStage() != storedPar->firstStage() ) {
406  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;
407  }
408  if ( param->secondStage() != storedPar->secondStage() ) {
409  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;
410  }
411 }
412 
414 {
415  this->initSockets();
416 }
417 
418 void LauSimFitMaster::runSimFit( const TString& fitNtupleFileName, UInt_t nExpt, UInt_t firstExpt, Bool_t useAsymmErrors, Bool_t twoStageFit )
419 {
420  // Routine to perform the total fit.
421 
422  // First, initialise
423  this->initialise();
424 
425  std::cout << "INFO in LauSimFitMaster::runSimFit : First experiment = " << firstExpt << std::endl;
426  std::cout << "INFO in LauSimFitMaster::runSimFit : Number of experiments = " << nExpt << std::endl;
427 
428  // Start the cumulative timer
429  cumulTimer_.Start();
430 
431  numberOKFits_ = 0, numberBadFits_ = 0;
432  fitStatus_ = -1;
433 
434  // Create and setup the fit results ntuple
435  std::cout << "INFO in LauSimFitMaster::runSimFit : Creating fit ntuple." << std::endl;
436  if (fitNtuple_ != 0) {delete fitNtuple_; fitNtuple_ = 0;}
437  fitNtuple_ = new LauFitNtuple(fitNtupleFileName, useAsymmErrors);
438 
439  // Loop over the number of experiments
440  for (iExpt_ = firstExpt; iExpt_ < (firstExpt+nExpt); ++iExpt_) {
441 
442  // Start the timer to see how long each fit takes
443  timer_.Start();
444 
445  // Instruct the slaves to read the data for this experiment
446  Bool_t readOK = this->readData();
447  if ( ! readOK ) {
448  std::cerr << "ERROR in LauSimFitMaster::runSimFit : One or more slaves reported problems with reading data for experiment " << iExpt_ << ", skipping..." << std::endl;
449  timer_.Stop();
450  continue;
451  }
452 
453  // Instruct the slaves to perform the caching
454  this->cacheInputData();
455 
456  // Do the fit
457  this->fitExpt( useAsymmErrors, twoStageFit );
458 
459  // Stop the timer and see how long the program took so far
460  timer_.Stop();
461  timer_.Print();
462 
463  // Instruct the slaves to finalise the results
464  this->finalise();
465 
466  // Keep track of how many fits succeeded or failed
467  if (fitStatus_ == 3) {
468  ++numberOKFits_;
469  } else {
470  ++numberBadFits_;
471  }
472  }
473 
474  // Print out total timing info.
475  std::cout << "INFO in LauSimFitMaster::runSimFit : Cumulative timing:" << std::endl;
476  cumulTimer_.Stop();
477  cumulTimer_.Print();
478 
479  // Print out stats on OK fits.
480  std::cout << "INFO in LauSimFitMaster::runSimFit : Number of OK Fits = " << numberOKFits_ << std::endl;
481  std::cout << "INFO in LauSimFitMaster::runSimFit : Number of Failed Fits = " << numberBadFits_ << std::endl;
482  Double_t fitEff(0.0);
483  if (nExpt != 0) {fitEff = numberOKFits_/(1.0*nExpt);}
484  std::cout << "INFO in LauSimFitMaster::runSimFit : Fit efficiency = " << fitEff*100.0 << "%." << std::endl;
485 
486  // Instruct the slaves to write out any fit results (ntuples etc...).
487  this->writeOutResults();
488 }
489 
491 {
492  if ( socketMonitor_ == 0 ) {
493  std::cerr << "ERROR in LauSimFitMaster::readData : Sockets not initialised." << std::endl;
494  return kFALSE;
495  }
496 
497  // Construct a message, requesting to read the data for the given experiment
498  TString msgStr("Read Expt");
499  TMessage message( kMESS_STRING );
500  message.WriteTString( msgStr );
501  message.WriteUInt( iExpt_ );
502 
503  // Send the message to the slaves
504  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
505  sSlaves_[iSlave]->Send(message);
506  }
507 
508  TSocket* sActive(0);
509  UInt_t responsesReceived(0);
510  Bool_t ok(kTRUE);
511  while ( responsesReceived != nSlaves_ ) {
512 
513  // Get the next queued response
514  sActive = socketMonitor_->Select();
515 
516  // Extract from the message the ID of the slave and the number of events read
517  sActive->Recv( messageFromSlave_ );
518  UInt_t iSlave(0);
519  UInt_t nEvents(0);
520  messageFromSlave_->ReadUInt( iSlave );
521  messageFromSlave_->ReadUInt( nEvents );
522 
523  if ( nEvents <= 0 ) {
524  std::cerr << "ERROR in LauSimFitMaster::readData : Slave " << iSlave << " reports no events found for experiment " << iExpt_ << std::endl;
525  ok = kFALSE;
526  } else {
527  std::cerr << "INFO in LauSimFitMaster::readData : Slave " << iSlave << " reports " << nEvents << " events found for experiment " << iExpt_ << std::endl;
528  }
529 
530  ++responsesReceived;
531  }
532 
533  return ok;
534 }
535 
537 {
538  if ( socketMonitor_ == 0 ) {
539  std::cerr << "ERROR in LauSimFitMaster::cacheInputData : Sockets not initialised." << std::endl;
540  return kFALSE;
541  }
542 
543  // Construct a message, requesting it to read the data for the given experiment
544  TString msgStr("Cache");
545  TMessage message( kMESS_STRING );
546  message.WriteTString( msgStr );
547 
548  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
549  // Send the message to the slave
550  sSlaves_[iSlave]->Send(message);
551  }
552 
553  TSocket* sActive(0);
554  UInt_t responsesReceived(0);
555  Bool_t allOK(kTRUE);
556  while ( responsesReceived != nSlaves_ ) {
557 
558  // Get the next queued response
559  sActive = socketMonitor_->Select();
560 
561  // Extract from the message the ID of the slave and the success/failure flag
562  sActive->Recv( messageFromSlave_ );
563  UInt_t iSlave(0);
564  Bool_t ok(kTRUE);
565  messageFromSlave_->ReadUInt( iSlave );
566  messageFromSlave_->ReadBool( ok );
567 
568  if ( ! ok ) {
569  std::cerr << "ERROR in LauSimFitMaster::cacheInputData : Slave " << iSlave << " reports an error performing caching" << std::endl;
570  allOK = kFALSE;
571  }
572 
573  ++responsesReceived;
574  }
575 
576  return allOK;
577 }
578 
580 {
581  this->getParametersFromSlaves();
582  this->printParInfo();
583 }
584 
585 void LauSimFitMaster::fitExpt( Bool_t useAsymmErrors, Bool_t twoStageFit )
586 {
587  // Routine to perform the actual fit for the given experiment
588 
589  // Instruct the salves to update initial fit parameters if required (e.g. if using random numbers).
590  this->checkInitFitParams();
591 
592  // Initialise the fitter
593  LauFitter::fitter()->useAsymmFitErrors( useAsymmErrors );
594  LauFitter::fitter()->twoStageFit( twoStageFit );
596 
599 
600  // Now ready for minimisation step
601  std::cout << "\nINFO in LauSimFitMaster::fitExpt : Start minimisation...\n";
602  std::pair<Int_t,Double_t> fitResult = LauFitter::fitter()->minimise();
603 
604  fitStatus_ = fitResult.first;
605  NLL_ = fitResult.second;
606 
607  // If we're doing a two stage fit we can now release (i.e. float)
608  // the 2nd stage parameters and re-fit
609  if ( twoStageFit ) {
610 
611  if ( fitStatus_ != 3 ) {
612  std::cerr << "ERROR in LauSimFitMaster:fitExpt : Not running second stage fit since first stage failed." << std::endl;
614  } else {
619  fitResult = LauFitter::fitter()->minimise();
620  }
621  }
622 
623  fitStatus_ = fitResult.first;
624  NLL_ = fitResult.second;
625  const TMatrixD& covMat = LauFitter::fitter()->covarianceMatrix();
626  covMatrix_.Clear();
627  covMatrix_.ResizeTo( covMat.GetNrows(), covMat.GetNcols() );
628  covMatrix_.SetMatrixArray( covMat.GetMatrixArray() );
629 
630  // Store the final fit results and errors into protected internal vectors that
631  // all sub-classes can use within their own finalFitResults implementation
632  // used below (e.g. putting them into an ntuple in a root file)
635 }
636 
637 void LauSimFitMaster::setParsFromMinuit(Double_t* par, Int_t npar)
638 {
639  // This function sets the internal parameters based on the values
640  // that Minuit is using when trying to minimise the total likelihood function.
641 
642  // MINOS reports different numbers of free parameters depending on the
643  // situation, so disable this check
644  if ( ! withinAsymErrorCalc_ ) {
645  if (static_cast<UInt_t>(npar) != nFreeParams_) {
646  std::cerr << "ERROR in LauSimFitMaster::setParsFromMinuit : Unexpected number of free parameters: " << npar << ".\n";
647  std::cerr << " Expected: " << nFreeParams_ << ".\n" << std::endl;
648  gSystem->Exit(EXIT_FAILURE);
649  }
650  }
651 
652  // Despite npar being the number of free parameters
653  // the par array actually contains all the parameters,
654  // free and floating...
655  // Update all the parameters with their new values.
656  // Change the value in the array to be sent out to the slaves and the
657  // parameters themselves (so that constraints are correctly calculated)
658  for (UInt_t i(0); i<nParams_; ++i) {
659  if (!params_[i]->fixed()) {
660  parValues_[i] = par[i];
661  params_[i]->value(par[i]);
662  }
663  }
664 }
665 
667 {
668  if ( socketMonitor_ == 0 ) {
669  std::cerr << "ERROR in LauSimFitMaster::getTotNegLogLikelihood : Sockets not initialised." << std::endl;
670  return 0.0;
671  }
672 
673  // Send current values of the parameters to the slaves.
674  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
675 
676  std::vector<UInt_t>& indices = slaveIndices_[iSlave];
677  UInt_t nPars = indices.size();
678  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
679  vectorPar_[iSlave][iPar] = parValues_[ indices[iPar] ];
680  }
681 
682  TMessage* message = messagesToSlaves_[iSlave];
683  message->Reset( kMESS_ANY );
684  message->WriteUInt( nPars );
685  message->WriteFastArray( vectorPar_[iSlave], nPars );
686 
687  sSlaves_[iSlave]->Send(*message);
688  }
689 
690  Double_t negLogLike(0.0);
691  TSocket *sActive(0);
692  UInt_t responsesReceived(0);
693  while ( responsesReceived != nSlaves_ ) {
694 
695  sActive = socketMonitor_->Select();
696  sActive->Recv(messageFromSlave_);
697 
698  messageFromSlave_->ReadDouble( vectorRes_[responsesReceived] );
699 
700  negLogLike += vectorRes_[responsesReceived];
701 
702  ++responsesReceived;
703  }
704 
705  // Calculate any penalty terms from Gaussian constrained variables
706  if ( ! conVars_.empty() ){
707  negLogLike += this->getLogLikelihoodPenalty();
708  }
709 
710  return negLogLike;
711 }
712 
714 {
715  Double_t penalty(0.0);
716 
717  for ( std::vector<LauAbsRValue*>::const_iterator iter = conVars_.begin(); iter != conVars_.end(); ++iter ) {
718  Double_t val = (*iter)->value();
719  Double_t mean = (*iter)->constraintMean();
720  Double_t width = (*iter)->constraintWidth();
721 
722  Double_t term = ( val - mean )*( val - mean );
723  penalty += term/( 2*width*width );
724  }
725 
726  return penalty;
727 }
728 
729 void LauSimFitMaster::addConstraint(const TString& formula, const std::vector<TString>& pars, const Double_t mean, const Double_t width)
730 {
731  StoreConstraints newCon;
732  newCon.formula_ = formula;
733  newCon.conPars_ = pars;
734  newCon.mean_ = mean;
735  newCon.width_ = width;
736  storeCon_.push_back(newCon);
737 }
738 
740 {
741  // Add penalties from the constraints to fit parameters
742 
743  // First, constraints on the fit parameters themselves
744  for ( std::vector<LauParameter*>::const_iterator iter = params_.begin(); iter != params_.end(); ++iter ) {
745  if ( (*iter)->gaussConstraint() ) {
746  conVars_.push_back( *iter );
747  std::cout << "INFO in LauSimFitMaster::addConParameters : Added Gaussian constraint to parameter "<< (*iter)->name() << std::endl;
748  }
749  }
750 
751  // Second, constraints on arbitrary combinations
752  for ( std::vector<StoreConstraints>::iterator iter = storeCon_.begin(); iter != storeCon_.end(); ++iter ) {
753  std::vector<TString> names = (*iter).conPars_;
754  std::vector<LauParameter*> params;
755  for ( std::vector<TString>::iterator iternames = names.begin(); iternames != names.end(); ++iternames ) {
756  for ( std::vector<LauParameter*>::const_iterator iterfit = params_.begin(); iterfit != params_.end(); ++iterfit ) {
757  if ( (*iternames) == (*iterfit)->name() ){
758  params.push_back(*iterfit);
759  }
760  }
761  }
762 
763  // If the parameters are not found, skip it
764  if ( params.size() != (*iter).conPars_.size() ) {
765  std::cerr << "WARNING in LauSimFitMaster::addConParameters: Could not find parameters to constrain in the formula... skipping" << std::endl;
766  continue;
767  }
768 
769  LauFormulaPar* formPar = new LauFormulaPar( (*iter).formula_, (*iter).formula_, params );
770  formPar->addGaussianConstraint( (*iter).mean_, (*iter).width_ );
771  conVars_.push_back(formPar);
772 
773  std::cout << "INFO in LauSimFitMaster::addConParameters : Added Gaussian constraint to formula\n";
774  std::cout << " : Formula: " << (*iter).formula_ << std::endl;
775  for ( std::vector<LauParameter*>::iterator iterparam = params.begin(); iterparam != params.end(); ++iterparam ) {
776  std::cout << " : Parameter: " << (*iterparam)->name() << std::endl;
777  }
778  }
779 
780 }
781 
783 {
784  if ( socketMonitor_ == 0 ) {
785  std::cerr << "ERROR in LauSimFitMaster::finalise : Sockets not initialised." << std::endl;
786  return kFALSE;
787  }
788 
789  // Prepare the covariance matrices
790  covMatrices_.resize( nSlaves_ );
791 
792  LauParamFixed pred;
793 
794  std::map<UInt_t,UInt_t> freeParIndices;
795 
796  UInt_t counter(0);
797  for ( UInt_t iPar(0); iPar < nParams_; ++iPar ) {
798  const LauParameter* par = params_[iPar];
799  if ( ! pred(par) ) {
800  freeParIndices.insert( std::make_pair(iPar,counter) );
801  ++counter;
802  }
803  }
804 
805  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
806  const UInt_t nPar = slaveIndices_[iSlave].size();
807 
808  std::vector<UInt_t> freeIndices;
809  freeIndices.reserve( nPar );
810 
811  for ( UInt_t iPar(0); iPar < nPar; ++iPar ) {
812  UInt_t index = slaveIndices_[iSlave][iPar];
813  std::map<UInt_t,UInt_t>::iterator freeIter = freeParIndices.find(index);
814  if ( freeIter == freeParIndices.end() ) {
815  continue;
816  }
817  UInt_t freeIndex = freeIter->second;
818  freeIndices.push_back( freeIndex );
819  }
820 
821  const UInt_t nFreePars = freeIndices.size();
822  TMatrixD& covMat = covMatrices_[iSlave];
823  covMat.ResizeTo( nFreePars, nFreePars );
824 
825  for ( UInt_t iPar(0); iPar < nFreePars; ++iPar ) {
826  for ( UInt_t jPar(0); jPar < nFreePars; ++jPar ) {
827  UInt_t i = freeIndices[iPar];
828  UInt_t j = freeIndices[jPar];
829  covMat( iPar, jPar ) = covMatrix_( i, j );
830  }
831  }
832  }
833 
834  // The array to hold the parameters
835  TObjArray array;
836 
837  // Send messages to all slaves containing the final parameters and fit status, NLL
838  // TODO - at present we lose the information on the correlations between the parameters that are unique to each slave
839  // - so should we store the full correlation matrix in an ntuple? along with all the parameters?
840  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
841 
842  array.Clear();
843 
844  std::vector<UInt_t>& indices = slaveIndices_[iSlave];
845  UInt_t nPars = indices.size();
846  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
847  array.Add( params_[ indices[iPar] ] );
848  }
849 
850  TMatrixD& covMat = covMatrices_[iSlave];
851 
852  TMessage* message = messagesToSlaves_[iSlave];
853  message->Reset( kMESS_OBJECT );
854  message->WriteInt( fitStatus_ );
855  message->WriteDouble( NLL_ );
856  message->WriteObject( &array );
857  message->WriteObject( &covMat );
858 
859  sSlaves_[iSlave]->Send(*message);
860  }
861 
862  TSocket *sActive(0);
863  UInt_t responsesReceived(0);
864  Bool_t allOK(kTRUE);
865  while ( responsesReceived != nSlaves_ ) {
866 
867  // Get the next queued response
868  sActive = socketMonitor_->Select();
869 
870  // Extract from the message the ID of the slave and the number of events read
871  sActive->Recv( messageFromSlave_ );
872  UInt_t iSlave(0);
873  Bool_t ok(kTRUE);
874  messageFromSlave_->ReadUInt( iSlave );
875  messageFromSlave_->ReadBool( ok );
876 
877  if ( ok ) {
878  TObjArray * objarray = dynamic_cast<TObjArray*>( messageFromSlave_->ReadObject( messageFromSlave_->GetClass() ) );
879  if ( ! objarray ) {
880  std::cerr << "ERROR in LauSimFitMaster::finalise : Error reading finalised parameters from slave" << std::endl;
881  allOK = kFALSE;
882  } else {
883  // We want to auto-delete the supplied parameters since we only copy their values in this case
884  objarray->SetOwner(kTRUE);
885 
886  const UInt_t nPars = objarray->GetEntries();
887  if ( nPars != slaveIndices_[iSlave].size() ) {
888  std::cerr << "ERROR in LauSimFitMaster::finalise : Unexpected number of finalised parameters received from slave" << std::endl;
889  allOK = kFALSE;
890  } else {
891  for ( UInt_t iPar(0); iPar < nPars; ++iPar ) {
892  LauParameter* parameter = dynamic_cast<LauParameter*>( (*objarray)[iPar] );
893  if ( ! parameter ) {
894  std::cerr << "ERROR in LauSimFitMaster::finalise : Error reading parameter from slave" << std::endl;
895  allOK = kFALSE;
896  continue;
897  }
898 
899  TString parname = parameter->name();
900 
901  std::map< TString, UInt_t >::iterator iter = parIndices_.find( parname );
902  if ( iter == parIndices_.end() ) {
903  std::cerr << "ERROR in LauSimFitMaster::finalise : Unexpected parameter name received from slave" << std::endl;
904  allOK = kFALSE;
905  continue;
906  }
907 
908  const UInt_t index = iter->second;
909  if ( slaveIndices_[iSlave][iPar] != index ) {
910  std::cerr << "ERROR in LauSimFitMaster::finalise : Unexpected parameter received from slave" << std::endl;
911  allOK = kFALSE;
912  continue;
913  }
914 
915  Double_t parvalue = parameter->value();
916  params_[index]->value( parvalue );
917  parValues_[index] = parvalue;
918  vectorPar_[iSlave][iPar] = parvalue;
919  }
920  }
921  delete objarray;
922  }
923  } else {
924  std::cerr << "ERROR in LauSimFitMaster::finalise : Slave " << iSlave << " reports an error performing finalisation" << std::endl;
925  allOK = kFALSE;
926  }
927 
928  ++responsesReceived;
929  }
930 
931  // Fill our ntuple as well
932  if ( fitNtuple_ != 0 ) {
933  for ( std::vector<LauParameter*>::iterator iter = params_.begin(); iter != params_.end(); ++iter ) {
934  if (!(*iter)->fixed()) {
935  (*iter)->updatePull();
936  }
937  }
938  std::vector<LauParameter> extraVars;
939  fitNtuple_->storeParsAndErrors( params_, extraVars );
942  }
943 
944  return allOK;
945 }
946 
948 {
949  if ( socketMonitor_ == 0 ) {
950  std::cerr << "ERROR in LauSimFitMaster::writeOutResults : Sockets not initialised." << std::endl;
951  return kFALSE;
952  }
953 
954  // Construct a message, requesting to write out the fit results
955  TString msgStr("Write Results");
956  TMessage message( kMESS_STRING );
957  message.WriteTString( msgStr );
958 
959  // Send the message to the slaves
960  for ( UInt_t iSlave(0); iSlave<nSlaves_; ++iSlave ) {
961  sSlaves_[iSlave]->Send(message);
962  }
963 
964  TSocket *sActive(0);
965  UInt_t responsesReceived(0);
966  Bool_t allOK(kTRUE);
967  while ( responsesReceived != nSlaves_ ) {
968 
969  // Get the next queued response
970  sActive = socketMonitor_->Select();
971 
972  // Extract from the message the ID of the slave and the number of events read
973  sActive->Recv( messageFromSlave_ );
974  UInt_t iSlave(0);
975  Bool_t ok(kTRUE);
976  messageFromSlave_->ReadUInt( iSlave );
977  messageFromSlave_->ReadBool( ok );
978 
979  if ( ! ok ) {
980  std::cerr << "ERROR in LauSimFitMaster::writeOutResults : Slave " << iSlave << " reports an error performing finalisation" << std::endl;
981  allOK = kFALSE;
982  }
983 
984  ++responsesReceived;
985  }
986 
987  // Write out our ntuple as well
988  if (fitNtuple_ != 0) {
990  }
991 
992  return allOK;
993 }
994 
TMonitor * socketMonitor_
Parallel setup monitor.
std::vector< std::vector< UInt_t > > slaveIndices_
Lists of indices for each slave.
std::vector< LauAbsRValue * > conVars_
Gaussian constraints.
Bool_t readData()
Instruct the slaves to read the input data for the given experiment.
File containing declaration of LauFormulaPar class.
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.
void addConParameters()
Add parameters to the list of Gaussian constrained parameters.
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.
ClassImp(LauAbsCoeffSet)
std::vector< LauParameter * > params_
Parameters.
const TString & name() const
The parameter name.
Class for defining combinations of fit parameter objects.
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.
Double_t mean_
The mean value of the Gaussian constraint to be applied.
std::vector< TSocket * > sSlaves_
Sockets for each of the slaves.
const UInt_t reqPort_
The requested port.
std::vector< StoreConstraints > storeCon_
Store the constraints for fit parameters until initialisation is complete.
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.
void addGaussianConstraint(Double_t newGaussMean, Double_t newGaussWidth)
Add a Gaussian constraint (or modify an existing one)
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 class.
void initSockets()
Initialise socket connections for the slaves.
Struct to store constraint information until the fit is run.
std::vector< TString > conPars_
The list of LauParameter names to be used in the LauFormulaPar.
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.
Double_t getLogLikelihoodPenalty()
Calculate the penalty terms to the log likelihood from Gaussian constraints.
void checkInitFitParams()
Instruct the slaves to update the initial fit parameter values, if required.
virtual void addConstraint(const TString &formula, const std::vector< TString > &pars, const Double_t mean, const Double_t width)
Store constraint information for fit parameters.
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:64
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:33
void fitExpt(Bool_t useAsymmErrors, Bool_t twoStageFit)
Perform the fit for the current experiment.
TString formula_
The formula to be used in the LauFormulaPar.
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.
Double_t width_
The width of the Gaussian constraint to be applied.
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.