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