35 model_name=
"DsTopipi0pi0";
85 width[0] = 6.7400e-02;
86 width[1] = 2.6500e-01;
87 width[2] = 1.8500e-01;
88 width[3] = 1.7400e-01;
89 width[4] = 1.7400e-01;
106 mass_Pi0 = 0.1349766;
117 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
118 for (
int i=0; i<4; i++) {
119 for (
int j=0; j<4; j++) {
170 double P1[4], P2[4], P3[4];
171 P1[0] = D1.
get(0); P1[1] = D1.
get(1); P1[2] = D1.
get(2); P1[3] = D1.
get(3);
172 P2[0] = D2.
get(0); P2[1] = D2.
get(1); P2[2] = D2.
get(2); P2[3] = D2.
get(3);
173 P3[0] = D3.
get(0); P3[1] = D3.
get(1); P3[2] = D3.
get(2); P3[3] = D3.
get(3);
180 int g0[5]={3,1,1,0,0};
181 int spin[5]={0,0,2,2,2};
184 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
191void EvtDsTopipi0pi0::Com_Multi(
double a1[2],
double a2[2],
double res[2])
193 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
194 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
196void EvtDsTopipi0pi0::Com_Divide(
double a1[2],
double a2[2],
double res[2])
198 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
199 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
200 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
203double EvtDsTopipi0pi0::SCADot(
double a1[4],
double a2[4])
205 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
208double EvtDsTopipi0pi0::barrier(
int l,
double sa,
double sb,
double sc,
double r,
double mass)
210 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
216 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
217 if(q0 < 0) q0 = 1e-16;
221 if(l == 1) F = sqrt((1+z0)/(1+z));
222 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
226void EvtDsTopipi0pi0::calt1(
double daug1[4],
double daug2[4],
double t1[4])
230 for(
int i=0; i<4; i++) {
231 pa[i] = daug1[i] + daug2[i];
232 qa[i] = daug1[i] - daug2[i];
237 for(
int i=0; i<4; i++) {
238 t1[i] = qa[i] - tmp*pa[i];
241void EvtDsTopipi0pi0::calt2(
double daug1[4],
double daug2[4],
double t2[4][4])
245 calt1(daug1,daug2,t1);
246 r = SCADot(t1,t1)/3.0;
247 for(
int i=0; i<4; i++) {
248 pa[i] = daug1[i] + daug2[i];
251 for(
int i=0; i<4; i++) {
252 for(
int j=0; j<4; j++) {
253 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
259double EvtDsTopipi0pi0::wid(
double mass2,
double mass,
double sa,
double sb,
double sc,
double r2,
int l)
264 double tmp1 = sa+tmp;
265 double q = 0.25*tmp1*tmp1/sa-sb;
267 double tmp2 = mass2+tmp;
268 double q0 = 0.25*tmp2*tmp2/mass2-sb;
273 if(l == 0) {widm = sqrt(
t)*
mass/m;}
274 else if(l == 1) {widm =
t*sqrt(
t)*
mass/m*(1+z0)/(1+z);}
275 else if(l == 2) {widm =
t*
t*sqrt(
t)*
mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
278double EvtDsTopipi0pi0::widl1(
double mass2,
double mass,
double sa,
double sb,
double sc,
double r2)
283 double tmp1 = sa+tmp;
284 double q = 0.25*tmp1*tmp1/sa-sb;
286 double tmp2 = mass2+tmp;
287 double q0 = 0.25*tmp2*tmp2/mass2-sb;
291 double F = (1+z0)/(1+z);
293 widm =
t*sqrt(
t)*
mass/m*F;
296void EvtDsTopipi0pi0::propagatorRBW(
double mass2,
double mass,
double width,
double sa,
double sb,
double sc,
double r2,
int l,
double prop[2])
302 b[1] = -
mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
303 Com_Divide(a,
b,prop);
306void EvtDsTopipi0pi0::propagatorGS(
double mass2,
double mass,
double width,
double sa,
double sb,
double sc,
double r2,
double prop[2])
310 double tmp1 = sa+tmp;
311 double q2 = 0.25*tmp1*tmp1/sa-sb;
314 double tmp2 = mass2+tmp;
315 double q02 = 0.25*tmp2*tmp2/mass2-sb;
316 if(q02<0) q02 = 1e-16;
319 double q0 = sqrt(q02);
322 double tmp3 = log(mass+2*q0)+1.2926305904;
324 double h = GS1*
q/m*(log(m+2*
q)+1.2926305904);
325 double h0 = GS1*q0/
mass*tmp3;
326 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
327 double d = GS2/q02*tmp3+GS3*
mass/q0-GS4*
mass/q03;
328 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
330 a[0] = 1.0+d*width/
mass;
332 b[0] = mass2-sa+width*
f;
333 b[1] = -
mass*width*widl1(mass2,mass,sa,sb,sc,r2);
334 Com_Divide(a,
b,prop);
337void EvtDsTopipi0pi0::PiPiSWASS(
double sa,
double sb,
double sc,
double prop[2]) {
342void EvtDsTopipi0pi0::Flatte_rhoab(
double sa,
double sb,
double rho[2])
344 double q = 1.0-(4*sb/sa);
356void EvtDsTopipi0pi0:: propagator980(
double mass,
double sx,
double *sb,
double prop[2])
359 double gpipi1 = 2.0/3.0*0.165;
360 double gpipi2 = 1.0/3.0*0.165;
361 double gKK1 = 0.5*0.69465;
362 double gKK2 = 0.5*0.69465;
364 double unit[2]={1.0};
367 Flatte_rhoab(sx,sb[0],rho1);
369 Flatte_rhoab(sx,sb[1],rho2);
371 Flatte_rhoab(sx,sb[2],rho3);
373 Flatte_rhoab(sx,sb[3],rho4);
375 double tmp1[2]={gpipi1,0};
377 double tmp2[2]={gpipi2,0};
379 double tmp3[2]={gKK1,0};
381 double tmp4[2]={gKK2,0};
384 Com_Multi(tmp1,rho1,tmp11);
385 Com_Multi(tmp2,rho2,tmp22);
386 Com_Multi(tmp3,rho3,tmp33);
387 Com_Multi(tmp4,rho4,tmp44);
389 double tmp5[2]={tmp11[0]+tmp22[0]+tmp33[0]+tmp44[0],tmp11[1]+tmp22[1]+tmp33[1]+tmp44[1]};
391 Com_Multi(tmp5, ci,tmp51);
393 double tmp6[2]={
mass*
mass-sx-tmp51[0], -1.0*tmp51[1]};
394 Com_Divide(
unit,tmp6, prop);
397double EvtDsTopipi0pi0::DDalitz(
double P1[4],
double P2[4],
double P3[4],
int Ang,
double mass){
399 double temp_PDF, v_re;
402 double B[2], s1, s2, s3, sR, sD;
403 for(
int i=0; i<4; i++){
404 pR[i] = P1[i] + P2[i];
405 pD[i] = pR[i] + P3[i];
428 B[0] = barrier(1,sR,s1,s2,3.0,mass);
429 B[1] = barrier(1,sD,sR,s3,5.0,1.9683);
434 for(
int i=0; i<4; i++){
435 temp_PDF += t1[i]*T1[i]*G[i][i];
439 B[0] = barrier(2,sR,s1,s2,3.0,mass);
440 B[1] = barrier(2,sD,sR,s3,5.0,1.9683);
441 double t2[4][4], T2[4][4];
445 for(
int i=0; i<4; i++){
446 for(
int j=0; j<4; j++){
447 temp_PDF += t2[i][j]*T2[j][i]*G[i][i]*G[j][j];
451 v_re = temp_PDF*
B[0]*
B[1];
456void EvtDsTopipi0pi0::calEva(
double* Pic,
double* Pi01,
double* Pi02,
double *mass1,
double *width1,
double *amp,
double *phase,
int* g0,
int* spin,
int* modetype,
int nstates,
double & Result)
458 double P12[4], P23[4], P13[4];
459 double cof[2], amp_PDF[2], PDF[2];
460 double scpi, spi01, spi02;
462 for(
int i=0; i<4; i++){
463 P23[i] = Pi01[i] + Pi02[i];
464 P12[i] = Pic[i] + Pi01[i];
465 P13[i] = Pic[i] + Pi02[i];
467 scpi = SCADot(Pic,Pic);
468 spi01 = SCADot(Pi01,Pi01);
469 spi02 = SCADot(Pi02,Pi02);
470 s12 = SCADot(P12,P12);
471 s13 = SCADot(P13,P13);
472 s23 = SCADot(P23,P23);
473 double spi012[4]={mass_Pion*mass_Pion,spi02,mass_Kaon*mass_Kaon,mk0*mk0};
475 double pro[2], temp_PDF, amp_tmp[2],temp_PDF1 ,temp_PDF2,pro1[2],pro2[2];
483 for(
int i=0; i<nstates; i++) {
486 mass1sq = mass1[i]*mass1[i];
487 cof[0] = amp[i]*
cos(phase[i]);
488 cof[1] = amp[i]*
sin(phase[i]);
491 if(modetype[i] == 23){
492 temp_PDF = DDalitz(Pi01, Pi02, Pic, spin[i], mass1[i]);
493 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],
s23,spi01,spi02,rRes,spin[i],pro);
494 if(g0[i]==12) PiPiSWASS(
s23,spi01,spi02,pro);
495 if(g0[i]==3) propagator980(mass1[i],
s23,spi012,pro);
502 amp_tmp[0] = temp_PDF*pro[0];
503 amp_tmp[1] = temp_PDF*pro[1];
506 if(modetype[i] == 12){
507 temp_PDF1 = DDalitz(Pic, Pi01, Pi02, spin[i], mass1[i]);
508 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],
s12,scpi,spi01,rRes,spin[i],pro1);
509 if(g0[i]==12) PiPiSWASS(
s12,scpi,spi01,pro1);
515 temp_PDF2 = DDalitz(Pic, Pi02, Pi01, spin[i], mass1[i]);
516 if(g0[i]==1) propagatorRBW(mass1sq, mass1[i],width1[i],s13,scpi,spi02,rRes,spin[i],pro2);
517 if(g0[i]==12) PiPiSWASS(s13,scpi,spi02,pro2);
522 amp_tmp[0] = temp_PDF1*pro1[0] + temp_PDF2*pro2[0];
523 amp_tmp[1] = temp_PDF1*pro1[1] + temp_PDF2*pro2[1];
526 Com_Multi(amp_tmp,cof,amp_PDF);
527 PDF[0] += amp_PDF[0];
528 PDF[1] += amp_PDF[1];
533 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
double sin(const BesAngle a)
double cos(const BesAngle a)
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in !Latex Output unit
****INTEGER imax DOUBLE PRECISION m_pi *DOUBLE PRECISION m_amfin DOUBLE PRECISION m_Chfin DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_sinw2 DOUBLE PRECISION m_GFermi DOUBLE PRECISION m_MfinMin DOUBLE PRECISION m_ta2 INTEGER m_out INTEGER m_KeyFSR INTEGER m_KeyQCD *COMMON c_Semalib $ !copy of input $ !CMS energy $ !beam mass $ !final mass $ !beam charge $ !final charge $ !smallest final mass $ !Z mass $ !Z width $ !EW mixing angle $ !Gmu Fermi $ alphaQED at q
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
void checkSpinParent(EvtSpinType::spintype sp)
void setProbMax(double prbmx)
void checkNDaug(int d1, int d2=-1)
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
void setProb(double prob)
virtual ~EvtDsTopipi0pi0()
void decay(EvtParticle *p)
void getName(std::string &name)
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
double double double double * s12
double double double double double * s23