[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenBase / EvtPto3PAmpFactory.cpp

//-----------------------------------------------------------------------
// File and Version Information: 
//      $Id: EvtPto3PAmpFactory.cpp,v 1.3 2009-03-16 15:44:04 robbep Exp $
// 
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations. If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information:
//      Copyright (C) 1998 Caltech, UCSB
//
// Module creator:
//      Alexei Dvoretskii, Caltech, 2001-2002.
//-----------------------------------------------------------------------
#include "EvtGenBase/EvtPatches.hh"

// AmpFactory for building a P -> 3P decay
// (pseudoscalar to three pseudoscalars)

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtComplex.hh"
#include "EvtGenBase/EvtCyclic3.hh"
#include "EvtGenBase/EvtSpinType.hh"
#include "EvtGenBase/EvtPto3PAmp.hh"
#include "EvtGenBase/EvtNonresonantAmp.hh"
#include "EvtGenBase/EvtFlatAmp.hh"
#include "EvtGenBase/EvtLASSAmp.hh"
#include "EvtGenBase/EvtPto3PAmpFactory.hh"
#include "EvtGenBase/EvtPropBreitWigner.hh"
#include "EvtGenBase/EvtPropFlatte.hh"
#include "EvtGenBase/EvtPropBreitWignerRel.hh"
#include "EvtGenBase/EvtDalitzResPdf.hh"
#include "EvtGenBase/EvtDalitzFlatPdf.hh"

using namespace EvtCyclic3;
#include <iostream>

void EvtPto3PAmpFactory::processAmp(EvtComplex c, std::vector<std::string> vv, bool conj)
{
  if(_verbose) {
    
    printf("Make %samplitude\n",conj ? "CP conjugate" : "");
    unsigned i;
    for(i=0;i<vv.size();i++) printf("%s\n",vv[i].c_str());
    printf("\n");
  }
  
  EvtAmplitude<EvtDalitzPoint>* amp = 0;
  EvtPdf<EvtDalitzPoint>* pdf = 0;
  std::string name;
  Pair pairRes=AB;

  size_t i;
  /*
         Experimental amplitudes
  */
  if(vv[0] == "PHASESPACE") {
    
    pdf = new EvtDalitzFlatPdf(_dp);
    amp = new EvtFlatAmp<EvtDalitzPoint>();
    name = "NR";
  }
  else if (!vv[0].find("NONRES")) {
    double alpha=0;
    EvtPto3PAmp::NumType typeNRes=EvtPto3PAmp::NONRES;
    if      (vv[0]=="NONRES_LIN") {
      typeNRes=EvtPto3PAmp::NONRES_LIN;
      pairRes=strToPair(vv[1].c_str());
    }
    else if (vv[0]=="NONRES_EXP") {
      typeNRes=EvtPto3PAmp::NONRES_EXP;
      pairRes = strToPair(vv[1].c_str());
      alpha   = strtod(vv[2].c_str(),0);
    }
    else assert(0);
    pdf = new EvtDalitzFlatPdf(_dp);
    amp = new EvtNonresonantAmp( &_dp, typeNRes, pairRes, alpha);  
  }
  else if (vv[0]=="LASS" || vv[0]=="LASS_ELASTIC" || vv[0]=="LASS_RESONANT") {
    pairRes = strToPair(vv[1].c_str());
    double m0 = strtod(vv[2].c_str(),0);
    double g0 = strtod(vv[3].c_str(),0);
    double a  = strtod(vv[4].c_str(),0);
    double r  = strtod(vv[5].c_str(),0);
    double cutoff  = strtod(vv[6].c_str(),0);
    pdf = new EvtDalitzResPdf(_dp,m0,g0,pairRes);
    amp = new EvtLASSAmp( &_dp, pairRes, m0, g0, a, r, cutoff, vv[0]);
  }

  /*
      Resonant amplitudes
  */
  else if(vv[0] == "RESONANCE") {
    EvtPto3PAmp* partAmp = 0;
      
    // RESONANCE stanza
    
    pairRes = strToPair(vv[1].c_str());
    EvtSpinType::spintype spinR = EvtSpinType::SCALAR;
    double mR, gR;      
    name = vv[2];
    EvtId resId = EvtPDL::getId(vv[2]);
    if(_verbose) printf("Particles %s form %sresonance %s\n",
			vv[1].c_str(),vv[2].c_str(), conj ? "(conj) " : "");

    // If no valid particle name is given, assume that 
    // it is the spin, the mass and the width of the particle.
      
    if(resId.getId() == -1) {
	
      switch(atoi(vv[2].c_str())) {
	
      case 0: { spinR = EvtSpinType::SCALAR; break; }
      case 1: { spinR = EvtSpinType::VECTOR; break; }
      case 2: { spinR = EvtSpinType::TENSOR; break; }
      case 3: { spinR = EvtSpinType::SPIN3; break; }
      case 4: { spinR = EvtSpinType::SPIN4; break; }
      default: { assert(0); break; }
      }
	
      mR = strtod(vv[3].c_str(),0);
      gR = strtod(vv[4].c_str(),0);
      i = 4;
    }
    else {
      
      // For a valid particle get spin, mass and width
      
      spinR = EvtPDL::getSpinType(resId);
      mR = EvtPDL::getMeanMass(resId);
      gR = EvtPDL::getWidth(resId);
      i = 2;
      
      // It's possible to specify mass and width of a particle 
      // explicitly
      
      if(vv[3] != "ANGULAR") {
	
	if(_verbose) 
	  printf("Setting m(%s)=%s g(%s)=%s\n",
		 vv[2].c_str(),vv[3].c_str(),vv[2].c_str(),vv[4].c_str());

	mR = strtod(vv[3].c_str(),0);
	gR = strtod(vv[4].c_str(),0);
	i = 4;
      }
    }
    
    // ANGULAR stanza
    
    if(vv[++i] != "ANGULAR") {

      printf("%s instead of ANGULAR\n",vv[i].c_str());
      exit(0);
    }
    Pair pairAng = strToPair(vv[++i].c_str());
    if(_verbose) printf("Angle is measured between particles %s\n",vv[i].c_str());
      
    // TYPE stanza
    
    std::string typeName = vv[++i];
    assert(typeName == "TYPE");
    std::string type = vv[++i];
    if(_verbose) printf("Propagator type %s\n",vv[i].c_str());
    
    if(type == "NBW") {      

      EvtPropBreitWigner prop(mR,gR);
      partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::NBW);
    }
    else if(type == "RBW_ZEMACH") {
      
      EvtPropBreitWignerRel prop(mR,gR);
      partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_ZEMACH);
    }
    else if(type == "RBW_KUEHN") {
      
      EvtPropBreitWignerRel prop(mR,gR);
      partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_KUEHN);
    }
    else if(type == "RBW_CLEO") {
      
      EvtPropBreitWignerRel prop(mR,gR);
      partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_CLEO);
    }     
    else if(type == "FLATTE") {
      
      double m1a = _dp.m( first(pairRes) );
      double m1b = _dp.m( second(pairRes) );    
      // 2nd channel
      double g2  = strtod(vv[++i].c_str(),0);
      double m2a = strtod(vv[++i].c_str(),0);
      double m2b = strtod(vv[++i].c_str(),0);
      EvtPropFlatte  prop( mR, gR, m1a, m1b, g2, m2a, m2b );
      partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::FLATTE);
    }
    else assert(0);
      
    // Optional DVFF, BVFF stanzas
    
    if(i < vv.size() - 1) {
      if(vv[i+1] == "DVFF") {	
	i++;
	if(vv[++i] == "BLATTWEISSKOPF") {
	  
	  double R = strtod(vv[++i].c_str(),0);
	  partAmp->set_fd(R);
	}
	else assert(0);
      }
    }
      
    if(i < vv.size() - 1) {
      if(vv[i+1] == "BVFF") {	
	i++;
	if(vv[++i] == "BLATTWEISSKOPF") {

	  if(_verbose) printf("BVFF=%s\n",vv[i].c_str());
	  double R = strtod(vv[++i].c_str(),0);
	  partAmp->set_fb(R);
	}
	else assert(0);
      }
    }

    const int minwidths=5;
    //Optional resonance minimum and maximum
    if(i < vv.size() - 1) {
      if(vv[i+1] == "CUTOFF") {	
	i++;
	if(vv[i+1] == "MIN") {
	  i++;
	  double min = strtod(vv[++i].c_str(),0);
	  if(_verbose) std::cout<<"CUTOFF MIN = "<<min<<" "<<minwidths<<std::endl;
	  //ensure against cutting off too close to the resonance
	  assert( min<(mR-minwidths*gR) );
	  partAmp->setmin(min);
	}
	else if (vv[i+1] == "MAX") {
	  i++;
	  double max = strtod(vv[++i].c_str(),0);
	  if(_verbose) std::cout<<"CUTOFF MAX = "<<max<<" "<<minwidths<<std::endl;
	  //ensure against cutting off too close to the resonance
	  assert( max>(mR+minwidths*gR) );
	  partAmp->setmax(max);
	}
	else assert(0);
      }
    }

    //2nd iteration in case min and max are both specified
    if(i < vv.size() - 1) {
      if(vv[i+1] == "CUTOFF") {	
	i++;
	if(vv[i+1] == "MIN") {
	  i++;
	  double min = strtod(vv[++i].c_str(),0);
	  if(_verbose) std::cout<<"CUTOFF MIN = "<<min<<std::endl;
	  //ensure against cutting off too close to the resonance
	  assert( min<(mR-minwidths*gR) );
	  partAmp->setmin(min);
	}
	else if (vv[i+1] == "MAX") {
	  i++;
	  double max = strtod(vv[++i].c_str(),0);
	  if(_verbose) std::cout<<"CUTOFF MAX = "<<max<<std::endl;
	  //ensure against cutting off too close to the resonance
	  assert( max>(mR+minwidths*gR) );
	  partAmp->setmax(max);
	}
	else assert(0);
      }
    }


    i++;
    
    pdf = new EvtDalitzResPdf(_dp,mR,gR,pairRes);
    amp = partAmp;
  }

  assert(amp);
  assert(pdf);

  if(!conj) {
    _amp->addOwnedTerm(c,amp);
  }
  else {
    _ampConj->addOwnedTerm(c,amp);
  }

  double scale = matchIsobarCoef(_amp, pdf, pairRes);
  _pc->addOwnedTerm(abs2(c)*scale,pdf);

  _names.push_back(name);
}
  
double EvtPto3PAmpFactory::matchIsobarCoef(EvtAmplitude<EvtDalitzPoint>* amp,
					   EvtPdf<EvtDalitzPoint>* pdf, 
					   EvtCyclic3::Pair ipair) {

  // account for differences in the definition of amplitudes by matching 
  //        Integral( c'*pdf ) = Integral( c*|A|^2 ) 
  // to improve generation efficiency ...

  double Ipdf  = pdf->compute_integral(10000).value();
  double Iamp2 = 0;


  EvtCyclic3::Pair jpair = EvtCyclic3::next(ipair);
  EvtCyclic3::Pair kpair = EvtCyclic3::next(jpair);

  // Trapezoidal integral
  int N=10000;
  
  double di = (_dp.qAbsMax(ipair) - _dp.qAbsMin(ipair))/((double) N);
  
  double siMin = _dp.qAbsMin(ipair);
  
  double s[3]; // playing with fire
  for(int i=1; i<N; i++) {
    
    s[ipair] = siMin + di*i;
    s[jpair] = _dp.q(jpair, 0.9999, ipair, s[ipair]);    
    s[kpair] = _dp.bigM()*_dp.bigM() - s[ipair] - s[jpair]
      + _dp.mA()*_dp.mA() + _dp.mB()*_dp.mB() + _dp.mC()*_dp.mC();
    
    EvtDalitzPoint point( _dp.mA(), _dp.mB(), _dp.mC(), 
			  s[EvtCyclic3::AB], s[EvtCyclic3::BC], s[EvtCyclic3::CA]);
    
    if (!point.isValid()) continue;
    
    double p = point.p(other(ipair), ipair);
    double q = point.p(first(ipair), ipair);
    
    double itg = abs2( amp->evaluate(point) )*di*4*q*p;
    Iamp2 += itg;
    
  }
  if (_verbose) std::cout << "integral = " << Iamp2 << "  pdf="<<Ipdf << std::endl;
  
  assert(Ipdf>0 && Iamp2>0);
  
  return Iamp2/Ipdf;
}
Commit	Line	Data
da0e9ce3	1	//-----------------------------------------------------------------------
da0e9ce3	2	// File and Version Information:
0ca57c2f	3	// $Id: EvtPto3PAmpFactory.cpp,v 1.3 2009-03-16 15:44:04 robbep Exp $
da0e9ce3	4	//
	5	// Environment:
	6	// This software is part of the EvtGen package developed jointly
	7	// for the BaBar and CLEO collaborations. If you use all or part
	8	// of it, please give an appropriate acknowledgement.
	9	//
	10	// Copyright Information:
	11	// Copyright (C) 1998 Caltech, UCSB
	12	//
	13	// Module creator:
	14	// Alexei Dvoretskii, Caltech, 2001-2002.
	15	//-----------------------------------------------------------------------
	16	#include "EvtGenBase/EvtPatches.hh"
	17
	18	// AmpFactory for building a P -> 3P decay
	19	// (pseudoscalar to three pseudoscalars)
	20
	21	#include <assert.h>
	22	#include <math.h>
	23	#include <stdio.h>
	24	#include <stdlib.h>
	25
	26	#include "EvtGenBase/EvtId.hh"
	27	#include "EvtGenBase/EvtPDL.hh"
	28	#include "EvtGenBase/EvtConst.hh"
	29	#include "EvtGenBase/EvtComplex.hh"
	30	#include "EvtGenBase/EvtCyclic3.hh"
	31	#include "EvtGenBase/EvtSpinType.hh"
	32	#include "EvtGenBase/EvtPto3PAmp.hh"
	33	#include "EvtGenBase/EvtNonresonantAmp.hh"
	34	#include "EvtGenBase/EvtFlatAmp.hh"
	35	#include "EvtGenBase/EvtLASSAmp.hh"
	36	#include "EvtGenBase/EvtPto3PAmpFactory.hh"
	37	#include "EvtGenBase/EvtPropBreitWigner.hh"
	38	#include "EvtGenBase/EvtPropFlatte.hh"
	39	#include "EvtGenBase/EvtPropBreitWignerRel.hh"
	40	#include "EvtGenBase/EvtDalitzResPdf.hh"
	41	#include "EvtGenBase/EvtDalitzFlatPdf.hh"
	42
	43	using namespace EvtCyclic3;
	44	#include <iostream>
	45
	46	void EvtPto3PAmpFactory::processAmp(EvtComplex c, std::vector<std::string> vv, bool conj)
	47	{
	48	if(_verbose) {
	49
	50	printf("Make %samplitude\n",conj ? "CP conjugate" : "");
	51	unsigned i;
	52	for(i=0;i<vv.size();i++) printf("%s\n",vv[i].c_str());
	53	printf("\n");
	54	}
	55
	56	EvtAmplitude<EvtDalitzPoint>* amp = 0;
	57	EvtPdf<EvtDalitzPoint>* pdf = 0;
	58	std::string name;
	59	Pair pairRes=AB;
	60
	61	size_t i;
	62	/*
	63	Experimental amplitudes
	64	*/
	65	if(vv[0] == "PHASESPACE") {
	66
	67	pdf = new EvtDalitzFlatPdf(_dp);
68	amp = new EvtFlatAmp<EvtDalitzPoint>();
69	name = "NR";
70	}
71	else if (!vv[0].find("NONRES")) {
72	double alpha=0;
73	EvtPto3PAmp::NumType typeNRes=EvtPto3PAmp::NONRES;
74	if (vv[0]=="NONRES_LIN") {
75	typeNRes=EvtPto3PAmp::NONRES_LIN;
76	pairRes=strToPair(vv[1].c_str());
77	}
78	else if (vv[0]=="NONRES_EXP") {
79	typeNRes=EvtPto3PAmp::NONRES_EXP;
80	pairRes = strToPair(vv[1].c_str());
81	alpha = strtod(vv[2].c_str(),0);
82	}
83	else assert(0);
84	pdf = new EvtDalitzFlatPdf(_dp);
85	amp = new EvtNonresonantAmp( &_dp, typeNRes, pairRes, alpha);
86	}
0ca57c2f	87	else if (vv[0]=="LASS" \|\| vv[0]=="LASS_ELASTIC" \|\| vv[0]=="LASS_RESONANT") {
da0e9ce3	88	pairRes = strToPair(vv[1].c_str());
	89	double m0 = strtod(vv[2].c_str(),0);
	90	double g0 = strtod(vv[3].c_str(),0);
	91	double a = strtod(vv[4].c_str(),0);
	92	double r = strtod(vv[5].c_str(),0);
	93	double cutoff = strtod(vv[6].c_str(),0);
	94	pdf = new EvtDalitzResPdf(_dp,m0,g0,pairRes);
0ca57c2f	95	amp = new EvtLASSAmp( &_dp, pairRes, m0, g0, a, r, cutoff, vv[0]);
da0e9ce3	96	}
	97
	98	/*
	99	Resonant amplitudes
	100	*/
	101	else if(vv[0] == "RESONANCE") {
	102	EvtPto3PAmp* partAmp = 0;
	103
	104	// RESONANCE stanza
	105
	106	pairRes = strToPair(vv[1].c_str());
0ca57c2f	107	EvtSpinType::spintype spinR = EvtSpinType::SCALAR;
da0e9ce3	108	double mR, gR;
	109	name = vv[2];
	110	EvtId resId = EvtPDL::getId(vv[2]);
	111	if(_verbose) printf("Particles %s form %sresonance %s\n",
	112	vv[1].c_str(),vv[2].c_str(), conj ? "(conj) " : "");
	113
	114	// If no valid particle name is given, assume that
	115	// it is the spin, the mass and the width of the particle.
	116
	117	if(resId.getId() == -1) {
	118
	119	switch(atoi(vv[2].c_str())) {
	120
	121	case 0: { spinR = EvtSpinType::SCALAR; break; }
	122	case 1: { spinR = EvtSpinType::VECTOR; break; }
	123	case 2: { spinR = EvtSpinType::TENSOR; break; }
	124	case 3: { spinR = EvtSpinType::SPIN3; break; }
	125	case 4: { spinR = EvtSpinType::SPIN4; break; }
	126	default: { assert(0); break; }
	127	}
	128
	129	mR = strtod(vv[3].c_str(),0);
	130	gR = strtod(vv[4].c_str(),0);
	131	i = 4;
	132	}
	133	else {
	134
	135	// For a valid particle get spin, mass and width
	136
	137	spinR = EvtPDL::getSpinType(resId);
	138	mR = EvtPDL::getMeanMass(resId);
	139	gR = EvtPDL::getWidth(resId);
	140	i = 2;
	141
	142	// It's possible to specify mass and width of a particle
	143	// explicitly
	144
	145	if(vv[3] != "ANGULAR") {
	146
	147	if(_verbose)
	148	printf("Setting m(%s)=%s g(%s)=%s\n",
	149	vv[2].c_str(),vv[3].c_str(),vv[2].c_str(),vv[4].c_str());
	150
	151	mR = strtod(vv[3].c_str(),0);
	152	gR = strtod(vv[4].c_str(),0);
	153	i = 4;
	154	}
	155	}
	156
	157	// ANGULAR stanza
	158
	159	if(vv[++i] != "ANGULAR") {
	160
	161	printf("%s instead of ANGULAR\n",vv[i].c_str());
	162	exit(0);
	163	}
	164	Pair pairAng = strToPair(vv[++i].c_str());
	165	if(_verbose) printf("Angle is measured between particles %s\n",vv[i].c_str());
	166
	167	// TYPE stanza
	168
0ca57c2f	169	std::string typeName = vv[++i];
0ca57c2f	170	assert(typeName == "TYPE");
da0e9ce3	171	std::string type = vv[++i];
	172	if(_verbose) printf("Propagator type %s\n",vv[i].c_str());
	173
	174	if(type == "NBW") {
	175
	176	EvtPropBreitWigner prop(mR,gR);
	177	partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::NBW);
	178	}
	179	else if(type == "RBW_ZEMACH") {
	180
	181	EvtPropBreitWignerRel prop(mR,gR);
	182	partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_ZEMACH);
	183	}
	184	else if(type == "RBW_KUEHN") {
	185
	186	EvtPropBreitWignerRel prop(mR,gR);
	187	partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_KUEHN);
	188	}
	189	else if(type == "RBW_CLEO") {
	190
	191	EvtPropBreitWignerRel prop(mR,gR);
	192	partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_CLEO);
	193	}
	194	else if(type == "FLATTE") {
	195
	196	double m1a = _dp.m( first(pairRes) );
	197	double m1b = _dp.m( second(pairRes) );
	198	// 2nd channel
	199	double g2 = strtod(vv[++i].c_str(),0);
	200	double m2a = strtod(vv[++i].c_str(),0);
	201	double m2b = strtod(vv[++i].c_str(),0);
	202	EvtPropFlatte prop( mR, gR, m1a, m1b, g2, m2a, m2b );
	203	partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::FLATTE);
	204	}
	205	else assert(0);
	206
	207	// Optional DVFF, BVFF stanzas
	208
	209	if(i < vv.size() - 1) {
	210	if(vv[i+1] == "DVFF") {
	211	i++;
	212	if(vv[++i] == "BLATTWEISSKOPF") {
	213
	214	double R = strtod(vv[++i].c_str(),0);
	215	partAmp->set_fd(R);
	216	}
	217	else assert(0);
	218	}
	219	}
	220
	221	if(i < vv.size() - 1) {
	222	if(vv[i+1] == "BVFF") {
	223	i++;
	224	if(vv[++i] == "BLATTWEISSKOPF") {
	225
	226	if(_verbose) printf("BVFF=%s\n",vv[i].c_str());
	227	double R = strtod(vv[++i].c_str(),0);
	228	partAmp->set_fb(R);
	229	}
	230	else assert(0);
	231	}
	232	}
	233
	234	const int minwidths=5;
235	//Optional resonance minimum and maximum
236	if(i < vv.size() - 1) {
237	if(vv[i+1] == "CUTOFF") {
238	i++;
239	if(vv[i+1] == "MIN") {
240	i++;
241	double min = strtod(vv[++i].c_str(),0);
0ca57c2f	242	if(_verbose) std::cout<<"CUTOFF MIN = "<<min<<" "<<minwidths<<std::endl;
da0e9ce3	243	//ensure against cutting off too close to the resonance
	244	assert( min<(mR-minwidths*gR) );
	245	partAmp->setmin(min);
	246	}
	247	else if (vv[i+1] == "MAX") {
	248	i++;
	249	double max = strtod(vv[++i].c_str(),0);
0ca57c2f	250	if(_verbose) std::cout<<"CUTOFF MAX = "<<max<<" "<<minwidths<<std::endl;
da0e9ce3	251	//ensure against cutting off too close to the resonance
	252	assert( max>(mR+minwidths*gR) );
	253	partAmp->setmax(max);
	254	}
	255	else assert(0);
	256	}
	257	}
	258
	259	//2nd iteration in case min and max are both specified
	260	if(i < vv.size() - 1) {
	261	if(vv[i+1] == "CUTOFF") {
	262	i++;
	263	if(vv[i+1] == "MIN") {
	264	i++;
	265	double min = strtod(vv[++i].c_str(),0);
	266	if(_verbose) std::cout<<"CUTOFF MIN = "<<min<<std::endl;
	267	//ensure against cutting off too close to the resonance
	268	assert( min<(mR-minwidths*gR) );
	269	partAmp->setmin(min);
	270	}
	271	else if (vv[i+1] == "MAX") {
	272	i++;
	273	double max = strtod(vv[++i].c_str(),0);
	274	if(_verbose) std::cout<<"CUTOFF MAX = "<<max<<std::endl;
	275	//ensure against cutting off too close to the resonance
	276	assert( max>(mR+minwidths*gR) );
	277	partAmp->setmax(max);
	278	}
	279	else assert(0);
	280	}
	281	}
	282
	283
	284	i++;
	285
	286	pdf = new EvtDalitzResPdf(_dp,mR,gR,pairRes);
	287	amp = partAmp;
	288	}
	289
	290	assert(amp);
	291	assert(pdf);
	292
	293	if(!conj) {
	294	_amp->addOwnedTerm(c,amp);
	295	}
	296	else {
	297	_ampConj->addOwnedTerm(c,amp);
	298	}
	299
	300	double scale = matchIsobarCoef(_amp, pdf, pairRes);
	301	_pc->addOwnedTerm(abs2(c)*scale,pdf);
	302
	303	_names.push_back(name);
	304	}
	305
	306	double EvtPto3PAmpFactory::matchIsobarCoef(EvtAmplitude<EvtDalitzPoint>* amp,
	307	EvtPdf<EvtDalitzPoint>* pdf,
	308	EvtCyclic3::Pair ipair) {
	309
	310	// account for differences in the definition of amplitudes by matching
	311	// Integral( c'pdf ) = Integral( c\|A\|^2 )
	312	// to improve generation efficiency ...
	313
	314	double Ipdf = pdf->compute_integral(10000).value();
315	double Iamp2 = 0;
316
317
318	EvtCyclic3::Pair jpair = EvtCyclic3::next(ipair);
319	EvtCyclic3::Pair kpair = EvtCyclic3::next(jpair);
320
321	// Trapezoidal integral
322	int N=10000;
323
324	double di = (_dp.qAbsMax(ipair) - _dp.qAbsMin(ipair))/((double) N);
325
326	double siMin = _dp.qAbsMin(ipair);
327
328	double s[3]; // playing with fire
329	for(int i=1; i<N; i++) {
330
331	s[ipair] = siMin + di*i;
332	s[jpair] = _dp.q(jpair, 0.9999, ipair, s[ipair]);
333	s[kpair] = _dp.bigM()*_dp.bigM() - s[ipair] - s[jpair]
334	+ _dp.mA()_dp.mA() + _dp.mB()_dp.mB() + _dp.mC()*_dp.mC();
335
336	EvtDalitzPoint point( _dp.mA(), _dp.mB(), _dp.mC(),
337	s[EvtCyclic3::AB], s[EvtCyclic3::BC], s[EvtCyclic3::CA]);
338
339	if (!point.isValid()) continue;
340
341	double p = point.p(other(ipair), ipair);
342	double q = point.p(first(ipair), ipair);
343
344	double itg = abs2( amp->evaluate(point) )di4qp;
345	Iamp2 += itg;
346
347	}
348	if (_verbose) std::cout << "integral = " << Iamp2 << " pdf="<<Ipdf << std::endl;
349
350	assert(Ipdf>0 && Iamp2>0);
351
352	return Iamp2/Ipdf;
353	}