#include <DalitzBox.h>

Public Member Functions
void	getReady ()

void	setFlatPhaseSpace (bool flat=true)

double	amps (IDalitzEvent &ep)

double	phaseSpace (IDalitzEvent &ep)

double	ampsWithPhaseSpace (IDalitzEvent &ep)

	DalitzBox ()

	DalitzBox (const DalitzEventPattern &pat, MINT::IReturnRealForEvent< IDalitzEvent > amps=0, TRandom rnd=gRandom)

	DalitzBox (const DalitzEventPattern &pat, const DalitzCoordinate &limit, MINT::IReturnRealForEvent< IDalitzEvent > amps=0, TRandom rnd=gRandom)

	DalitzBox (const DalitzEventPattern &pat, const std::vector< DalitzCoordinate > &limits, MINT::IReturnRealForEvent< IDalitzEvent > amps=0, TRandom rnd=gRandom)

	DalitzBox (const DalitzBox &other)

	~DalitzBox ()

void	setNumber (int num)

int	number () const

void	setName (const std::string &name)

const std::string &	name () const

void	setGuessedHeight (double h)

double	guessedHeight () const

const MappedDalitzArea &	area () const

MappedDalitzArea &	area ()

double	height () const

double	setHeight (double h)

DalitzEventList &	eventList ()

const DalitzEventList &	eventList () const

double	volume () const

void	setDaddy (const DalitzBox *daddy)

bool	setAmps (MINT::IReturnRealForEvent< IDalitzEvent > *amps)

bool	insideArea (const DalitzEvent &evt) const

bool	insideDaddysArea (const DalitzEvent &evt) const

bool	insideMyOrDaddysArea (const DalitzEvent &evt) const

MINT::counted_ptr< DalitzEvent >	tryEventForOwner ()

MINT::counted_ptr< DalitzEvent >	tryWeightedEventForOwner ()

void	encloseInPhaseSpaceArea ()

std::vector< DalitzBox >	split (unsigned int nWays) const

std::vector< DalitzBox >	splitIfWiderThan (double maxWidth) const

bool	setRnd (TRandom *rnd=gRandom)

void	print (std::ostream &os=std::cout) const

Private Member Functions
int	makeFlatEventsKeepAll (int N)

bool	estimateHeight (std::vector< double > &vals)

bool	estimateHeight_old (std::vector< double > &vals)

bool	estimateHeightMinuit ()

bool	makeStarterSet ()

int	throwAwayData (const std::vector< double > &vals)

MINT::counted_ptr< DalitzEvent >	tryNewEvent ()

MINT::counted_ptr< DalitzEvent >	tryEventFromList ()

MINT::counted_ptr< DalitzEvent >	popEventFromList ()

double	getEventsPdf (DalitzEvent &evt)

Private Attributes
double	_max_s0

double	_max_s1

double	_max_s2

double	_max_s3

double	_max_s4

double	_maxPhaseSpace

int	_number

std::string	_name

bool	_ready

MappedDalitzArea	_area

DalitzEventPattern	_pat

MINT::IReturnRealForEvent< IDalitzEvent > *	_amps

TRandom *	_rnd

double	_height

DalitzEventList	_eventList

int	_nTries

int	_nEventsForTries

const DalitzBox *	_daddy

double	_guessedHeight

unsigned int	_heightProblems

bool	_flatPhaseSpace

Detailed Description

Definition at line 22 of file DalitzBox.h.

Constructor & Destructor Documentation

◆ DalitzBox() [1/5]

DalitzBox::DalitzBox ( )

Definition at line 76 of file DalitzBox.cpp.

   :  _max_s0(-1), _max_s1(-1), _max_s2(-1), _max_s3(-1), _max_s4(-1)
   , _maxPhaseSpace(-9999)
   , _number(-1)
   , _name("noName")
   , _ready(false)
   , _area()
   , _pat()
   , _amps(0)
   , _rnd(0)
   , _height(-9999)
   , _eventList()
   , _nTries(0)
   , _nEventsForTries(0)
   , _daddy(0)
   , _guessedHeight(-9999)
   , _heightProblems(0)
   , _flatPhaseSpace(0)
 {}

◆ DalitzBox() [2/5]

DalitzBox::DalitzBox	(	const DalitzEventPattern &	pat,
		MINT::IReturnRealForEvent< IDalitzEvent > *	amps = `0`,
		TRandom *	rnd = `gRandom`
	)

Definition at line 97 of file DalitzBox.cpp.

   :  _max_s0(-1), _max_s1(-1), _max_s2(-1), _max_s3(-1), _max_s4(-1)
   , _maxPhaseSpace(-9999)
   , _number(-1)
   , _name("noName")
   , _ready(false)
   , _area(pat, rnd)
   , _pat(pat)
   , _amps(0)
   , _rnd(rnd)
   , _height(-9999)
   , _eventList()
   , _nTries(0)
   , _nEventsForTries(0)
   , _daddy(0)
   , _guessedHeight(-9999)
   , _heightProblems(0)
   , _flatPhaseSpace(0)
 {
   setAmps(amps);
 }

◆ DalitzBox() [3/5]

DalitzBox::DalitzBox	(	const DalitzEventPattern &	pat,
		const DalitzCoordinate &	limit,
		MINT::IReturnRealForEvent< IDalitzEvent > *	amps = `0`,
		TRandom *	rnd = `gRandom`
	)

Definition at line 121 of file DalitzBox.cpp.

   :  _max_s0(-1), _max_s1(-1), _max_s2(-1), _max_s3(-1), _max_s4(-1)
   , _maxPhaseSpace(-9999)
   , _number(-1)
   , _name("noName")
   , _ready(false)
   , _area(pat, limit, rnd)
   , _pat(pat)
   , _amps(0)
   , _rnd(rnd)
   , _height(-9999)
   , _eventList()
   , _nTries(0)
   , _nEventsForTries(0)
   , _daddy(0)
   , _guessedHeight(-9999)
   , _heightProblems(0)
   , _flatPhaseSpace(0)
 {
   setAmps(amps);
 }

◆ DalitzBox() [4/5]

DalitzBox::DalitzBox	(	const DalitzEventPattern &	pat,
		const std::vector< DalitzCoordinate > &	limits,
		MINT::IReturnRealForEvent< IDalitzEvent > *	amps = `0`,
		TRandom *	rnd = `gRandom`
	)

Definition at line 146 of file DalitzBox.cpp.

   :  _max_s0(-1), _max_s1(-1), _max_s2(-1), _max_s3(-1), _max_s4(-1)
   , _maxPhaseSpace(-9999)
   , _number(-1)
   , _name("noName")
   , _ready(false)
   , _area(pat, limits, rnd)
   , _pat(pat)
   , _amps(0)
   , _rnd(rnd)
   , _height(-9999)
   , _eventList()
   , _nTries(0)
   , _nEventsForTries(0)
   , _daddy(0)
   , _guessedHeight(-9999)
   , _heightProblems(0)
   , _flatPhaseSpace(0)
 {
   setAmps(amps);
 }

◆ DalitzBox() [5/5]

DalitzBox::DalitzBox ( const DalitzBox & other )

Definition at line 171 of file DalitzBox.cpp.

   :  _max_s0(other._max_s0)
   , _max_s1(other._max_s1)
   , _max_s2(other._max_s2)
   , _max_s3(other._max_s3)
   , _max_s4(other._max_s4)
   , _maxPhaseSpace(other._maxPhaseSpace)
   , _number(other._number)
   ,  _name(other._name)
   , _ready(other._ready)
   , _area(other._area)
   , _pat(other._pat)
   , _amps(other._amps)
   , _rnd(other._rnd)
   , _height(other._height)
   , _eventList(other._eventList) // expensive!
   , _nTries(other._nTries)
   , _nEventsForTries(other._nEventsForTries)
   , _daddy(other._daddy)
   , _guessedHeight(other._guessedHeight)
   , _heightProblems(other._heightProblems)
   , _flatPhaseSpace(other._flatPhaseSpace)
 {
   //  _ready = false;
   //  cout << "copying box" << endl;
 }

◆ ~DalitzBox()

DalitzBox::~DalitzBox ( )

Definition at line 206 of file DalitzBox.cpp.

206 {

207 }

Member Function Documentation

◆ amps()

double DalitzBox::amps ( IDalitzEvent & ep )

Definition at line 268 of file DalitzBox.cpp.

                                        {
   bool dbThis=false;
 
   if(0 == _amps) return 1;
   if(dbThis) cout << " the event: " << endl;
   if(dbThis) evt.print();
   if(dbThis) cout << " ... now getting RealVal() " << endl;
   double val = _amps->RealVal(evt);
   if(dbThis) cout << "  all done, returning " << val << endl;
   return val;
 }

◆ ampsWithPhaseSpace()

double DalitzBox::ampsWithPhaseSpace ( IDalitzEvent & ep )

Definition at line 285 of file DalitzBox.cpp.

                                                         {
   bool dbThis=false;
 
   double a = amps(evt_in);
   if(dbThis) cout << "amps " << a << endl;
   double ps = phaseSpace(evt_in);
   if(dbThis) cout << "ps " << ps << endl;
 
   return a*ps;
 }

◆ area() [1/2]

const MappedDalitzArea& DalitzBox::area ( ) const

inline

Definition at line 103 of file DalitzBox.h.

103 {return _area;}

DalitzBox::_area

MappedDalitzArea _area

Definition: DalitzBox.h:31

◆ area() [2/2]

MappedDalitzArea& DalitzBox::area ( )

inline

Definition at line 104 of file DalitzBox.h.

104 {return _area;}

DalitzBox::_area

MappedDalitzArea _area

Definition: DalitzBox.h:31

◆ encloseInPhaseSpaceArea()

void DalitzBox::encloseInPhaseSpaceArea ( )

inline

Definition at line 128 of file DalitzBox.h.

                                 {
     area().encloseInPhaseSpaceArea();
   }

◆ estimateHeight()

bool DalitzBox::estimateHeight ( std::vector< double > & vals )

private

Definition at line 406 of file DalitzBox.cpp.

                     {
   bool dbThis=false;
   int maxEvents =   100000;
   int maxTries  = 10000000;
   double meanPhaseSpaceSum = 0;
   if(fastTest){
     maxEvents =   10000;
     maxTries  = 1000000;
   }
   if(veryFastTest || ultraFast){
     maxEvents =    1000;
     maxTries  = 1000000;
   }
 
   if(guessedHeight() > 0){
     _height = guessedHeight();
     return 0;
   }
 
   //  double saveFactor =  2.0;
 
   cout << "DalitzBox::estimateHeight" << endl;
   //time_t startTime = time(0);
   
   double max = -9999;
 
   int numEvents=0, counter=0;
   while(numEvents < maxEvents && counter < maxTries){
     counter++;
     double tmp_s0 = _rnd->Rndm();
     double tmp_s1 = _rnd->Rndm();
     double tmp_s2 = _rnd->Rndm();
     double tmp_s3 = _rnd->Rndm();
     double tmp_s4 = _rnd->Rndm();
     counted_ptr<DalitzEvent> tmpEvtPtr(area().makeEventForOwner(tmp_s0, tmp_s1, tmp_s2, tmp_s3, tmp_s4));
     if(0 == tmpEvtPtr) continue;
     if(dbThis) cout << "estimate height got event "
             << " printing it: " << endl;
     if(dbThis) cout << *tmpEvtPtr << endl;
 
     //cout << " found a non-zero event!!!" << endl;
     if(this->insideDaddysArea(*tmpEvtPtr)) continue;
 
     numEvents++;
 
 
     DalitzEvent tmpEvt(*tmpEvtPtr);
     double val = ampsWithPhaseSpace(tmpEvt);
     vals.push_back(val);
     if(val > max){
       max=val;
       _max_s0 = tmp_s0; 
       _max_s1 = tmp_s1; 
       _max_s2 = tmp_s2;
       _max_s3 = tmp_s3;
       _max_s4 = tmp_s4;
     }
     double p =  phaseSpace(tmpEvt);
 
     meanPhaseSpaceSum += p;
     if(p > _maxPhaseSpace){
       _maxPhaseSpace = p;
     }
   };
 
   double meanPhaseSpace = meanPhaseSpaceSum/((double)numEvents);
   cout << "mean phase space " << meanPhaseSpace
        << ", max: " << _maxPhaseSpace << endl;
 
   //double epsilon = 0.25;
   
   double CL = 1.0 - 1.e-7;
   std::cout << "actual max = " << max << std::endl;
 
   if(numEvents < 50){
     cout << " no Paret of only " << _eventList.size()
      << " events. Taking 5*max = 5*" << max 
      << " = " << 5*max
      << " and return false." << endl;
     _height = 5.0 * max;
     return false;
   }
 
   double actualMax=-1, paretoMax=-1;
 
   max = generalisedPareto_estimateMaximum(vals
                       , CL
                       , actualMax
                       , paretoMax
                       );
   
   std::cout << " Pareto max " << paretoMax << std::endl;
   /*
   max *= 1.0 + epsilon;
   std::cout << "Now added " << epsilon * 100 << "% for safety "
         << max << std::endl;
 
   max *= saveFactor;
 
   std::cout << "and multiplied by safety factor " 
         << saveFactor
         << " giving " << max 
         << endl;
   */
   cout << "compare to guessed height: " 
        << guessedHeight()*_maxPhaseSpace << endl;
   _height = max;
 
   if(fabs( (paretoMax - actualMax)/actualMax) > 0.3 ){
     return false;
   }
 
   return true;
 }

◆ estimateHeight_old()

bool DalitzBox::estimateHeight_old ( std::vector< double > & vals )

private

Definition at line 522 of file DalitzBox.cpp.

                     {
   double saveFactor = 2.0;
 
   cout << "DalitzBox::estimateHeight" << endl;
   if(false && bigDebug){
     vals.resize(_eventList.size(), 1);
     _height=1;
     return true;
   }
   if(_eventList.empty()){
     _height = 0;
     return false;
   }
   time_t startTime = time(0);
   
   unsigned int counter=0;
 
   vals.clear();
 
   vals.resize(_eventList.size());
 
   double max = -9999;
 
   double maxPhaseSpace = -9999;
 
   for(unsigned int i=0; i < _eventList.size(); i++){
     IDalitzEvent& evt = _eventList[i];
      double p = evt.phaseSpace();
     if(p > maxPhaseSpace || 0 == counter){
       maxPhaseSpace = p;
     }
     double d=ampsWithPhaseSpace(evt);
 
     if(d > max || 0 == counter){
       max=d;
       //      _maxEvent = *(_eventList.getEvent());
     }
     vals[i] = d;
     if(i < 5 || 0 == i%(_eventList.size()/20 + 1)){
       std::cout << " calculated amps for event " 
         << i 
         << ", its value is " << d
         << std::endl;
       double deltaT = difftime(time(0), startTime);
       std::cout << " this took " << deltaT << " s";
       if(deltaT > 0.5){
     std::cout << "; rate (before throwing away) = " 
           << i/deltaT
           << " evts/s";
       }
       std::cout << std::endl;
     }
   };
 
   cout << " max phase space: " << maxPhaseSpace << endl;
 
   double epsilon = 0.25;
   
   double CL = 1.0 - 1./_eventList.size();
   std::cout << "actual max = " << max << std::endl;
 
   if(_eventList.size() < 50){
     cout << " no Paret of only " << _eventList.size()
      << " events. Taking 5*max = 5*" << max 
      << " = " << 5*max
      << " and return false." << endl;
     _height = 5.0 * max;
     return false;
   }
 
   double actualMax=-1, paretoMax=-1;
 
   max = generalisedPareto_estimateMaximum(vals
                       , CL
                       , actualMax
                       , paretoMax
                       );
   
   std::cout << " Pareto max " << paretoMax << std::endl;
   max *= 1.0 + epsilon;
   std::cout << "Now added " << epsilon * 100 << "% for safety "
         << max << std::endl;
 
   max *= saveFactor;
 
   std::cout << "and multiplied by safety factor " 
         << saveFactor
         << " giving " << max 
         << endl;
   cout << "compare to guessed height: " 
        << guessedHeight()*maxPhaseSpace << endl;
   _height = max;
 
   if(fabs( (paretoMax - actualMax)/actualMax) > 0.3 ){
     return false;
   }
 
   return true;
 }

◆ estimateHeightMinuit()

bool DalitzBox::estimateHeightMinuit ( )

private

Definition at line 296 of file DalitzBox.cpp.

                                     {
   MinuitParameterSet parset;
 
   /*
   // step 1: Find starting values:
   int minEvents=   10000000;
   int maxEvents=   10000001;
   int counter = 0;
   double max = -1;
   double s0=0.5, s1=0.5, s2=0.5, s3=0.5, s4=0.5;
   bool foundOne=false;
   cout << " starting loop" << endl;
   while((! (foundOne && counter > minEvents)) && counter < maxEvents){
     counter++;
     //    cout << " counter " << counter << endl;
     double tmp_s0 = _rnd->Rndm();
     double tmp_s1 = _rnd->Rndm();
     double tmp_s2 = _rnd->Rndm();
     double tmp_s3 = _rnd->Rndm();
     double tmp_s4 = _rnd->Rndm();
     counted_ptr<DalitzEvent> tmpEvtPtr(area().makeEventForOwner(tmp_s0, tmp_s1, tmp_s2, tmp_s3, tmp_s4));
     if(0 == tmpEvtPtr) continue;
     //cout << " found a non-zero event!!!" << endl;
     if(this->insideDaddysArea(*tmpEvtPtr)) return 0;
 
     foundOne=true;
     _amps->setEvent(&(*tmpEvtPtr));
     double val = _amps->RealVal();
     _amps->resetEventRecord();
     if(val > max){
       max=val;
       s0 = tmp_s0; s1=tmp_s1, s2=tmp_s2, s3=tmp_s3, s4=tmp_s4;
     }
   }
 
   if(max < 0){
     cout << " didn't find starting point "
      << " giving up; "
      << endl;
     _height = -1;
     return true;
   }
 
   cout << " max = " << max << " at "
        << s0 << ", " << s1 << ", " << s2 << ", " << s3 << ", " << s4 << endl;
 
   */
   FitParameter fit_scale0("scale0"
               , 0
               , _max_s0
               , 0.2
               , -0.05
               , 1.05
               , &parset
               );
   FitParameter fit_scale1("scale1"
               , 0
               , _max_s1
               , 0.2
               , -0.05
               , 1.05
               , &parset
               );
   FitParameter fit_scale2("scale2"
               , 0
               , _max_s2
               , 0.2
               , -0.05
               , 1.05
               , &parset
               );
   FitParameter fit_scale3("scale3"
               , 0
               , _max_s3
               , 0.2
               , -0.05
               , 1.05
               , &parset
               );
   FitParameter fit_scale4("scale4"
               , 0
               , _max_s4
               , 0.2
               , -0.05
               , 1.05
               , &parset
               );
 
   FindMaxFCN fcn(fit_scale0
          , fit_scale1
          , fit_scale2
          , fit_scale3
          , fit_scale4
          , this, &parset);
   
   Minimiser mini(&fcn);
   //  mini.Command("Set Strategy 2");
   mini.prepFit();
   mini.Command("Simplex 10000000");
   mini.Command("MIGRAD");
 
   double safetyMargin = 0;
 
   _height = -(mini.getFCNVal() * (1.0 + safetyMargin));
 
   if(_height <= 0.0) _height = -1.0;
   cout << " setting height to " << _height << endl;
   return true;
 }

◆ eventList() [1/2]

DalitzEventList& DalitzBox::eventList ( )

inline

Definition at line 112 of file DalitzBox.h.

112 { return _eventList;}

DalitzBox::_eventList

DalitzEventList _eventList

Definition: DalitzBox.h:38

◆ eventList() [2/2]

const DalitzEventList& DalitzBox::eventList ( ) const

inline

Definition at line 113 of file DalitzBox.h.

113 { return _eventList;}

DalitzBox::_eventList

DalitzEventList _eventList

Definition: DalitzBox.h:38

◆ getEventsPdf()

double DalitzBox::getEventsPdf ( DalitzEvent & evt )

private

Definition at line 830 of file DalitzBox.cpp.

                                               {
 
   if(0 == _amps || bigDebug) return evt.phaseSpace();
 
   double val = _amps->RealVal(evt) * evt.phaseSpace();
 
   return val;
 }

◆ getReady()

void DalitzBox::getReady ( )

Definition at line 796 of file DalitzBox.cpp.

                         {
   time_t startTime = time(0);
 
   cout << " DalitzBox::getReady() " << *this << endl;
   
   _nTries = 0;
   _nEventsForTries=0;
   _eventList.clear();
 
   makeStarterSet();
   _nEventsForTries = _eventList.size();
   
   double deltaT = difftime(time(0), startTime);
   cout << " DalitzBox::getRead(): am ready. This took " 
        << deltaT/60 << " min" << endl;
 
   if(no_reuse) _eventList.clear();
 
   _ready = true;
 }

◆ guessedHeight()

double DalitzBox::guessedHeight ( ) const

inline

Definition at line 101 of file DalitzBox.h.

101 {return _guessedHeight;}

DalitzBox::_guessedHeight

double _guessedHeight

Definition: DalitzBox.h:44

◆ height()

double DalitzBox::height ( ) const

inline

Definition at line 106 of file DalitzBox.h.

                        {
     return _height;
   }

◆ insideArea()

bool DalitzBox::insideArea ( const DalitzEvent & evt ) const

Definition at line 214 of file DalitzBox.cpp.

                                                       {
   /*
   cout << "Hello from DalitzBox::insideArea()" << endl;
   cout << " my area is " << _area << endl;
   cout << ", or the same with (): " << area() << endl;
   cout << " and the event: " << evt << endl;
   cout << " Let's check if the event is inside the area: " << endl;
   */
   bool inside = area().isInside(evt);
 
   /*
   cout << "Goodbye from DalitzBox::insideArea(); returning " 
        << inside << endl;
   */
 
   return inside;
 }

◆ insideDaddysArea()

bool DalitzBox::insideDaddysArea ( const DalitzEvent & evt ) const

Definition at line 232 of file DalitzBox.cpp.

                                                             {
   if(0 == _daddy) return false;
   return _daddy->insideMyOrDaddysArea(evt);
 }

◆ insideMyOrDaddysArea()

bool DalitzBox::insideMyOrDaddysArea ( const DalitzEvent & evt ) const

Definition at line 237 of file DalitzBox.cpp.

                                                                 {
   return insideArea(evt) || insideDaddysArea(evt);
 }

◆ makeFlatEventsKeepAll()

int DalitzBox::makeFlatEventsKeepAll ( int N )

private

Definition at line 241 of file DalitzBox.cpp.

                                          {
   cout << "making " << N << " flat event" << endl;
   int counter = 0;
   int guessedWasteFactor = 110*5;
   if(fastTest) guessedWasteFactor = 150;
 
   for(int i=0; i< N*guessedWasteFactor; i++){
     _nTries ++;
     bool printdbg = (_nTries < 5 || 0 == _nTries%100000);
     if(printdbg) cout << " nTries " << _nTries << endl;
       
     counted_ptr<DalitzEvent> evt(area().makeEventForOwner());
     if(printdbg) cout << " made event" << endl;
 
     if(0 != evt && (! insideDaddysArea(*evt))){ // this should work
       if(printdbg)cout << " area OK " << endl;
       _eventList.Add(*evt);
       if(printdbg)cout << " added event" << endl;
 
       if(++counter >= N) break;
     }
   }
 
   cout << " got " << _eventList.size() << endl;
   return _nTries;
 }

◆ makeStarterSet()

bool DalitzBox::makeStarterSet ( )

private

Definition at line 623 of file DalitzBox.cpp.

                               {
   double safetyFactor = 2.5;
 
   time_t startTime = time(0);
   cout << "making starter set" << endl;
   // unsigned int maxEvents = 400000;
   int NperLoop           =  20000;
   //  int maxTries           = maxEvents * 110;
   int maxLoop            =  5;
   if(fastTest){
     NperLoop = 10000;
     maxLoop  = 3;
   }
   if(veryFastTest){
     NperLoop = 1000;
     maxLoop  = 1;
   }
 
   bool successfulHeightEstimate = false;
   std::vector<double> vals;
 
   int counter = 0;
   double prevHeight = -0.1;
   double heightErrorEstimate = 0.1;
 
   bool doOldHeight=(guessedHeight() < 10);
   bool doMinuitHeight=false;
   if(doOldHeight){
     do{    
       //  makeFlatEventsKeepAll(NperLoop);
       successfulHeightEstimate = estimateHeight(vals);
       heightErrorEstimate 
     = fabs(prevHeight - _height)/( 0.5*(fabs(prevHeight) + fabs(_height) ));
 
       /*
       successfulHeightEstimate |= heightErrorEstimate < 0.25;
       // the above line means: I believe it if pareto and actual
       // are close - or if at least pareto is consistent.
       */
 
       successfulHeightEstimate &= heightErrorEstimate < 0.25;
       // the above line means: I believe it if pareto and actual
       // are close AND if Pareto is consistent.
       counter++;
       double deltaT = difftime(time(0), startTime);
       cout << " making " << counter 
        << " " << NperLoop << "-events starter set(s)"
        << " took " << deltaT/60 << " min" << endl;
       cout << " current heightErrorEstimate: " << heightErrorEstimate << endl;
       prevHeight = _height;
       // no event inside box - box must have 0 volume
       // (most likely overlap 100% with other volume)
     }while(counter < maxLoop
        && (! successfulHeightEstimate)
        );
     
     cout << " DONE \"old\" height guess, now MINUIT " << endl;
   }
   double search_height = _height;
   double minuit_height = -1;
   if(doMinuitHeight){
     successfulHeightEstimate |= estimateHeightMinuit();
     minuit_height = _height;
   }
 
 
   _height = ( search_height > minuit_height ? search_height : minuit_height);
   _height = ( guessedHeight() > _height     ? guessedHeight() : _height);
 
   cout << "Different Heights:"
        << "\n\t  search_height " << search_height
        << "\n\t  minuit_height " << minuit_height
        << "\n\t guessed height " << guessedHeight()
        << "\n\t  chosen height " << _height
        << endl;
 
   if(! successfulHeightEstimate){
     cout << "DalitzBox::makeStarterSet "
      << "WARNING - no reliable estimate of box-height"
      << endl;
     cout << " I'll multiply it by " << 1.0 + 2*heightErrorEstimate
      << ", and continue..." << endl;
     //    _height *= 1.0 + 2.0*heightErrorEstimate;
   }else{
     cout << "DalitzBox::makeStarterSet "
      << "I think I found a reliable Height Estimate of "
      << _height
      << endl;
     /*    _height *= 1.0 + 2.0*heightErrorEstimate;
       cout << "For safety, I'll add 2* the Error estimate of "
       << heightErrorEstimate << " which gives a final height of "
       << _height
       << endl;
     */
   }
   cout << " appliying safety factor of " << safetyFactor << endl;
   _height *= safetyFactor;
   cout << "final height " << _height << endl;
 
   //throwAwayData(vals);
 
   return true;
 }

◆ name()

const std::string& DalitzBox::name ( ) const

inline

Definition at line 98 of file DalitzBox.h.

98 {return _name;}

DalitzBox::_name

std::string _name

Definition: DalitzBox.h:27

◆ number()

int DalitzBox::number ( ) const

inline

Definition at line 95 of file DalitzBox.h.

95 {return _number;}

DalitzBox::_number

int _number

Definition: DalitzBox.h:26

◆ phaseSpace()

double DalitzBox::phaseSpace ( IDalitzEvent & ep )

Definition at line 280 of file DalitzBox.cpp.

                                                 {
   if(_flatPhaseSpace) return 1;
   return evt_in.phaseSpace();
 }

◆ popEventFromList()

counted_ptr< DalitzEvent > DalitzBox::popEventFromList ( )

private

Definition at line 908 of file DalitzBox.cpp.

                                                     {
   if(_eventList.empty()){
     return counted_ptr<DalitzEvent>(0);
   }
   counted_ptr<DalitzEvent> evt(new DalitzEvent(_eventList[_eventList.size()-1]));
   _eventList.resize(_eventList.size()-1);
   return evt;
 }

◆ print()

void DalitzBox::print ( std::ostream & os = std::cout ) const

Definition at line 961 of file DalitzBox.cpp.

                                          {
   os << "DalitzBox: " << name() << " number " << number()
      << "\n area " << area()
      << "\n guessed height " << guessedHeight()
      << "\n height " << height()
      << "\n daddy: ";
   if(0 == _daddy){
     os << " none";
   }else{
     os << " number " << _daddy->number();
   }
 }

◆ setAmps()

bool DalitzBox::setAmps ( MINT::IReturnRealForEvent< IDalitzEvent > * amps )

Definition at line 198 of file DalitzBox.cpp.

                                                               {
   _ready = false;
   _amps = amps;
   _eventList.clear();
   if(0 == _amps) return false;
   return true;
 }

◆ setDaddy()

void DalitzBox::setDaddy ( const DalitzBox * daddy )

Definition at line 209 of file DalitzBox.cpp.

                                               {
   _ready = false;
   _daddy = daddy;
 }

◆ setFlatPhaseSpace()

void DalitzBox::setFlatPhaseSpace ( bool flat = true )

inline

Definition at line 68 of file DalitzBox.h.

                                         {
     _flatPhaseSpace = flat;
   }

◆ setGuessedHeight()

void DalitzBox::setGuessedHeight ( double h )

inline

Definition at line 100 of file DalitzBox.h.

100 {_guessedHeight = h;}

DalitzBox::_guessedHeight

double _guessedHeight

Definition: DalitzBox.h:44

◆ setHeight()

double DalitzBox::setHeight ( double h )

inline

Definition at line 109 of file DalitzBox.h.

                             {
     return _height = h;
   }

◆ setName()

void DalitzBox::setName ( const std::string & name )

inline

Definition at line 97 of file DalitzBox.h.

97 { _name = name;}

DalitzBox::name

const std::string & name() const

Definition: DalitzBox.h:98

DalitzBox::_name

std::string _name

Definition: DalitzBox.h:27

◆ setNumber()

void DalitzBox::setNumber ( int num )

inline

Definition at line 94 of file DalitzBox.h.

94 {_number = num;}

DalitzBox::_number

int _number

Definition: DalitzBox.h:26

◆ setRnd()

bool DalitzBox::setRnd ( TRandom * rnd = gRandom )

Definition at line 974 of file DalitzBox.cpp.

                                   {
   _rnd = rnd;
   area().setRnd(rnd);
   return true;
 }

◆ split()

std::vector< DalitzBox > DalitzBox::split ( unsigned int nWays ) const

Definition at line 917 of file DalitzBox.cpp.

                            {
   std::vector<MappedDalitzArea> areaList(_area.split(nWays));
 
   std::vector<DalitzBox> newSet;
   if(areaList.empty()) return newSet;
   bool foundCtr=false;
   for(unsigned int i=0; i < areaList.size(); i++){
     DalitzBox newBox(*this);
     newBox.area() = areaList[i];
     if( ! (newBox.area().isInside((this->area().centre())))){
       //      newBox.setGuessedHeight(0.6*guessedHeight());
       newBox.setGuessedHeight(-1);
     }else{
       foundCtr=true;
     }
     newSet.push_back(newBox);
   }
   if(false && (! foundCtr)){
     cout << "ERROR split box " << *this << " " << nWays << " ways,"
      << " but center: is gone." << endl;
     cout << "newSet\n: " << newSet << endl;
   }
   return newSet;
 }

◆ splitIfWiderThan()

std::vector< DalitzBox > DalitzBox::splitIfWiderThan ( double maxWidth ) const

Definition at line 943 of file DalitzBox.cpp.

                                   {
   std::vector<MappedDalitzArea> areaList(_area.splitIfWiderThan(maxWidth));
 
   std::vector<DalitzBox> newSet;
   if(areaList.empty()) return newSet;
   for(unsigned int i=0; i < areaList.size(); i++){
     DalitzBox newBox(*this);
     newBox.area() = areaList[i];
     if( ! (newBox.area().isInside((this->area().centre())))){
       //      newBox.setGuessedHeight(0.6*guessedHeight());
       newBox.setGuessedHeight(-1);
     }
     newSet.push_back(newBox);
   }
   return newSet;
 }

◆ throwAwayData()

int DalitzBox::throwAwayData ( const std::vector< double > & vals )

private

Definition at line 727 of file DalitzBox.cpp.

                                                          {
   if(false && bigDebug) return 0;
   time_t startTime = time(0);
   
   if(vals.size() != _eventList.size()){
     cout << " ERROR in  DalitzBox::throwAwayData: "
      << " vals.size()=" << vals.size()
      << ", but _eventList.size()=" 
      << _eventList.size()
      << endl;
     throw "this is bad.";
   }
   
   DalitzEventList newList;
 
   int rememberSize = _eventList.size();
 
   unsigned int counter=0;
 
   for(unsigned int i=0; i < _eventList.size(); i++){
     IDalitzEvent& evt(_eventList[i]);
 
     double d=ampsWithPhaseSpace(evt);
     //    double d = vals[counter]; 
 
     if(_rnd->Rndm()*_height < d){
       newList.Add( _eventList[i] );
     }
     if(counter < 10 || 0 == counter%(_eventList.size()/20 + 1)){
       std::cout << " remembering amps for event " 
         << counter 
         << ", its value is " << d
         << std::endl;
     }
     counter++;
   }
 
   _eventList = newList;
   std::cout << "now my size has changed to " << _eventList.size()
         << std::endl;
   std::cout << " So the waste factor for \n" << this->name() << " is ";
   if(_eventList.size() > 0) std::cout << rememberSize/_eventList.size();
   else std::cout << " infinity! - that's big!";
   std::cout << std::endl;
   
   
   double deltaTFinal = difftime(time(0), startTime);
   std::cout << " this took " << deltaTFinal/60.0 << " min";
   if(deltaTFinal > 0){
     std::cout << " rate = " << (_eventList.size()/deltaTFinal) << " evts/s"
           << " or " << (_eventList.size()/deltaTFinal) *60 << " evts/m"
           << " or " << (_eventList.size()/deltaTFinal) *60.*60./1000. 
           << "k evts/h";
   }
   std::cout << std::endl;
   
   std::cout << " ---------------\n " << std::endl;
 
   return _eventList.size();
 }

◆ tryEventForOwner()

counted_ptr< DalitzEvent > DalitzBox::tryEventForOwner ( )

Definition at line 817 of file DalitzBox.cpp.

                                                     {
   if(! _ready) getReady();
   
   counted_ptr<DalitzEvent> evt(0);
 
   if(_eventList.empty()){
     evt = tryNewEvent();
   }else{
     evt = tryEventFromList();
   }
   return evt;
 }

◆ tryEventFromList()

counted_ptr< DalitzEvent > DalitzBox::tryEventFromList ( )

private

Definition at line 894 of file DalitzBox.cpp.

                                                     {
   // the list has already been subjected to
   // data throwing away.
   // _nTries were made, and _nEventsForTries are
   // left over. So each event in the list represents
   // on average _nTries/_nEventsForTries tries;
 
   double rndNumber = _rnd->Rndm()*_nTries;
   if(rndNumber > (double) _nEventsForTries){
     return counted_ptr<DalitzEvent>(0);
   }
   return popEventFromList();  
 }

◆ tryNewEvent()

counted_ptr< DalitzEvent > DalitzBox::tryNewEvent ( )

private

Definition at line 839 of file DalitzBox.cpp.

                                                {
 
   counted_ptr<DalitzEvent> evt( area().makeEventForOwner() );
 
   if( 0 == evt){                         // bad ptr
     return counted_ptr<DalitzEvent>(0);
   }
   if(insideDaddysArea(*evt)){            // out of area
     return counted_ptr<DalitzEvent>(0);
   }
 
   double val = getEventsPdf(*evt);
   if(val > _height){
     _heightProblems++;
     cout << "ERROR in DalitzBox::tryNewEvent()\n" << (*this)
      << "\n   val > _height: " << val << " > " << _height
      << "\n  (for phase space = " << evt->phaseSpace() << ")."
      << "\n  if _height had been correctly estimated,"
      << " this should not be possible." << endl;
     cout << " This has happened " << _heightProblems 
      << " times in this box " << endl;
   }
   if(_rnd->Rndm()*_height > val){ // then we don't take it.
     evt = counted_ptr<DalitzEvent>(0);
   }
   return evt;
 }

◆ tryWeightedEventForOwner()

counted_ptr< DalitzEvent > DalitzBox::tryWeightedEventForOwner ( )

Definition at line 867 of file DalitzBox.cpp.

                                                             {
   bool dbThis=false;
   counted_ptr<DalitzEvent> evt( area().makeEventForOwner() );
 
   if( 0 == evt){                         // bad ptr
     return counted_ptr<DalitzEvent>(0);
   }
   if(insideDaddysArea(*evt)){            // out of area
     return counted_ptr<DalitzEvent>(0);
   }
 
   if(dbThis)cout << " DalitzBox::tryWeightedEventForOwner() "
          << " got event with ptr "
          << evt
          << ", setting val: " << endl;
   if(dbThis)cout << " now printing event: "
          << *evt
          << endl;
   double val = ampsWithPhaseSpace(*evt);
   if(dbThis)cout << " val = " << val << ", now setting weight " << endl;
   evt->setWeight(val/_height);
   if(dbThis)cout << " done weight, now generator..." << endl;
   evt->setGeneratorPdfRelativeToPhaseSpace(amps(*evt));
   if(dbThis)cout << " all OK, returning " << evt << endl;
   return evt;
 }

◆ volume()

double DalitzBox::volume ( ) const

Definition at line 788 of file DalitzBox.cpp.

                               {
   double a = area().size();
   cout << " volume: returning area " << a
        << " times height " << _height
        << " = " << a*_height
        << endl;
   return a * _height;
 }

Member Data Documentation

◆ _amps

MINT::IReturnRealForEvent<IDalitzEvent>* DalitzBox::_amps

private

Definition at line 34 of file DalitzBox.h.

◆ _area

MappedDalitzArea DalitzBox::_area

private

Definition at line 31 of file DalitzBox.h.

◆ _daddy

const DalitzBox* DalitzBox::_daddy

private

Definition at line 42 of file DalitzBox.h.

◆ _eventList

DalitzEventList DalitzBox::_eventList

private

Definition at line 38 of file DalitzBox.h.

◆ _flatPhaseSpace

bool DalitzBox::_flatPhaseSpace

private

Definition at line 47 of file DalitzBox.h.

◆ _guessedHeight

double DalitzBox::_guessedHeight

private

Definition at line 44 of file DalitzBox.h.

◆ _height

double DalitzBox::_height

private

Definition at line 37 of file DalitzBox.h.

◆ _heightProblems

unsigned int DalitzBox::_heightProblems

private

Definition at line 45 of file DalitzBox.h.

◆ _max_s0

double DalitzBox::_max_s0

private

Definition at line 23 of file DalitzBox.h.

◆ _max_s1

double DalitzBox::_max_s1

private

Definition at line 23 of file DalitzBox.h.

◆ _max_s2

double DalitzBox::_max_s2

private

Definition at line 23 of file DalitzBox.h.

◆ _max_s3

double DalitzBox::_max_s3

private

Definition at line 23 of file DalitzBox.h.

◆ _max_s4

double DalitzBox::_max_s4

private

Definition at line 23 of file DalitzBox.h.

◆ _maxPhaseSpace

double DalitzBox::_maxPhaseSpace

private

Definition at line 24 of file DalitzBox.h.

◆ _name

std::string DalitzBox::_name

private

Definition at line 27 of file DalitzBox.h.

◆ _nEventsForTries

int DalitzBox::_nEventsForTries

private

Definition at line 41 of file DalitzBox.h.

◆ _nTries

int DalitzBox::_nTries

private

Definition at line 40 of file DalitzBox.h.

◆ _number

int DalitzBox::_number

private

Definition at line 26 of file DalitzBox.h.

◆ _pat

DalitzEventPattern DalitzBox::_pat

private

Definition at line 33 of file DalitzBox.h.

◆ _ready

bool DalitzBox::_ready

private

Definition at line 29 of file DalitzBox.h.

◆ _rnd

TRandom* DalitzBox::_rnd

private

Definition at line 36 of file DalitzBox.h.

The documentation for this class was generated from the following files:

Mint/DalitzBox.h
src/Mojito/BetterMC/DalitzBox.cpp

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Constructor & Destructor Documentation

◆ DalitzBox() [1/5]

◆ DalitzBox() [2/5]

◆ DalitzBox() [3/5]

◆ DalitzBox() [4/5]

◆ DalitzBox() [5/5]

◆ ~DalitzBox()

Member Function Documentation

◆ amps()

◆ ampsWithPhaseSpace()

◆ area() [1/2]

◆ area() [2/2]

◆ encloseInPhaseSpaceArea()

◆ estimateHeight()

◆ estimateHeight_old()

◆ estimateHeightMinuit()

◆ eventList() [1/2]

◆ eventList() [2/2]

◆ getEventsPdf()

◆ getReady()

◆ guessedHeight()

◆ height()

◆ insideArea()

◆ insideDaddysArea()

◆ insideMyOrDaddysArea()

◆ makeFlatEventsKeepAll()

◆ makeStarterSet()

◆ name()

◆ number()

◆ phaseSpace()

◆ popEventFromList()

◆ print()

◆ setAmps()

◆ setDaddy()

◆ setFlatPhaseSpace()

◆ setGuessedHeight()

◆ setHeight()

◆ setName()

◆ setNumber()

◆ setRnd()

◆ split()

◆ splitIfWiderThan()

◆ throwAwayData()

◆ tryEventForOwner()

◆ tryEventFromList()

◆ tryNewEvent()

◆ tryWeightedEventForOwner()

◆ volume()

Member Data Documentation

◆ _amps

◆ _area

◆ _daddy

◆ _eventList

◆ _flatPhaseSpace

◆ _guessedHeight

◆ _height

◆ _heightProblems

◆ _max_s0

◆ _max_s1

◆ _max_s2

◆ _max_s3

◆ _max_s4

◆ _maxPhaseSpace

◆ _name

◆ _nEventsForTries

◆ _nTries

◆ _number

◆ _pat

◆ _ready

◆ _rnd