MINT2
Classes | Macros | Functions
ampFit.cpp File Reference
#include "Mint/FitParameter.h"
#include "Mint/NamedParameter.h"
#include "Mint/Minimiser.h"
#include "Mint/Neg2LL.h"
#include "Mint/Neg2LLSum.h"
#include "Mint/DalitzEventList.h"
#include "Mint/CLHEPPhysicalConstants.h"
#include "Mint/PdfBase.h"
#include "Mint/DalitzPdfBase.h"
#include "Mint/DalitzPdfBaseFastInteg.h"
#include "Mint/DalitzSumPdf.h"
#include "Mint/FitAmplitude.h"
#include "Mint/FitAmpSum.h"
#include "Mint/FitAmpIncoherentSum.h"
#include "Mint/DalitzEvent.h"
#include "Mint/AmpRatios.h"
#include "Mint/IEventGenerator.h"
#include "Mint/DalitzBWBoxSet.h"
#include "Mint/DalitzBoxSet.h"
#include "Mint/SignalGenerator.h"
#include "Mint/FromFileGenerator.h"
#include "Mint/LASSO.h"
#include "TGraph.h"
#include "TFile.h"
#include "TCanvas.h"
#include "TRandom2.h"
#include "TRandom3.h"
#include <ctime>
#include <iostream>
#include "Mint/DalitzPdfNormChecker.h"
#include "Mint/DalitzPdfBaseFlexiFastInteg.h"
#include "Mint/Chi2Binning.h"

Go to the source code of this file.

Classes

class  AmpsPdf
 

Macros

#define MULTIAMPS
 

Functions

int ampFit ()
 
int main ()
 

Macro Definition Documentation

◆ MULTIAMPS

#define MULTIAMPS

Definition at line 54 of file ampFit.cpp.

Function Documentation

◆ ampFit()

int ampFit ( )

Definition at line 139 of file ampFit.cpp.

139  {
140  time_t startTime = time(0);
141  //MinuitParameterSet& fitMPS(*MinuitParameterSet::getDefaultSet());
142  MinuitParameterSet fitMPS;//(*MinuitParameterSet::getDefaultSet());
143  MinuitParameterSet dummy;
144 
145  cout << "pset pointer in default " << MinuitParameterSet::getDefaultSet() << endl;
146  cout << "pset pointer in ampsFit " << &fitMPS << endl;
147 
148  cout << "sizes " << MinuitParameterSet::getDefaultSet()->size()
149  << ", " << fitMPS.size() << endl;
150 
151 
152  TRandom3 ranLux;
153  NamedParameter<int> RandomSeed("RandomSeed", 0);
154  ranLux.SetSeed((int)RandomSeed);
155  gRandom = &ranLux;
156 
157  FitAmplitude::AutogenerateFitFile();
158 
159  NamedParameter<double> lambda("lambda", 1.);
160 
161  NamedParameter<int> MakeIntegrators("MakeIntegrators", 1);
162  bool makeIntegrators = (int) MakeIntegrators;
163  cout << "MakeIntegrator= " << MakeIntegrators << endl;
164 
165  NamedParameter<string> InputFileName("InputFileName", (std::string) "");
166  bool generateNew = (std::string) InputFileName == "";
167 
168  string InputFileName1 = ((string)InputFileName) + "_1.root";
169  string InputFileName2 = ((string)InputFileName) + "_2.root";
170  string InputFileName3 = ((string)InputFileName) + "_3.root";
171  string InputFileName4 = ((string)InputFileName) + "_4.root";
172  string InputFileName5 = ((string)InputFileName) + "_5.root";
173 
174  NamedParameter<string> IntegratorEventFile("IntegratorEventFile"
175  , (std::string) "SignalIntegrationEvents"
176  , (char*) 0);
177  string SgIntegratorEventFile1 = "Sg" + ((string)IntegratorEventFile) + "_1.root";
178  string SgIntegratorEventFile2 = "Sg" + ((string)IntegratorEventFile) + "_2.root";
179  string SgIntegratorEventFile3 = "Sg" + ((string)IntegratorEventFile) + "_3.root";
180  string SgIntegratorEventFile4 = "Sg" + ((string)IntegratorEventFile) + "_4.root";
181  string SgIntegratorEventFile5 = "Sg" + ((string)IntegratorEventFile) + "_5.root";
182 
183  string BgIntegratorEventFile1 = "Bg" + ((string)IntegratorEventFile) + "_1.root";
184  string BgIntegratorEventFile2 = "Bg" + ((string)IntegratorEventFile) + "_2.root";
185  string BgIntegratorEventFile3 = "Bg" + ((string)IntegratorEventFile) + "_3.root";
186  string BgIntegratorEventFile4 = "Bg" + ((string)IntegratorEventFile) + "_4.root";
187  string BgIntegratorEventFile5 = "Bg" + ((string)IntegratorEventFile) + "_5.root";
188 
189 
190 
191  NamedParameter<string> IntegratorInputFile("IntegratorInputFile"
192  , (std::string) "sgIntegrator"
193  , (char*) 0);
194  NamedParameter<int> Nevents("Nevents", 10000);
195  NamedParameter<int> doScan("doScan", 0);
196  NamedParameter<std::string> integMethod("IntegMethod", (std::string) "efficient");
197  NamedParameter<double> integPrecision("IntegPrecision", 1.e-4);
198 
199  NamedParameter<int> EventPattern("Event Pattern", 421, -321, 211, 211, -211);
200  DalitzEventPattern pat(EventPattern.getVector());
201 
202  NamedParameter<int> doNormCheck("doNormCheck", 0);
203  NamedParameter<int> saveEvents("SaveEvents", 1);
204  NamedParameter<int> doFinalStats("DoFinalStats", 1);
205  cout << " got event pattern: " << pat << endl;
206 
207 
208  DalitzEventList eventList1, eventList2, eventList3, eventList4, eventList5;
209 
210  cout << "1 pset pointer in default " << MinuitParameterSet::getDefaultSet() << endl;
211  cout << "1 pset pointer in ampsFit " << &fitMPS << endl;
212 
213  cout << "1 sizes " << MinuitParameterSet::getDefaultSet()->size()
214  << ", " << fitMPS.size() << endl;
215 
216  if(! generateNew){
217  cout << "reading events from file " << InputFileName1 << endl;
218  eventList1.fromFile(InputFileName1);
219  cout << " I've got " << eventList1.size() << " events (1)." << endl;
220 
221  cout << "reading events from file " << InputFileName2 << endl;
222  eventList2.fromFile(InputFileName2);
223  cout << " I've got " << eventList2.size() << " events (2)." << endl;
224 
225  cout << "reading events from file " << InputFileName3 << endl;
226  eventList3.fromFile(InputFileName3);
227  cout << " I've got " << eventList3.size() << " events (3)." << endl;
228 
229  cout << "reading events from file " << InputFileName4 << endl;
230  eventList4.fromFile(InputFileName4);
231  cout << " I've got " << eventList4.size() << " events (4)." << endl;
232 
233  cout << "reading events from file " << InputFileName5 << endl;
234  eventList5.fromFile(InputFileName5);
235  cout << " I've got " << eventList5.size() << " events (5)." << endl;
236  }
237 
238  if(generateNew){
239  SignalGenerator sg(pat);
240 
241  cout << "Generating " << Nevents << " signal events (1)." << endl;
242  sg.FillEventList(eventList1, Nevents);
243  if((int) saveEvents)eventList1.save("KKpipi_1.root");
244 
245  cout << "Generating " << Nevents << " signal events (2)." << endl;
246  sg.FillEventList(eventList2, Nevents);
247  if((int) saveEvents)eventList2.save("KKpipi_2.root");
248 
249  cout << "Generating " << Nevents << " signal events (3)." << endl;
250  sg.FillEventList(eventList3, Nevents);
251  if((int) saveEvents)eventList3.save("KKpipi_3.root");
252 
253  cout << "Generating " << Nevents << " signal events (4)." << endl;
254  sg.FillEventList(eventList4, Nevents);
255  if((int) saveEvents)eventList4.save("KKpipi_4.root");
256 
257  cout << "Generating " << Nevents << " signal events (5)." << endl;
258  sg.FillEventList(eventList5, Nevents);
259  if((int) saveEvents)eventList4.save("KKpipi_5.root");
260  }
261 
262  //DalitzHistoSet datH = eventList.histoSet();
263  //datH.save("plotsFromEventList.root");
264 
265 
266  cout << "2 pset pointer in default " << MinuitParameterSet::getDefaultSet() << endl;
267  cout << "2 pset pointer in ampsFit " << &fitMPS << endl;
268 
269  cout << "2 sizes " << MinuitParameterSet::getDefaultSet()->size()
270  << ", " << fitMPS.size() << endl;
271 
272  /*
273  AmpsPdf amps(pat, &fitMPS
274  , integPrecision
275  , integMethod
276  , (std::string) IntegratorEventFile);
277  */
278 
279  FitAmpSum ampSum1(pat, &fitMPS);
280  FitAmpIncoherentSum ampSumBg1(pat, &dummy, "Inco1_");
281 
282 #ifdef MULTIAMPS
283  /*
284  FitAmpSum ampSum2(pat, &fitMPS);
285  FitAmpSum ampSum3(pat, &fitMPS);
286  FitAmpSum ampSum4(pat, &fitMPS);
287  FitAmpSum ampSum5(pat, &fitMPS);
288 
289  FitAmpIncoherentSum ampSumBg2(pat, &dummy, "Inco2_");
290  FitAmpIncoherentSum ampSumBg3(pat, &dummy, "Inco3_");
291  FitAmpIncoherentSum ampSumBg4(pat, &dummy, "Inco4_");
292  FitAmpIncoherentSum ampSumBg5(pat, &dummy, "Inco5_");
293  */
294 #endif
295 
296  counted_ptr<FitAmpSum> ampSum2Ptr = ampSum1.GetCloneSameFitParameters();
297  counted_ptr<FitAmpSum> ampSum3Ptr = ampSum1.GetCloneSameFitParameters();
298  counted_ptr<FitAmpSum> ampSum4Ptr = ampSum1.GetCloneSameFitParameters();
299  counted_ptr<FitAmpSum> ampSum5Ptr = ampSum1.GetCloneSameFitParameters();
300 
301  counted_ptr<FitAmpSum> ampSumBg2Ptr = ampSumBg1.GetCloneSameFitParameters();
302  counted_ptr<FitAmpSum> ampSumBg3Ptr = ampSumBg1.GetCloneSameFitParameters();
303  counted_ptr<FitAmpSum> ampSumBg4Ptr = ampSumBg1.GetCloneSameFitParameters();
304  counted_ptr<FitAmpSum> ampSumBg5Ptr = ampSumBg1.GetCloneSameFitParameters();
305 
306 
307  FitParameter sgFrac("SignalFraction");
308 
309  AmpsPdf ampsSg1(pat, &ampSum1
310  , integPrecision
311  , integMethod
312  , SgIntegratorEventFile1, true, &fitMPS);
313 
314 
315 #ifdef MULTIAMPS
316  // Note: Crash also occurs if we use the same ampSum in each AmpsPdf
317  AmpsPdf ampsSg2(pat, ampSum2Ptr.get() // <<< putting ampSum1 here also crashes.
318  , integPrecision
319  , integMethod
320  , SgIntegratorEventFile2, true, &fitMPS);
321 
322  AmpsPdf ampsSg3(pat, ampSum3Ptr.get()
323  , integPrecision
324  , integMethod
325  , SgIntegratorEventFile3, true, &fitMPS);
326 
327  AmpsPdf ampsSg4(pat, ampSum4Ptr.get()
328  , integPrecision
329  , integMethod
330  , SgIntegratorEventFile4, true, &fitMPS);
331 
332  AmpsPdf ampsSg5(pat, ampSum5Ptr.get()
333  , integPrecision
334  , integMethod
335  , SgIntegratorEventFile5, true, &fitMPS);
336 #endif
337  AmpsPdf ampsBg1(pat, &ampSumBg1
338  , integPrecision
339  , integMethod
340  , BgIntegratorEventFile1, true, &fitMPS);
341 
342 
343 #ifdef MULTIAMPS
344  // Note: Crash also occurs if we use the same ampSum in each AmpsPdf
345  AmpsPdf ampsBg2(pat, ampSumBg2Ptr.get()
346  , integPrecision
347  , integMethod
348  , BgIntegratorEventFile2, true, &fitMPS);
349 
350  AmpsPdf ampsBg3(pat, ampSumBg3Ptr.get()
351  , integPrecision
352  , integMethod
353  , BgIntegratorEventFile3, true, &fitMPS);
354 
355  AmpsPdf ampsBg4(pat, ampSumBg4Ptr.get()
356  , integPrecision
357  , integMethod
358  , BgIntegratorEventFile4, true, &fitMPS);
359 
360  AmpsPdf ampsBg5(pat, ampSumBg5Ptr.get()
361  , integPrecision
362  , integMethod
363  , BgIntegratorEventFile5, true, &fitMPS);
364 #endif
365 
366 
367  DalitzSumPdf amps1(sgFrac, ampsSg1, ampsBg1);
368 #ifdef MULTIAMPS
369  DalitzSumPdf amps2(sgFrac, ampsSg2, ampsBg2);
370  DalitzSumPdf amps3(sgFrac, ampsSg3, ampsBg3);
371  DalitzSumPdf amps4(sgFrac, ampsSg4, ampsBg4);
372  DalitzSumPdf amps5(sgFrac, ampsSg5, ampsBg5);
373 #endif
374 
375  cout << "3 pset pointer in default " << MinuitParameterSet::getDefaultSet() << endl;
376  cout << "3 pset pointer in ampsFit " << &fitMPS << endl;
377 
378  cout << "3 sizes " << MinuitParameterSet::getDefaultSet()->size()
379  << ", " << fitMPS.size() << endl;
380 
381  if(! makeIntegrators){
382  ampsSg1.setIntegratorFileName("sgIntegrator1");
383 #ifdef MULTIAMPS
384  ampsSg2.setIntegratorFileName("sgIntegrator2");
385  ampsSg3.setIntegratorFileName("sgIntegrator3");
386  ampsSg4.setIntegratorFileName("sgIntegrator4");
387  ampsSg5.setIntegratorFileName("sgIntegrator5");
388 #endif
389  }
390  if(! makeIntegrators){
391  ampsBg1.setIntegratorFileName("bgIntegrator1");
392 #ifdef MULTIAMPS
393  ampsBg2.setIntegratorFileName("bgIntegrator2");
394  ampsBg3.setIntegratorFileName("bgIntegrator3");
395  ampsBg4.setIntegratorFileName("bgIntegrator4");
396  ampsBg5.setIntegratorFileName("bgIntegrator5");
397 #endif
398  }
399 
400 
401 #ifdef MULTIAMPS
402  /* this also crashes
403  Neg2LL fcn1(amps1, eventList1);//, &fitMPS);
404  Neg2LL fcn2(amps1, eventList2);//, &fitMPS);
405  Neg2LL fcn3(amps1, eventList3);//, &fitMPS);
406  Neg2LL fcn4(amps1, eventList4);//, &fitMPS);
407  */
408 
409  Neg2LL fcn1(amps1, eventList1, &fitMPS);
410  Neg2LL fcn2(amps2, eventList2, &fitMPS);
411  Neg2LL fcn3(amps3, eventList3, &fitMPS);
412  Neg2LL fcn4(amps4, eventList4, &fitMPS);
413  Neg2LL fcn5(amps5, eventList5, &fitMPS);
414 
415 
416  LASSO lasso(&ampsSg3, lambda);
417 
418  //Neg2LLSum fcn(&fcn1, &fcn2, &fcn3, &fcn4, &fcn5, &lasso, &fitMPS);
419  //Neg2LLSum fcn(&fcn1, &fcn2, &fcn3, &fcn4, &lasso, &fitMPS);
420  Neg2LLSum fcn(&fitMPS, &fcn1, &fcn2, &fcn3, &fcn4, &lasso);
421 #else
422  Neg2LL fcn(amps1, eventList1);
423 #endif
424 
425  if((int) doNormCheck){
426  DalitzPdfNormChecker nc(&amps1, pat);
427  nc.checkNorm();
428  }
429 
430 
431 
432  Minimiser mini(&fcn);
433 
434  mini.doFit();
435 
436  mini.printResultVsInput();
437  ampsSg1.saveEachAmpsHistograms("singleAmpHistos");
438 
439  if(makeIntegrators){
440  cout << "now saving all integrators" <<endl;
441  ampsSg1.saveIntegrator("sgIntegrator1");
442 #ifdef MULTIAMPS
443  ampsSg2.saveIntegrator("sgIntegrator2");
444  ampsSg3.saveIntegrator("sgIntegrator3");
445  ampsSg4.saveIntegrator("sgIntegrator4");
446  ampsSg5.saveIntegrator("sgIntegrator5");
447 #endif
448  ampsBg1.saveIntegrator("sgIntegrator1");
449 #ifdef MULTIAMPS
450  ampsBg2.saveIntegrator("bgIntegrator2");
451  ampsBg3.saveIntegrator("bgIntegrator3");
452  ampsBg4.saveIntegrator("bgIntegrator4");
453  ampsBg5.saveIntegrator("bgIntegrator5");
454 #endif
455  }
456 
457  cout << "now calling doFinalStat" << endl;
458 
459  if((bool) doFinalStats){
460  ampsSg1.doFinalStats(&mini);
461 #ifdef MULTIAMPS
462  ampsSg2.doFinalStats(&mini);
463  ampsSg3.doFinalStats(&mini);
464  ampsSg4.doFinalStats(&mini);
465  ampsSg5.doFinalStats(&mini);
466 #endif
467  }
468 
469  cout<<"final fcn1 getVal() "<< fcn1.getVal()<<endl;
470 #ifdef MULTIAMPS
471  cout<<"final fcn2 getVal() "<< fcn2.getVal()<<endl;
472  cout<<"final fcn3 getVal() "<< fcn3.getVal()<<endl;
473  cout<<"final fcn4 getVal() "<< fcn4.getVal()<<endl;
474  cout<<"final fcn5 getVal() "<< fcn5.getVal()<<endl;
475 #endif
476  //return 0;
477 
478  cout << " ampFit done. Took " << (time(0) - startTime)/60.
479  << " min. Returning 0." << endl;
480  return 0;
481 }
DalitzEventList::fromFile
bool fromFile(const std::string &fname="DalitzEvents.root")
Definition: MyDalitzEventList.h:320
FitAmplitude::AutogenerateFitFile
static void AutogenerateFitFile(const std::string &fname="protoFitAmplitudeFile.txt", const DalitzEventPattern &pat=DalitzEventPattern::NoPattern)
Definition: FitAmplitude.cpp:50
DalitzEventList::save
bool save(const std::string &fname="DalitzEvents.root") const
Definition: MyDalitzEventList.h:317
MINT::MinuitParameterSet::size
unsigned int size() const
Definition: MinuitParameterSet.cpp:96
MINT::counted_ptr< FitAmpSum >
MINT::EventList::size
virtual unsigned int size() const
Definition: EventList.h:59
lambda
double lambda(double x, double y, double z)
Definition: lambda.h:8
MINT::counted_ptr::get
X * get() const
Definition: counted_ptr.h:123

◆ main()

int main ( )

Definition at line 484 of file ampFit.cpp.

484  {
485 
486  return ampFit();
487 
488 }
ampFit
int ampFit()
Definition: ampFit.cpp:139
Generated by   doxygen 1.8.15