MINT2
Public Member Functions | Public Attributes | Private Member Functions | Private Attributes | List of all members
DalitzBWArea Class Reference

#include <DalitzBWArea.h>

Public Member Functions

bool checkIntegration () const
 
bool unWeightPs () const
 
void setUnWeightPs (bool doSo=true)
 
DalitzCoordinate & t01 ()
 
DalitzCoordinate & s12 ()
 
DalitzCoordinate & s23 ()
 
DalitzCoordinate & s34 ()
 
DalitzCoordinate & t40 ()
 
const DalitzCoordinate & t01 () const
 
const DalitzCoordinate & s12 () const
 
const DalitzCoordinate & s23 () const
 
const DalitzCoordinate & s34 () const
 
const DalitzCoordinate & t40 () const
 
MINT::counted_ptr< IGenFct > f_t01 ()
 
MINT::counted_ptr< IGenFct > f_s12 ()
 
MINT::counted_ptr< IGenFct > f_s23 ()
 
MINT::counted_ptr< IGenFct > f_s34 ()
 
MINT::counted_ptr< IGenFct > f_t40 ()
 
 DalitzBWArea ()
 
 DalitzBWArea (const DalitzEventPattern &pat, TRandom *rnd=gRandom)
 
 DalitzBWArea (const DalitzBWArea &other)
 
 ~DalitzBWArea ()
 
DalitzBWArea & operator= (const DalitzBWArea &other)
 
void setFcn (const DalitzCoordinate &c, const MINT::counted_ptr< IGenFct > &fcn)
 
double genValue (const IDalitzEvent *evtPtr) const
 
double genValue (const IDalitzEvent *evtPtr, const Permutation &mapping) const
 
void setPattern (const DalitzEventPattern &pat)
 
bool isInside (const DalitzEvent &evt) const
 
bool isInside (const DalitzEvent &evt, const Permutation &mapping) const
 
bool isInside (const DalitzCoordinate &dc) const
 
bool isInside (const std::vector< DalitzCoordinate > &dcList) const
 
double size () const
 
double integral () const
 
double integral3 () const
 
double integral4 () const
 
MINT::counted_ptr< DalitzEvent > tryEventForOwner (const Permutation &mapping=Permutation::unity()) const
 
bool setRnd (TRandom *rnd=gRandom)
 
void print (std::ostream &os=std::cout) const
 

Public Attributes

std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
 
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _mappedCoords
 

Private Member Functions

void makeCoords ()
 
void makeCoord (int i, int j)
 
void makeCoord (int i, int j, int k)
 
void makeMiMa ()
 
int ResonanceConfigurationNumber () const
 
MINT::counted_ptr< DalitzEvent > try3Event (const Permutation &mapping=Permutation::unity()) const
 
MINT::counted_ptr< DalitzEvent > try4Event (const Permutation &mapping=Permutation::unity()) const
 
MINT::counted_ptr< DalitzEvent > make4EventWithPhaseSpace (const Permutation &mapping=Permutation::unity()) const
 
MINT::counted_ptr< DalitzEvent > make4EventWithPhaseSpace (double &nTries, const Permutation &mapping=Permutation::unity()) const
 
MINT::counted_ptr< DalitzEvent > try4EventWithPhaseSpace (double &maxWeight, const Permutation &mapping=Permutation::unity()) const
 
MINT::counted_ptr< DalitzEvent > try4EventWithPhaseSpace (const Permutation &mapping=Permutation::unity()) const
 
MINT::counted_ptr< DalitzEvent > try3EventWithPhaseSpace (double &maxWeight, const Permutation &mapping=Permutation::unity()) const
 
MINT::counted_ptr< DalitzEvent > try3EventWithPhaseSpace (const Permutation &mapping=Permutation::unity()) const
 
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf (const DalitzCoordinate &c)
 
const std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf (const DalitzCoordinate &c) const
 
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf (int i, int j, int k)
 
const std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf (int i, int j, int k) const
 
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf (int i, int j)
 
const std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf (int i, int j) const
 
double MC4Integral (double &prec) const
 
double MC4IntegralNoTransform (double &prec) const
 
MINT::counted_ptr< DalitzEvent > tryFlat4EventWithPhaseSpace (double &maxWeight, const Permutation &mapping=Permutation::unity()) const
 
MINT::counted_ptr< DalitzEvent > tryFlat4EventWithPhaseSpace (const Permutation &mapping=Permutation::unity()) const
 
double genValueRho (const IDalitzEvent *evtPtr) const
 
std::vector< TLorentzVector > & mapP4 (const std::vector< TLorentzVector > &v_in, const Permutation &mapping, std::vector< TLorentzVector > &v_out) const
 
std::vector< TLorentzVector > mapP4 (const std::vector< TLorentzVector > &v_in, const Permutation &mapping) const
 

Private Attributes

DalitzEventPattern _pat
 
TRandom * _rnd
 
bool _madeCMap
 
double _areaIntegral
 
bool _unWeightPs
 

Detailed Description

Definition at line 23 of file DalitzBWArea.h.

Constructor & Destructor Documentation

◆ DalitzBWArea() [1/3]

DalitzBWArea::DalitzBWArea ( )

Definition at line 24 of file DalitzBWArea.cpp.

25  : _pat()
26  , _rnd(gRandom)
27  , _madeCMap(false)
28  , _areaIntegral(-9999.0)
29  , _unWeightPs(false)
30  , _coords()
31  , _mappedCoords()
32 {
33  makeCoords();
34 }
DalitzBWArea::_mappedCoords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _mappedCoords
Definition: DalitzBWArea.h:121
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzBWArea::_rnd
TRandom * _rnd
Definition: DalitzBWArea.h:27
DalitzBWArea::_unWeightPs
bool _unWeightPs
Definition: DalitzBWArea.h:34
DalitzBWArea::_areaIntegral
double _areaIntegral
Definition: DalitzBWArea.h:33
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ DalitzBWArea() [2/3]

DalitzBWArea::DalitzBWArea ( const DalitzEventPattern &  pat,
TRandom *  rnd = gRandom 
)

Definition at line 58 of file DalitzBWArea.cpp.

59  : _pat(pat)
60  , _rnd(rnd)
61  , _madeCMap(false)
62  , _areaIntegral(-9999.0)
63  , _unWeightPs(false)
64  , _coords()
65  , _mappedCoords()
66 {
67  makeCoords();
68 
69  if(_pat.numDaughters() != 4
70  && _pat.numDaughters() != 3
71  ){
72  cout << "DalitzBWArea::DalitzBWArea:"
73  << "\n Sorry, can only deal with 3 and 4 body decays. "
74  << "\n Please improve me so I can deal with decays like this: "
75  << pat << endl;
76  cout << " I'll crash now." << endl;
77  throw "DalitzBWArea: \"I'm no superman\"";
78  }
79  makeMiMa();
80 }
DalitzBWArea::_mappedCoords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _mappedCoords
Definition: DalitzBWArea.h:121
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzBWArea::_rnd
TRandom * _rnd
Definition: DalitzBWArea.h:27
DalitzBWArea::_unWeightPs
bool _unWeightPs
Definition: DalitzBWArea.h:34
DalitzBWArea::_areaIntegral
double _areaIntegral
Definition: DalitzBWArea.h:33
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ DalitzBWArea() [3/3]

DalitzBWArea::DalitzBWArea ( const DalitzBWArea &  other)

Definition at line 36 of file DalitzBWArea.cpp.

37  : _pat(other._pat)
38  , _rnd(other._rnd)
39  , _madeCMap(false)
40  , _areaIntegral(other._areaIntegral)
41  , _unWeightPs(other._unWeightPs)
42  , _coords(other._coords)
43  , _mappedCoords(other._mappedCoords)
44 {
45 }
DalitzBWArea::_mappedCoords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _mappedCoords
Definition: DalitzBWArea.h:121
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzBWArea::_rnd
TRandom * _rnd
Definition: DalitzBWArea.h:27
DalitzBWArea::_unWeightPs
bool _unWeightPs
Definition: DalitzBWArea.h:34
DalitzBWArea::_areaIntegral
double _areaIntegral
Definition: DalitzBWArea.h:33
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ ~DalitzBWArea()

DalitzBWArea::~DalitzBWArea ( )

Definition at line 161 of file DalitzBWArea.cpp.

161  {
162 }

Member Function Documentation

◆ checkIntegration()

bool DalitzBWArea::checkIntegration ( ) const

Definition at line 611 of file DalitzBWArea.cpp.

611  {
612 
613  cout << " DalitzBWArea::checkIntegration() Checking integration."
614  << "\n ResonanceConfigurationNumber "
615  << ResonanceConfigurationNumber()
616  << "\n Analytical: ";
617  if(_pat.numDaughters() != 4){
618  cout << "DalitzBWArea::checkIntegration() only works for 4-body decays"
619  << " returning 0." << endl;
620  return 0;
621  }
622 
623  double analytical=0;
624  if(ResonanceConfigurationNumber() == 0){
625  PhaseSpaceIntegral4bodyWith_s12s34 psi(&(*sf(1,2).second)
626  , &(*sf(3,4).second), _pat);
627 
628  analytical = psi.getVal();
629  }else{
630  PhaseSpaceIntegral4bodyWith_s123s12 psi(&(*sf(1,2,3).second)
631  , &(*sf(1,2).second), _pat);
632  analytical = psi.getVal();
633  }
634  double prec_1=-9999, prec_2=-9999;
635 
636  double numerical_1 = MC4Integral(prec_1);
637  double numerical_2 = MC4IntegralNoTransform(prec_2);
638 
639  cout << "DalitzBWArea::checkIntegration(): Result:"
640  << "\n\t root's integral " << analytical
641  << "\n\t MC with trafo " << numerical_1
642  << " +/- " << prec_1 * numerical_1
643  << "\n\t MC w/o trafo " << numerical_2
644  << " +/- " << prec_2 * numerical_2
645  << endl;
646  return true;
647 
648 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165
DalitzBWArea::MC4IntegralNoTransform
double MC4IntegralNoTransform(double &prec) const
Definition: DalitzBWArea.cpp:555
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzBWArea::MC4Integral
double MC4Integral(double &prec) const
Definition: DalitzBWArea.cpp:505
second
static const double second
Definition: CLHEPSystemOfUnits.h:115
PhaseSpaceIntegral4bodyWith_s12s34
Definition: phaseSpaceIntegrals.h:66
PhaseSpaceIntegral4bodyWith_s123s12
Definition: phaseSpaceIntegrals.h:78
DalitzBWArea::ResonanceConfigurationNumber
int ResonanceConfigurationNumber() const
Definition: DalitzBWArea.cpp:736

◆ f_s12()

MINT::counted_ptr<IGenFct> DalitzBWArea::f_s12 ( )
inline

Definition at line 162 of file DalitzBWArea.h.

162  {
163  return sf(1,2).second;
164  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ f_s23()

MINT::counted_ptr<IGenFct> DalitzBWArea::f_s23 ( )
inline

Definition at line 165 of file DalitzBWArea.h.

165  {
166  return sf(2,3).second;
167  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ f_s34()

MINT::counted_ptr<IGenFct> DalitzBWArea::f_s34 ( )
inline

Definition at line 168 of file DalitzBWArea.h.

168  {
169  return sf(3,4).second;
170  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ f_t01()

MINT::counted_ptr<IGenFct> DalitzBWArea::f_t01 ( )
inline

Definition at line 159 of file DalitzBWArea.h.

159  {
160  return sf(2,3,4).second;
161  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ f_t40()

MINT::counted_ptr<IGenFct> DalitzBWArea::f_t40 ( )
inline

Definition at line 171 of file DalitzBWArea.h.

171  {
172  return sf(1,2,3).second;
173  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ genValue() [1/2]

double DalitzBWArea::genValue ( const IDalitzEvent *  evtPtr) const

Definition at line 297 of file DalitzBWArea.cpp.

297  {
298  //bool dbThis=false;
299  if(0 == evtPtr) return 0;
300  double p=1;
301 
302  for(map<DalitzCoordKey, pair<DalitzCoordinate, counted_ptr<IGenFct> > >::const_iterator
303  it = _coords.begin();
304  it != _coords.end(); it++){
305  //DalitzCoordinate c = it->second.first;
306  counted_ptr<IGenFct> fct(it->second.second);
307 
308  p *= fct->generatingFctValue( evtPtr->sij(it->second.first) );
309  }
310  return fabs(p);
311 }
IDalitzEvent::sij
virtual double sij(const MINT::PolymorphVector< int > &indices) const =0
MINT::counted_ptr< IGenFct >
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ genValue() [2/2]

double DalitzBWArea::genValue ( const IDalitzEvent *  evtPtr,
const Permutation &  mapping 
) const

Definition at line 312 of file DalitzBWArea.cpp.

313  {
314  bool dbThis=false;
315 
316  if(0 == evtPtr) return 0;
317 
318  /*
319  we save significant time by not doing this check (below) - so let's risk it.
320 
321  if(mapping.mapOrder(_pat) != evtPtr->eventPattern()){
322  cout << "ERROR in DalitzBWArea::genValue(const IDalitzEvent* evtPtr) const"
323  << " patterns don't match!"
324  << " compare: " << _pat << " and "
325  << evtPtr->eventPattern() << endl;
326  return 0;
327  }
328  */
329 
330 
331  double p=1;
332 
333  for(map<DalitzCoordKey, pair<DalitzCoordinate, counted_ptr<IGenFct> > >::const_iterator
334  it = _coords.begin();
335  it != _coords.end(); it++){
336  DalitzCoordinate c = it->second.first;
337  counted_ptr<IGenFct> fct(it->second.second);
338 
339 
340  if(dbThis){
341  cout << "DalitzBWArea::genValue: factor for variable "
342  << c
343  << " s with mapping = " << evtPtr->sij(mapping.mapValues(c))
344  << " s with new mapping = " << evtPtr->sij(c.mapMe(mapping))
345  << ", and without= " << evtPtr->sij(c)
346  << " fct->generatingPDFValue with mapping "
347  << fct->generatingFctValue( evtPtr->sij(mapping.mapValues(c)) )
348  << " new mapping "
349  << fct->generatingFctValue( evtPtr->sij(c.mapMe(mapping)) )
350  << ", and without: "
351  << fct->generatingFctValue( evtPtr->sij(c) )
352  << endl;
353  }
354 
355  p *= fct->generatingFctValue( evtPtr->sij(it->second.first.mapMe(mapping)));
356  if(dbThis) cout << " p = " << p << endl;
357  }
358  if(dbThis) cout << " returning p " << p << endl;
359 
360  return fabs(p);
361 }
IDalitzEvent::sij
virtual double sij(const MINT::PolymorphVector< int > &indices) const =0
MINT::counted_ptr< IGenFct >
DalitzCoordinate::mapMe
DalitzCoordinate mapMe(const Permutation &perm) const
Definition: DalitzCoordinate.cpp:81
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119
Permutation::mapValues
std::vector< int > & mapValues(const std::vector< int > &in, std::vector< int > &out) const
Definition: Permutation.h:70

◆ genValueRho()

double DalitzBWArea::genValueRho ( const IDalitzEvent *  evtPtr) const
private

Definition at line 363 of file DalitzBWArea.cpp.

363  {
364  // for debugging only.
365  bool dbThis=false;
366 
367  if(0 == evtPtr) return 0;
368  if(_pat != evtPtr->eventPattern()){
369  cout << "ERROR in DalitzBWArea::genValue(const IDalitzEvent* evtPtr) const"
370  << " patterns don't match!"
371  << " compare: " << _pat << " and "
372  << evtPtr->eventPattern() << endl;
373  return 0;
374  }
375  double p=1;
376 
377  for(map<DalitzCoordKey, pair<DalitzCoordinate, counted_ptr<IGenFct> > >::const_iterator
378  it = _coords.begin();
379  it != _coords.end(); it++){
380  DalitzCoordinate c = it->second.first;
381  counted_ptr<IGenFct> fct(it->second.second);
382  double rho = fct->coordTransformFromS(evtPtr->sij(c));
383  if(dbThis)cout << "DalitzBWArea::genValueRho: factor for variable "
384  << c
385  << " s = " << evtPtr->sij(c)
386  << " fct->generatingPDFValue "
387  << fct->transformedFctValue(rho )
388  << endl;
389  p *= fct->transformedFctValue(rho );
390  if(dbThis) cout << " p = " << p << endl;
391  }
392  if(dbThis) cout << " returning p " << p << endl;
393 
394  return fabs(p);
395 }
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
IDalitzEvent::sij
virtual double sij(const MINT::PolymorphVector< int > &indices) const =0
MINT::counted_ptr< IGenFct >
IDalitzEvent::eventPattern
virtual const DalitzEventPattern & eventPattern() const =0
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ integral()

double DalitzBWArea::integral ( ) const

Definition at line 397 of file DalitzBWArea.cpp.

397  {
398  //bool dbThis=false;
399  if(_areaIntegral < 0){
400 
401  if(_pat.numDaughters() == 3){
402  _areaIntegral = integral3();
403  }else if(_pat.numDaughters() == 4){
404  _areaIntegral = integral4();
405  }else{
406  cout << "ERROR in DalitzBWArea::integral() can only make events"
407  << " with 3 or 4 daughters. You want : " << _pat
408  << endl;
409  _areaIntegral = 0;
410  }
411  }
412  return _areaIntegral;
413 }
DalitzBWArea::integral4
double integral4() const
Definition: DalitzBWArea.cpp:650
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzBWArea::_areaIntegral
double _areaIntegral
Definition: DalitzBWArea.h:33
DalitzBWArea::integral3
double integral3() const
Definition: DalitzBWArea.cpp:721

◆ integral3()

double DalitzBWArea::integral3 ( ) const

Definition at line 721 of file DalitzBWArea.cpp.

721  {
722  if(unWeightPs()){
723  // See integral4() for details on the meaning of unWeightPs().
724  // Here the integration is only over one parameter, s12, otherwise
725  // it's the same idea.
726  cout << "Don't know how to calculate integral3WithPhaseSpace() "
727  << "with option \"unWeightPs\" - please implement this."
728  << "\n I will crash now."
729  << endl;
730  throw "no unweighted ps for 3-body, yet";
731  }else{
732  return sf(1,2).second->integral();
733  }
734 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165
DalitzBWArea::unWeightPs
bool unWeightPs() const
Definition: DalitzBWArea.h:124

◆ integral4()

double DalitzBWArea::integral4 ( ) const

Definition at line 650 of file DalitzBWArea.cpp.

650  {
651  bool dbThis=false;
652 
653  if(_pat.numDaughters() != 4){
654  cout << "DalitzBWArea::integral4() called, although this pattern: "
655  << _pat
656  << " is a " << _pat.numDaughters() << " body decay."
657  << " I'll crash now." << endl;
658  throw "wrong integral4";
659  }
660  if(dbThis){
661  if(sf(1,2).second->flat() && sf(3,4).second->flat()){
662  cout << "check this integral: " << endl;
663  PhaseSpaceIntegral4bodyWith_s12s34 psiA(&(*sf(1,2).second)
664  , &(*sf(3,4).second), _pat);
665  double vA = psiA.getVal();
666  cout << " 12, 34: " << vA << endl;
667  PhaseSpaceIntegral4bodyWith_s123s12 psiB(&(*sf(1,2,3).second)
668  , &(*sf(1,2).second), _pat);
669  double vB = psiB.getVal();
670  cout << " 123, 12: " << vB << endl;
671  PhaseSpaceIntegral4body psiC;
672  double vC = psiC.getVal(_pat);
673  cout << " tested " << vC << endl;
674  cout << " check integrals: "
675  << vA << " , " << vB << " , " << vC << endl;
676  }
677  }
678 
679  // below:
680  // unWeightPs() == true: We produce events of weight 1 according to
681  // PDF = BW * phase-space-density
682  // this means, to get the relative normalisations of the different
683  // Breit Wigners right (i.e. the probability that this particular
684  // instance of DalitzBWArea is called to produce an event), we
685  // need to calculate the integral WITH phase space, i.e.
686  // integral BreitWigner * phase-space-density dt01 ds12
687  // or
688  // integral BreitWigner * phase-space-density ds12 ds23
689  // depending on the resonance.
690  //
691  // unWeightPs() == false: we produce events of weight phase-space-density according to
692  // PDF = BW
693  // this means, to get the relative normalisations of the different
694  // Breit Wigners right (i.e. the probability that this particular
695  // instance of DalitzBWArea is called to produce an event), we
696  // need to calculate simply the integral WITHOUT phase space, i.e.
697  // integral BreitWigner * ds123 ds12
698  // or
699  // integral BreitWigner * ds12 ds34
700  // which is much easier.
701 
702  if(ResonanceConfigurationNumber() == 0){
703  if(unWeightPs()){
704  PhaseSpaceIntegral4bodyWith_s12s34 psi(&(*sf(1,2).second)
705  , &(*sf(3,4).second), _pat);
706  return psi.getVal();
707  }else{
708  return sf(1,2).second->integral() * sf(3,4).second->integral();
709  }
710  }else{
711  if(unWeightPs()){
712  PhaseSpaceIntegral4bodyWith_s123s12 psi(&(*sf(1,2,3).second)
713  , &(*sf(1,2).second), _pat);
714  return psi.getVal();
715  }else{
716  return sf(1,2,3).second->integral() * sf(1,2).second->integral();
717  }
718  }
719 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
second
static const double second
Definition: CLHEPSystemOfUnits.h:115
DalitzBWArea::unWeightPs
bool unWeightPs() const
Definition: DalitzBWArea.h:124
PhaseSpaceIntegral4bodyWith_s12s34
Definition: phaseSpaceIntegrals.h:66
PhaseSpaceIntegral4bodyWith_s123s12
Definition: phaseSpaceIntegrals.h:78
DalitzBWArea::ResonanceConfigurationNumber
int ResonanceConfigurationNumber() const
Definition: DalitzBWArea.cpp:736
PhaseSpaceIntegral4body::getVal
double getVal(const DalitzEventPattern &_pat)
Definition: phaseSpaceIntegrals.cpp:260

◆ isInside() [1/4]

bool DalitzBWArea::isInside ( const DalitzEvent &  evt) const

Definition at line 211 of file DalitzBWArea.cpp.

211  {
212  Permutation mapping(evt.eventPattern().size());
213  mapping.makeUnity();
214  return isInside(evt, mapping);
215 }
DalitzEvent::eventPattern
virtual const DalitzEventPattern & eventPattern() const
Definition: DalitzEvent.h:173
Permutation::makeUnity
void makeUnity()
Definition: Permutation.cpp:42
MINT::PolymorphVector::size
unsigned int size() const
Definition: PolymorphVector.h:46
DalitzBWArea::isInside
bool isInside(const DalitzEvent &evt) const
Definition: DalitzBWArea.cpp:211

◆ isInside() [2/4]

bool DalitzBWArea::isInside ( const DalitzEvent &  evt,
const Permutation &  mapping 
) const

Definition at line 217 of file DalitzBWArea.cpp.

218  {
219  bool dbThis=false;
220  //cout << "Hello from DalitzBWArea::isInside" << endl;
221  if(mapping.mapOrder(_pat) != (MINT::PolymorphVector<DecayTreeItem>) evt.eventPattern()){
222  if(dbThis) cout << "patterns don't match in isInside" << endl;
223  return false;
224  }
225  //cout << " Made Coordinate map " << endl;
226  std::vector<int> mappedValues;
227  for(map<DalitzCoordKey, pair<DalitzCoordinate, counted_ptr<IGenFct> > >::const_iterator
228  it = _coords.begin();
229  it != _coords.end(); it++){
230  if(dbThis) {
231  cout << "old: " << evt.sij(it->second.first)
232  << " new " << evt.sij(mapping.mapValues(it->second.first))
233  << " newer " << evt.sij(it->second.first.mapMe(mapping))
234  << endl;
235  }
236  // the following line is instead of the also possible, nicer looking,
237  // but malloc-intensive evt.sij(mapping.mapValues(it->second.first))
238  mapping.mapValues(it->second.first, mappedValues);
239  double val = evt.sij(mappedValues);
240  if(val < it->second.first.min() || val >= it->second.first.max()) return false;
241  }
242  //cout << "returning true" << endl;
243  return true;
244 }
DalitzEvent::sij
virtual double sij(const std::vector< int > &indices) const
Definition: DalitzEvent.cpp:709
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzEvent::eventPattern
virtual const DalitzEventPattern & eventPattern() const
Definition: DalitzEvent.h:173
MINT::counted_ptr< IGenFct >
second
static const double second
Definition: CLHEPSystemOfUnits.h:115
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119
Permutation::mapValues
std::vector< int > & mapValues(const std::vector< int > &in, std::vector< int > &out) const
Definition: Permutation.h:70
Permutation::mapOrder
std::vector< T > & mapOrder(const std::vector< T > &in, std::vector< T > &out) const
Definition: Permutation.h:91

◆ isInside() [3/4]

bool DalitzBWArea::isInside ( const DalitzCoordinate &  dc) const

Definition at line 245 of file DalitzBWArea.cpp.

245  {
246 
247  //cout << " Made Coordinate map " << endl;
248  double val = dc.val();
249  map<DalitzCoordKey, pair<DalitzCoordinate, counted_ptr<IGenFct> > >::const_iterator it = _coords.find(dc);
250  if(it == _coords.end()){
251  cout << "ERROR in DalitzBWArea::isInside - unknown coordinate "
252  << dc << endl;
253  return false;
254  }
255  if(it->second.first.min() <= val && it->second.first.max() > val){
256  return true;
257  }else{
258  return false;
259  }
260 }
DalitzCoordinate::val
double val() const
Definition: DalitzCoordinate.h:41
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ isInside() [4/4]

bool DalitzBWArea::isInside ( const std::vector< DalitzCoordinate > &  dcList) const

Definition at line 261 of file DalitzBWArea.cpp.

261  {
262  for(unsigned int i=0; i < dcList.size(); i++){
263  if(! isInside(dcList[i])){
264  return false;
265  }
266  }
267  return true;
268 }
DalitzBWArea::isInside
bool isInside(const DalitzEvent &evt) const
Definition: DalitzBWArea.cpp:211

◆ make4EventWithPhaseSpace() [1/2]

counted_ptr< DalitzEvent > DalitzBWArea::make4EventWithPhaseSpace ( const Permutation &  mapping = Permutation::unity()) const
private

Definition at line 480 of file DalitzBWArea.cpp.

480  {
481  double nt;
482  return make4EventWithPhaseSpace(nt, mapping);
483 }
DalitzBWArea::make4EventWithPhaseSpace
MINT::counted_ptr< DalitzEvent > make4EventWithPhaseSpace(const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:480

◆ make4EventWithPhaseSpace() [2/2]

counted_ptr< DalitzEvent > DalitzBWArea::make4EventWithPhaseSpace ( double &  nTries,
const Permutation &  mapping = Permutation::unity() 
) const
private

Definition at line 484 of file DalitzBWArea.cpp.

485  {
486  int nearInfinity = 100000000;
487  counted_ptr<DalitzEvent> evtPtr(0);
488  nTries=0;
489  double maxW=-9999;
490  do{
491  evtPtr = try4EventWithPhaseSpace(maxW, mapping);
492  nTries++;
493  if(nTries > nearInfinity){
494  cout << "DalitzBWArea::makeEventWithPhaseSpace() "
495  << " giving up after " << nTries
496  << " unsuccessful tries." << endl;
497  return evtPtr;
498  }
499  }while(0 == evtPtr || _rnd->Rndm()*maxW > evtPtr->getWeight());
500  evtPtr->setWeight(1);
501 
502  return evtPtr;
503 }
DalitzBWArea::_rnd
TRandom * _rnd
Definition: DalitzBWArea.h:27
MINT::counted_ptr
Definition: counted_ptr.h:39
DalitzBWArea::try4EventWithPhaseSpace
MINT::counted_ptr< DalitzEvent > try4EventWithPhaseSpace(double &maxWeight, const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:772

◆ makeCoord() [1/2]

void DalitzBWArea::makeCoord ( int  i,
int  j 
)
private

Definition at line 111 of file DalitzBWArea.cpp.

111  {
112  DalitzCoordinate c(i, j);
113  counted_ptr<IGenFct> fcn(new FlatFct(c));
114  std::pair<DalitzCoordinate, counted_ptr<IGenFct> > p (c, fcn);
115  _coords[c.myKey()] = p;
116 }
MINT::counted_ptr< IGenFct >
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ makeCoord() [2/2]

void DalitzBWArea::makeCoord ( int  i,
int  j,
int  k 
)
private

Definition at line 117 of file DalitzBWArea.cpp.

117  {
118  DalitzCoordinate c(i, j, k);
119  counted_ptr<IGenFct> fcn(new FlatFct(c));
120  std::pair<DalitzCoordinate, counted_ptr<IGenFct> > p (c, fcn);
121  _coords[c.myKey()] = p;
122 }
MINT::counted_ptr< IGenFct >
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ makeCoords()

void DalitzBWArea::makeCoords ( )
private

Definition at line 103 of file DalitzBWArea.cpp.

103  {
104  _coords.clear();
105  if(_pat.numDaughters() >= 4)makeCoord(2,3,4);
106  makeCoord(1,2);
107  makeCoord(2,3);
108  if(_pat.numDaughters() >= 4) makeCoord(3,4);
109  if(_pat.numDaughters() >= 4) makeCoord(1,2,3);
110 }
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzBWArea::makeCoord
void makeCoord(int i, int j)
Definition: DalitzBWArea.cpp:111
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ makeMiMa()

void DalitzBWArea::makeMiMa ( )
private

Definition at line 129 of file DalitzBWArea.cpp.

129  {
130  bool dbThis=false;
131  double safetyFactor=1; // for debug only - otherwise 1
132  // double safetyFactor=1.00001;
133  s12().setMinMax( _pat.sijMin(1,2)/safetyFactor
134  , _pat.sijMax(1,2)*safetyFactor);
135  s23().setMinMax( _pat.sijMin(2,3)/safetyFactor
136  , _pat.sijMax(2,3)*safetyFactor);
137  if(dbThis) cout << "_pat.numDaughters() " << _pat.numDaughters() << endl;
138  if(_pat.numDaughters() >= 4){
139  t01().setMinMax( _pat.sijMin(2,3,4)/safetyFactor
140  , _pat.sijMax(2,3,4)*safetyFactor);
141  s34().setMinMax( _pat.sijMin(3,4)/safetyFactor
142  , _pat.sijMax(3,4)*safetyFactor);
143  t40().setMinMax( _pat.sijMin(1,2,3)/safetyFactor
144  , _pat.sijMax(1,2,3)*safetyFactor);
145  }
146  f_s12()->setLimits( _pat.sijMin(1,2)/safetyFactor
147  , _pat.sijMax(1,2)*safetyFactor);
148  f_s23()->setLimits( _pat.sijMin(2,3)/safetyFactor
149  , _pat.sijMax(2,3)*safetyFactor);
150  if(_pat.numDaughters() >= 4){
151  f_t01()->setLimits( _pat.sijMin(2,3,4)/safetyFactor
152  , _pat.sijMax(2,3,4)*safetyFactor);
153  f_s34()->setLimits( _pat.sijMin(3,4)/safetyFactor
154  , _pat.sijMax(3,4)*safetyFactor);
155  f_t40()->setLimits( _pat.sijMin(1,2,3)/safetyFactor
156  , _pat.sijMax(1,2,3)*safetyFactor);
157  }
158 
159 }
DalitzEventPattern::sijMax
double sijMax(const MINT::PolymorphVector< int > &indices) const
Definition: DalitzEventPattern.cpp:254
DalitzBWArea::t40
DalitzCoordinate & t40()
Definition: DalitzBWArea.h:139
DalitzBWArea::f_t40
MINT::counted_ptr< IGenFct > f_t40()
Definition: DalitzBWArea.h:171
DalitzBWArea::f_s23
MINT::counted_ptr< IGenFct > f_s23()
Definition: DalitzBWArea.h:165
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzBWArea::s34
DalitzCoordinate & s34()
Definition: DalitzBWArea.h:136
IGenFct::setLimits
virtual void setLimits(double sMin, double sMax)=0
DalitzBWArea::f_t01
MINT::counted_ptr< IGenFct > f_t01()
Definition: DalitzBWArea.h:159
DalitzBWArea::s23
DalitzCoordinate & s23()
Definition: DalitzBWArea.h:133
DalitzBWArea::f_s12
MINT::counted_ptr< IGenFct > f_s12()
Definition: DalitzBWArea.h:162
DalitzEventPattern::sijMin
double sijMin(const MINT::PolymorphVector< int > &indices) const
Definition: DalitzEventPattern.cpp:242
DalitzBWArea::s12
DalitzCoordinate & s12()
Definition: DalitzBWArea.h:130
DalitzBWArea::f_s34
MINT::counted_ptr< IGenFct > f_s34()
Definition: DalitzBWArea.h:168
DalitzBWArea::t01
DalitzCoordinate & t01()
Definition: DalitzBWArea.h:127
DalitzCoordinate::setMinMax
void setMinMax(double min, double max)
Definition: DalitzCoordinate.h:35

◆ mapP4() [1/2]

std::vector<TLorentzVector>& DalitzBWArea::mapP4 ( const std::vector< TLorentzVector > &  v_in,
const Permutation &  mapping,
std::vector< TLorentzVector > &  v_out 
) const
inlineprivate

Definition at line 76 of file DalitzBWArea.h.

79  {
80  unsigned int n = v_in.size();
81  if(mapping.isUnity()){
82  v_out = v_in;
83  return v_out;
84  }
85 
86  v_out.resize(n);
87  for(unsigned int i=0; i < n; i++){
88  int mappedIndex = mapping[i];
89  if(mappedIndex < 0 || mappedIndex + 1 > (int) n){
90  std::cout << "ERROR in DalitzBWArea::mapP4()"
91  << "\n Index out of range: " << mappedIndex
92  << " not in [0, " << n-1
93  << std::endl;
94  throw "index out of range.";
95  }
96  v_out[mappedIndex] = v_in[i];
97  }
98  return v_out;
99  }
Permutation::isUnity
bool isUnity() const
Definition: Permutation.cpp:47

◆ mapP4() [2/2]

std::vector<TLorentzVector> DalitzBWArea::mapP4 ( const std::vector< TLorentzVector > &  v_in,
const Permutation &  mapping 
) const
inlineprivate

Definition at line 100 of file DalitzBWArea.h.

102  {
103  std::vector<TLorentzVector> v_out(v_in.size());
104  return mapP4(v_in, mapping, v_out);
105  }
DalitzBWArea::mapP4
std::vector< TLorentzVector > & mapP4(const std::vector< TLorentzVector > &v_in, const Permutation &mapping, std::vector< TLorentzVector > &v_out) const
Definition: DalitzBWArea.h:76

◆ MC4Integral()

double DalitzBWArea::MC4Integral ( double &  prec) const
private

Definition at line 505 of file DalitzBWArea.cpp.

505  {
506  // for cross checks only
507  int N=1000000;
508  double sum=0;
509  double sumsq=0;
510  for(int i=0; i < N; i++){
511  counted_ptr<DalitzEvent> evtPtr(try4EventWithPhaseSpace());
512  if(0 == evtPtr) continue;
513  double val = evtPtr->getWeight();// * genValueRho(&(*evtPtr));
514  sum += val;
515  sumsq += val*val;
516  }
517  double mean = sum/((double)N);
518  double meansq = sumsq/((double)N);
519  double variance = (meansq - mean*mean)/((double)N);
520  double sigma = sqrt(variance);
521  precision = sigma/mean;
522 
523  return mean * integral();
524 
525  // ============== not used ==================
526  double xmi=0, xma=0, ymi=0, yma=0;
527  if(ResonanceConfigurationNumber() == 0){
528  xmi = sf(1,2).second->getRhoMi();
529  xma = sf(1,2).second->getRhoMa();
530  ymi = sf(3,4).second->getRhoMi();
531  yma = sf(3,4).second->getRhoMa();
532  }else{
533  xmi = sf(1,2,3).second->getRhoMi();
534  xma = sf(1,2,3).second->getRhoMa();
535  ymi = sf(1,2).second->getRhoMi();
536  yma = sf(1,2).second->getRhoMa();
537  }
538 
539  double intArea = (xma - xmi)*(yma - ymi);
540  if(intArea < 0){
541  cout << "DalitzBWArea::MC4Integral() negative area? "
542  << "\n x: ( " << xma << " - " << xmi << ") = "
543  << xma - xmi
544  << "\n y: ( " << yma << " - " << ymi << ") = "
545  << yma - ymi
546  << endl;
547  }
548 
549  cout << "DalitzBWArea::MC4Integral: returning "
550  << mean << " * " << intArea << endl;
551  return mean * intArea;
552 
553 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165
MINT::counted_ptr
Definition: counted_ptr.h:39
DalitzBWArea::try4EventWithPhaseSpace
MINT::counted_ptr< DalitzEvent > try4EventWithPhaseSpace(double &maxWeight, const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:772
DalitzBWArea::ResonanceConfigurationNumber
int ResonanceConfigurationNumber() const
Definition: DalitzBWArea.cpp:736
DalitzBWArea::integral
double integral() const
Definition: DalitzBWArea.cpp:397

◆ MC4IntegralNoTransform()

double DalitzBWArea::MC4IntegralNoTransform ( double &  prec) const
private

Definition at line 555 of file DalitzBWArea.cpp.

555  {
556  // for cross checks only.
557  int N=1000000;
558  double sum=0;
559  double sumsq=0;
560  for(int i=0; i < N; i++){
561  counted_ptr<DalitzEvent> evtPtr(tryFlat4EventWithPhaseSpace());
562  if(0 == evtPtr) continue;
563  double val = evtPtr->getWeight() * genValue(&(*evtPtr));
564  sum += val;
565  sumsq += val*val;
566  }
567  double mean = sum/((double)N);
568  double meansq = sumsq/((double)N);
569  double variance = (meansq - mean*mean)/((double)N);
570  prec = sqrt(variance)/mean;
571 
572  double xmi=0, xma=0, ymi=0, yma=0;
573  if(ResonanceConfigurationNumber() == 0){
574  xmi = sf(1,2).second->getSMi();
575  xma = sf(1,2).second->getSMa();
576  if(_pat.numDaughters()>=4){
577  ymi = sf(3,4).second->getSMi();
578  yma = sf(3,4).second->getSMa();
579  }
580  }else{
581  xmi = sf(1,2,3).second->getSMi();
582  xma = sf(1,2,3).second->getSMa();
583  if(_pat.numDaughters()>=4){
584  ymi = sf(1,2).second->getSMi();
585  yma = sf(1,2).second->getSMa();
586  }
587  }
588 
589  double intArea = 0;
590 
591  if(_pat.numDaughters()>=4){
592  intArea = (xma - xmi)*(yma - ymi);
593  }else{
594  intArea = (xma - xmi);
595  }
596  if(intArea < 0){
597  cout << "DalitzBWArea::MC4IntegralNoTransform() negative area? "
598  << "\n x: ( " << xma << " - " << xmi << ") = "
599  << xma - xmi
600  << "\n y: ( " << yma << " - " << ymi << ") = "
601  << yma - ymi
602  << endl;
603  }
604 
605  cout << "DalitzBWArea::MC4IntegralNoTransform: returning "
606  << mean << " * " << intArea << endl;
607  return mean * intArea;
608 
609 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
MINT::counted_ptr
Definition: counted_ptr.h:39
DalitzBWArea::genValue
double genValue(const IDalitzEvent *evtPtr) const
Definition: DalitzBWArea.cpp:297
DalitzBWArea::ResonanceConfigurationNumber
int ResonanceConfigurationNumber() const
Definition: DalitzBWArea.cpp:736
DalitzBWArea::tryFlat4EventWithPhaseSpace
MINT::counted_ptr< DalitzEvent > tryFlat4EventWithPhaseSpace(double &maxWeight, const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:1058

◆ operator=()

DalitzBWArea & DalitzBWArea::operator= ( const DalitzBWArea &  other)

Definition at line 47 of file DalitzBWArea.cpp.

47  {
48  _pat = other._pat;
49  _rnd = other._rnd;
50  _madeCMap = false;
51  _coords = other._coords;
52  _mappedCoords = other._mappedCoords;
53  _areaIntegral = other._areaIntegral;
54  _unWeightPs = other._unWeightPs;
55  return (*this);
56 }
DalitzBWArea::_mappedCoords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _mappedCoords
Definition: DalitzBWArea.h:121
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzBWArea::_rnd
TRandom * _rnd
Definition: DalitzBWArea.h:27
DalitzBWArea::_unWeightPs
bool _unWeightPs
Definition: DalitzBWArea.h:34
DalitzBWArea::_areaIntegral
double _areaIntegral
Definition: DalitzBWArea.h:33
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ print()

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

Definition at line 1434 of file DalitzBWArea.cpp.

1434  {
1435  os << "Area size: " << size() << endl;
1436  for(std::map<DalitzCoordKey, std::pair<DalitzCoordinate, counted_ptr<IGenFct> > >::const_iterator
1437  it =_coords.begin();
1438  it != _coords.end();
1439  it++){
1440 
1441  os << "\n " << it->second.first;
1442  }
1443 }
MINT::counted_ptr< IGenFct >
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119
DalitzBWArea::size
double size() const
Definition: DalitzBWArea.cpp:270

◆ ResonanceConfigurationNumber()

int DalitzBWArea::ResonanceConfigurationNumber ( ) const
private

Definition at line 736 of file DalitzBWArea.cpp.

736  {
737  if(_pat.numDaughters() <= 3) return 0; // only one configurtion: D->X(d1, d2), d3
738  if(_pat.numDaughters() > 4) return 0; // no known configuration
739 
740  // so it's 4 daughters, m1, m2, m3, m4.
741  // 0: treat as D->X Y, with X->m1, m2; Y->m3,m4
742  // 1: treat as D->X m4, with X->Y m3, with Y->m1, m2
743 
744  if((! sf(1,2).second->flat())
745  && (! sf(3,4).second->flat())){
746  return 0; // D-> 2 resonances
747  }else if((! sf(1,2).second->flat())
748  && sf(3,4).second->flat()
749  && sf(1,2,3).second->flat()){
750  return 0; // D-> 1 resonance + 2 bodies
751  }else if( sf(1,2).second->flat()
752  && sf(3,4).second->flat()
753  && sf(1,2,3).second->flat()){
754  return 0; // D-> 4 body (non-resonant)
755  }else{
756  return 1; // D-> resonance + 1 body; resonance->another resonance + 1 body
757  }
758 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
second
static const double second
Definition: CLHEPSystemOfUnits.h:115

◆ s12() [1/2]

DalitzCoordinate& DalitzBWArea::s12 ( )
inline

Definition at line 130 of file DalitzBWArea.h.

130  {
131  return sf(1,2).first;
132  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ s12() [2/2]

const DalitzCoordinate& DalitzBWArea::s12 ( ) const
inline

Definition at line 146 of file DalitzBWArea.h.

146  {
147  return sf(1,2).first;
148  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ s23() [1/2]

DalitzCoordinate& DalitzBWArea::s23 ( )
inline

Definition at line 133 of file DalitzBWArea.h.

133  {
134  return sf(2,3).first;
135  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ s23() [2/2]

const DalitzCoordinate& DalitzBWArea::s23 ( ) const
inline

Definition at line 149 of file DalitzBWArea.h.

149  {
150  return sf(2,3).first;
151  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ s34() [1/2]

DalitzCoordinate& DalitzBWArea::s34 ( )
inline

Definition at line 136 of file DalitzBWArea.h.

136  {
137  return sf(3,4).first;
138  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ s34() [2/2]

const DalitzCoordinate& DalitzBWArea::s34 ( ) const
inline

Definition at line 152 of file DalitzBWArea.h.

152  {
153  return sf(3,4).first;
154  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ setFcn()

void DalitzBWArea::setFcn ( const DalitzCoordinate &  c,
const MINT::counted_ptr< IGenFct > &  fcn 
)

Definition at line 124 of file DalitzBWArea.cpp.

124  {
125  fct->setLimits(sf(c).first.min(), sf(c).first.max());
126  sf(c).second = fct;
127 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ setPattern()

void DalitzBWArea::setPattern ( const DalitzEventPattern &  pat)

Definition at line 82 of file DalitzBWArea.cpp.

82  {
83  _pat = pat;
84  makeCoords();
85  if(_pat.numDaughters() != 4
86  && _pat.numDaughters() != 3
87  ){
88  cout << "DalitzBWArea::DalitzBWArea:"
89  << "\n Sorry, can only deal with 3 or 4 body decays. "
90  << "\n Please improve me so I can deal with decays like this: "
91  << pat << endl;
92  cout << " I'll crash now." << endl;
93  throw "DalitzBWArea: \"I'm no superman\"";
94  }
95  makeMiMa();
96 }
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25

◆ setRnd()

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

Definition at line 98 of file DalitzBWArea.cpp.

98  {
99  _rnd = rnd;
100  return true;
101 }
DalitzBWArea::_rnd
TRandom * _rnd
Definition: DalitzBWArea.h:27

◆ setUnWeightPs()

void DalitzBWArea::setUnWeightPs ( bool  doSo = true)
inline

Definition at line 125 of file DalitzBWArea.h.

125 {_unWeightPs = doSo;}
DalitzBWArea::_unWeightPs
bool _unWeightPs
Definition: DalitzBWArea.h:34

◆ sf() [1/6]

std::pair< DalitzCoordinate, counted_ptr< IGenFct > > & DalitzBWArea::sf ( const DalitzCoordinate &  c)
private

Definition at line 165 of file DalitzBWArea.cpp.

165  {
166  map<DalitzCoordKey, pair<DalitzCoordinate, counted_ptr<IGenFct> > >::iterator it = _coords.find(c.myKey());
167  if(it == _coords.end()){
168  cout << "ERROR in DalitzBWArea::sf"
169  << " unknown coordinate: " << c
170  << endl;
171  throw "CAN'T HANDLE THAT...";
172  }
173  return it->second;
174 }
DalitzCoordinate::myKey
const DalitzCoordKey & myKey() const
Definition: DalitzCoordinate.h:57
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ sf() [2/6]

const std::pair< DalitzCoordinate, counted_ptr< IGenFct > > & DalitzBWArea::sf ( const DalitzCoordinate &  c) const
private

Definition at line 177 of file DalitzBWArea.cpp.

177  {
178  map<DalitzCoordKey, pair<DalitzCoordinate, counted_ptr<IGenFct> > >::iterator it =
179  _coords.find(c.myKey());
180  if(it == _coords.end()){
181  cout << "ERROR in DalitzBWArea::sf"
182  << " unknown coordinate: " << c
183  << endl;
184  throw "CAN'T HANDLE THAT...";
185  }
186  return it->second;
187 }
DalitzCoordinate::myKey
const DalitzCoordKey & myKey() const
Definition: DalitzCoordinate.h:57
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ sf() [3/6]

std::pair< DalitzCoordinate, counted_ptr< IGenFct > > & DalitzBWArea::sf ( int  i,
int  j,
int  k 
)
private

Definition at line 190 of file DalitzBWArea.cpp.

190  {
191  DalitzCoordinate c(i,j,k);
192  return sf(c);
193 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ sf() [4/6]

const std::pair< DalitzCoordinate, counted_ptr< IGenFct > > & DalitzBWArea::sf ( int  i,
int  j,
int  k 
) const
private

Definition at line 196 of file DalitzBWArea.cpp.

196  {
197  DalitzCoordinate c(i,j,k);
198  return sf(c);
199 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ sf() [5/6]

std::pair< DalitzCoordinate, counted_ptr< IGenFct > > & DalitzBWArea::sf ( int  i,
int  j 
)
private

Definition at line 201 of file DalitzBWArea.cpp.

201  {
202  DalitzCoordinate c(i,j);
203  return sf(c);
204 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ sf() [6/6]

const std::pair< DalitzCoordinate, counted_ptr< IGenFct > > & DalitzBWArea::sf ( int  i,
int  j 
) const
private

Definition at line 206 of file DalitzBWArea.cpp.

206  {
207  DalitzCoordinate c(i,j);
208  return sf(c);
209 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ size()

double DalitzBWArea::size ( ) const

Definition at line 270 of file DalitzBWArea.cpp.

270  {
271  if(_pat.size() < 3) return 0;
272 
273  double p=1;
274 
275  /* cout << " DalitzBWArea::size(): returning..." << endl;
276  cout << " e.g. _t01 " << _t01 << endl;
277  cout << " e.g. _s12 " << _s12 << endl;
278  cout << "_madeCMap = " << _madeCMap << endl;
279  */
280  for(map<DalitzCoordKey, pair<DalitzCoordinate, counted_ptr<IGenFct> > >::const_iterator it = _coords.begin();
281  it != _coords.end(); it++){
282 
283  p *= (it->second.first.max() - it->second.first.min());
284  }
285  return fabs(p);
286 }
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
MINT::counted_ptr< IGenFct >
MINT::PolymorphVector::size
unsigned int size() const
Definition: PolymorphVector.h:46
DalitzBWArea::_coords
std::map< DalitzCoordKey, std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > > _coords
Definition: DalitzBWArea.h:119

◆ t01() [1/2]

DalitzCoordinate& DalitzBWArea::t01 ( )
inline

Definition at line 127 of file DalitzBWArea.h.

127  {
128  return sf(2,3,4).first;
129  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ t01() [2/2]

const DalitzCoordinate& DalitzBWArea::t01 ( ) const
inline

Definition at line 143 of file DalitzBWArea.h.

143  {
144  return sf(2,3,4).first;
145  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ t40() [1/2]

DalitzCoordinate& DalitzBWArea::t40 ( )
inline

Definition at line 139 of file DalitzBWArea.h.

139  {
140  return sf(1,2,3).first;
141  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ t40() [2/2]

const DalitzCoordinate& DalitzBWArea::t40 ( ) const
inline

Definition at line 155 of file DalitzBWArea.h.

155  {
156  return sf(1,2,3).first;
157  }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165

◆ try3Event()

counted_ptr< DalitzEvent > DalitzBWArea::try3Event ( const Permutation &  mapping = Permutation::unity()) const
private

Definition at line 451 of file DalitzBWArea.cpp.

451  {
452  bool dbThis=false;
453 
454  counted_ptr<DalitzEvent> evtPtr(try3EventWithPhaseSpace(mapping));
455  if(dbThis && 0 != evtPtr){
456  cout << "weight in DalitzBWArea::make3Event() "
457  << evtPtr->getWeight() << endl;
458  }
459  return evtPtr;
460 }
MINT::counted_ptr
Definition: counted_ptr.h:39
DalitzBWArea::try3EventWithPhaseSpace
MINT::counted_ptr< DalitzEvent > try3EventWithPhaseSpace(double &maxWeight, const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:1316

◆ try3EventWithPhaseSpace() [1/2]

counted_ptr< DalitzEvent > DalitzBWArea::try3EventWithPhaseSpace ( double &  maxWeight,
const Permutation &  mapping = Permutation::unity() 
) const
private

Definition at line 1316 of file DalitzBWArea.cpp.

1317  {
1318 
1319  bool dbThis=false;
1320 
1321  if(dbThis){
1322  cout << " trying to make event with phase space for "
1323  << *this
1324  << endl;
1325  }
1326  static TGenPhaseSpaceWithRnd gps;
1327  gps.setRandom(_rnd); // should not be necessary to repeat ch time, but no time to check and make sure
1328 
1329  TLorentzVector mumsP4;
1330  mumsP4.SetXYZM(0, 0, 0, _pat[0].mass());
1331 
1332  static vector<TLorentzVector> p4_final(4), p4_finalMapped(4);
1333  static TLorentzVector p4_inter[1];
1334  static DalitzEventPattern patMapped;
1335 
1336  counted_ptr<DalitzEvent> nullEvtPtr(0);
1337  counted_ptr<DalitzEvent> returnEvent(0);
1338 
1339  // all set up such that we know that only
1340  // s12 can be filled;
1341 
1342  maxWeight = 1.0;
1343 
1344  double s12 = sf(1,2).second->generate(_rnd);
1345  if(s12 < 0) return nullEvtPtr;
1346  if(sqrt(s12) < _pat[1].mass() + _pat[2].mass()) return nullEvtPtr;
1347 
1348  if(sqrt(s12) + _pat[3].mass() > _pat[0].mass())return nullEvtPtr;
1349 
1350  if(dbThis){
1351  cout << "made s12 " << s12 << endl;
1352  }
1353  double s12Min = sf(1,2).second->getCoordinate().min();
1354  double s12Max = sf(1,2).second->getCoordinate().max();
1355 
1356  Double_t dgt[2];
1357 
1358  // D->s12, m3
1359 
1360  dgt[0] = sqrt(s12);
1361  dgt[1] = _pat[3].mass();
1362  mumsP4.SetXYZM(0, 0, 0, _pat[0].mass());
1363  p4_final[0] = mumsP4;
1364 
1365  if(! gps.SetDecay(mumsP4, 2, dgt)) return nullEvtPtr;
1366  if(dbThis){
1367  cout << " set decay to " << mumsP4 << " -> "
1368  << dgt[0] << ", " << dgt[1]
1369  << endl;
1370  }
1371  gps.Generate();
1372  /* Note:
1373  usually, gps.Generate() produces weighted events, and
1374  to get un-weighted ones (i.e.events with weight 1) one
1375  needs to do:
1376  while(rnd->Rndm() > gps.Generate());
1377  However, 2-body decays all have weight 1.
1378 
1379  */
1380  p4_inter[0] = *(gps.GetDecay(0));
1381  p4_final[3] = *(gps.GetDecay(1));
1382 
1383  double ps = phaseSpaceIntegral2body(_pat[0].mass()
1384  , sqrt(s12)
1385  , _pat[3].mass());
1386  maxWeight *= phaseSpaceIntegral2body(_pat[0].mass()
1387  ,sqrt(s12Min)
1388  ,_pat[3].mass());
1389 
1390  // s12 -> m1, m2
1391  dgt[0] = _pat[1].mass();
1392  dgt[1] = _pat[2].mass();
1393  if(! gps.SetDecay(p4_inter[0], 2, dgt)) return nullEvtPtr;
1394  if(dbThis){
1395  cout << " set decay to " << p4_inter[0] << " -> "
1396  << dgt[0] << ", " << dgt[1]
1397  << endl;
1398  }
1399  gps.Generate();
1400  p4_final[1] = *(gps.GetDecay(0));
1401  p4_final[2] = *(gps.GetDecay(1));
1402 
1403  ps *= phaseSpaceIntegral2body(sqrt(s12)
1404  , _pat[1].mass()
1405  , _pat[2].mass());
1406  maxWeight *= phaseSpaceIntegral2body(sqrt(s12Max)
1407  , _pat[1].mass()
1408  , _pat[2].mass());
1409 
1410  mapP4(p4_final, mapping, p4_finalMapped);
1411  mapping.mapOrder(_pat, patMapped);
1412  counted_ptr<DalitzEvent> thisEvent(new DalitzEvent(patMapped
1413  , p4_finalMapped));
1414  if(thisEvent->phaseSpace() <= 0.0){
1415  cout << "WARNING in try3EventWithPhaseSpace: "
1416  << " made 'good' event with bad phase space"
1417  << " , sqrt(s12) " << sqrt(s12)
1418  << " , sqrt(s12) + m3 " << sqrt(s12) + _pat[3].mass()
1419  << " and m(D) " << _pat[0].mass()
1420  << endl;
1421  cout << " compare " << s12 << ", " << thisEvent->s(1,2)
1422  << endl;
1423  return nullEvtPtr;
1424  }
1425  thisEvent->setWeight(ps);
1426  returnEvent = thisEvent;
1427  if(dbThis){
1428  cout << " got new event with weight "
1429  << returnEvent->getWeight() << endl;
1430  }
1431  return returnEvent;
1432 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
phaseSpaceIntegral2body
double phaseSpaceIntegral2body(double mum, double d1, double d2)
Definition: phaseSpaceIntegrals.h:20
TGenPhaseSpaceWithRnd::setRandom
void setRandom(TRandom *rnd=gRandom)
Definition: TGenPhaseSpaceWithRnd.h:11
DalitzBWArea::_rnd
TRandom * _rnd
Definition: DalitzBWArea.h:27
DalitzBWArea::mapP4
std::vector< TLorentzVector > & mapP4(const std::vector< TLorentzVector > &v_in, const Permutation &mapping, std::vector< TLorentzVector > &v_out) const
Definition: DalitzBWArea.h:76
MINT::counted_ptr
Definition: counted_ptr.h:39
DalitzBWArea::s12
DalitzCoordinate & s12()
Definition: DalitzBWArea.h:130
Permutation::mapOrder
std::vector< T > & mapOrder(const std::vector< T > &in, std::vector< T > &out) const
Definition: Permutation.h:91

◆ try3EventWithPhaseSpace() [2/2]

counted_ptr< DalitzEvent > DalitzBWArea::try3EventWithPhaseSpace ( const Permutation &  mapping = Permutation::unity()) const
private

Definition at line 1312 of file DalitzBWArea.cpp.

1312  {
1313  double maxWeight;
1314  return try3EventWithPhaseSpace(maxWeight, mapping);
1315 }
DalitzBWArea::try3EventWithPhaseSpace
MINT::counted_ptr< DalitzEvent > try3EventWithPhaseSpace(double &maxWeight, const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:1316

◆ try4Event()

counted_ptr< DalitzEvent > DalitzBWArea::try4Event ( const Permutation &  mapping = Permutation::unity()) const
private

Definition at line 462 of file DalitzBWArea.cpp.

462  {
463  bool dbThis=false;
464 
465  counted_ptr<DalitzEvent> evtPtr(0);
466  if(unWeightPs()){
467  evtPtr = make4EventWithPhaseSpace(mapping);
468  }else{
469  evtPtr = try4EventWithPhaseSpace(mapping);
470  }
471  if(dbThis && 0 != evtPtr){
472  cout << "weight in DalitzBWArea::make4Event() "
473  << evtPtr->getWeight() << endl;
474  }
475  return evtPtr;
476  // return make4EventWithPhaseSpace();
477 
478 }
MINT::counted_ptr
Definition: counted_ptr.h:39
DalitzBWArea::unWeightPs
bool unWeightPs() const
Definition: DalitzBWArea.h:124
DalitzBWArea::try4EventWithPhaseSpace
MINT::counted_ptr< DalitzEvent > try4EventWithPhaseSpace(double &maxWeight, const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:772
DalitzBWArea::make4EventWithPhaseSpace
MINT::counted_ptr< DalitzEvent > make4EventWithPhaseSpace(const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:480

◆ try4EventWithPhaseSpace() [1/2]

counted_ptr< DalitzEvent > DalitzBWArea::try4EventWithPhaseSpace ( double &  maxWeight,
const Permutation &  mapping = Permutation::unity() 
) const
private

Definition at line 772 of file DalitzBWArea.cpp.

773  {
774  // return counted_ptr<DalitzEvent>(0);
775  //}
776 
777  bool dbThis=false;
778 
779  if(dbThis){
780  cout << " trying to make event with phase space for "
781  << *this
782  << endl;
783  }
784  // hoping to save some allocation time
785  // through use of static:
786  static TGenPhaseSpaceWithRnd gps(_rnd);
787  static TLorentzVector mumsP4;
788  mumsP4.SetXYZM(0, 0, 0, _pat[0].mass());
789 
790  static vector<TLorentzVector> p4_final(5);
791  static vector<TLorentzVector> p4_finalMapped(5);
792  static DalitzEventPattern patMapped;
793 
794  static TLorentzVector p4_inter[4];
795 
796  counted_ptr<DalitzEvent> nullEvtPtr(0);
797  counted_ptr<DalitzEvent> returnEvent(0);
798 
799  // all set up such that we know that only
800  // s123, s12 can be filled, or s12, s34
801  // (or single limits: s123 or s12)
802 
803  maxWeight = 1.0;
804 
805 
806  if(ResonanceConfigurationNumber() == 0){
807  // like D->K* rho, K*->Kpi, rho->pipi
808  if(dbThis) cout << " making s12, s34 configuration " << endl;
809 
810  double rho12 = sf(1,2).second->generateRho(_rnd);
811  double s12 = sf(1,2).second->coordTransformToS(rho12);
812 
813  if(s12 < 0) return nullEvtPtr;
814  double m12 = sqrt(s12);
815  if(m12 < _pat[1].mass() + _pat[2].mass()) return nullEvtPtr;
816 
817  double rho34 = sf(3,4).second->generateRho(_rnd);
818  double s34 = sf(3,4).second->coordTransformToS(rho34);
819 
820  if(s34 < 0) return nullEvtPtr;
821  double m34 = sqrt(s34);
822  if(m34 < _pat[3].mass() + _pat[4].mass()) return nullEvtPtr;
823  if(m12 + m34 > _pat[0].mass())return nullEvtPtr;
824 
825  double s12Min = sf(1,2).second->getCoordinate().min();
826  double s12Max = sf(1,2).second->getCoordinate().max();
827  double s34Min = sf(3,4).second->getCoordinate().min();
828  double s34Max = sf(3,4).second->getCoordinate().max();
829 
830 
831  double fct_12 = 1;//sf(1,2).second->transformedFctValue(rho12);
832  double fct_34 = 1;//sf(3,4).second->transformedFctValue(rho34);
833 
834  double fct_12Max = 1;//sf(1,2).second->transformedFctMax();
835  double fct_34Max = 1;//sf(3,4).second->transformedFctMax();
836 
837  maxWeight *= fct_12Max;
838  maxWeight *= fct_34Max;
839 
840  Double_t dgt[2];
841 
842  // D->s12, s34
843 
844  dgt[0] = m12;
845  dgt[1] = m34;
846  mumsP4.SetXYZM(0, 0, 0, _pat[0].mass());
847  p4_final[0] = mumsP4;
848 
849  if(! gps.SetDecay(mumsP4, 2, dgt)) return nullEvtPtr;
850  gps.Generate();
851  /* Note:
852  usually, gps.Generate() produces weighted events, and
853  to get un-weighted ones (i.e.events with weight 1) one
854  needs to do:
855  while(rnd->Rndm() > gps.Generate());
856  However, 2-body decays all have weight 1.
857 
858  */
859  p4_inter[0] = *(gps.GetDecay(0));
860  p4_inter[1] = *(gps.GetDecay(1));
861 
862  double ps1 = phaseSpaceIntegral2body(_pat[0].mass()
863  , m12
864  , m34);
865  maxWeight *= phaseSpaceIntegral2body(_pat[0].mass()
866  ,sqrt(s12Min)
867  ,sqrt(s34Min));
868  if(ps1 <=0 ) return nullEvtPtr;
869 
870  // s12 -> m1, m2
871  dgt[0] = _pat[1].mass();
872  dgt[1] = _pat[2].mass();
873  if(! gps.SetDecay(p4_inter[0], 2, dgt)) return nullEvtPtr;
874  gps.Generate();
875  p4_final[1] = *(gps.GetDecay(0));
876  p4_final[2] = *(gps.GetDecay(1));
877 
878  double ps2 = phaseSpaceIntegral2body(m12
879  , _pat[1].mass()
880  , _pat[2].mass());
881  maxWeight *= phaseSpaceIntegral2body(sqrt(s12Max)
882  , _pat[1].mass()
883  , _pat[2].mass());
884 
885  if(ps2 <=0)return nullEvtPtr;
886 
887  // s34 -> m3, m4
888  dgt[0] = _pat[3].mass();
889  dgt[1] = _pat[4].mass();
890  if(! gps.SetDecay(p4_inter[1], 2, dgt)) return nullEvtPtr;
891  gps.Generate();
892  p4_final[3] = *(gps.GetDecay(0));
893  p4_final[4] = *(gps.GetDecay(1));
894 
895  double ps3 = phaseSpaceIntegral2body(m34
896  , _pat[3].mass()
897  , _pat[4].mass());
898  maxWeight *= phaseSpaceIntegral2body(sqrt(s34Max)
899  , _pat[3].mass()
900  , _pat[4].mass());
901 
902  if(ps3 <= 0) return nullEvtPtr;
903 
904  double ps = ps1 * ps2 * ps3;
905 
906  double w = fct_12 * fct_34 * ps;
907 
908  mapP4(p4_final, mapping, p4_finalMapped);
909  mapping.mapOrder(_pat, patMapped);
910  counted_ptr<DalitzEvent> thisEvent(new DalitzEvent(patMapped, p4_finalMapped));
911  if(thisEvent->makeCalculate4BodyProps().phaseSpaceFactor() <= 0.0){
912  Calculate4BodyProps c4bp = thisEvent->makeCalculate4BodyProps();
913  cout << "WARNING in DalitzBWArea::try4EventWithPhaseSpace (0):"
914  << "\n made 'good' event with bad phase space "
915  << c4bp.showPhaseSpaceFactorCalculation()
916  << "\n" << *thisEvent
917  << "masses:\n";
918  for(int i=0; i <= 4; i++) cout << thisEvent->p(i).M() << ", ";
919  cout << "\n should be; \n";
920  for(int i=0; i <= 4; i++) cout << _pat[i].mass() << ", ";
921  cout <<"\n momentum sum ";
922  TLorentzVector psum(.0,.0,.0,.0);
923  for(int i=1; i <=4; i++) psum += thisEvent->p(i);
924  cout << psum;
925  cout << "\n ps1 " << ps1 << ", ps2 " << ps2 << ", ps3 " << ps3
926  << "; ps " << ps;
927  cout << endl;
928  return nullEvtPtr;
929  }
930  thisEvent->setWeight(w);
931  returnEvent = thisEvent;
932  }else{
933  if(dbThis) cout << " making s123, s12 configuration " << endl;
934  // like D->K1 pi, K1->K* pi, K*->Kpi
935 
936  double rho123 = sf(1,2,3).second->generateRho(_rnd);
937  double s123 = sf(1,2,3).second->coordTransformToS(rho123);
938  if(s123 < 0) return nullEvtPtr;
939  double m123 = sqrt(s123);
940  if(m123 + _pat[4].mass() > _pat[0].mass()) return nullEvtPtr;
941  if(m123 < _pat[1].mass() + _pat[2].mass() + _pat[3].mass()){
942  return nullEvtPtr;
943  }
944 
945  double rho12 = sf(1,2).second->generateRho(_rnd);
946  double s12 = sf(1,2).second->coordTransformToS(rho12);
947  if(s12 < 0) return nullEvtPtr;
948  double m12 = sqrt(s12);
949  if(m12 + _pat[3].mass() > m123)return nullEvtPtr;
950  if(m12 + _pat[3].mass() + _pat[4].mass() > _pat[0].mass() ){
951  return nullEvtPtr;
952  }
953  if(m12 < _pat[1].mass() + _pat[2].mass()) return nullEvtPtr;
954 
955  double s123Min = sf(1,2,3).second->getCoordinate().min();
956  double s123Max = sf(1,2,3).second->getCoordinate().max();
957  double s12Min = sf(1,2).second->getCoordinate().min();
958  double s12Max = sf(1,2).second->getCoordinate().max();
959 
960  Double_t dgt[2];
961 
962 
963  double fct_123 = 1;//sf(1,2,3).second->transformedFctValue(rho123);
964  double fct_12 = 1;//sf(1,2).second->transformedFctValue(rho12);
965  double fct_123Max = 1;//sf(1,2,3).second->transformedFctMax();
966  double fct_12Max = 1;//sf(1,2).second->transformedFctMax();
967 
968  maxWeight *= fct_123Max;
969  maxWeight *= fct_12Max;
970 
971  // D->s123, m4
972 
973  dgt[0] = m123;
974  dgt[1] = _pat[4].mass();
975  mumsP4.SetXYZM(0, 0, 0, _pat[0].mass());
976  p4_final[0] = mumsP4;
977 
978  if(! gps.SetDecay(mumsP4, 2, dgt)) return nullEvtPtr;
979  gps.Generate();
980  p4_inter[0] = *(gps.GetDecay(0));
981  p4_final[4] = *(gps.GetDecay(1));
982 
983  double ps1 = phaseSpaceIntegral2body(_pat[0].mass()
984  , m123
985  , _pat[4].mass());
986  maxWeight *= phaseSpaceIntegral2body(_pat[0].mass()
987  ,sqrt(s123Min)
988  , _pat[4].mass());
989 
990  if(ps1 <=0) return nullEvtPtr;
991  // s123 -> s12, m3
992  dgt[0] = m12;
993  dgt[1] = _pat[3].mass();
994  if(! gps.SetDecay(p4_inter[0], 2, dgt)) return nullEvtPtr;
995  gps.Generate();
996  p4_inter[1] = *(gps.GetDecay(0));
997  p4_final[3] = *(gps.GetDecay(1));
998 
999  double ps2 = phaseSpaceIntegral2body(m123
1000  , m12
1001  , _pat[3].mass());
1002  maxWeight *= phaseSpaceIntegral2body(sqrt(s123Max)
1003  ,sqrt(s12Min)
1004  , _pat[3].mass());
1005  if(ps2 <=0) return nullEvtPtr;
1006 
1007  // s12 -> m1, m2
1008  dgt[0] = _pat[1].mass();
1009  dgt[1] = _pat[2].mass();
1010  if(! gps.SetDecay(p4_inter[1], 2, dgt)) return nullEvtPtr;
1011  gps.Generate();
1012  p4_final[1] = *(gps.GetDecay(0));
1013  p4_final[2] = *(gps.GetDecay(1));
1014 
1015  double ps3 = phaseSpaceIntegral2body(m12
1016  , _pat[1].mass()
1017  , _pat[2].mass());
1018  maxWeight *= phaseSpaceIntegral2body(sqrt(s12Max)
1019  , _pat[1].mass()
1020  , _pat[2].mass());
1021 
1022  if(ps3 <=0) return nullEvtPtr;
1023 
1024  double ps = ps1 * ps2 * ps3;
1025 
1026  double w = fct_123 * fct_12 * ps;
1027 
1028  mapP4(p4_final, mapping, p4_finalMapped);
1029  mapping.mapOrder(_pat, patMapped);
1030  counted_ptr<DalitzEvent> thisEvent(new DalitzEvent(patMapped
1031  , p4_finalMapped));
1032  if(thisEvent->makeCalculate4BodyProps().phaseSpaceFactor() <= 0.0){
1033  cout << "WARNING in DalitzBWArea::try4EventWithPhaseSpace (2): "
1034  << " made 'good' event with bad phase space"
1035  << " for sqrt(s123) " << sqrt(s123) << " = " << m123
1036  << " , sqrt(s12) " << sqrt(s12) << " = " << m12
1037  << " , sqrt(s12) + m3 " << sqrt(s12) + _pat[3].mass()
1038  << " , sqrt(s123) + m4 " << sqrt(s123) + _pat[4].mass()
1039  << " and m(D) " << _pat[0].mass()
1040  << endl;
1041  cout << " compare " << s123 << ", " << thisEvent->t(4,0)
1042  << " compare " << s12 << ", " << thisEvent->s(1,2)
1043  << endl;
1044  cout << "compare mapping: " << endl;
1045  for(unsigned int i=0; i < p4_final.size() && i < p4_finalMapped.size(); i++){
1046  cout << "premapping " << p4_final[i] << " post: " << p4_finalMapped[i] << endl;
1047  }
1048  return nullEvtPtr;
1049  }
1050  thisEvent->setWeight(w);
1051  returnEvent = thisEvent;
1052  }
1053  if(dbThis)cout << " got new event with weight "
1054  << returnEvent->getWeight() << endl;
1055  return returnEvent;
1056 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzBWArea::s34
DalitzCoordinate & s34()
Definition: DalitzBWArea.h:136
phaseSpaceIntegral2body
double phaseSpaceIntegral2body(double mum, double d1, double d2)
Definition: phaseSpaceIntegrals.h:20
DalitzBWArea::_rnd
TRandom * _rnd
Definition: DalitzBWArea.h:27
DalitzBWArea::mapP4
std::vector< TLorentzVector > & mapP4(const std::vector< TLorentzVector > &v_in, const Permutation &mapping, std::vector< TLorentzVector > &v_out) const
Definition: DalitzBWArea.h:76
MINT::counted_ptr
Definition: counted_ptr.h:39
DalitzBWArea::s12
DalitzCoordinate & s12()
Definition: DalitzBWArea.h:130
DalitzBWArea::ResonanceConfigurationNumber
int ResonanceConfigurationNumber() const
Definition: DalitzBWArea.cpp:736
Permutation::mapOrder
std::vector< T > & mapOrder(const std::vector< T > &in, std::vector< T > &out) const
Definition: Permutation.h:91

◆ try4EventWithPhaseSpace() [2/2]

counted_ptr< DalitzEvent > DalitzBWArea::try4EventWithPhaseSpace ( const Permutation &  mapping = Permutation::unity()) const
private

Definition at line 760 of file DalitzBWArea.cpp.

760  {
761  double maxWeight;
762  return try4EventWithPhaseSpace(maxWeight, mapping);
763 }
DalitzBWArea::try4EventWithPhaseSpace
MINT::counted_ptr< DalitzEvent > try4EventWithPhaseSpace(double &maxWeight, const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:772

◆ tryEventForOwner()

counted_ptr< DalitzEvent > DalitzBWArea::tryEventForOwner ( const Permutation &  mapping = Permutation::unity()) const

Definition at line 415 of file DalitzBWArea.cpp.

415  {
416  bool dbThis=false;
417  counted_ptr<DalitzEvent> evtPtr(0);
418  if(_pat.numDaughters() == 3){
419  evtPtr = try3Event(mapping);
420  }else if(_pat.numDaughters() == 4){
421  evtPtr = try4Event(mapping);
422  if(dbThis && 0 != evtPtr){
423  cout << " DalitzBWArea::makeEventForOwner() "
424  << " returning event with weight "
425  << evtPtr->getWeight()
426  << endl;
427  }
428  }else{
429  cout << "ERROR in DalitzBWArea::tryEventForOwner() can only make events"
430  << " with 3 or 4 daughters. You want : " << _pat
431  << endl;
432  return counted_ptr<DalitzEvent>(0);
433  }
434 
435  if(dbThis && 0 != evtPtr) cout << "Event before P-con " << *evtPtr << endl;
436  if(0 != evtPtr && _pat[0] < 0) evtPtr->P_conjugateYourself();
437  // the above ensures that, for the same random seed,
438  // identical but CP conjugate events are generated
439  // for D->f and Dbar->fbar.
440  // Note that this step of the event generation is
441  // completely P-even, and the event generation would
442  // still be correct without this P-conjugation. The
443  // crucial P-senstive step is the reweighting applied
444  // later, which will then take into account the full
445  // amplitude model. The P-conjugation here is just to keep the
446  // random numbers in sync between CP conjugate event generations.
447  if(dbThis && 0 != evtPtr) cout << "Event after P-con " << *evtPtr << endl;
448  return evtPtr;
449 }
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzBWArea::try3Event
MINT::counted_ptr< DalitzEvent > try3Event(const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:451
MINT::counted_ptr
Definition: counted_ptr.h:39
DalitzBWArea::try4Event
MINT::counted_ptr< DalitzEvent > try4Event(const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:462

◆ tryFlat4EventWithPhaseSpace() [1/2]

counted_ptr< DalitzEvent > DalitzBWArea::tryFlat4EventWithPhaseSpace ( double &  maxWeight,
const Permutation &  mapping = Permutation::unity() 
) const
private

Definition at line 1058 of file DalitzBWArea.cpp.

1058  {
1059  // return counted_ptr<DalitzEvent>(0);
1060  //}
1061 
1062  bool dbThis=false;
1063 
1064  if(dbThis){
1065  cout << " trying to make event with phase space for "
1066  << *this
1067  << endl;
1068  }
1069  TGenPhaseSpaceWithRnd gps(_rnd);
1070 
1071  TLorentzVector mumsP4;
1072  mumsP4.SetXYZM(0, 0, 0, _pat[0].mass());
1073 
1074  static vector<TLorentzVector> p4_final(5), p4_finalMapped(5);
1075  static TLorentzVector p4_inter[4];
1076  static DalitzEventPattern patMapped;
1077 
1078  counted_ptr<DalitzEvent> nullEvtPtr(0);
1079  counted_ptr<DalitzEvent> returnEvent(0);
1080 
1081  // all set up such that we know that only
1082  // s123, s12 can be filled, or s12, s34
1083  // (or single limits: s123 or s12)
1084 
1085  maxWeight = 1.0;
1086 
1087 
1088  if(ResonanceConfigurationNumber() == 0){
1089  if(dbThis) cout << " making s12, s34 configuration " << endl;
1090 
1091  double s12Min = sf(1,2).second->getCoordinate().min();
1092  double s12Max = sf(1,2).second->getCoordinate().max();
1093  double s34Min = sf(3,4).second->getCoordinate().min();
1094  double s34Max = sf(3,4).second->getCoordinate().max();
1095 
1096  double s12 = s12Min + _rnd->Rndm()*(s12Max - s12Min);
1097 
1098  if(s12 < 0) return nullEvtPtr;
1099  if(sqrt(s12) < _pat[1].mass() + _pat[2].mass()) return nullEvtPtr;
1100 
1101  double s34 = s34Min + _rnd->Rndm()*(s34Max - s34Min);
1102 
1103  if(s34 < 0) return nullEvtPtr;
1104  if(sqrt(s34) < _pat[3].mass() + _pat[4].mass()) return nullEvtPtr;
1105  if(sqrt(s12) + sqrt(s34) > _pat[0].mass())return nullEvtPtr;
1106 
1107  double fct_12 = 1;//sf(1,2).second->transformedFctValue(rho12);
1108  double fct_34 = 1;//sf(3,4).second->transformedFctValue(rho34);
1109 
1110  double fct_12Max = 1;//sf(1,2).second->transformedFctMax();
1111  double fct_34Max = 1;//sf(3,4).second->transformedFctMax();
1112 
1113  maxWeight *= fct_12Max;
1114  maxWeight *= fct_34Max;
1115 
1116  Double_t dgt[2];
1117 
1118  // D->s12, s34
1119 
1120  dgt[0] = sqrt(s12);
1121  dgt[1] = sqrt(s34);
1122  mumsP4.SetXYZM(0, 0, 0, _pat[0].mass());
1123  p4_final[0] = mumsP4;
1124 
1125  if(! gps.SetDecay(mumsP4, 2, dgt)) return nullEvtPtr;
1126  gps.Generate();
1127  /* Note:
1128  usually, gps.Generate() produces weighted events, and
1129  to get un-weighted ones (i.e.events with weight 1) one
1130  needs to do:
1131  while(rnd->Rndm() > gps.Generate());
1132  However, 2-body decays all have weight 1.
1133 
1134  */
1135  p4_inter[0] = *(gps.GetDecay(0));
1136  p4_inter[1] = *(gps.GetDecay(1));
1137 
1138  double ps1 = phaseSpaceIntegral2body(_pat[0].mass()
1139  , sqrt(s12)
1140  , sqrt(s34));
1141  maxWeight *= phaseSpaceIntegral2body(_pat[0].mass()
1142  ,sqrt(s12Min)
1143  ,sqrt(s34Min));
1144  if(ps1 <=0 ) return nullEvtPtr;
1145 
1146  // s12 -> m1, m2
1147  dgt[0] = _pat[1].mass();
1148  dgt[1] = _pat[2].mass();
1149  if(! gps.SetDecay(p4_inter[0], 2, dgt)) return nullEvtPtr;
1150  gps.Generate();
1151  p4_final[1] = *(gps.GetDecay(0));
1152  p4_final[2] = *(gps.GetDecay(1));
1153 
1154  double ps2 = phaseSpaceIntegral2body(sqrt(s12)
1155  , _pat[1].mass()
1156  , _pat[2].mass());
1157  maxWeight *= phaseSpaceIntegral2body(sqrt(s12Max)
1158  , _pat[1].mass()
1159  , _pat[2].mass());
1160 
1161  if(ps2 <=0)return nullEvtPtr;
1162 
1163  // s34 -> m3, m4
1164  dgt[0] = _pat[3].mass();
1165  dgt[1] = _pat[4].mass();
1166  if(! gps.SetDecay(p4_inter[1], 2, dgt)) return nullEvtPtr;
1167  gps.Generate();
1168  p4_final[3] = *(gps.GetDecay(0));
1169  p4_final[4] = *(gps.GetDecay(1));
1170 
1171  double ps3 = phaseSpaceIntegral2body(sqrt(s34)
1172  , _pat[3].mass()
1173  , _pat[4].mass());
1174  maxWeight *= phaseSpaceIntegral2body(sqrt(s34Max)
1175  , _pat[3].mass()
1176  , _pat[4].mass());
1177 
1178  if(ps3 <= 0) return nullEvtPtr;
1179 
1180  double ps = ps1 * ps2 * ps3;
1181 
1182  double w = fct_12 * fct_34 * ps;
1183 
1184  mapP4(p4_final, mapping, p4_finalMapped);
1185  mapping.mapOrder(_pat, patMapped);
1186  counted_ptr<DalitzEvent> thisEvent(new DalitzEvent(patMapped, p4_finalMapped));
1187  if(thisEvent->makeCalculate4BodyProps().phaseSpaceFactor() <= 0.0){
1188  cout << "WARNING in tryFlat4EventWithPhaseSpace: "
1189  << " made 'good' event with bad phase space"
1190  << endl;
1191  return nullEvtPtr;
1192  }
1193  thisEvent->setWeight(w);
1194  returnEvent = thisEvent;
1195  }else{
1196  if(dbThis) cout << " making s123, s12 configuration " << endl;
1197 
1198  double s123Min = sf(1,2,3).second->getCoordinate().min();
1199  double s123Max = sf(1,2,3).second->getCoordinate().max();
1200  double s12Min = sf(1,2).second->getCoordinate().min();
1201  double s12Max = sf(1,2).second->getCoordinate().max();
1202 
1203  double s123 = s123Min + _rnd->Rndm()*(s123Max - s123Min);
1204 
1205  if(s123 < 0) return nullEvtPtr;
1206  if(sqrt(s123) + _pat[4].mass() > _pat[0].mass()) return nullEvtPtr;
1207  if(sqrt(s123) <
1208  _pat[1].mass() + _pat[2].mass() + _pat[3].mass()) return nullEvtPtr;
1209 
1210 
1211  double s12 = s12Min + _rnd->Rndm()*(s12Max - s12Min);
1212 
1213  if(s12 < 0) return nullEvtPtr;
1214  if(sqrt(s12) + _pat[3].mass() > sqrt(s123))return nullEvtPtr;
1215  if(sqrt(s12) + _pat[3].mass() + _pat[4].mass() > _pat[0].mass() )return nullEvtPtr;
1216  if(sqrt(s12) < _pat[1].mass() + _pat[2].mass()) return nullEvtPtr;
1217 
1218 
1219  Double_t dgt[2];
1220 
1221 
1222  double fct_123 = 1;//sf(1,2,3).second->transformedFctValue(rho123);
1223  double fct_12 = 1;//sf(1,2).second->transformedFctValue(rho12);
1224  double fct_123Max = 1;//sf(1,2,3).second->transformedFctMax();
1225  double fct_12Max = 1;//sf(1,2).second->transformedFctMax();
1226 
1227  maxWeight *= fct_123Max;
1228  maxWeight *= fct_12Max;
1229 
1230  // D->s123, m4
1231 
1232  dgt[0] = sqrt(s123);
1233  dgt[1] = _pat[4].mass();
1234  mumsP4.SetXYZM(0, 0, 0, _pat[0].mass());
1235  p4_final[0] = mumsP4;
1236 
1237  if(! gps.SetDecay(mumsP4, 2, dgt)) return nullEvtPtr;
1238  gps.Generate();
1239  p4_inter[0] = *(gps.GetDecay(0));
1240  p4_final[4] = *(gps.GetDecay(1));
1241 
1242  double ps1 = phaseSpaceIntegral2body(_pat[0].mass()
1243  , sqrt(s123)
1244  , _pat[4].mass());
1245  maxWeight *= phaseSpaceIntegral2body(_pat[0].mass()
1246  ,sqrt(s123Min)
1247  , _pat[4].mass());
1248 
1249  if(ps1 <=0) return nullEvtPtr;
1250  // s123 -> s12, m3
1251  dgt[0] = sqrt(s12);
1252  dgt[1] = _pat[3].mass();
1253  if(! gps.SetDecay(p4_inter[0], 2, dgt)) return nullEvtPtr;
1254  gps.Generate();
1255  p4_inter[1] = *(gps.GetDecay(0));
1256  p4_final[3] = *(gps.GetDecay(1));
1257 
1258  double ps2 = phaseSpaceIntegral2body(sqrt(s123)
1259  , sqrt(s12)
1260  , _pat[3].mass());
1261  maxWeight *= phaseSpaceIntegral2body(sqrt(s123Max)
1262  ,sqrt(s12Min)
1263  , _pat[3].mass());
1264  if(ps2 <=0) return nullEvtPtr;
1265 
1266  // s12 -> m1, m2
1267  dgt[0] = _pat[1].mass();
1268  dgt[1] = _pat[2].mass();
1269  if(! gps.SetDecay(p4_inter[1], 2, dgt)) return nullEvtPtr;
1270  gps.Generate();
1271  p4_final[1] = *(gps.GetDecay(0));
1272  p4_final[2] = *(gps.GetDecay(1));
1273 
1274  double ps3 = phaseSpaceIntegral2body(sqrt(s12)
1275  , _pat[1].mass()
1276  , _pat[2].mass());
1277  maxWeight *= phaseSpaceIntegral2body(sqrt(s12Max)
1278  , _pat[1].mass()
1279  , _pat[2].mass());
1280 
1281  if(ps3 <=0) return nullEvtPtr;
1282 
1283  double ps = ps1 * ps2 * ps3;
1284 
1285  double w = fct_123 * fct_12 * ps;
1286 
1287  mapP4(p4_final, mapping, p4_finalMapped);
1288  mapping.mapOrder(_pat, patMapped);
1289 
1290  counted_ptr<DalitzEvent> thisEvent(new DalitzEvent(patMapped, p4_finalMapped));
1291  if(thisEvent->makeCalculate4BodyProps().phaseSpaceFactor() <= 0.0){
1292  cout << "Warning in tryFlat4EventWithPhaseSpace "
1293  << " made 'good' event with bad phase space"
1294  << " for sqrt(s123) " << sqrt(s123)
1295  << " , sqrt(s12) " << sqrt(s12)
1296  << " , sqrt(s12) + m3 " << sqrt(s12) + _pat[3].mass()
1297  << " , sqrt(s123) + m4 " << sqrt(s123) + _pat[4].mass()
1298  << " and m(D) " << _pat[0].mass()
1299  << endl;
1300  cout << " compare " << s123 << ", " << thisEvent->t(4,0)
1301  << " compare " << s12 << ", " << thisEvent->s(1,2)
1302  << endl;
1303  return nullEvtPtr;
1304  }
1305  thisEvent->setWeight(w);
1306  returnEvent = thisEvent;
1307  }
1308  if(dbThis)cout << " got new event with weight " << returnEvent->getWeight() << endl;
1309  return returnEvent;
1310 }
DalitzBWArea::sf
std::pair< DalitzCoordinate, MINT::counted_ptr< IGenFct > > & sf(const DalitzCoordinate &c)
Definition: DalitzBWArea.cpp:165
DalitzBWArea::_pat
DalitzEventPattern _pat
Definition: DalitzBWArea.h:25
DalitzBWArea::s34
DalitzCoordinate & s34()
Definition: DalitzBWArea.h:136
phaseSpaceIntegral2body
double phaseSpaceIntegral2body(double mum, double d1, double d2)
Definition: phaseSpaceIntegrals.h:20
DalitzBWArea::_rnd
TRandom * _rnd
Definition: DalitzBWArea.h:27
DalitzBWArea::mapP4
std::vector< TLorentzVector > & mapP4(const std::vector< TLorentzVector > &v_in, const Permutation &mapping, std::vector< TLorentzVector > &v_out) const
Definition: DalitzBWArea.h:76
MINT::counted_ptr
Definition: counted_ptr.h:39
DalitzBWArea::s12
DalitzCoordinate & s12()
Definition: DalitzBWArea.h:130
DalitzBWArea::ResonanceConfigurationNumber
int ResonanceConfigurationNumber() const
Definition: DalitzBWArea.cpp:736
Permutation::mapOrder
std::vector< T > & mapOrder(const std::vector< T > &in, std::vector< T > &out) const
Definition: Permutation.h:91

◆ tryFlat4EventWithPhaseSpace() [2/2]

counted_ptr< DalitzEvent > DalitzBWArea::tryFlat4EventWithPhaseSpace ( const Permutation &  mapping = Permutation::unity()) const
private

Definition at line 765 of file DalitzBWArea.cpp.

765  {
766  double maxWeight;
767  return tryFlat4EventWithPhaseSpace(maxWeight, mapping);
768 }
DalitzBWArea::tryFlat4EventWithPhaseSpace
MINT::counted_ptr< DalitzEvent > tryFlat4EventWithPhaseSpace(double &maxWeight, const Permutation &mapping=Permutation::unity()) const
Definition: DalitzBWArea.cpp:1058

◆ unWeightPs()

bool DalitzBWArea::unWeightPs ( ) const
inline

Definition at line 124 of file DalitzBWArea.h.

124 { return _unWeightPs;}
DalitzBWArea::_unWeightPs
bool _unWeightPs
Definition: DalitzBWArea.h:34

Member Data Documentation

◆ _areaIntegral

double DalitzBWArea::_areaIntegral
mutableprivate

Definition at line 33 of file DalitzBWArea.h.

◆ _coords

std::map<DalitzCoordKey , std::pair<DalitzCoordinate, MINT::counted_ptr<IGenFct> > > DalitzBWArea::_coords
mutable

Definition at line 119 of file DalitzBWArea.h.

◆ _madeCMap

bool DalitzBWArea::_madeCMap
mutableprivate

Definition at line 31 of file DalitzBWArea.h.

◆ _mappedCoords

std::map<DalitzCoordKey , std::pair<DalitzCoordinate, MINT::counted_ptr<IGenFct> > > DalitzBWArea::_mappedCoords
mutable

Definition at line 121 of file DalitzBWArea.h.

◆ _pat

DalitzEventPattern DalitzBWArea::_pat
private

Definition at line 25 of file DalitzBWArea.h.

◆ _rnd

TRandom* DalitzBWArea::_rnd
private

Definition at line 27 of file DalitzBWArea.h.

◆ _unWeightPs

bool DalitzBWArea::_unWeightPs
private

Definition at line 34 of file DalitzBWArea.h.

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.15