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