[u/mrichter/AliRoot.git] / TEvtGen / EvtGenBase / EvtGenKine.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: EvtGenKine.cc
//
// Description: Tools for generating distributions of four vectors in
//              phasespace
//
// Modification history:
//
//    DJL/RYD     September 25, 1996         Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <iostream>
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtConst.hh"
#include <math.h>
using std::endl;


double EvtPawt(double a,double b,double c)
{
  double temp=(a*a-(b+c)*(b+c))*(a*a-(b-c)*(b-c));

  if (temp<=0) {
     return 0.0;
  }

  return sqrt(temp)/(2.0*a);
}


double EvtGenKine::PhaseSpace( int ndaug, double mass[30], EvtVector4R p4[30], 
			       double mp )

//  N body phase space routine.  Send parent with
//  daughters already defined ( Number and masses )
//  Returns four vectors in parent frame.

{

  double energy, p3, alpha, beta;

  if ( ndaug == 1 ) {
     p4[0].set(mass[0],0.0,0.0,0.0);
     return 1.0;
  }


  if ( ndaug == 2 ) {

     //Two body phase space

     energy = ( mp*mp + mass[0]*mass[0] -
              mass[1]*mass[1] ) / ( 2.0 * mp );

     p3 = sqrt( energy*energy - mass[0]*mass[0] );

     p4[0].set( energy, 0.0, 0.0, p3 );

     energy = mp - energy;
     p3 = -1.0*p3;
     p4[1].set( energy, 0.0, 0.0, p3 );

     //Now rotate four vectors.

     alpha = EvtRandom::Flat( EvtConst::twoPi );
     beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
   
     p4[0].applyRotateEuler( alpha, beta, -alpha );
     p4[1].applyRotateEuler( alpha, beta, -alpha );

     return 1.0;
  }

  if ( ndaug != 2 ) {

    double wtmax=0.0;
    double pm[5][30],to[4],pmin,pmax,psum,rnd[30];
    double ran,wt,pa,costh,sinth,phi,p[4][30],be[4],bep,temp;
    int i,il,ilr,i1,il1u,il1,il2r,ilu;
    int il2=0;
    
     for(i=0;i<ndaug;i++){
       pm[4][i]=0.0;
       rnd[i]=0.0;
     }     
    
     pm[0][0]=mp;
     pm[1][0]=0.0;
     pm[2][0]=0.0;
     pm[3][0]=0.0;
     pm[4][0]=mp;

     to[0]=mp;
     to[1]=0.0;
     to[2]=0.0;
     to[3]=0.0;

     psum=0.0;
     for(i=1;i<ndaug+1;i++){
       psum=psum+mass[i-1];
     }

     pm[4][ndaug-1]=mass[ndaug-1];

     switch (ndaug) {

     case 1:
       wtmax=1.0/16.0;
       break;
     case 2:
       wtmax=1.0/150.0;
       break;
     case 3:
       wtmax=1.0/2.0;
       break;
     case 4:
       wtmax=1.0/5.0;
       break;
     case 5:
       wtmax=1.0/15.0;
       break;
     case 6:
       wtmax=1.0/15.0;
       break;
     case 7:
       wtmax=1.0/15.0;
       break;
     case 8:
       wtmax=1.0/15.0;
       break;
     case 9:
       wtmax=1.0/15.0;
       break;
     case 10:
       wtmax=1.0/15.0;
       break;
     case 11:
       wtmax=1.0/15.0;
       break;
     case 12:
       wtmax=1.0/15.0;
       break;
     case 13:
       wtmax=1.0/15.0;
       break;
     case 14:
       wtmax=1.0/15.0;
       break;
     case 15:
       wtmax=1.0/15.0;
       break;
     default:
       report(ERROR,"EvtGen") << "too many daughters for phase space..." << ndaug << " "<< mp <<endl;;
       break;
     }

     pmax=mp-psum+mass[ndaug-1];

     pmin=0.0;

     for(ilr=2;ilr<ndaug+1;ilr++){
       il=ndaug+1-ilr;
       pmax=pmax+mass[il-1];
       pmin=pmin+mass[il+1-1];
       wtmax=wtmax*EvtPawt(pmax,pmin,mass[il-1]);
     }

     do{

       rnd[0]=1.0;
       il1u=ndaug-1;

       for (il1=2;il1<il1u+1;il1++){
	 ran=EvtRandom::Flat();
	 for (il2r=2;il2r<il1+1;il2r++){
	   il2=il1+1-il2r;
	   if (ran<=rnd[il2-1]) goto two39;
	   rnd[il2+1-1]=rnd[il2-1];
	 }
       two39:
	 rnd[il2+1-1]=ran;
       }

       rnd[ndaug-1]=0.0;
       wt=1.0;
       for(ilr=2;ilr<ndaug+1;ilr++){
	 il=ndaug+1-ilr;
	 pm[4][il-1]=pm[4][il+1-1]+mass[il-1]+
           (rnd[il-1]-rnd[il+1-1])*(mp-psum);
	 wt=wt*EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
       }
       if (wt>wtmax) {
	 report(ERROR,"EvtGen") << "wtmax to small in EvtPhaseSpace with "
				<< ndaug <<" daughters"<<endl;;
       }
     } while (wt<EvtRandom::Flat(wtmax));
     
     ilu=ndaug-1;

     for (il=1;il<ilu+1;il++){
       pa=EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
       costh=EvtRandom::Flat(-1.0,1.0);
       sinth=sqrt(1.0-costh*costh);
       phi=EvtRandom::Flat(EvtConst::twoPi);
       p[1][il-1]=pa*sinth*cos(phi);
       p[2][il-1]=pa*sinth*sin(phi);
       p[3][il-1]=pa*costh;
       pm[1][il+1-1]=-p[1][il-1];
       pm[2][il+1-1]=-p[2][il-1];
       pm[3][il+1-1]=-p[3][il-1];
       p[0][il-1]=sqrt(pa*pa+mass[il-1]*mass[il-1]);
       pm[0][il+1-1]=sqrt(pa*pa+pm[4][il+1-1]*pm[4][il+1-1]);
     }

      p[0][ndaug-1]=pm[0][ndaug-1];
      p[1][ndaug-1]=pm[1][ndaug-1];
      p[2][ndaug-1]=pm[2][ndaug-1];
      p[3][ndaug-1]=pm[3][ndaug-1];

      for (ilr=2;ilr<ndaug+1;ilr++){
        il=ndaug+1-ilr;
        be[0]=pm[0][il-1]/pm[4][il-1];
        be[1]=pm[1][il-1]/pm[4][il-1];
        be[2]=pm[2][il-1]/pm[4][il-1];
        be[3]=pm[3][il-1]/pm[4][il-1];

        for (i1=il;i1<ndaug+1;i1++){
          bep=be[1]*p[1][i1-1]+be[2]*p[2][i1-1]+
              be[3]*p[3][i1-1]+be[0]*p[0][i1-1];
          temp=(p[0][i1-1]+bep)/(be[0]+1.0);
          p[1][i1-1]=p[1][i1-1]+temp*be[1];
          p[2][i1-1]=p[2][i1-1]+temp*be[2];
          p[3][i1-1]=p[3][i1-1]+temp*be[3];
          p[0][i1-1]=bep;

        }
      }

     for (ilr=0;ilr<ndaug;ilr++){
       p4[ilr].set(p[0][ilr],p[1][ilr],p[2][ilr],p[3][ilr]);
     }

     return 1.0;
  }

  return 1.0;

}


double EvtGenKine::PhaseSpacePole(double M, double m1, double m2, double m3, 
				  double a,EvtVector4R p4[10])

//  generate kinematics for 3 body decays, pole for the m1,m2 mass.

{

  //f1   = 1  (phasespace)
  //f2   = a*(1/m12sq)^2

  double m12sqmax=(M-m3)*(M-m3);
  double m12sqmin=(m1+m2)*(m1+m2);

  double m13sqmax=(M-m2)*(M-m2);
  double m13sqmin=(m1+m3)*(m1+m3);

  double v1=(m12sqmax-m12sqmin)*(m13sqmax-m13sqmin);
  double v2= a*(1.0/m12sqmin-1.0/m12sqmax)*(m13sqmax-m13sqmin);

  double m12sq,m13sq;

  double r=v1/(v1+v2);

  double m13min,m13max;

  do{

    m13sq=EvtRandom::Flat(m13sqmin,m13sqmax);
  
    if (r>EvtRandom::Flat()){
      m12sq=EvtRandom::Flat(m12sqmin,m12sqmax);
    }
    else{
      m12sq=1.0/(1.0/m12sqmin-EvtRandom::Flat()*(1.0/m12sqmin-1.0/m12sqmax));
    }

    //kinematically allowed?
    double E3star=(M*M-m12sq-m3*m3)/sqrt(4*m12sq);
    double E1star=(m12sq+m1*m1-m2*m2)/sqrt(4*m12sq);
    double p3star=sqrt(E3star*E3star-m3*m3);
    double p1star=sqrt(E1star*E1star-m1*m1);
    m13max=(E3star+E1star)*(E3star+E1star)-
      (p3star-p1star)*(p3star-p1star);
    m13min=(E3star+E1star)*(E3star+E1star)-
      (p3star+p1star)*(p3star+p1star);

  }while(m13sq<m13min||m13sq>m13max);

  
  double E2=(M*M+m2*m2-m13sq)/(2.0*M);
  double E3=(M*M+m3*m3-m12sq)/(2.0*M);
  double E1=M-E2-E3;
  double p1mom=sqrt(E1*E1-m1*m1);      
  double p3mom=sqrt(E3*E3-m3*m3);
  double cost13=(2.0*E1*E3+m1*m1+m3*m3-m13sq)/(2.0*p1mom*p3mom);

  //report(INFO,"EvtGen") << m13sq << endl;
  //report(INFO,"EvtGen") << m12sq << endl;
  //report(INFO,"EvtGen") << E1 << endl;
  //report(INFO,"EvtGen") << E2 << endl;
  //report(INFO,"EvtGen") << E3 << endl;
  //report(INFO,"EvtGen") << p1mom << endl;
  //report(INFO,"EvtGen") << p3mom << endl;
  //report(INFO,"EvtGen") << cost13 << endl;
  

  p4[2].set(E3,0.0,0.0,p3mom);
  p4[0].set(E1,p1mom*sqrt(1.0-cost13*cost13),0.0,p1mom*cost13);
  p4[1].set(E2,-p1mom*sqrt(1.0-cost13*cost13),0.0,-p1mom*cost13-p3mom);

  //report(INFO,"EvtGen") << "p4:"<<p4[0]<<p4[1]<<p4[2]<<endl;

  double alpha = EvtRandom::Flat( EvtConst::twoPi );
  double beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
  double gamma = EvtRandom::Flat( EvtConst::twoPi );
  
  p4[0].applyRotateEuler( alpha, beta, gamma );
  p4[1].applyRotateEuler( alpha, beta, gamma );
  p4[2].applyRotateEuler( alpha, beta, gamma );
  
  return 1.0+a/(m12sq*m12sq);

}
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: EvtGenKine.cc
	12	//
	13	// Description: Tools for generating distributions of four vectors in
	14	// phasespace
	15	//
	16	// Modification history:
	17	//
	18	// DJL/RYD September 25, 1996 Module created
	19	//
	20	//------------------------------------------------------------------------
	21	//
	22	#include "EvtGenBase/EvtPatches.hh"
	23	#include <iostream>
	24	#include "EvtGenBase/EvtGenKine.hh"
	25	#include "EvtGenBase/EvtRandom.hh"
	26	#include "EvtGenBase/EvtVector4R.hh"
	27	#include "EvtGenBase/EvtReport.hh"
	28	#include "EvtGenBase/EvtConst.hh"
	29	#include <math.h>
	30	using std::endl;
	31
	32
	33	double EvtPawt(double a,double b,double c)
	34	{
	35	double temp=(aa-(b+c)(b+c))(aa-(b-c)*(b-c));
	36
	37	if (temp<=0) {
	38	return 0.0;
	39	}
	40
	41	return sqrt(temp)/(2.0*a);
	42	}
	43
	44
	45	double EvtGenKine::PhaseSpace( int ndaug, double mass[30], EvtVector4R p4[30],
	46	double mp )
	47
	48	// N body phase space routine. Send parent with
	49	// daughters already defined ( Number and masses )
	50	// Returns four vectors in parent frame.
	51
	52	{
	53
	54	double energy, p3, alpha, beta;
	55
	56	if ( ndaug == 1 ) {
	57	p4[0].set(mass[0],0.0,0.0,0.0);
	58	return 1.0;
	59	}
	60
	61
	62	if ( ndaug == 2 ) {
	63
	64	//Two body phase space
65
66	energy = ( mpmp + mass[0]mass[0] -
67	mass[1]mass[1] ) / ( 2.0 mp );
68
69	p3 = sqrt( energyenergy - mass[0]mass[0] );
70
71	p4[0].set( energy, 0.0, 0.0, p3 );
72
73	energy = mp - energy;
74	p3 = -1.0*p3;
75	p4[1].set( energy, 0.0, 0.0, p3 );
76
77	//Now rotate four vectors.
78
79	alpha = EvtRandom::Flat( EvtConst::twoPi );
80	beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
81
82	p4[0].applyRotateEuler( alpha, beta, -alpha );
83	p4[1].applyRotateEuler( alpha, beta, -alpha );
84
85	return 1.0;
86	}
87
88	if ( ndaug != 2 ) {
89
90	double wtmax=0.0;
91	double pm[5][30],to[4],pmin,pmax,psum,rnd[30];
92	double ran,wt,pa,costh,sinth,phi,p[4][30],be[4],bep,temp;
93	int i,il,ilr,i1,il1u,il1,il2r,ilu;
94	int il2=0;
95
96	for(i=0;i<ndaug;i++){
97	pm[4][i]=0.0;
98	rnd[i]=0.0;
99	}
100
101	pm[0][0]=mp;
102	pm[1][0]=0.0;
103	pm[2][0]=0.0;
104	pm[3][0]=0.0;
105	pm[4][0]=mp;
106
107	to[0]=mp;
108	to[1]=0.0;
109	to[2]=0.0;
110	to[3]=0.0;
111
112	psum=0.0;
113	for(i=1;i<ndaug+1;i++){
114	psum=psum+mass[i-1];
115	}
116
117	pm[4][ndaug-1]=mass[ndaug-1];
118
119	switch (ndaug) {
120
121	case 1:
122	wtmax=1.0/16.0;
123	break;
124	case 2:
125	wtmax=1.0/150.0;
126	break;
127	case 3:
128	wtmax=1.0/2.0;
129	break;
130	case 4:
131	wtmax=1.0/5.0;
132	break;
133	case 5:
134	wtmax=1.0/15.0;
135	break;
136	case 6:
137	wtmax=1.0/15.0;
138	break;
139	case 7:
140	wtmax=1.0/15.0;
141	break;
142	case 8:
143	wtmax=1.0/15.0;
144	break;
145	case 9:
146	wtmax=1.0/15.0;
147	break;
148	case 10:
149	wtmax=1.0/15.0;
150	break;
151	case 11:
152	wtmax=1.0/15.0;
153	break;
154	case 12:
155	wtmax=1.0/15.0;
156	break;
157	case 13:
158	wtmax=1.0/15.0;
159	break;
160	case 14:
161	wtmax=1.0/15.0;
162	break;
163	case 15:
164	wtmax=1.0/15.0;
165	break;
166	default:
167	report(ERROR,"EvtGen") << "too many daughters for phase space..." << ndaug << " "<< mp <<endl;;
168	break;
169	}
170
171	pmax=mp-psum+mass[ndaug-1];
172
173	pmin=0.0;
174
175	for(ilr=2;ilr<ndaug+1;ilr++){
176	il=ndaug+1-ilr;
177	pmax=pmax+mass[il-1];
178	pmin=pmin+mass[il+1-1];
179	wtmax=wtmax*EvtPawt(pmax,pmin,mass[il-1]);
180	}
181
182	do{
183
184	rnd[0]=1.0;
185	il1u=ndaug-1;
186
187	for (il1=2;il1<il1u+1;il1++){
188	ran=EvtRandom::Flat();
189	for (il2r=2;il2r<il1+1;il2r++){
190	il2=il1+1-il2r;
191	if (ran<=rnd[il2-1]) goto two39;
192	rnd[il2+1-1]=rnd[il2-1];
193	}
194	two39:
195	rnd[il2+1-1]=ran;
196	}
197
198	rnd[ndaug-1]=0.0;
199	wt=1.0;
200	for(ilr=2;ilr<ndaug+1;ilr++){
201	il=ndaug+1-ilr;
202	pm[4][il-1]=pm[4][il+1-1]+mass[il-1]+
203	(rnd[il-1]-rnd[il+1-1])*(mp-psum);
204	wt=wt*EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
205	}
206	if (wt>wtmax) {
207	report(ERROR,"EvtGen") << "wtmax to small in EvtPhaseSpace with "
208	<< ndaug <<" daughters"<<endl;;
209	}
210	} while (wt<EvtRandom::Flat(wtmax));
211
212	ilu=ndaug-1;
213
214	for (il=1;il<ilu+1;il++){
215	pa=EvtPawt(pm[4][il-1],pm[4][il+1-1],mass[il-1]);
216	costh=EvtRandom::Flat(-1.0,1.0);
217	sinth=sqrt(1.0-costh*costh);
218	phi=EvtRandom::Flat(EvtConst::twoPi);
219	p[1][il-1]=pasinthcos(phi);
220	p[2][il-1]=pasinthsin(phi);
221	p[3][il-1]=pa*costh;
222	pm[1][il+1-1]=-p[1][il-1];
223	pm[2][il+1-1]=-p[2][il-1];
224	pm[3][il+1-1]=-p[3][il-1];
225	p[0][il-1]=sqrt(papa+mass[il-1]mass[il-1]);
226	pm[0][il+1-1]=sqrt(papa+pm[4][il+1-1]pm[4][il+1-1]);
227	}
228
229	p[0][ndaug-1]=pm[0][ndaug-1];
230	p[1][ndaug-1]=pm[1][ndaug-1];
231	p[2][ndaug-1]=pm[2][ndaug-1];
232	p[3][ndaug-1]=pm[3][ndaug-1];
233
234	for (ilr=2;ilr<ndaug+1;ilr++){
235	il=ndaug+1-ilr;
236	be[0]=pm[0][il-1]/pm[4][il-1];
237	be[1]=pm[1][il-1]/pm[4][il-1];
238	be[2]=pm[2][il-1]/pm[4][il-1];
239	be[3]=pm[3][il-1]/pm[4][il-1];
240
241	for (i1=il;i1<ndaug+1;i1++){
242	bep=be[1]p[1][i1-1]+be[2]p[2][i1-1]+
243	be[3]p[3][i1-1]+be[0]p[0][i1-1];
244	temp=(p[0][i1-1]+bep)/(be[0]+1.0);
245	p[1][i1-1]=p[1][i1-1]+temp*be[1];
246	p[2][i1-1]=p[2][i1-1]+temp*be[2];
247	p[3][i1-1]=p[3][i1-1]+temp*be[3];
248	p[0][i1-1]=bep;
249
250	}
251	}
252
253	for (ilr=0;ilr<ndaug;ilr++){
254	p4[ilr].set(p[0][ilr],p[1][ilr],p[2][ilr],p[3][ilr]);
255	}
256
257	return 1.0;
258	}
259
260	return 1.0;
261
262	}
263
264
265	double EvtGenKine::PhaseSpacePole(double M, double m1, double m2, double m3,
266	double a,EvtVector4R p4[10])
267
268	// generate kinematics for 3 body decays, pole for the m1,m2 mass.
269
270	{
271
272	//f1 = 1 (phasespace)
273	//f2 = a*(1/m12sq)^2
274
275	double m12sqmax=(M-m3)*(M-m3);
276	double m12sqmin=(m1+m2)*(m1+m2);
277
278	double m13sqmax=(M-m2)*(M-m2);
279	double m13sqmin=(m1+m3)*(m1+m3);
280
281	double v1=(m12sqmax-m12sqmin)*(m13sqmax-m13sqmin);
282	double v2= a(1.0/m12sqmin-1.0/m12sqmax)(m13sqmax-m13sqmin);
283
284	double m12sq,m13sq;
285
286	double r=v1/(v1+v2);
287
288	double m13min,m13max;
289
290	do{
291
292	m13sq=EvtRandom::Flat(m13sqmin,m13sqmax);
293
294	if (r>EvtRandom::Flat()){
295	m12sq=EvtRandom::Flat(m12sqmin,m12sqmax);
296	}
297	else{
298	m12sq=1.0/(1.0/m12sqmin-EvtRandom::Flat()*(1.0/m12sqmin-1.0/m12sqmax));
299	}
300
301	//kinematically allowed?
302	double E3star=(MM-m12sq-m3m3)/sqrt(4*m12sq);
303	double E1star=(m12sq+m1m1-m2m2)/sqrt(4*m12sq);
304	double p3star=sqrt(E3starE3star-m3m3);
305	double p1star=sqrt(E1starE1star-m1m1);
306	m13max=(E3star+E1star)*(E3star+E1star)-
307	(p3star-p1star)*(p3star-p1star);
308	m13min=(E3star+E1star)*(E3star+E1star)-
309	(p3star+p1star)*(p3star+p1star);
310
311	}while(m13sq<m13min\|\|m13sq>m13max);
312
313
314	double E2=(MM+m2m2-m13sq)/(2.0*M);
315	double E3=(MM+m3m3-m12sq)/(2.0*M);
316	double E1=M-E2-E3;
317	double p1mom=sqrt(E1E1-m1m1);
318	double p3mom=sqrt(E3E3-m3m3);
319	double cost13=(2.0E1E3+m1m1+m3m3-m13sq)/(2.0p1momp3mom);
320
321	//report(INFO,"EvtGen") << m13sq << endl;
322	//report(INFO,"EvtGen") << m12sq << endl;
323	//report(INFO,"EvtGen") << E1 << endl;
324	//report(INFO,"EvtGen") << E2 << endl;
325	//report(INFO,"EvtGen") << E3 << endl;
326	//report(INFO,"EvtGen") << p1mom << endl;
327	//report(INFO,"EvtGen") << p3mom << endl;
328	//report(INFO,"EvtGen") << cost13 << endl;
329
330
331	p4[2].set(E3,0.0,0.0,p3mom);
332	p4[0].set(E1,p1momsqrt(1.0-cost13cost13),0.0,p1mom*cost13);
333	p4[1].set(E2,-p1momsqrt(1.0-cost13cost13),0.0,-p1mom*cost13-p3mom);
334
335	//report(INFO,"EvtGen") << "p4:"<<p4[0]<<p4[1]<<p4[2]<<endl;
336
337	double alpha = EvtRandom::Flat( EvtConst::twoPi );
338	double beta = acos(EvtRandom::Flat( -1.0, 1.0 ));
339	double gamma = EvtRandom::Flat( EvtConst::twoPi );
340
341	p4[0].applyRotateEuler( alpha, beta, gamma );
342	p4[1].applyRotateEuler( alpha, beta, gamma );
343	p4[2].applyRotateEuler( alpha, beta, gamma );
344
345	return 1.0+a/(m12sq*m12sq);
346
347	}
348