[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtRelBreitWignerBarrierFact.cxx

//--------------------------------------------------------------------------
//
// 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: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtLineShape.cc
//
// Description: Store particle properties for one particle.
//
// Modification history:
//
//    Lange      March 10, 2001        Module created
//    Dvoretskii June  03, 2002        Reimplemented rollMass()
//
//------------------------------------------------------------------------
#include "EvtGenBase/EvtPatches.hh"

#include "EvtGenBase/EvtPatches.hh"

#include "EvtGenBase/EvtPredGen.hh"

#include "EvtGenBase/EvtRelBreitWignerBarrierFact.hh"
#include "EvtGenBase/EvtTwoBodyVertex.hh"
#include "EvtGenBase/EvtPropBreitWignerRel.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtAmpPdf.hh"
#include "EvtGenBase/EvtMassAmp.hh"
#include "EvtGenBase/EvtSpinType.hh"
#include "EvtGenBase/EvtIntervalFlatPdf.hh"
#include <algorithm>
EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact() {

}

EvtRelBreitWignerBarrierFact::~EvtRelBreitWignerBarrierFact() {
}

EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(double mass, double width, double maxRange, EvtSpinType::spintype sp) :
  EvtAbsLineShape(mass,width,maxRange,sp)
{ // double mDaug1, double mDaug2, int l) {

  _includeDecayFact=true;
  _includeBirthFact=true;
  _mass=mass;
  _width=width;
  _spin=sp;
  _blatt=3.0;
  _maxRange=maxRange;
  _errorCond=false;

  double maxdelta = 15.0*width;

  if ( maxRange > 0.00001 ) {
    _massMax=mass+maxdelta;
    _massMin=mass-maxRange;
  }
  else{
    _massMax=mass+maxdelta;
    _massMin=mass-15.0*width;
  }
 
  _massMax=mass+maxdelta;
  if ( _massMin< 0. ) _massMin=0.;
}

EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(const EvtRelBreitWignerBarrierFact& x) :
  EvtAbsLineShape(x)
{
  _massMax=x._massMax;
  _massMin=x._massMin;
  _blatt=x._blatt;
  _maxRange=x._maxRange;
  _includeDecayFact=x._includeDecayFact;
  _includeBirthFact=x._includeBirthFact;
  _errorCond=x._errorCond;

}

EvtRelBreitWignerBarrierFact& EvtRelBreitWignerBarrierFact::operator=(const EvtRelBreitWignerBarrierFact& x) {
  _mass=x._mass;
  _width=x._width;
  _spin=x._spin;
  _massMax=x._massMax;
  _massMin=x._massMin;
  _blatt=x._blatt;
  _maxRange=x._maxRange;
  _includeDecayFact=x._includeDecayFact;
  _includeBirthFact=x._includeBirthFact;
  _errorCond=x._errorCond;
  
  return *this;
}

EvtAbsLineShape* EvtRelBreitWignerBarrierFact::clone() {

  return new EvtRelBreitWignerBarrierFact(*this);
}


double EvtRelBreitWignerBarrierFact::getMassProb(double mass, double massPar,int nDaug, double *massDau) {

  _errorCond=false;
  //return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);
  if (nDaug!=2) return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);

  double dTotMass=0.;

  int i;
  for (i=0; i<nDaug; i++) {
    dTotMass+=massDau[i];
  }
  //report(INFO,"EvtGen") << mass << " " << massPar << " " << dTotMass << " "<< endl;
  //    if ( (mass-dTotMass)<0.0001 ) return 0.;
  //report(INFO,"EvtGen") << mass << " " << dTotMass << endl;
  if ( (mass<dTotMass) ) return 0.;

  if ( _width< 0.0001) return 1.;

  if ( massPar>0.0000000001 ) {
    if ( mass > massPar) return 0.;
  }

  if ( _errorCond ) return 0.;

  // we did all the work in getRandMass
  return 1.;
}

double EvtRelBreitWignerBarrierFact::getRandMass(EvtId *parId,int nDaug, EvtId *dauId, EvtId *othDaugId, double maxMass, double *dauMasses) {
  if ( nDaug!=2) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);

  if ( _width< 0.00001) return _mass;

  //first figure out L - take the lowest allowed.

  EvtSpinType::spintype spinD1=EvtPDL::getSpinType(dauId[0]);
  EvtSpinType::spintype spinD2=EvtPDL::getSpinType(dauId[1]);

  int t1=EvtSpinType::getSpin2(spinD1);
  int t2=EvtSpinType::getSpin2(spinD2);
  int t3=EvtSpinType::getSpin2(_spin);

  int Lmin=-10;


  // the user has overridden the partial wave to use.
  for (unsigned int vC=0; vC<_userSetPW.size(); vC++) {
    if ( dauId[0]==_userSetPWD1[vC] &&  dauId[1]==_userSetPWD2[vC] ) {
      Lmin=2*_userSetPW[vC]; 
    }
    if ( dauId[0]==_userSetPWD2[vC] &&  dauId[1]==_userSetPWD1[vC] ) {
      Lmin=2*_userSetPW[vC]; 
    }
  }
  
  // allow for special cases.
  if (Lmin<-1 ) {
    
    //There are some things I don't know how to deal with
    if ( t3>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
    if ( t1>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
    if ( t2>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
    
    //figure the min and max allowwed "spins" for the daughters state
    Lmin=std::max(t3-t2-t1,std::max(t2-t3-t1,t1-t3-t2));
    if (Lmin<0) Lmin=0;
    assert(Lmin==0||Lmin==2||Lmin==4);
  }

  //double massD1=EvtPDL::getMeanMass(dauId[0]);
  //double massD2=EvtPDL::getMeanMass(dauId[1]);
  double massD1=dauMasses[0];
  double massD2=dauMasses[1];

  // I'm not sure how to define the vertex factor here - so retreat to nonRel code.
  if ( (massD1+massD2)> _mass ) return  EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);

  //parent vertex factor not yet implemented
  double massOthD=-10.;
  double massParent=-10.;
  int birthl=-10;
  if ( othDaugId) {
    EvtSpinType::spintype spinOth=EvtPDL::getSpinType(*othDaugId);
    EvtSpinType::spintype spinPar=EvtPDL::getSpinType(*parId);
    
    int tt1=EvtSpinType::getSpin2(spinOth);
    int tt2=EvtSpinType::getSpin2(spinPar);
    int tt3=EvtSpinType::getSpin2(_spin);
    
    
    //figure the min and max allowwed "spins" for the daughters state
    if ( (tt1<=4) && ( tt2<=4) ) {
      birthl=std::max(tt3-tt2-tt1,std::max(tt2-tt3-tt1,tt1-tt3-tt2));
      if (birthl<0) birthl=0;
    
      massOthD=EvtPDL::getMeanMass(*othDaugId);
      massParent=EvtPDL::getMeanMass(*parId);
    
    }

    // allow user to override
    for (size_t vC=0; vC<_userSetBirthPW.size(); vC++) {
      if ( *othDaugId==_userSetBirthOthD[vC] && *parId==_userSetBirthPar[vC]){
	birthl=2*_userSetBirthPW[vC];
      } 
    }

  }
  double massM=_massMax;
  if ( (maxMass > -0.5) && (maxMass < massM) ) massM=maxMass;

  //special case... if the parent mass is _fixed_ we can do a little better
  //and only for a two body decay as that seems to be where we have problems

  // Define relativistic propagator amplitude

  EvtTwoBodyVertex vd(massD1,massD2,_mass,Lmin/2);
  vd.set_f(_blatt);
  EvtPropBreitWignerRel bw(_mass,_width);
  EvtMassAmp amp(bw,vd);

  if ( _includeDecayFact) {
    amp.addDeathFact();
    amp.addDeathFactFF();
  }
  if ( massParent>-1.) {
    if ( _includeBirthFact ) {

      EvtTwoBodyVertex vb(_mass,massOthD,massParent,birthl/2);
      //whoops 060116
      //needed to actually tell the vertex about the form factor!
      //otherwise its just 1
      if ( _applyFixForSP8 ) vb.set_f(_blatt);
      amp.setBirthVtx(vb);
      amp.addBirthFact();
      amp.addBirthFactFF();
    }
  }


  EvtAmpPdf<EvtPoint1D> pdf(amp);

  // Estimate maximum and create predicate for accept reject


  double tempMaxLoc=_mass;
  if ( maxMass>-0.5 && maxMass<_mass) tempMaxLoc=maxMass;
  double tempMax=_massMax;
  if ( maxMass>-0.5 && maxMass<_massMax) tempMax=maxMass;
  double tempMinMass=_massMin;
  if ( massD1+massD2 > _massMin) tempMinMass=massD1+massD2;

  //redo sanity check - is there a solution to our problem.
  //if not return an error condition that is caught by the
  //mass prob calculation above.
  if ( tempMinMass > tempMax ) {
    _errorCond=true;
    return tempMinMass;
  }
  
  if ( tempMaxLoc < tempMinMass) tempMaxLoc=tempMinMass;

  double safetyFactor=1.2;
  if ( _applyFixForSP8 ) safetyFactor=1.4;

  EvtPdfMax<EvtPoint1D> max(safetyFactor*pdf.evaluate(EvtPoint1D(tempMinMass,tempMax,tempMaxLoc)));

  EvtPdfPred<EvtPoint1D> pred(pdf);
  pred.setMax(max);

  EvtIntervalFlatPdf flat(tempMinMass,tempMax);
  EvtPdfGen<EvtPoint1D> gen(flat);
  EvtPredGen<EvtPdfGen<EvtPoint1D>,EvtPdfPred<EvtPoint1D> > predgen(gen,pred);

  EvtPoint1D point = predgen();
  return point.value();

}
Commit	Line	Data
da0e9ce3	1	//--------------------------------------------------------------------------
	2	//
	3	// Environment:
	4	// This software is part of the EvtGen package developed jointly
	5	// for the BaBar and CLEO collaborations. If you use all or part
	6	// of it, please give an appropriate acknowledgement.
	7	//
	8	// Copyright Information: See EvtGen/COPYRIGHT
	9	// Copyright (C) 1998 Caltech, UCSB
	10	//
	11	// Module: EvtLineShape.cc
	12	//
	13	// Description: Store particle properties for one particle.
	14	//
	15	// Modification history:
	16	//
	17	// Lange March 10, 2001 Module created
	18	// Dvoretskii June 03, 2002 Reimplemented rollMass()
	19	//
	20	//------------------------------------------------------------------------
	21	#include "EvtGenBase/EvtPatches.hh"
	22
	23	#include "EvtGenBase/EvtPatches.hh"
	24
	25	#include "EvtGenBase/EvtPredGen.hh"
	26
	27	#include "EvtGenBase/EvtRelBreitWignerBarrierFact.hh"
	28	#include "EvtGenBase/EvtTwoBodyVertex.hh"
	29	#include "EvtGenBase/EvtPropBreitWignerRel.hh"
	30	#include "EvtGenBase/EvtPDL.hh"
	31	#include "EvtGenBase/EvtAmpPdf.hh"
	32	#include "EvtGenBase/EvtMassAmp.hh"
	33	#include "EvtGenBase/EvtSpinType.hh"
	34	#include "EvtGenBase/EvtIntervalFlatPdf.hh"
	35	#include <algorithm>
	36	EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact() {
	37
	38	}
	39
	40	EvtRelBreitWignerBarrierFact::~EvtRelBreitWignerBarrierFact() {
	41	}
	42
	43	EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(double mass, double width, double maxRange, EvtSpinType::spintype sp) :
	44	EvtAbsLineShape(mass,width,maxRange,sp)
	45	{ // double mDaug1, double mDaug2, int l) {
	46
	47	_includeDecayFact=true;
	48	_includeBirthFact=true;
	49	_mass=mass;
	50	_width=width;
	51	_spin=sp;
	52	_blatt=3.0;
	53	_maxRange=maxRange;
	54	_errorCond=false;
	55
	56	double maxdelta = 15.0*width;
	57
	58	if ( maxRange > 0.00001 ) {
	59	_massMax=mass+maxdelta;
	60	_massMin=mass-maxRange;
	61	}
	62	else{
	63	_massMax=mass+maxdelta;
	64	_massMin=mass-15.0*width;
65	}
66
67	_massMax=mass+maxdelta;
68	if ( _massMin< 0. ) _massMin=0.;
69	}
70
71	EvtRelBreitWignerBarrierFact::EvtRelBreitWignerBarrierFact(const EvtRelBreitWignerBarrierFact& x) :
72	EvtAbsLineShape(x)
73	{
74	_massMax=x._massMax;
75	_massMin=x._massMin;
76	_blatt=x._blatt;
77	_maxRange=x._maxRange;
78	_includeDecayFact=x._includeDecayFact;
79	_includeBirthFact=x._includeBirthFact;
80	_errorCond=x._errorCond;
81
82	}
83
84	EvtRelBreitWignerBarrierFact& EvtRelBreitWignerBarrierFact::operator=(const EvtRelBreitWignerBarrierFact& x) {
85	_mass=x._mass;
86	_width=x._width;
87	_spin=x._spin;
88	_massMax=x._massMax;
89	_massMin=x._massMin;
90	_blatt=x._blatt;
91	_maxRange=x._maxRange;
92	_includeDecayFact=x._includeDecayFact;
93	_includeBirthFact=x._includeBirthFact;
94	_errorCond=x._errorCond;
95
96	return *this;
97	}
98
99	EvtAbsLineShape* EvtRelBreitWignerBarrierFact::clone() {
100
101	return new EvtRelBreitWignerBarrierFact(*this);
102	}
103
104
105	double EvtRelBreitWignerBarrierFact::getMassProb(double mass, double massPar,int nDaug, double *massDau) {
106
107	_errorCond=false;
108	//return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);
109	if (nDaug!=2) return EvtAbsLineShape::getMassProb(mass,massPar,nDaug,massDau);
110
111	double dTotMass=0.;
112
113	int i;
114	for (i=0; i<nDaug; i++) {
115	dTotMass+=massDau[i];
116	}
117	//report(INFO,"EvtGen") << mass << " " << massPar << " " << dTotMass << " "<< endl;
118	// if ( (mass-dTotMass)<0.0001 ) return 0.;
119	//report(INFO,"EvtGen") << mass << " " << dTotMass << endl;
120	if ( (mass<dTotMass) ) return 0.;
121
122	if ( _width< 0.0001) return 1.;
123
124	if ( massPar>0.0000000001 ) {
125	if ( mass > massPar) return 0.;
126	}
127
128	if ( _errorCond ) return 0.;
129
130	// we did all the work in getRandMass
131	return 1.;
132	}
133
134	double EvtRelBreitWignerBarrierFact::getRandMass(EvtId parId,int nDaug, EvtId dauId, EvtId othDaugId, double maxMass, double dauMasses) {
135	if ( nDaug!=2) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
136
137	if ( _width< 0.00001) return _mass;
138
139	//first figure out L - take the lowest allowed.
140
141	EvtSpinType::spintype spinD1=EvtPDL::getSpinType(dauId[0]);
142	EvtSpinType::spintype spinD2=EvtPDL::getSpinType(dauId[1]);
143
144	int t1=EvtSpinType::getSpin2(spinD1);
145	int t2=EvtSpinType::getSpin2(spinD2);
146	int t3=EvtSpinType::getSpin2(_spin);
147
148	int Lmin=-10;
149
150
151	// the user has overridden the partial wave to use.
152	for (unsigned int vC=0; vC<_userSetPW.size(); vC++) {
153	if ( dauId[0]==_userSetPWD1[vC] && dauId[1]==_userSetPWD2[vC] ) {
154	Lmin=2*_userSetPW[vC];
155	}
156	if ( dauId[0]==_userSetPWD2[vC] && dauId[1]==_userSetPWD1[vC] ) {
157	Lmin=2*_userSetPW[vC];
158	}
159	}
160
161	// allow for special cases.
162	if (Lmin<-1 ) {
163
164	//There are some things I don't know how to deal with
165	if ( t3>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
166	if ( t1>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
167	if ( t2>4) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
168
169	//figure the min and max allowwed "spins" for the daughters state
170	Lmin=std::max(t3-t2-t1,std::max(t2-t3-t1,t1-t3-t2));
171	if (Lmin<0) Lmin=0;
172	assert(Lmin==0\|\|Lmin==2\|\|Lmin==4);
173	}
174
175	//double massD1=EvtPDL::getMeanMass(dauId[0]);
176	//double massD2=EvtPDL::getMeanMass(dauId[1]);
177	double massD1=dauMasses[0];
178	double massD2=dauMasses[1];
179
180	// I'm not sure how to define the vertex factor here - so retreat to nonRel code.
181	if ( (massD1+massD2)> _mass ) return EvtAbsLineShape::getRandMass(parId,nDaug,dauId,othDaugId,maxMass,dauMasses);
182
183	//parent vertex factor not yet implemented
184	double massOthD=-10.;
185	double massParent=-10.;
186	int birthl=-10;
187	if ( othDaugId) {
188	EvtSpinType::spintype spinOth=EvtPDL::getSpinType(*othDaugId);
189	EvtSpinType::spintype spinPar=EvtPDL::getSpinType(*parId);
190
191	int tt1=EvtSpinType::getSpin2(spinOth);
192	int tt2=EvtSpinType::getSpin2(spinPar);
193	int tt3=EvtSpinType::getSpin2(_spin);
194
195
196	//figure the min and max allowwed "spins" for the daughters state
197	if ( (tt1<=4) && ( tt2<=4) ) {
198	birthl=std::max(tt3-tt2-tt1,std::max(tt2-tt3-tt1,tt1-tt3-tt2));
199	if (birthl<0) birthl=0;
200
201	massOthD=EvtPDL::getMeanMass(*othDaugId);
202	massParent=EvtPDL::getMeanMass(*parId);
203
204	}
205
206	// allow user to override
207	for (size_t vC=0; vC<_userSetBirthPW.size(); vC++) {
208	if ( othDaugId==_userSetBirthOthD[vC] && parId==_userSetBirthPar[vC]){
209	birthl=2*_userSetBirthPW[vC];
210	}
211	}
212
213	}
214	double massM=_massMax;
215	if ( (maxMass > -0.5) && (maxMass < massM) ) massM=maxMass;
216
217	//special case... if the parent mass is _fixed_ we can do a little better
218	//and only for a two body decay as that seems to be where we have problems
219
220	// Define relativistic propagator amplitude
221
222	EvtTwoBodyVertex vd(massD1,massD2,_mass,Lmin/2);
223	vd.set_f(_blatt);
224	EvtPropBreitWignerRel bw(_mass,_width);
225	EvtMassAmp amp(bw,vd);
226
227	if ( _includeDecayFact) {
228	amp.addDeathFact();
229	amp.addDeathFactFF();
230	}
231	if ( massParent>-1.) {
232	if ( _includeBirthFact ) {
233
234	EvtTwoBodyVertex vb(_mass,massOthD,massParent,birthl/2);
235	//whoops 060116
236	//needed to actually tell the vertex about the form factor!
237	//otherwise its just 1
238	if ( _applyFixForSP8 ) vb.set_f(_blatt);
239	amp.setBirthVtx(vb);
240	amp.addBirthFact();
241	amp.addBirthFactFF();
242	}
243	}
244
245
246	EvtAmpPdf<EvtPoint1D> pdf(amp);
247
248	// Estimate maximum and create predicate for accept reject
249
250
251	double tempMaxLoc=_mass;
252	if ( maxMass>-0.5 && maxMass<_mass) tempMaxLoc=maxMass;
253	double tempMax=_massMax;
254	if ( maxMass>-0.5 && maxMass<_massMax) tempMax=maxMass;
255	double tempMinMass=_massMin;
256	if ( massD1+massD2 > _massMin) tempMinMass=massD1+massD2;
257
258	//redo sanity check - is there a solution to our problem.
259	//if not return an error condition that is caught by the
260	//mass prob calculation above.
261	if ( tempMinMass > tempMax ) {
262	_errorCond=true;
263	return tempMinMass;
264	}
265
266	if ( tempMaxLoc < tempMinMass) tempMaxLoc=tempMinMass;
267
268	double safetyFactor=1.2;
269	if ( _applyFixForSP8 ) safetyFactor=1.4;
270
271	EvtPdfMax<EvtPoint1D> max(safetyFactor*pdf.evaluate(EvtPoint1D(tempMinMass,tempMax,tempMaxLoc)));
272
273	EvtPdfPred<EvtPoint1D> pred(pdf);
274	pred.setMax(max);
275
276	EvtIntervalFlatPdf flat(tempMinMass,tempMax);
277	EvtPdfGen<EvtPoint1D> gen(flat);
278	EvtPredGen<EvtPdfGen<EvtPoint1D>,EvtPdfPred<EvtPoint1D> > predgen(gen,pred);
279
280	EvtPoint1D point = predgen();
281	return point.value();
282
283	}
284
285
286
287
288
289
290
291
292