BOSS 7.0.4
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtConExc.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of models developed at BES collaboration
4// based on the EvtGen framework. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/BesCopyright
8// Copyright (A) 2006 Ping Rong-Gang @IHEP, Wang Weiping @USTC
9//
10// Module: EvtXsection.hh
11//
12// Description: To define cross section for the continuum exclusive process
13// Experimental cross section taken from PRD73,012005, PRD76,092006, also
14// see a review: Rev. Mod. Phys. 83,1545
15// Modification history:
16//
17// Ping R.-G. Nov., 2012 Module created
18// Wang Weiping: Nov., 2017 added the mode==85: e+e- --> D_s*+D_s*-
19// Wang Weiping: Nov., 2017 fixed final states of mode==28,33,92
20// Wang Weiping: Nov., 2017 fixed the angle sampling of VP states
21//
22 /*******************--- mode definition: also see EvtXsection.cc
23 mode==0: ppbar
24 mode==1: nnbar
25 mode==2: Lambda0 anti-Lambda0
26 mode==3: Sigma0 anti-Sigma0
27 mode==4: Lambda0 anti-Sigma0
28 mode==5: Sigma0 anti-Lambda0
29 mode==6: pi+ pi-
30 mode==7: pi+ pi- pi0
31 mode==8: K+K- pi0
32 mode==9: KsK+pi-
33 mode==10: KsK-pi+
34 mode==11: K+K-eta
35 mode==12: 2(pi+pi-)
36 mode==13: pi+pi-2pi0
37 mode==14: K+K-pi+pi-
38 mode==15: K+K-2pi0
39 mode==16: 2(K+K-)
40 mode==17: 2(pi+pi-)pi0
41 mode==18: 2(pi+pi-)eta
42 mode==19: K+K-pi+pi-pi0
43 mode==20: K+K-pi+pi-eta
44 mode==21: 3(pi+pi-)
45 mode==22: 2(pi+pi-pi0)
46 mode==23: phi eta
47 mode==24: phi pi0
48 mode==25: K+K*-
49 mode==26: K-K*+
50 mode==27: K_SK*0-bar
51 mode==28: K*0(892)K+pi-
52 mode==29: K*0(892)K-pi+
53 mode==30: K*+K-pi0
54 mode==31: K*-K+pi0
55 mode==32: K_2*(1430)0 K+pi-
56 mode==33: K_2*(1430)0 K-pi+
57 mode==34: K+K-rho
58 mode==35: phi pi+pi-
59 mode==36: phi f0(980)
60 mode==37: eta pi+pi-
61 mode==38: omega pi+ pi-
62 mode==39: omega f0(980)
63 mode==40: eta' pi+ pi-
64 mode==41: f_1(1285)pi+pi-
65 mode==42: omega K+K-
66 mode==43: omega pi+pi-pi0
67 mode==44: Sigma+ Sigma- (Sigma0 Sigma0-bar SU(3) extention )
68 mode==45: K+K-
69 mode==46: K_S K_L
70 mode==47: omega eta
71 mode==48: phi eta'
72 mode==49: phi K+K-
73 mode==50: ppbar pi0
74 mode==51: ppbar eta
75 mode==52: omega pi0
76
77 mode==70: D0 anti-D0
78 mode==71: D+D-
79 mode==72: D+D*-
80 mode==73: D-D*+
81 mode==74: D*+D*-
82 mode==68: D0 anti-D*0
83 mode==69: anti-D0 D*0
84 mode==67: D*0 anti-D*0
85
86 mode==59: D0 anti-D0 pi0
87 mode==63: D+D-pi0
88 mode==75: D0D-pi+
89 mode==76: anti-D0 D+pi-
90
91 mode==77: D0D*-pi+
92 mode==78: anti-D0 D*+pi-
93 mode==65: D-D*0pi+
94 mode==66: D+ anti-D*0pi-
95 mode==60: anti-D0 D*0pi0
96 mode==61: D0 anti-D*0pi0
97 mode==62: D+D*-pi0
98 mode==64: D-D*+pi0
99
100 mode==90: pi+pi- J/psi
101 mode==92: pi0pi0 J/psi
102 mode==91: pi+pi- psi(2S)
103 mode==79: pi0pi0 psi(2S)
104 mode==80: eta J/psi
105 mode==81: pi+pi- h_c
106 mode==82: pi0pi0 h_c
107 mode==83: K+K- J/psi
108 mode==84: KsKs J/psi
109
110 mode==93: D_s+ D_s-
111 mode==94: D_s^{*+}D_s^-
112 mode==95: D_s^{*-}D_s^+
113 mode==85: D_s^{*+}D_s^{*-}
114
115 mode==96: Lambda_c+ anti-Lamba_c-
116 mode==111: e+e- ->Zc(3900)+ pi-
117 mode==112: e+e- ->Zc(3900)- pi+
118//---
119 For cross section, see Ref. Rev.Mod. Phys.83,1545
120 *************************************/
121//------------------------------------------------------------------------
122#include <math.h>
123#include "EvtGenBase/EvtPatches.hh"
124#include <stdlib.h>
125#include "EvtGenBase/EvtParticle.hh"
126#include "EvtGenBase/EvtGenKine.hh"
127#include "EvtGenBase/EvtPDL.hh"
128#include "EvtGenBase/EvtReport.hh"
129#include "EvtGenModels/EvtPhsp.hh"
130#include "EvtGenModels/EvtConExc.hh"
131#include "EvtGenBase/EvtVector4R.hh"
132#include "EvtGenBase/EvtParticleFactory.hh"
133#include "EvtGenBase/EvtRandom.hh"
134#include "EvtGenBase/EvtHelSys.hh"
135#include "EvtGenBase/EvtDecayTable.hh"
136#include "EvtGenBase/EvtDecayBase.hh"
137#include "EvtGenBase/EvtPDL.hh"
138#include "EvtGenModels/EvtGlobalSet.hh"
139#include <string>
140#include <iostream>
141#include <fstream>
142#include <istream>
143#include <strstream>
144#include <stdio.h>
145#include <fstream>
146#include <strstream>
147#include "TGraphErrors.h"
148#include "TCanvas.h"
149#include "TPostScript.h"
150#include "TStyle.h"
151#include "TMultiGraph.h"
152using namespace std;
153
154int EvtConExc::nconexcdecays=0;
155EvtDecayBasePtr* EvtConExc::conexcdecays=0;
156int EvtConExc::ntable=0;
157
158int EvtConExc::ncommand=0;
159int EvtConExc::lcommand=0;
160std::string* EvtConExc::commands=0;
161double EvtConExc::AF[600];
162double EvtConExc::AA[600];
163double EvtConExc::MH[600];
164
165double EvtConExc::mlow=0;
166double EvtConExc::mup=0;
167
168extern "C" {
169 extern double dgauss_( double (*fun)(double*), double *,double *, double *);
170}
171
172extern "C" {
173 extern double divdif_( float*, float *,int *, float *,int*);
174}
175
176extern "C" {
177 extern void polint_( float*, float *,int *, float *,float*);
178}
179
180extern "C" {
181 extern void hadr5n12_( float*, float *,float *, float *,float *, float *);
182}
183
184
185double Rad2difXs(double *m);
186double Rad2difXs_er(double *m);
187
188EvtXsection *EvtConExc::myxsection;
189EvtXsection *EvtConExc::staxsection;
190double EvtConExc::_cms;
191double EvtConExc::XS_max;
192
193double EvtConExc::_xs0=0;
194double EvtConExc::_xs1=0;
195double EvtConExc::_er0=0;
196double EvtConExc::_er1=0;
197int EvtConExc::_nevt=0;
198
199std::ofstream oa;
200int EvtConExc::getMode;
201int EvtConExc::conexcmode=-1;
202int EvtConExc::multimode=0;
203std::vector<std::vector <double> > EvtConExc::VXS;
204
205EvtConExc::~EvtConExc() {
206// if(_mode==74110)checkEvtRatio();
207 if(myFisr.size())OutStaISR();
208 if(mydbg){
209 // xs->Write();
210 myfile->Write();
211 }
212
213 delete myxsection;
214 delete staxsection;
215 gamH->deleteTree();
216
217 //the deletion of commands is really uggly!
218
219 int i;
220 if (nconexcdecays==0) {
221 delete [] commands;
222 commands=0;
223 return;
224 }
225
226 for(i=0;i<nconexcdecays;i++){
227 if (conexcdecays[i]==this){
228 conexcdecays[i]=conexcdecays[nconexcdecays-1];
229 nconexcdecays--;
230 if (nconexcdecays==0) {
231 delete [] commands;
232 commands=0;
233 }
234 return;
235 }
236 }
237
238
239}
240
241void EvtConExc::getName(std::string& model_name){
242
243 model_name="ConExc";
244
245}
246
247EvtDecayBase* EvtConExc::clone(){
248
249 return new EvtConExc;
250
251}
252
253
254void EvtConExc::init(){
255 mydbg=false;
256 myFisr.clear();
257 ReadVP();
258 VXS.resize(120);
259 for(int i=0;i<120;i++){
260 VXS[i].resize(600);
261 }
262 _mode = getArg(0);
263 for(int i=0;i<97;i++){
264 if(53<=i && i<=58) continue;
265 if(86<=i && i<=89) continue;
266 if(_mode==74110 ||_mode==-100) vmode.push_back(i);
267 }
268 if(_mode==74110||_mode==-100){
269 vmode.push_back(74100);
270 vmode.push_back(74101);
271 vmode.push_back(74102);
272 vmode.push_back(74103);
273 vmode.push_back(74104);
274 vmode.push_back(74105);
275 vmode.push_back(74110);
276 }
277
278 std::cout<<"==== Parameters set by users====="<<std::endl;
279 for(int i=0;i<ncommand;i++){
280 std::cout<<commands[i]<<std::endl;
281 }
282 std::cout<<"==================================="<<std::endl;
283 InitPars();
284 std::cout<<threshold<<" "<<beamEnergySpread<<std::endl;
285 //----------------
286 // check that there are 0 arguments
287 // checkNArg(1);
288 //Vector ISR process
289 VISR=0;
290 if(getNDaug()==2){
291 if(getDaugs()[0]==EvtPDL::getId("gamma")&&getDaugs()[1]==EvtPDL::getId("gamma*") || getDaugs()[0]==EvtPDL::getId("gamma*")&&getDaugs()[1]==EvtPDL::getId("gamma")) VISR=1;
292 }else if(getNDaug()>2){std::cout<<"Bad daughter specified"<<std::endl;abort();
293 }
294
295 // cmsspread=0.0015; //CMC energy spread
296 testflag=0;
297 getResMass();
298 if(getArg(0) == -1){
299 radflag=0;
300 mydbg=false;
301 _mode = getArg(0);
302 pdgcode = getArg(1);
303 pid = EvtPDL::evtIdFromStdHep(pdgcode );
304 nson = getNArg()-2;
305 std::cout<<"The decay daughters are "<<std::endl;
306 for(int i=2;i<getNArg();i++){son[i-2]=EvtPDL::evtIdFromStdHep(getArg(i));std::cout<<EvtPDL::name(son[i-2])<<" ";}
307 std::cout<<std::endl;
308 }else if(getArg(0)==-2){
309 radflag=0;mydbg=false;
310 _mode = getArg(0);
311 nson = getNArg()-1;
312 for(int i=1;i<getNArg();i++){son[i-1]=EvtPDL::evtIdFromStdHep(getArg(i));}
313 }else if(getArg(0)==-100){
314 _mode = getArg(0);
315 multimode=1;
316 vmd.clear();
317 for(int i=1;i<getNArg();i++){vmd.push_back(getArg(i));}
318 std::cout<<" multi-exclusive mode "<<std::endl;
319 for(int i=1;i<getNArg();i++){std::cout<<getArg(i)<<" ";}
320 std::cout<<std::endl;
321 if(vmd[vmd.size()-1]==74110){std::cout<<"74110 is not allowd for multimode simulation"<<std::endl;abort();}
322 }else if(getNArg()==1){ _mode = getArg(0); radflag=0;mydbg=false;
323 }else if(getNArg()==2){ _mode = getArg(0); radflag=getArg(1);mydbg=false;
324 }else if(getNArg()==3){ _mode = getArg(0); radflag=getArg(1);mydbg=getArg(2);
325 }else{std::cout<<"ConExc:number of parameter should be 1,2 or 3, but you set "<<getNArg()<<std::endl;::abort(); }
326 //--- debugging
327 //std::cout<<_mode<<" "<<pdgcode<<" "<<nson<<" "<<EvtPDL::name(son[0])<<" "<<EvtPDL::name(son[1])<<std::endl;
328
329 gamId = EvtPDL::getId(std::string("gamma"));
330 init_mode(_mode);
331 XS_max=-1;
332 //-----debugging to make root file
333 mydbg = false; //designer need to remove this line
334 if(mydbg){
335 myfile = new TFile("mydbg.root","RECREATE");
336 xs=new TTree ("xs","momentum of rad. gamma and hadrons");
337 xs->Branch("imode",&imode,"imode/I");
338 xs->Branch("pgam",pgam,"pgam[4]/D");
339 xs->Branch("phds",phds,"phds[4]/D");
340 xs->Branch("ph1",ph1,"ph1[4]/D");
341 xs->Branch("ph2",ph2,"ph2[4]/D");
342 xs->Branch("mhds",&mhds,"mhds/D");
343 xs->Branch("mass1",&mass1,"mass1/D");
344 xs->Branch("mass2",&mass2,"mass2/D");
345 xs->Branch("costheta",&costheta,"costheta/D");
346 xs->Branch("cosp",&cosp,"cosp/D");
347 xs->Branch("cosk",&cosk,"cosk/D");
348 xs->Branch("sumxs",&sumxs,"sumxs/D");
349 xs->Branch("selectmode",&selectmode,"selectmode/D");
350 }
351 //--- prepare the output information
352 EvtId parentId =EvtPDL::getId(std::string("vpho"));
353 EvtPDL::reSetWidth(parentId,0.0);
354 double parentMass = EvtPDL::getMass(parentId);
355 std::cout<<"parent mass = "<<parentMass<<std::endl;
356
357
358 EvtVector4R p_init(parentMass,0.0,0.0,0.0);
359 EvtParticle *p=EvtParticleFactory::particleFactory(parentId,p_init);
360 theparent = p;
361 myxsection = new EvtXsection (_mode);
362 staxsection = new EvtXsection (_mode);
363 if(_mode==-100) {
364 myxsection->setModes(vmd);
365 myxsection->ini_data_multimode();
366 staxsection->setModes(vmd);
367 staxsection->ini_data_multimode();
368 }
369 double mth=0;
370 double mx = EvtPDL::getMeanMass(parentId); //p->getP4().mass();
371 if(_mode ==-1){
372 myxsection->setBW(pdgcode);
373 staxsection->setBW(pdgcode);
374 for(int i=0;i<nson;i++){mth +=EvtPDL::getMeanMass(son[i]);}
375 if(mx<mth){std::cout<<"The CMS energy is lower than the threshold of the final states"<<std::endl;abort();}
376 }else if(_mode == -2){
377 mth=myxsection->getXlw();
378 }else{ mth=myxsection->getXlw();}
379 _cms = mx;
380 _unit = myxsection->getUnit();
381
382 std::cout<<"The specified mode is "<<_mode<<std::endl;
383 EvtConExc::getMode = _mode;
384 //std::cout<<"getXlw= "<<mth<<std::endl;
385
386 Egamcut = 0.001; //set photon energy cut according to the BESIII detector
387 double Esig = 0.0001; //sigularity energy
388 double x = 2*Egamcut/_cms;
389 double s = _cms*_cms;
390 double mhscut = sqrt(s*(1-x));
391
392 //for vaccum polarization
393 float fe,fst2,fder,ferrder,fdeg,ferrdeg;
394 fe=_cms;
395 //using the vacuum pol. value as given by http://www-com.physik.hu-berlin.de/~fjeger/software.html
396 vph=getVP(_cms);
397 if(3.0943<_cms && _cms<3.102) vph=1;// For J/psi, the vacuum factor is included in the resonance
398 std::cout<<"sum of leptonic, hadronic contributions(u,d,s,c,b) to vacuum polarization: "<<vph<<" @"<<fe<<"GeV"<<std::endl;
399
400 //caculate the xs for ISR to psiprime, J/psi and phi narrow resonance.
401 double totxsIRS=0;
402 init_Br_ee();
403 double mthrld = EvtPDL::getMeanMass(daugs[0]); //threshold mass of hadrons
404 for(int jj=1;jj<_ndaugs;jj++){
405 mthrld += EvtPDL::getMeanMass(daugs[jj]);
406 }
407
408 ISRXS.clear();ISRM.clear();ISRFLAG.clear();
409 for(int ii=0;ii<3;ii++){
410 double mR = EvtPDL::getMeanMass(ResId[ii]);
411 if(mx<mR || _mode !=74110) continue;
412 double narRxs=Ros_xs(mx,BR_ee[ii],ResId[ii]);
413 std::cout<<"The corss section for gamma "<<EvtPDL::name(ResId[ii]).c_str()<<" is: "<< narRxs<<"*Br (nb)."<<std::endl;
414 ISRXS.push_back(narRxs);
415 ISRM.push_back(mR);
416 ISRFLAG.push_back(1.);
417 ISRID.push_back(ResId[ii]);
418 totxsIRS += narRxs;
419 }
420 std::cout<<"==========================================================================="<<std::endl;
421
422 //-- calculated with MC integration method
423 double mhdL=myxsection->getXlw();
424 if(mhdL<SetMthr) mhdL=SetMthr;
425std::cout<<"SetMthr= "<<SetMthr<<std::endl;
426 EgamH = (s-mhdL*mhdL)/2/sqrt(s);
427 int nmc=1000000;
428 _xs0 = gamHXSection(s,Esig,Egamcut,nmc);
429 _er0 = _er1; // stored when calculate _xs0
430 std::cout<<"_er0= "<<_er0<<std::endl;
431 _xs1 = gamHXSection(s,Egamcut,EgamH,nmc);
432 std::cout<<"_xs1= "<<_xs0<<std::endl;
433 double xs1_err = _er1;
434
435
436 //--- sigularity point x from 0 to 0.0001GeV
437 double xb= 2*Esig/_cms;
438 double sig_xs = SoftPhoton_xs(s,xb)*(myxsection->getXsection(mx));
439 _xs0 += sig_xs;
440 std::cout<<"_xs0= "<<_xs0<<std::endl;
441
442 //prepare the observed cross section with interpotation function divdif in CERNLIB
443 testflag=1;
444 int Nstp=598;
445 double stp=(EgamH - Egamcut)/Nstp;
446 for(int i=0;i<Nstp;i++){//points within Egam(max) ~ Egamcut=25MeV
447 double Eg0=EgamH - i*stp;
448 double Eg1=EgamH - (i+1)*stp;
449 int nmc=20000;
450 int nint=100;
451 //double xsi= gamHXSection(s,Eg1,Eg0,nmc);
452 double xsi= trapezoid(s,Eg1,Eg0,nint,myxsection);
453 AA[i]=(Eg1+Eg0)/2;
454 double mhi=sqrt(_cms*_cms-2*_cms*AA[i]);
455 double mh2=sqrt(_cms*_cms-2*_cms*Eg1);
456 double binwidth = mh2-mhi;
457 //if(_mode==74110) xsi += addNarrowRXS(mhi,binwidth);
458 if(i==0) AF[0]=xsi;
459 if(i>=1) AF[i]=AF[i-1]+xsi;
460 RadXS[i]=xsi/stp;
461 }
462 //add the interval 0~Esig
463 AA[598]=Egamcut; AA[599]=0; //AA[i]: E_gamma
464 AF[598]=AF[597];
465 AF[599]=AF[598]+ _xs0;
466 RadXS[599]=_xs0;
467
468
469 //prepare observed cross section for each mode
470std::cout<<"mhdL = "<<mhdL<<", SetMthr= "<<SetMthr<<", EgamH = "<<EgamH<<std::endl;
471 for(int i=0;i<vmode.size();i++){
472 if(_mode==74110||_mode==-100) mk_VXS(Esig,Egamcut,EgamH,i);
473 }
474 if(_mode==74110||_mode==-100) writeDecayTabel();
475 //after mk_VXS, restore the myxsection to _mode
476 if(_mode !=-100){
477 delete myxsection;
478 myxsection = new EvtXsection (_mode);
479 }
480 //debugging VXS
481 /*
482 double xtmp=0;
483 for(int i=0;i<vmode.size();i++){
484 xtmp+=VXS[i][599];
485 for(int j=0;j<600;j++){
486 std::cout<<VXS[i][j]<<" ";
487 }
488 std::cout<<std::endl;
489 }
490 */
491 //output observed cross section for the gamma Resonance mode
492 //std::cout<<"vmode.size======="<<vmode.size()<<std::endl;
493 for(int i=0;i<vmode.size();i++){
494 int md=vmode[i];
495 if(md<74100 || md>74106) continue;
496 std::string partname="";
497 if(md==74100) {partname="J/psi";}
498 else if(md==74101) {partname="psi(2S)";}
499 else if(md==74102) {partname="psi(3770)";}
500 else if(md==74103) {partname="phi";}
501 else if(md==74104) {partname="omega";}
502 else if(md==74105) {partname="rho0";}
503 else if(md==74106) {partname="rho(3S)0";}
504 else{;}
505 std::cout<<"The observed cross section for gamma "<<partname<<": "<<VXS[i][599]<<" nb"<<std::endl;
506 }
507 //--
508
509 for(int i=0;i<600;i++){
510 MH[i]=sqrt(_cms*_cms-2*_cms*AA[i]);
511 //std::cout<<i<<" Egamma "<<AA[i]<<" Mhadons "<<MH[i]<<std::endl;
512 }
513 //for debugging
514 //for(int i=0;i<600;i++){double s=_cms*_cms; double mhds=sqrt(s-2*_cms*AA[i]);std::cout<<"Mhds="<<mhds<<" Egam="<<AA[i]<<" "<<AF[i]<<std::endl;}
515 //std::cout<<"VXS[79][599]: "<<VXS[79][599]<<" VXS[79][598]= "<<VXS[79][598]<<std::endl;
516 std::cout<<"EgamH = "<<EgamH<<", obsvXS = "<<_xs0+_xs1<<", _xs1 = "<<_xs1 << ", _xs0 = "<<_xs0<<std::endl;
517 std::cout<<"EgamH = "<<EgamH<<", AF[599] = "<< AF[599]<< ", AF[598] = "<<AF[598]<<", _xs0 = "<<_xs0<<std::endl;
518 //double Etst=0.0241838;
519 //std::cout<<"Etst="<<Etst<<" "<<gamHXSection(s,Etst,EgamH,nmc)<<" "<< lgr(AA,AF,600,Etst)<<std::endl; abort();
520
521 //for information output
522 double bornXS = myxsection->getXsection(mx); // std::cout<<"BornXS= "<<bornXS<<std::endl;
523 double bornXS_er= myxsection->getErr(mx);
524 double obsvXS = _xs0 + _xs1; //gamHXSection(mth ,mx);
525 double obsvXS_er= _er0 + xs1_err;
526 double corr = obsvXS/bornXS;
527 double corr_er =corr*sqrt(bornXS_er*bornXS_er/bornXS/bornXS + obsvXS_er*obsvXS_er/obsvXS/obsvXS);
528
529
530 if(bornXS==0){std::cout<<"EvtConExc::init : The Born cross section at this point is zero!"<<std::endl;abort();}
531 if(obsvXS==0){std::cout<<"EvtConExc::init : The observed cross section at this point is zero!"<<std::endl;abort();}
532
533 //_xs0 += bornXS;// events with very soft photon (Egam<0.025) are taken as the born process
534 //_er0 = sqrt(_er0*_er0 + bornXS_er*bornXS_er);
535
536 std::cout<<""<<std::endl;
537 std::cout<<"========= Generation using cross section (Egamcut="<<Egamcut<<" GeV) =============="<<std::endl;
538 std::cout<<" sqrt(s)= "<<mx<< " GeV "<<std::endl;
539 if(_mode>=0 || _mode ==-2 || _mode ==-100){
540 std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<"+/-"<<_er0<<" in Unit "<<_unit.c_str()<<std::endl;
541 std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"+/-"<<xs1_err<<" in Unit "<<_unit.c_str()<<std::endl;
542 std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
543 std::cout<<"The radiative correction factor f_ISR*f_vacuum= sigma_obs/sigma_born(s) = "<<corr<<"+/-"<< fabs(corr_er)<<","<<std::endl;
544 std::cout<<"which is calcualted with these quantities:"<<std::endl;
545 std::cout<<"the observed cross section is "<<obsvXS<<"+/-"<< obsvXS_er<<_unit.c_str()<<std::endl;
546 std::cout<<"and the Born cross section (s) is "<<bornXS<<"+/-"<<bornXS_er<<_unit.c_str()<<std::endl;
547 std::cout<<"and the vaccum polarziation factor (lep. + hadr.) 1/|1-Pi|^2= "<<vph<<std::endl;
548 std::cout<<"Within the range "<<myxsection->getXlw()<<" ~ "<<myxsection->getXup()<<" GeV, "<<myxsection->getMsg().c_str()<<std::endl;
549 std::cout<<"==========================================================================="<<std::endl;
550 }else if(_mode==-1){
551 std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<" *Br_ee"<<" in Unit "<<_unit.c_str()<<std::endl;
552 std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"*Br_ee in Unit "<<_unit.c_str()<<std::endl;
553 std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
554 std::cout<<"The radiative correction factor sigma_obs/sigma_born(s) = (1+delta) = "<<corr<<"+/-"<< fabs(corr_er)<<std::endl;
555 std::cout<<"==========================================================================="<<std::endl;
556 }
557 std::cout<<std::endl;
558 std::cout<<std::endl;
559
560 findMaxXS(p);
561 _mhdL=myxsection->getXlw();
562 //--debugging
563 //std::cout<<"maxXS= "<<maxXS<<std::endl;
564
565 if(_xs0 == 0 && _xs1==0){std::cout<<"EvtConExc::ini() has zero cross section"<<std::endl;abort();}
566
567 std::cout<<std::endl;
568 std::cout<<std::endl;
569 mlow=myxsection->getXlw();
570 mup =myxsection->getXup();
571 //--- for debugging
572 if(mydbg && _mode!=74110){
573 int nd = myxsection->getYY().size();
574 double xx[10000],yy[10000],xer[10000],yer[10000];
575 for(int i=0;i<nd;i++){
576 xx[i]=myxsection->getXX()[i];
577 yy[i]=myxsection->getYY()[i];
578 yer[i]=myxsection->getEr()[i];
579 xer[i]=0.;
580 //std::cout<<"yy "<<yy[i]<<std::endl;
581 }
582 myth=new TH1F("myth","Exp.data",200,xx[0],xx[nd]);
583 for(int i=0;i<nd;i++){
584 myth->Fill(xx[i],yy[i]);
585 }
586 myth->SetError(yer);
587 myth->Write();
588
589 //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
590 }else if(mydbg && _mode==74110){
591 int nd = myxsection->getYY().size();
592 double xx[10000],yy[10000],xer[10000],yer[10000],ysum[10000],yysum[10000];
593 for(int i=0;i<nd;i++){
594 xx[i]=myxsection->getXX()[i];
595 yy[i]=myxsection->getYY()[i];
596 yer[i]=myxsection->getEr()[i];
597 xer[i]=0.;
598 std::cout<<"yy "<<yy[i]<<std::endl;
599 }
600
601#include "sty.C"
602 double xhigh=5.0;
603 double xlow = myxsection->getXlw();
604 TGraphErrors *Gdt = new TGraphErrors(nd,xx,yy,xer,yer);
605
606 myth=new TH1F("myth","Exp.data",600,xlow,xhigh);
607 Xsum=new TH1F("sumxs","sum of exclusive xs",600,xlow,xhigh);
608 double mstp=(xhigh-xlow)/600;
609 for(int i=0;i<600;i++){
610 double mx=i*mstp+xlow;
611 double xsi = myxsection->getXsection(mx);
612 if(xsi==0) continue;
613 myth->Fill(mx,xsi);
614 //std::cout<<mx<<" "<<xsi<<std::endl;
615 }
616
617 for(int i=0;i<600;i++){
618 double mx=i*mstp+xlow;
619 ysum[i]=sumExc(mx);
620 if(ysum[i]==0) continue;
621 Xsum->Fill(mx,ysum[i]);
622 //std::cout<<mx<<" "<<ysum[i]<<std::endl;
623 }
624
625 for(int i=0;i<nd;i++){
626 yysum[i]=sumExc(xx[i]);
627 }
628
629 myth->SetError(yer);
630 myth->Write();
631 Xsum->Write();
632
633 TGraphErrors *Gsum = new TGraphErrors(nd,xx,yysum,xer,yer);
634 //draw graph
635 double ex[600]={600*0};
636 double ey[600],ta[600];
637 double exz[600]={600*0};
638 double eyz[600]={600*0};
639 for(int i=0;i<600;i++){
640 ey[i]=AF[i]*0.001;
641 }
642 TGraphErrors *gr0 = new TGraphErrors(600,MH,AF,ex,ey);
643 TGraphErrors *gr1 = new TGraphErrors(600,MH,RadXS,exz,eyz);
644 TPostScript *ps = new TPostScript("xsobs.ps",111);
645 TCanvas *c1 = new TCanvas("c1","TGraphs for data",200,10,600,400);
646 gPad-> SetLogy(1);
647 // c1->SetLogy(1);
648 gStyle->SetCanvasColor(0);
649 gStyle->SetStatBorderSize(1);
650 c1->Divide(1,1);
651
652 c1->Update();
653 ps->NewPage();
654 c1->Draw();
655 c1->cd(1);
656 c1->SetLogy(1);
657 gr0->SetMarkerStyle(10);
658 gr0->Draw("AP");
659 gr0->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
660 gr0->GetYaxis()->SetTitle("observed cross section (nb)");
661 gr0->Write();
662
663 c1->Update();
664 ps->NewPage();
665 c1->Draw();
666 c1->cd(1);
667 c1->SetLogy(1);
668 gr1->SetMarkerStyle(10);
669 gr1->Draw("AP");
670 gr1->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
671 gr1->GetYaxis()->SetTitle("Rad*BornXS");
672 gr1->Write();
673
674 //--------- imposing graphs to a pad
675 TMultiGraph *mg = new TMultiGraph();
676 mg->SetTitle("Exclusion graphs");
677 Gdt->SetMarkerStyle(2);
678 Gdt->SetMarkerSize(1);
679 Gsum->SetLineColor(2);
680 Gsum->SetLineWidth(1);
681 mg->Add(Gdt);
682 mg->Add(Gsum);
683
684 c1->Update();
685 ps->NewPage();
686 c1->Draw();
687 c1->cd(1);
688 gPad-> SetLogy(1);
689 mg->Draw("APL");
690 mg->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
691 mg->GetYaxis()->SetTitle("observed cross section (nb)");
692 mg->Write();
693 //-------
694
695 c1->Update();
696 ps->NewPage();
697 c1->Draw();
698 c1->cd(1);
699 Gdt->SetMarkerStyle(2);
700 Gdt->Draw("AP");
701 Gdt->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
702 Gdt->GetYaxis()->SetTitle("observed cross section (nb)");
703 Gdt->Write();
704
705 c1->Update();
706 ps->NewPage();
707 c1->Draw();
708 c1->cd(1);
709 Gsum->SetMarkerStyle(2);
710 Gsum->SetMarkerSize(1);
711 Gsum->Draw("AP");
712 Gsum->SetLineWidth(0);
713 Gsum->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
714 Gsum->GetYaxis()->SetTitle("observed cross section (nb)");
715 Gsum->Write();
716
717 c1->Update();
718 ps->NewPage();
719 c1->Draw();
720 c1->cd(1);
721 // gPad-> SetLogx(1);
722 Gdt->SetMarkerStyle(2);
723 Gdt->SetMarkerSize(1);
724 Gdt->SetMaximum(8000.0);
725 Gsum->SetLineColor(2);
726 Gsum->SetLineWidth(1.5);
727 Gsum->Draw("AL"); //A: draw axis
728 Gdt->Draw("P"); // superpose to the Gsum, without A-option
729 Gsum->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
730 Gsum->GetYaxis()->SetTitle("cross section (nb)");
731 Gsum->Write();
732
733 ps->Close();
734 } //if(mydbg)_block
735 //-------------------------
736
737}
738
739
740//--
741void EvtConExc::init_mode(int mode){
742 if(mode==0){
743 _ndaugs =2;
744 daugs[0]=EvtPDL::getId(std::string("p+"));
745 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
746 }else if(mode==1){
747 _ndaugs =2;
748 daugs[0]=EvtPDL::getId(std::string("n0"));
749 daugs[1]=EvtPDL::getId(std::string("anti-n0"));
750 }else if(mode==2){
751 _ndaugs =2;
752 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
753 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
754 }else if(mode==3){
755 _ndaugs =2;
756 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
757 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
758 }else if(mode==4){
759 _ndaugs =2;
760 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
761 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
762 }else if(mode==5){
763 _ndaugs =2;
764 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
765 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
766 }else if(mode==6){
767 _ndaugs =2;
768 daugs[0]=EvtPDL::getId(std::string("pi+"));
769 daugs[1]=EvtPDL::getId(std::string("pi-"));
770 }else if(mode==7){
771 _ndaugs =3;
772 daugs[0]=EvtPDL::getId(std::string("pi+"));
773 daugs[1]=EvtPDL::getId(std::string("pi-"));
774 daugs[2]=EvtPDL::getId(std::string("pi0"));
775 }else if(mode==8){
776 _ndaugs =3;
777 daugs[0]=EvtPDL::getId(std::string("K+"));
778 daugs[1]=EvtPDL::getId(std::string("K-"));
779 daugs[2]=EvtPDL::getId(std::string("pi0"));
780 }else if(mode==9){
781 _ndaugs =3;
782 daugs[0]=EvtPDL::getId(std::string("K_S0"));
783 daugs[1]=EvtPDL::getId(std::string("K+"));
784 daugs[2]=EvtPDL::getId(std::string("pi-"));
785 }else if(mode==10){
786 _ndaugs =3;
787 daugs[0]=EvtPDL::getId(std::string("K_S0"));
788 daugs[1]=EvtPDL::getId(std::string("K-"));
789 daugs[2]=EvtPDL::getId(std::string("pi+"));
790 }else if(mode==11){
791 _ndaugs =3;
792 daugs[0]=EvtPDL::getId(std::string("K+"));
793 daugs[1]=EvtPDL::getId(std::string("K+"));
794 daugs[2]=EvtPDL::getId(std::string("eta"));
795 }else if(mode==12){
796 _ndaugs =4;
797 daugs[0]=EvtPDL::getId(std::string("pi+"));
798 daugs[1]=EvtPDL::getId(std::string("pi-"));
799 daugs[2]=EvtPDL::getId(std::string("pi+"));
800 daugs[3]=EvtPDL::getId(std::string("pi-"));
801 }else if(mode==13){
802 _ndaugs =4;
803 daugs[0]=EvtPDL::getId(std::string("pi+"));
804 daugs[1]=EvtPDL::getId(std::string("pi-"));
805 daugs[2]=EvtPDL::getId(std::string("pi0"));
806 daugs[3]=EvtPDL::getId(std::string("pi0"));
807 }else if(mode==14){
808 _ndaugs =4;
809 daugs[0]=EvtPDL::getId(std::string("K+"));
810 daugs[1]=EvtPDL::getId(std::string("K-"));
811 daugs[2]=EvtPDL::getId(std::string("pi+"));
812 daugs[3]=EvtPDL::getId(std::string("pi-"));
813 }else if(mode==15){
814 _ndaugs =4;
815 daugs[0]=EvtPDL::getId(std::string("K+"));
816 daugs[1]=EvtPDL::getId(std::string("K-"));
817 daugs[2]=EvtPDL::getId(std::string("pi0"));
818 daugs[3]=EvtPDL::getId(std::string("pi0"));
819 }else if(mode==16){
820 _ndaugs =4;
821 daugs[0]=EvtPDL::getId(std::string("K+"));
822 daugs[1]=EvtPDL::getId(std::string("K-"));
823 daugs[2]=EvtPDL::getId(std::string("K+"));
824 daugs[3]=EvtPDL::getId(std::string("K-"));
825 }else if(mode==17){
826 _ndaugs =5;
827 daugs[0]=EvtPDL::getId(std::string("pi+"));
828 daugs[1]=EvtPDL::getId(std::string("pi-"));
829 daugs[2]=EvtPDL::getId(std::string("pi+"));
830 daugs[3]=EvtPDL::getId(std::string("pi-"));
831 daugs[4]=EvtPDL::getId(std::string("pi0"));
832 }else if(mode==18){
833 _ndaugs =5;
834 daugs[0]=EvtPDL::getId(std::string("pi+"));
835 daugs[1]=EvtPDL::getId(std::string("pi-"));
836 daugs[2]=EvtPDL::getId(std::string("pi+"));
837 daugs[3]=EvtPDL::getId(std::string("pi-"));
838 daugs[4]=EvtPDL::getId(std::string("eta"));
839 }else if(mode==19){
840 _ndaugs =5;
841 daugs[0]=EvtPDL::getId(std::string("K+"));
842 daugs[1]=EvtPDL::getId(std::string("K-"));
843 daugs[2]=EvtPDL::getId(std::string("pi+"));
844 daugs[3]=EvtPDL::getId(std::string("pi-"));
845 daugs[4]=EvtPDL::getId(std::string("pi0"));
846 }else if(mode==20){
847 _ndaugs =5;
848 daugs[0]=EvtPDL::getId(std::string("K+"));
849 daugs[1]=EvtPDL::getId(std::string("K-"));
850 daugs[2]=EvtPDL::getId(std::string("pi+"));
851 daugs[3]=EvtPDL::getId(std::string("pi-"));
852 daugs[4]=EvtPDL::getId(std::string("eta"));
853 }else if(mode==21){
854 _ndaugs =6;
855 daugs[0]=EvtPDL::getId(std::string("pi+"));
856 daugs[1]=EvtPDL::getId(std::string("pi-"));
857 daugs[2]=EvtPDL::getId(std::string("pi+"));
858 daugs[3]=EvtPDL::getId(std::string("pi-"));
859 daugs[4]=EvtPDL::getId(std::string("pi+"));
860 daugs[5]=EvtPDL::getId(std::string("pi-"));
861 }else if(mode==22){
862 _ndaugs =6;
863 daugs[0]=EvtPDL::getId(std::string("pi+"));
864 daugs[1]=EvtPDL::getId(std::string("pi-"));
865 daugs[2]=EvtPDL::getId(std::string("pi+"));
866 daugs[3]=EvtPDL::getId(std::string("pi-"));
867 daugs[4]=EvtPDL::getId(std::string("pi0"));
868 daugs[5]=EvtPDL::getId(std::string("pi0"));
869 }else if(mode==23){
870 _ndaugs =2;
871 daugs[0]=EvtPDL::getId(std::string("eta"));
872 daugs[1]=EvtPDL::getId(std::string("phi"));
873 }else if(mode==24){
874 _ndaugs =2;
875 daugs[0]=EvtPDL::getId(std::string("phi"));
876 daugs[1]=EvtPDL::getId(std::string("pi0"));
877 }else if(mode==25){
878 _ndaugs =2;
879 daugs[0]=EvtPDL::getId(std::string("K+"));
880 daugs[1]=EvtPDL::getId(std::string("K*-"));
881 }else if(mode==26){
882 _ndaugs =2;
883 daugs[0]=EvtPDL::getId(std::string("K-"));
884 daugs[1]=EvtPDL::getId(std::string("K*+"));
885 }else if(mode==27){
886 _ndaugs =2;
887 daugs[0]=EvtPDL::getId(std::string("K_S0"));
888 daugs[1]=EvtPDL::getId(std::string("anti-K*0"));
889 }else if(mode==28){
890 _ndaugs =3;
891 daugs[0]=EvtPDL::getId(std::string("anti-K*0"));
892 daugs[1]=EvtPDL::getId(std::string("K+"));
893 daugs[2]=EvtPDL::getId(std::string("pi-"));
894 }else if(mode==29){
895 _ndaugs =3;
896 daugs[0]=EvtPDL::getId(std::string("K*0"));
897 daugs[1]=EvtPDL::getId(std::string("K-"));
898 daugs[2]=EvtPDL::getId(std::string("pi+"));
899 }else if(mode==30){
900 _ndaugs =3;
901 daugs[0]=EvtPDL::getId(std::string("K*+"));
902 daugs[1]=EvtPDL::getId(std::string("K-"));
903 daugs[2]=EvtPDL::getId(std::string("pi0"));
904 }else if(mode==31){
905 _ndaugs =3;
906 daugs[0]=EvtPDL::getId(std::string("K*-"));
907 daugs[1]=EvtPDL::getId(std::string("K+"));
908 daugs[2]=EvtPDL::getId(std::string("pi0"));
909 }else if(mode==32){
910 _ndaugs =3;
911 daugs[0]=EvtPDL::getId(std::string("anti-K_2*0"));
912 daugs[1]=EvtPDL::getId(std::string("K+"));
913 daugs[2]=EvtPDL::getId(std::string("pi-"));
914 }else if(mode==33){
915 _ndaugs =3;
916 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
917 daugs[1]=EvtPDL::getId(std::string("K-"));
918 daugs[2]=EvtPDL::getId(std::string("pi+"));
919 }else if(mode==34){
920 _ndaugs =3;
921 daugs[0]=EvtPDL::getId(std::string("K+"));
922 daugs[1]=EvtPDL::getId(std::string("K-"));
923 daugs[2]=EvtPDL::getId(std::string("rho0"));
924 }else if(mode==35){
925 _ndaugs =3;
926 daugs[0]=EvtPDL::getId(std::string("phi"));
927 daugs[1]=EvtPDL::getId(std::string("pi-"));
928 daugs[2]=EvtPDL::getId(std::string("pi+"));
929 }else if(mode==36){
930 _ndaugs =2;
931 daugs[0]=EvtPDL::getId(std::string("phi"));
932 daugs[1]=EvtPDL::getId(std::string("f_0"));
933 }else if(mode==37){
934 _ndaugs =3;
935 daugs[0]=EvtPDL::getId(std::string("eta"));
936 daugs[1]=EvtPDL::getId(std::string("pi+"));
937 daugs[2]=EvtPDL::getId(std::string("pi-"));
938 }else if(mode==38){
939 _ndaugs =3;
940 daugs[0]=EvtPDL::getId(std::string("omega"));
941 daugs[1]=EvtPDL::getId(std::string("pi+"));
942 daugs[2]=EvtPDL::getId(std::string("pi-"));
943 }else if(mode==39){
944 _ndaugs =2;
945 daugs[0]=EvtPDL::getId(std::string("omega"));
946 daugs[1]=EvtPDL::getId(std::string("f_0"));
947 }else if(mode==40){
948 _ndaugs =3;
949 daugs[0]=EvtPDL::getId(std::string("eta'"));
950 daugs[1]=EvtPDL::getId(std::string("pi+"));
951 daugs[2]=EvtPDL::getId(std::string("pi-"));
952 }else if(mode==41){
953 _ndaugs =3;
954 daugs[0]=EvtPDL::getId(std::string("f_1"));
955 daugs[1]=EvtPDL::getId(std::string("pi+"));
956 daugs[2]=EvtPDL::getId(std::string("pi-"));
957 }else if(mode==42){
958 _ndaugs =3;
959 daugs[0]=EvtPDL::getId(std::string("omega"));
960 daugs[1]=EvtPDL::getId(std::string("K+"));
961 daugs[2]=EvtPDL::getId(std::string("K-"));
962 }else if(mode==43){
963 _ndaugs =4;
964 daugs[0]=EvtPDL::getId(std::string("omega"));
965 daugs[1]=EvtPDL::getId(std::string("pi+"));
966 daugs[2]=EvtPDL::getId(std::string("pi-"));
967 daugs[3]=EvtPDL::getId(std::string("pi0"));
968 }else if(mode==44){
969 _ndaugs =2;
970 daugs[0]=EvtPDL::getId(std::string("Sigma-"));
971 daugs[1]=EvtPDL::getId(std::string("anti-Sigma+"));
972 }else if(mode==45){
973 _ndaugs =2;
974 daugs[0]=EvtPDL::getId(std::string("K+"));
975 daugs[1]=EvtPDL::getId(std::string("K-"));
976 }else if(mode==46){
977 _ndaugs =2;
978 daugs[0]=EvtPDL::getId(std::string("K_S0"));
979 daugs[1]=EvtPDL::getId(std::string("K_L0"));
980 }else if(mode==47){
981 _ndaugs =2;
982 daugs[0]=EvtPDL::getId(std::string("omega"));
983 daugs[1]=EvtPDL::getId(std::string("eta"));
984 }else if(mode==48){
985 _ndaugs =2;
986 daugs[0]=EvtPDL::getId(std::string("phi"));
987 daugs[1]=EvtPDL::getId(std::string("eta'"));
988 }else if(mode==49){
989 _ndaugs =3;
990 daugs[0]=EvtPDL::getId(std::string("phi"));
991 daugs[1]=EvtPDL::getId(std::string("K+"));
992 daugs[2]=EvtPDL::getId(std::string("K-"));
993 }else if(mode==50){
994 _ndaugs =3;
995 daugs[0]=EvtPDL::getId(std::string("p+"));
996 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
997 daugs[2]=EvtPDL::getId(std::string("pi0"));
998 }else if(mode==51){
999 _ndaugs =3;
1000 daugs[0]=EvtPDL::getId(std::string("p+"));
1001 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
1002 daugs[2]=EvtPDL::getId(std::string("eta"));
1003 }else if(mode==52){
1004 _ndaugs =2;
1005 daugs[0]=EvtPDL::getId(std::string("omega"));
1006 daugs[1]=EvtPDL::getId(std::string("pi0"));
1007 }else if(mode==59){
1008 _ndaugs =3;
1009 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1010 daugs[1]=EvtPDL::getId(std::string("D0"));
1011 daugs[2]=EvtPDL::getId(std::string("pi0"));
1012 }else if(mode==60){
1013 _ndaugs =3;
1014 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1015 daugs[1]=EvtPDL::getId(std::string("D*0"));
1016 daugs[2]=EvtPDL::getId(std::string("pi0"));
1017 }else if(mode==61){
1018 _ndaugs =3;
1019 daugs[0]=EvtPDL::getId(std::string("D0"));
1020 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1021 daugs[2]=EvtPDL::getId(std::string("pi0"));
1022 }else if(mode==62){
1023 _ndaugs =3;
1024 daugs[0]=EvtPDL::getId(std::string("D+"));
1025 daugs[1]=EvtPDL::getId(std::string("D*-"));
1026 daugs[2]=EvtPDL::getId(std::string("pi0"));
1027 }else if(mode==63){
1028 _ndaugs =3;
1029 daugs[0]=EvtPDL::getId(std::string("D-"));
1030 daugs[1]=EvtPDL::getId(std::string("D+"));
1031 daugs[2]=EvtPDL::getId(std::string("pi0"));
1032 }else if(mode==64){
1033 _ndaugs =3;
1034 daugs[0]=EvtPDL::getId(std::string("D-"));
1035 daugs[1]=EvtPDL::getId(std::string("D*+"));
1036 daugs[2]=EvtPDL::getId(std::string("pi0"));
1037 }else if(mode==65){
1038 _ndaugs =3;
1039 daugs[0]=EvtPDL::getId(std::string("D-"));
1040 daugs[1]=EvtPDL::getId(std::string("D*0"));
1041 daugs[2]=EvtPDL::getId(std::string("pi+"));
1042 }else if(mode==66){
1043 _ndaugs =3;
1044 daugs[0]=EvtPDL::getId(std::string("D+"));
1045 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1046 daugs[2]=EvtPDL::getId(std::string("pi-"));
1047 }else if(mode==67){
1048 _ndaugs =2;
1049 daugs[0]=EvtPDL::getId(std::string("D*0"));
1050 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1051 }else if(mode==68){
1052 _ndaugs =2;
1053 daugs[0]=EvtPDL::getId(std::string("D0"));
1054 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1055 }else if(mode==69){
1056 _ndaugs =2;
1057 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1058 daugs[1]=EvtPDL::getId(std::string("D*0"));
1059 }else if(mode==70){
1060 _ndaugs =2;
1061 daugs[0]=EvtPDL::getId(std::string("D0"));
1062 daugs[1]=EvtPDL::getId(std::string("anti-D0"));
1063 }else if(mode==71){
1064 _ndaugs =2;
1065 daugs[0]=EvtPDL::getId(std::string("D+"));
1066 daugs[1]=EvtPDL::getId(std::string("D-"));
1067 }else if(mode==72){
1068 _ndaugs =2;
1069 daugs[0]=EvtPDL::getId(std::string("D+"));
1070 daugs[1]=EvtPDL::getId(std::string("D*-"));
1071 }else if(mode==73){
1072 _ndaugs =2;
1073 daugs[0]=EvtPDL::getId(std::string("D-"));
1074 daugs[1]=EvtPDL::getId(std::string("D*+"));
1075 }else if(mode==74){
1076 _ndaugs =2;
1077 daugs[0]=EvtPDL::getId(std::string("D*+"));
1078 daugs[1]=EvtPDL::getId(std::string("D*-"));
1079 }else if(mode==75){
1080 _ndaugs =3;
1081 daugs[0]=EvtPDL::getId(std::string("D0"));
1082 daugs[1]=EvtPDL::getId(std::string("D-"));
1083 daugs[2]=EvtPDL::getId(std::string("pi+"));
1084 }else if(mode==76){
1085 _ndaugs =3;
1086 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1087 daugs[1]=EvtPDL::getId(std::string("D+"));
1088 daugs[2]=EvtPDL::getId(std::string("pi-"));
1089 }else if(mode==77){
1090 _ndaugs =3;
1091 daugs[0]=EvtPDL::getId(std::string("D0"));
1092 daugs[1]=EvtPDL::getId(std::string("D*-"));
1093 daugs[2]=EvtPDL::getId(std::string("pi+"));
1094 }else if(mode==78){
1095 _ndaugs =3;
1096 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1097 daugs[1]=EvtPDL::getId(std::string("D*+"));
1098 daugs[2]=EvtPDL::getId(std::string("pi-"));
1099 }else if(mode==79){
1100 _ndaugs =3;
1101 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1102 daugs[1]=EvtPDL::getId(std::string("pi0"));
1103 daugs[2]=EvtPDL::getId(std::string("pi0"));
1104 }else if(mode==80){
1105 _ndaugs =2;
1106 daugs[0]=EvtPDL::getId(std::string("eta"));
1107 daugs[1]=EvtPDL::getId(std::string("J/psi"));
1108 }else if(mode==81){
1109 _ndaugs =3;
1110 daugs[0]=EvtPDL::getId(std::string("h_c"));
1111 daugs[1]=EvtPDL::getId(std::string("pi+"));
1112 daugs[2]=EvtPDL::getId(std::string("pi-"));
1113 }else if(mode==82){
1114 _ndaugs =3;
1115 daugs[0]=EvtPDL::getId(std::string("h_c"));
1116 daugs[1]=EvtPDL::getId(std::string("pi0"));
1117 daugs[2]=EvtPDL::getId(std::string("pi0"));
1118 }else if(mode==83){
1119 _ndaugs =3;
1120 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1121 daugs[1]=EvtPDL::getId(std::string("K+"));
1122 daugs[2]=EvtPDL::getId(std::string("K-"));
1123 }else if(mode==84){
1124 _ndaugs =3;
1125 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1126 daugs[1]=EvtPDL::getId(std::string("K_S0"));
1127 daugs[2]=EvtPDL::getId(std::string("K_S0"));
1128 }else if(mode==85){
1129 _ndaugs =2;
1130 daugs[0]=EvtPDL::getId(std::string("D_s*+"));
1131 daugs[1]=EvtPDL::getId(std::string("D_s*-"));
1132 }else if(mode==90){
1133 _ndaugs =3;
1134 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1135 daugs[1]=EvtPDL::getId(std::string("pi+"));
1136 daugs[2]=EvtPDL::getId(std::string("pi-"));
1137 }else if(mode==91){
1138 _ndaugs =3;
1139 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1140 daugs[1]=EvtPDL::getId(std::string("pi+"));
1141 daugs[2]=EvtPDL::getId(std::string("pi-"));
1142 }else if(mode==92){
1143 _ndaugs =3;
1144 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1145 daugs[1]=EvtPDL::getId(std::string("pi0"));
1146 daugs[2]=EvtPDL::getId(std::string("pi0"));
1147 }else if(mode==93){
1148 _ndaugs =2;
1149 daugs[0]=EvtPDL::getId(std::string("D_s+"));
1150 daugs[1]=EvtPDL::getId(std::string("D_s-"));
1151 }else if(mode==94){
1152 _ndaugs =2;
1153 daugs[0]=EvtPDL::getId(std::string("D_s*+"));
1154 daugs[1]=EvtPDL::getId(std::string("D_s-"));
1155 }else if(mode==95){
1156 _ndaugs =2;
1157 daugs[0]=EvtPDL::getId(std::string("D_s*-"));
1158 daugs[1]=EvtPDL::getId(std::string("D_s+"));
1159 }else if(mode==96){
1160 _ndaugs =2;
1161 daugs[0]=EvtPDL::getId(std::string("Lambda_c+"));
1162 daugs[1]=EvtPDL::getId(std::string("anti-Lambda_c-"));
1163 }else if(mode==111){
1164 _ndaugs =2;
1165 daugs[0]=EvtPDL::getId(std::string("Zc(4020)+"));
1166 daugs[1]=EvtPDL::getId(std::string("pi-"));
1167 }else if(mode==112){
1168 _ndaugs =2;
1169 daugs[0]=EvtPDL::getId(std::string("Zc(4020)-"));
1170 daugs[1]=EvtPDL::getId(std::string("pi+"));
1171 }else if(mode==74100){
1172 _ndaugs =1;
1173 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1174 }else if(mode==74101){
1175 _ndaugs =1;
1176 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1177 }else if(mode==74102){
1178 _ndaugs =1;
1179 daugs[0]=EvtPDL::getId(std::string("psi(3770)"));
1180 }else if(mode==74103){
1181 _ndaugs =1;
1182 daugs[0]=EvtPDL::getId(std::string("phi"));
1183 }else if(mode==74104){
1184 _ndaugs =1;
1185 daugs[0]=EvtPDL::getId(std::string("omega"));
1186 }else if(mode==74105){
1187 _ndaugs =1;
1188 daugs[0]=EvtPDL::getId(std::string("rho0"));
1189 }else if(mode==74106){
1190 _ndaugs =1;
1191 daugs[0]=EvtPDL::getId(std::string("rho(3S)0"));
1192 }else if(mode==74107){
1193 _ndaugs =1;
1194 daugs[0]=EvtPDL::getId(std::string("omega(2S)"));
1195 }else if(mode==74110){
1196 _ndaugs =1;
1197 daugs[0]=EvtPDL::getId(std::string("gamma*"));
1198 EvtId myvpho = EvtPDL::getId(std::string("gamma*"));
1199 EvtPDL::reSetMass(myvpho,mhds*0.9); //for calculating maxXS, it will be resent when this mode is selected
1200 _modeFlag.clear();
1201 _modeFlag.push_back(74110); //R-value generator tag
1202 _modeFlag.push_back(0); //to contain the mode selected
1203 }else if(mode==-1){
1204 _ndaugs = nson;
1205 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1206 std::cout<<"The paraent particle: "<<EvtPDL::name(pid)<<std::endl;
1207 }else if(mode==-2){
1208 _ndaugs = nson;
1209 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1210 }else if(mode==-100){
1211 _ndaugs =1;
1212 daugs[0]=EvtPDL::getId(std::string("gamma*"));
1213 _modeFlag.clear();
1214 _modeFlag.push_back(-100); //R-value generator tag
1215 _modeFlag.push_back(0); //to contain the mode selected
1216 return;
1217 }else if(mode==10000){//use for check ISR
1218 _ndaugs =2;
1219 daugs[0]=EvtPDL::getId(std::string("pi+"));
1220 daugs[1]=EvtPDL::getId(std::string("pi-"));
1221 }else {
1222 std::cout<<"Bad mode_index number " <<mode<<", please refer to the manual."<<std::endl;
1223 ::abort();
1224 }
1225 // std::cout<<_ndaugs<<", "<<daugs[0]<<daugs[1]<<std::endl;
1226
1227 if(VISR){
1228 _ndaugs=2;
1229 daugs[0]=getDaugs()[0];//EvtPDL::getId(std::string("gamma"));
1230 daugs[1]=getDaugs()[1];//EvtPDL::getId(std::string("gamma*"));
1231 }
1232
1233 double fmth=0;
1234 for(int i=0;i<_ndaugs;i++){
1235 double xmi=EvtPDL::getMinMass(daugs[i]);
1236 fmth +=xmi;
1237 }
1238 myxsection = new EvtXsection (mode);
1239 Mthr=myxsection->getXlw();
1240 if(!(_mode==74110 || _mode==-100)){
1241 if(Mthr<fmth) {std::cout<<"Low end of cross section: ("<<Mthr<<" < (mass of final state)"<<fmth<<") GeV."<< std::endl; abort();}
1242 }
1243}
1244
1245//--
1246std::vector<EvtId> EvtConExc::get_mode(int mode){
1247 int _ndaugs;
1248 EvtId daugs[100];
1249 if(mode==0){
1250 _ndaugs =2;
1251 daugs[0]=EvtPDL::getId(std::string("p+"));
1252 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
1253 }else if(mode==1){
1254 _ndaugs =2;
1255 daugs[0]=EvtPDL::getId(std::string("n0"));
1256 daugs[1]=EvtPDL::getId(std::string("anti-n0"));
1257 }else if(mode==2){
1258 _ndaugs =2;
1259 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
1260 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
1261 }else if(mode==3){
1262 _ndaugs =2;
1263 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
1264 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
1265 }else if(mode==4){
1266 _ndaugs =2;
1267 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
1268 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
1269 }else if(mode==5){
1270 _ndaugs =2;
1271 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
1272 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
1273 }else if(mode ==6){
1274 _ndaugs =2;
1275 daugs[0]=EvtPDL::getId(std::string("pi+"));
1276 daugs[1]=EvtPDL::getId(std::string("pi-"));
1277 }else if(mode==7){
1278 _ndaugs =3;
1279 daugs[0]=EvtPDL::getId(std::string("pi+"));
1280 daugs[1]=EvtPDL::getId(std::string("pi-"));
1281 daugs[2]=EvtPDL::getId(std::string("pi0"));
1282 }else if(mode==8){
1283 _ndaugs =3;
1284 daugs[0]=EvtPDL::getId(std::string("K+"));
1285 daugs[1]=EvtPDL::getId(std::string("K-"));
1286 daugs[2]=EvtPDL::getId(std::string("pi0"));
1287 }else if(mode==9){
1288 _ndaugs =3;
1289 daugs[0]=EvtPDL::getId(std::string("K_S0"));
1290 daugs[1]=EvtPDL::getId(std::string("K+"));
1291 daugs[2]=EvtPDL::getId(std::string("pi-"));
1292 }else if(mode==10){
1293 _ndaugs =3;
1294 daugs[0]=EvtPDL::getId(std::string("K_S0"));
1295 daugs[1]=EvtPDL::getId(std::string("K-"));
1296 daugs[2]=EvtPDL::getId(std::string("pi+"));
1297 }else if(mode==11){
1298 _ndaugs =3;
1299 daugs[0]=EvtPDL::getId(std::string("K+"));
1300 daugs[1]=EvtPDL::getId(std::string("K-"));
1301 daugs[2]=EvtPDL::getId(std::string("eta"));
1302 }else if(mode==12){
1303 _ndaugs =4;
1304 daugs[0]=EvtPDL::getId(std::string("pi+"));
1305 daugs[1]=EvtPDL::getId(std::string("pi-"));
1306 daugs[2]=EvtPDL::getId(std::string("pi+"));
1307 daugs[3]=EvtPDL::getId(std::string("pi-"));
1308 }else if(mode==13){
1309 _ndaugs =4;
1310 daugs[0]=EvtPDL::getId(std::string("pi+"));
1311 daugs[1]=EvtPDL::getId(std::string("pi-"));
1312 daugs[2]=EvtPDL::getId(std::string("pi0"));
1313 daugs[3]=EvtPDL::getId(std::string("pi0"));
1314 }else if(mode==14){
1315 _ndaugs =4;
1316 daugs[0]=EvtPDL::getId(std::string("K+"));
1317 daugs[1]=EvtPDL::getId(std::string("K-"));
1318 daugs[2]=EvtPDL::getId(std::string("pi+"));
1319 daugs[3]=EvtPDL::getId(std::string("pi-"));
1320 }else if(mode==15){
1321 _ndaugs =4;
1322 daugs[0]=EvtPDL::getId(std::string("K+"));
1323 daugs[1]=EvtPDL::getId(std::string("K-"));
1324 daugs[2]=EvtPDL::getId(std::string("pi0"));
1325 daugs[3]=EvtPDL::getId(std::string("pi0"));
1326 }else if(mode==16){
1327 _ndaugs =4;
1328 daugs[0]=EvtPDL::getId(std::string("K+"));
1329 daugs[1]=EvtPDL::getId(std::string("K-"));
1330 daugs[2]=EvtPDL::getId(std::string("K+"));
1331 daugs[3]=EvtPDL::getId(std::string("K-"));
1332 }else if(mode==17){
1333 _ndaugs =5;
1334 daugs[0]=EvtPDL::getId(std::string("pi+"));
1335 daugs[1]=EvtPDL::getId(std::string("pi-"));
1336 daugs[2]=EvtPDL::getId(std::string("pi+"));
1337 daugs[3]=EvtPDL::getId(std::string("pi-"));
1338 daugs[4]=EvtPDL::getId(std::string("pi0"));
1339 }else if(mode==18){
1340 _ndaugs =5;
1341 daugs[0]=EvtPDL::getId(std::string("pi+"));
1342 daugs[1]=EvtPDL::getId(std::string("pi-"));
1343 daugs[2]=EvtPDL::getId(std::string("pi+"));
1344 daugs[3]=EvtPDL::getId(std::string("pi-"));
1345 daugs[4]=EvtPDL::getId(std::string("eta"));
1346 }else if(mode==19){
1347 _ndaugs =5;
1348 daugs[0]=EvtPDL::getId(std::string("K+"));
1349 daugs[1]=EvtPDL::getId(std::string("K-"));
1350 daugs[2]=EvtPDL::getId(std::string("pi+"));
1351 daugs[3]=EvtPDL::getId(std::string("pi-"));
1352 daugs[4]=EvtPDL::getId(std::string("pi0"));
1353 }else if(mode==20){
1354 _ndaugs =5;
1355 daugs[0]=EvtPDL::getId(std::string("K+"));
1356 daugs[1]=EvtPDL::getId(std::string("K-"));
1357 daugs[2]=EvtPDL::getId(std::string("pi+"));
1358 daugs[3]=EvtPDL::getId(std::string("pi-"));
1359 daugs[4]=EvtPDL::getId(std::string("eta"));
1360 }else if(mode==21){
1361 _ndaugs =6;
1362 daugs[0]=EvtPDL::getId(std::string("pi+"));
1363 daugs[1]=EvtPDL::getId(std::string("pi-"));
1364 daugs[2]=EvtPDL::getId(std::string("pi+"));
1365 daugs[3]=EvtPDL::getId(std::string("pi-"));
1366 daugs[4]=EvtPDL::getId(std::string("pi+"));
1367 daugs[5]=EvtPDL::getId(std::string("pi-"));
1368 }else if(mode==22){
1369 _ndaugs =6;
1370 daugs[0]=EvtPDL::getId(std::string("pi+"));
1371 daugs[1]=EvtPDL::getId(std::string("pi-"));
1372 daugs[2]=EvtPDL::getId(std::string("pi+"));
1373 daugs[3]=EvtPDL::getId(std::string("pi-"));
1374 daugs[4]=EvtPDL::getId(std::string("pi0"));
1375 daugs[5]=EvtPDL::getId(std::string("pi0"));
1376 }else if(mode==23){
1377 _ndaugs =2;
1378 daugs[0]=EvtPDL::getId(std::string("eta"));
1379 daugs[1]=EvtPDL::getId(std::string("phi"));
1380 }else if(mode==24){
1381 _ndaugs =2;
1382 daugs[0]=EvtPDL::getId(std::string("phi"));
1383 daugs[1]=EvtPDL::getId(std::string("pi0"));
1384 }else if(mode==25){
1385 _ndaugs =2;
1386 daugs[0]=EvtPDL::getId(std::string("K+"));
1387 daugs[1]=EvtPDL::getId(std::string("K*-"));
1388 }else if(mode==26){
1389 _ndaugs =2;
1390 daugs[0]=EvtPDL::getId(std::string("K-"));
1391 daugs[1]=EvtPDL::getId(std::string("K*+"));
1392 }else if(mode==27){
1393 _ndaugs =2;
1394 daugs[0]=EvtPDL::getId(std::string("K_S0"));
1395 daugs[1]=EvtPDL::getId(std::string("anti-K*0"));
1396 }else if(mode==28){
1397 _ndaugs =3;
1398 daugs[0]=EvtPDL::getId(std::string("anti-K*0"));
1399 daugs[1]=EvtPDL::getId(std::string("K+"));
1400 daugs[2]=EvtPDL::getId(std::string("pi-"));
1401 }else if(mode==29){
1402 _ndaugs =3;
1403 daugs[0]=EvtPDL::getId(std::string("K*0"));
1404 daugs[1]=EvtPDL::getId(std::string("K-"));
1405 daugs[2]=EvtPDL::getId(std::string("pi+"));
1406 }else if(mode==30){
1407 _ndaugs =3;
1408 daugs[0]=EvtPDL::getId(std::string("K*+"));
1409 daugs[1]=EvtPDL::getId(std::string("K-"));
1410 daugs[2]=EvtPDL::getId(std::string("pi0"));
1411 }else if(mode==31){
1412 _ndaugs =3;
1413 daugs[0]=EvtPDL::getId(std::string("K*-"));
1414 daugs[1]=EvtPDL::getId(std::string("K+"));
1415 daugs[2]=EvtPDL::getId(std::string("pi0"));
1416 }else if(mode==32){
1417 _ndaugs =3;
1418 daugs[0]=EvtPDL::getId(std::string("anti-K_2*0"));
1419 daugs[1]=EvtPDL::getId(std::string("K+"));
1420 daugs[2]=EvtPDL::getId(std::string("pi-"));
1421 }else if(mode==33){
1422 _ndaugs =3;
1423 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
1424 daugs[1]=EvtPDL::getId(std::string("K-"));
1425 daugs[2]=EvtPDL::getId(std::string("pi+"));
1426 }else if(mode==34){
1427 _ndaugs =3;
1428 daugs[0]=EvtPDL::getId(std::string("K+"));
1429 daugs[1]=EvtPDL::getId(std::string("K-"));
1430 daugs[2]=EvtPDL::getId(std::string("rho0"));
1431 }else if(mode==35){
1432 _ndaugs =3;
1433 daugs[0]=EvtPDL::getId(std::string("phi"));
1434 daugs[1]=EvtPDL::getId(std::string("pi-"));
1435 daugs[2]=EvtPDL::getId(std::string("pi+"));
1436 }else if(mode==36){
1437 _ndaugs =2;
1438 daugs[0]=EvtPDL::getId(std::string("phi"));
1439 daugs[1]=EvtPDL::getId(std::string("f_0"));
1440 }else if(mode==37){
1441 _ndaugs =3;
1442 daugs[0]=EvtPDL::getId(std::string("eta"));
1443 daugs[1]=EvtPDL::getId(std::string("pi+"));
1444 daugs[2]=EvtPDL::getId(std::string("pi-"));
1445 }else if(mode==38){
1446 _ndaugs =3;
1447 daugs[0]=EvtPDL::getId(std::string("omega"));
1448 daugs[1]=EvtPDL::getId(std::string("pi+"));
1449 daugs[2]=EvtPDL::getId(std::string("pi-"));
1450 }else if(mode==39){
1451 _ndaugs =2;
1452 daugs[0]=EvtPDL::getId(std::string("omega"));
1453 daugs[1]=EvtPDL::getId(std::string("f_0"));
1454 }else if(mode==40){
1455 _ndaugs =3;
1456 daugs[0]=EvtPDL::getId(std::string("eta'"));
1457 daugs[1]=EvtPDL::getId(std::string("pi+"));
1458 daugs[2]=EvtPDL::getId(std::string("pi-"));
1459 }else if(mode==41){
1460 _ndaugs =3;
1461 daugs[0]=EvtPDL::getId(std::string("f_1"));
1462 daugs[1]=EvtPDL::getId(std::string("pi+"));
1463 daugs[2]=EvtPDL::getId(std::string("pi-"));
1464 }else if(mode==42){
1465 _ndaugs =3;
1466 daugs[0]=EvtPDL::getId(std::string("omega"));
1467 daugs[1]=EvtPDL::getId(std::string("K+"));
1468 daugs[2]=EvtPDL::getId(std::string("K-"));
1469 }else if(mode==43){
1470 _ndaugs =4;
1471 daugs[0]=EvtPDL::getId(std::string("omega"));
1472 daugs[1]=EvtPDL::getId(std::string("pi+"));
1473 daugs[2]=EvtPDL::getId(std::string("pi-"));
1474 daugs[3]=EvtPDL::getId(std::string("pi0"));
1475 }else if(mode==44){
1476 _ndaugs =2;
1477 daugs[0]=EvtPDL::getId(std::string("Sigma-"));
1478 daugs[1]=EvtPDL::getId(std::string("anti-Sigma+"));
1479 }else if(mode==45){
1480 _ndaugs =2;
1481 daugs[0]=EvtPDL::getId(std::string("K+"));
1482 daugs[1]=EvtPDL::getId(std::string("K-"));
1483 }else if(mode==46){
1484 _ndaugs =2;
1485 daugs[0]=EvtPDL::getId(std::string("K_S0"));
1486 daugs[1]=EvtPDL::getId(std::string("K_L0"));
1487 }else if(mode==47){
1488 _ndaugs =2;
1489 daugs[0]=EvtPDL::getId(std::string("omega"));
1490 daugs[1]=EvtPDL::getId(std::string("eta"));
1491 }else if(mode==48){
1492 _ndaugs =2;
1493 daugs[0]=EvtPDL::getId(std::string("phi"));
1494 daugs[1]=EvtPDL::getId(std::string("eta'"));
1495 }else if(mode==49){
1496 _ndaugs =3;
1497 daugs[0]=EvtPDL::getId(std::string("phi"));
1498 daugs[1]=EvtPDL::getId(std::string("K+"));
1499 daugs[2]=EvtPDL::getId(std::string("K-"));
1500 }else if(mode==50){
1501 _ndaugs =3;
1502 daugs[0]=EvtPDL::getId(std::string("p+"));
1503 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
1504 daugs[2]=EvtPDL::getId(std::string("pi0"));
1505 }else if(mode==51){
1506 _ndaugs =3;
1507 daugs[0]=EvtPDL::getId(std::string("p+"));
1508 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
1509 daugs[2]=EvtPDL::getId(std::string("eta"));
1510 }else if(mode==52){
1511 _ndaugs =2;
1512 daugs[0]=EvtPDL::getId(std::string("omega"));
1513 daugs[1]=EvtPDL::getId(std::string("pi0"));
1514 }else if(mode==59){
1515 _ndaugs =3;
1516 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1517 daugs[1]=EvtPDL::getId(std::string("D0"));
1518 daugs[2]=EvtPDL::getId(std::string("pi0"));
1519 }else if(mode==60){
1520 _ndaugs =3;
1521 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1522 daugs[1]=EvtPDL::getId(std::string("D*0"));
1523 daugs[2]=EvtPDL::getId(std::string("pi0"));
1524 }else if(mode==61){
1525 _ndaugs =3;
1526 daugs[0]=EvtPDL::getId(std::string("D0"));
1527 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1528 daugs[2]=EvtPDL::getId(std::string("pi0"));
1529 }else if(mode==62){
1530 _ndaugs =3;
1531 daugs[0]=EvtPDL::getId(std::string("D+"));
1532 daugs[1]=EvtPDL::getId(std::string("D*-"));
1533 daugs[2]=EvtPDL::getId(std::string("pi0"));
1534 }else if(mode==63){
1535 _ndaugs =3;
1536 daugs[0]=EvtPDL::getId(std::string("D-"));
1537 daugs[1]=EvtPDL::getId(std::string("D+"));
1538 daugs[2]=EvtPDL::getId(std::string("pi0"));
1539 }else if(mode==64){
1540 _ndaugs =3;
1541 daugs[0]=EvtPDL::getId(std::string("D-"));
1542 daugs[1]=EvtPDL::getId(std::string("D*+"));
1543 daugs[2]=EvtPDL::getId(std::string("pi0"));
1544 }else if(mode==65){
1545 _ndaugs =3;
1546 daugs[0]=EvtPDL::getId(std::string("D-"));
1547 daugs[1]=EvtPDL::getId(std::string("D*0"));
1548 daugs[2]=EvtPDL::getId(std::string("pi+"));
1549 }else if(mode==66){
1550 _ndaugs =3;
1551 daugs[0]=EvtPDL::getId(std::string("D+"));
1552 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1553 daugs[2]=EvtPDL::getId(std::string("pi-"));
1554 }else if(mode==67){
1555 _ndaugs =2;
1556 daugs[0]=EvtPDL::getId(std::string("D*0"));
1557 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1558 }else if(mode==68){
1559 _ndaugs =2;
1560 daugs[0]=EvtPDL::getId(std::string("D0"));
1561 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
1562 }else if(mode==69){
1563 _ndaugs =2;
1564 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1565 daugs[1]=EvtPDL::getId(std::string("D*0"));
1566 }else if(mode==70){
1567 _ndaugs =2;
1568 daugs[0]=EvtPDL::getId(std::string("D0"));
1569 daugs[1]=EvtPDL::getId(std::string("anti-D0"));
1570 }else if(mode==71){
1571 _ndaugs =2;
1572 daugs[0]=EvtPDL::getId(std::string("D+"));
1573 daugs[1]=EvtPDL::getId(std::string("D-"));
1574 }else if(mode==72){
1575 _ndaugs =2;
1576 daugs[0]=EvtPDL::getId(std::string("D+"));
1577 daugs[1]=EvtPDL::getId(std::string("D*-"));
1578 }else if(mode==73){
1579 _ndaugs =2;
1580 daugs[0]=EvtPDL::getId(std::string("D-"));
1581 daugs[1]=EvtPDL::getId(std::string("D*+"));
1582 }else if(mode==74){
1583 _ndaugs =2;
1584 daugs[0]=EvtPDL::getId(std::string("D*+"));
1585 daugs[1]=EvtPDL::getId(std::string("D*-"));
1586 }else if(mode==75){
1587 _ndaugs =3;
1588 daugs[0]=EvtPDL::getId(std::string("D0"));
1589 daugs[1]=EvtPDL::getId(std::string("D-"));
1590 daugs[2]=EvtPDL::getId(std::string("pi+"));
1591 }else if(mode==76){
1592 _ndaugs =3;
1593 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1594 daugs[1]=EvtPDL::getId(std::string("D+"));
1595 daugs[2]=EvtPDL::getId(std::string("pi-"));
1596 }else if(mode==77){
1597 _ndaugs =3;
1598 daugs[0]=EvtPDL::getId(std::string("D0"));
1599 daugs[1]=EvtPDL::getId(std::string("D*-"));
1600 daugs[2]=EvtPDL::getId(std::string("pi+"));
1601 }else if(mode==78){
1602 _ndaugs =3;
1603 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
1604 daugs[1]=EvtPDL::getId(std::string("D*+"));
1605 daugs[2]=EvtPDL::getId(std::string("pi-"));
1606 }else if(mode==79){
1607 _ndaugs =3;
1608 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1609 daugs[1]=EvtPDL::getId(std::string("pi0"));
1610 daugs[2]=EvtPDL::getId(std::string("pi0"));
1611 }else if(mode==80){
1612 _ndaugs =2;
1613 daugs[0]=EvtPDL::getId(std::string("eta"));
1614 daugs[1]=EvtPDL::getId(std::string("J/psi"));
1615 }else if(mode==81){
1616 _ndaugs =3;
1617 daugs[0]=EvtPDL::getId(std::string("h_c"));
1618 daugs[1]=EvtPDL::getId(std::string("pi+"));
1619 daugs[2]=EvtPDL::getId(std::string("pi-"));
1620 }else if(mode==82){
1621 _ndaugs =3;
1622 daugs[0]=EvtPDL::getId(std::string("h_c"));
1623 daugs[1]=EvtPDL::getId(std::string("pi0"));
1624 daugs[2]=EvtPDL::getId(std::string("pi0"));
1625 }else if(mode==83){
1626 _ndaugs =3;
1627 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1628 daugs[1]=EvtPDL::getId(std::string("K+"));
1629 daugs[2]=EvtPDL::getId(std::string("K-"));
1630 }else if(mode==84){
1631 _ndaugs =3;
1632 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1633 daugs[1]=EvtPDL::getId(std::string("K_S0"));
1634 daugs[2]=EvtPDL::getId(std::string("K_S0"));
1635 }else if(mode==85){
1636 _ndaugs =2;
1637 daugs[0]=EvtPDL::getId(std::string("D_s*+"));
1638 daugs[1]=EvtPDL::getId(std::string("D_s*-"));
1639 }else if(mode==90){
1640 _ndaugs =3;
1641 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1642 daugs[1]=EvtPDL::getId(std::string("pi+"));
1643 daugs[2]=EvtPDL::getId(std::string("pi-"));
1644 }else if(mode==91){
1645 _ndaugs =3;
1646 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1647 daugs[1]=EvtPDL::getId(std::string("pi+"));
1648 daugs[2]=EvtPDL::getId(std::string("pi-"));
1649 }else if(mode==92){
1650 _ndaugs =3;
1651 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1652 daugs[1]=EvtPDL::getId(std::string("pi0"));
1653 daugs[2]=EvtPDL::getId(std::string("pi0"));
1654 }else if(mode==93){
1655 _ndaugs =2;
1656 daugs[0]=EvtPDL::getId(std::string("D_s+"));
1657 daugs[1]=EvtPDL::getId(std::string("D_s-"));
1658 }else if(mode==94){
1659 _ndaugs =2;
1660 daugs[0]=EvtPDL::getId(std::string("D_s*+"));
1661 daugs[1]=EvtPDL::getId(std::string("D_s-"));
1662 }else if(mode==95){
1663 _ndaugs =2;
1664 daugs[0]=EvtPDL::getId(std::string("D_s*-"));
1665 daugs[1]=EvtPDL::getId(std::string("D_s+"));
1666 }else if(mode==96){
1667 _ndaugs =2;
1668 daugs[0]=EvtPDL::getId(std::string("Lambda_c+"));
1669 daugs[1]=EvtPDL::getId(std::string("anti-Lambda_c-"));
1670 }else if(mode==111){
1671 _ndaugs =2;
1672 daugs[0]=EvtPDL::getId(std::string("Zc(4020)+"));
1673 daugs[1]=EvtPDL::getId(std::string("pi-"));
1674 }else if(mode==112){
1675 _ndaugs =2;
1676 daugs[0]=EvtPDL::getId(std::string("Zc(4020)-"));
1677 daugs[1]=EvtPDL::getId(std::string("pi+"));
1678 }else if(mode==74100){
1679 _ndaugs =2;
1680 daugs[0]=EvtPDL::getId(std::string("J/psi"));
1681 daugs[1]=EvtPDL::getId(std::string("gamma"));
1682 }else if(mode==74101){
1683 _ndaugs =2;
1684 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
1685 daugs[1]=EvtPDL::getId(std::string("gamma"));
1686 }else if(mode==74102){
1687 _ndaugs =2;
1688 daugs[0]=EvtPDL::getId(std::string("psi(3770)"));
1689 daugs[1]=EvtPDL::getId(std::string("gamma"));
1690 }else if(mode==74103){
1691 _ndaugs =2;
1692 daugs[0]=EvtPDL::getId(std::string("phi"));
1693 daugs[1]=EvtPDL::getId(std::string("gamma"));
1694 }else if(mode==74104){
1695 _ndaugs =2;
1696 daugs[0]=EvtPDL::getId(std::string("omega"));
1697 daugs[1]=EvtPDL::getId(std::string("gamma"));
1698 }else if(mode==74105){
1699 _ndaugs =2;
1700 daugs[0]=EvtPDL::getId(std::string("rho0"));
1701 daugs[1]=EvtPDL::getId(std::string("gamma"));
1702 }else if(mode==74106){
1703 _ndaugs =2;
1704 daugs[0]=EvtPDL::getId(std::string("rho(3S)0"));
1705 daugs[1]=EvtPDL::getId(std::string("gamma"));
1706 }else if(mode==74107){
1707 _ndaugs =2;
1708 daugs[0]=EvtPDL::getId(std::string("omega(2S)"));
1709 daugs[1]=EvtPDL::getId(std::string("gamma"));
1710 }else if(mode==74110){
1711 _ndaugs =2;
1712 daugs[0]=EvtPDL::getId(std::string("gamma*"));
1713 daugs[1]=EvtPDL::getId(std::string("gamma"));
1714 EvtId myvpho = EvtPDL::getId(std::string("gamma*"));
1715 EvtPDL::reSetMass(myvpho,mhds*0.9); //for calculating maxXS, it will be resent when this mode is selected
1716 _modeFlag.clear();
1717 _modeFlag.push_back(74110); //R-value generator tag
1718 _modeFlag.push_back(0); //to contain the mode selected
1719 }else if(mode==-1){
1720 _ndaugs = nson;
1721 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1722 std::cout<<"The paraent particle: "<<EvtPDL::name(pid)<<std::endl;
1723 }else if(mode==-2){
1724 _ndaugs = nson;
1725 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1726 }else if(mode==-100){
1727 _ndaugs =1;
1728 daugs[0]=EvtPDL::getId(std::string("gamma*"));
1729 _modeFlag.clear();
1730 _modeFlag.push_back(-100); //R-value generator tag
1731 _modeFlag.push_back(0); //to contain the mode selected
1732
1733 }else if(mode==10000){//use for check ISR
1734 _ndaugs =2;
1735 daugs[0]=EvtPDL::getId(std::string("pi+"));
1736 daugs[1]=EvtPDL::getId(std::string("pi-"));
1737 }else {
1738 std::cout<<"Bad mode_index number " <<mode<<", please refer to the manual."<<std::endl;
1739 ::abort();
1740 }
1741 // std::cout<<_ndaugs<<", "<<daugs[0]<<daugs[1]<<std::endl;
1742
1743 if(VISR){
1744 _ndaugs=2;
1745 daugs[0]=getDaugs()[0];//EvtPDL::getId(std::string("gamma"));
1746 daugs[1]=getDaugs()[1];//EvtPDL::getId(std::string("gamma*"));
1747 }
1748
1749 std::vector<EvtId> theDaugs;
1750 for(int i=0;i<_ndaugs;i++){
1751 theDaugs.push_back(daugs[i]);
1752 }
1753 if(theDaugs.size()) {return theDaugs;} else {std::cout<<"EvtConExc::get_mode: zero size"<<std::endl;abort();}
1754}
1755
1756
1757void EvtConExc::initProbMax(){
1758
1759 noProbMax();
1760
1761}
1762
1763void EvtConExc::decay( EvtParticle *p){
1764 //std::cout<<"_cms= "<<_cms<<" mode= "<<_mode<<std::endl;
1765 EvtId myvpho = EvtPDL::getId(std::string("vpho"));
1766 if(myvpho != p->getId()){
1767 std::cout<<"Parent needs to be vpho, but found "<<EvtPDL::name(p->getId())<<std::endl;
1768 abort();
1769 }
1770 //-- to read _cms from database or to consider beam energy spread
1771 double mvpho=EvtPDL::getMeanMass(myvpho);
1772 bool beamflag=0;
1773 if(mvpho!=EvtConExc::_cms){ //if read cms energy from database
1774 EvtConExc::_cms=mvpho;
1775 beamflag=1;
1776 }
1777
1778 /*
1779 if(beamEnergySpread!=0){
1780 cmsspread=sqrt(2.)*beamEnergySpread;
1781 _cms = energySpread(_cms,cmsspread);
1782 beamflag=1;
1783 std::cout<<"test_cms "<<_cms<<std::endl;
1784 }
1785 */
1786
1787 if(beamflag) calAF(_cms);
1788
1789
1790 //debugging
1791 //double mvpho=EvtPDL::getMeanMass(myvpho);
1792 //std::cout<<"myvpho= "<<mvpho<<std::endl;
1793
1794 //for fill root tree
1795 EvtVector4R vgam,hd1,hd2,vhds;
1796
1797 //first for Rscan generator with _mode=74110
1798 if(_mode==74110){
1799 //preparation of mode sampling
1800 std::vector<int> vmod; vmod.clear();
1801 int mn[83]={0,1,2,3,4,5,6,7,8,9,10,11,13,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,32,33,34,35,36,37,38,40,41,44,45,46, // 12,14, 30,31,39,42 and 43 is removed
1802 50,51,
1803 59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,
1804 90,91,92,93,94,95,96,
1805 74100,74101,74102,74103,74104,74105,74110};
1806 int mn2[84]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,32,33,34,35,36,37,38,40,41,44,45,46,//mode 14,30,31,42, 43 are removed
1807 50,51,
1808 59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,
1809 90,91,92,93,94,95,96,
1810 74100,74101,74102,74103,74104,74105,74110};
1811 int mn3[77]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,44,45,46,
1812 50,51,
1813 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,
1814 90,91,92,93,94,95,96,
1815 74100,74101,74102,74110}; //remove 43, 74103,74104,74105, this is included in PHOKHARA
1816 double mx = p->getP4().mass();
1817
1818 if(_cms>2.5){
1819 for(int i=0;i<83;i++){vmod.push_back(mn[i]);}
1820 }else{
1821 for(int i=0;i<84;i++){vmod.push_back(mn2[i]);}
1822 }
1823
1824 //for(int i=0;i<76;i++){vmod.push_back(mn3[i]);}
1825
1826 int mymode;
1827 double pm= EvtRandom::Flat(0.,1);
1828
1829 //std::cout<<"AF[598], AF[598] "<<AF[598]<<" "<<AF[599]<<std::endl;
1830 //--
1831 //-- where AF[599]: cross section with full region (0-Egamma^max GeV), AF[598]: cross section within Egam=(0.025-Egamma^max GeV)
1832 if(pm <_xs0/(_xs0 + _xs1)){//without ISR photon for mode=74110
1833
1834 mhds=_cms;
1835 mymode=selectMode(vmod,mhds);
1836 _selectedMode = mymode;
1837 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
1838 delete myxsection;
1839 myxsection =new EvtXsection (mymode);
1840 init_mode(mymode);
1841 resetResMass();
1842
1843 if(_ndaugs>1){//for e+e- -> exclusive decays
1844 checkA:
1845 p->makeDaughters(_ndaugs,daugs);
1846 double totMass=0;
1847 for(int i=0;i< _ndaugs;i++){
1848 EvtParticle* di=p->getDaug(i);
1849 double xmi=EvtPDL::getMass(di->getId());
1850 di->setMass(xmi);
1851 totMass+=xmi;
1852 }
1853 if(totMass>p->mass()) goto checkA;
1854 p->initializePhaseSpace( _ndaugs , daugs);
1855 if(!checkdecay(p)) goto checkA;
1856 vhds = p->getP4();
1857 if(_cms>2.5 && !angularSampling(p)) goto checkA;
1858 if(_cms<2.5 && !photonSampling(p)) goto checkA;
1859 }else{// for e+e- -> vector resonace, add a very soft photon
1860 EvtId mydaugs[2];
1861 mydaugs[0]=daugs[0];
1862 EvtPDL::reSetMass(mydaugs[0],mhds-0.001);
1863 //EvtPDL::reSetWidth(mydaugs[0],0);
1864 mydaugs[1]=gamId; //ISR photon
1865 p->makeDaughters(2,mydaugs);
1866 checkB:
1867 double totMass=0;
1868 for(int i=0;i< 2;i++){
1869 EvtParticle* di=p->getDaug(i);
1870 double xmi=EvtPDL::getMass(di->getId());
1871 di->setMass(xmi);
1872 totMass+=xmi;
1873 }
1874 //std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<" "<<checkdecay(p)<<std::endl;
1875 if(totMass>p->mass()) goto checkB;
1876 p->initializePhaseSpace(2,mydaugs);
1877
1878 if(!checkdecay(p)) goto checkB;
1879 vhds = p->getDaug(0)->getP4();;
1880 vgam = p->getDaug(1)->getP4();
1881 }
1882 //-----
1883 }else{// with ISR photon for mode=74110
1884Sampling_mhds:
1885 mhds=Mhad_sampling(MH,AF);
1886 //std::cout<<"SetMthr= "<<SetMthr<<std::endl;
1887 if(mhds<SetMthr) goto Sampling_mhds;
1888 double xeng=1-mhds*mhds/(_cms*_cms); //photon energy ratio
1889 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
1890
1891 if(mydbg) mass2=mhds;
1892
1893 //generating events
1894 ModeSelection:
1895 EvtGlobalSet::ConExcPythia=1;
1896 mymode=selectMode(vmod,mhds);
1897 if(mymode==-10) goto Sampling_mhds;
1898 conexcmode = mymode; //save for model REXC ( see EvtRexc.cc) decay
1899 if(mhds<2.3 && mymode==74110) {goto ModeSelection;} // E<2.3 GeV, fully exclusive production
1900 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
1901 _selectedMode = mymode;
1902 delete myxsection;
1903 myxsection =new EvtXsection (mymode);
1904 init_mode(mymode);
1905 EvtId myvpho = EvtPDL::getId(std::string("vgam"));
1906 EvtPDL::reSetMass(myvpho,mhds); //set to continum cms energy
1907 EvtPDL::reSetWidth(myvpho,0);
1908
1909 //debugging
1910 //for(int i=0;i<_ndaugs+1;i++){std::cout<<_ndaugs<<" "<<EvtPDL::name(daugs[i])<<std::endl;}
1911
1912 //decay e+e- ->gamma_ISR gamma*
1913 EvtId mydaugs[2];
1914 if(_ndaugs>1) { //gamma* -> exclusive decays
1915 resetResMass();
1916 mydaugs[0]=myvpho;
1917 }else{// vector resonance decays
1918 resetResMass();
1919 mydaugs[0]=daugs[0];
1920 EvtPDL::reSetMass(mydaugs[0],mhds);
1921 //EvtPDL::reSetWidth(mydaugs[0],0);
1922 }
1923 mydaugs[1]=gamId; //ISR photon
1924 int maxflag=0;
1925 int trycount=0;
1926 ISRphotonAng_sampling:
1927 double totMass=0;
1928 p->makeDaughters(2,mydaugs);
1929 for(int i=0;i< 2;i++){
1930 EvtParticle* di=p->getDaug(i);
1931 double xmi=EvtPDL::getMass(di->getId());
1932 di->setMass(xmi);
1933 totMass+=xmi;
1934 }
1935 if(totMass>p->mass()) goto ISRphotonAng_sampling;
1936 //std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
1937 double weight1 = p->initializePhaseSpace(2,mydaugs);
1938 if(!checkdecay(p)) goto ISRphotonAng_sampling;
1939 //set the ISR photon angular distribution
1940 SetP4Rvalue(p,mhds,xeng,theta); //end of decay e+e- -> vpho gamma_ISR
1941
1942 if(maxflag==0) {findMaxXS(mhds); maxflag=1;}
1943 vhds = p->getDaug(0)->getP4();
1944 vgam=p->getDaug(1)->getP4();
1945 double gx=vgam.get(1);
1946 double gy=vgam.get(2);
1947 double sintheta= sqrt(gx*gx+gy*gy)/vgam.d3mag();
1948 bool gam_ang=gam_sampling(mhds,sintheta);
1949 trycount++;
1950
1951 } // with and without ISR sampling block
1952 // finish event generation
1953 // for debugging
1954 if(mydbg){
1955 pgam[0]=vgam.get(0);
1956 pgam[1]=vgam.get(1);
1957 pgam[2]=vgam.get(2);
1958 pgam[3]=vgam.get(3);
1959
1960 phds[0]=vhds.get(0);
1961 phds[1]=vhds.get(1);
1962 phds[2]=vhds.get(2);
1963 phds[3]=vhds.get(3);
1964 costheta = vgam.get(3)/vgam.d3mag();
1965 selectmode = mymode;
1966 xs->Fill();
1967 //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
1968 }
1969 _modeFlag[1]= _selectedMode;
1970 p->setIntFlag(_modeFlag);
1971 return;
1972 }//end block if(_mode==74110)
1973
1974 // ======================================================
1975 // multi-exclusive mode
1976 //=======================================================
1977 if(_mode==-100){
1978 int mymode;
1979 double pm= EvtRandom::Flat(0.,1);
1980
1981 //std::cout<<"AF[598], AF[598] "<<AF[598]<<" "<<AF[599]<<std::endl;
1982 //--
1983 //-- where AF[599]: cross section with full region (0-Egamma^max GeV), AF[598]: cross section within Egam=(0.025-Egamma^max GeV)
1984 if(pm <_xs0/(_xs0 + _xs1)){//without ISR photon for mode=74110
1985
1986 mhds=_cms;
1987 mymode=selectMode(vmd,mhds);
1988 _selectedMode = mymode;
1989 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
1990 delete myxsection;
1991 myxsection =new EvtXsection (mymode);
1992 init_mode(mymode);
1993 resetResMass();
1994
1995 if(_ndaugs>1){//for e+e- -> exclusive decays
1996 checkAA:
1997 p->makeDaughters(_ndaugs,daugs);
1998 double totMass=0;
1999 for(int i=0;i< _ndaugs;i++){
2000 EvtParticle* di=p->getDaug(i);
2001 double xmi=EvtPDL::getMass(di->getId());
2002 di->setMass(xmi);
2003 totMass+=xmi;
2004 }
2005 if(totMass>p->mass()) goto checkAA;
2006 p->initializePhaseSpace( _ndaugs , daugs);
2007 if(!checkdecay(p)) goto checkAA;
2008 vhds = p->getP4();
2009 bool byn_ang;
2010 EvtVector4R pd1 = p->getDaug(0)->getP4();
2011 EvtVector4R pcm(_cms,0,0,0);
2012 //p->printTree();
2013 if(_ndaugs==2){//angular distribution for the hadrons only for two-body decays
2014 byn_ang=hadron_angle_sampling(pd1, pcm);
2015 if(!byn_ang) goto checkAA;
2016 }
2017
2018 }else{// for e+e- -> vector resonace, add a very soft photon
2019 EvtId mydaugs[2];
2020 mydaugs[0]=daugs[0];
2021 EvtPDL::reSetMass(mydaugs[0],mhds-0.001);
2022 //EvtPDL::reSetWidth(mydaugs[0],0);
2023 mydaugs[1]=gamId; //ISR photon
2024 p->makeDaughters(2,mydaugs);
2025 checkBA:
2026 double totMass=0;
2027 for(int i=0;i< 2;i++){
2028 EvtParticle* di=p->getDaug(i);
2029 double xmi=EvtPDL::getMass(di->getId());
2030 di->setMass(xmi);
2031 totMass+=xmi;
2032 }
2033 //std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<" "<<checkdecay(p)<<std::endl;
2034 if(totMass>p->mass()) goto checkBA;
2035 p->initializePhaseSpace(2,mydaugs);
2036
2037 if(!checkdecay(p)) goto checkBA;
2038 vhds = p->getDaug(0)->getP4();;
2039 vgam = p->getDaug(1)->getP4();
2040 }
2041 //-----
2042 }else{// with ISR photon for mode=-100
2043Sampling_mhds_A:
2044 mhds=Mhad_sampling(MH,AF);
2045 //std::cout<<"SetMthr= "<<SetMthr<<std::endl;
2046 if(mhds<SetMthr) goto Sampling_mhds_A;
2047 double xeng=1-mhds*mhds/(_cms*_cms); //photon energy ratio
2048 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
2049
2050 if(mydbg) mass2=mhds;
2051
2052 //generating events
2053 ModeSelection_A:
2054 mymode=selectMode(vmd,mhds);
2055 if(mymode==-10) goto Sampling_mhds_A;
2056 conexcmode = mymode; //save for model REXC ( see EvtRexc.cc) decay
2057 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
2058 _selectedMode = mymode;
2059 delete myxsection;
2060 myxsection =new EvtXsection (mymode);
2061 init_mode(mymode);
2062 EvtId myvpho = EvtPDL::getId(std::string("vgam"));
2063 EvtPDL::reSetMass(myvpho,mhds); //set to continum cms energy
2064 EvtPDL::reSetWidth(myvpho,0);
2065
2066 //debugging
2067 // for(int i=0;i<_ndaugs+1;i++){std::cout<<_ndaugs<<" "<<EvtPDL::name(daugs[i])<<std::endl;}
2068
2069 //decay e+e- ->gamma_ISR gamma*
2070 EvtId mydaugs[2];
2071 if(_ndaugs>1) { //gamma* -> exclusive decays
2072 resetResMass();
2073 mydaugs[0]=myvpho;
2074 }else{// vector resonance decays
2075 resetResMass();
2076 mydaugs[0]=daugs[0];
2077 EvtPDL::reSetMass(mydaugs[0],mhds);
2078 //EvtPDL::reSetWidth(mydaugs[0],0);
2079 }
2080 mydaugs[1]=gamId; //ISR photon
2081 int maxflag=0;
2082 int trycount=0;
2083 ISRphotonAng_sampling_A:
2084 double totMass=0;
2085 p->makeDaughters(2,mydaugs);
2086 for(int i=0;i< 2;i++){
2087 EvtParticle* di=p->getDaug(i);
2088 double xmi=EvtPDL::getMass(di->getId());
2089 di->setMass(xmi);
2090 totMass+=xmi;
2091 }
2092 if(totMass>p->mass()) goto ISRphotonAng_sampling_A;
2093 //std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
2094 double weight1 = p->initializePhaseSpace(2,mydaugs);
2095 if(!checkdecay(p)) goto ISRphotonAng_sampling_A;
2096 //set the ISR photon angular distribution
2097 SetP4Rvalue(p,mhds,xeng,theta); //end of decay e+e- -> vpho gamma_ISR
2098
2099 //--debugging
2100 //p->printTree();
2101
2102 if(maxflag==0) {findMaxXS(mhds); maxflag=1;}
2103 vhds = p->getDaug(0)->getP4();
2104 vgam=p->getDaug(1)->getP4();
2105 double gx=vgam.get(1);
2106 double gy=vgam.get(2);
2107 double sintheta= sqrt(gx*gx+gy*gy)/vgam.d3mag();
2108 bool gam_ang=gam_sampling(mhds,sintheta);
2109 trycount++;
2110
2111 } // with and without ISR sampling block
2112 _modeFlag[0]=-100;
2113 _modeFlag[1]= _selectedMode;
2114 p->setIntFlag(_modeFlag);
2115 // finish event generation
2116 // for debugging
2117 if(mydbg){
2118 pgam[0]=vgam.get(0);
2119 pgam[1]=vgam.get(1);
2120 pgam[2]=vgam.get(2);
2121 pgam[3]=vgam.get(3);
2122
2123 phds[0]=vhds.get(0);
2124 phds[1]=vhds.get(1);
2125 phds[2]=vhds.get(2);
2126 phds[3]=vhds.get(3);
2127 costheta = vgam.get(3)/vgam.d3mag();
2128 selectmode = mymode;
2129 xs->Fill();
2130 //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
2131 }
2132 return;
2133 }//end block if(_mode==-100)
2134
2135 //=======================================================
2136 // e+e- -> gamma_ISR gamma*
2137 //=======================================================
2138 if(VISR){
2139 EvtId gid=EvtPDL::getId("gamma*");
2140 double pm= EvtRandom::Flat(0.,1);
2141
2142 if(pm <_xs0/(_xs0 + _xs1)){//with a very soft ISR photon
2143 EvtPDL::reSetMass(gid,(_cms-0.0001) );
2144 mhds = _cms-0.0001;
2145 }else{
2146 mhds=Mhad_sampling(MH,AF);
2147 EvtPDL::reSetMass(gid,mhds);
2148 }
2149 //debugging
2150 std::cout<<"generatedMass "<<EvtPDL::getMeanMass(gid)<<std::endl;
2151 double xeng=1-mhds*mhds/(_cms*_cms); //photon energy ratio
2152 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
2153 p->makeDaughters(2,daugs);
2154 for(int i=0;i< 2;i++){
2155 EvtParticle* di=p->getDaug(i);
2156 //if(EvtPDL::name(di->getId())=="gamma*") di->setVectorSpinDensity();
2157 double xmi=EvtPDL::getMeanMass(di->getId());
2158 di->setMass(xmi);
2159 }
2160 p->initializePhaseSpace(2,daugs);
2161 SetP4(p,mhds,xeng,theta);
2162 return;
2163 }
2164
2165
2166 //========================================================
2167 //=== for other mode generation with _mode != 74110
2168 //========================================================
2169
2170 if((_xs0 + _xs1)==0) {std::cout<<"EvtConExc::zero total xsection"<<std::endl;::abort();}
2171 double pm= EvtRandom::Flat(0.,1);
2172 double xsratio = _xs0/(_xs0 + _xs1);
2173 int iflag=2; //flag = 0 for ee->hadrons, 1 for ee->gamma_ISR hadrons, 2: mix 0 and 1 case
2174 if(getArg(0)!= -2){// exclude DIY case
2175 if(getNArg()==3 && radflag==1){iflag=1;xsratio=0;} // only generating gamma hadrons mode
2176 else if(getNArg()==3 && radflag==0) {iflag=0;xsratio=1;}
2177 }
2178
2179 // std::cout<<"xsratio= "<<xsratio<<std::endl;
2180
2181 if(pm<xsratio){// no ISR photon
2182 masscheck:
2183 double summassx=0;
2184 p->makeDaughters(_ndaugs,daugs);
2185 for(int i=0;i< _ndaugs;i++){
2186 EvtParticle* di=p->getDaug(i);
2187 double xmi=EvtPDL::getMass(di->getId());
2188 di->setMass(xmi);
2189 summassx += xmi;
2190 //std::cout<<"PID and mass "<<di->getId()<<" "<<xmi<<std::endl;
2191 }
2192 if(summassx>p->mass()) goto masscheck;
2193 angular_hadron:
2194 p->initializePhaseSpace(_ndaugs,daugs);
2195 bool byn_ang;
2196 EvtVector4R pd1 = p->getDaug(0)->getP4();
2197 EvtVector4R pcm(_cms,0,0,0);
2198 if(_ndaugs==2){//angular distribution for the hadrons only for two-body decays
2199 byn_ang=hadron_angle_sampling(pd1, pcm);
2200 if(!byn_ang) goto angular_hadron;
2201 }
2202
2203 //for histogram
2204 cosp = pd1.get(3)/pd1.d3mag();
2205 mhds = _cms;
2206 //std::cout<<"cosp, mhds "<<cosp<<" "<<mhds<<std::endl;
2207 //p->printTree();
2208
2209 }else{// sampling m_hadrons and decay e+e- ->gamma gamma*
2210 double mhdr=Mhad_sampling(MH,AF);
2211 double xeng=1-mhdr*mhdr/(_cms*_cms); //photon energy ratio
2212 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
2213 EvtId gid =EvtPDL::getId("vhdr");
2214 EvtPDL::reSetMass(gid,mhdr);
2215 int ndaugs =2;
2216 EvtId mydaugs[2];
2217 mydaugs[0] =EvtPDL::getId(std::string("gamma"));
2218 mydaugs[1]=EvtPDL::getId(std::string("vhdr"));
2219
2220 //for histogram
2221 mhds = mhdr;
2222 costheta = cos(theta);
2223 //--
2224
2225 p->makeDaughters(2,mydaugs);
2226 for(int i=0;i< 2;i++){
2227 EvtParticle* di=p->getDaug(i);
2228 double xmi=EvtPDL::getMass(di->getId());
2229 //if(EvtPDL::name(di->getId())=="vhdr") di->setVectorSpinDensity();
2230 di->setMass(xmi);
2231 }
2232 p->initializePhaseSpace(2,mydaugs);
2233 SetP4(p,mhdr,xeng,theta); //end of decay e+e- -> gamma_ISR gamma*
2234 //p->printTree();
2235
2236 //----
2237 }//end of gamma gamma*, gamma*->hadrons generation
2238 // p->printTree();
2239
2240 //-----------------------------------------
2241 // End of decays for all topology
2242 //----------------------------------------
2243 //================================= event analysis
2244 if(_nevt ==0){
2245 std::cout<<"The decay chain: "<<std::endl;
2246 p->printTree();
2247 }
2248 _nevt++;
2249 //--- for debuggint to make root file
2250 if(mydbg){
2251 xs->Fill();
2252 }
2253
2254 //----
2255 return ;
2256}
2257
2258bool EvtConExc::hadron_angle_sampling(EvtVector4R ppi, EvtVector4R pcm){
2259 EvtVector4R pbst=-1*pcm;
2260 pbst.set(0,pcm.get(0));
2261 EvtVector4R p4i = boostTo(ppi,pbst);
2262 if((_mode>=0 && _mode<=5) || _mode==44 || _mode==96){//ee->two baryon mode, VP, SP, mode=-2 is excluded
2263 bool byn_ang = baryon_sampling(pcm, p4i);
2264 return byn_ang;
2265 }else if(_mode==6 || _mode==45 || _mode==46 || _mode==70 || _mode==71|| _mode==93){// ee->PP (pi+pi-,..) mode
2266 bool msn_ang = meson_sampling(pcm,p4i);
2267 return msn_ang;
2268 }else if(_mode==23 || _mode==24 || _mode==25 || _mode==26 || _mode==27 || _mode==36 || _mode==39 || _mode==47 || _mode==48 || _mode==52 || _mode==68 || _mode==69 || _mode==72 || _mode==73 || _mode==80 || _mode==94 || _mode==95){
2269 bool msn_ang = VP_sampling(pcm,p4i);
2270 return msn_ang;
2271 }else if(_mode==-2){
2272 EvtSpinType::spintype type0 = EvtPDL::getSpinType(daugs[0]);
2273 EvtSpinType::spintype type1 = EvtPDL::getSpinType(daugs[1]);
2274 if(type0==EvtSpinType::SCALAR &&type1==EvtSpinType::SCALAR){
2275 bool msn_ang = meson_sampling(pcm,p4i);
2276 return msn_ang;
2277 }else if(type0==EvtSpinType::VECTOR &&type1==EvtSpinType::SCALAR || type1==EvtSpinType::VECTOR &&type0==EvtSpinType::SCALAR){
2278 bool msn_ang = VP_sampling(pcm,p4i);
2279 return msn_ang;
2280 }
2281 }
2282 return true;
2283}
2284
2285
2286double EvtConExc::gamHXSection(EvtParticle *p,double El,double Eh,int nmc){//El, Eh : the lower and higher energy of radiative photons
2287 //std::cout<<"nmc= "<<nmc<<std::endl;
2288 gamId = EvtPDL::getId(std::string("gamma"));
2289 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
2290 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
2291 double totxs = 0;
2292 maxXS=-1;
2293 _er1=0;
2294 Rad2Xs =0;
2295 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
2296 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
2297 int gamdaugs = _ndaugs+1;
2298
2299 EvtId theDaugs[10];
2300 for(int i=0;i<_ndaugs;i++){
2301 theDaugs[i] = daugs[i];
2302 }
2303 theDaugs[_ndaugs]=gamId;
2304
2305 gamH->makeDaughters(gamdaugs,theDaugs);
2306
2307 for(int i=0;i<gamdaugs;i++){
2308 EvtParticle* di=gamH->getDaug(i);
2309 double xmi=EvtPDL::getMass(di->getId());
2310 di->setMass(xmi);
2311 }
2312 loop:
2313 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
2314 //-- slection:theta_gamma > 1 degree
2315 EvtVector4R pgam = gamH -> getDaug(_ndaugs)->getP4();
2316 double pmag = pgam.d3mag();
2317 double pz = pgam.get(3);
2318 double sintheta = sqrt( pmag*pmag - pz*pz)/pmag;
2319 double onedegree = 1./180.*3.1415926;
2320 //if(sintheta < onedegree) {goto loop;}
2321 if(pmag <El || pmag >Eh) {goto loop;}
2322 //--------
2323 // std::cout<<"pmag= "<<pmag<<std::endl;
2324
2325 double thexs = difgamXs(gamH); //prepare the photon angular distribution
2326 Rad2Xs += Rad2difXs( gamH );
2327 if(thexs>maxXS) {maxXS=thexs;}
2328 double s = p_init.mass2();
2329 double x = 2*pgam.get(0)/sqrt(s);
2330 double rad1xs = Rad1difXs(gamH); //first order xs by Eq(4)hep-ph/9910523
2331 totxs += rad1xs;
2332 _er1 += differ;
2333 gamH->deleteDaughters();
2334 } //for int_i block
2335 _er1/=nmc;
2336
2337 Rad2Xs/=nmc; // the leading second order correction
2338 totxs/=nmc; // first order correction XS
2339
2340 // return totxs; // first order correction XS
2341 return Rad2Xs; // the leading second order correction
2342}
2343
2344
2345double EvtConExc::gamHXSection(double s, double El,double Eh,int nmc){//El, Eh : the lower and higher energy of radiative photons
2346 //--for narrow resonance psi(2S),J/psi, phi, which can not calculated with the MC integration
2347 //double mxL = sqrt( s-2*Eh*sqrt(s) ); //low mass
2348 //double mxH = sqrt( s-2*El*sqrt(s) ); //high mass
2349 //double sigv = narrowRXS(mxL,mxH);
2350 //--
2351
2352 double totxs = 0;
2353 maxXS=-1;
2354 _er1=0;
2355 Rad2Xs =0;
2356 double xL=2*El/sqrt(s);
2357 double xH=2*Eh/sqrt(s);
2358 for(int i=0;i<nmc;i++){//It is found that the narrow resonance can not included in this MC integartion
2359 double rdm = EvtRandom::Flat(0.,1.);// set angular cut 1^o
2360 double x= xL+ (xH-xL)*rdm;
2361 Rad2Xs += Rad2difXs(s,x);
2362 _er1 += differ2; //stored when calculate Rad2Xs
2363 // std::cout<<"i= "<<i<<" Rad2Xs= "<<Rad2Xs<<std::endl;
2364 }
2365 _er1/=nmc;
2366 _er1*=(xH-xL);
2367 //std::cout<<"_er1= "<<_er1<<std::endl;
2368 Rad2Xs/=nmc; // the leading second order correction
2369 double thexs= Rad2Xs*(xH-xL); //xH-xL is the Jacobi factor
2370 //std::cout<<"thexs= "<<thexs<<std::endl;
2371 //if(sigv != -1) thexs += sigv; //add narrow resonance ISR cross section
2372 return thexs;
2373}
2374
2375
2376
2377double EvtConExc::gamHXSection(double El,double Eh){// using Gaussian integration subroutine from Cern lib
2378 std::cout<<" "<<std::endl;
2379 extern double Rad2difXs(double *x);
2380 extern double Rad2difXs2(double *x);
2381 double eps = 0.01;
2382 double Rad2Xs;
2383 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2384 if(Rad2Xs==0){
2385 for(int iter=0;iter<10;iter++){
2386 eps += 0.01;
2387 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2388 if(!Rad2Xs) return Rad2Xs;
2389 }
2390 }
2391 return Rad2Xs; // the leading second order correction
2392}
2393
2394
2395double EvtConExc::gamHXSection(double El,double Eh, int mode){// using Gaussian integration subroutine from Cern lib
2396 std::cout<<" "<<std::endl;
2397 extern double Rad2difXs(double *x);
2398 extern double Rad2difXs2(double *x);
2399 double eps = 0.01;
2400 double Rad2Xs;
2401 if(mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2402 if(Rad2Xs==0){
2403 for(int iter=0;iter<10;iter++){
2404 eps += 0.01;
2405 if(mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
2406 if(!Rad2Xs) return Rad2Xs;
2407 }
2408 }
2409 return Rad2Xs; // the leading second order correction
2410}
2411
2412
2413double EvtConExc::gamHXSection_er(double El,double Eh){// using Gaussian integration subroutine from Cern lib
2414 std::cout<<" "<<std::endl;
2415 extern double Rad2difXs_er(double *x);
2416 extern double Rad2difXs_er2(double *x);
2417 double eps = 0.01;
2418 double Rad2Xs;
2419 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
2420 if(Rad2Xs==0){
2421 for(int iter=0;iter<10;iter++){
2422 eps += 0.01;
2423 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
2424 if(!Rad2Xs) return Rad2Xs;;
2425 }
2426 }
2427 return Rad2Xs; // the leading second order correction
2428}
2429
2430
2431void EvtConExc::findMaxXS(EvtParticle *p){
2432 //std::cout<<"nmc= "<<nmc<<std::endl;
2433 gamId = EvtPDL::getId(std::string("gamma"));
2434 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
2435 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
2436 double totxs = 0;
2437 maxXS=-1;
2438 _er1=0;
2439 Rad2Xs =0;
2440 int nmc = 50000;
2441 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
2442 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
2443 int gamdaugs = _ndaugs+1;
2444
2445 EvtId theDaugs[10];
2446 for(int i=0;i<_ndaugs;i++){
2447 theDaugs[i] = daugs[i];
2448 }
2449 theDaugs[_ndaugs]=gamId;
2450
2451 gamH->makeDaughters(gamdaugs,theDaugs);
2452 loop:
2453 double totm=0;
2454 for(int i=0;i<gamdaugs;i++){
2455 EvtParticle* di=gamH->getDaug(i);
2456 double xmi=EvtPDL::getMass(di->getId());
2457 di->setMass(xmi);
2458 totm += xmi;
2459 }
2460
2461 //std::cout<<totm<<" "<<p_init.get(0)<<std::endl;
2462 if(totm >= p_init.get(0) ) goto loop;
2463
2464 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
2465 double thexs = difgamXs(gamH); //prepare the photon angular distribution
2466 EvtVector4R p4gam = gamH->getDaug(_ndaugs)->getP4();
2467 double costheta = p4gam.get(3)/p4gam.d3mag();
2468 double sintheta = sqrt(1-costheta*costheta);
2469 bool acut=(sintheta>0.04) && (sintheta<0.96); //set photon auglar cut 2^o
2470 if(thexs>maxXS && acut ) {maxXS=thexs;}
2471 //gamH->deleteTree();
2472 }
2473
2474}
2475
2476double EvtConExc::difgamXs(EvtParticle *p){// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
2477 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
2478 EvtVector4R p0 = p->getDaug(0)->getP4();
2479 for(int i=1;i<_ndaugs;i++){
2480 p0 += p->getDaug(i)->getP4();
2481 }
2482 EvtVector4R pgam = p->getDaug(_ndaugs)->getP4();
2483 double mhs = p0.mass();
2484 double egam = pgam.get(0);
2485 double sin2theta = pgam.get(3)/pgam.d3mag();
2486 sin2theta = 1-sin2theta*sin2theta;
2487
2488 double cms = p->getP4().mass();
2489 double alpha = 1./137.;
2490 double pi = 3.1415926;
2491 double x = 2*egam/cms;
2492 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
2493 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
2494
2495 double xs = myxsection->getXsection(mhs);
2496 double difxs = 2*mhs/cms/cms * wsx *xs;
2497 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
2498 return difxs;
2499
2500}
2501
2502
2503bool EvtConExc::gam_sampling(EvtParticle *p){//photon angular distribution sampling
2504 double pm= EvtRandom::Flat(0.,1);
2505 double xs = difgamXs( p );
2506 double xsratio = xs/maxXS;
2507 if(pm<xsratio){return true;}else{return false;}
2508}
2509
2510bool EvtConExc::gam_sampling(double mhds,double sintheta){
2511 double pm= EvtRandom::Flat(0.,1);
2512 double xs = difgamXs( mhds,sintheta );
2513 double xsratio = xs/maxXS;
2514 if(pm<xsratio){return true;}else{return false;}
2515}
2516
2517bool EvtConExc::xs_sampling(double xs){
2518 double pm= EvtRandom::Flat(0.,1);
2519 //std::cout<<"Rdm= "<<pm<<std::endl;
2520 if(pm <xs/XS_max) {return true;} else {return false;}
2521}
2522
2523bool EvtConExc::xs_sampling(double xs,double xs0){
2524 double pm= EvtRandom::Flat(0.,1);
2525 // std::cout<<"Rdm= "<<pm<<std::endl;
2526 if(pm <xs/(xs0*1.1)) {return true;} else {return false;}
2527}
2528
2529bool EvtConExc::baryon_sampling(EvtVector4R pcm, EvtVector4R pi){
2530 EvtHelSys angles(pcm,pi); //using helicity sys.angles
2531 double theta=angles.getHelAng(1);
2532 double phi =angles.getHelAng(2);
2533 double costheta=cos(theta); //using helicity angles in parent system
2534
2535 // double costh = pcm.dot(pi)/pcm.d3mag()/pi.d3mag();
2536 //std::cout<<"costheta: "<<costheta<<", "<<costh<<std::endl;
2537 double alpha = baryonAng(_cms);
2538 if(_mode ==96){alpha=-0.34;}
2539 double pm= EvtRandom::Flat(0.,1);
2540 double ang = 1 + alpha*costheta*costheta;
2541 double myang;
2542 if(alpha>=0){myang=1+alpha;}else{myang=1;}
2543 if(pm< ang/myang) {return true;}else{return false;}
2544}
2545
2546bool EvtConExc::meson_sampling(EvtVector4R pcm, EvtVector4R pi){
2547 EvtHelSys angles(pcm,pi); //using helicity sys.angles
2548 double theta=angles.getHelAng(1);
2549 double phi =angles.getHelAng(2);
2550 double costheta=cos(theta); //using helicity angles in parent system
2551
2552 double pm= EvtRandom::Flat(0.,1);
2553 double ang = 1 - costheta*costheta;
2554 if(pm< ang/1.) {return true;}else{return false;}
2555}
2556
2557bool EvtConExc::VP_sampling(EvtVector4R pcm, EvtVector4R pi){
2558 EvtHelSys angles(pcm,pi); //using helicity sys.angles
2559 double theta=angles.getHelAng(1);
2560 double phi =angles.getHelAng(2);
2561 double costheta=cos(theta); //using helicity angles in parent system
2562
2563 double pm= EvtRandom::Flat(0.,1);
2564 double ang = 1 + costheta*costheta;
2565 if(pm< ang/2.) {return true;}else{return false;}
2566}
2567
2568double EvtConExc::Rad1(double s, double x){ //first order correction
2569 //radiator correction upto second order, see Int. J. of Moder.Phys. A
2570 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
2571 double alpha = 1./137.;
2572 double pi=3.1415926;
2573 double me = 0.5*0.001; //e mass
2574 double sme = sqrt(s)/me;
2575 double L = 2*log(sme);
2576 double wsx = 2*alpha/pi/x*(L-1)*(1 - x + x*x/2 );
2577 return wsx;
2578}
2579
2580double EvtConExc::Rad2(double s, double x){
2581 //radiator correction upto second order, see Int. J. of Moder.Phys. A, hep-ph/9910523, Eq(8)
2582 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
2583 double alpha = 1./137.;
2584 double pi=3.1415926;
2585 double me = 0.5*0.001; //e mass
2586 double xi2 = 1.64493407;
2587 double xi3=1.2020569;
2588 double sme = sqrt(s)/me;
2589 double L = 2*log(sme);
2590 double beta = 2*alpha/pi*(L-1);
2591 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.;
2592 double ap = alpha/pi;
2593 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
2594 double wsx = Delta * beta * pow(x,beta-1)-0.5*beta*(2-x);
2595 double wsx2 = (2-x)*( 3*log(1-x)-4*log(x) ) -4*log(1-x)/x-6+x;
2596 wsx = wsx + beta*beta/8. *wsx2;
2597 double mymx = sqrt(s*(1-x));
2598 // return wsx*vph; //getVP(mymx);//vph is vaccum polarization factor
2599 return wsx*getVP(mymx);//vph is vaccum polarization factor
2600}
2601
2602
2603
2604double EvtConExc::Rad2difXs(EvtParticle *p){// leading second order xsection
2605 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
2606 double summass = p->getDaug(0)->getP4().mass();
2607 EvtVector4R v4p=p->getDaug(0)->getP4();
2608 for(int i=1;i<_ndaugs;i++){
2609 summass += p->getDaug(i)->getP4().mass();
2610 v4p += p->getDaug(i)->getP4();
2611 }
2612
2613 double Egam = p->getDaug(_ndaugs)->getP4().d3mag();
2614 double cms = p->getP4().mass();
2615 double mrg = cms - summass;
2616 double pm= EvtRandom::Flat(0.,1);
2617 double mhs = pm*mrg + summass;
2618
2619 double s = cms * cms;
2620 double x = 2*Egam/cms;
2621 //double mhs = sqrt( s*(1-x) );
2622 double wsx = Rad2(s,x);
2623
2624 // std::cout<<"random : "<<pm<<std::endl;
2625
2626 double xs = myxsection->getXsection(mhs);
2627 if(xs>XS_max){XS_max = xs;}
2628 double difxs = 2*mhs/cms/cms * wsx *xs;
2629 differ2 = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
2630 differ *= mrg; //Jacobi factor for variable m_hds
2631 difxs *= mrg;
2632 return difxs;
2633}
2634
2635double EvtConExc::Rad2difXs(double s, double x){// leading second order xsection
2636 double wsx = Rad2(s,x);
2637 double mhs = sqrt(s*(1-x));
2638 double xs = myxsection->getXsection(mhs);
2639 //if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
2640 if(xs>XS_max){XS_max = xs;}
2641 double difxs = wsx *xs;
2642 differ2 = wsx *(myxsection->getErr(mhs));
2643 return difxs;
2644}
2645
2646double EvtConExc::Rad2difXs(double s, double x, EvtXsection* myxsection){// leading second order xsection
2647 double wsx = Rad2(s,x);
2648 double mhs = sqrt(s*(1-x));
2649 double xs = myxsection->getXsection(mhs);
2650 //if(testflag==1) std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
2651 if(xs>XS_max){XS_max = xs;}
2652 double difxs = wsx *xs;
2653 differ2 = wsx *(myxsection->getErr(mhs));
2654 return difxs;
2655}
2656
2657
2658double EvtConExc::Rad1difXs(EvtParticle *p){// // first order xsection
2659 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
2660 double summass = p->getDaug(0)->getP4().mass();
2661 for(int i=1;i<_ndaugs;i++){
2662 summass += p->getDaug(i)->getP4().mass();
2663 }
2664
2665 double cms = p->getP4().mass();
2666 double mrg = cms - summass;
2667 double pm= EvtRandom::Flat(0.,1);
2668 double mhs = pm*mrg + summass;
2669
2670 double s = cms * cms;
2671 double x = 1 - mhs*mhs/s;
2672 double wsx = Rad1(s,x);
2673
2674 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
2675
2676 double xs = myxsection->getXsection(mhs);
2677 if(xs>XS_max){XS_max = xs;}
2678 double difxs = 2*mhs/cms/cms * wsx *xs;
2679
2680 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
2681 differ *= mrg; //Jacobi factor for variable m_hds
2682 difxs *= mrg;
2683 return difxs;
2684}
2685
2686void EvtConExc::init_Br_ee(){
2687 // double psipp_ee=9.6E-06;
2688 double psip_ee =7.73E-03;
2689 double jsi_ee =5.94*0.01;
2690 double phi_ee =2.954E-04;
2691 // double omega_ee=7.28E-05;
2692 // double rho_ee = 4.72E-05;
2693 EvtId psppId=EvtPDL::getId(std::string("psi(3770)"));
2694 EvtId pspId=EvtPDL::getId(std::string("psi(2S)"));
2695 EvtId jsiId=EvtPDL::getId(std::string("J/psi"));
2696 EvtId phiId=EvtPDL::getId(std::string("phi"));
2697 EvtId omegaId=EvtPDL::getId(std::string("omega"));
2698 EvtId rhoId=EvtPDL::getId(std::string("rho0"));
2699 BR_ee.clear();
2700 ResId.clear();
2701
2702 //BR_ee.push_back(rho_ee);
2703 //BR_ee.push_back(omega_ee);
2704 BR_ee.push_back(phi_ee);
2705 BR_ee.push_back(jsi_ee);
2706 BR_ee.push_back(psip_ee);
2707 //BR_ee.push_back(psipp_ee);
2708
2709 //ResId.push_back(rhoId);
2710 //ResId.push_back(omegaId);
2711 ResId.push_back(phiId);
2712 ResId.push_back(jsiId);
2713 ResId.push_back(pspId);
2714 //ResId.push_back(psppId);
2715
2716}
2717
2718double EvtConExc::Ros_xs(double mx,double bree, EvtId pid){// in unit nb
2719 double pi=3.1415926;
2720 double s= mx*mx;
2721 double width = EvtPDL::getWidth(pid);
2722 double mass = EvtPDL::getMeanMass(pid);
2723 double xv = 1-mass*mass/s;
2724 double sigma = 12*pi*pi*bree*width/mass/s;
2725 sigma *= Rad2(s, xv);
2726 double nbar = 3.89379304*100000; // ! GeV-2 = 3.89379304*10^5 nbar
2727 return sigma*nbar;
2728}
2729
2730
2731 double Rad2difXs(double *mhs){// leading second order xsection, mhs: the mass of final states
2732 double cms = EvtConExc::_cms;
2733 double s = cms * cms;
2734 double x = 1 - (*mhs)*(*mhs)/s;
2735 EvtConExc myconexc;
2736 double wsx;
2737 double dhs=(*mhs);
2738 double xs = EvtConExc::myxsection->getXsection(dhs);
2739 wsx=myconexc.Rad2(s,x);
2740 if(xs>EvtConExc::XS_max){EvtConExc::XS_max = xs;}
2741 double difxs = 2*dhs/cms/cms * wsx *xs;
2742 std::ofstream oa;oa<<x<<std::endl; //without this line, the running will breaking, I don't know why
2743 return difxs;
2744}
2745double Rad2difXs_er(double *mhs){// leading second order xsection, mhs: the mass of final states
2746 double cms = EvtConExc::_cms;
2747 double s = cms * cms;
2748 double x = 1 - (*mhs)*(*mhs)/s;
2749 EvtConExc myconexc;
2750 double wsx;
2751 double xs = EvtConExc::myxsection->getErr(*mhs);
2752 wsx=myconexc.Rad2(s,x);
2753 double differ2 = 2*(*mhs)/cms/cms * wsx *xs;
2754 std::ofstream oa;oa<<x<<std::endl;
2755 return differ2;
2756}
2757
2758double Rad2difXs2(double *mhs){// leading second order xsection, mhs: the mass of final states
2759 double cms = EvtConExc::_cms;
2760 double s = cms * cms;
2761 double x = 1 - (*mhs)*(*mhs)/s;
2762 EvtConExc myconexc;
2763 double wsx;
2764 double dhs=(*mhs);
2765 double xs = EvtConExc::myxsection->getXsection(dhs);
2766 wsx=myconexc.Rad2(s,x);
2767 if(xs>EvtConExc::XS_max){EvtConExc::XS_max = xs;}
2768 double difxs = 2*dhs/cms/cms * wsx *xs;
2769 oa<<x<<std::endl; //without this line, the running will breaking, I don't know why
2770 return difxs;
2771}
2772
2773double Rad2difXs_er2(double *mhs){// leading second order xsection, mhs: the mass of final states
2774 double cms = EvtConExc::_cms;
2775 double s = cms * cms;
2776 double x = 1 - (*mhs)*(*mhs)/s;
2777 EvtConExc myconexc;
2778 double wsx;
2779 double xs = EvtConExc::myxsection->getErr(*mhs);
2780 wsx=myconexc.Rad2(s,x);
2781 double differ2 = 2*(*mhs)/cms/cms * wsx *xs;
2782 oa<<x<<std::endl;
2783 return differ2;
2784}
2785
2786
2787double EvtConExc::SoftPhoton_xs(double s, double b){
2788 double alpha = 1./137.;
2789 double pi=3.1415926;
2790 double me = 0.5*0.001; //e mass
2791 double xi2 = 1.64493407;
2792 double xi3=1.2020569;
2793 double sme = sqrt(s)/me;
2794 double L = 2*log(sme);
2795 double beta = 2*alpha/pi*(L-1);
2796 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.;
2797 double ap = alpha/pi;
2798 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
2799
2800 double beta2 = beta*beta;
2801 double b2 = b*b;
2802
2803 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 -
2804 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) +
2805 16*pow(beta,2)*Li2(b))/32. ;
2806 double mymx = sqrt(s*(1-b));
2807 //return xs*vph; //getVP(mymx); //vph :the vaccum polarzation factor
2808 return xs*getVP(mymx); //vph :the vaccum polarzation factor
2809
2810}
2811
2812double EvtConExc::Li2(double x){
2813// double li2= -x +x*x/4. - x*x*x/9; // wangwp: may be not correct!
2814 double li2= +x +x*x/4. + x*x*x/9; // corrected by wangwp
2815 return li2;
2816}
2817
2818
2819double EvtConExc::lgr(double *x,double *y,int n,double t)
2820{ int n0=-1;
2821 double z;
2822 for(int i=0;i<n-1;i++){
2823 if(x[i]>=t && t> x[i+1]) {n0=i;break;}
2824 }
2825 if(n0==-1) {return 0.0;}else{
2826 double k=(y[n0]-y[n0+1])/(x[n0]-x[n0+1]);
2827 z= y[n0+1]+k*(t-x[n0+1]);
2828 }
2829 //std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
2830 //std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
2831 return z;
2832}
2833
2834bool EvtConExc::islgr(double *x,double *y,int n,double t)
2835{ int n0=-1;
2836 double z;
2837 for(int i=0;i<n-1;i++){
2838 if(x[i]>=t && t> x[i+1]) {n0=i;break;}
2839 }
2840 double xstotal=y[599];
2841 if(n0==-1) {return 0;}else{
2842 double p1=y[n0]/xstotal;
2843 double p2=y[n0+1]/xstotal;
2844 double pm= EvtRandom::Flat(0.,1);
2845 if(p1<pm && pm<=p2) {return 1;}else{return 0;}
2846 }
2847}
2848
2849double EvtConExc::LLr(double *x,double *y,int n,double t)
2850{ int n0=-1;
2851 double z;
2852 if( t==x[n-1] ) return y[n-1];
2853 for(int i=0;i<n-1;i++){
2854 if(x[i]<=t && t< x[i+1]) {n0=i;break;}
2855 }
2856 if(n0==-1) {return 0.0;}else{
2857 double k=(y[n0]-y[n0+1])/(x[n0]-x[n0+1]);
2858 z= y[n0+1]+k*(t-x[n0+1]);
2859 }
2860 //std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
2861 //std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
2862 return z;
2863}
2864
2865double EvtConExc::Mhad_sampling(double *x,double *y){//sample ISR hadrons, return Mhadron
2866 //x=MH: array contains the Mhad
2867 //y=AF: array contains the Mhad*Xsection
2868 //n: specify how many bins in the hadron sampling
2869 int n=598; //AF[599] is the total xs including Ecms point
2870 double pm= EvtRandom::Flat(0.,1);
2871 double xrt=1;
2872 int mybin=1;
2873 double xst=y[n];
2874 for(int i=0;i<n;i++){
2875 if((y[i]/xst)<pm && pm<=(y[i+1]/xst)){
2876 mybin=i+1;
2877 break;
2878 }
2879 }
2880 pm= EvtRandom::Flat(0.,1);
2881 if(pm>1){std::cout<<"random larger than 1: "<<pm<<std::endl;}
2882 double mhd=x[mybin-1]+(x[mybin]-x[mybin-1])*pm;
2883
2884 if((mhd - _cms)>0){std::cout<<"selected mhd larger than Ecms: "<<mhd<<" > "<<x[mybin] <<std::endl;abort();}
2885 if(mhd<_mhdL){
2886 std::cout<<"the sample mhassrons is less than the low end of XS"<<mhd<<" <"<<_mhdL<<std::endl;
2887 for(int i=0;i<598;i++){std::cout<<i<<" "<<x[i]<<" "<<y[i]<<std::endl;}
2888 abort();
2889 }
2890 return mhd;
2891}
2892
2893double EvtConExc::ISR_ang_integrate(double x,double theta){
2894 //see Eq(11) in arXif:hep-ph/9910523, return h(theta)/h(pi)
2895 double costheta=cos(theta);
2896 double eb=_cms/2;
2897 double cos2=costheta*costheta;
2898 double sin2=1-cos2;
2899 double me=0.51*0.001;
2900 double L=2*log(_cms/me);
2901 double meE2= me*me/eb/eb;
2902 double hpi=L-1;
2903 double hq1= meE2/2*costheta/(sin2+meE2*cos2);
2904 double hq2= -0.5*log((1+sqrt(1-meE2)*costheta)/(1-sqrt(1-meE2)*costheta));
2905 double hq3=x*x*costheta/2/(x*x-2*x+2)*(1-meE2/(sin2+meE2*cos2));
2906 double hq=(L-1)/2.+hq1+hq2+hq3;
2907 hq /= hpi; //Eq (11) in arXif:hep-ph/9910523
2908 return hq;
2909}
2910
2911double EvtConExc::ISR_ang_sampling(double x){
2912 double xx[180],yy[180];
2913 double dgr2rad=1./180.*3.1415926;
2914 xx[0]=0;
2915 xx[1]=5*dgr2rad; //first bin is 5 degree
2916 int nc=2;
2917 for(int i=6;i<=175;i++){ //last bin is 5 degree
2918 xx[nc]=i*dgr2rad;
2919 nc++;
2920 }
2921 xx[nc]=180*dgr2rad;
2922 for(int j=0;j<=nc;j++){
2923 yy[j]=ISR_ang_integrate(x,xx[j]);
2924 }
2925 double pm= EvtRandom::Flat(0.,1);
2926 int mypt=1;
2927 for(int j=1;j<=nc;j++){
2928 if(yy[j-1]<pm && pm<=yy[j]){mypt=j;break;}
2929 }
2930 pm= EvtRandom::Flat(0.,1);
2931 double mytheta= xx[mypt-1]+ (xx[mypt]-xx[mypt-1])*pm;
2932 return mytheta; //theta in rad unit
2933}
2934
2935void EvtConExc::SetP4(EvtParticle *p, double mhdr, double xeng,double theta){ //set the gamma and gamma* momentum according sampled results
2936 double pm= EvtRandom::Flat(0.,1);
2937 double phi=2*3.1415926*pm;
2938 double gamE = _cms/2*xeng;
2939 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
2940 double px = gamE*sin(theta)*cos(phi);
2941 double py = gamE*sin(theta)*sin(phi);
2942 double pz = gamE*cos(theta);
2943 EvtVector4R p4[2];
2944 p4[0].set(gamE,px,py,pz);//ISR photon
2945 p4[1].set(hdrE,-px,-py,-pz); //mhdraons
2946
2947 for(int i=0;i<2;i++){
2948 EvtId myid = p->getDaug(i)->getId();
2949 p->getDaug(i)->init(myid,p4[i]);
2950 if(EvtPDL::name(myid)=="gamma*"){
2951 if( (_cms-mhdr)<0.002){
2952 EvtVector4R PG(_cms,0,0,0);
2953 p->getDaug(i)->init(myid,PG);
2954 }
2955 }
2956 }
2957}
2958
2959void EvtConExc::SetP4Rvalue(EvtParticle *p, double mhdr, double xeng,double theta){ //set the gamma and gamma* momentum according sampled results
2960 double pm= EvtRandom::Flat(0.,1);
2961 double phi=2*3.1415926*pm;
2962 double gamE = _cms/2*xeng;
2963 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
2964 double px = gamE*sin(theta)*cos(phi);
2965 double py = gamE*sin(theta)*sin(phi);
2966 double pz = gamE*cos(theta);
2967 EvtVector4R p4[2];
2968 p4[0].set(hdrE,px,py,pz);//mhdraons
2969 p4[1].set(gamE,-px,-py,-pz); //ISR photon
2970 for(int i=0;i<2;i++){
2971 EvtId myid = p->getDaug(i)->getId();
2972 p->getDaug(i)->init(myid,p4[i]);
2973 }
2974}
2975
2976
2977void EvtConExc::findMaxXS(double mhds){// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
2978 maxXS=-1;
2979 for(int i=0;i<90000;i++){
2980 double x = 1-(mhds/_cms)*(mhds/_cms);
2981 double cos = EvtRandom::Flat(0.0006,0.9994); //set cut on ISR gamma 2^o
2982 double sin2theta =sqrt(1-cos*cos);
2983 double alpha = 1./137.;
2984 double pi = 3.1415926;
2985 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
2986 double xs = myxsection->getXsection(mhds);
2987 double difxs = 2*mhds/_cms/_cms * wsx *xs;
2988 if(difxs>maxXS)maxXS=difxs;
2989 }
2990}
2991
2992double EvtConExc::difgamXs(double mhds,double sintheta){// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
2993 double x = 1-(mhds/_cms)*(mhds/_cms);
2994 double sin2theta = sintheta*sintheta;
2995 double alpha = 1./137.;
2996 double pi = 3.1415926;
2997 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
2998 double xs = myxsection->getXsection(mhds);
2999 double difxs = 2*mhds/_cms/_cms * wsx *xs;
3000 return difxs;
3001}
3002
3003int EvtConExc::selectMode(std::vector<int> vmod, double mhds){
3004 //first get xsection for each mode in vmod array
3005 //std::cout<<"vmod.size, mhds "<<vmod.size()<<" "<<mhds<<std::endl;
3006 std::vector<int>excmod;
3007 excmod.push_back(0);
3008 excmod.push_back(1);
3009 excmod.push_back(2);
3010 excmod.push_back(6);
3011 excmod.push_back(7);
3012 excmod.push_back(12);
3013 excmod.push_back(13);
3014 excmod.push_back(45);
3015 excmod.push_back(46);
3016 double uscale = 1;
3017 std::vector<double> vxs;vxs.clear();
3018 for (int i=0;i<vmod.size();i++){
3019 int imod = vmod[i];
3020 delete myxsection;
3021 myxsection =new EvtXsection (imod);
3022 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
3023 //if(uscale!=1) std::cout<<"mode="<<imod<<" uscale= "<<uscale<<std::endl;
3024
3025 double myxs=uscale*myxsection->getXsection(mhds);
3026
3027 if(i==0) {vxs.push_back(myxs);}
3028 else if(imod==74110){//for continuum process
3029 double xcon = myxs - vxs[i-1];
3030 if(xcon<0) {xcon=vxs[i-1];}else{xcon=myxs;}
3031 if(mhds<2.0) xcon=vxs[i-1]; //continuum cut: _cms energy less than 1.1, sampling phi,rho0 and omega resonance, see parp(2)=1.08 at Pythia.F
3032 vxs.push_back(xcon);
3033 }else{
3034 vxs.push_back(myxs+vxs[i-1]);
3035 }
3036 }//for_i_block
3037 //degugging
3038 // for(int i=0;i<vxs.size();i++){std::cout<<"Mhds="<<mhds<<" Mode="<<vmod[i]<<" xs_i= "<<vxs[i]<<" totalxs "<< vxs[vxs.size()-1]<<std::endl;}
3039
3040
3041 double totxs = vxs[vxs.size()-1];
3042 //if(totxs==0){std::cout<<"totalXs=0, vxs.size()= "<<vxs.size()<<std::endl;}
3043 int icount=0;
3044 int themode;
3045 mode_iter:
3046 double pm= EvtRandom::Flat(0.,1); //std::cout<<"pm= "<<pm<<" mhds= "<<mhds<<std::endl;
3047 if(pm <= vxs[0]/totxs) {
3048 themode= vmod[0];
3049 std::vector<EvtId> theDaug=get_mode(themode);
3050 EvtId daugs[100];
3051 for(int i=0;i<theDaug.size();i++){daugs[i]=theDaug[i];}
3052 int exmode=EvtDecayTable::inChannelList(theparent->getId(),theDaug.size(),daugs);
3053 if(_mode!=74110){return themode;}
3054 else if(exmode==-1){ return themode;}else{goto mycount;}
3055 }
3056
3057 for(int j=1;j<vxs.size();j++){
3058 double x0 = vxs[j-1]/totxs;
3059 double x1 = vxs[j]/totxs; //std::cout<<"j,x0,x1 "<<j<<" "<<x0<<" "<<x1<<std::endl;
3060 if(x0<pm && pm<=x1){
3061 themode=vmod[j];
3062 std::vector<EvtId> theDaug=get_mode(themode);
3063 EvtId daugs[100];
3064 for(int i=0;i<theDaug.size();i++){daugs[i]=theDaug[i];}
3065 int exmode=EvtDecayTable::inChannelList(theparent->getId(),theDaug.size(),daugs);
3066 if(_mode!=74110){return themode;}
3067 else if(exmode==-1){ return themode;} else{ break;}
3068 }
3069 }
3070
3071 mycount:
3072 icount++;
3073 if(icount<20000 ) goto mode_iter;
3074 //debugging
3075 //for(int i=0;i<vxs.size();i++){std::cout<<"Random="<<pm<<" Mode,Mhad "<<vmod[i]<<", "<<mhds<<" xs_i "<<vxs[i]<<" xsi/totalxs "<<vxs[i]/totxs<<std::endl;}
3076 std::cout<<"EvtConExc::selectMode can not find a mode with 20000 iteration for Mhadrons="<<mhds<<std::endl;
3077 return -10;
3078 //abort();
3079
3080}
3081
3082
3083void EvtConExc::resetResMass(){
3084 // EvtGen base class resetwidth has bugs, it make the mass of particle changed.
3085 EvtId myres = EvtPDL::getId(std::string("J/psi"));
3086 EvtPDL::reSetMass(myres,mjsi);
3087 //EvtPDL::reSetWidth(myres,wjsi);
3088
3089 myres = EvtPDL::getId(std::string("psi(2S)"));
3090 EvtPDL::reSetMass(myres,mpsip);
3091 //EvtPDL::reSetWidth(myres,wpsip);
3092
3093 myres = EvtPDL::getId(std::string("psi(3770)"));
3094 EvtPDL::reSetMass(myres,mpsipp);
3095 //EvtPDL::reSetWidth(myres,wpsipp);
3096
3097 myres = EvtPDL::getId(std::string("phi"));
3098 EvtPDL::reSetMass(myres,mphi);
3099 //EvtPDL::reSetWidth(myres,wphi);
3100
3101 myres = EvtPDL::getId(std::string("omega"));
3102 EvtPDL::reSetMass(myres,momega);
3103 //EvtPDL::reSetWidth(myres,womega);
3104
3105 myres = EvtPDL::getId(std::string("rho0"));
3106 EvtPDL::reSetMass(myres,mrho0);
3107 //EvtPDL::reSetWidth(myres,wrho0);
3108
3109 myres = EvtPDL::getId(std::string("rho(3S)0"));
3110 EvtPDL::reSetMass(myres,mrho3s);
3111 //EvtPDL::reSetWidth(myres,wrho3s);
3112
3113 myres = EvtPDL::getId(std::string("omega(2S)"));
3114 EvtPDL::reSetMass(myres,momega2s);
3115 //EvtPDL::reSetWidth(myres,womega2s);
3116}
3117
3118void EvtConExc::getResMass(){
3119 EvtId myres = EvtPDL::getId(std::string("J/psi"));
3120 mjsi = EvtPDL::getMeanMass(myres);
3121 wjsi = EvtPDL::getWidth(myres);
3122
3123 myres = EvtPDL::getId(std::string("psi(2S)"));
3124 mpsip= EvtPDL::getMeanMass(myres);
3125 wpsip= EvtPDL::getWidth(myres);
3126
3127 myres = EvtPDL::getId(std::string("psi(3770)"));
3128 mpsipp= EvtPDL::getMeanMass(myres);
3129 wpsipp= EvtPDL::getWidth(myres);
3130
3131 myres = EvtPDL::getId(std::string("phi"));
3132 mphi = EvtPDL::getMeanMass(myres);
3133 wphi = EvtPDL::getWidth(myres);
3134
3135 myres = EvtPDL::getId(std::string("omega"));
3136 momega= EvtPDL::getMeanMass(myres);
3137 womega= EvtPDL::getWidth(myres);
3138
3139 myres = EvtPDL::getId(std::string("rho0"));
3140 mrho0 = EvtPDL::getMeanMass(myres);
3141 wrho0 = EvtPDL::getWidth(myres);
3142
3143 myres = EvtPDL::getId(std::string("rho(3S)0"));
3144 mrho3s= EvtPDL::getMeanMass(myres);
3145 wrho3s= EvtPDL::getWidth(myres);
3146
3147
3148 myres = EvtPDL::getId(std::string("omega(2S)"));
3149 momega2s= EvtPDL::getMeanMass(myres);
3150 womega2s= EvtPDL::getWidth(myres);
3151
3152}
3153
3154void EvtConExc::showResMass(){
3155 std::cout<<"J/psi: "<<mjsi<<" "<<wjsi<<std::endl;
3156 std::cout<<"psipr: "<<mpsip<<" "<<wpsip<<std::endl;
3157 std::cout<<"psipp: "<<mpsipp<<" "<<wpsipp<<std::endl;
3158 std::cout<<"phi : "<<mphi<<" "<<wphi<<std::endl;
3159 std::cout<<"omega: "<<momega<<" "<<womega<<std::endl;
3160 std::cout<<"rho0 : "<<mrho0<<" "<<wrho0<<std::endl;
3161 std::cout<<"rho(3S)0: "<<mrho3s<<" "<<wrho3s<<std::endl;
3162 std::cout<<"omega(2S): "<<momega2s<<" "<<womega2s<<std::endl;
3163}
3164
3165bool EvtConExc::checkdecay(EvtParticle* p){
3166 int nson=p->getNDaug();
3167 double massflag=1;
3168 for(int i=0;i<nson;i++){
3169 if( EvtPDL::getStdHep(p->getDaug(i)->getId())==22 ) continue;
3170 massflag *= p->getDaug(i)->mass();
3171 }
3172 if(massflag==0){
3173 std::cout<<"Zero mass!"<<std::endl;
3174 return 0;
3175 }else{return 1;}
3176}
3177
3178
3179double EvtConExc::sumExc(double mx){//summation all cross section of exclusive decays
3180 std::vector<int> vmod; vmod.clear();
3181 int mn[83]={0,1,2,3,4,5,6,7,8,9,10,11,13,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,32,33,34,35,36,37,38,40,41,44,45,46, // 12,14, 30,31,39,42 and 43 is removed
3182 50,51,
3183 59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,
3184 90,91,92,93,94,95,96,
3185 74100,74101,74102,74103,74104,74105,74110};
3186 int mn2[84]={0,1,2,3,4,5,6,7,8,9,10,11,12,13,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,32,33,34,35,36,37,38,40,41,44,45,46,//mode 14,30,31,42, 43 are removed
3187 50,51,
3188 59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,
3189 90,91,92,93,94,95,96,
3190 74100,74101,74102,74103,74104,74105,74110};
3191
3192 if(_cms > 2.5){
3193 for(int i=0;i<83;i++){vmod.push_back(mn[i]);}
3194 }else{
3195 for(int i=0;i<84;i++){vmod.push_back(mn2[i]);}
3196 }
3197
3198 // for(int i=0;i<79;i++){vmod.push_back(mn2[i]);}
3199
3200 double xsum=0;
3201 double uscale = 1;
3202 for(int i=0;i<vmod.size();i++){
3203 int mymode = vmod[i];
3204 if(mymode ==74110) continue;
3205 delete myxsection;
3206 myxsection =new EvtXsection (mymode);
3207 init_mode(mymode);
3208 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
3209 xsum += uscale*myxsection->getXsection(mx);
3210 //if(mx==3.0967) {std::cout<<vmod[i]<<" "<<uscale<<" "<<xsum<<std::endl;}
3211 }
3212 return xsum;
3213}
3214
3215bool EvtConExc::angularSampling(EvtParticle* par){
3216 bool tagp,tagk;
3217 tagk=0;
3218 tagp=0;
3219 int nds = par->getNDaug();
3220 for(int i=0;i<par->getNDaug();i++){
3221 EvtId di=par->getDaug(i)->getId();
3222 EvtVector4R p4i = par->getDaug(i)->getP4Lab();
3223 int pdgcode =EvtPDL::getStdHep( di );
3224 double alpha=1;
3225 /*
3226 if(fabs(pdgcode)==2212){//proton and anti-proton
3227 alpha = baryonAng(p4i.mass());
3228 cosp = cos(p4i.theta());
3229 tagp=1;
3230 }else if(fabs(pdgcode)==321||fabs(pdgcode)==211){
3231 alpha=1;
3232 cosk = cos(p4i.theta());
3233 tagk=1;
3234 }else continue;
3235 */
3236 double angmax = 1+alpha;
3237 double costheta = cos(p4i.theta());
3238 double ang=1+alpha*costheta*costheta;
3239 double xratio = ang/angmax;
3240 double xi=EvtRandom::Flat(0.,1);
3241 //std::cout<<"pdgcode "<<pdgcode<<std::endl;
3242 //std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
3243 if(xi>xratio) return false;
3244 }//loop over duaghters
3245 //if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
3246 //if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
3247
3248 return true;
3249}
3250
3251double EvtConExc::baryonAng(double mx){
3252 double mp=0.938;
3253 double u = 0.938/mx;
3254 u = u*u;
3255 double u2 = (1+u)*(1+u);
3256 double uu = u*(1+6*u);
3257 double alpha = (u2-uu)/(u2+uu);
3258 return alpha;
3259}
3260
3261bool EvtConExc::photonSampling(EvtParticle* par){
3262 bool tagp,tagk;
3263 tagk=0;
3264 tagp=0;
3265 int nds = par->getNDaug();
3266 for(int i=0;i<par->getNDaug();i++){
3267 EvtId di=par->getDaug(i)->getId();
3268 EvtVector4R p4i = par->getDaug(i)->getP4Lab();
3269 int pdgcode =EvtPDL::getStdHep( di );
3270 double alpha=0;
3271
3272 if(pdgcode==111 ||pdgcode==221 || pdgcode==331){//for photon
3273 alpha = 1;
3274 }else continue;
3275
3276 double angmax = 1+alpha;
3277 double costheta = cos(p4i.theta());
3278 double ang=1+alpha*costheta*costheta;
3279 double xratio = ang/angmax;
3280 double xi=EvtRandom::Flat(0.,1);
3281 //std::cout<<"pdgcode "<<pdgcode<<std::endl;
3282 //std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
3283 if(xi>xratio) return false;
3284 }//loop over duaghters
3285 //if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
3286 //if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
3287
3288 return true;
3289}
3290
3291
3292double EvtConExc::narrowRXS(double mxL,double mxH){
3293 //br for ee
3294 double psippee,psipee,jsiee,phiee,omegaee,rhoee;
3295 double kev2Gev=0.000001;
3296 psippee=0.262*kev2Gev;
3297 psipee =2.36*kev2Gev;
3298 jsiee =5.55*kev2Gev;
3299 phiee=4.266*0.001*2.954*0.0001;
3300 omegaee=0.6*kev2Gev;
3301 rhoee=7.04*kev2Gev;
3302 double s=_cms*_cms;
3303 double sigv=0;
3304 double mv=0;
3305 double widee=0;
3306 double xpi=12*3.1415926*3.1415926;
3307 if(mxL<=3.0957 && 3.0957<=mxH || mxL<=3.0983 && 3.0983<=mxH) {
3308 widee = jsiee;
3309 mv = 3.096916;
3310 }else if(mxL<=3.6847 && 3.6847<=mxH || mxL<=3.6873 && 3.6873<=mxH){
3311 widee = psipee;
3312 mv = 3.686109;
3313 }else if(mxL<=1.01906 && 1.01906<=mxH || mxL<= 1.01986 && 1.01986<=mxH){
3314 widee = phiee;
3315 mv = 1.01946;
3316 }else{return -1.0;}
3317
3318 if(_cms<=mv) return -1.;
3319 double xv = 1-mv*mv/s;
3320 sigv += xpi*widee/mv/s*Rad2(s,xv);
3321 double unic = 3.89379304*100000; //in unit nb
3322 return sigv*unic; //vaccum factor has included in the Rad2
3323}
3324
3325
3326double EvtConExc::addNarrowRXS(double mhi,double binwidth){
3327 double myxs = 0;
3328 for(int i=0;i<ISRXS.size();i++){ // std::cout<<"ISRM "<<ISRM[i]<<std::endl;
3329 if( (ISRM[i]-mhi)<binwidth && ISRFLAG[i]){ISRFLAG[i]=0; myxs = ISRXS[i];}
3330 }
3331 //std::cout<<mhi<<" "<<binwidth<<" "<<myxs<<std::endl;
3332 return myxs;
3333}
3334
3335double EvtConExc::selectMass(){
3336 double pm,mhdL,mhds;
3337 pm = EvtRandom::Flat(0.,1);
3338 mhdL =_mhdL;
3339 mhds = pm*(_cms - mhdL)+mhdL; //non narrow resonance
3340 std::vector<double> sxs;
3341 for(int i=0;i<ISRID.size();i++){
3342 double ri=ISRXS[i]/AF[598];
3343 sxs.push_back(ri);
3344 }
3345 for(int j=0;j<sxs.size();j++){
3346 if(j>0) sxs[j] += sxs[j-1];
3347 }
3348 pm = EvtRandom::Flat(0.,1);
3349 if(pm<sxs[0]) {
3350 mhds = EvtPDL::getMass(ISRID[0]); //narrow resonance
3351 }
3352 for(int j=1;j<sxs.size();j++){
3353 double x0 = sxs[j-1];
3354 double x1 = sxs[j];
3355 if(x0<pm && pm<=x1) mhds = EvtPDL::getMass(ISRID[j]); //narrow resonance
3356 }
3357 return mhds;
3358}
3359
3360double EvtConExc::getVP(double mx){
3361 double xx,r1,i1;
3362 double x1,y1,y2;
3363 xx=0;
3364 int mytag=1;
3365 for(int i=0;i<4001;i++){
3366 x1=vpx[i];
3367 y1=vpr[i];
3368 y2=vpi[i];
3369 if(xx<=mx && mx <x1){xx=x1; r1=y1; i1=y2; mytag=0; break;}
3370 xx=x1; r1=y1; i1=y2;
3371 }
3372 if(mytag==1){std::cout<<"No vacuum polarization value found"<<std::endl;abort();}
3373 EvtComplex cvp(r1,i1);
3374 double thevp=abs(1./(1-cvp));
3375 if(3.0933<mx && mx<3.1035) return 1.; //J/psi peak excluded
3376 if(3.6810<mx && mx<3.6913) return 1.; //psi(2S) peak excluded
3377 return thevp*thevp;
3378}
3379
3380
3381void EvtConExc::ReadVP(){
3382 vpx.clear();vpr.clear();vpi.clear();
3383 int vpsize=4001;
3384 vpx.resize(vpsize);
3385 vpr.resize(vpsize);
3386 vpi.resize(vpsize);
3387 std::string locvp=getenv("BESEVTGENROOT");
3388 locvp +="/share/vp.dat";
3389 ifstream m_inputFile;
3390 m_inputFile.open(locvp.c_str());
3391 double xx,r1,i1;
3392 double x1,y1,y2;
3393 xx=0;
3394 int mytag=1;
3395 for(int i=0;i<4001;i++){
3396 m_inputFile>>vpx[i]>>vpr[i]>>vpi[i];
3397 }
3398
3399}
3400
3401void EvtConExc::mk_VXS(double Esig, double Egamcut, double EgamH, int midx){
3402 //the mode index is put into vmode
3403 //initial
3404 //midx: mode index in vmode[midx] and VXS[midx][bin]
3405 int mode=vmode[midx];
3406 double uscale;
3407
3408 EvtXsection* myxsection = new EvtXsection (mode);
3409 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
3410 double s = _cms*_cms;
3411 double mlow=myxsection->getXlw();
3412 if(_cms <= mlow){
3413 for(int i=0;i<600;i++){VXS[midx][i]=0;}
3414 return;
3415 }
3416
3417 int nint=50;
3418 int nmc=800;
3419 //double myxs0 =uscale*gamHXSection(s,Esig,Egamcut,nmc);
3420 double myxs0 = uscale*trapezoid(s,Esig,Egamcut,nint,myxsection);
3421 double xb= 2*Esig/_cms;
3422 double sig_xs = uscale*SoftPhoton_xs(s,xb)*(myxsection->getXsection(_cms));
3423 myxs0 += sig_xs;
3424
3425 int Nstp=598;
3426 double stp=(EgamH - Egamcut)/Nstp;
3427 for(int i=0;i<Nstp;i++){//points within Egam(max) ~ Egamcut=25MeV
3428 double Eg0=EgamH - i*stp;
3429 double Eg1=EgamH - (i+1)*stp;
3430 //double xsi= uscale*gamHXSection(s,Eg1,Eg0,nmc);
3431 double xsi= uscale*trapezoid(s,Eg1,Eg0,nint,myxsection);
3432 if(i==0) VXS[midx][0]=xsi;
3433 if(i>=1) VXS[midx][i]=VXS[midx][i-1]+xsi;
3434 }
3435 VXS[midx][598]=VXS[midx][597];
3436 VXS[midx][599]=VXS[midx][598]+ myxs0;
3437 //std::cout<<"mode "<<mode<<" "<<uscale<<" "<<VXS[midx][599]<<std::endl;
3438}
3439
3440int EvtConExc::get_mode_index(int mode){
3441 int i=0;
3442 for(int j=0;i<vmode.size();j++){
3443 if(mode==vmode[j]) return j;
3444 }
3445 std::cout<<" EvtConExc::get_mode_index: no index is found in vmode"<<std::endl;
3446 abort();
3447}
3448
3449double EvtConExc::getObsXsection(double mhds,int mode){
3450 double hwid=(AA[0]-AA[1])/2.;
3451 double s=_cms*_cms;
3452 int idx=get_mode_index(mode);
3453 double Egam=(s-mhds*mhds)/2/_cms;
3454 for(int i=0;i<600;i++){
3455 if( (AA[i]-hwid)<=Egam && Egam<(AA[i]+hwid) ) return VXS[idx][i];
3456 }
3457
3458 std:cout<<"EvtConExc::getObsXsection : no observed xs is found "<<endl;
3459 abort();
3460}
3461
3462double EvtConExc::Egam2Mhds(double Egam){
3463 double s=_cms*_cms;
3464 double mhds = sqrt( s - 2*Egam*_cms);
3465 return mhds;
3466}
3467
3468void EvtConExc::writeDecayTabel(){
3469 ofstream oa,ob;
3470 oa.open("_pkhr.dec");
3471 ob.open("obsxs.dat");
3472 //
3473 int im=getModeIndex(74110);
3474
3475 double xscon=VXS[im][599];
3476 //std::cout<<"tsize "<<tsize<<" "<<xscon<<" "<<VXS[70][599]<<std::endl;
3477 std::vector<int> Vmode;
3478 Vmode.push_back(6);//1:pi+pi-
3479 Vmode.push_back(13);//2:pi+pi-2pi0
3480 Vmode.push_back(12);//3:2(pi+pi-)
3481 Vmode.push_back(0);//4:ppbar
3482 Vmode.push_back(1);//5:nnbar
3483 Vmode.push_back(45);//6:K+K-
3484 Vmode.push_back(46);//7:K0K0bar
3485 Vmode.push_back(7);//8:pi+pi-pi0
3486 Vmode.push_back(2);//9:Lambda anti-Lambda
3487 Vmode.push_back(37);//10: pi+pi- eta
3488 std::vector<std::string> vmdl;
3489 vmdl.push_back("PHOKHARA_pipi");
3490 vmdl.push_back("PHOKHARA_pi0pi0pipi");
3491 vmdl.push_back("PHOKHARA_4pi");
3492 vmdl.push_back("PHOKHARA_ppbar");
3493 vmdl.push_back("PHOKHARA_nnbar");
3494 vmdl.push_back("PHOKHARA_KK");
3495 vmdl.push_back("PHOKHARA_K0K0");
3496 vmdl.push_back("PHOKHARA_pipipi0");
3497 vmdl.push_back("PHOKHARA_LLB");
3498 vmdl.push_back("PHOKHARA_pipieta");
3499
3500 ob<<"mode_index observed cross /nb"<<std::endl;
3501 for(int i=0;i<vmode.size();i++){
3502 ob<<vmode[i]<<" "<<VXS[getModeIndex(vmode[i])][599]<<std::endl;
3503//cout<<"vmode["<<i<<"] = "<<vmode[i]<<", VXS["<<getModeIndex(vmode[i])<<"][599] = "<<XS[getModeIndex(vmode[i])][599]<<std::endl; // wangwp
3504 }
3505
3506 for(int i=0;i<Vmode.size();i++){
3507 //std::cout<<"i: "<<xscon<<" "<<VXS[Vmode[i]][599]<<std::endl;
3508 xscon -= VXS[ getModeIndex(Vmode[i]) ][599];
3509 }
3510 //debugging
3511 for(int i=0;i<Vmode.size();i++){
3512 // std::cout<<Vmode[i]<<"-th mode: "<<VXS[getModeIndex(Vmode[i])][599]<<" "<<std::endl;
3513 }
3514
3515 oa<<"LundaPar PARJ(11)=0.611798"<<std::endl;
3516 oa<<"LundaPar PARJ(12)=7.92527E-12"<<std::endl;
3517 oa<<"LundaPar PARJ(14)=0.244495"<<std::endl;
3518 oa<<"LundaPar PARJ(15)=1.16573E-15"<<std::endl;
3519 oa<<"LundaPar PARJ(16)=0.436516"<<std::endl;
3520 oa<<"LundaPar PARJ(17)=0.530517"<<std::endl;
3521 oa<<"LundaPar PARJ(1)=0.0651577"<<std::endl;
3522 oa<<"LundaPar PARJ(2)=0.260378"<<std::endl;
3523 oa<<"LundaPar PARJ(21)=0.0664835"<<std::endl;
3524 oa<<"LundaPar RALPA(15)=0.576687"<<std::endl;
3525 oa<<"LundaPar RALPA(16)=0.364796"<<std::endl;
3526 oa<<"LundaPar RALPA(17)=3.19486E-06"<<std::endl;
3527 oa<<"noPhotos"<<std::endl;
3528 oa<<"Particle vpho "<<_cms<<" 0"<<std::endl;
3529 oa<<"Decay vpho"<<std::endl;
3530 oa<<"0 pi+ pi- PHSP ;"<<std::endl;
3531 oa<<"0 pi0 pi0 pi- pi+ PHSP ;"<<std::endl;
3532 oa<<"0 pi+ pi- pi- pi+ PHSP ;"<<std::endl;
3533 oa<<"0 anti-p- p+ PHSP ;"<<std::endl;
3534 oa<<"0 anti-n0 n0 PHSP ;"<<std::endl;
3535 oa<<"0 K+ K- PHSP ;"<<std::endl;
3536 oa<<"0 K_S0 K_L0 PHSP ;"<<std::endl;
3537 oa<<"0 pi+ pi- pi0 PHSP ;"<<std::endl;
3538 oa<<"0 Lambda0 anti-Lambda0 PHSP ;"<<std::endl;
3539 oa<<"0 pi+ pi- eta PHSP ;"<<std::endl;
3540 oa<<"0 gamma phi PHSP;"<<std::endl;
3541 oa<<"0 gamma rho0 PHSP;"<<std::endl;
3542 oa<<"0 gamma omega PHSP;"<<std::endl;
3543 oa<<xscon<<" ConExc 74110;"<<std::endl;
3544 for(int i=0;i<Vmode.size();i++){
3545 oa<<VXS[ getModeIndex(Vmode[i]) ][599]<<" "<<vmdl[i]<<" ;"<<std::endl;
3546 }
3547 oa<<"Enddecay"<<std::endl;
3548 oa<<"End"<<std::endl;
3549}
3550
3551
3552
3553void EvtConExc::checkEvtRatio(){
3554 ofstream oa;
3555 oa.open("_evtR.dat");
3556 //
3557 int im=getModeIndex(74110);
3558 double xscon=VXS[im][599];
3559 double xscon0=xscon;
3560 oa<<"Ecms= "<<_cms<<" GeV"<<std::endl;
3561 oa<<"The total observed xs: "<<xscon<<" nb"<<std::endl;
3562 oa<<"=== mode =========== ratio ==========="<<std::endl;
3563 for(int i=0;i<vmode.size();i++){
3564 //std::cout<<"i: "<<xscon<<" "<<VXS[Vmode[i]][599]<<std::endl;
3565 int themode=getModeIndex(vmode[i]);
3566 if(vmode[i]==74110) continue;
3567 xscon -= VXS[themode ][599];
3568 }
3569 if(xscon<0) xscon=0;
3570 //debugging
3571 for(int i=0;i<vmode.size();i++){
3572 int themode=getModeIndex(vmode[i]);
3573 if(vmode[i]==74110) continue;
3574 oa<<vmode[i]<<"-th mode: "<<100*VXS[themode][599]/xscon0<<" % "<<std::endl;
3575 }
3576 oa<<"74110-th mode: "<<100*xscon/xscon0<<" % "<<std::endl;
3577
3578
3579}
3580
3581int EvtConExc::getModeIndex(int m){
3582 for (int i=0;i<vmode.size();i++){
3583 if(m==vmode[i]) return i;
3584 }
3585 std::cout<<"EvtConExc::getModeIndex: no index in vmode found "<<std::endl;
3586 abort();
3587}
3588
3589std::string EvtConExc::commandName(){
3590
3591 return std::string("ConExcPar");
3592
3593}
3594
3595void EvtConExc::command(std::string cmd){
3596 if (ncommand==lcommand){
3597 lcommand=10+2*lcommand;
3598 std::string* newcommands=new std::string[lcommand];
3599 int i;
3600 for(i=0;i<ncommand;i++){
3601 newcommands[i]=commands[i];
3602 }
3603 delete [] commands;
3604 commands=newcommands;
3605 }
3606 commands[ncommand]=cmd;
3607 ncommand++;
3608}
3609
3610
3611void EvtConExc::store(EvtDecayBase* jsdecay){
3612
3613 if (nconexcdecays==ntable){
3614
3615 EvtDecayBasePtr* newconexcdecays=new EvtDecayBasePtr[2*ntable+10];
3616 int i;
3617 for(i=0;i<ntable;i++){
3618 newconexcdecays[i]=conexcdecays[i];
3619 }
3620 ntable=2*ntable+10;
3621 delete [] conexcdecays;
3622 conexcdecays=newconexcdecays;
3623 }
3624
3625 conexcdecays[nconexcdecays++]=jsdecay;
3626
3627}
3628
3629
3630std::vector<std::string> EvtConExc::split(std::string str,std::string pattern)
3631{
3632 std::string::size_type pos;
3633 std::vector<std::string> result;
3634 str+=pattern;
3635 int size=str.size();
3636
3637 for(int i=0; i<size; i++)
3638 {
3639 pos=str.find(pattern,i);
3640 if(pos<size)
3641 {
3642 std::string s=str.substr(i,pos-i);
3643 result.push_back(s);
3644 i=pos+pattern.size()-1;
3645 }
3646 }
3647 return result;
3648}
3649
3650
3651void EvtConExc::InitPars(){
3652 threshold=0;
3653 beamEnergySpread=0;
3654 std::string pattern="=";
3655 for(int i=0;i<ncommand;i++){
3656 std::vector<std::string> result=split(commands[i],pattern);
3657 if(result[0]=="threshold") { threshold=atof(result[1].c_str());}else
3658 if(result[0]=="beamEnergySpread"){ beamEnergySpread=atof(result[1].c_str());}
3659 else{
3660 std::cout<<"Your parameter in decay card \""<<result[0]<<"\" incorect"<<std::endl;
3661 std::cout<<"Possible words: threshold, beamEnergySpread"<<std::endl;
3662 abort();
3663 }
3664 }
3665
3666}
3667
3668
3669double EvtConExc::energySpread(double mu,double sigma){
3670 //mu: mean value in Gaussian
3671 //sigma: variance in Gaussian
3672 rloop:
3673 int n=12;
3674 double ri=0;
3675 for(int i=0;i<n;i++){
3676 double pm= EvtRandom::Flat(0.,1);
3677 ri += pm;
3678 }
3679 double eta=sqrt(n*12.0)*(ri/12-0.5);
3680 double xsig= sigma*eta+mu;//smapling Gaussion value
3681 double mx0=staxsection->getXlw();
3682 double mx1=staxsection->getXup();
3683 if(xsig<mx0 || xsig>mx1) goto rloop;
3684 return xsig;
3685}
3686
3687void EvtConExc::calAF(double myecms){
3688
3689
3690 double _cms=myecms;
3691 double Esig = 0.0001; //sigularity engy
3692 double x = 2*Egamcut/_cms;
3693 double s = _cms*_cms;
3694
3695 double mhdL=staxsection->getXlw();
3696 double EgamH = (s-mhdL*mhdL)/2/sqrt(s);
3697
3698 int nint=50;
3699 int nmc= 1000;
3700 _xs0 = trapezoid(s,Esig,Egamcut,nint); // std::cout<<"Int: _xs0= "<<_xs0<<std::endl;
3701 //--- sigularity point x from 0 to 0.0001GeV
3702 double xb= 2*Esig/_cms;
3703 double sig_xs = SoftPhoton_xs(s,xb)*(staxsection->getXsection(_cms));
3704 _xs0 += sig_xs;
3705
3706 int Nstp=598;
3707 double stp=(EgamH - Egamcut)/Nstp;
3708 for(int i=0;i<Nstp;i++){//points within Egam(max) ~ Egamcut=25MeV
3709 double Eg0=EgamH - i*stp;
3710 double Eg1=EgamH - (i+1)*stp;
3711 double xsi= trapezoid(s,Eg1,Eg0,nint); // std::cout<<"Int xsi= " <<xsi<<std::endl;
3712
3713 AA[i]=(Eg1+Eg0)/2;
3714 double mhi=sqrt(_cms*_cms-2*_cms*AA[i]);
3715 double mh2=sqrt(_cms*_cms-2*_cms*Eg1);
3716 double binwidth = mh2-mhi;
3717 if(i==0) AF[0]=xsi;
3718 if(i>=1) AF[i]=AF[i-1]+xsi;
3719 RadXS[i]=xsi/stp;
3720 }
3721 //add the interval 0~Esig
3722 AA[598]=Egamcut; AA[599]=0; //AA[i]: E_gamma
3723 AF[598]=AF[597];
3724 AF[599]=AF[598]+ _xs0;
3725 //--
3726 for(int i=0;i<600;i++){
3727 MH[i]=sqrt(_cms*_cms-2*_cms*AA[i]);
3728 }
3729 //--debugging
3730 double bornXS = staxsection->getXsection(_cms);
3731 if(bornXS==0){std::cout<<"EvtConExc::calAF: bad BornXS at "<<_cms<<" GeV"<<std::endl;abort();}
3732 double fisr=AF[599]/bornXS;
3733 myFisr.push_back(fisr);
3734 //std::cout<<"fisr= "<<fisr<<std::endl;
3735}
3736
3737
3738void EvtConExc::OutStaISR(){
3739 int ntot=myFisr.size();
3740 double mymu=0;
3741 for(int i=0;i<ntot;i++){mymu += myFisr[i];}
3742 mymu /= ntot;
3743 double mysig=0;
3744 for(int i=0;i<ntot;i++){ mysig += (myFisr[i]-mymu)*(myFisr[i]-mymu);}
3745 mysig /=ntot;
3746 mysig = sqrt(mysig);
3747 std::cout<<"========= Iteration times over ISR factor: "<<ntot<<std::endl;
3748 std::cout<<" ISR factor * VP factor= observedXS/BornXS: "<<mymu<<" + "<<mysig<<std::endl;
3749}
3750
3751
3752double EvtConExc::trapezoid(double s, double El,double Eh,int n)// integration with trapezoid method
3753{
3754 double xL=2*El/sqrt(s);
3755 double xH=2*Eh/sqrt(s);
3756 double sum = 0.0;
3757 double gaps = (xH-xL)/double(n); //interval
3758 for (int i = 0; i < n; i++)
3759 {
3760 sum += (gaps/2.0) * (Rad2difXs(s,xL + i*gaps, staxsection) + Rad2difXs(s,xL + (i+1)*gaps,staxsection) );
3761 }
3762 return sum;
3763}
3764
3765double EvtConExc::trapezoid(double s, double El,double Eh,int n, EvtXsection* myxs)// integration with trapezoid method
3766{
3767 double xL=2*El/sqrt(s);
3768 double xH=2*Eh/sqrt(s);
3769 double sum = 0.0;
3770 double gaps = (xH-xL)/double(n); //interval
3771 for (int i = 0; i < n; i++)
3772 {
3773 sum += (gaps/2.0) * (Rad2difXs(s,xL + i*gaps, myxs) + Rad2difXs(s,xL + (i+1)*gaps,myxs) );
3774 }
3775 return sum;
3776}
mass
double mass
Definition: CosmicGenerator.cxx:138
n
const Int_t n
Definition: DataBase/tau_mode.c:65
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
Rad2difXs2
double Rad2difXs2(double *mhs)
Definition: EvtConExc.cc:2758
hadr5n12_
void hadr5n12_(float *, float *, float *, float *, float *, float *)
divdif_
double divdif_(float *, float *, int *, float *, int *)
oa
std::ofstream oa
Definition: EvtConExc.cc:199
Rad2difXs
double Rad2difXs(double *m)
Definition: EvtConExc.cc:2731
Rad2difXs_er2
double Rad2difXs_er2(double *mhs)
Definition: EvtConExc.cc:2773
dgauss_
double dgauss_(double(*fun)(double *), double *, double *, double *)
Rad2difXs_er
double Rad2difXs_er(double *m)
Definition: EvtConExc.cc:2745
polint_
void polint_(float *, float *, int *, float *, float *)
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition: EvtDiracSpinor.cc:75
xlow
DOUBLE_PRECISION xlow[20]
Definition: EvtPhokharaDef.hh:56
mstp
int mstp[200]
Definition: EvtPycont.cc:54
iter
EvtStreamInputIterator< typename Generator::result_type > iter(Generator gen, int N=0)
Definition: EvtStreamInputIterator.hh:100
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
sin
double sin(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: InstallArea/include/MdcGeom/MdcGeom/BesAngle.h:213
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
phiee
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rh2ee DOUBLE PRECISION omepi DOUBLE PRECISION om3ee DOUBLE PRECISION phiee
Definition: RRes.h:37
rhoee
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations xmb DOUBLE PRECISION xcb DOUBLE PRECISION xmDcs DOUBLE PRECISION xmBbc DOUBLE PRECISION rhoee
Definition: RRes.h:34
EvtConExc::checkdecay
bool checkdecay(EvtParticle *p)
Definition: EvtConExc.cc:3165
EvtConExc::Li2
double Li2(double x)
Definition: EvtConExc.cc:2812
EvtConExc::findMaxXS
void findMaxXS(EvtParticle *p)
Definition: EvtConExc.cc:2431
EvtConExc::getMode
static int getMode
Definition: EvtConExc.hh:194
EvtConExc::conexcmode
static int conexcmode
Definition: EvtConExc.hh:211
EvtConExc::addNarrowRXS
double addNarrowRXS(double mhi, double binwidth)
Definition: EvtConExc.cc:3326
EvtConExc::init
void init()
Definition: EvtConExc.cc:254
EvtConExc::narrowRXS
double narrowRXS(double mxL, double mxH)
Definition: EvtConExc.cc:3292
EvtConExc::VP_sampling
bool VP_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:2557
EvtConExc::command
void command(std::string cmd)
Definition: EvtConExc.cc:3595
EvtConExc::init_Br_ee
void init_Br_ee()
Definition: EvtConExc.cc:2686
EvtConExc::gamHXSection_er
double gamHXSection_er(double El, double Eh)
Definition: EvtConExc.cc:2413
EvtConExc::islgr
bool islgr(double *x, double *y, int n, double t)
Definition: EvtConExc.cc:2834
EvtConExc::showResMass
void showResMass()
Definition: EvtConExc.cc:3154
EvtConExc::lgr
double lgr(double *x, double *y, int n, double t)
Definition: EvtConExc.cc:2819
EvtConExc::ReadVP
void ReadVP()
Definition: EvtConExc.cc:3381
EvtConExc::baryonAng
double baryonAng(double mx)
Definition: EvtConExc.cc:3251
EvtConExc::InitPars
void InitPars()
Definition: EvtConExc.cc:3651
EvtConExc::Ros_xs
double Ros_xs(double mx, double bree, EvtId pid)
Definition: EvtConExc.cc:2718
EvtConExc::OutStaISR
void OutStaISR()
Definition: EvtConExc.cc:3738
EvtConExc::Rad1
double Rad1(double s, double x)
Definition: EvtConExc.cc:2568
EvtConExc::staxsection
static EvtXsection * staxsection
Definition: EvtConExc.hh:189
EvtConExc::myxsection
static EvtXsection * myxsection
Definition: EvtConExc.hh:189
EvtConExc::SetP4Rvalue
void SetP4Rvalue(EvtParticle *part, double mhdr, double xeng, double theta)
Definition: EvtConExc.cc:2959
EvtConExc::getModeIndex
int getModeIndex(int m)
Definition: EvtConExc.cc:3581
EvtConExc::meson_sampling
bool meson_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:2546
EvtConExc::multimode
static int multimode
Definition: EvtConExc.hh:211
EvtConExc::photonSampling
bool photonSampling(EvtParticle *part)
Definition: EvtConExc.cc:3261
EvtConExc::ISR_ang_integrate
double ISR_ang_integrate(double x, double theta)
Definition: EvtConExc.cc:2893
EvtConExc::SetMthr
static double SetMthr
Definition: EvtConExc.hh:192
EvtConExc::xs_sampling
bool xs_sampling(double xs)
Definition: EvtConExc.cc:2517
EvtConExc::selectMode
int selectMode(std::vector< int > vmod, double mhds)
Definition: EvtConExc.cc:3003
EvtConExc::gamHXSection
double gamHXSection(EvtParticle *p, double El, double Eh, int nmc=100000)
Definition: EvtConExc.cc:2286
EvtConExc::Rad1difXs
double Rad1difXs(EvtParticle *p)
Definition: EvtConExc.cc:2658
EvtConExc::sumExc
double sumExc(double mx)
Definition: EvtConExc.cc:3179
EvtConExc::init_mode
void init_mode(int mode)
Definition: EvtConExc.cc:741
EvtConExc::baryon_sampling
bool baryon_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:2529
EvtConExc::get_mode_index
int get_mode_index(int mode)
Definition: EvtConExc.cc:3440
EvtConExc::LLr
double LLr(double *x, double *y, int n, double t)
Definition: EvtConExc.cc:2849
EvtConExc::Rad2difXs
double Rad2difXs(EvtParticle *p)
Definition: EvtConExc.cc:2604
EvtConExc::SoftPhoton_xs
double SoftPhoton_xs(double s, double b)
Definition: EvtConExc.cc:2787
EvtConExc::angularSampling
bool angularSampling(EvtParticle *part)
Definition: EvtConExc.cc:3215
EvtConExc::difgamXs
double difgamXs(EvtParticle *p)
Definition: EvtConExc.cc:2476
EvtConExc::ISR_ang_sampling
double ISR_ang_sampling(double x)
Definition: EvtConExc.cc:2911
EvtConExc::checkEvtRatio
void checkEvtRatio()
Definition: EvtConExc.cc:3553
EvtConExc::initProbMax
void initProbMax()
Definition: EvtConExc.cc:1757
EvtConExc::Mhad_sampling
double Mhad_sampling(double *x, double *y)
Definition: EvtConExc.cc:2865
EvtConExc::_cms
static double _cms
Definition: EvtConExc.hh:190
EvtConExc::writeDecayTabel
void writeDecayTabel()
Definition: EvtConExc.cc:3468
EvtConExc::energySpread
double energySpread(double mu, double sigma)
Definition: EvtConExc.cc:3669
EvtConExc::get_mode
std::vector< EvtId > get_mode(int mode)
Definition: EvtConExc.cc:1246
EvtConExc::mlow
static double mlow
Definition: EvtConExc.hh:234
EvtConExc::resetResMass
void resetResMass()
Definition: EvtConExc.cc:3083
EvtConExc::getObsXsection
double getObsXsection(double mhds, int mode)
Definition: EvtConExc.cc:3449
EvtConExc::mk_VXS
void mk_VXS(double Esig, double Egamcut, double EgamH, int midx)
Definition: EvtConExc.cc:3401
EvtConExc::calAF
void calAF(double myecms)
Definition: EvtConExc.cc:3687
EvtConExc::hadron_angle_sampling
bool hadron_angle_sampling(EvtVector4R ppi, EvtVector4R pcm)
Definition: EvtConExc.cc:2258
EvtConExc::~EvtConExc
virtual ~EvtConExc()
Definition: EvtConExc.cc:205
EvtConExc::getResMass
void getResMass()
Definition: EvtConExc.cc:3118
EvtConExc::gam_sampling
bool gam_sampling(EvtParticle *p)
Definition: EvtConExc.cc:2503
EvtConExc::Rad2
double Rad2(double s, double x)
Definition: EvtConExc.cc:2580
EvtConExc::selectMass
double selectMass()
Definition: EvtConExc.cc:3335
EvtConExc::getVP
double getVP(double cms)
Definition: EvtConExc.cc:3360
EvtConExc::trapezoid
double trapezoid(double s, double a, double b, int n)
Definition: EvtConExc.cc:3752
EvtConExc::getName
void getName(std::string &name)
Definition: EvtConExc.cc:241
EvtConExc::decay
void decay(EvtParticle *p)
Definition: EvtConExc.cc:1763
EvtConExc::split
std::vector< std::string > split(std::string str, std::string pattern)
Definition: EvtConExc.cc:3630
EvtConExc::commandName
std::string commandName()
Definition: EvtConExc.cc:3589
EvtConExc::XS_max
static double XS_max
Definition: EvtConExc.hh:191
EvtConExc::Egam2Mhds
double Egam2Mhds(double Egam)
Definition: EvtConExc.cc:3462
EvtConExc::mup
static double mup
Definition: EvtConExc.hh:234
EvtConExc::clone
EvtDecayBase * clone()
Definition: EvtConExc.cc:247
EvtConExc::SetP4
void SetP4(EvtParticle *part, double mhdr, double xeng, double theta)
Definition: EvtConExc.cc:2935
EvtDecayBase::getArg
double getArg(int j)
Definition: EvtDecayBase.cc:564
EvtDecayBase::noProbMax
void noProbMax()
Definition: EvtDecayBase.cc:304
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:65
EvtDecayBase::getDaug
EvtId getDaug(int i)
Definition: EvtDecayBase.hh:66
EvtDecayTable::inChannelList
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
Definition: EvtDecayTable.cc:795
EvtGlobalSet::ConExcPythia
static bool ConExcPythia
Definition: EvtGlobalSet.hh:20
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::getStdHep
static int getStdHep(EvtId id)
Definition: EvtPDL.hh:56
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition: EvtPDL.hh:45
EvtPDL::reSetWidth
static void reSetWidth(EvtId i, double width)
Definition: EvtPDL.hh:69
EvtPDL::reSetMass
static void reSetMass(EvtId i, double mass)
Definition: EvtPDL.hh:68
EvtPDL::evtIdFromStdHep
static EvtId evtIdFromStdHep(int stdhep)
Definition: EvtPDL.cc:244
EvtPDL::name
static std::string name(EvtId i)
Definition: EvtPDL.hh:64
EvtPDL::getMinMass
static double getMinMass(EvtId i)
Definition: EvtPDL.hh:51
EvtPDL::getSpinType
static EvtSpinType::spintype getSpinType(EvtId i)
Definition: EvtPDL.hh:61
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:372
EvtParticle::makeDaughters
void makeDaughters(int ndaug, EvtId *id)
Definition: EvtParticle.cc:1177
EvtParticle::init
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtParticle::getP4Lab
EvtVector4R getP4Lab()
Definition: EvtParticle.cc:685
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cc:113
EvtParticle::setIntFlag
void setIntFlag(std::vector< int > vi)
Definition: EvtParticle.hh:151
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:121
EvtParticle::printTree
void printTree() const
Definition: EvtParticle.cc:897
EvtParticle::getNDaug
int getNDaug() const
Definition: EvtParticle.cc:125
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:85
EvtParticle::deleteDaughters
void deleteDaughters(bool keepChannel=false)
Definition: EvtParticle.cc:540
EvtParticle::mass
double mass() const
Definition: EvtParticle.cc:127
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cc:1071
EvtParticle::deleteTree
void deleteTree()
Definition: EvtParticle.cc:557
EvtRandom::Flat
static double Flat()
Definition: EvtRandom.cc:74
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
EvtVector4R::theta
double theta()
Definition: EvtVector4R.cc:249
EvtXsection::getXlw
double getXlw()
Definition: EvtXsection.hh:181
EvtXsection::getErr
double getErr(double mx)
Definition: EvtXsection.cc:260
EvtXsection::getMsg
std::string getMsg()
Definition: EvtXsection.hh:182
EvtXsection::setModes
void setModes(std::vector< int > vmd)
Definition: EvtXsection.cc:378
EvtXsection::getEr
std::vector< double > getEr()
Definition: EvtXsection.hh:179
EvtXsection::ini_data_multimode
void ini_data_multimode()
Definition: EvtXsection.cc:164
EvtXsection::getXup
double getXup()
Definition: EvtXsection.hh:180
EvtXsection::getUnit
std::string getUnit()
Definition: EvtXsection.hh:175
EvtXsection::getXsection
double getXsection(double mx)
Definition: EvtXsection.cc:254
EvtXsection::getYY
std::vector< double > getYY()
Definition: EvtXsection.hh:178
EvtXsection::setBW
void setBW(int pdg)
Definition: EvtXsection.cc:358
EvtXsection::getXX
std::vector< double > getXX()
Definition: EvtXsection.hh:177
mp
const double mp
Definition: incllambda.cxx:45
t
int t()
Definition: t.c:1