BOSS 7.0.1
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;
164
165EvtConExc::~EvtConExc() {
166 if(mydbg){
167 // xs->Write();
168 myfile->Write();
169 }
170 delete myxsection;
171 gamH->deleteTree();
172}
173
174void EvtConExc::getName(std::string& model_name){
175
176 model_name="ConExc";
177
178}
179
180EvtDecayBase* EvtConExc::clone(){
181
182 return new EvtConExc;
183
184}
185
186
187void EvtConExc::init(){
188 ReadVP();
189 //----------------
190 // check that there are 0 arguments
191 // checkNArg(1);
192 //Vector ISR process
193 VISR=0;
194 if(getNDaug()==2){
195 if(getDaugs()[0]==EvtPDL::getId("gamma") &&getDaugs()[1]==EvtPDL::getId("gamma*") || getDaugs()[0]==EvtPDL::getId("gamma*") &&getDaugs()[1]==EvtPDL::getId("gamma") ) VISR=1;
196 }else if(getNDaug()>2){std::cout<<"Bad daughter specified"<<std::endl;abort();}
197
198 cmsspread=0.0015; //CMC energy spread
199 testflag=0;
200 getResMass();
201 if(getArg(0) == -1){
202 radflag=0;mydbg=false;
203 _mode = getArg(0);
204 pdgcode = getArg(1);
205 pid = EvtPDL::evtIdFromStdHep(pdgcode );
206 nson = getNArg()-2;
207 std::cout<<"The decay daughters are "<<std::endl;
208 for(int i=2;i<getNArg();i++ ){son[i-2]=EvtPDL::evtIdFromStdHep(getArg(i));std::cout<<EvtPDL::name(son[i-2])<<" ";}
209 std::cout<<std::endl;
210 }else if(getArg(0)==-2){
211 radflag=0;mydbg=false;
212 _mode = getArg(0);
213 nson = getNArg()-1;
214 for(int i=1;i<getNArg();i++ ){son[i-1]=EvtPDL::evtIdFromStdHep(getArg(i));}
215 }
216 else if(getNArg()==1){ _mode = getArg(0);radflag=0;mydbg=false;}
217 else if(getNArg()==2){ _mode = getArg(0); radflag=getArg(1);mydbg=false;}
218 else if(getNArg()==3){ _mode = getArg(0); radflag=getArg(1);mydbg=getArg(2);}
219 else{std::cout<<"ConExc:umber of parameter should be 1,2 or 3, but you set "<<getNArg()<<std::endl;::abort(); }
220 //--- debugging
221 //std::cout<<_mode<<" "<<pdgcode<<" "<<nson<<" "<<EvtPDL::name(son[0])<<" "<<EvtPDL::name(son[1])<<std::endl;
222
223 gamId = EvtPDL::getId(std::string("gamma"));
224 init_mode(_mode);
225 XS_max=-1;
226 //-----debugging to make root file
227 if(mydbg){
228 myfile = new TFile("mydbg.root","RECREATE");
229 xs=new TTree ("xs","momentum of rad. gamma and hadrons");
230 xs->Branch("imode",&imode,"imode/I");
231 xs->Branch("pgam",pgam,"pgam[4]/D");
232 xs->Branch("phds",phds,"phds[4]/D");
233 xs->Branch("ph1",ph1,"ph1[4]/D");
234 xs->Branch("ph2",ph2,"ph2[4]/D");
235 xs->Branch("mhds",&mhds,"mhds/D");
236 xs->Branch("mass1",&mass1,"mass1/D");
237 xs->Branch("mass2",&mass2,"mass2/D");
238 xs->Branch("costheta",&costheta,"costheta/D");
239 xs->Branch("cosp",&cosp,"cosp/D");
240 xs->Branch("cosk",&cosk,"cosk/D");
241 xs->Branch("sumxs",&sumxs,"sumxs/D");
242 xs->Branch("selectmode",&selectmode,"selectmode/D");
243 }
244 //--- prepare the output information
245 EvtId parentId =EvtPDL::getId(std::string("vpho"));
246 EvtPDL::reSetWidth(parentId,0.0);
247 double parentMass = EvtPDL::getMass(parentId);
248 std::cout<<"parent mass = "<<parentMass<<std::endl;
249
250
251 EvtVector4R p_init(parentMass,0.0,0.0,0.0);
252 EvtParticle *p=EvtParticleFactory::particleFactory(parentId,p_init);
253 myxsection =new EvtXsection (_mode);
254 double mth=0;
255 double mx = EvtPDL::getMeanMass(parentId); //p->getP4().mass();
256 if(_mode ==-1){
257 myxsection->setBW(pdgcode);
258 for(int i=0;i<nson;i++){mth +=EvtPDL::getMeanMass(son[i]);}
259 if(mx<mth){std::cout<<"The CMS energy is lower than the threshold of the final states"<<std::endl;abort();}
260 }else if(_mode == -2){
261 mth=myxsection->getXlw();
262 }else{ mth=myxsection->getXlw();}
263 _cms = mx;
264 _unit = myxsection->getUnit();
265
266 std::cout<<"The specified mode is "<<_mode<<std::endl;
267 EvtConExc::getMode = _mode;
268 //std::cout<<"getXlw= "<<mth<<std::endl;
269
270 Egamcut = 0.001; //set photon energy cut according to the BESIII detector
271 double Esig = 0.0001; //sigularity engy
272 double x = 2*Egamcut/_cms;
273 double s = _cms*_cms;
274 double mhscut = sqrt(s*(1-x));
275
276 //for vaccum polarization
277 float fe,fst2,fder,ferrder,fdeg,ferrdeg;
278 fe=_cms;
279 //using the vacuum pol. value as given by http://www-com.physik.hu-berlin.de/~fjeger/software.html
280 vph=getVP(_cms);
281 if(3.0943<_cms && _cms<3.102) vph=1;// For J/psi, the vacuum factor is included in the resonance
282 std::cout<<"sum of leptonic, hadronic contributions(u,d,s,c,b) to vacuum polarization: "<<vph<<" @"<<fe<<"GeV"<<std::endl;
283
284 //caculate the xs for ISR to psiprime, J/psi and phi narrow resonance.
285 double totxsIRS=0;
286 init_Br_ee();
287 double mthrld = EvtPDL::getMeanMass(daugs[0]); //threshold mass of hadrons
288 for(int jj=1;jj<_ndaugs;jj++){
289 mthrld += EvtPDL::getMeanMass(daugs[jj]);
290 }
291
292 ISRXS.clear();ISRM.clear();ISRFLAG.clear();
293 for(int ii=0;ii<3;ii++){
294 double mR = EvtPDL::getMeanMass(ResId[ii]);
295 if(mx<mR || _mode !=74110) continue;
296 double narRxs=Ros_xs(mx,BR_ee[ii],ResId[ii]);
297 std::cout<<"The corss section for gamma "<<EvtPDL::name(ResId[ii]).c_str()<<" is: "<< narRxs<<"*Br (nb)."<<std::endl;
298 ISRXS.push_back(narRxs);
299 ISRM.push_back(mR);
300 ISRFLAG.push_back(1.);
301 ISRID.push_back(ResId[ii]);
302 totxsIRS += narRxs;
303 }
304 std::cout<<"==========================================================================="<<std::endl;
305
306 //-- calculated with MC integration method
307 double mhdL=myxsection->getXlw();
308 EgamH = (s-mhdL*mhdL)/2/sqrt(s);
309 int nmc=1000000;
310 _xs0 = gamHXSection(s,Esig,Egamcut,nmc);
311 _er0 = _er1; // stored when calculate _xs0
312 std::cout<<"_er0= "<<_er0<<std::endl;
313 _xs1 = gamHXSection(s,Egamcut,EgamH,nmc);
314 double xs1_err = _er1;
315
316
317 //--- sigularity point x from 0 to 0.0001GeV
318 double xb= 2*Esig/_cms;
319 double sig_xs = SoftPhoton_xs(s,xb)*(myxsection->getXsection(mx));
320 _xs0 += sig_xs;
321
322 //prepare the observed cross section with interpotation function divdif in CERNLIB
323 testflag=1;
324 int Nstp=598;
325 double stp=(EgamH - Egamcut)/Nstp;
326 for(int i=0;i<Nstp;i++){//points within Egam(max) ~ Egamcut=25MeV
327 double Eg0=EgamH - i*stp;
328 double Eg1=EgamH - (i+1)*stp;
329 int nmc=20000;
330 double xsi= gamHXSection(s,Eg1,Eg0,nmc);
331 AA[i]=(Eg1+Eg0)/2;
332 double mhi=sqrt(_cms*_cms-2*_cms*AA[i]);
333 double mh2=sqrt(_cms*_cms-2*_cms*Eg1);
334 double binwidth = mh2-mhi;
335 if(_mode==74110) xsi += addNarrowRXS(mhi,binwidth);
336 if(i==0) AF[0]=xsi;
337 if(i>=1) AF[i]=AF[i-1]+xsi;
338 RadXS[i]=xsi/stp;
339 }
340 //add the interval 0~Esig
341 AA[598]=Egamcut; AA[599]=0; //AA[i]: E_gamma
342 AF[598]=AF[597];
343 AF[599]=AF[598]+ _xs0;
344 RadXS[599]=_xs0;
345 std::cout<<"mhadL= "<<mhdL<<" ecms "<<_cms<<std::endl;
346 for(int i=0;i<600;i++){
347 MH[i]=sqrt(_cms*_cms-2*_cms*AA[i]);
348 //std::cout<<i<<" Egamma "<<AA[i]<<" Mhadons "<<MH[i]<<std::endl;
349 }
350 //for debugging
351 //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;}
352 std::cout<<"EgamH="<<EgamH<<" "<<_xs0+_xs1<<" AF[599]="<<AF[599]<<" AF[598] "<<AF[598]<<std::endl;
353 //double Etst=0.0241838;
354 //std::cout<<"Etst="<<Etst<<" "<<gamHXSection(s,Etst,EgamH,nmc)<<" "<< lgr(AA,AF,600,Etst)<<std::endl; abort();
355
356 //for information output
357 double bornXS = myxsection->getXsection(mx); // std::cout<<"BornXS= "<<bornXS<<std::endl;
358 double bornXS_er=myxsection->getErr(mx);
359 double obsvXS = _xs0 + _xs1; //gamHXSection(mth ,mx);
360 double obsvXS_er= _er0 + xs1_err;
361 double corr = obsvXS/bornXS;
362 double corr_er =corr*sqrt(bornXS_er*bornXS_er/bornXS/bornXS + obsvXS_er*obsvXS_er/obsvXS/obsvXS);
363
364
365 if(bornXS==0){std::cout<<"EvtConExc::init : The Born cross section at this point is zero!"<<std::endl;abort();}
366 if(obsvXS==0){std::cout<<"EvtConExc::init : The observed cross section at this point is zero!"<<std::endl;abort();}
367
368 //_xs0 += bornXS;// events with very soft photon (Egam<0.025) are taken as the born process
369 //_er0 = sqrt(_er0*_er0 + bornXS_er*bornXS_er);
370
371 std::cout<<""<<std::endl;
372 std::cout<<"========= Generation using cross section (Egamcut="<<Egamcut<<" GeV) =============="<<std::endl;
373 std::cout<<" sqrt(s)= "<<mx<< " GeV "<<std::endl;
374 if(_mode>=0 || _mode ==-2 ){
375 std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<"+/-"<<_er0<<" in Unit "<<_unit.c_str()<<std::endl;
376 std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"+/-"<<xs1_err<<" in Unit "<<_unit.c_str()<<std::endl;
377 std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
378 std::cout<<"The radiative correction factor f_ISR*f_vacuum= sigma_obs/sigma_born(s) = "<<corr<<"+/-"<< fabs(corr_er)<<","<<std::endl;
379 std::cout<<"which is calcualted with these quantities:"<<std::endl;
380 std::cout<<"the observed cross section is "<<obsvXS<<"+/-"<< obsvXS_er<<_unit.c_str()<<std::endl;
381 std::cout<<"and the Born cross section (s) is "<<bornXS<<"+/-"<<bornXS_er<<_unit.c_str()<<std::endl;
382 std::cout<<"and the vaccum polarziation factor (lep. + hadr.) 1/|1-Pi|^2= "<<vph<<std::endl;
383 std::cout<<"Within the range "<<myxsection->getXlw()<<" ~"<<myxsection->getXup()<<" GeV, "<<myxsection->getMsg().c_str()<<std::endl;
384 std::cout<<"==========================================================================="<<std::endl;
385 }else if(_mode==-1){
386 std::cout<<"The generated born cross section (sigma_b): "<<_xs0<<" *Br_ee"<<" in Unit "<<_unit.c_str()<<std::endl;
387 std::cout<<"The generated radiative correction cross section (sigma_isr): "<<_xs1<<"*Br_ee in Unit "<<_unit.c_str()<<std::endl;
388 std::cout<<"The Born process is sampled according to the ratio sigma_b/(sigma_b+sigma_isr)"<<std::endl;
389 std::cout<<"The radiative correction factor sigma_obs/sigma_born(s) = (1+delta) = "<<corr<<"+/-"<< fabs(corr_er)<<std::endl;
390 std::cout<<"==========================================================================="<<std::endl;
391 }
392 std::cout<<std::endl;
393 std::cout<<std::endl;
394
395 findMaxXS(p);
396 _mhdL=myxsection->getXlw();
397 //--debugging
398 //std::cout<<"maxXS= "<<maxXS<<std::endl;
399
400 if(_xs0 == 0 && _xs1==0){std::cout<<"EvtConExc::ini() has zero cross section"<<std::endl;abort();}
401
402 std::cout<<std::endl;
403 std::cout<<std::endl;
404
405 //--- for debugging
406 if(mydbg && _mode!=74110){
407 int nd = myxsection->getYY().size();
408 double xx[10000],yy[10000],xer[10000],yer[10000];
409 for(int i=0;i<nd;i++){
410 xx[i]=myxsection->getXX()[i];
411 yy[i]=myxsection->getYY()[i];
412 yer[i]=myxsection->getEr()[i];
413 xer[i]=0.;
414 //std::cout<<"yy "<<yy[i]<<std::endl;
415 }
416 myth=new TH1F("myth","Exp.data",200,xx[0],xx[nd]);
417 for(int i=0;i<nd;i++){
418 myth->Fill(xx[i],yy[i]);
419 }
420 myth->SetError(yer);
421 myth->Write();
422
423 //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
424 }else
425
426 if(mydbg && _mode==74110 ){
427 int nd = myxsection->getYY().size();
428 double xx[10000],yy[10000],xer[10000],yer[10000],ysum[10000],yysum[10000];
429 for(int i=0;i<nd;i++){
430 xx[i]=myxsection->getXX()[i];
431 yy[i]=myxsection->getYY()[i];
432 yer[i]=myxsection->getEr()[i];
433 xer[i]=0.;
434 //std::cout<<"yy "<<yy[i]<<std::endl;
435 }
436#include "sty.C"
437 double xhigh=5.0;
438 double xlow = myxsection->getXlw();
439 TGraphErrors *Gdt = new TGraphErrors(nd,xx,yy,xer,yer);
440
441 myth=new TH1F("myth","Exp.data",600,xlow,xhigh);
442 Xsum=new TH1F("sumxs","sum of exclusive xs",600,xlow,xhigh);
443 double mstp=(xhigh-xlow)/600;
444 for(int i=0;i<600;i++){
445 double mx=i*mstp+xlow;
446 double xsi = myxsection->getXsection(mx);
447 if(xsi==0) continue;
448 myth->Fill(mx,xsi);
449 //std::cout<<mx<<" "<<xsi<<std::endl;
450 }
451
452 for(int i=0;i<600;i++){
453 double mx=i*mstp+xlow;
454 ysum[i]=sumExc(mx);
455 if(ysum[i]==0) continue;
456 Xsum->Fill(mx,ysum[i]);
457 //std::cout<<mx<<" "<<ysum[i]<<std::endl;
458 }
459
460 for(int i=0;i<nd;i++){
461 yysum[i]=sumExc(xx[i]);
462 }
463
464 myth->SetError(yer);
465 myth->Write();
466 Xsum->Write();
467
468 TGraphErrors *Gsum = new TGraphErrors(nd,xx,yysum,xer,yer);
469 //draw graph
470 double ex[600]={600*0};
471 double ey[600],ta[600];
472 double exz[600]={600*0};
473 double eyz[600]={600*0};
474 for(int i=0;i<600;i++){
475 ey[i]=AF[i]*0.001;
476 }
477 TGraphErrors *gr0 = new TGraphErrors(600,MH,AF,ex,ey);
478 TGraphErrors *gr1 = new TGraphErrors(600,MH,RadXS,exz,eyz);
479 TPostScript *ps = new TPostScript("xsobs.ps",111);
480 TCanvas *c1 = new TCanvas("c1","TGraphs for data",200,10,600,400);
481 gPad-> SetLogy(1);
482 // c1->SetLogy(1);
483 gStyle->SetCanvasColor(0);
484 gStyle->SetStatBorderSize(1);
485 c1->Divide(1,1);
486
487 c1->Update();
488 ps->NewPage();
489 c1->Draw();
490 c1->cd(1);
491 c1->SetLogy(1);
492 gr0->SetMarkerStyle(10);
493 gr0->Draw("AP");
494 gr0->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
495 gr0->GetYaxis()->SetTitle("observed cross section (nb)");
496 gr0->Write();
497
498 c1->Update();
499 ps->NewPage();
500 c1->Draw();
501 c1->cd(1);
502 c1->SetLogy(1);
503 gr1->SetMarkerStyle(10);
504 gr1->Draw("AP");
505 gr1->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
506 gr1->GetYaxis()->SetTitle("Rad*BornXS");
507 gr1->Write();
508
509 //--------- imposing graphs to a pad
510 TMultiGraph *mg = new TMultiGraph();
511 mg->SetTitle("Exclusion graphs");
512 Gdt->SetMarkerStyle(2);
513 Gdt->SetMarkerSize(1);
514 Gsum->SetLineColor(2);
515 Gsum->SetLineWidth(1);
516 mg->Add(Gdt);
517 mg->Add(Gsum);
518
519 c1->Update();
520 ps->NewPage();
521 c1->Draw();
522 c1->cd(1);
523 gPad-> SetLogy(1);
524 mg->Draw("APL");
525 mg->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
526 mg->GetYaxis()->SetTitle("observed cross section (nb)");
527 mg->Write();
528 //-------
529
530 c1->Update();
531 ps->NewPage();
532 c1->Draw();
533 c1->cd(1);
534 Gdt->SetMarkerStyle(2);
535 Gdt->Draw("AP");
536 Gdt->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
537 Gdt->GetYaxis()->SetTitle("observed cross section (nb)");
538 Gdt->Write();
539
540 c1->Update();
541 ps->NewPage();
542 c1->Draw();
543 c1->cd(1);
544 Gsum->SetMarkerStyle(2);
545 Gsum->SetMarkerSize(1);
546 Gsum->Draw("AP");
547 Gsum->SetLineWidth(0);
548 Gsum->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
549 Gsum->GetYaxis()->SetTitle("observed cross section (nb)");
550 Gsum->Write();
551
552 c1->Update();
553 ps->NewPage();
554 c1->Draw();
555 c1->cd(1);
556 // gPad-> SetLogx(1);
557 Gdt->SetMarkerStyle(2);
558 Gdt->SetMarkerSize(1);
559 Gdt->SetMaximum(8000.0);
560 Gsum->SetLineColor(2);
561 Gsum->SetLineWidth(1.5);
562 Gsum->Draw("AL"); //A: draw axis
563 Gdt->Draw("P"); // superpose to the Gsum, without A-option
564 Gsum->GetXaxis()->SetTitle("M_{hadrons} (GeV)");
565 Gsum->GetYaxis()->SetTitle("cross section (nb)");
566 Gsum->Write();
567
568 ps->Close();
569 } //if(mydbg)_block
570 //-------------------------
571
572}
573
574
575//--
576void EvtConExc::init_mode(int mode){
577 if(mode ==0 ){
578 _ndaugs =2;
579 daugs[0]=EvtPDL::getId(std::string("p+"));
580 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
581 }else if(mode ==1 ){
582 _ndaugs =2;
583 daugs[0]=EvtPDL::getId(std::string("n0"));
584 daugs[1]=EvtPDL::getId(std::string("anti-n0"));
585 }else if(mode ==2 ){
586 _ndaugs =2;
587 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
588 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
589 }else if(mode ==3 ){
590 _ndaugs =2;
591 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
592 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
593 }else if(mode ==4 ){
594 _ndaugs =2;
595 daugs[0]=EvtPDL::getId(std::string("Lambda0"));
596 daugs[1]=EvtPDL::getId(std::string("anti-Sigma0"));
597 }else if(mode ==5 ){
598 _ndaugs =2;
599 daugs[0]=EvtPDL::getId(std::string("Sigma0"));
600 daugs[1]=EvtPDL::getId(std::string("anti-Lambda0"));
601 } else if(mode ==6 ){
602 _ndaugs =2;
603 daugs[0]=EvtPDL::getId(std::string("pi+"));
604 daugs[1]=EvtPDL::getId(std::string("pi-"));
605 }else if(mode ==7 ){
606 _ndaugs =3;
607 daugs[0]=EvtPDL::getId(std::string("pi+"));
608 daugs[1]=EvtPDL::getId(std::string("pi-"));
609 daugs[2]=EvtPDL::getId(std::string("pi0"));
610 }else if(mode ==8 ){
611 _ndaugs =3;
612 daugs[0]=EvtPDL::getId(std::string("K+"));
613 daugs[1]=EvtPDL::getId(std::string("K-"));
614 daugs[2]=EvtPDL::getId(std::string("pi0"));
615 }else if(mode ==9 ){
616 _ndaugs =3;
617 daugs[0]=EvtPDL::getId(std::string("K_S0"));
618 daugs[1]=EvtPDL::getId(std::string("K+"));
619 daugs[2]=EvtPDL::getId(std::string("pi-"));
620 }else if(mode ==10 ){
621 _ndaugs =3;
622 daugs[0]=EvtPDL::getId(std::string("K_S0"));
623 daugs[1]=EvtPDL::getId(std::string("K-"));
624 daugs[2]=EvtPDL::getId(std::string("pi+"));
625 }else if(mode ==11 ){
626 _ndaugs =3;
627 daugs[0]=EvtPDL::getId(std::string("K+"));
628 daugs[1]=EvtPDL::getId(std::string("K+"));
629 daugs[2]=EvtPDL::getId(std::string("eta"));
630 }else if(mode ==12 ){
631 _ndaugs =4;
632 daugs[0]=EvtPDL::getId(std::string("pi+"));
633 daugs[1]=EvtPDL::getId(std::string("pi-"));
634 daugs[2]=EvtPDL::getId(std::string("pi+"));
635 daugs[3]=EvtPDL::getId(std::string("pi-"));
636 }else if(mode ==13 ){
637 _ndaugs =4;
638 daugs[0]=EvtPDL::getId(std::string("pi+"));
639 daugs[1]=EvtPDL::getId(std::string("pi-"));
640 daugs[2]=EvtPDL::getId(std::string("pi0"));
641 daugs[3]=EvtPDL::getId(std::string("pi0"));
642 }else if(mode ==14 ){
643 _ndaugs =4;
644 daugs[0]=EvtPDL::getId(std::string("K+"));
645 daugs[1]=EvtPDL::getId(std::string("K-"));
646 daugs[2]=EvtPDL::getId(std::string("pi+"));
647 daugs[3]=EvtPDL::getId(std::string("pi-"));
648 }else if(mode ==15 ){
649 _ndaugs =4;
650 daugs[0]=EvtPDL::getId(std::string("K+"));
651 daugs[1]=EvtPDL::getId(std::string("K-"));
652 daugs[2]=EvtPDL::getId(std::string("pi0"));
653 daugs[3]=EvtPDL::getId(std::string("pi0"));
654 }else if(mode ==16 ){
655 _ndaugs =4;
656 daugs[0]=EvtPDL::getId(std::string("K+"));
657 daugs[1]=EvtPDL::getId(std::string("K-"));
658 daugs[2]=EvtPDL::getId(std::string("K+"));
659 daugs[3]=EvtPDL::getId(std::string("K-"));
660 }else if(mode ==17 ){
661 _ndaugs =5;
662 daugs[0]=EvtPDL::getId(std::string("pi+"));
663 daugs[1]=EvtPDL::getId(std::string("pi-"));
664 daugs[2]=EvtPDL::getId(std::string("pi+"));
665 daugs[3]=EvtPDL::getId(std::string("pi-"));
666 daugs[4]=EvtPDL::getId(std::string("pi0"));
667 }else if(mode ==18 ){
668 _ndaugs =5;
669 daugs[0]=EvtPDL::getId(std::string("pi+"));
670 daugs[1]=EvtPDL::getId(std::string("pi-"));
671 daugs[2]=EvtPDL::getId(std::string("pi+"));
672 daugs[3]=EvtPDL::getId(std::string("pi-"));
673 daugs[4]=EvtPDL::getId(std::string("eta"));
674 }else if(mode ==19 ){
675 _ndaugs =5;
676 daugs[0]=EvtPDL::getId(std::string("K+"));
677 daugs[1]=EvtPDL::getId(std::string("K-"));
678 daugs[2]=EvtPDL::getId(std::string("pi+"));
679 daugs[3]=EvtPDL::getId(std::string("pi-"));
680 daugs[4]=EvtPDL::getId(std::string("pi0"));
681 }else if(mode ==20 ){
682 _ndaugs =5;
683 daugs[0]=EvtPDL::getId(std::string("K+"));
684 daugs[1]=EvtPDL::getId(std::string("K-"));
685 daugs[2]=EvtPDL::getId(std::string("pi+"));
686 daugs[3]=EvtPDL::getId(std::string("pi-"));
687 daugs[4]=EvtPDL::getId(std::string("eta"));
688 }else if(mode ==21 ){
689 _ndaugs =6;
690 daugs[0]=EvtPDL::getId(std::string("pi+"));
691 daugs[1]=EvtPDL::getId(std::string("pi-"));
692 daugs[2]=EvtPDL::getId(std::string("pi+"));
693 daugs[3]=EvtPDL::getId(std::string("pi-"));
694 daugs[4]=EvtPDL::getId(std::string("pi+"));
695 daugs[5]=EvtPDL::getId(std::string("pi-"));
696 }else if(mode ==22 ){
697 _ndaugs =6;
698 daugs[0]=EvtPDL::getId(std::string("pi+"));
699 daugs[1]=EvtPDL::getId(std::string("pi-"));
700 daugs[2]=EvtPDL::getId(std::string("pi+"));
701 daugs[3]=EvtPDL::getId(std::string("pi-"));
702 daugs[4]=EvtPDL::getId(std::string("pi0"));
703 daugs[5]=EvtPDL::getId(std::string("pi0"));
704 }else if(mode == 23){
705 _ndaugs =2;
706 daugs[0]=EvtPDL::getId(std::string("eta"));
707 daugs[1]=EvtPDL::getId(std::string("phi"));
708 }else if(mode == 24){
709 _ndaugs =2;
710 daugs[0]=EvtPDL::getId(std::string("phi"));
711 daugs[1]=EvtPDL::getId(std::string("pi0"));
712 }else if(mode == 25){
713 _ndaugs =2;
714 daugs[0]=EvtPDL::getId(std::string("K+"));
715 daugs[1]=EvtPDL::getId(std::string("K*-"));
716 }else if(mode == 26){
717 _ndaugs =2;
718 daugs[0]=EvtPDL::getId(std::string("K-"));
719 daugs[1]=EvtPDL::getId(std::string("K*+"));
720 }else if(mode == 27){
721 _ndaugs =2;
722 daugs[0]=EvtPDL::getId(std::string("K_S0"));
723 daugs[1]=EvtPDL::getId(std::string("anti-K*0"));
724 }else if(mode == 28){
725 _ndaugs =3;
726 daugs[0]=EvtPDL::getId(std::string("K*0"));
727 daugs[1]=EvtPDL::getId(std::string("K+"));
728 daugs[2]=EvtPDL::getId(std::string("pi-"));
729 }else if(mode == 29){
730 _ndaugs =3;
731 daugs[0]=EvtPDL::getId(std::string("K*0"));
732 daugs[1]=EvtPDL::getId(std::string("K-"));
733 daugs[2]=EvtPDL::getId(std::string("pi+"));
734 }else if(mode == 30){
735 _ndaugs =3;
736 daugs[0]=EvtPDL::getId(std::string("K*+"));
737 daugs[1]=EvtPDL::getId(std::string("K-"));
738 daugs[2]=EvtPDL::getId(std::string("pi0"));
739 }else if(mode == 31){
740 _ndaugs =3;
741 daugs[0]=EvtPDL::getId(std::string("K*-"));
742 daugs[1]=EvtPDL::getId(std::string("K+"));
743 daugs[2]=EvtPDL::getId(std::string("pi0"));
744 }else if(mode == 32){
745 _ndaugs =3;
746 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
747 daugs[1]=EvtPDL::getId(std::string("K+"));
748 daugs[2]=EvtPDL::getId(std::string("pi-"));
749 }else if(mode == 33){
750 _ndaugs =3;
751 daugs[0]=EvtPDL::getId(std::string("K_2*0"));
752 daugs[1]=EvtPDL::getId(std::string("K-"));
753 daugs[2]=EvtPDL::getId(std::string("pi+"));
754 }else if(mode == 34){
755 _ndaugs =3;
756 daugs[0]=EvtPDL::getId(std::string("K+"));
757 daugs[1]=EvtPDL::getId(std::string("K-"));
758 daugs[2]=EvtPDL::getId(std::string("rho0"));
759 }else if(mode == 35){
760 _ndaugs =3;
761 daugs[0]=EvtPDL::getId(std::string("phi"));
762 daugs[1]=EvtPDL::getId(std::string("pi-"));
763 daugs[2]=EvtPDL::getId(std::string("pi+"));
764 }else if(mode == 36){
765 _ndaugs =2;
766 daugs[0]=EvtPDL::getId(std::string("phi"));
767 daugs[1]=EvtPDL::getId(std::string("f_0"));
768 }else if(mode == 37){
769 _ndaugs =3;
770 daugs[0]=EvtPDL::getId(std::string("eta"));
771 daugs[1]=EvtPDL::getId(std::string("pi+"));
772 daugs[2]=EvtPDL::getId(std::string("pi-"));
773 }else if(mode == 38){
774 _ndaugs =3;
775 daugs[0]=EvtPDL::getId(std::string("omega"));
776 daugs[1]=EvtPDL::getId(std::string("pi+"));
777 daugs[2]=EvtPDL::getId(std::string("pi-"));
778 }else if(mode == 39){
779 _ndaugs =2;
780 daugs[0]=EvtPDL::getId(std::string("omega"));
781 daugs[1]=EvtPDL::getId(std::string("f_0"));
782 }else if(mode == 40){
783 _ndaugs =3;
784 daugs[0]=EvtPDL::getId(std::string("eta'"));
785 daugs[1]=EvtPDL::getId(std::string("pi+"));
786 daugs[2]=EvtPDL::getId(std::string("pi-"));
787 }else if(mode == 41){
788 _ndaugs =3;
789 daugs[0]=EvtPDL::getId(std::string("f_1"));
790 daugs[1]=EvtPDL::getId(std::string("pi+"));
791 daugs[2]=EvtPDL::getId(std::string("pi-"));
792 }else if(mode == 42){
793 _ndaugs =3;
794 daugs[0]=EvtPDL::getId(std::string("omega"));
795 daugs[1]=EvtPDL::getId(std::string("K+"));
796 daugs[2]=EvtPDL::getId(std::string("K-"));
797 }else if(mode == 43){
798 _ndaugs =4;
799 daugs[0]=EvtPDL::getId(std::string("omega"));
800 daugs[1]=EvtPDL::getId(std::string("pi+"));
801 daugs[2]=EvtPDL::getId(std::string("pi-"));
802 daugs[3]=EvtPDL::getId(std::string("pi0"));
803 }else if(mode == 44){
804 _ndaugs =2;
805 daugs[0]=EvtPDL::getId(std::string("Sigma-"));
806 daugs[1]=EvtPDL::getId(std::string("anti-Sigma+"));
807 }else if(mode == 45){
808 _ndaugs =2;
809 daugs[0]=EvtPDL::getId(std::string("K+"));
810 daugs[1]=EvtPDL::getId(std::string("K-"));
811 }else if(mode == 46){
812 _ndaugs =2;
813 daugs[0]=EvtPDL::getId(std::string("K_S0"));
814 daugs[1]=EvtPDL::getId(std::string("K_L0"));
815 }else if(mode == 47){
816 _ndaugs =2;
817 daugs[0]=EvtPDL::getId(std::string("omega"));
818 daugs[1]=EvtPDL::getId(std::string("eta"));
819 }else if(mode == 48){
820 _ndaugs =2;
821 daugs[0]=EvtPDL::getId(std::string("phi"));
822 daugs[1]=EvtPDL::getId(std::string("eta'"));
823 }else if(mode == 49){
824 _ndaugs =3;
825 daugs[0]=EvtPDL::getId(std::string("phi"));
826 daugs[1]=EvtPDL::getId(std::string("K+"));
827 daugs[2]=EvtPDL::getId(std::string("K-"));
828 }else if(mode == 50){
829 _ndaugs =3;
830 daugs[0]=EvtPDL::getId(std::string("p+"));
831 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
832 daugs[2]=EvtPDL::getId(std::string("pi0"));
833 }else if(mode == 51){
834 _ndaugs =3;
835 daugs[0]=EvtPDL::getId(std::string("p+"));
836 daugs[1]=EvtPDL::getId(std::string("anti-p-"));
837 daugs[2]=EvtPDL::getId(std::string("eta"));
838 }else if(mode == 65){
839 _ndaugs =3;
840 daugs[0]=EvtPDL::getId(std::string("D-"));
841 daugs[1]=EvtPDL::getId(std::string("D*0"));
842 daugs[2]=EvtPDL::getId(std::string("pi+"));
843 }else if(mode == 66){
844 _ndaugs =3;
845 daugs[0]=EvtPDL::getId(std::string("D+"));
846 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
847 daugs[2]=EvtPDL::getId(std::string("pi-"));
848 }else if(mode == 67){
849 _ndaugs =2;
850 daugs[0]=EvtPDL::getId(std::string("D*0"));
851 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
852 }else if(mode == 68){
853 _ndaugs =2;
854 daugs[0]=EvtPDL::getId(std::string("D0"));
855 daugs[1]=EvtPDL::getId(std::string("anti-D*0"));
856
857 }else if(mode == 69){
858 _ndaugs =2;
859 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
860 daugs[1]=EvtPDL::getId(std::string("D*0"));
861
862 }else if(mode == 70){
863 _ndaugs =2;
864 daugs[0]=EvtPDL::getId(std::string("D0"));
865 daugs[1]=EvtPDL::getId(std::string("anti-D0"));
866
867}else if(mode == 71){
868 _ndaugs =2;
869 daugs[0]=EvtPDL::getId(std::string("D+"));
870 daugs[1]=EvtPDL::getId(std::string("D-"));
871 }else if(mode == 72){
872 _ndaugs =2;
873 daugs[0]=EvtPDL::getId(std::string("D+"));
874 daugs[1]=EvtPDL::getId(std::string("D*-"));
875
876 }else if(mode == 73){
877 _ndaugs =2;
878 daugs[0]=EvtPDL::getId(std::string("D-"));
879 daugs[1]=EvtPDL::getId(std::string("D*+"));
880
881 }else if(mode == 74){
882 _ndaugs =2;
883 daugs[0]=EvtPDL::getId(std::string("D*+"));
884 daugs[1]=EvtPDL::getId(std::string("D*-"));
885
886 }else if(mode == 75){
887 _ndaugs =3;
888 daugs[0]=EvtPDL::getId(std::string("D0"));
889 daugs[1]=EvtPDL::getId(std::string("D-"));
890 daugs[2]=EvtPDL::getId(std::string("pi+"));
891 }else if(mode == 76){
892 _ndaugs =3;
893 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
894 daugs[1]=EvtPDL::getId(std::string("D+"));
895 daugs[2]=EvtPDL::getId(std::string("pi-"));
896 }else if(mode == 77){
897 _ndaugs =3;
898 daugs[0]=EvtPDL::getId(std::string("D0"));
899 daugs[1]=EvtPDL::getId(std::string("D*-"));
900 daugs[2]=EvtPDL::getId(std::string("pi+"));
901 }else if(mode == 78){
902 _ndaugs =3;
903 daugs[0]=EvtPDL::getId(std::string("anti-D0"));
904 daugs[1]=EvtPDL::getId(std::string("D*+"));
905 daugs[2]=EvtPDL::getId(std::string("pi-"));
906 }else if(mode == 79){
907 _ndaugs =3;
908 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
909 daugs[1]=EvtPDL::getId(std::string("pi0"));
910 daugs[2]=EvtPDL::getId(std::string("pi0"));
911 }else if(mode == 80){
912 _ndaugs =2;
913 daugs[0]=EvtPDL::getId(std::string("eta"));
914 daugs[1]=EvtPDL::getId(std::string("J/psi"));
915 }else if(mode == 81){
916 _ndaugs =3;
917 daugs[0]=EvtPDL::getId(std::string("h_c"));
918 daugs[1]=EvtPDL::getId(std::string("pi+"));
919 daugs[2]=EvtPDL::getId(std::string("pi-"));
920 }else if(mode == 82){
921 _ndaugs =3;
922 daugs[0]=EvtPDL::getId(std::string("h_c"));
923 daugs[1]=EvtPDL::getId(std::string("pi0"));
924 daugs[2]=EvtPDL::getId(std::string("pi0"));
925 }else if(mode == 83){
926 _ndaugs =3;
927 daugs[0]=EvtPDL::getId(std::string("J/psi"));
928 daugs[1]=EvtPDL::getId(std::string("K+"));
929 daugs[2]=EvtPDL::getId(std::string("K-"));
930 }else if(mode == 84){
931 _ndaugs =3;
932 daugs[0]=EvtPDL::getId(std::string("J/psi"));
933 daugs[1]=EvtPDL::getId(std::string("K_S0"));
934 daugs[2]=EvtPDL::getId(std::string("K_S0"));
935 }else if(mode == 90){
936 _ndaugs =3;
937 daugs[0]=EvtPDL::getId(std::string("J/psi"));
938 daugs[1]=EvtPDL::getId(std::string("pi+"));
939 daugs[2]=EvtPDL::getId(std::string("pi-"));
940 }else if(mode == 91){
941 _ndaugs =3;
942 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
943 daugs[1]=EvtPDL::getId(std::string("pi+"));
944 daugs[2]=EvtPDL::getId(std::string("pi-"));
945 }else if(mode == 92){
946 _ndaugs =3;
947 daugs[0]=EvtPDL::getId(std::string("J/psi"));
948 daugs[1]=EvtPDL::getId(std::string("K+"));
949 daugs[2]=EvtPDL::getId(std::string("K-"));
950 }else if(mode == 93){
951 _ndaugs =2;
952 daugs[0]=EvtPDL::getId(std::string("D_s+"));
953 daugs[1]=EvtPDL::getId(std::string("D_s-"));
954 }else if(mode == 94){
955 _ndaugs =2;
956 daugs[0]=EvtPDL::getId(std::string("D_s*+"));
957 daugs[1]=EvtPDL::getId(std::string("D_s-"));
958 }else if(mode == 95){
959 _ndaugs =2;
960 daugs[0]=EvtPDL::getId(std::string("D_s*-"));
961 daugs[1]=EvtPDL::getId(std::string("D_s+"));
962 }else if(mode == 96){
963 _ndaugs =2;
964 daugs[0]=EvtPDL::getId(std::string("Lambda_c+"));
965 daugs[1]=EvtPDL::getId(std::string("anti-Lambda_c-"));
966 }else if(mode == 74100){
967 _ndaugs =1;
968 daugs[0]=EvtPDL::getId(std::string("J/psi"));
969 }else if(mode == 74101){
970 _ndaugs =1;
971 daugs[0]=EvtPDL::getId(std::string("psi(2S)"));
972 }else if(mode == 74102){
973 _ndaugs =1;
974 daugs[0]=EvtPDL::getId(std::string("psi(3770)"));
975 }else if(mode == 74103){
976 _ndaugs =1;
977 daugs[0]=EvtPDL::getId(std::string("phi"));
978 }else if(mode == 74104){
979 _ndaugs =1;
980 daugs[0]=EvtPDL::getId(std::string("omega"));
981 }else if(mode == 74105){
982 _ndaugs =1;
983 daugs[0]=EvtPDL::getId(std::string("rho0"));
984 }else if(mode == 74106){
985 _ndaugs =1;
986 daugs[0]=EvtPDL::getId(std::string("rho(3S)0"));
987 }else if(mode == 74107){
988 _ndaugs =1;
989 daugs[0]=EvtPDL::getId(std::string("omega(2S)"));
990 }else if(mode == 74110){
991 _ndaugs =1;
992 daugs[0]=EvtPDL::getId(std::string("gamma*"));
993 EvtId myvpho = EvtPDL::getId(std::string("gamma*"));
994 EvtPDL::reSetMass(myvpho,mhds*0.9); //for calculating maxXS, it will be resent when this mode is selected
995 _modeFlag.clear();
996 _modeFlag.push_back(74110); //R-value generator tag
997 _modeFlag.push_back(0); //to contain the mode selected
998 }else if(mode == -1){
999 _ndaugs = nson;
1000 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1001 std::cout<<"The paraent particle: "<<EvtPDL::name(pid)<<std::endl;
1002 }else if(mode == -2){
1003 _ndaugs = nson;
1004 for(int i=0;i<nson;i++){ daugs[i]=son[i]; }
1005 }else if(mode == 10000){//use for check ISR
1006 _ndaugs =2;
1007 daugs[0]=EvtPDL::getId(std::string("pi+"));
1008 daugs[1]=EvtPDL::getId(std::string("pi-"));
1009 }else {
1010 std::cout<<"Bad mode_index number " <<mode<<", please refer to the manual."<<std::endl;
1011 ::abort();
1012 }
1013 // std::cout<<_ndaugs<<", "<<daugs[0]<<daugs[1]<<std::endl;
1014
1015 if(VISR){
1016 _ndaugs=2;
1017 daugs[0]=getDaugs()[0];//EvtPDL::getId(std::string("gamma"));
1018 daugs[1]=getDaugs()[1];//EvtPDL::getId(std::string("gamma*"));
1019 }
1020
1021}
1022
1023
1024void EvtConExc::initProbMax(){
1025
1026 noProbMax();
1027
1028}
1029
1030void EvtConExc::decay( EvtParticle *p ){
1031 //std::cout<<"_cms= "<<_cms<<std::endl;
1032 EvtId myvpho = EvtPDL::getId(std::string("vpho"));
1033 if(myvpho != p->getId()){
1034 std::cout<<"Parent needs to be vpho, but found "<<EvtPDL::name(p->getId())<<std::endl;
1035 abort();
1036 }
1037 //for fill root tree
1038 EvtVector4R vgam,hd1,hd2,vhds;
1039
1040 //first for Rscan generator with _mode=74110
1041 if(_mode==74110){
1042 //preparation of mode sampling
1043 std::vector<int> vmod; vmod.clear();
1044 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
1045 50,51,
1046 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
1047 90,91,92,93,94,95,96,
1048 74100,74101,74102,74103,74104,74105,74110};
1049 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
1050 50,51,
1051 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
1052 90,91,92,93,94,95,96,
1053 74100,74101,74102,74103,74104,74105,74110};
1054 double mx = p->getP4().mass();
1055 if(_cms>2.5){
1056 for(int i=0;i<78;i++){vmod.push_back(mn[i]);}
1057 }else{
1058 for(int i=0;i<79;i++){vmod.push_back(mn2[i]);}
1059 }
1060 int mymode;
1061 double pm= EvtRandom::Flat(0.,1);
1062
1063 //std::cout<<"AF[598], AF[598] "<<AF[598]<<" "<<AF[599]<<std::endl;
1064 //--
1065 //-- where AF[599]: cross section with full region (0-Egamma^max GeV), AF[598]: cross section within Egam=(0.025-Egamma^max GeV)
1066 if(pm <_xs0/(_xs0 + _xs1) ){//without ISR photon for mode=74110
1067
1068 mhds=_cms;
1069 mymode=selectMode(vmod,mhds);
1070 _selectedMode = mymode;
1071 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
1072 delete myxsection;
1073 myxsection =new EvtXsection (mymode);
1074 init_mode(mymode);
1075 resetResMass();
1076
1077
1078 if(_ndaugs>1){//for e+e- -> exclusive decays
1079 checkA:
1080 p->makeDaughters(_ndaugs,daugs);
1081 double totMass=0;
1082 for(int i=0;i< _ndaugs;i++){
1083 EvtParticle* di=p->getDaug(i);
1084 double xmi=EvtPDL::getMass(di->getId());
1085 di->setMass(xmi);
1086 totMass+=xmi;
1087 }
1088 if(totMass>p->mass()) goto checkA;
1089 p->initializePhaseSpace( _ndaugs , daugs);
1090 if(!checkdecay(p)) goto checkA;
1091 vhds = p->getP4();
1092 if(_cms>2.5 && !angularSampling(p)) goto checkA;
1093 if(_cms<2.5 && !photonSampling(p)) goto checkA;
1094 }else{// for e+e- -> vector resonace, add a very soft photon
1095 EvtId mydaugs[2];
1096 mydaugs[0]=daugs[0];
1097 EvtPDL::reSetMass(mydaugs[0],mhds*0.9999);
1098 //EvtPDL::reSetWidth(mydaugs[0],0);
1099 mydaugs[1]=gamId; //ISR photon
1100 p->makeDaughters(2,mydaugs);
1101 checkB:
1102 double totMass=0;
1103 for(int i=0;i< 2;i++){
1104 EvtParticle* di=p->getDaug(i);
1105 double xmi=EvtPDL::getMass(di->getId());
1106 di->setMass(xmi);
1107 totMass+=xmi;
1108 }
1109 //std::cout<<"checkB "<<totMass<<" p->mass() "<<p->mass()<<" "<<checkdecay(p)<<std::endl;
1110 if(totMass>p->mass()) goto checkB;
1111 p->initializePhaseSpace(2,mydaugs);
1112
1113 if(!checkdecay(p)) goto checkB;
1114 vhds = p->getDaug(0)->getP4();;
1115 vgam = p->getDaug(1)->getP4();
1116 }
1117 //-----
1118 }else{// with ISR photon for mode=74110
1119 mhds=Mhad_sampling(MH,AF);
1120 double xeng=1-mhds*mhds/(_cms*_cms); //photon energy ratio
1121 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
1122
1123 if(mydbg) mass2=mhds;
1124
1125 //generating events
1126 ModeSelection:
1127 EvtGlobalSet::ConExcPythia=1;
1128 mymode=selectMode(vmod,mhds);
1129 conexcmode = mymode; //save for model REXC ( see EvtRexc.cc) decay
1130 if(mhds<2.3 && mymode==74110) {goto ModeSelection;} // E<2.3 GeV, fully exclusive production
1131 std::cout<<"A selected mode is "<<mymode<<" with Mhds= "<<mhds<<std::endl; //debugging
1132 _selectedMode = mymode;
1133 delete myxsection;
1134 myxsection =new EvtXsection (mymode);
1135 init_mode(mymode);
1136 EvtId myvpho = EvtPDL::getId(std::string("vgam"));
1137 EvtPDL::reSetMass(myvpho,mhds); //set to continum cms energy
1138 EvtPDL::reSetWidth(myvpho,0);
1139
1140 //debugging
1141 //for(int i=0;i<_ndaugs+1;i++){std::cout<<_ndaugs<<" "<<EvtPDL::name(daugs[i])<<std::endl;}
1142
1143 //decay e+e- ->gamma_ISR gamma*
1144 EvtId mydaugs[2];
1145 if(_ndaugs>1) { //gamma* -> exclusive decays
1146 resetResMass();
1147 mydaugs[0]=myvpho;
1148 }else{// vector resonance decays
1149 resetResMass();
1150 mydaugs[0]=daugs[0];
1151 EvtPDL::reSetMass(mydaugs[0],mhds);
1152 //EvtPDL::reSetWidth(mydaugs[0],0);
1153 }
1154 mydaugs[1]=gamId; //ISR photon
1155 int maxflag=0;
1156 int trycount=0;
1157 ISRphotonAng_sampling:
1158 double totMass=0;
1159 p->makeDaughters(2,mydaugs);
1160 for(int i=0;i< 2;i++){
1161 EvtParticle* di=p->getDaug(i);
1162 double xmi=EvtPDL::getMass(di->getId());
1163 di->setMass(xmi);
1164 totMass+=xmi;
1165 }
1166 if(totMass>p->mass()) goto ISRphotonAng_sampling;
1167 //std::cout<<"Ndaugs= "<<_ndaugs<<" Mhds= "<<mhds<<std::endl;
1168 double weight1 = p->initializePhaseSpace(2,mydaugs);
1169 if(!checkdecay(p)) goto ISRphotonAng_sampling;
1170 //set the ISR photon angular distribution
1171 SetP4Rvalue(p,mhds,xeng,theta); //end of decay e+e- -> vpho gamma_ISR
1172
1173 if(maxflag==0) {findMaxXS(mhds); maxflag=1;}
1174 vhds = p->getDaug(0)->getP4();
1175 vgam=p->getDaug(1)->getP4();
1176 double gx=vgam.get(1);
1177 double gy=vgam.get(2);
1178 double sintheta= sqrt(gx*gx+gy*gy)/vgam.d3mag();
1179 bool gam_ang=gam_sampling(mhds,sintheta);
1180 trycount++;
1181
1182 } // with and without ISR sampling block
1183 // finish event generation
1184 // for debugging
1185 if(mydbg){
1186 pgam[0]=vgam.get(0);
1187 pgam[1]=vgam.get(1);
1188 pgam[2]=vgam.get(2);
1189 pgam[3]=vgam.get(3);
1190
1191 phds[0]=vhds.get(0);
1192 phds[1]=vhds.get(1);
1193 phds[2]=vhds.get(2);
1194 phds[3]=vhds.get(3);
1195 costheta = vgam.get(3)/vgam.d3mag();
1196 selectmode = mymode;
1197 xs->Fill();
1198 //std::cout<<"mass1= "<<mass1<<" mass2= "<<mass2<<std::endl;
1199 }
1200 _modeFlag[1]= _selectedMode;
1201 p->setIntFlag(_modeFlag);
1202 return;
1203 }//end block if(_mode==74110)
1204
1205 //=======================================================
1206 // e+e- -> gamma_ISR gamma*
1207 //=======================================================
1208 if(VISR){
1209 EvtId gid=EvtPDL::getId("gamma*");
1210 double pm= EvtRandom::Flat(0.,1);
1211
1212 if(pm <_xs0/(_xs0 + _xs1) ){//with a very soft ISR photon
1213 EvtPDL::reSetMass(gid,_cms*0.995);
1214 mhds = _cms*0.995;
1215 }else{
1216 mhds=Mhad_sampling(MH,AF);
1217 EvtPDL::reSetMass(gid,mhds);
1218 }
1219 //debugging
1220 std::cout<<"generatedMass "<<EvtPDL::getMeanMass(gid)<<std::endl;
1221 double xeng=1-mhds*mhds/(_cms*_cms); //photon energy ratio
1222 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
1223 p->makeDaughters(2,daugs);
1224 for(int i=0;i< 2;i++){
1225 EvtParticle* di=p->getDaug(i);
1226 double xmi=EvtPDL::getMeanMass(di->getId());
1227 di->setMass(xmi);
1228 }
1229 p->initializePhaseSpace(2,daugs);
1230 SetP4(p,mhds,xeng,theta);
1231 return;
1232 }
1233
1234
1235 //========================================================
1236 //=== for other mode generation with _mode != 74110
1237 //========================================================
1238
1239 if((_xs0 + _xs1)==0) {std::cout<<"EvtConExc::zero total xsection"<<std::endl;::abort();}
1240 double pm= EvtRandom::Flat(0.,1);
1241 double xsratio = _xs0/(_xs0 + _xs1);
1242 int iflag=2; //flag = 0 for ee->hadrons, 1 for ee->gamma_ISR hadrons, 2: mix 0 and 1 case
1243 if(getArg(0)!= -2){// exclude DIY case
1244 if(getNArg()==3 && radflag==1){iflag=1;xsratio=0;} // only generating gamma hadrons mode
1245 else if(getNArg()==3 && radflag==0) {iflag=0;xsratio=1;}
1246 }
1247
1248 // std::cout<<"xsratio= "<<xsratio<<std::endl;
1249
1250 if(pm<xsratio ){// no ISR photon
1251
1252 p->makeDaughters(_ndaugs,daugs);
1253 for(int i=0;i< _ndaugs;i++){
1254 EvtParticle* di=p->getDaug(i);
1255 double xmi=EvtPDL::getMass(di->getId());
1256 di->setMass(xmi);
1257 }
1258 angular_hadron:
1259 p->initializePhaseSpace(_ndaugs,daugs);
1260 bool byn_ang;
1261 EvtVector4R pd1 = p->getDaug(0)->getP4();
1262 EvtVector4R pcm(_cms,0,0,0);
1263 if(_ndaugs==2){//angular distribution for the hadrons only for two-body decays
1264 byn_ang=hadron_angle_sampling(pd1, pcm);
1265 if(!byn_ang) goto angular_hadron;
1266 }
1267
1268 //for histogram
1269 cosp = pd1.get(3)/pd1.d3mag();
1270 mhds = _cms;
1271 //std::cout<<"cosp, mhds "<<cosp<<" "<<mhds<<std::endl;
1272 //p->printTree();
1273
1274 }else{// sampling m_hadrons and decay e+e- ->gamma gamma*
1275 double mhdr=Mhad_sampling(MH,AF);
1276 double xeng=1-mhdr*mhdr/(_cms*_cms); //photon energy ratio
1277 double theta=ISR_ang_sampling(xeng); //ISR photon polar angle in rad unit
1278 EvtId gid =EvtPDL::getId("vhdr");
1279 EvtPDL::reSetMass(gid,mhdr);
1280 int ndaugs =2;
1281 EvtId mydaugs[2];
1282 mydaugs[0] =EvtPDL::getId(std::string("gamma"));
1283 mydaugs[1]=EvtPDL::getId(std::string("vhdr"));
1284
1285 //for histogram
1286 mhds = mhdr;
1287 costheta = cos(theta);
1288 //--
1289
1290 p->makeDaughters(2,mydaugs);
1291 for(int i=0;i< 2;i++){
1292 EvtParticle* di=p->getDaug(i);
1293 double xmi=EvtPDL::getMass(di->getId());
1294 di->setMass(xmi);
1295 }
1296 p->initializePhaseSpace(2,mydaugs);
1297 SetP4(p,mhdr,xeng,theta); //end of decay e+e- -> gamma_ISR gamma*
1298 //p->printTree();
1299
1300 //----
1301 }//end of gamma gamma*, gamma*->hadrons generation
1302
1303
1304 //-----------------------------------------
1305 // End of decays for all topology
1306 //----------------------------------------
1307 //================================= event analysis
1308 if(_nevt ==0 ){
1309 std::cout<<"The decay chain: "<<std::endl;
1310 p->printTree();
1311 }
1312 _nevt++;
1313 //--- for debuggint to make root file
1314 if(mydbg){
1315 xs->Fill();
1316 }
1317
1318 //----
1319 return ;
1320}
1321
1322bool EvtConExc::hadron_angle_sampling(EvtVector4R ppi, EvtVector4R pcm){
1323 EvtVector4R pbst=-1*pcm;
1324 pbst.set(0,pcm.get(0));
1325 EvtVector4R p4i = boostTo(ppi,pbst);
1326 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
1327 bool byn_ang = baryon_sampling(pcm, p4i);
1328 return byn_ang;
1329 }else if(_mode ==6||_mode==45 || _mode==46 || _mode==70 || _mode==71|| _mode==93){// ee->PP (pi+pi-,..) mode
1330 bool msn_ang = meson_sampling(pcm,p4i);
1331 return msn_ang;
1332 }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){
1333 bool msn_ang = VP_sampling(pcm,p4i);
1334 return msn_ang;
1335 }
1336 return true;
1337}
1338
1339
1340double EvtConExc::gamHXSection(EvtParticle *p,double El,double Eh,int nmc){//El, Eh : the lower and higher energy of radiative photons
1341 //std::cout<<"nmc= "<<nmc<<std::endl;
1342 gamId = EvtPDL::getId(std::string("gamma"));
1343 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
1344 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
1345 double totxs = 0;
1346 maxXS=-1;
1347 _er1=0;
1348 Rad2Xs =0;
1349 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
1350 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
1351 int gamdaugs = _ndaugs+1;
1352
1353 EvtId theDaugs[10];
1354 for(int i=0;i<_ndaugs;i++){
1355 theDaugs[i] = daugs[i];
1356 }
1357 theDaugs[_ndaugs]=gamId;
1358
1359 gamH->makeDaughters(gamdaugs,theDaugs);
1360
1361 for(int i=0;i<gamdaugs;i++){
1362 EvtParticle* di=gamH->getDaug(i);
1363 double xmi=EvtPDL::getMass(di->getId());
1364 di->setMass(xmi);
1365 }
1366 loop:
1367 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
1368 //-- slection:theta_gamma > 1 degree
1369 EvtVector4R pgam = gamH -> getDaug(_ndaugs)->getP4();
1370 double pmag = pgam.d3mag();
1371 double pz = pgam.get(3);
1372 double sintheta = sqrt( pmag*pmag - pz*pz)/pmag;
1373 double onedegree = 1./180.*3.1415926;
1374 //if(sintheta < onedegree) {goto loop;}
1375 if(pmag <El || pmag >Eh) {goto loop;}
1376 //--------
1377 // std::cout<<"pmag= "<<pmag<<std::endl;
1378
1379 double thexs = difgamXs(gamH); //prepare the photon angular distribution
1380 Rad2Xs += Rad2difXs( gamH );
1381 if(thexs>maxXS) {maxXS=thexs;}
1382 double s = p_init.mass2();
1383 double x = 2*pgam.get(0)/sqrt(s);
1384 double rad1xs = Rad1difXs(gamH); //first order xs by Eq(4)hep-ph/9910523
1385 totxs += rad1xs;
1386 _er1 += differ;
1387 gamH->deleteDaughters();
1388 } //for int_i block
1389 _er1/=nmc;
1390
1391 Rad2Xs/=nmc; // the leading second order correction
1392 totxs/=nmc; // first order correction XS
1393
1394 // return totxs; // first order correction XS
1395 return Rad2Xs; // the leading second order correction
1396}
1397
1398
1399double EvtConExc::gamHXSection(double s, double El,double Eh,int nmc){//El, Eh : the lower and higher energy of radiative photons
1400 //--for narrow resonance psi(2S),J/psi, phi, which can not calculated with the MC integration
1401 //double mxL = sqrt( s-2*Eh*sqrt(s) ); //low mass
1402 //double mxH = sqrt( s-2*El*sqrt(s) ); //high mass
1403 //double sigv = narrowRXS(mxL,mxH);
1404 //--
1405
1406 double totxs = 0;
1407 maxXS=-1;
1408 _er1=0;
1409 Rad2Xs =0;
1410 double xL=2*El/sqrt(s);
1411 double xH=2*Eh/sqrt(s);
1412 for(int i=0;i<nmc;i++){//It is found that the narrow resonance can not included in this MC integartion
1413 double rdm = EvtRandom::Flat(0.,1.);// set angular cut 1^o
1414 double x= xL+ (xH-xL)*rdm;
1415 Rad2Xs += Rad2difXs(s,x);
1416 _er1 += differ2; //stored when calculate Rad2Xs
1417 // std::cout<<"i= "<<i<<" Rad2Xs= "<<Rad2Xs<<std::endl;
1418 }
1419 _er1/=nmc;
1420 _er1*=(xH-xL);
1421 //std::cout<<"_er1= "<<_er1<<std::endl;
1422 Rad2Xs/=nmc; // the leading second order correction
1423 double thexs= Rad2Xs*(xH-xL); //xH-xL is the Jacobi factor
1424 //std::cout<<"thexs= "<<thexs<<std::endl;
1425 //if(sigv != -1) thexs += sigv; //add narrow resonance ISR cross section
1426 return thexs;
1427}
1428
1429
1430
1431double EvtConExc::gamHXSection(double El,double Eh){// using Gaussian integration subroutine from Cern lib
1432 std::cout<<" "<<std::endl;
1433 extern double Rad2difXs(double *x);
1434 extern double Rad2difXs2(double *x);
1435 double eps = 0.01;
1436 double Rad2Xs;
1437 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
1438 if(Rad2Xs==0){
1439 for(int iter=0;iter<10;iter++){
1440 eps += 0.01;
1441 if(_mode!=-2) {Rad2Xs= dgauss_(&Rad2difXs,&El,&Eh,&eps);}else{Rad2Xs = dgauss_(&Rad2difXs2,&El,&Eh,&eps);}
1442 if(!Rad2Xs) return Rad2Xs;
1443 }
1444 }
1445 return Rad2Xs; // the leading second order correction
1446}
1447
1448double EvtConExc::gamHXSection_er(double El,double Eh){// using Gaussian integration subroutine from Cern lib
1449 std::cout<<" "<<std::endl;
1450 extern double Rad2difXs_er(double *x);
1451 extern double Rad2difXs_er2(double *x);
1452 double eps = 0.01;
1453 double Rad2Xs;
1454 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
1455 if(Rad2Xs==0){
1456 for(int iter=0;iter<10;iter++){
1457 eps += 0.01;
1458 if(_mode !=-2){Rad2Xs = dgauss_(&Rad2difXs_er,&El,&Eh,&eps);}else{ Rad2Xs = dgauss_(&Rad2difXs_er2,&El,&Eh,&eps);}
1459 if(!Rad2Xs) return Rad2Xs;;
1460 }
1461 }
1462 return Rad2Xs; // the leading second order correction
1463}
1464
1465
1466void EvtConExc::findMaxXS(EvtParticle *p){
1467 //std::cout<<"nmc= "<<nmc<<std::endl;
1468 gamId = EvtPDL::getId(std::string("gamma"));
1469 EvtId xvp = EvtPDL::getId(std::string("xvpho"));
1470 EvtVector4R p_init( p->getP4().mass(),0.,0.,0.);
1471 double totxs = 0;
1472 maxXS=-1;
1473 _er1=0;
1474 Rad2Xs =0;
1475 int nmc = 50000;
1476 for(int ii=0;ii<nmc;ii++){// integarted the gamma hadrons xsection
1477 gamH=EvtParticleFactory::particleFactory(xvp, p_init);
1478 int gamdaugs = _ndaugs+1;
1479
1480 EvtId theDaugs[10];
1481 for(int i=0;i<_ndaugs;i++){
1482 theDaugs[i] = daugs[i];
1483 }
1484 theDaugs[_ndaugs]=gamId;
1485
1486 gamH->makeDaughters(gamdaugs,theDaugs);
1487 loop:
1488 double totm=0;
1489 for(int i=0;i<gamdaugs;i++){
1490 EvtParticle* di=gamH->getDaug(i);
1491 double xmi=EvtPDL::getMass(di->getId());
1492 di->setMass(xmi);
1493 totm += xmi;
1494 }
1495
1496 //std::cout<<totm<<" "<<p_init.get(0)<<std::endl;
1497 if(totm >= p_init.get(0) ) goto loop;
1498
1499 double weight = gamH->initializePhaseSpace( gamdaugs , theDaugs);
1500 double thexs = difgamXs(gamH); //prepare the photon angular distribution
1501 EvtVector4R p4gam = gamH->getDaug(_ndaugs)->getP4();
1502 double costheta = p4gam.get(3)/p4gam.d3mag();
1503 double sintheta = sqrt(1-costheta*costheta);
1504 bool acut=(sintheta>0.04) && (sintheta<0.96); //set photon auglar cut 2^o
1505 if(thexs>maxXS && acut ) {maxXS=thexs;}
1506 //gamH->deleteTree();
1507 }
1508
1509}
1510
1511double EvtConExc::difgamXs(EvtParticle *p){// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
1512 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
1513 EvtVector4R p0 = p->getDaug(0)->getP4();
1514 for(int i=1;i<_ndaugs;i++){
1515 p0 += p->getDaug(i)->getP4();
1516 }
1517 EvtVector4R pgam = p->getDaug(_ndaugs)->getP4();
1518 double mhs = p0.mass();
1519 double egam = pgam.get(0);
1520 double sin2theta = pgam.get(3)/pgam.d3mag();
1521 sin2theta = 1-sin2theta*sin2theta;
1522
1523 double cms = p->getP4().mass();
1524 double alpha = 1./137.;
1525 double pi = 3.1415926;
1526 double x = 2*egam/cms;
1527 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
1528 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
1529
1530 double xs = myxsection->getXsection(mhs);
1531 double difxs = 2*mhs/cms/cms * wsx *xs;
1532 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
1533 return difxs;
1534
1535}
1536
1537
1538bool EvtConExc::gam_sampling(EvtParticle *p){//photon angular distribution sampling
1539 double pm= EvtRandom::Flat(0.,1);
1540 double xs = difgamXs( p );
1541 double xsratio = xs/maxXS;
1542 if(pm<xsratio){return true;}else{return false;}
1543}
1544
1545bool EvtConExc::gam_sampling(double mhds,double sintheta){
1546 double pm= EvtRandom::Flat(0.,1);
1547 double xs = difgamXs( mhds,sintheta );
1548 double xsratio = xs/maxXS;
1549 if(pm<xsratio){return true;}else{return false;}
1550}
1551
1552bool EvtConExc::xs_sampling(double xs){
1553 double pm= EvtRandom::Flat(0.,1);
1554 //std::cout<<"Rdm= "<<pm<<std::endl;
1555 if(pm <xs/XS_max) {return true;} else {return false;}
1556}
1557
1558bool EvtConExc::xs_sampling(double xs,double xs0){
1559 double pm= EvtRandom::Flat(0.,1);
1560 // std::cout<<"Rdm= "<<pm<<std::endl;
1561 if(pm <xs/(xs0*1.1)) {return true;} else {return false;}
1562}
1563
1564bool EvtConExc::baryon_sampling(EvtVector4R pcm, EvtVector4R pi){
1565 EvtHelSys angles(pcm,pi); //using helicity sys.angles
1566 double theta=angles.getHelAng(1);
1567 double phi =angles.getHelAng(2);
1568 double costheta=cos(theta); //using helicity angles in parent system
1569
1570 // double costh = pcm.dot(pi)/pcm.d3mag()/pi.d3mag();
1571 //std::cout<<"costheta: "<<costheta<<", "<<costh<<std::endl;
1572 double alpha = baryonAng(_cms);
1573 if(_mode ==96){alpha=-0.34;}
1574 double pm= EvtRandom::Flat(0.,1);
1575 double ang = 1 + alpha*costheta*costheta;
1576 double myang;
1577 if(alpha>=0){myang=1+alpha;}else{myang=1;}
1578 if(pm< ang/myang) {return true;}else{return false;}
1579}
1580
1581bool EvtConExc::meson_sampling(EvtVector4R pcm, EvtVector4R pi){
1582 EvtHelSys angles(pcm,pi); //using helicity sys.angles
1583 double theta=angles.getHelAng(1);
1584 double phi =angles.getHelAng(2);
1585 double costheta=cos(theta); //using helicity angles in parent system
1586
1587 double pm= EvtRandom::Flat(0.,1);
1588 double ang = 1 - costheta*costheta;
1589 if(pm< ang/1.) {return true;}else{return false;}
1590}
1591
1592bool EvtConExc::VP_sampling(EvtVector4R pcm, EvtVector4R pi){
1593 EvtHelSys angles(pcm,pi); //using helicity sys.angles
1594 double theta=angles.getHelAng(1);
1595 double phi =angles.getHelAng(2);
1596 double costheta=cos(theta); //using helicity angles in parent system
1597
1598 double pm= EvtRandom::Flat(0.,1);
1599 double ang = 1 + costheta*costheta;
1600 if(pm< ang/2.) {return true;}else{return false;}
1601}
1602
1603double EvtConExc::Rad1(double s, double x){ //first order correction
1604 //radiator correction upto second order, see Int. J. of Moder.Phys. A
1605 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
1606 double alpha = 1./137.;
1607 double pi=3.1415926;
1608 double me = 0.5*0.001; //e mass
1609 double sme = sqrt(s)/me;
1610 double L = 2*log(sme);
1611 double wsx = 2*alpha/pi/x*(L-1)*(1 - x + x*x/2 );
1612 return wsx;
1613}
1614
1615double EvtConExc::Rad2(double s, double x){
1616 //radiator correction upto second order, see Int. J. of Moder.Phys. A, hep-ph/9910523, Eq(8)
1617 // Spectroscopy at B-factory using hard emision photon, by M. Benayoun, et al
1618 double alpha = 1./137.;
1619 double pi=3.1415926;
1620 double me = 0.5*0.001; //e mass
1621 double xi2 = 1.64493407;
1622 double xi3=1.2020569;
1623 double sme = sqrt(s)/me;
1624 double L = 2*log(sme);
1625 double beta = 2*alpha/pi*(L-1);
1626 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.;
1627 double ap = alpha/pi;
1628 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
1629 double wsx = Delta * beta * pow(x,beta-1)-0.5*beta*(2-x);
1630 double wsx2 = (2-x)*( 3*log(1-x)-4*log(x) ) -4*log(1-x)/x-6+x;
1631 wsx = wsx + beta*beta/8. *wsx2;
1632 double mymx = sqrt(s*(1-x));
1633 // return wsx*vph; //getVP(mymx);//vph is vaccum polarization factor
1634 return wsx*getVP(mymx);//vph is vaccum polarization factor
1635}
1636
1637
1638
1639double EvtConExc::Rad2difXs(EvtParticle *p){// leading second order xsection
1640 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
1641 double summass = p->getDaug(0)->getP4().mass();
1642 EvtVector4R v4p=p->getDaug(0)->getP4();
1643 for(int i=1;i<_ndaugs;i++){
1644 summass += p->getDaug(i)->getP4().mass();
1645 v4p += p->getDaug(i)->getP4();
1646 }
1647
1648 double Egam = p->getDaug(_ndaugs)->getP4().d3mag();
1649 double cms = p->getP4().mass();
1650 double mrg = cms - summass;
1651 double pm= EvtRandom::Flat(0.,1);
1652 double mhs = pm*mrg + summass;
1653
1654 double s = cms * cms;
1655 double x = 2*Egam/cms;
1656 //double mhs = sqrt( s*(1-x) );
1657 double wsx = Rad2(s,x);
1658
1659 // std::cout<<"random : "<<pm<<std::endl;
1660
1661 double xs = myxsection->getXsection(mhs);
1662 if(xs>XS_max){XS_max = xs;}
1663 double difxs = 2*mhs/cms/cms * wsx *xs;
1664 differ2 = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
1665 differ *= mrg; //Jacobi factor for variable m_hds
1666 difxs *= mrg;
1667 return difxs;
1668}
1669
1670double EvtConExc::Rad2difXs(double s, double x ){// leading second order xsection
1671 double wsx = Rad2(s,x);
1672 double mhs = sqrt(s*(1-x));
1673 double xs = myxsection->getXsection(mhs);
1674 //if(testflag==1)std::cout<<"s= "<<s<<" mhs= "<<mhs<<" x= "<<x<<" xs= "<<xs<<std::endl;
1675 if(xs>XS_max){XS_max = xs;}
1676 double difxs = wsx *xs;
1677 differ2 = wsx *(myxsection->getErr(mhs));
1678 return difxs;
1679}
1680
1681
1682double EvtConExc::Rad1difXs(EvtParticle *p){// // first order xsection
1683 // where p -> hdarons(0 ~ _ndaugs-1) + gamma
1684 double summass = p->getDaug(0)->getP4().mass();
1685 for(int i=1;i<_ndaugs;i++){
1686 summass += p->getDaug(i)->getP4().mass();
1687 }
1688
1689 double cms = p->getP4().mass();
1690 double mrg = cms - summass;
1691 double pm= EvtRandom::Flat(0.,1);
1692 double mhs = pm*mrg + summass;
1693
1694 double s = cms * cms;
1695 double x = 1 - mhs*mhs/s;
1696 double wsx = Rad1(s,x);
1697
1698 // std::cout<<"pgam"<<pgam<<" mhs,egam,sin2theta,cms,x,wsx"<< mhs<<", "<<egam<<", "<<sin2theta<<", "<<cms<<", "<<x<<", "<<wsx<<std::endl;
1699
1700 double xs = myxsection->getXsection(mhs);
1701 if(xs>XS_max){XS_max = xs;}
1702 double difxs = 2*mhs/cms/cms * wsx *xs;
1703
1704 differ = 2*mhs/cms/cms * wsx *(myxsection->getErr(mhs));
1705 differ *= mrg; //Jacobi factor for variable m_hds
1706 difxs *= mrg;
1707 return difxs;
1708}
1709
1710void EvtConExc::init_Br_ee(){
1711 // double psipp_ee=9.6E-06;
1712 double psip_ee =7.73E-03;
1713 double jsi_ee =5.94*0.01;
1714 double phi_ee =2.954E-04;
1715 // double omega_ee=7.28E-05;
1716 // double rho_ee = 4.72E-05;
1717 EvtId psppId=EvtPDL::getId(std::string("psi(3770)"));
1718 EvtId pspId=EvtPDL::getId(std::string("psi(2S)"));
1719 EvtId jsiId=EvtPDL::getId(std::string("J/psi"));
1720 EvtId phiId=EvtPDL::getId(std::string("phi"));
1721 EvtId omegaId=EvtPDL::getId(std::string("omega"));
1722 EvtId rhoId=EvtPDL::getId(std::string("rho0"));
1723 BR_ee.clear();
1724 ResId.clear();
1725
1726 //BR_ee.push_back(rho_ee);
1727 //BR_ee.push_back(omega_ee);
1728 BR_ee.push_back(phi_ee);
1729 BR_ee.push_back(jsi_ee);
1730 BR_ee.push_back(psip_ee);
1731 //BR_ee.push_back(psipp_ee);
1732
1733 //ResId.push_back(rhoId);
1734 //ResId.push_back(omegaId);
1735 ResId.push_back(phiId);
1736 ResId.push_back(jsiId);
1737 ResId.push_back(pspId);
1738 //ResId.push_back(psppId);
1739
1740}
1741
1742double EvtConExc::Ros_xs(double mx,double bree, EvtId pid){// in unit nb
1743 double pi=3.1415926;
1744 double s= mx*mx;
1745 double width = EvtPDL::getWidth(pid);
1746 double mass = EvtPDL::getMeanMass(pid);
1747 double xv = 1-mass*mass/s;
1748 double sigma = 12*pi*pi*bree*width/mass/s;
1749 sigma *= Rad2(s, xv);
1750 double nbar = 3.89379304*100000; // ! GeV-2 = 3.89379304*10^5 nbar
1751 return sigma*nbar;
1752}
1753
1754
1755double Rad2difXs(double *mhs){// leading second order xsection, mhs: the mass of final states
1756 double cms = EvtConExc::_cms;
1757 double s = cms * cms;
1758 double x = 1 - (*mhs)*(*mhs)/s;
1759 EvtConExc myconexc;
1760 double wsx;
1761 double dhs=(*mhs);
1762 double xs = EvtConExc::myxsection->getXsection(dhs);
1763 wsx=myconexc.Rad2(s,x);
1764 if(xs>EvtConExc::XS_max){EvtConExc::XS_max = xs;}
1765 double difxs = 2*dhs/cms/cms * wsx *xs;
1766 std::ofstream oa;oa<<x<<std::endl; //without this line, the running will breaking, I don't know why
1767 return difxs;
1768}
1769double Rad2difXs_er(double *mhs){// leading second order xsection, mhs: the mass of final states
1770 double cms = EvtConExc::_cms;
1771 double s = cms * cms;
1772 double x = 1 - (*mhs)*(*mhs)/s;
1773 EvtConExc myconexc;
1774 double wsx;
1775 double xs = EvtConExc::myxsection->getErr(*mhs);
1776 wsx=myconexc.Rad2(s,x);
1777 double differ2 = 2*(*mhs)/cms/cms * wsx *xs;
1778 std::ofstream oa;oa<<x<<std::endl;
1779 return differ2;
1780}
1781
1782double Rad2difXs2(double *mhs){// leading second order xsection, mhs: the mass of final states
1783 double cms = EvtConExc::_cms;
1784 double s = cms * cms;
1785 double x = 1 - (*mhs)*(*mhs)/s;
1786 EvtConExc myconexc;
1787 double wsx;
1788 double dhs=(*mhs);
1789 double xs = EvtConExc::myxsection->getXsection(dhs);
1790 wsx=myconexc.Rad2(s,x);
1791 if(xs>EvtConExc::XS_max){EvtConExc::XS_max = xs;}
1792 double difxs = 2*dhs/cms/cms * wsx *xs;
1793 oa<<x<<std::endl; //without this line, the running will breaking, I don't know why
1794 return difxs;
1795}
1796
1797double Rad2difXs_er2(double *mhs){// leading second order xsection, mhs: the mass of final states
1798 double cms = EvtConExc::_cms;
1799 double s = cms * cms;
1800 double x = 1 - (*mhs)*(*mhs)/s;
1801 EvtConExc myconexc;
1802 double wsx;
1803 double xs = EvtConExc::myxsection->getErr(*mhs);
1804 wsx=myconexc.Rad2(s,x);
1805 double differ2 = 2*(*mhs)/cms/cms * wsx *xs;
1806 oa<<x<<std::endl;
1807 return differ2;
1808}
1809
1810
1811double EvtConExc::SoftPhoton_xs(double s, double b){
1812 double alpha = 1./137.;
1813 double pi=3.1415926;
1814 double me = 0.5*0.001; //e mass
1815 double xi2 = 1.64493407;
1816 double xi3=1.2020569;
1817 double sme = sqrt(s)/me;
1818 double L = 2*log(sme);
1819 double beta = 2*alpha/pi*(L-1);
1820 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.;
1821 double ap = alpha/pi;
1822 double Delta = 1 + ap *(1.5*L + 1./3.*pi*pi-2) + ap*ap *delta2;
1823
1824 double beta2 = beta*beta;
1825 double b2 = b*b;
1826
1827 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 -
1828 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) +
1829 16*pow(beta,2)*Li2(b))/32. ;
1830 double mymx = sqrt(s*(1-b));
1831 //return xs*vph; //getVP(mymx); //vph :the vaccum polarzation factor
1832 return xs*getVP(mymx); //vph :the vaccum polarzation factor
1833
1834}
1835
1836double EvtConExc::Li2(double x){
1837 double li2= -x +x*x/4. - x*x*x/9;
1838 return li2;
1839}
1840
1841
1842double EvtConExc::lgr(double *x,double *y,int n,double t)
1843{ int n0=-1;
1844 double z;
1845 for(int i=0;i<n-1;i++){
1846 if(x[i]>=t && t> x[i+1]) {n0=i;break;}
1847 }
1848 if(n0==-1) {return 0.0;}else{
1849 double k=(y[n0]-y[n0+1])/(x[n0]-x[n0+1]);
1850 z= y[n0+1]+k*(t-x[n0+1]);
1851 }
1852 //std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
1853 //std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
1854 return z;
1855}
1856
1857bool EvtConExc::islgr(double *x,double *y,int n,double t)
1858{ int n0=-1;
1859 double z;
1860 for(int i=0;i<n-1;i++){
1861 if(x[i]>=t && t> x[i+1]) {n0=i;break;}
1862 }
1863 double xstotal=y[599];
1864 if(n0==-1) {return 0;}else{
1865 double p1=y[n0]/xstotal;
1866 double p2=y[n0+1]/xstotal;
1867 double pm= EvtRandom::Flat(0.,1);
1868 if(p1<pm && pm<=p2) {return 1;}else{return 0;}
1869 }
1870}
1871
1872double EvtConExc::LLr(double *x,double *y,int n,double t)
1873{ int n0=-1;
1874 double z;
1875 if( t==x[n-1] ) return y[n-1];
1876 for(int i=0;i<n-1;i++){
1877 if(x[i]<=t && t< x[i+1]) {n0=i;break;}
1878 }
1879 if(n0==-1) {return 0.0;}else{
1880 double k=(y[n0]-y[n0+1])/(x[n0]-x[n0+1]);
1881 z= y[n0+1]+k*(t-x[n0+1]);
1882 }
1883 //std::cout<<"t="<<t<<" "<<x[n0+1]<<" "<<x[n0]<<std::endl;
1884 //std::cout<<"t="<<t<<" "<<y[n0+1]<<" "<<y[n0]<<std::endl;
1885 return z;
1886}
1887
1888double EvtConExc::Mhad_sampling(double *x,double *y){//sample ISR hadrons, return Mhadron
1889 //x=MH: array contains the Mhad
1890 //y=AF: array contains the Mhad*Xsection
1891 //n: specify how many bins in the hadron sampling
1892 int n=598; //AF[599] is the total xs including Ecms point
1893 double pm= EvtRandom::Flat(0.,1);
1894 int mybin=1;
1895 double xst=y[n];
1896 for(int i=0;i<n;i++){
1897 if((y[i]/xst)<pm && pm<=(y[i+1]/xst)){mybin=i+1;break;}
1898 }
1899 pm= EvtRandom::Flat(0.,1);
1900 double mhd=x[mybin-1]+(x[mybin]-x[mybin-1])*pm;
1901
1902 if(mhd>_cms) {std::cout<<"selected mhd larger than Ecms "<<mhd<<" > "<<x[mybin]<<std::endl;abort();}
1903 if(mhd<_mhdL){
1904 std::cout<<"the sample mhassrons is less than the low end of XS"<<mhd<<" <"<<_mhdL<<std::endl;
1905 for(int i=0;i<598;i++){std::cout<<i<<" "<<x[i]<<" "<<y[i]<<std::endl;}
1906 abort();
1907 }
1908 return mhd;
1909}
1910
1911double EvtConExc::ISR_ang_integrate(double x,double theta){
1912 //see Eq(11) in arXif:hep-ph/9910523, return h(theta)/h(pi)
1913 double costheta=cos(theta);
1914 double eb=_cms/2;
1915 double cos2=costheta*costheta;
1916 double sin2=1-cos2;
1917 double me=0.51*0.001;
1918 double L=2*log(_cms/me);
1919 double meE2= me*me/eb/eb;
1920 double hpi=L-1;
1921 double hq1= meE2/2*costheta/(sin2+meE2*cos2);
1922 double hq2= -0.5*log((1+sqrt(1-meE2)*costheta)/(1-sqrt(1-meE2)*costheta));
1923 double hq3=x*x*costheta/2/(x*x-2*x+2)*(1-meE2/(sin2+meE2*cos2));
1924 double hq=(L-1)/2.+hq1+hq2+hq3;
1925 hq /= hpi; //Eq (11) in arXif:hep-ph/9910523
1926 return hq;
1927}
1928
1929double EvtConExc::ISR_ang_sampling(double x){
1930 double xx[180],yy[180];
1931 double dgr2rad=1./180.*3.1415926;
1932 xx[0]=0;
1933 xx[1]=5*dgr2rad; //first bin is 5 degree
1934 int nc=2;
1935 for(int i=6;i<=175;i++){ //last bin is 5 degree
1936 xx[nc]=i*dgr2rad;
1937 nc++;
1938 }
1939 xx[nc]=180*dgr2rad;
1940 for(int j=0;j<=nc;j++){
1941 yy[j]=ISR_ang_integrate(x,xx[j]);
1942 }
1943 double pm= EvtRandom::Flat(0.,1);
1944 int mypt=1;
1945 for(int j=1;j<=nc;j++){
1946 if(yy[j-1]<pm && pm<=yy[j]){mypt=j;break;}
1947 }
1948 pm= EvtRandom::Flat(0.,1);
1949 double mytheta= xx[mypt-1]+ (xx[mypt]-xx[mypt-1])*pm;
1950 return mytheta; //theta in rad unit
1951}
1952
1953void EvtConExc::SetP4(EvtParticle *p, double mhdr, double xeng,double theta){ //set the gamma and gamma* momentum according sampled results
1954 double pm= EvtRandom::Flat(0.,1);
1955 double phi=2*3.1415926*pm;
1956 double gamE = _cms/2*xeng;
1957 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
1958 double px = gamE*sin(theta)*cos(phi);
1959 double py = gamE*sin(theta)*sin(phi);
1960 double pz = gamE*cos(theta);
1961 EvtVector4R p4[2];
1962 p4[0].set(gamE,px,py,pz);//ISR photon
1963 p4[1].set(hdrE,-px,-py,-pz); //mhdraons
1964 for(int i=0;i<2;i++){
1965 EvtId myid = p->getDaug(i)->getId();
1966 p->getDaug(i)->init(myid,p4[i]);
1967 }
1968}
1969
1970void EvtConExc::SetP4Rvalue(EvtParticle *p, double mhdr, double xeng,double theta){ //set the gamma and gamma* momentum according sampled results
1971 double pm= EvtRandom::Flat(0.,1);
1972 double phi=2*3.1415926*pm;
1973 double gamE = _cms/2*xeng;
1974 double hdrE = sqrt(mhdr*mhdr+gamE*gamE);
1975 double px = gamE*sin(theta)*cos(phi);
1976 double py = gamE*sin(theta)*sin(phi);
1977 double pz = gamE*cos(theta);
1978 EvtVector4R p4[2];
1979 p4[0].set(hdrE,px,py,pz);//mhdraons
1980 p4[1].set(gamE,-px,-py,-pz); //ISR photon
1981 for(int i=0;i<2;i++){
1982 EvtId myid = p->getDaug(i)->getId();
1983 p->getDaug(i)->init(myid,p4[i]);
1984 }
1985}
1986
1987
1988void EvtConExc::findMaxXS(double mhds ){// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
1989 maxXS=-1;
1990 for(int i=0;i<90000;i++){
1991 double x = 1-(mhds/_cms)*(mhds/_cms);
1992 double cos = EvtRandom::Flat(0.0006,0.9994); //set cut on ISR gamma 2^o
1993 double sin2theta =sqrt(1-cos*cos);
1994 double alpha = 1./137.;
1995 double pi = 3.1415926;
1996 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
1997 double xs = myxsection->getXsection(mhds);
1998 double difxs = 2*mhds/_cms/_cms * wsx *xs;
1999 if(difxs>maxXS)maxXS=difxs;
2000 }
2001}
2002
2003double EvtConExc::difgamXs(double mhds,double sintheta ){// dsigma/dm/dcos\theta* (Eq(1)@PRD73,012005)
2004 double x = 1-(mhds/_cms)*(mhds/_cms);
2005 double sin2theta = sintheta*sintheta;
2006 double alpha = 1./137.;
2007 double pi = 3.1415926;
2008 double wsx = alpha/pi/x*( (2-2*x+x*x)/sin2theta-x*x/2. );
2009 double xs = myxsection->getXsection(mhds);
2010 double difxs = 2*mhds/_cms/_cms * wsx *xs;
2011 return difxs;
2012}
2013
2014int EvtConExc::selectMode(std::vector<int> vmod, double mhds){
2015 //first get xsection for each mode in vmod array
2016 double uscale = 1;
2017 std::vector<double> vxs;vxs.clear();
2018 for (int i=0;i<vmod.size();i++){
2019 int imod = vmod[i];
2020 delete myxsection;
2021 myxsection =new EvtXsection (imod);
2022 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
2023 //if(uscale!=1) std::cout<<"mode="<<imod<<" uscale= "<<uscale<<std::endl;
2024
2025 double myxs=uscale*myxsection->getXsection(mhds);
2026 if(imod==0) {vxs.push_back(myxs);}
2027 else if(imod==74110){//for continuum process
2028 double xcon = myxs - vxs[i-1];
2029 if(xcon<0) {xcon=vxs[i-1];}else{xcon=myxs;}
2030 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
2031 vxs.push_back(xcon);
2032 }else{
2033 vxs.push_back(myxs+vxs[i-1]);
2034 }
2035 }//for_i_block
2036 //degugging
2037 // 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;}
2038
2039 double totxs = vxs[vxs.size()-1];
2040 int icount=0;
2041 mode_iter:
2042 double pm= EvtRandom::Flat(0.,1);
2043 if(pm <= vxs[0]/totxs) return vmod[0];
2044 for(int j=1;j<vxs.size();j++){
2045 double x0 = vxs[j-1]/totxs;
2046 double x1 = vxs[j]/totxs;
2047 if(x0<pm && pm<=x1) return vmod[j];
2048 }
2049
2050 icount++;
2051 if(icount<10000) goto mode_iter;
2052
2053 //debugging
2054 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;}
2055
2056 std::cout<<"EvtConExc::selectMode can not find a mode with 100 iteration"<<std::endl;
2057 abort();
2058}
2059
2060
2061void EvtConExc::resetResMass(){
2062 // EvtGen base class resetwidth has bugs, it make the mass of particle changed.
2063 EvtId myres = EvtPDL::getId(std::string("J/psi"));
2064 EvtPDL::reSetMass(myres,mjsi);
2065 //EvtPDL::reSetWidth(myres,wjsi);
2066
2067 myres = EvtPDL::getId(std::string("psi(2S)"));
2068 EvtPDL::reSetMass(myres,mpsip);
2069 //EvtPDL::reSetWidth(myres,wpsip);
2070
2071 myres = EvtPDL::getId(std::string("psi(3770)"));
2072 EvtPDL::reSetMass(myres,mpsipp);
2073 //EvtPDL::reSetWidth(myres,wpsipp);
2074
2075 myres = EvtPDL::getId(std::string("phi"));
2076 EvtPDL::reSetMass(myres,mphi);
2077 //EvtPDL::reSetWidth(myres,wphi);
2078
2079 myres = EvtPDL::getId(std::string("omega"));
2080 EvtPDL::reSetMass(myres,momega);
2081 //EvtPDL::reSetWidth(myres,womega);
2082
2083 myres = EvtPDL::getId(std::string("rho0"));
2084 EvtPDL::reSetMass(myres,mrho0);
2085 //EvtPDL::reSetWidth(myres,wrho0);
2086
2087 myres = EvtPDL::getId(std::string("rho(3S)0"));
2088 EvtPDL::reSetMass(myres,mrho3s);
2089 //EvtPDL::reSetWidth(myres,wrho3s);
2090
2091 myres = EvtPDL::getId(std::string("omega(2S)"));
2092 EvtPDL::reSetMass(myres,momega2s);
2093 //EvtPDL::reSetWidth(myres,womega2s);
2094}
2095
2096void EvtConExc::getResMass(){
2097 EvtId myres = EvtPDL::getId(std::string("J/psi"));
2098 mjsi = EvtPDL::getMeanMass(myres);
2099 wjsi = EvtPDL::getWidth(myres);
2100
2101 myres = EvtPDL::getId(std::string("psi(2S)"));
2102 mpsip= EvtPDL::getMeanMass(myres);
2103 wpsip= EvtPDL::getWidth(myres);
2104
2105 myres = EvtPDL::getId(std::string("psi(3770)"));
2106 mpsipp= EvtPDL::getMeanMass(myres);
2107 wpsipp= EvtPDL::getWidth(myres);
2108
2109 myres = EvtPDL::getId(std::string("phi"));
2110 mphi = EvtPDL::getMeanMass(myres);
2111 wphi = EvtPDL::getWidth(myres);
2112
2113 myres = EvtPDL::getId(std::string("omega"));
2114 momega= EvtPDL::getMeanMass(myres);
2115 womega= EvtPDL::getWidth(myres);
2116
2117 myres = EvtPDL::getId(std::string("rho0"));
2118 mrho0 = EvtPDL::getMeanMass(myres);
2119 wrho0 = EvtPDL::getWidth(myres);
2120
2121 myres = EvtPDL::getId(std::string("rho(3S)0"));
2122 mrho3s= EvtPDL::getMeanMass(myres);
2123 wrho3s= EvtPDL::getWidth(myres);
2124
2125
2126 myres = EvtPDL::getId(std::string("omega(2S)"));
2127 momega2s= EvtPDL::getMeanMass(myres);
2128 womega2s= EvtPDL::getWidth(myres);
2129
2130}
2131
2132void EvtConExc::showResMass(){
2133 std::cout<<"J/psi: "<<mjsi<<" "<<wjsi<<std::endl;
2134 std::cout<<"psipr: "<<mpsip<<" "<<wpsip<<std::endl;
2135 std::cout<<"psipp: "<<mpsipp<<" "<<wpsipp<<std::endl;
2136 std::cout<<"phi : "<<mphi<<" "<<wphi<<std::endl;
2137 std::cout<<"omega: "<<momega<<" "<<womega<<std::endl;
2138 std::cout<<"rho0 : "<<mrho0<<" "<<wrho0<<std::endl;
2139 std::cout<<"rho(3S)0: "<<mrho3s<<" "<<wrho3s<<std::endl;
2140 std::cout<<"omega(2S): "<<momega2s<<" "<<womega2s<<std::endl;
2141}
2142
2143bool EvtConExc::checkdecay(EvtParticle* p){
2144 int nson=p->getNDaug();
2145 double massflag=1;
2146 for(int i=0;i<nson;i++){
2147 if( EvtPDL::getStdHep(p->getDaug(i)->getId())==22 ) continue;
2148 massflag *= p->getDaug(i)->mass();
2149 }
2150 if(massflag==0){
2151 std::cout<<"Zero mass!"<<std::endl;
2152 return 0;
2153 }else{return 1;}
2154}
2155
2156
2157double EvtConExc::sumExc(double mx){//summation all cross section of exclusive decays
2158 std::vector<int> vmod; vmod.clear();
2159 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
2160 50,51,
2161 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
2162 90,91,92,93,94,95,96,
2163 74100,74101,74102,74103,74104,74105,74110};
2164 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
2165 50,51,
2166 68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,
2167 90,91,92,93,94,95,96,
2168 74100,74101,74102,74103,74104,74105,74110};
2169 if(_cms > 2.5){
2170 for(int i=0;i<78;i++){vmod.push_back(mn[i]);}
2171 }else{
2172 for(int i=0;i<79;i++){vmod.push_back(mn2[i]);}
2173 }
2174 double xsum=0;
2175 double uscale = 1;
2176 for(int i=0;i<vmod.size();i++){
2177 int mymode = vmod[i];
2178 if(mymode ==74110) continue;
2179 delete myxsection;
2180 myxsection =new EvtXsection (mymode);
2181 init_mode(mymode);
2182 if(myxsection->getUnit() =="pb") {uscale =0.001;}else{uscale=1;}
2183 xsum += uscale*myxsection->getXsection(mx);
2184 //if(mx==3.0967) {std::cout<<vmod[i]<<" "<<uscale<<" "<<xsum<<std::endl;}
2185 }
2186 return xsum;
2187}
2188
2189bool EvtConExc::angularSampling(EvtParticle* par){
2190 bool tagp,tagk;
2191 tagk=0;
2192 tagp=0;
2193 int nds = par->getNDaug();
2194 for(int i=0;i<par->getNDaug();i++){
2195 EvtId di=par->getDaug(i)->getId();
2196 EvtVector4R p4i = par->getDaug(i)->getP4Lab();
2197 int pdgcode =EvtPDL::getStdHep( di );
2198 double alpha=1;
2199 /*
2200 if(fabs(pdgcode)==2212){//proton and anti-proton
2201 alpha = baryonAng(p4i.mass());
2202 cosp = cos(p4i.theta());
2203 tagp=1;
2204 }else if(fabs(pdgcode)==321||fabs(pdgcode)==211 ){
2205 alpha=1;
2206 cosk = cos(p4i.theta());
2207 tagk=1;
2208 }else continue;
2209 */
2210 double angmax = 1+alpha;
2211 double costheta = cos(p4i.theta());
2212 double ang=1+alpha*costheta*costheta;
2213 double xratio = ang/angmax;
2214 double xi=EvtRandom::Flat(0.,1);
2215 //std::cout<<"pdgcode "<<pdgcode<<std::endl;
2216 //std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
2217 if(xi>xratio) return false;
2218 }//loop over duaghters
2219 //if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
2220 //if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
2221
2222 return true;
2223}
2224
2225double EvtConExc::baryonAng(double mx){
2226 double mp=0.938;
2227 double u = 0.938/mx;
2228 u = u*u;
2229 double u2 = (1+u)*(1+u);
2230 double uu = u*(1+6*u);
2231 double alpha = (u2-uu)/(u2+uu);
2232 return alpha;
2233}
2234
2235bool EvtConExc::photonSampling(EvtParticle* par){
2236 bool tagp,tagk;
2237 tagk=0;
2238 tagp=0;
2239 int nds = par->getNDaug();
2240 for(int i=0;i<par->getNDaug();i++){
2241 EvtId di=par->getDaug(i)->getId();
2242 EvtVector4R p4i = par->getDaug(i)->getP4Lab();
2243 int pdgcode =EvtPDL::getStdHep( di );
2244 double alpha=0;
2245
2246 if(pdgcode==111 ||pdgcode==221 || pdgcode==331 ){//for photon
2247 alpha = 1;
2248 }else continue;
2249
2250 double angmax = 1+alpha;
2251 double costheta = cos(p4i.theta());
2252 double ang=1+alpha*costheta*costheta;
2253 double xratio = ang/angmax;
2254 double xi=EvtRandom::Flat(0.,1);
2255 //std::cout<<"pdgcode "<<pdgcode<<std::endl;
2256 //std::cout<<ang<<" "<<angmax<<" "<<xi<<" "<<xratio<<std::endl;
2257 if(xi>xratio) return false;
2258 }//loop over duaghters
2259 //if(tagp) std::cout<<" cosp "<<cosp<<std::endl;
2260 //if(tagk) std::cout<<" cosk "<<cosk<<std::endl;
2261
2262 return true;
2263}
2264
2265
2266double EvtConExc::narrowRXS(double mxL,double mxH){
2267 //br for ee
2268 double psippee,psipee,jsiee,phiee,omegaee,rhoee;
2269 double kev2Gev=0.000001;
2270 psippee=0.262*kev2Gev;
2271 psipee =2.36*kev2Gev;
2272 jsiee =5.55*kev2Gev;
2273 phiee=4.266*0.001*2.954*0.0001;
2274 omegaee=0.6*kev2Gev;
2275 rhoee=7.04*kev2Gev;
2276 double s=_cms*_cms;
2277 double sigv=0;
2278 double mv=0;
2279 double widee=0;
2280 double xpi=12*3.1415926*3.1415926;
2281 if(mxL<=3.0957 && 3.0957<=mxH || mxL<=3.0983 && 3.0983<=mxH) {
2282 widee = jsiee;
2283 mv = 3.096916;
2284 }else if(mxL<=3.6847 && 3.6847<=mxH || mxL<=3.6873 && 3.6873<=mxH){
2285 widee = psipee;
2286 mv = 3.686109;
2287 }else if(mxL<=1.01906 && 1.01906<=mxH || mxL<= 1.01986 && 1.01986<=mxH){
2288 widee = phiee;
2289 mv = 1.01946;
2290 }else{return -1.0;}
2291
2292 if(_cms<=mv) return -1.;
2293 double xv = 1-mv*mv/s;
2294 sigv += xpi*widee/mv/s*Rad2(s,xv);
2295 double unic = 3.89379304*100000; //in unit nb
2296 return sigv*unic; //vaccum factor has included in the Rad2
2297}
2298
2299
2300double EvtConExc::addNarrowRXS(double mhi,double binwidth){
2301 double myxs = 0;
2302 for(int i=0;i<ISRXS.size();i++){ // std::cout<<"ISRM "<<ISRM[i]<<std::endl;
2303 if( (ISRM[i]-mhi)<binwidth && ISRFLAG[i]){ISRFLAG[i]=0; myxs = ISRXS[i];}
2304 }
2305 //std::cout<<mhi<<" "<<binwidth<<" "<<myxs<<std::endl;
2306 return myxs;
2307}
2308
2309double EvtConExc::selectMass(){
2310 double pm,mhdL,mhds;
2311 pm = EvtRandom::Flat(0.,1);
2312 mhdL =_mhdL;
2313 mhds = pm*(_cms - mhdL)+mhdL; //non narrow resonance
2314 std::vector<double> sxs;
2315 for(int i=0;i<ISRID.size();i++){
2316 double ri=ISRXS[i]/AF[598];
2317 sxs.push_back(ri);
2318 }
2319 for(int j=0;j<sxs.size();j++){
2320 if(j>0) sxs[j] += sxs[j-1];
2321 }
2322 pm = EvtRandom::Flat(0.,1);
2323 if(pm<sxs[0]) {
2324 mhds = EvtPDL::getMass(ISRID[0]); //narrow resonance
2325 }
2326 for(int j=1;j<sxs.size();j++){
2327 double x0 = sxs[j-1];
2328 double x1 = sxs[j];
2329 if(x0<pm && pm<=x1) mhds = EvtPDL::getMass(ISRID[j]); //narrow resonance
2330 }
2331 return mhds;
2332}
2333
2334double EvtConExc::getVP(double mx){
2335 double xx,r1,i1;
2336 double x1,y1,y2;
2337 xx=0;
2338 int mytag=1;
2339 for(int i=0;i<4001;i++){
2340 x1=vpx[i];
2341 y1=vpr[i];
2342 y2=vpi[i];
2343 if(xx<=mx && mx <x1){xx=x1; r1=y1; i1=y2; mytag=0; break;}
2344 xx=x1; r1=y1; i1=y2;
2345 }
2346 if(mytag==1){std::cout<<"No vacuum polarization value found"<<std::endl;abort();}
2347 EvtComplex cvp(r1,i1);
2348 double thevp=abs(1./(1-cvp));
2349 if(3.0933<mx && mx<3.1035) return 1.; //J/psi peak excluded
2350 if(3.6810<mx && mx<3.6913) return 1.; //psi(2S) peak excluded
2351 return thevp*thevp;
2352}
2353
2354
2355void EvtConExc::ReadVP(){
2356 vpx.clear();vpr.clear();vpi.clear();
2357 int vpsize=4001;
2358 vpx.resize(vpsize);
2359 vpr.resize(vpsize);
2360 vpi.resize(vpsize);
2361 std::string locvp=getenv("BESEVTGENROOT");
2362 locvp +="/share/vp.dat";
2363 ifstream m_inputFile;
2364 m_inputFile.open(locvp.c_str());
2365 double xx,r1,i1;
2366 double x1,y1,y2;
2367 xx=0;
2368 int mytag=1;
2369 for(int i=0;i<4001;i++){
2370 m_inputFile>>vpx[i]>>vpr[i]>>vpi[i];
2371 }
2372
2373}
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:1782
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:1755
Rad2difXs_er2
double Rad2difXs_er2(double *mhs)
Definition: EvtConExc.cc:1797
dgauss_
double dgauss_(double(*fun)(double *), double *, double *, double *)
Rad2difXs_er
double Rad2difXs_er(double *m)
Definition: EvtConExc.cc:1769
polint_
void polint_(float *, float *, int *, float *, float *)
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition: EvtDiracSpinor.cc:75
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:2143
EvtConExc::Li2
double Li2(double x)
Definition: EvtConExc.cc:1836
EvtConExc::findMaxXS
void findMaxXS(EvtParticle *p)
Definition: EvtConExc.cc:1466
EvtConExc::getMode
static int getMode
Definition: EvtConExc.hh:140
EvtConExc::conexcmode
static int conexcmode
Definition: EvtConExc.hh:157
EvtConExc::addNarrowRXS
double addNarrowRXS(double mhi, double binwidth)
Definition: EvtConExc.cc:2300
EvtConExc::init
void init()
Definition: EvtConExc.cc:187
EvtConExc::narrowRXS
double narrowRXS(double mxL, double mxH)
Definition: EvtConExc.cc:2266
EvtConExc::VP_sampling
bool VP_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:1592
EvtConExc::init_Br_ee
void init_Br_ee()
Definition: EvtConExc.cc:1710
EvtConExc::gamHXSection_er
double gamHXSection_er(double El, double Eh)
Definition: EvtConExc.cc:1448
EvtConExc::islgr
bool islgr(double *x, double *y, int n, double t)
Definition: EvtConExc.cc:1857
EvtConExc::showResMass
void showResMass()
Definition: EvtConExc.cc:2132
EvtConExc::lgr
double lgr(double *x, double *y, int n, double t)
Definition: EvtConExc.cc:1842
EvtConExc::ReadVP
void ReadVP()
Definition: EvtConExc.cc:2355
EvtConExc::baryonAng
double baryonAng(double mx)
Definition: EvtConExc.cc:2225
EvtConExc::Ros_xs
double Ros_xs(double mx, double bree, EvtId pid)
Definition: EvtConExc.cc:1742
EvtConExc::Rad1
double Rad1(double s, double x)
Definition: EvtConExc.cc:1603
EvtConExc::myxsection
static EvtXsection * myxsection
Definition: EvtConExc.hh:136
EvtConExc::SetP4Rvalue
void SetP4Rvalue(EvtParticle *part, double mhdr, double xeng, double theta)
Definition: EvtConExc.cc:1970
EvtConExc::meson_sampling
bool meson_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:1581
EvtConExc::photonSampling
bool photonSampling(EvtParticle *part)
Definition: EvtConExc.cc:2235
EvtConExc::ISR_ang_integrate
double ISR_ang_integrate(double x, double theta)
Definition: EvtConExc.cc:1911
EvtConExc::xs_sampling
bool xs_sampling(double xs)
Definition: EvtConExc.cc:1552
EvtConExc::selectMode
int selectMode(std::vector< int > vmod, double mhds)
Definition: EvtConExc.cc:2014
EvtConExc::gamHXSection
double gamHXSection(EvtParticle *p, double El, double Eh, int nmc=100000)
Definition: EvtConExc.cc:1340
EvtConExc::Rad1difXs
double Rad1difXs(EvtParticle *p)
Definition: EvtConExc.cc:1682
EvtConExc::sumExc
double sumExc(double mx)
Definition: EvtConExc.cc:2157
EvtConExc::init_mode
void init_mode(int mode)
Definition: EvtConExc.cc:576
EvtConExc::baryon_sampling
bool baryon_sampling(EvtVector4R pcm, EvtVector4R pi)
Definition: EvtConExc.cc:1564
EvtConExc::LLr
double LLr(double *x, double *y, int n, double t)
Definition: EvtConExc.cc:1872
EvtConExc::Rad2difXs
double Rad2difXs(EvtParticle *p)
Definition: EvtConExc.cc:1639
EvtConExc::SoftPhoton_xs
double SoftPhoton_xs(double s, double b)
Definition: EvtConExc.cc:1811
EvtConExc::angularSampling
bool angularSampling(EvtParticle *part)
Definition: EvtConExc.cc:2189
EvtConExc::difgamXs
double difgamXs(EvtParticle *p)
Definition: EvtConExc.cc:1511
EvtConExc::ISR_ang_sampling
double ISR_ang_sampling(double x)
Definition: EvtConExc.cc:1929
EvtConExc::initProbMax
void initProbMax()
Definition: EvtConExc.cc:1024
EvtConExc::Mhad_sampling
double Mhad_sampling(double *x, double *y)
Definition: EvtConExc.cc:1888
EvtConExc::_cms
static double _cms
Definition: EvtConExc.hh:137
EvtConExc::resetResMass
void resetResMass()
Definition: EvtConExc.cc:2061
EvtConExc::hadron_angle_sampling
bool hadron_angle_sampling(EvtVector4R ppi, EvtVector4R pcm)
Definition: EvtConExc.cc:1322
EvtConExc::~EvtConExc
virtual ~EvtConExc()
Definition: EvtConExc.cc:165
EvtConExc::getResMass
void getResMass()
Definition: EvtConExc.cc:2096
EvtConExc::gam_sampling
bool gam_sampling(EvtParticle *p)
Definition: EvtConExc.cc:1538
EvtConExc::Rad2
double Rad2(double s, double x)
Definition: EvtConExc.cc:1615
EvtConExc::selectMass
double selectMass()
Definition: EvtConExc.cc:2309
EvtConExc::getVP
double getVP(double cms)
Definition: EvtConExc.cc:2334
EvtConExc::getName
void getName(std::string &name)
Definition: EvtConExc.cc:174
EvtConExc::decay
void decay(EvtParticle *p)
Definition: EvtConExc.cc:1030
EvtConExc::XS_max
static double XS_max
Definition: EvtConExc.hh:138
EvtConExc::clone
EvtDecayBase * clone()
Definition: EvtConExc.cc:180
EvtConExc::SetP4
void SetP4(EvtParticle *part, double mhdr, double xeng, double theta)
Definition: EvtConExc.cc:1953
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
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::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