BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsTopipi0pi0.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of the EvtGen package developed jointly
4// for the BaBar and CLEO collaborations. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/COPYRIGHT
8// Copyright (C) 1998 Caltech, UCSB
9//
10// Module: EvtDsTopipi0pi0.cc
11//
12// Description: Routine to handle three-body decays of Ds+/Ds-
13//
14// Modification history:
15// Liaoyuan Dong Sun Jun 02 01:23:25 2024 Module created
16//------------------------------------------------------------------------
18#include <fstream>
19#include <stdlib.h>
22#include "EvtGenBase/EvtPDL.hh"
27#include <string>
30using std::endl;
31
33
34void EvtDsTopipi0pi0::getName(std::string& model_name){
35 model_name="DsTopipi0pi0";
36}
37
41
43 // check that there are 0 arguments
44 checkNArg(0);
45 checkNDaug(3);
50
51 phi[0] = 0;
52 rho[0] = 1.0;
53 phi[1] = 0;
54 rho[1] = 0;
55 phi[2] = 0;
56 rho[2] = 0;
57 phi[3] = 0;
58 rho[3] = 0;
59 phi[4] = 0;
60 rho[4] = 0;
61
62 phi[0] = 0.0;
63 phi[1] = -7.3409e-01;
64 phi[2] = -9.4467e-01;
65 phi[3] = -4.3801e+00;
66 phi[4] = -2.3138e+00;
67
68 rho[0] = 1.0000e+00;
69 rho[1] = 8.5416e-01;
70 rho[2] = 6.8573e-01;
71 rho[3] = 1.7286e+00;
72 rho[4] = 7.5781e-01;
73
74 modetype[0]= 23;
75 modetype[1]= 23;
76 modetype[2]= 23;
77 modetype[3]= 23;
78 modetype[4]= 12;
79
80 //cout << "DsTopipi0pi0 :" << endl;
81 //for (int i=0; i<5; i++) {
82 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
83 //}
84
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;
90
91 mass[0] = 9.5500e-01;
92 mass[1] = 1.3500e+00;
93 mass[2] = 1.2755e+00;
94 mass[3] = 1.0000e+00;
95 mass[4] = 1.0000e+00;
96
97 mD = 1.86486;
98 mDs = 1.9683;
99 rRes = 9.0;
100 rD = 5.0;
101 metap = 0.95778;
102 mkstr = 0.89594;
103 mk0 = 0.497614;
104 mass_Kaon = 0.49368;
105 mass_Pion = 0.13957;
106 mass_Pi0 = 0.1349766;
107 math_pi = 3.1415926;
108 ma0 = 0.99;
109 Ga0 = 0.0756;
110 meta = 0.547862;
111
112 GS1 = 0.636619783;
113 GS2 = 0.01860182466;
114 GS3 = 0.1591549458; // 1/(2*math_2pi)
115 GS4 = 0.00620060822; // mass_Pion2/math_pi
116
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++) {
120 G[i][j] = GG[i][j];
121 }
122 }
123}
124
128
130/*
131 double maxprob = 0.0;
132 for(int ir=0;ir<=60000000;ir++){
133 p->initializePhaseSpace(getNDaug(),getDaugs());
134 EvtVector4R D1 = p->getDaug(0)->getP4();
135 EvtVector4R D2 = p->getDaug(1)->getP4();
136 EvtVector4R D3 = p->getDaug(2)->getP4();
137
138 double P1[4], P2[4], P3[4];
139 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
140 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
141 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
142
143//P1[0] = 0.420940262140625; P1[1] = 0.186467651715440; P1[2] =-0.006370318623730;P1[3] =0.350571109833864;
144//P2[0] = 0.882122776747018; P2[1] =0.132998408143823; P2[2] =-0.666777627124688; P2[3] =-0.545566661456526;
145//P3[0] = 0.712403566511058; P3[1] =-0.587918056029798; P3[2] =0.378201754131628; P3[3] =0.024818375264997;
146
147 double value;
148 //value = calDalEva(P1, P2, P3);
149 int g0[5]={3,1,1,0,0};
150 int spin[5]={0,0,2,2,2};
151 int nstates=5;
152 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
153 if (value<0) continue;
154 if(value>maxprob) {
155 maxprob=value;
156 cout << "ir= " << ir << endl;
157// cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<< endl;
158// cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< endl;
159// cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< endl;
160 cout << "MAX====> " << maxprob << endl;
161 }
162 }
163 printf("MAXprob = %.10f\n",maxprob);
164*/
166 EvtVector4R D1 = p->getDaug(0)->getP4();
167 EvtVector4R D2 = p->getDaug(1)->getP4();
168 EvtVector4R D3 = p->getDaug(2)->getP4();
169
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);
174
175//P1[0] = 0.420940262140625; P1[1] = 0.186467651715440; P1[2] =-0.006370318623730;P1[3] =0.350571109833864;
176//P2[0] = 0.882122776747018; P2[1] =0.132998408143823; P2[2] =-0.666777627124688; P2[3] =-0.545566661456526;
177//P3[0] = 0.712403566511058; P3[1] =-0.587918056029798; P3[2] =0.378201754131628; P3[3] =0.024818375264997;
178
179 double value;
180 int g0[5]={3,1,1,0,0};
181 int spin[5]={0,0,2,2,2};
182 int nstates=5;
183
184 calEva(P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value);
185
186 setProb(value);
187 return ;
188
189}
190
191void EvtDsTopipi0pi0::Com_Multi(double a1[2], double a2[2], double res[2])
192{
193 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
194 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
195}
196void EvtDsTopipi0pi0::Com_Divide(double a1[2], double a2[2], double res[2])
197{
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;
201}
202//------------base---------------------------------
203double EvtDsTopipi0pi0::SCADot(double a1[4], double a2[4])
204{
205 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
206 return _cal;
207}
208double EvtDsTopipi0pi0::barrier(int l, double sa, double sb, double sc, double r, double mass)
209{
210 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
211 if(q < 0) q = 1e-16;
212 double z;
213 z = q*r*r;
214 double sa0;
215 sa0 = mass*mass;
216 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
217 if(q0 < 0) q0 = 1e-16;
218 double z0 = q0*r*r;
219 double F = 0.0;
220 if(l == 0) F = 1;
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));
223 return F;
224}
225
226void EvtDsTopipi0pi0::calt1(double daug1[4], double daug2[4], double t1[4])
227{
228 double p, pq, tmp;
229 double pa[4], qa[4];
230 for(int i=0; i<4; i++) {
231 pa[i] = daug1[i] + daug2[i];
232 qa[i] = daug1[i] - daug2[i];
233 }
234 p = SCADot(pa,pa);
235 pq = SCADot(pa,qa);
236 tmp = pq/p;
237 for(int i=0; i<4; i++) {
238 t1[i] = qa[i] - tmp*pa[i];
239 }
240}
241void EvtDsTopipi0pi0::calt2(double daug1[4], double daug2[4], double t2[4][4])
242{
243 double p, r;
244 double pa[4], t1[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];
249 }
250 p = SCADot(pa,pa);
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);
254 }
255 }
256}
257//-------------------prop--------------------------------------------
258
259double EvtDsTopipi0pi0::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
260{
261 double widm = 0.;
262 double m = sqrt(sa);
263 double tmp = sb-sc;
264 double tmp1 = sa+tmp;
265 double q = 0.25*tmp1*tmp1/sa-sb;
266 if(q<0) q = 1e-16;
267 double tmp2 = mass2+tmp;
268 double q0 = 0.25*tmp2*tmp2/mass2-sb;
269 if(q0<0) q0 = 1e-16;
270 double z = q*r2;
271 double z0 = q0*r2;
272 double t = q/q0;
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);}
276 return widm;
277}
278double EvtDsTopipi0pi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
279{
280 double widm = 0.;
281 double m = sqrt(sa);
282 double tmp = sb-sc;
283 double tmp1 = sa+tmp;
284 double q = 0.25*tmp1*tmp1/sa-sb;
285 if(q<0) q = 1e-16;
286 double tmp2 = mass2+tmp;
287 double q0 = 0.25*tmp2*tmp2/mass2-sb;
288 if(q0<0) q0 = 1e-16;
289 double z = q*r2;
290 double z0 = q0*r2;
291 double F = (1+z0)/(1+z);
292 double t = q/q0;
293 widm = t*sqrt(t)*mass/m*F;
294 return widm;
295}
296void EvtDsTopipi0pi0::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
297{
298 double a[2], b[2];
299 a[0] = 1;
300 a[1] = 0;
301 b[0] = mass2-sa;
302 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
303 Com_Divide(a,b,prop);
304}
305
306void EvtDsTopipi0pi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
307{
308 double a[2], b[2];
309 double tmp = sb-sc;
310 double tmp1 = sa+tmp;
311 double q2 = 0.25*tmp1*tmp1/sa-sb;
312 if(q2<0) q2 = 1e-16;
313
314 double tmp2 = mass2+tmp;
315 double q02 = 0.25*tmp2*tmp2/mass2-sb;
316 if(q02<0) q02 = 1e-16;
317
318 double q = sqrt(q2);
319 double q0 = sqrt(q02);
320 double m = sqrt(sa);
321 double q03 = q0*q02;
322 double tmp3 = log(mass+2*q0)+1.2926305904; // log(mpi0+mpip) = -1.2926305904;
323
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);
329
330 a[0] = 1.0+d*width/mass;
331 a[1] = 0.0;
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);
335}
336
337void EvtDsTopipi0pi0::PiPiSWASS(double sa, double sb, double sc, double prop[2]) {
338 prop[0] = 1;
339 prop[1] = 0;
340}
341
342void EvtDsTopipi0pi0::Flatte_rhoab(double sa, double sb, double rho[2])
343{
344 double q = 1.0-(4*sb/sa);
345
346 if(q>0) {
347 rho[0]=sqrt(q);
348 rho[1]=0;
349 }
350 else if(q<0){
351 rho[0]=0;
352 rho[1]=sqrt(-q);
353 }
354}
355
356void EvtDsTopipi0pi0:: propagator980(double mass, double sx, double *sb, double prop[2])
357{
358
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;
363
364 double unit[2]={1.0};
365 double ci[2]={0,1};
366 double rho1[2];
367 Flatte_rhoab(sx,sb[0],rho1);
368 double rho2[2];
369 Flatte_rhoab(sx,sb[1],rho2);
370 double rho3[2];
371 Flatte_rhoab(sx,sb[2],rho3);
372 double rho4[2];
373 Flatte_rhoab(sx,sb[3],rho4);
374
375 double tmp1[2]={gpipi1,0};
376 double tmp11[2];
377 double tmp2[2]={gpipi2,0};
378 double tmp22[2];
379 double tmp3[2]={gKK1,0};
380 double tmp33[2];
381 double tmp4[2]={gKK2,0};
382 double tmp44[2];
383
384 Com_Multi(tmp1,rho1,tmp11);
385 Com_Multi(tmp2,rho2,tmp22);
386 Com_Multi(tmp3,rho3,tmp33);
387 Com_Multi(tmp4,rho4,tmp44);
388
389 double tmp5[2]={tmp11[0]+tmp22[0]+tmp33[0]+tmp44[0],tmp11[1]+tmp22[1]+tmp33[1]+tmp44[1]};
390 double tmp51[2];
391 Com_Multi(tmp5, ci,tmp51);
392// double tmp6[2]={mass*mass-sx-mass*tmp51[0], -1.0*mass*tmp51[1]};
393 double tmp6[2]={mass*mass-sx-tmp51[0], -1.0*tmp51[1]};
394 Com_Divide( unit,tmp6, prop);
395}
396
397double EvtDsTopipi0pi0::DDalitz(double P1[4], double P2[4], double P3[4], int Ang, double mass){
398 double pR[4], pD[4];
399 double temp_PDF, v_re;
400 temp_PDF = 0.0;
401 v_re = 0.0;
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];
406 }
407 s1 = SCADot(P1,P1);
408 s2 = SCADot(P2,P2);
409 s3 = SCADot(P3,P3);
410 sR = SCADot(pR,pR);
411 sD = SCADot(pD,pD);
412 //int G[4][4];
413 //for(int i=0; i!=4; i++){
414 // for(int j=0; j!=4; j++){
415 // if(i==j){
416 // if(i==0) G[i][j] = 1;
417 // else G[i][j] = -1;
418 // }
419 // else G[i][j] = 0;
420 // }
421 //}
422 if(Ang == 0){
423 B[0] = 1;
424 B[1] = 1;
425 temp_PDF = 1;
426 }
427 if(Ang == 1){
428 B[0] = barrier(1,sR,s1,s2,3.0,mass);
429 B[1] = barrier(1,sD,sR,s3,5.0,1.9683);
430 double t1[4], T1[4];
431 calt1(P1,P2,t1);
432 calt1(pR,P3,T1);
433 temp_PDF = 0;
434 for(int i=0; i<4; i++){
435 temp_PDF += t1[i]*T1[i]*G[i][i];
436 }
437 }
438 if(Ang == 2){
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];
442 calt2(P1,P2,t2);
443 calt2(pR,P3,T2);
444 temp_PDF = 0;
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];
448 }
449 }
450 }
451 v_re = temp_PDF*B[0]*B[1];
452 return v_re;
453}
454
455
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)
457{
458 double P12[4], P23[4], P13[4];
459 double cof[2], amp_PDF[2], PDF[2];
460 double scpi, spi01, spi02;
461 double s12, s13, s23;
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];
466 }
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};
474
475 double pro[2], temp_PDF, amp_tmp[2],temp_PDF1 ,temp_PDF2,pro1[2],pro2[2];
476 double mass1sq;
477 amp_PDF[0] = 0;
478 amp_PDF[1] = 0;
479 PDF[0] = 0;
480 PDF[1] = 0;
481 amp_tmp[0] = 0;
482 amp_tmp[1] = 0;
483 for(int i=0; i<nstates; i++) {
484 amp_tmp[0] = 0;
485 amp_tmp[1] = 0;
486 mass1sq = mass1[i]*mass1[i];
487 cof[0] = amp[i]*cos(phase[i]);
488 cof[1] = amp[i]*sin(phase[i]);
489 temp_PDF = 0;
490
491 if(modetype[i] == 23){//pi0pi0
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);
496 if(g0[i]==0){
497 pro[0] = 1;
498 pro[1] = 0;
499 }
500// printf("cof1 = %.15f \n",cof[1]);
501
502 amp_tmp[0] = temp_PDF*pro[0];
503 amp_tmp[1] = temp_PDF*pro[1];
504
505 }
506 if(modetype[i] == 12){//pipi0
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);
510 if(g0[i]==0){
511 pro1[0] = 1;
512 pro1[1] = 0;
513 }
514
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);
518 if(g0[i]==0){
519 pro2[0] = 1;
520 pro2[1] = 0;
521 }
522 amp_tmp[0] = temp_PDF1*pro1[0] + temp_PDF2*pro2[0];
523 amp_tmp[1] = temp_PDF1*pro1[1] + temp_PDF2*pro2[1];
524
525 }
526 Com_Multi(amp_tmp,cof,amp_PDF);
527 PDF[0] += amp_PDF[0];
528 PDF[1] += amp_PDF[1];
529
530
531}
532
533 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
534 //printf("value = %.15f\n",value);
535 Result = value;
536 }
double sin(const BesAngle a)
Definition BesAngle.h:210
double cos(const BesAngle a)
Definition BesAngle.h:213
double mass
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
Definition FoamA.h:90
****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
Definition KKsem.h:33
TTree * t
Definition binning.cxx:23
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
void checkSpinParent(EvtSpinType::spintype sp)
void setProbMax(double prbmx)
void checkNDaug(int d1, int d2=-1)
EvtId * getDaugs()
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)
EvtDecayBase * clone()
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
double get(int i) const
double double double double * s12
Definition qcdloop1.h:77
double double double double double * s23
Definition qcdloop1.h:77
const double b
Definition slope.cxx:9