CGEM BOSS 6.6.5.i
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtConExc.cc
Go to the documentation of this file.
1
2//--------------------------------------------------------------------------
3//
4// Environment:
5// This software is part of models developed at BES collaboration
6// based on the EvtGen framework. If you use all or part
7// of it, please give an appropriate acknowledgement.
8//
9// Copyright Information: See EvtGen/BesCopyright
10// Copyright (A) 2006 Ping Rong-Gang @IHEP
11//
12// Module: VVS.hh
13//
14// Description: To define cross section for the continuum exclusive process
15// Experimental cross section taken from PRD73,012005, PRD76,092006, also
16// see a review: Rev. Mod. Phys. 83,1545
17 /*******************--- mode definition:also see EvtXsection.cc
18 0: ppbar
19 1: nnbar
20 2: Lambda0 anti-Lambda0
21 3: Sigma0 anti-Sigma0
22 4: Lambda0 anti-Sigma0
23 5: Sigma0 anti-Lambda0
24 6: pi+ pi-
25 7: pi+ pi- pi0
26 8: K+K- pi0
27 9: KsK+pi-
28 10: KsK-pi+
29 11: K+K-eta
30 12: 2(pi+pi-)
31 13: pi+pi-2pi0
32 14: K+K-pi+pi-
33 15: K+K-2pi0
34 16: 2(K+K-)
35 17: 2(pi+pi-)pi0
36 18: 2(pi+pi-)eta
37 19: K+K-pi+pi-pi0
38 20: K+K-pi+pi-eta
39 21: 3(pi+pi-)
40 22: 2(pi+pi-pi0)
41 23: phi eta
42 24: phi pi0
43 25: K+K*-
44 26: K-K*+
45 27: K_SK*0-bar
46 28: K*0(892)K+pi-
47 29: K*0(892)K-pi+
48 30: K*+K-pi0
49 31: K*-K+pi0
50 32: K_2*(1430)0 K+pi-
51 33: K_2*(1430)0 K-pi+
52 34: K+K-rho
53 35: phi pi+pi-
54 36: phi f0(980)
55 37: eta pi+pi-
56 38: omega pi+ pi-
57 39: omega f0(980)
58 40: eta' pi+ pi-
59 41: f_1(1285)pi+pi-
60 42: omega K+K-
61 43: omega pi+pi-pi0
62 44: Sigma+ Sigma- (Sigma0 Sigma0-bar SU(3) extention )
63 45: K+K-
64 46: K_S K_L
65
66 70: D0D0-bar
67 71: D+D-
68 72: D+D*-
69 73: D-D*+
70 74: D*+D*-
71 75: D0D-pi+
72 76: D0D+pi-
73 77: D0D*-pi+
74 78: D0D*+pi-
75 79: Lambda_c+ Lambda_c-
76
77
78 90: J/psi pi+ pi-
79 91: psi(2S)pi+pi-
80 92: J/psi K+K-
81 93: D_s+ D_s-
82 94: D_s^{*+}D_s^-
83 95: D_s^{*-}D_s^+
84 *************************************/
85// Modification history:
86//
87// Ping R.-G. Nov., 2012 Module created
88//
89//------------------------------------------------------------------------
90//
91#include <math.h>
92#include "EvtGenBase/EvtPatches.hh"
93#include <stdlib.h>
94#include "EvtGenBase/EvtParticle.hh"
95#include "EvtGenBase/EvtGenKine.hh"
96#include "EvtGenBase/EvtPDL.hh"
97#include "EvtGenBase/EvtReport.hh"
98#include "EvtGenModels/EvtPhsp.hh"
99#include "EvtGenModels/EvtConExc.hh"
100#include "EvtGenBase/EvtVector4R.hh"
101#include "EvtGenBase/EvtParticleFactory.hh"
102#include "EvtGenBase/EvtRandom.hh"
103#include "EvtGenBase/EvtHelSys.hh"
104#include "EvtGenBase/EvtDecayTable.hh"
105#include "EvtGenBase/EvtDecayBase.hh"
106#include "EvtGenBase/EvtPDL.hh"
107#include <string>
108#include <iostream>
109#include <fstream>
110#include <istream>
111#include <strstream>
112#include <stdio.h>
113#include <fstream>
114#include <strstream>
115using namespace std;
116
117extern "C" {
118 extern double dgauss_( double (*fun)(double*), double *,double *, double *);
119}
120double Rad2difXs(double *m);
121double Rad2difXs_er(double *m);
122
123EvtXsection *EvtConExc::myxsection;
124double EvtConExc::_cms;
125double EvtConExc::XS_max;
126
127double EvtConExc::_xs0=0;
128double EvtConExc::_xs1=0;
129double EvtConExc::_er0=0;
130double EvtConExc::_er1=0;
131int EvtConExc::_nevt=0;
132
133std::ofstream oa;
134int EvtConExc::getMode;
135
136EvtConExc::~EvtConExc() {
137 if(mydbg){
138 // xs->Write();
139 myfile->Write();
140 }
141 delete myxsection;
142 gamH->deleteTree();
143}
144
145void EvtConExc::getName(std::string& model_name){
146
147 model_name="ConExc";
148
149}
150
151EvtDecayBase* EvtConExc::clone(){
152
153 return new EvtConExc;
154
155}
156
157
158void EvtConExc::init(){
159
160 // check that there are 0 arguments
161 // checkNArg(1);
162 if(getArg(0) == -1){
163 radflag=0;mydbg=false;
164 _mode = getArg(0);
165 pdgcode = getArg(1);
166 pid = EvtPDL::evtIdFromStdHep(pdgcode );
167 nson = getNArg()-2;
168 std::cout<<"The decay daughters are "<<std::endl;
169 for(int i=2;i<getNArg();i++ ){son[i-2]=EvtPDL::evtIdFromStdHep(getArg(i));std::cout<<EvtPDL::name(son[i-2])<<" ";}
170 std::cout<<std::endl;
171 }else if(getArg(0)==-2){
172 radflag=0;mydbg=false;
173 _mode = getArg(0);
174 nson = getNArg()-1;
175 for(int i=1;i<getNArg();i++ ){son[i-1]=EvtPDL::evtIdFromStdHep(getArg(i));}
176 }
177 else if(getNArg()==1){ _mode = getArg(0);radflag=0;mydbg=false;}
178 else if(getNArg()==2){_mode = getArg(0); radflag=getArg(1);mydbg=false;}
179 else if(getNArg()==3){ _mode = getArg(0); radflag=getArg(1);mydbg=getArg(2);}
180 else{std::cout<<"ConExc:umber of parameter should be 1,2 or 3, but you set "<<getNArg()<<std::endl;::abort(); }
181 //--- debugging
182 //std::cout<<_mode<<" "<<pdgcode<<" "<<nson<<" "<<EvtPDL::name(son[0])<<" "<<EvtPDL::name(son[1])<<std::endl;
183
184 gamId = EvtPDL::getId(std::string("gamma"));
185 init_mode(_mode);
186 XS_max=-1;
187 //-----debugging to make root file
188 if(mydbg){
189 myfile = new TFile("mydbg.root","RECREATE");
190 xs=new TTree ("xs","momentum of rad. gamma and hadrons");
191 xs->Branch("imode",&imode,"imode/I");
192 xs->Branch("pgam",pgam,"pgam[4]/D");
193 xs->Branch("phds",phds,"phds[4]/D");
194 xs->Branch("ph1",ph1,"ph1[4]/D");
195 xs->Branch("ph2",ph2,"ph2[4]/D");
196 xs->Branch("mhds",&mhds,"mhds/D");
197 xs->Branch("mass1",&mass1,"mass1/D");
198 xs->Branch("mass2",&mass2,"mass2/D");
199 }
200 //--- prepare the output information
201 EvtId parentId =EvtPDL::getId(std::string("vpho"));
202 double parentMass = EvtPDL::getMass(parentId);
203 // std::cout<<"parent mass = "<<parentMass<<std::endl;
204 EvtVector4R p_init(parentMass,0.0,0.0,0.0);
205 EvtParticle *p=EvtParticleFactory::particleFactory(parentId,p_init);
206 myxsection =new EvtXsection (_mode);
207 double mth=0;
208 double mx = p->getP4().mass();
209 if(_mode ==-1){
210 myxsection->setBW(pdgcode);
211 for(int i=0;i<nson;i++){mth +=EvtPDL::getMeanMass(son[i]);}
212 if(mx<mth){std::cout<<"The CMS energy is lower than the threshold of the final states"<<std::endl;abort();}
213 }else if(_mode == -2){
214 mth=myxsection->getXlw();
215 }else{ mth=myxsection->getXlw();}
216 _cms = mx;
217 _unit = myxsection->getUnit();
218
219 std::cout<<"The specified mode is "<<_mode<<std::endl;
220 EvtConExc::getMode = _mode;
221 //std::cout<<"getXlw= "<<mth<<std::endl;
222
223 double Egamcut = 0.025; //set photon energy cut according to the BESIII detector
224 double Esig = 0.0001; //sigularity engy
225 double x = 2*Egamcut/_cms;
226 double s = _cms*_cms;
227 double mhscut = sqrt(s*(1-x));
228
229 // std::cout<<"mhscut= "<<mhscut<<" _cms= "<<_cms<<" mth= "<<mth<<std::endl;
230 //--- calculated with Gaussian integration
231 /*
232 _xs0 = gamHXSection(mhscut,_cms);
233 _er0 = gamHXSection_er(mhscut,_cms);
234 _xs1 = gamHXSection(mth,mhscut);
235 double xs1_err = gamHXSection_er(mth,mhscut);
236 */
237 //-- calculated with MC integration method
238
239 double mhdL=myxsection->getXlw();
240 double EgamH = (s-mhdL*mhdL)/2/sqrt(s);
241 int nmc=100000;
242 _xs0 = gamHXSection(s,Esig,Egamcut,nmc);
243 _er0 = _er1; // stored when calculate _xs0
244 std::cout<<"_xs0= "<<_er0<<std::endl;
245 _xs1 = gamHXSection(s,Egamcut,EgamH,nmc);
246 double xs1_err = _er1;
247
248 //--- sigularity point x from 0 to 0.0001GeV
249 double xb= 2*Esig/_cms;
250 double sig_xs = SoftPhoton_xs(s,xb)*(myxsection->getXsection(mx));
251 _xs0 += sig_xs;
252
253 //for information output
254 double bornXS = myxsection->getXsection(mx); // std::cout<<"BornXS= "<<bornXS<<std::endl;
255 double bornXS_er=myxsection->getErr(mx);
256 double obsvXS = _xs0 + _xs1; //gamHXSection(mth ,mx);
257 double obsvXS_er= _er0 + xs1_err;
258 double corr = obsvXS/bornXS;
259 double corr_er =corr*sqrt(bornXS_er*bornXS_er/bornXS/bornXS + obsvXS_er*obsvXS_er/obsvXS/obsvXS);
260
261
262 if(bornXS==0){std::cout<<"EvtConExc::init : The Born cross section at this point is zero!"<<std::endl;abort();}
263 if(obsvXS==0){std::cout<<"EvtConExc::init : The observed cross section at this point is zero!"<<std::endl;abort();}
264
265 //_xs0 += bornXS;// events with very soft photon (Egam<0.025) are taken as the born process
266 //_er0 = sqrt(_er0*_er0 + bornXS_er*bornXS_er);
267
268 std::cout<<""<<std::endl;
269 std::cout<<"========= Generation using cross section (Egamcut=0.025 GeV) =============="<<std::endl;
270 std::cout<<" sqrt(s)= "<<mx<< " GeV "<<std::endl;
271 if(_mode>=0 || _mode ==-2 ){
272 std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<"+/-"<<_er0<<" in Unit "<<_unit.c_str()<<std::endl;
273 std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"+/-"<<xs1_err<<" in Unit "<<_unit.c_str()<<std::endl;
274 std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
275 std::cout<<"The radiative correction factor sigma_obs/sigma_born(s) = (1+delta) = "<<corr<<"+/-"<< fabs(corr_er)<<","<<std::endl;
276 std::cout<<"which is calcualted with these quantities:"<<std::endl;
277 std::cout<<"the observed cross section is "<<obsvXS<<"+/-"<< obsvXS_er<<_unit.c_str()<<std::endl;
278 std::cout<<"and the Born cross section (s) is "<<bornXS<<"+/-"<<bornXS_er<<_unit.c_str()<<std::endl;
279 std::cout<<"Within the range "<<myxsection->getXlw()<<" ~"<<myxsection->getXup()<<" GeV, "<<myxsection->getMsg().c_str()<<std::endl;
280 std::cout<<"==========================================================================="<<std::endl;
281 }else if(_mode==-1){
282 std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<" *Br_ee"<<" in Unit "<<_unit.c_str()<<std::endl;
283 std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"*Br_ee in Unit "<<_unit.c_str()<<std::endl;
284 std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
285 std::cout<<"The radiative correction factor sigma_obs/sigma_born(s) = (1+delta) = "<<corr<<"+/-"<< fabs(corr_er)<<std::endl;
286 std::cout<<"==========================================================================="<<std::endl;
287 }
288 std::cout<<std::endl;
289 std::cout<<std::endl;
290
291 findMaxXS(p);
292 //--debugging
293 //std::cout<<"maxXS= "<<maxXS<<std::endl;
294
295 if(_xs0 == 0 && _xs1==0){std::cout<<"EvtConExc::ini() has zero cross section"<<std::endl;abort();}
296
297 init_Br_ee();
298 double mthrld = EvtPDL::getMeanMass(daugs[0]); //threshold mass of hadrons
299 for(int jj=1;jj<_ndaugs;jj++){
300 mthrld += EvtPDL::getMeanMass(daugs[jj]);
301 }
302
303 for(int ii=0;ii<6;ii++){
304 double mR = EvtPDL::getMeanMass(ResId[ii]);
305 if(mx<mR || mR<mthrld) continue;
306 double narRxs=Ros_xs(mx,BR_ee[ii],ResId[ii]);
307 std::cout<<"The corss section for gamma "<<EvtPDL::name(ResId[ii]).c_str()<<" is: "<< narRxs<<"*Br (nb)."<<std::endl;
308 }
309 std::cout<<"==========================================================================="<<std::endl;
310 std::cout<<std::endl;
311 std::cout<<std::endl;
312
313 //--- for debugging
314 if(mydbg){
315 int nd = myxsection->getYY().size();
316 double xx[10000],yy[10000],xer[10000],yer[10000];
317 for(int i=0;i<nd;i++){
318 xx[i]=myxsection->getXX()[i];
319 yy[i]=myxsection->getYY()[i];
320 yer[i]=myxsection->getEr()[i];
321 xer[i]=0.;
322 //std::cout<<"yy "<<yy[i]<<std::endl;
323 }
324 myth=new TH1F("myth","Exp.data",200,xx[0],xx[nd]);
325 for(int i=0;i<nd;i++){
326 myth->Fill(xx[i],yy[i]);
327 }
328 myth->SetError(yer);
329 myth->Write();
330 } //if(mydbg)_block
331 //-------------------------
332
333}
334
335
336//--
337void EvtConExc::init_mode(int mode){
338 if(mode ==0 ){
339 _ndaugs =2;
340 daugs[0]=EvtPDL::getId(std::string("p+"));
341 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
342 }else if(mode ==1 ){
343 _ndaugs =2;
344 daugs[0]=EvtPDL::getId(std::string("n0"));
345 daugs[1]=EvtPDL::getId(std::string("anti-n0"));
346 }else if(mode ==2 ){
347 _ndaugs =2;
348 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
349 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
350 }else if(mode ==3 ){
351 _ndaugs =2;
352 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
353 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
354 }else if(mode ==4 ){
355 _ndaugs =2;
356 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
357 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
358 }else if(mode ==5 ){
359 _ndaugs =2;
360 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
361 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
362 } else if(mode ==6 ){
363 _ndaugs =2;
364 daugs[0]=EvtPDL::getId(std::string("pi+"));
365 daugs[1]=EvtPDL::getId(std::string("pi-"));
366 }else if(mode ==7 ){
367 _ndaugs =3;
368 daugs[0]=EvtPDL::getId(std::string("pi+"));
369 daugs[1]=EvtPDL::getId(std::string("pi-"));
370 daugs[2]=EvtPDL::getId(std::string("pi0"));
371 }else if(mode ==8 ){
372 _ndaugs =3;
373 daugs[0]=EvtPDL::getId(std::string("K+"));
374 daugs[1]=EvtPDL::getId(std::string("K-"));
375 daugs[2]=EvtPDL::getId(std::string("pi0"));
376 }else if(mode ==9 ){
377 _ndaugs =3;
378 daugs[0]=EvtPDL::getId(std::string("K_S0"));
379 daugs[1]=EvtPDL::getId(std::string("K+"));
380 daugs[2]=EvtPDL::getId(std::string("pi-"));
381 }else if(mode ==10 ){
382 _ndaugs =3;
383 daugs[0]=EvtPDL::getId(std::string("K_S0"));
384 daugs[1]=EvtPDL::getId(std::string("K-"));
385 daugs[2]=EvtPDL::getId(std::string("pi+"));
386 }else if(mode ==11 ){
387 _ndaugs =3;
388 daugs[0]=EvtPDL::getId(std::string("K+"));
389 daugs[1]=EvtPDL::getId(std::string("K+"));
390 daugs[2]=EvtPDL::getId(std::string("eta"));
391 }else if(mode ==12 ){
392 _ndaugs =4;
393 daugs[0]=EvtPDL::getId(std::string("pi+"));
394 daugs[1]=EvtPDL::getId(std::string("pi-"));
395 daugs[2]=EvtPDL::getId(std::string("pi+"));
396 daugs[3]=EvtPDL::getId(std::string("pi-"));
397 }else if(mode ==13 ){
398 _ndaugs =4;
399 daugs[0]=EvtPDL::getId(std::string("pi+"));
400 daugs[1]=EvtPDL::getId(std::string("pi-"));
401 daugs[2]=EvtPDL::getId(std::string("pi0"));
402 daugs[3]=EvtPDL::getId(std::string("pi0"));
403 }else if(mode ==14 ){
404 _ndaugs =4;
405 daugs[0]=EvtPDL::getId(std::string("K+"));
406 daugs[1]=EvtPDL::getId(std::string("K-"));
407 daugs[2]=EvtPDL::getId(std::string("pi+"));
408 daugs[3]=EvtPDL::getId(std::string("pi-"));
409 }else if(mode ==15 ){
410 _ndaugs =4;
411 daugs[0]=EvtPDL::getId(std::string("K+"));
412 daugs[1]=EvtPDL::getId(std::string("K-"));
413 daugs[2]=EvtPDL::getId(std::string("pi0"));
414 daugs[3]=EvtPDL::getId(std::string("pi0"));
415 }else if(mode ==16 ){
416 _ndaugs =4;
417 daugs[0]=EvtPDL::getId(std::string("K+"));
418 daugs[1]=EvtPDL::getId(std::string("K-"));
419 daugs[2]=EvtPDL::getId(std::string("K+"));
420 daugs[3]=EvtPDL::getId(std::string("K-"));
421 }else if(mode ==17 ){
422 _ndaugs =5;
423 daugs[0]=EvtPDL::getId(std::string("pi+"));
424 daugs[1]=EvtPDL::getId(std::string("pi-"));
425 daugs[2]=EvtPDL::getId(std::string("pi+"));
426 daugs[3]=EvtPDL::getId(std::string("pi-"));
427 daugs[4]=EvtPDL::getId(std::string("pi0"));
428 }else if(mode ==18 ){
429 _ndaugs =5;
430 daugs[0]=EvtPDL::getId(std::string("pi+"));
431 daugs[1]=EvtPDL::getId(std::string("pi-"));
432 daugs[2]=EvtPDL::getId(std::string("pi+"));
433 daugs[3]=EvtPDL::getId(std::string("pi-"));
434 daugs[4]=EvtPDL::getId(std::string("eta"));
435 }else if(mode ==19 ){
436 _ndaugs =5;
437 daugs[0]=EvtPDL::getId(std::string("K+"));
438 daugs[1]=EvtPDL::getId(std::string("K-"));
439 daugs[2]=EvtPDL::getId(std::string("pi+"));
440 daugs[3]=EvtPDL::getId(std::string("pi-"));
441 daugs[4]=EvtPDL::getId(std::string("pi0"));
442 }else if(mode ==20 ){
443 _ndaugs =5;
444 daugs[0]=EvtPDL::getId(std::string("K+"));
445 daugs[1]=EvtPDL::getId(std::string("K-"));
446 daugs[2]=EvtPDL::getId(std::string("pi+"));
447 daugs[3]=EvtPDL::getId(std::string("pi-"));
448 daugs[4]=EvtPDL::getId(std::string("eta"));
449 }else if(mode ==21 ){
450 _ndaugs =6;
451 daugs[0]=EvtPDL::getId(std::string("pi+"));
452 daugs[1]=EvtPDL::getId(std::string("pi-"));
453 daugs[2]=EvtPDL::getId(std::string("pi+"));
454 daugs[3]=EvtPDL::getId(std::string("pi-"));
455 daugs[4]=EvtPDL::getId(std::string("pi+"));
456 daugs[5]=EvtPDL::getId(std::string("pi-"));
457 }else if(mode ==22 ){
458 _ndaugs =6;
459 daugs[0]=EvtPDL::getId(std::string("pi+"));
460 daugs[1]=EvtPDL::getId(std::string("pi-"));
461 daugs[2]=EvtPDL::getId(std::string("pi+"));
462 daugs[3]=EvtPDL::getId(std::string("pi-"));
463 daugs[4]=EvtPDL::getId(std::string("pi0"));
464 daugs[5]=EvtPDL::getId(std::string("pi0"));
465 }else if(mode == 23){
466 _ndaugs =2;
467 daugs[0]=EvtPDL::getId(std::string("eta"));
468 daugs[1]=EvtPDL::getId(std::string("phi"));
469 }else if(mode == 24){
470 _ndaugs =2;
471 daugs[0]=EvtPDL::getId(std::string("phi"));
472 daugs[1]=EvtPDL::getId(std::string("pi0"));
473 }else if(mode == 25){
474 _ndaugs =2;
475 daugs[0]=EvtPDL::getId(std::string("K+"));
476 daugs[1]=EvtPDL::getId(std::string("K*-"));
477 }else if(mode == 26){
478 _ndaugs =2;
479 daugs[0]=EvtPDL::getId(std::string("K-"));
480 daugs[1]=EvtPDL::getId(std::string("K*+"));
481 }else if(mode == 27){
482 _ndaugs =2;
483 daugs[0]=EvtPDL::getId(std::string("K_S0"));
484 daugs[1]=EvtPDL::getId(std::string("anti-K*0"));
485 }else if(mode == 28){
486 _ndaugs =3;
487 daugs[0]=EvtPDL::getId(std::string("K*0"));
488 daugs[1]=EvtPDL::getId(std::string("K+"));
489 daugs[2]=EvtPDL::getId(std::string("pi-"));
490 }else if(mode == 29){
491 _ndaugs =3;
492 daugs[0]=EvtPDL::getId(std::string("K*0"));
493 daugs[1]=EvtPDL::getId(std::string("K-"));
494 daugs[2]=EvtPDL::getId(std::string("pi+"));
495 }else if(mode == 30){
496 _ndaugs =3;
497 daugs[0]=EvtPDL::getId(std::string("K*+"));
498 daugs[1]=EvtPDL::getId(std::string("K-"));
499 daugs[2]=EvtPDL::getId(std::string("pi0"));
500 }else if(mode == 31){
501 _ndaugs =3;
502 daugs[0]=EvtPDL::getId(std::string("K*-"));
503 daugs[1]=EvtPDL::getId(std::string("K+"));
504 daugs[2]=EvtPDL::getId(std::string("pi0"));
505 }else if(mode == 32){
506 _ndaugs =3;
507 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
508 daugs[1]=EvtPDL::getId(std::string("K+"));
509 daugs[2]=EvtPDL::getId(std::string("pi-"));
510 }else if(mode == 33){
511 _ndaugs =3;
512 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
513 daugs[1]=EvtPDL::getId(std::string("K-"));
514 daugs[2]=EvtPDL::getId(std::string("pi+"));
515 }else if(mode == 34){
516 _ndaugs =3;
517 daugs[0]=EvtPDL::getId(std::string("K+"));
518 daugs[1]=EvtPDL::getId(std::string("K-"));
519 daugs[2]=EvtPDL::getId(std::string("rho0"));
520 }else if(mode == 35){
521 _ndaugs =3;
522 daugs[0]=EvtPDL::getId(std::string("phi"));
523 daugs[1]=EvtPDL::getId(std::string("pi-"));
524 daugs[2]=EvtPDL::getId(std::string("pi+"));
525 }else if(mode == 36){
526 _ndaugs =2;
527 daugs[0]=EvtPDL::getId(std::string("phi"));
528 daugs[1]=EvtPDL::getId(std::string("f_0"));
529 }else if(mode == 37){
530 _ndaugs =3;
531 daugs[0]=EvtPDL::getId(std::string("eta"));
532 daugs[1]=EvtPDL::getId(std::string("pi+"));
533 daugs[2]=EvtPDL::getId(std::string("pi-"));
534 }else if(mode == 38){
535 _ndaugs =3;
536 daugs[0]=EvtPDL::getId(std::string("omega"));
537 daugs[1]=EvtPDL::getId(std::string("pi+"));
538 daugs[2]=EvtPDL::getId(std::string("pi-"));
539 }else if(mode == 39){
540 _ndaugs =2;
541 daugs[0]=EvtPDL::getId(std::string("omega"));
542 daugs[1]=EvtPDL::getId(std::string("f_0"));
543 }else if(mode == 40){
544 _ndaugs =3;
545 daugs[0]=EvtPDL::getId(std::string("eta'"));
546 daugs[1]=EvtPDL::getId(std::string("pi+"));
547 daugs[2]=EvtPDL::getId(std::string("pi-"));
548 }else if(mode == 41){
549 _ndaugs =3;
550 daugs[0]=EvtPDL::getId(std::string("f_1'"));
551 daugs[1]=EvtPDL::getId(std::string("pi+"));
552 daugs[2]=EvtPDL::getId(std::string("pi-"));
553 }else if(mode == 42){
554 _ndaugs =3;
555 daugs[0]=EvtPDL::getId(std::string("omega'"));
556 daugs[1]=EvtPDL::getId(std::string("K+"));
557 daugs[2]=EvtPDL::getId(std::string("K-"));
558 }else if(mode == 43){
559 _ndaugs =4;
560 daugs[0]=EvtPDL::getId(std::string("omega'"));
561 daugs[1]=EvtPDL::getId(std::string("pi+"));
562 daugs[2]=EvtPDL::getId(std::string("pi-"));
563 daugs[3]=EvtPDL::getId(std::string("pi0"));
564 }else if(mode == 44){
565 _ndaugs =2;
566 daugs[0]=EvtPDL::getId(std::string("Sigma-"));
567 daugs[1]=EvtPDL::getId(std::string("anti-Sigma+"));
568 }else if(mode == 45){
569 _ndaugs =2;
570 daugs[0]=EvtPDL::getId(std::string("K+"));
571 daugs[1]=EvtPDL::getId(std::string("K-"));
572 }else if(mode == 46){
573 _ndaugs =2;
574 daugs[0]=EvtPDL::getId(std::string("K_S0"));
575 daugs[1]=EvtPDL::getId(std::string("K_L0"));
576 }else if(mode == 50){
577 _ndaugs =3;
578 daugs[0]=EvtPDL::getId(std::string("p+"));
579 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
580 daugs[2]=EvtPDL::getId(std::string("pi0"));
581 }else if(mode == 51){
582 _ndaugs =3;
583 daugs[0]=EvtPDL::getId(std::string("p+"));
584 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
585 daugs[2]=EvtPDL::getId(std::string("eta"));
586 }else if(mode == 68){
587 _ndaugs =2;
588 daugs[0]=EvtPDL::getId(std::string("D0"));
589 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
590
591 }else if(mode == 69){
592 _ndaugs =2;
593 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
594 daugs[1]=EvtPDL::getId(std::string("D*0"));
595
596 }else if(mode == 70){
597 _ndaugs =2;
598 daugs[0]=EvtPDL::getId(std::string("D0"));
599 daugs[1]=EvtPDL::getId(std::string("anti-D0"));
600
601}else if(mode == 71){
602 _ndaugs =2;
603 daugs[0]=EvtPDL::getId(std::string("D+"));
604 daugs[1]=EvtPDL::getId(std::string("D-"));
605 }else if(mode == 72){
606 _ndaugs =2;
607 daugs[0]=EvtPDL::getId(std::string("D+"));
608 daugs[1]=EvtPDL::getId(std::string("D*-"));
609
610 }else if(mode == 73){
611 _ndaugs =2;
612 daugs[0]=EvtPDL::getId(std::string("D-"));
613 daugs[1]=EvtPDL::getId(std::string("D*+"));
614
615 }else if(mode == 74){
616 _ndaugs =2;
617 daugs[0]=EvtPDL::getId(std::string("D*+"));
618 daugs[1]=EvtPDL::getId(std::string("D*-"));
619
620 }else if(mode == 75){
621 _ndaugs =3;
622 daugs[0]=EvtPDL::getId(std::string("D0"));
623 daugs[1]=EvtPDL::getId(std::string("D-"));
624 daugs[2]=EvtPDL::getId(std::string("pi+"));
625 }else if(mode == 76){
626 _ndaugs =3;
627 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
628 daugs[1]=EvtPDL::getId(std::string("D+"));
629 daugs[2]=EvtPDL::getId(std::string("pi-"));
630 }else if(mode == 77){
631 _ndaugs =3;
632 daugs[0]=EvtPDL::getId(std::string("D0"));
633 daugs[1]=EvtPDL::getId(std::string("D*-"));
634 daugs[2]=EvtPDL::getId(std::string("pi+"));
635 }else if(mode == 78){
636 _ndaugs =3;
637 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
638 daugs[1]=EvtPDL::getId(std::string("D*+"));
639 daugs[2]=EvtPDL::getId(std::string("pi-"));
640 }else if(mode == 79){
641 _ndaugs =2;
642 daugs[0]=EvtPDL::getId(std::string("Lambda_c+"));
643 daugs[1]=EvtPDL::getId(std::string("anti-Lambda_c-"));
644 }else if(mode == 80){
645 _ndaugs =2;
646 daugs[0]=EvtPDL::getId(std::string("eta"));
647 daugs[1]=EvtPDL::getId(std::string("J/psi"));
648 }else if(mode == 81){
649 _ndaugs =3;
650 daugs[0]=EvtPDL::getId(std::string("h_c"));
651 daugs[1]=EvtPDL::getId(std::string("pi+"));
652 daugs[2]=EvtPDL::getId(std::string("pi-"));
653 }else if(mode == 82){
654 _ndaugs =3;
655 daugs[0]=EvtPDL::getId(std::string("h_c"));
656 daugs[1]=EvtPDL::getId(std::string("pi0"));
657 daugs[2]=EvtPDL::getId(std::string("pi0"));
658 }else if(mode == 83){
659 _ndaugs =3;
660 daugs[0]=EvtPDL::getId(std::string("J/psi"));
661 daugs[1]=EvtPDL::getId(std::string("K+"));
662 daugs[2]=EvtPDL::getId(std::string("K-"));
663 }else if(mode == 84){
664 _ndaugs =3;
665 daugs[0]=EvtPDL::getId(std::string("J/psi"));
666 daugs[1]=EvtPDL::getId(std::string("K_S0"));
667 daugs[2]=EvtPDL::getId(std::string("K_S0"));
668 }else if(mode == 90){
669 _ndaugs =3;
670 daugs[0]=EvtPDL::getId(std::string("J/psi"));
671 daugs[1]=EvtPDL::getId(std::string("pi+"));
672 daugs[2]=EvtPDL::getId(std::string("pi-"));
673 }else if(mode == 91){
674 _ndaugs =3;
675 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
676 daugs[1]=EvtPDL::getId(std::string("pi+"));
677 daugs[2]=EvtPDL::getId(std::string("pi-"));
678 }else if(mode == 92){
679 _ndaugs =3;
680 daugs[0]=EvtPDL::getId(std::string("J/psi"));
681 daugs[1]=EvtPDL::getId(std::string("K+"));
682 daugs[2]=EvtPDL::getId(std::string("K-"));
683 }else if(mode == 93){
684 _ndaugs =2;
685 daugs[0]=EvtPDL::getId(std::string("D_s+"));
686 daugs[1]=EvtPDL::getId(std::string("D_s-"));
687 }else if(mode == 94){
688 _ndaugs =2;
689 daugs[0]=EvtPDL::getId(std::string("D_s*+"));
690 daugs[1]=EvtPDL::getId(std::string("D_s-"));
691 }else if(mode == 95){
692 _ndaugs =2;
693 daugs[0]=EvtPDL::getId(std::string("D_s*-"));
694 daugs[1]=EvtPDL::getId(std::string("D_s+"));
695 }else if(mode == -1){
696 _ndaugs = nson;
697 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
698 std::cout<<"The paraent particle: "<<EvtPDL::name(pid)<<std::endl;
699 }else if(mode == -2){
700 _ndaugs = nson;
701 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
702 }else {
703 std::cout<<"Bad mode_index number, please refer to the manual."<<std::endl;
704 ::abort();
705 }
706 // std::cout<<_ndaugs<<", "<<daugs[0]<<daugs[1]<<std::endl;
707
708}
709
710
711void EvtConExc::initProbMax(){
712
713 noProbMax();
714
715}
716
717void EvtConExc::decay( EvtParticle *p ){
718 //--- debugging
719 // std::cout<<"_mode= "<<_mode<<", XS_max= "<<XS_max<<std::endl;
720 //------
721
722 if((_xs0 + _xs1)==0) {std::cout<<"EvtConExc::zero total xsection"<<std::endl;::abort();}
723 double pm= EvtRandom::Flat(0.,1);
724 double xsratio = _xs0/(_xs0 + _xs1);
725 int iflag; //flag = 0 for ee->hadrons, 1 for ee->gamma hadrons
726 if(pm<xsratio){iflag = 0;}else{iflag=1;}
727
728 if(radflag==1){xsratio=1;iflag=1;} // only generating gamma hadrons mode
729
730 daugs[_ndaugs]=gamId;
731 p->makeDaughters(_ndaugs+iflag,daugs);
732
733 for(int i=0;i< _ndaugs+iflag;i++){
734 EvtParticle* di=p->getDaug(i);
735 double xmi=EvtPDL::getMass(di->getId());
736 di->setMass(xmi);
737 }
738 angular_sampling:
739 double weight1 = p->initializePhaseSpace( _ndaugs+iflag , daugs);
740 bool gam_ang,byn_ang,msn_ang,xs_ang;
741 EvtVector4R vhds=p->getDaug(0)->getP4();
742 EvtVector4R hd1=vhds;
743 EvtVector4R hd2 =p->getDaug(1)->getP4();;
744 mass1 = hd1.mass();
745 mass2 = hd2.mass();
746 for(int i=1;i<_ndaugs;i++){
747 vhds += p->getDaug(i)->getP4();
748 }
749 double xm = vhds.mass();
750 double mxs=myxsection->getXsection(xm);
751 xs_ang = xs_sampling(mxs);
752 //---debugging
753 //std::cout<<"mhds= "<<xm<<", xs= "<<mxs<<" XS_max="<<XS_max<<std::endl;
754
755 // if(!xs_ang && iflag==1){ goto angular_sampling;}
756
757 if(iflag ==1 ){
758 gam_ang = gam_sampling(p);
759 if(!gam_ang || !xs_ang ){ goto angular_sampling;}
760 }
761
762 EvtVector4R ppi = p->getDaug(0)->getP4();
763 EvtVector4R pcm = p->getDaug(1)->getP4() + ppi;
764 EvtVector4R pbst=-1*pcm;
765 pbst.set(0,pcm.get(0));
766
767 EvtVector4R pi = boostTo(ppi,pbst);
768 //--debugging
769 //EvtVector4R d2 = boostTo(p->getDaug(1)->getP4(),pbst);
770 //std::cout<<"pi,d2"<<pi<<d2<<std::endl;
771
772 if(_mode<=5 || _mode ==44 ||_mode ==72 ||_mode == 73||_mode==94||_mode==95 ||_mode ==79 || _mode>=23 &&_mode<=27 || _mode==36|| _mode ==39 || _mode == 80){//ee->two baryon mode, VP, SP
773 byn_ang = baryon_sampling(pcm, pi);
774 if(!byn_ang) { goto angular_sampling;}
775 }else if(_mode ==6||_mode==45 || _mode==46 || _mode==70 || _mode==71|| _mode==93){// ee->PP (pi+pi-,..) mode
776 msn_ang = meson_sampling(pcm,pi);
777 if(!msn_ang){goto angular_sampling;}
778 }
779
780 if(_nevt ==0 ){
781 std::cout<<"The decay chain: "<<std::endl;
782 p->printTree();
783 }
784 _nevt++;
785 //--- for debuggint to make root file
786 if(mydbg){
787 imode = iflag;
788 if(iflag==1){
789 EvtVector4R vgam = p->getDaug(_ndaugs)->getP4();
790 pgam[0]=vgam.get(0);
791 pgam[1]=vgam.get(1);
792 pgam[2]=vgam.get(2);
793 pgam[3]=vgam.get(3);
794 }
795
796 ph1[0]=hd1.get(0);
797 ph1[1]=hd1.get(1);
798 ph1[2]=hd1.get(2);
799 ph1[3]=hd1.get(3);
800
801 ph2[0]=hd2.get(0);
802 ph2[1]=hd2.get(1);
803 ph2[2]=hd2.get(2);
804 ph2[3]=hd2.get(3);
805
806 phds[0]=vhds.get(0);
807 phds[1]=vhds.get(1);
808 phds[2]=vhds.get(2);
809 phds[3]=vhds.get(3);
810 mhds = vhds.mass();
811 xs->Fill();
812 //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
813 }
814 //----
815 return ;
816}
817
818
819double EvtConExc::gamHXSection(EvtParticle *p,double El,double Eh,int nmc){//El, Eh : the lower and higher energy of radiative photons
820 //std::cout<<"nmc= "<<nmc<<std::endl;
821 gamId = EvtPDL::getId(std::string("gamma"));
822 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
823 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
824 double totxs = 0;
825 maxXS=-1;
826 _er1=0;
827 Rad2Xs =0;
828 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
829 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
830 int gamdaugs = _ndaugs+1;
831
832 EvtId theDaugs[10];
833 for(int i=0;i<_ndaugs;i++){
834 theDaugs[i] = daugs[i];
835 }
836 theDaugs[_ndaugs]=gamId;
837
838 gamH->makeDaughters(gamdaugs,theDaugs);
839
840 for(int i=0;i<gamdaugs;i++){
841 EvtParticle* di=gamH->getDaug(i);
842 double xmi=EvtPDL::getMass(di->getId());
843 di->setMass(xmi);
844 }
845 loop:
846 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
847 //-- slection:theta_gamma > 1 degree
848 EvtVector4R pgam = gamH -> getDaug(_ndaugs)->getP4();
849 double pmag = pgam.d3mag();
850 double pz = pgam.get(3);
851 double sintheta = sqrt( pmag*pmag - pz*pz)/pmag;
852 double onedegree = 1./180.*3.1415926;
853 //if(sintheta < onedegree) {goto loop;}
854 if(pmag <El || pmag >Eh) {goto loop;}
855 //--------
856 // std::cout<<"pmag= "<<pmag<<std::endl;
857
858 double thexs = difgamXs(gamH); //prepare the photon angular distribution
859 Rad2Xs += Rad2difXs( gamH );
860 if(thexs>maxXS) {maxXS=thexs;}
861 double s = p_init.mass2();
862 double x = 2*pgam.get(0)/sqrt(s);
863 double rad1xs = Rad1difXs(gamH); //first order xs by Eq(4)hep-ph/9910523
864 totxs += rad1xs;
865 _er1 += differ;
866 gamH->deleteDaughters();
867 } //for int_i block
868 _er1/=nmc;
869
870 Rad2Xs/=nmc; // the leading second order correction
871 totxs/=nmc; // first order correction XS
872
873 // return totxs; // first order correction XS
874 return Rad2Xs; // the leading second order correction
875}
876
877
878double EvtConExc::gamHXSection(double s, double El,double Eh,int nmc){//El, Eh : the lower and higher energy of radiative photons
879 double totxs = 0;
880 maxXS=-1;
881 _er1=0;
882 Rad2Xs =0;
883 double xL=2*El/sqrt(s);
884 double xH=2*Eh/sqrt(s);
885 for(int i=0;i<nmc;i++){
886 double rdm = EvtRandom::Flat(0.,1);// std::cout<<"rdm= "<<rdm<<std::endl;
887 double x= xL+ (xH-xL)*rdm;
888 Rad2Xs += Rad2difXs(s,x);
889 _er1 += differ2; //stored when calculate Rad2Xs
890 // std::cout<<"i= "<<i<<" Rad2Xs= "<<Rad2Xs<<std::endl;
891 }
892 _er1/=nmc;
893 _er1*=(xH-xL);
894 //std::cout<<"_er1= "<<_er1<<std::endl;
895 Rad2Xs/=nmc; // the leading second order correction
896 double thexs= Rad2Xs*(xH-xL); //xH-xL is the Jacobi factor
897 //std::cout<<"thexs= "<<thexs<<std::endl;
898 return thexs;
899}
900
901
902
903double EvtConExc::gamHXSection(double El,double Eh){// using Gaussian integration subroutine from Cern lib
904 std::cout<<" "<<std::endl;
905 extern double Rad2difXs(double *x);
906 extern double Rad2difXs2(double *x);
907 double eps = 0.01;
908 double Rad2Xs;
909 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
910 if(Rad2Xs==0){
911 for(int iter=0;iter<10;iter++){
912 eps += 0.01;
913 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
914 if(!Rad2Xs) return Rad2Xs;
915 }
916 }
917 return Rad2Xs; // the leading second order correction
918}
919
920double EvtConExc::gamHXSection_er(double El,double Eh){// using Gaussian integration subroutine from Cern lib
921 std::cout<<" "<<std::endl;
922 extern double Rad2difXs_er(double *x);
923 extern double Rad2difXs_er2(double *x);
924 double eps = 0.01;
925 double Rad2Xs;
926 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
927 if(Rad2Xs==0){
928 for(int iter=0;iter<10;iter++){
929 eps += 0.01;
930 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
931 if(!Rad2Xs) return Rad2Xs;;
932 }
933 }
934 return Rad2Xs; // the leading second order correction
935}
936
937
938void EvtConExc::findMaxXS(EvtParticle *p){
939 //std::cout<<"nmc= "<<nmc<<std::endl;
940 gamId = EvtPDL::getId(std::string("gamma"));
941 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
942 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
943 double totxs = 0;
944 maxXS=-1;
945 _er1=0;
946 Rad2Xs =0;
947 int nmc = 50000;
948 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
949 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
950 int gamdaugs = _ndaugs+1;
951
952 EvtId theDaugs[10];
953 for(int i=0;i<_ndaugs;i++){
954 theDaugs[i] = daugs[i];
955 }
956 theDaugs[_ndaugs]=gamId;
957
958 gamH->makeDaughters(gamdaugs,theDaugs);
959
960 for(int i=0;i<gamdaugs;i++){
961 EvtParticle* di=gamH->getDaug(i);
962 double xmi=EvtPDL::getMass(di->getId());
963 di->setMass(xmi);
964 }
965
966 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
967 double thexs = difgamXs(gamH); //prepare the photon angular distribution
968 if(thexs>maxXS) {maxXS=thexs;}
969 }
970
971}
972
973double EvtConExc::difgamXs(EvtParticle *p){// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
974 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
975 EvtVector4R p0 = p->getDaug(0)->getP4();
976 for(int i=1;i<_ndaugs;i++){
977 p0 += p->getDaug(i)->getP4();
978 }
979 EvtVector4R pgam = p->getDaug(_ndaugs)->getP4();
980 double mhs = p0.mass();
981 double egam = pgam.get(0);
982 double sin2theta = pgam.get(3)/pgam.d3mag();
983 sin2theta = 1-sin2theta*sin2theta;
984
985 double cms = p->getP4().mass();
986 double alpha = 1./137.;
987 double pi = 3.1415926;
988 double x = 2*egam/cms;
989 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
990 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
991
992 double xs = myxsection->getXsection(mhs);
993 double difxs = 2*mhs/cms/cms * wsx *xs;
994 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
995 return difxs;
996}
997
998
999bool EvtConExc::gam_sampling(EvtParticle *p){
1000 double pm= EvtRandom::Flat(0.,1);
1001 double xs = difgamXs( p );
1002 double xsratio = xs/maxXS;
1003 if(pm<xsratio){return true;}else{return false;}
1004}
1005
1006bool EvtConExc::xs_sampling(double xs){
1007 double pm= EvtRandom::Flat(0.,1);
1008 // std::cout<<"Rdm= "<<pm<<std::endl;
1009 if(pm <xs/XS_max) {return true;} else {return false;}
1010}
1011bool EvtConExc::baryon_sampling(EvtVector4R pcm, EvtVector4R pi){
1012 EvtHelSys angles(pcm,pi); //using helicity sys.angles
1013 double theta=angles.getHelAng(1);
1014 double phi =angles.getHelAng(2);
1015 double costheta=cos(theta); //using helicity angles in parent system
1016
1017 // double costh = pcm.dot(pi)/pcm.d3mag()/pi.d3mag();
1018 //std::cout<<"costheta: "<<costheta<<", "<<costh<<std::endl;
1019
1020 double pm= EvtRandom::Flat(0.,1);
1021 double ang = 1 + costheta*costheta;
1022 if(pm< ang/2.) {return true;}else{return false;}
1023}
1024
1025bool EvtConExc::meson_sampling(EvtVector4R pcm, EvtVector4R pi){
1026 EvtHelSys angles(pcm,pi); //using helicity sys.angles
1027 double theta=angles.getHelAng(1);
1028 double phi =angles.getHelAng(2);
1029 double costheta=cos(theta); //using helicity angles in parent system
1030
1031 double pm= EvtRandom::Flat(0.,1);
1032 double ang = 1 - costheta*costheta;
1033 if(pm< ang/1.) {return true;}else{return false;}
1034}
1035
1036double EvtConExc::Rad1(double s, double x){ //first order correction
1037 //radiator correction upto second order, see Int. J. of Moder.Phys. A
1038 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
1039 double alpha = 1./137.;
1040 double pi=3.1415926;
1041 double me = 0.5*0.001; //e mass
1042 double sme = sqrt(s)/me;
1043 double L = 2*log(sme);
1044 double wsx = 2*alpha/pi/x*(L-1)*(1 - x + x*x/2 );
1045 return wsx;
1046}
1047
1048double EvtConExc::Rad2(double s, double x){
1049 //radiator correction upto second order, see Int. J. of Moder.Phys. A
1050 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
1051 double alpha = 1./137.;
1052 double pi=3.1415926;
1053 double me = 0.5*0.001; //e mass
1054 double xi2 = 1.64493407;
1055 double xi3=1.2020569;
1056 double sme = sqrt(s)/me;
1057 double L = 2*log(sme);
1058 double beta = 2*alpha/pi*(L-1);
1059 double delta2 = (9./8.-2*xi2)*L*L - (45./16.-5.5*xi2-3*xi3)*L-6./5.*xi2*xi2-4.5*xi3-6*xi2*log(2.)+3./8.*xi2+57./12.;
1060 double ap = alpha/pi;
1061 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
1062 double wsx = Delta * beta * pow(x,beta-1)-0.5*beta*(2-x);
1063 double wsx2 = (2-x)*( 3*log(1-x)-4*log(x) ) -4*log(1-x)/x-6+x;
1064 wsx = wsx + beta*beta/8. *wsx2;
1065 return wsx;
1066}
1067
1068
1069
1070double EvtConExc::Rad2difXs(EvtParticle *p){// leading second order xsection
1071 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
1072 double summass = p->getDaug(0)->getP4().mass();
1073 EvtVector4R v4p=p->getDaug(0)->getP4();
1074 for(int i=1;i<_ndaugs;i++){
1075 summass += p->getDaug(i)->getP4().mass();
1076 v4p += p->getDaug(i)->getP4();
1077 }
1078
1079 double Egam = p->getDaug(_ndaugs)->getP4().d3mag();
1080 double cms = p->getP4().mass();
1081 double mrg = cms - summass;
1082 double pm= EvtRandom::Flat(0.,1);
1083 double mhs = pm*mrg + summass;
1084
1085 double s = cms * cms;
1086 double x = 2*Egam/cms;
1087 //double mhs = sqrt( s*(1-x) );
1088 double wsx = Rad2(s,x);
1089
1090 // std::cout<<"random : "<<pm<<std::endl;
1091
1092 double xs = myxsection->getXsection(mhs);
1093 if(xs>XS_max){XS_max = xs;}
1094 double difxs = 2*mhs/cms/cms * wsx *xs;
1095 differ2 = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
1096 differ *= mrg; //Jacobi factor for variable m_hds
1097 difxs *= mrg;
1098 return difxs;
1099}
1100
1101double EvtConExc::Rad2difXs(double s, double x ){// leading second order xsection
1102 double wsx = Rad2(s,x);
1103 double mhs = sqrt(s*(1-x));
1104 double xs = myxsection->getXsection(mhs);
1105 //std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
1106 if(xs>XS_max){XS_max = xs;}
1107 double difxs = wsx *xs;
1108 differ2 = wsx *(myxsection->getErr(mhs));
1109 return difxs;
1110}
1111
1112
1113double EvtConExc::Rad1difXs(EvtParticle *p){// // first order xsection
1114 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
1115 double summass = p->getDaug(0)->getP4().mass();
1116 for(int i=1;i<_ndaugs;i++){
1117 summass += p->getDaug(i)->getP4().mass();
1118 }
1119
1120 double cms = p->getP4().mass();
1121 double mrg = cms - summass;
1122 double pm= EvtRandom::Flat(0.,1);
1123 double mhs = pm*mrg + summass;
1124
1125 double s = cms * cms;
1126 double x = 1 - mhs*mhs/s;
1127 double wsx = Rad1(s,x);
1128
1129 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
1130
1131 double xs = myxsection->getXsection(mhs);
1132 if(xs>XS_max){XS_max = xs;}
1133 double difxs = 2*mhs/cms/cms * wsx *xs;
1134
1135 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
1136 differ *= mrg; //Jacobi factor for variable m_hds
1137 difxs *= mrg;
1138 return difxs;
1139}
1140
1141void EvtConExc::init_Br_ee(){
1142 double psipp_ee=9.6E-06;
1143 double psip_ee =7.73E-03;
1144 double jsi_ee =5.94*0.01;
1145 double phi_ee =2.954E-04;
1146 double omega_ee=7.28E-05;
1147 double rho_ee = 4.72E-05;
1148 EvtId psppId=EvtPDL::getId(std::string("psi(3770)"));
1149 EvtId pspId=EvtPDL::getId(std::string("psi(2S)"));
1150 EvtId jsiId=EvtPDL::getId(std::string("J/psi"));
1151 EvtId phiId=EvtPDL::getId(std::string("phi"));
1152 EvtId omegaId=EvtPDL::getId(std::string("omega"));
1153 EvtId rhoId=EvtPDL::getId(std::string("rho0"));
1154 BR_ee.clear();
1155 ResId.clear();
1156
1157 BR_ee.push_back(rho_ee);
1158 BR_ee.push_back(omega_ee);
1159 BR_ee.push_back(phi_ee);
1160 BR_ee.push_back(jsi_ee);
1161 BR_ee.push_back(psip_ee);
1162 BR_ee.push_back(psipp_ee);
1163
1164 ResId.push_back(rhoId);
1165 ResId.push_back(omegaId);
1166 ResId.push_back(phiId);
1167 ResId.push_back(jsiId);
1168 ResId.push_back(pspId);
1169 ResId.push_back(psppId);
1170
1171}
1172
1173double EvtConExc::Ros_xs(double mx,double bree, EvtId pid){// in unit nb
1174 double pi=3.1415926;
1175 double s= mx*mx;
1176 double width = EvtPDL::getWidth(pid);
1177 double mass = EvtPDL::getMeanMass(pid);
1178 double xv = 1-mass*mass/s;
1179 double sigma = 12*pi*pi*bree*width/mass/s;
1180 sigma *= Rad2(s, xv);
1181 double nbar = 4E05; // ! GeV-2 = 4*10^5 nbar
1182 return sigma*nbar;
1183}
1184
1185
1186double Rad2difXs(double *mhs){// leading second order xsection, mhs: the mass of final states
1187 double cms = EvtConExc::_cms;
1188 double s = cms * cms;
1189 double x = 1 - (*mhs)*(*mhs)/s;
1190 EvtConExc myconexc;
1191 double wsx;
1192 double dhs=(*mhs);
1193 double xs = EvtConExc::myxsection->getXsection(dhs);
1194 wsx=myconexc.Rad2(s,x);
1195 if(xs>EvtConExc::XS_max){EvtConExc::XS_max = xs;}
1196 double difxs = 2*dhs/cms/cms * wsx *xs;
1197 std::ofstream oa;oa<<x<<std::endl; //without this line, the running will breaking, I don't know why
1198 return difxs;
1199}
1200double Rad2difXs_er(double *mhs){// leading second order xsection, mhs: the mass of final states
1201 double cms = EvtConExc::_cms;
1202 double s = cms * cms;
1203 double x = 1 - (*mhs)*(*mhs)/s;
1204 EvtConExc myconexc;
1205 double wsx;
1206 double xs = EvtConExc::myxsection->getErr(*mhs);
1207 wsx=myconexc.Rad2(s,x);
1208 double differ2 = 2*(*mhs)/cms/cms * wsx *xs;
1209 std::ofstream oa;oa<<x<<std::endl;
1210 return differ2;
1211}
1212
1213double Rad2difXs2(double *mhs){// leading second order xsection, mhs: the mass of final states
1214 double cms = EvtConExc::_cms;
1215 double s = cms * cms;
1216 double x = 1 - (*mhs)*(*mhs)/s;
1217 EvtConExc myconexc;
1218 double wsx;
1219 double dhs=(*mhs);
1220 double xs = EvtConExc::myxsection->getXsection(dhs);
1221 wsx=myconexc.Rad2(s,x);
1222 if(xs>EvtConExc::XS_max){EvtConExc::XS_max = xs;}
1223 double difxs = 2*dhs/cms/cms * wsx *xs;
1224 oa<<x<<std::endl; //without this line, the running will breaking, I don't know why
1225 return difxs;
1226}
1227
1228double Rad2difXs_er2(double *mhs){// leading second order xsection, mhs: the mass of final states
1229 double cms = EvtConExc::_cms;
1230 double s = cms * cms;
1231 double x = 1 - (*mhs)*(*mhs)/s;
1232 EvtConExc myconexc;
1233 double wsx;
1234 double xs = EvtConExc::myxsection->getErr(*mhs);
1235 wsx=myconexc.Rad2(s,x);
1236 double differ2 = 2*(*mhs)/cms/cms * wsx *xs;
1237 oa<<x<<std::endl;
1238 return differ2;
1239}
1240
1241
1242double EvtConExc::SoftPhoton_xs(double s, double b){
1243 double alpha = 1./137.;
1244 double pi=3.1415926;
1245 double me = 0.5*0.001; //e mass
1246 double xi2 = 1.64493407;
1247 double xi3=1.2020569;
1248 double sme = sqrt(s)/me;
1249 double L = 2*log(sme);
1250 double beta = 2*alpha/pi*(L-1);
1251 double delta2 = (9./8.-2*xi2)*L*L - (45./16.-5.5*xi2-3*xi3)*L-6./5.*xi2*xi2-4.5*xi3-6*xi2*log(2.)+3./8.*xi2+57./12.;
1252 double ap = alpha/pi;
1253 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
1254
1255 double beta2 = beta*beta;
1256 double b2 = b*b;
1257
1258 double xs=(-32*b*beta + 8*pow(b,2)*beta - 10*b*pow(beta,2) + pow(b,2)*pow(beta,2) + 32*pow(b,beta)*Delta -
1259 6*(3 - 4*b + pow(b,2))*pow(beta,2)*log(1 - b) - 32*b*pow(beta,2)*log(b) + 8*pow(b,2)*pow(beta,2)*log(b) +
1260 16*pow(beta,2)*Li2(b))/32. ;
1261 return xs;
1262
1263}
1264
1265double EvtConExc::Li2(double x){
1266 double li2= -x +x*x/4. - x*x*x/9;
1267 return li2;
1268}
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
mass
double mass
Definition CosmicGenerator.cxx:138
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
Rad2difXs2
double Rad2difXs2(double *mhs)
Definition EvtConExc.cc:1213
oa
std::ofstream oa
Definition EvtConExc.cc:133
Rad2difXs
double Rad2difXs(double *m)
Definition EvtConExc.cc:1186
Rad2difXs_er2
double Rad2difXs_er2(double *mhs)
Definition EvtConExc.cc:1228
dgauss_
double dgauss_(double(*fun)(double *), double *, double *, double *)
Rad2difXs_er
double Rad2difXs_er(double *m)
Definition EvtConExc.cc:1200
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition EvtDiracSpinor.cc:75
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition EvtStreamInputIterator.hh:101
eps
EvtTensor3C eps(const EvtVector3R &v)
Definition EvtTensor3C.cc:307
alpha
const double alpha
Definition FastVertexFit.cxx:4
s
XmlRpcServer s
Definition HelloServer.cpp:11
weight
*********Class see also m_nmax DOUBLE PRECISION m_MasPhot DOUBLE PRECISION m_phsu DOUBLE PRECISION m_Xenph DOUBLE PRECISION m_r2 DOUBLE PRECISION m_WtMass INTEGER m_nmax INTEGER m_Nevgen INTEGER m_IsFSR INTEGER m_MarTot *COMMON c_KarFin $ !Output file $ !Event serial number $ !alpha QED at Thomson limit $ !minimum energy at CMS for remooval $ !infrared dimensionless $ !dummy photon IR regulator $ !crude photon multiplicity enhancement factor *EVENT $ !MC crude volume of PhhSpace *Sfactors $ !YFS formfactor IR part only $ !YFS formfactor non IR finite part $ !mass weight
Definition KarFin.h:34
me
const double me
Definition PipiJpsi.cxx:48
EvtConExc::Li2
double Li2(double x)
Definition EvtConExc.cc:1265
EvtConExc::findMaxXS
void findMaxXS(EvtParticle *p)
Definition EvtConExc.cc:938
EvtConExc::getMode
static int getMode
Definition EvtConExc.hh:128
EvtConExc::init
void init()
Definition EvtConExc.cc:158
EvtConExc::init_Br_ee
void init_Br_ee()
Definition EvtConExc.cc:1141
EvtConExc::gamHXSection_er
double gamHXSection_er(double El, double Eh)
Definition EvtConExc.cc:920
EvtConExc::Ros_xs
double Ros_xs(double mx, double bree, EvtId pid)
Definition EvtConExc.cc:1173
EvtConExc::Rad1
double Rad1(double s, double x)
Definition EvtConExc.cc:1036
EvtConExc::myxsection
static EvtXsection * myxsection
Definition EvtConExc.hh:124
EvtConExc::meson_sampling
bool meson_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition EvtConExc.cc:1025
EvtConExc::xs_sampling
bool xs_sampling(double xs)
Definition EvtConExc.cc:1006
EvtConExc::gamHXSection
double gamHXSection(EvtParticle *p, double El, double Eh, int nmc=100000)
Definition EvtConExc.cc:819
EvtConExc::Rad1difXs
double Rad1difXs(EvtParticle *p)
Definition EvtConExc.cc:1113
EvtConExc::init_mode
void init_mode(int mode)
Definition EvtConExc.cc:337
EvtConExc::baryon_sampling
bool baryon_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition EvtConExc.cc:1011
EvtConExc::Rad2difXs
double Rad2difXs(EvtParticle *p)
Definition EvtConExc.cc:1070
EvtConExc::SoftPhoton_xs
double SoftPhoton_xs(double s, double b)
Definition EvtConExc.cc:1242
EvtConExc::difgamXs
double difgamXs(EvtParticle *p)
Definition EvtConExc.cc:973
EvtConExc::initProbMax
void initProbMax()
Definition EvtConExc.cc:711
EvtConExc::_cms
static double _cms
Definition EvtConExc.hh:125
EvtConExc::~EvtConExc
virtual ~EvtConExc()
Definition EvtConExc.cc:136
EvtConExc::gam_sampling
bool gam_sampling(EvtParticle *p)
Definition EvtConExc.cc:999
EvtConExc::Rad2
double Rad2(double s, double x)
Definition EvtConExc.cc:1048
EvtConExc::getName
void getName(std::string &name)
Definition EvtConExc.cc:145
EvtConExc::decay
void decay(EvtParticle *p)
Definition EvtConExc.cc:717
EvtConExc::XS_max
static double XS_max
Definition EvtConExc.hh:126
EvtConExc::clone
EvtDecayBase * clone()
Definition EvtConExc.cc:151
EvtDecayBase::getArg
double getArg(int j)
Definition EvtDecayBase.cc:564
EvtDecayBase::getDaug
EvtId getDaug(int i)
Definition EvtDecayBase.hh:66
EvtHelSys::getHelAng
double getHelAng(int i)
Definition EvtHelSys.cc:54
EvtId
Definition EvtId.hh:27
EvtPDL::getWidth
static double getWidth(EvtId i)
Definition EvtPDL.hh:54
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition EvtPDL.hh:45
EvtPDL::evtIdFromStdHep
static EvtId evtIdFromStdHep(int stdhep)
Definition EvtPDL.cc:244
EvtPDL::name
static std::string name(EvtId i)
Definition EvtPDL.hh:64
EvtPDL::getMass
static double getMass(EvtId i)
Definition EvtPDL.hh:46
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition EvtPDL.cc:287
EvtParticleFactory::particleFactory
static EvtParticle * particleFactory(EvtSpinType::spintype spinType)
Definition EvtParticleFactory.cc:46
EvtParticle::setMass
void setMass(double m)
Definition EvtParticle.hh:358
EvtParticle::makeDaughters
void makeDaughters(int ndaug, EvtId *id)
Definition EvtParticle.cc:1175
EvtParticle::getId
EvtId getId() const
Definition EvtParticle.cc:113
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition EvtParticle.cc:121
EvtParticle::printTree
void printTree() const
Definition EvtParticle.cc:895
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition EvtParticle.cc:85
EvtParticle::deleteDaughters
void deleteDaughters(bool keepChannel=false)
Definition EvtParticle.cc:538
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition EvtParticle.cc:1069
EvtParticle::deleteTree
void deleteTree()
Definition EvtParticle.cc:555
EvtRandom::Flat
static double Flat()
Definition EvtRandom.cc:73
EvtVector4R::mass
double mass() const
Definition EvtVector4R.cc:39
EvtVector4R::get
double get(int i) const
Definition EvtVector4R.hh:179
EvtVector4R::d3mag
double d3mag() const
Definition EvtVector4R.cc:186
EvtVector4R::mass2
double mass2() const
Definition EvtVector4R.hh:116
EvtVector4R::set
void set(int i, double d)
Definition EvtVector4R.hh:183
EvtXsection::getXlw
double getXlw()
Definition EvtXsection.hh:150
EvtXsection::getErr
double getErr(double mx)
Definition EvtXsection.cc:11076
EvtXsection::getMsg
std::string getMsg()
Definition EvtXsection.hh:151
EvtXsection::getEr
std::vector< double > getEr()
Definition EvtXsection.hh:148
EvtXsection::getXup
double getXup()
Definition EvtXsection.hh:149
EvtXsection::getUnit
std::string getUnit()
Definition EvtXsection.hh:144
EvtXsection::getXsection
double getXsection(double mx)
Definition EvtXsection.cc:11070
EvtXsection::getYY
std::vector< double > getYY()
Definition EvtXsection.hh:147
EvtXsection::setBW
void setBW(int pdg)
Definition EvtXsection.cc:11170
EvtXsection::getXX
std::vector< double > getXX()
Definition EvtXsection.hh:146
pi
const float pi
Definition vector3.h:133