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