00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include "EvtGenBase/EvtPatches.hh"
00024 #include "EvtGenBase/EvtParticle.hh"
00025 #include "EvtGenBase/EvtStringParticle.hh"
00026 #include "EvtGenBase/EvtDecayTable.hh"
00027 #include "EvtGenBase/EvtPDL.hh"
00028 #include "EvtGenModels/EvtPhokhara.hh"
00029 #include "EvtGenModels/EvtPhokharaDef.hh"
00030 #include "EvtGenBase/EvtReport.hh"
00031 #include <string>
00032 #include "EvtGenBase/EvtId.hh"
00033 #include <iostream>
00034 #include <iomanip>
00035 #include <fstream>
00036 #include <string.h>
00037 #include <stdlib.h>
00038 #include <unistd.h>
00039 #include <stdio.h>
00040
00041 #include "GeneratorObject/McGenEvent.h"
00042 #include "BesKernel/IBesRndmGenSvc.h"
00043 #include "CLHEP/Random/RandomEngine.h"
00044 #include "cfortran/cfortran.h"
00045
00046 using namespace std;
00047 using namespace CLHEP;
00048
00049 using std::endl;
00050 using std::fstream;
00051 using std::ios;
00052 using std::ofstream;
00053 using std::resetiosflags;
00054 using std::setiosflags;
00055 using std::setw;
00056
00057 int EvtPhokhara::nphokharadecays=0;
00058 EvtDecayBasePtr* EvtPhokhara::phokharadecays=0;
00059 int EvtPhokhara::ntable=0;
00060
00061 int EvtPhokhara::ncommand=0;
00062 int EvtPhokhara::lcommand=0;
00063 std::string* EvtPhokhara::commands=0;
00064 int EvtPhokhara::nevt=0;
00065
00066 EvtPhokhara::EvtPhokhara(){}
00067 EvtPhokhara::~EvtPhokhara(){
00068
00069
00070
00071 if(flags_.pion == 9)
00072 maxima_.Mmax[1] = maxima_.Mmax[1] * (1.0 + lambda_par_.alpha_lamb)*(1.0 + lambda_par_.alpha_lamb)
00073 * lambda_par_.ratio_lamb*lambda_par_.ratio_lamb;
00074
00075
00076 if (flags_.nlo == 0) {
00077 sigma = maxima_.Mmax[0]/maxima_.count[0]*maxima_.tr[0];
00078 dsigm = maxima_.Mmax[0]*sqrt((maxima_.tr[0]/maxima_.count[0]-(maxima_.tr[0]/maxima_.count[0])*(maxima_.tr[0]/maxima_.count[0]))/maxima_.count[0]);
00079 } else {
00080 sigma1 = maxima_.Mmax[0]/maxima_.count[0]*maxima_.tr[0];
00081 sigma2 = maxima_.Mmax[1]/maxima_.count[1]*maxima_.tr[1];
00082 dsigm1 = maxima_.Mmax[0]*sqrt((maxima_.tr[0]/maxima_.count[0]-(maxima_.tr[0]/maxima_.count[0])*(maxima_.tr[0]/maxima_.count[0]))/maxima_.count[0]);
00083 dsigm2 = maxima_.Mmax[1]*sqrt((maxima_.tr[1]/maxima_.count[1]-(maxima_.tr[1]/maxima_.count[1])*(maxima_.tr[1]/maxima_.count[1]))/maxima_.count[1]);
00084
00085 sigma = sigma1+sigma2;
00086 dsigm = sqrt(dsigm1*dsigm1 + dsigm2*dsigm2);
00087 }
00088
00089 cout << "-------------------------------------------------------------" << endl;
00090 cout << " PHOKHARA 7.0 Final Statistics " << endl;
00091 cout << "-------------------------------------------------------------" << endl;
00092 cout << int(maxima_.tr[0]+maxima_.tr[1]) << " total events accepted of " << endl;
00093 cout << int(ievent) << " total events generated" << endl;
00094 cout << int(maxima_.tr[0]) << " one photon events accepted of " << endl;
00095 cout << int(maxima_.count[0]) << " events generated" << endl;
00096 cout << int(maxima_.tr[1]) << " two photon events accepted of " << endl;
00097 cout << int(maxima_.count[1]) << " events generated" << endl;
00098 cout << endl;
00099 if (flags_.nlo != 0) cout << "sigma1(nbarn) = " << sigma1 << " +- " << dsigm1 << endl;
00100 if (flags_.nlo != 0) cout << "sigma2(nbarn) = " << sigma2 << " +- " << dsigm2 << endl;
00101 cout << "sigma (nbarn) = " << sigma << " +- " << dsigm << endl;
00102 cout << endl;
00103 cout << "maximum1 = " << maxima_.gross[0] << " minimum1 = " << maxima_.klein[0] << endl;
00104 cout << "Mmax1 = " << maxima_.Mmax[0] << endl;
00105 cout << "maximum2 = " << maxima_.gross[1] << " minimum2 = " << maxima_.klein[1] << endl;
00106 cout << "Mmax2 = " << maxima_.Mmax[1] << endl;
00107 cout << "-------------------------------------------------------------" << endl;
00108
00109
00110 int i;
00111
00112
00113 if (nphokharadecays==0) {
00114 delete [] commands;
00115 commands=0;
00116 return;
00117 }
00118
00119 for(i=0;i<nphokharadecays;i++){
00120 if (phokharadecays[i]==this){
00121 phokharadecays[i]=phokharadecays[nphokharadecays-1];
00122 nphokharadecays--;
00123 if (nphokharadecays==0) {
00124 delete [] commands;
00125 commands=0;
00126 }
00127 return;
00128 }
00129 }
00130
00131 report(ERROR,"EvtGen") << "Error in destroying Phokhara model!"<<endl;
00132
00133 }
00134
00135
00136 void EvtPhokhara::getName(std::string& model_name){
00137
00138 model_name="PHOKHARA";
00139
00140 }
00141
00142 EvtDecayBase* EvtPhokhara::clone(){
00143
00144 return new EvtPhokhara;
00145
00146 }
00147
00148
00149 void EvtPhokhara::initProbMax(){
00150
00151 noProbMax();
00152
00153 }
00154
00155
00156 void EvtPhokhara::init_mode(EvtParticle* p){
00157 m_pion=getArg(0);
00158
00159
00160
00161
00162 if(m_pion<0 || m_pion>9){std::cout<<"mode index for phokhar 0~9, but you give "<<m_pion<<std::endl;abort();}
00163 EvtId myvpho=EvtPDL::getId("vpho");
00164 m_E = p->mass();
00166 m_tagged = 0;
00167 m_nm = 50000 ;
00168 m_nlo = 1;
00169 m_w = 0.0001;
00170 m_fsr = 1;
00171 m_fsrnlo = 1 ;
00172 m_NarrowRes = 1 ;
00173 m_FF_Kaon = 1 ;
00174 m_ivac = 0;
00175 m_FF_Pion = 0 ;
00176 m_f0_model = 0 ;
00177 m_q2min = 0.0;
00178 m_q2_min_c = 0.0447 ;
00179 m_q2_max_c = m_E*m_E;
00180 m_gmin = 0.001;
00181 m_phot1cut = 0.0;
00182 m_phot2cut = 180.0;
00183 m_pi1cut = 0.0 ;
00184 m_pi2cut = 180.0;
00185
00186 if(!(m_pion==0 || m_pion==1 || m_pion==6)){m_fsr = 0; m_fsrnlo = 0 ;}
00187 if( m_pion==9 ){m_nlo = 0 ;}
00188
00189 m_initSeed = 123456789;
00190 RLXDINIT(1, m_initSeed);
00191 maxima_.iprint = 0;
00192 flags_.nlo = m_nlo;
00193 flags_.pion = m_pion;
00194 flags_.fsr = m_fsr;
00195 flags_.fsrnlo = m_fsrnlo;
00196 flags_.ivac = m_ivac;
00197 flags_.FF_pion = m_FF_Pion;
00198 flags_.f0_model = m_f0_model;
00199 flags_.FF_kaon = m_FF_Kaon;
00200 flags_.narr_res = m_NarrowRes;
00201
00202 ctes_.Sp = m_E*m_E; ;
00203
00204 cuts_.w = m_w;
00205 cuts_.q2min = m_q2min;
00206 cuts_.q2_min_c = m_q2_min_c;
00207 cuts_.q2_max_c = m_q2_max_c;
00208 cuts_.gmin = m_gmin;
00209 cuts_.phot1cut = m_phot1cut;
00210 cuts_.phot2cut = m_phot2cut;
00211 cuts_.pi1cut = m_pi1cut;
00212 cuts_.pi2cut = m_pi2cut;
00213
00214 INPUT();
00215
00216
00217 cout << "-------------------------------------------------------------" << endl;
00218 if (flags_.pion == 0)
00219 cout << " PHOKHARA 7.0 : e^+ e^- -> mu^+ mu^- gamma" << endl;
00220 else if (flags_.pion == 1)
00221 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- gamma" << endl;
00222 else if (flags_.pion == 2)
00223 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- 2pi^0 gamma" << endl;
00224 else if (flags_.pion == 3)
00225 cout << " PHOKHARA 7.0: e^+ e^- -> 2pi^+ 2pi^- gamma" << endl;
00226 else if (flags_.pion == 4)
00227 cout << " PHOKHARA 7.0: e^+ e^- -> p pbar gamma" << endl;
00228 else if (flags_.pion == 5)
00229 cout << " PHOKHARA 7.0: e^+ e^- -> n nbar gamma" << endl;
00230 else if (flags_.pion == 6)
00231 cout << " PHOKHARA 7.0: e^+ e^- -> K^+ K^- gamma" << endl;
00232 else if (flags_.pion == 7)
00233 cout << " PHOKHARA 7.0: e^+ e^- -> K_0 K_0bar gamma" << endl;
00234 else if (flags_.pion == 8)
00235 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- pi^0 gamma" << endl;
00236 else if (flags_.pion == 9) {
00237 cout << "PHOKHARA 7.0 : e^+ e^- ->" << endl;
00238 cout << " Lambda (-> pi^- p) Lambda bar (-> pi^+ pbar) gamma" << endl;
00239 } else
00240 cout << " PHOKHARA 7.0: not yet implemented" << endl;
00241
00242
00243 cout << "--------------------------------------------------------------" << endl;
00244 printf(" %s %f %s\n","cms total energy = ",sqrt(ctes_.Sp)," GeV ");
00245
00246
00247 if(cuts_.gmin<0.001){
00248 cerr << " minimal missing energy set too small" << endl;
00249 abort();
00250 }
00251 printf(" %s %f %s\n","minimal tagged photon energy = ",cuts_.gmin," GeV ");
00252 printf(" %s %f,%f\n","angular cuts on tagged photon = ",cuts_.phot1cut,cuts_.phot2cut);
00253
00254
00255 if (flags_.pion == 0)
00256 printf(" %s %f,%f\n","angular cuts on muons = ",cuts_.pi1cut,cuts_.pi2cut);
00257 else if (flags_.pion == 4)
00258 printf(" %s %f,%f\n","angular cuts on protons = ",cuts_.pi1cut,cuts_.pi2cut);
00259 else if (flags_.pion == 5)
00260 printf(" %s %f,%f\n","angular cuts on neutrons = ", cuts_.pi1cut,cuts_.pi2cut);
00261 else if ((flags_.pion == 6)||(flags_.pion == 7))
00262 printf(" %s %f,%f\n","angular cuts on kaons = ", cuts_.pi1cut,cuts_.pi2cut);
00263 else if (flags_.pion == 9)
00264 printf(" %s %f,%f\n","angular cuts on pions and protons = ", cuts_.pi1cut,cuts_.pi2cut);
00265 else
00266 printf(" %s %f,%f\n","angular cuts on pions = ", cuts_.pi1cut,cuts_.pi2cut);
00267 if (flags_.pion == 0)
00268 printf(" %s %f %s\n","min. muons-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2");
00269 else if (flags_.pion == 4)
00270 printf(" %s %f %s\n","min. protons-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2" );
00271 else if (flags_.pion == 5)
00272 printf(" %s %f %s\n","min. neutrons-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2" );
00273 else if ((flags_.pion == 6)||(flags_.pion == 7))
00274 printf(" %s %f %s\n","min. kaons-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2" );
00275 else if (flags_.pion == 9)
00276 printf(" %s %f %s\n"," min. lambdas-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2" );
00277 else
00278 printf(" %s %f %s\n","min. pions-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2" );
00279 cos1min = cos(cuts_.phot2cut*ctes_.pi/180.0);
00280 cos1max = cos(cuts_.phot1cut*ctes_.pi/180.0);
00281 cos2min = -1.0;
00282 cos2max = 1.0;
00283 cos3min = -1.0;
00284 cos3max = 1.0;
00285 if (flags_.pion == 0)
00286 qqmin = 4.0*ctes_.mmu*ctes_.mmu;
00287 else if (flags_.pion == 1)
00288 qqmin = 4.0*ctes_.mpi*ctes_.mpi;
00289 else if (flags_.pion == 2)
00290 qqmin = 4.0*(ctes_.mpi+ctes_.mpi0)*(ctes_.mpi+ctes_.mpi0);
00291 else if (flags_.pion == 3)
00292 qqmin = 16.0*ctes_.mpi*ctes_.mpi;
00293 else if (flags_.pion == 4)
00294 qqmin = 4.0*ctes_.mp*ctes_.mp;
00295 else if (flags_.pion == 5)
00296 qqmin = 4.0*ctes_.mnt*ctes_.mnt;
00297 else if (flags_.pion == 6)
00298 qqmin = 4.0*ctes_.mKp*ctes_.mKp;
00299 else if (flags_.pion == 7)
00300 qqmin = 4.0*ctes_.mKn*ctes_.mKn;
00301 else if (flags_.pion == 8)
00302 qqmin = (2.0*ctes_.mpi+ctes_.mpi0)*(2.0*ctes_.mpi+ctes_.mpi0);
00303 else if (flags_.pion == 9)
00304 qqmin = 4.0*ctes_.mlamb*ctes_.mlamb;
00305 qqmax = ctes_.Sp-2.0*sqrt(ctes_.Sp)*cuts_.gmin;
00306 if (cuts_.q2_max_c < qqmax)
00307 qqmax=cuts_.q2_max_c;
00308
00309
00310 if ( (cuts_.q2_min_c > qqmin) && (cuts_.q2_min_c < (ctes_.Sp*(1.0-2.0*(cuts_.gmin/sqrt(ctes_.Sp)+cuts_.w))) ) )
00311 qqmin = cuts_.q2_min_c;
00312 else {
00313 cerr << "------------------------------" << endl;
00314 cerr << " Q^2_min TOO SMALL" << endl;
00315 cerr << " Q^2_min CHANGED BY PHOKHARA = " << qqmin << " GeV^2" << endl;
00316 cerr << "------------------------------" << endl;
00317 }
00318
00319 if(qqmax <= qqmin){
00320 cerr << " Q^2_max to small " << endl;
00321 cerr << " Q^2_max = " << qqmax << endl;
00322 cerr << " Q^2_min = " << qqmin << endl;
00323 abort();
00324 }
00325
00326
00327 if (flags_.pion == 0) {
00328 printf(" %s %f %s\n", "minimal muon-pair invariant mass^2 = ", qqmin," GeV^2");
00329 printf(" %s %f %s\n", "maximal muon-pair invariant mass^2 = ", qqmax," GeV^2");
00330 } else if (flags_.pion == 1) {
00331 printf(" %s %f %s\n", "minimal pion-pair invariant mass^2 = ", qqmin," GeV^2");
00332 printf(" %s %f %s\n", "maximal pion-pair invariant mass^2 = ", qqmax," GeV^2");
00333 } else if (flags_.pion == 4) {
00334 printf(" %s %f %s\n", "minimal proton-pair invariant mass^2 = ", qqmin," GeV^2");
00335 printf(" %s %f %s\n", "maximal proton-pair invariant mass^2 = ", qqmax," GeV^2");
00336 } else if (flags_.pion == 5) {
00337 printf(" %s %f %s\n", "minimal neutron-pair invariant mass^2 = ", qqmin," GeV^2");
00338 printf(" %s %f %s\n", "maximal neutron-pair invariant mass^2 = ", qqmax," GeV^2");
00339 } else if ((flags_.pion == 6) || (flags_.pion == 7)) {
00340 printf(" %s %f %s\n", "minimal kaon-pair invariant mass^2 = ", qqmin," GeV^2");
00341 printf(" %s %f %s\n", "maximal kaon-pair invariant mass^2 = ", qqmax," GeV^2");
00342 } else if(flags_.pion == 8){
00343 printf(" %s %f %s\n", "minimal three-pion invariant mass^2 = ", qqmin," GeV^2");
00344 printf(" %s %f %s\n", "maximal three-pion invariant mass^2 = ", qqmax," GeV^2");
00345 } else if(flags_.pion == 9){
00346 printf(" %s %f %s\n", "minimal lambda-pair invariant mass^2 = ", qqmin," GeV^2");
00347 printf(" %s %f %s\n", "maximal lambda-pair invariant mass^2 = ", qqmax," GeV^2");
00348 } else {
00349 printf(" %s %f %s\n", "minimal four-pion invariant mass^2 = ", qqmin," GeV^2" );
00350 printf(" %s %f %s\n", "maximal four-pion invariant mass^2 = ", qqmax," GeV^2");
00351 }
00352
00353 if (flags_.nlo == 0) {
00354 cout << "Born" << endl;
00355 if(flags_.fsrnlo != 0){
00356 cerr << "WRONG FSRNLO flag - only fsrnlo=0 allowed for Born" << endl;
00357 abort();
00358 }
00359 }
00360
00361 if((flags_.pion == 9) && (flags_.nlo != 0)) {
00362 cerr << "WRONG NLO flag - only Born allowed for Lambdas"<< endl;
00363 cerr << "If you feel that you need NLO"<< endl;
00364 cerr << "please contact the authors"<< endl;
00365 abort();
00366 }
00367
00368 if (flags_.nlo == 1) printf(" %s %f\n", "NLO: soft photon cutoff w = ",cuts_.w);
00369 if ((flags_.pion <= 1) || (flags_.pion == 6))
00370 {
00371
00372 if( ! ((flags_.fsr == 1) || (flags_.fsr == 2) || (flags_.fsrnlo == 0)
00373 || (flags_.fsr == 1) || (flags_.fsrnlo == 1)
00374 || (flags_.fsr == 0) || (flags_.fsrnlo == 0))) {
00375 cerr << "WRONG combination of FSR, FSRNLO flags" <<endl;
00376 abort();
00377 }
00378
00379 if (flags_.fsr == 0)
00380 cout << "ISR only" << endl;
00381 else if (flags_.fsr == 1)
00382 cout << "ISR+FSR" << endl;
00383 else if (flags_.fsr == 2) {
00384 if (flags_.nlo == 0)
00385 cout << "ISR+INT+FSR" << endl;
00386 else {
00387 cerr << "WRONG FSR flag: interference is included only for nlo=0" << endl;
00388 abort();
00389 }
00390 }
00391 else {
00392 cerr << "WRONG FSR flag" << flags_.fsr << endl;
00393 abort();
00394 }
00395 if(flags_.fsrnlo == 1)
00396 cout << "IFSNLO included" << endl;
00397 }
00398 else
00399 {
00400 if((flags_.fsr == 0) && (flags_.fsrnlo == 0))
00401 cout << "ISR only" << endl;
00402 else {
00403 cerr << "FSR is implemented only for pi+pi-, mu+mu- and K+K- modes" << endl;
00404 abort();
00405 }
00406 }
00407
00408 if(flags_.ivac == 0){
00409 cout << "Vacuum polarization is NOT included" << endl;}
00410 else if(flags_.ivac == 1){
00411 cout << "Vacuum polarization by Fred Jegerlehner (http://www-com.physik.hu-berlin.de/fjeger/alphaQEDn.uu)" << endl;}
00412 else if(flags_.ivac == 2){
00413 cout << "Vacuum polarization (VP_HLMNT_v1_3nonr) by Daisuke Nomura and Thomas Teubner" << endl;}
00414 else {
00415 cout << "WRONG vacuum polarization switch" << endl;
00416 abort();
00417 }
00418
00419
00420 if(flags_.pion == 1){
00421 if(flags_.FF_pion == 0)
00422 cout << "Kuhn-Santamaria PionFormFactor" << endl;
00423 else if(flags_.FF_pion == 1)
00424 cout << "Gounaris-Sakurai PionFormFactor old" << endl;
00425 else if(flags_.FF_pion == 2)
00426 cout << "Gounaris-Sakurai PionFormFactor new" << endl;
00427 else {
00428 cout << "WRONG PionFormFactor switch" << endl;
00429 abort();
00430 }
00431 if(flags_.fsr != 0){
00432 if(flags_.f0_model == 0)
00433 cout << "f0+f0(600): K+K- model" << endl;
00434 else if(flags_.f0_model == 1)
00435 cout << "f0+f0(600): \"no structure\" model" << endl;
00436 else if(flags_.f0_model == 2)
00437 cout << "NO f0+f0(600)" << endl;
00438 else if(flags_.f0_model == 3)
00439 cout << "only f0, KLOE: Cesare Bini-private communication" << endl;
00440 else {
00441 cout << "WRONG f0+f0(600) switch" << endl;
00442 abort();
00443 }
00444 }
00445 }
00446
00447 if((flags_.pion == 6) || (flags_.pion==7)){
00448 if(flags_.FF_kaon == 0)
00449 cout << "constrained KaonFormFactor" << endl;
00450 else if(flags_.FF_kaon == 1) {
00451 cout << "unconstrained KaonFormFactor" << endl;}
00452 else if(flags_.FF_kaon == 2) {
00453 cout << "Kuhn-Khodjamirian-Bruch KaonFormFactor" << endl;}
00454 else{
00455 cout << "WRONG KaonFormFactor switch" << endl;
00456 abort();
00457 }
00458 }
00459
00460 if((flags_.pion == 0) || (flags_.pion==1) || (flags_.pion == 6) || (flags_.pion == 7)){
00461 if(flags_.narr_res == 1){
00462 cout << "THE NARROW RESONANCE J/PSI INCLUDED" << endl;}
00463 else if(flags_.narr_res == 2){
00464 cout << "THE NARROW RESONANCE PSI(2S) INCLUDED" << endl;}
00465 else if(flags_.narr_res != 0){
00466 cout << "wrong flag narr_res: only 0, 1 or 2 allowed" << endl;
00467 abort();
00468 }
00469 }
00470
00471 cout << "--------------------------------------------------------------" << endl;
00472
00473
00474
00475 for(int i = 0; i<2; i++)
00476 {
00477 maxima_.Mmax[i] = 1.0;
00478 maxima_.gross[i] = 0.0;
00479 maxima_.klein[i] = 0.0;
00480 }
00481
00482 if (flags_.nlo == 0)
00483 maxima_.Mmax[1]=0.0;
00484
00485 maxima_.tr[0] = 0.0;
00486 maxima_.tr[1] = 0.0;
00487 maxima_.count[0] = 0.0;
00488 maxima_.count[1] = 0.0;
00489
00490
00491
00492 for(int j = 1; j <= m_nm; j++)
00493 {
00494 RANLXDF(Ar_r,1);
00495 Ar[1] = Ar_r[0];
00496 if (Ar[1] <= (maxima_.Mmax[0]/(maxima_.Mmax[0]+maxima_.Mmax[1]))) {
00497 maxima_.count[0] = maxima_.count[0]+1.0;
00498 GEN_1PH(1,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
00499 }
00500 else {
00501 maxima_.count[1] = maxima_.count[1]+1.0;
00502 GEN_2PH(1,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
00503 }
00504 }
00505
00506
00507
00508 maxima_.Mmax[0] = maxima_.gross[0]+.05*sqrt(maxima_.gross[0]*maxima_.gross[0]);
00509 maxima_.Mmax[1] = maxima_.gross[1]+(.03+.02*ctes_.Sp)*sqrt(maxima_.gross[1]*maxima_.gross[1]);
00510 theMmax[m_pion][0]=maxima_.Mmax[0];
00511 theMmax[m_pion][1]=maxima_.Mmax[1];
00512
00513 if((flags_.pion == 1) && (flags_.fsrnlo == 1))
00514 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
00515 if((flags_.pion == 0) && (flags_.fsrnlo == 1))
00516 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
00517
00518 if((flags_.pion == 0) && (flags_.fsr == 1) && (flags_.fsrnlo == 0))
00519 maxima_.Mmax[1]=maxima_.Mmax[1]*1.2;
00520
00521 if((flags_.pion == 2) || (flags_.pion == 3)){
00522 maxima_.Mmax[0]=maxima_.Mmax[0]*1.1;
00523 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
00524 }
00525
00526 if(flags_.pion == 8){
00527 maxima_.Mmax[0]=maxima_.Mmax[0]*1.08;
00528 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
00529 }
00530
00531
00532 maxima_.tr[0] = 0.0;
00533 maxima_.tr[1] = 0.0;
00534 maxima_.count[0] = 0.0;
00535 maxima_.count[1] = 0.0;
00536 }
00537
00538
00539
00540 void EvtPhokhara::init(){
00541 checkNDaug(0);
00542 checkNArg(1);
00543 nevtgen.resize(10);
00544 theMmax.resize(10);
00545 for(int i=0;i<=10;i++){ theMmax[i].resize(2); nevtgen[i]=0;}
00546
00547 Vfs.push_back(" mu+mu-");
00548 Vfs.push_back(" pi+pi-");
00549 Vfs.push_back(" 2pi0pi+pi-");
00550 Vfs.push_back(" 2pi+2pi-");
00551 Vfs.push_back(" ppbar");
00552 Vfs.push_back(" nnbar");
00553 Vfs.push_back(" K+K-");
00554 Vfs.push_back(" K0K0bar");
00555 Vfs.push_back(" pi+pi-pi0");
00556 Vfs.push_back(" Lambda Lambdabar");
00557
00558 std::string locvp=getenv("BESEVTGENROOT");
00559 system("cat $BESEVTGENROOT/share/phokhara.param>phokhara.param");
00560
00561
00562 if (getParentId().isAlias()){
00563
00564 report(ERROR,"EvtGen") << "EvtPhokhara finds that you are decaying the"<<endl
00565 << " aliased particle "
00566 << EvtPDL::name(getParentId()).c_str()
00567 << " with the Phokhara model"<<endl
00568 << " this does not work, please modify decay table."
00569 << endl;
00570 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00571 ::abort();
00572
00573 }
00574
00575
00576
00577 }
00578
00579
00580 std::string EvtPhokhara::commandName(){
00581
00582 return std::string("PhokharaPar");
00583
00584 }
00585
00586 void EvtPhokhara::command(std::string cmd){
00587
00588 if (ncommand==lcommand){
00589
00590 lcommand=10+2*lcommand;
00591
00592 std::string* newcommands=new std::string[lcommand];
00593
00594 int i;
00595
00596 for(i=0;i<ncommand;i++){
00597 newcommands[i]=commands[i];
00598 }
00599
00600 delete [] commands;
00601
00602 commands=newcommands;
00603
00604 }
00605
00606 commands[ncommand]=cmd;
00607
00608 ncommand++;
00609
00610 }
00611
00612
00613
00614 void EvtPhokhara::decay( EvtParticle *p){
00615 EvtId myvpho=EvtPDL::getId("vpho");
00616 if(p->getId()!=myvpho) {std::cout<<"Parent particle is required to be vpho for Phokhara model"<<std::endl;abort();}
00617 m_pion=getArg(0);
00618 if(nevtgen[m_pion]==0) {init_mode(p);}
00619 else{init_evt(p);}
00620 std::cout<<"PHOKHARA : "<<Vfs[m_pion]<<" mode "<<std::endl;
00621
00622 int istdheppar=EvtPDL::getStdHep(p->getId());
00623 int ntrials = 0;
00624 int tr_old[2];
00625 tr_old[0] = (int)maxima_.tr[0];
00626 tr_old[1] = (int)maxima_.tr[1];
00627
00628 while( ntrials < 10000)
00629 {
00630 ievent++;
00631 RANLXDF(Ar_r,1);
00632 Ar[1] = Ar_r[0];
00633
00634 if (Ar[1] <= (maxima_.Mmax[0]/(maxima_.Mmax[0]+maxima_.Mmax[1]))) {
00635 maxima_.count[0] = maxima_.count[0]+1.0;
00636 GEN_1PH(2,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
00637 }
00638 else {
00639 maxima_.count[1] = maxima_.count[1]+1.0;
00640 GEN_2PH(2,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
00641 }
00642
00643 if( ((int)maxima_.tr[0]+(int)maxima_.tr[1]) > (tr_old[0]+tr_old[1]) )
00644 {
00645 goto storedEvents;
00646 }
00647 ntrials ++;
00648 }
00649
00650 std::cout <<"FATAL: Could not satisfy cuts after " << ntrials << "trials. Terminate." <<std::endl;
00651
00652 storedEvents:
00653 int more=0;
00654 int numstable=0;
00655 int numparton=0;
00656 EvtId evtnumstable[100];
00657 EvtVector4R p4[20];
00658
00659
00660 if (flags_.pion == 0) {
00661
00662 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-13);
00663 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00664 numstable++;
00665
00666 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(13);
00667 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00668 numstable++;
00669 }
00670 if (flags_.pion == 1) {
00671
00672 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
00673 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00674 numstable++;
00675
00676 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
00677 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00678 numstable++;
00679 }
00680 if (flags_.pion == 2) {
00681
00682 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
00683 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00684 numstable++;
00685
00686 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
00687 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00688 numstable++;
00689
00690 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
00691 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
00692 numstable++;
00693
00694 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
00695 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
00696 numstable++;
00697 }
00698 if (flags_.pion == 3) {
00699
00700 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
00701 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00702 numstable++;
00703
00704 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
00705 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00706 numstable++;
00707
00708 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
00709 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
00710 numstable++;
00711
00712 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
00713 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
00714 numstable++;
00715 }
00716 if (flags_.pion == 4) {
00717
00718 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2212);
00719 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00720 numstable++;
00721
00722 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2212);
00723 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00724 numstable++;
00725 }
00726 if (flags_.pion == 5) {
00727
00728 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2112);
00729 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00730 numstable++;
00731
00732 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2112);
00733 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00734 numstable++;
00735 }
00736 if (flags_.pion == 6) {
00737
00738 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(321);
00739 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00740 numstable++;
00741
00742 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-321);
00743 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00744 numstable++;
00745 }
00746 if (flags_.pion == 7) {
00747
00748 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(310);
00749 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00750 numstable++;
00751
00752 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(130);
00753 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00754 numstable++;
00755 }
00756 if (flags_.pion == 8) {
00757
00758 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
00759 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00760 numstable++;
00761
00762 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
00763 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00764 numstable++;
00765
00766 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
00767 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
00768 numstable++;
00769 }
00770 if (flags_.pion == 9) {
00771
00772 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
00773 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
00774 numstable++;
00775
00776 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2212);
00777 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
00778 numstable++;
00779
00780 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
00781 p4[numstable].set(ctes_.momenta[0][9],ctes_.momenta[1][9], ctes_.momenta[2][9], ctes_.momenta[3][9]);
00782 numstable++;
00783
00784 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2212);
00785 p4[numstable].set(ctes_.momenta[0][10],ctes_.momenta[1][10], ctes_.momenta[2][10], ctes_.momenta[3][10]);
00786 numstable++;
00787 }
00788
00789
00790 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
00791 p4[numstable].set(ctes_.momenta[0][2],ctes_.momenta[1][2], ctes_.momenta[2][2], ctes_.momenta[3][2]);
00792 numstable++;
00793 if( ctes_.momenta[0][3] != 0 )
00794 {
00795 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
00796 p4[numstable].set(ctes_.momenta[0][3],ctes_.momenta[1][3], ctes_.momenta[2][3], ctes_.momenta[3][3]);
00797 numstable++;
00798 }
00799
00800 int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
00801 more=(channel!=-1);
00802 if(more) {std::cout<<"Existence of mode "<<channel<<" in exclusive decay list has the same final state as this one"<<std::endl;abort(); }
00803
00804 p->makeDaughters(numstable,evtnumstable);
00805
00806 int ndaugFound=0;
00807 EvtVector4R SUMP4(0,0,0,0);
00808 for(int i=0;i<numstable;i++){
00809 p->getDaug(i)->init(evtnumstable[i],p4[i]);
00810 ndaugFound++;
00811 }
00812 if ( ndaugFound == 0 ) {
00813 report(ERROR,"EvtGen") << "Phokhara has failed to do a decay ";
00814 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass()<<endl;
00815 assert(0);
00816 }
00817
00818 nevtgen[m_pion]++;
00819 return ;
00820
00821 }
00822
00823
00824
00825 void EvtPhokhara::store(EvtDecayBase* jsdecay){
00826
00827 if (nphokharadecays==ntable){
00828
00829 EvtDecayBasePtr* newphokharadecays=new EvtDecayBasePtr[2*ntable+10];
00830 int i;
00831 for(i=0;i<ntable;i++){
00832 newphokharadecays[i]=phokharadecays[i];
00833 }
00834 ntable=2*ntable+10;
00835 delete [] phokharadecays;
00836 phokharadecays=newphokharadecays;
00837 }
00838
00839 phokharadecays[nphokharadecays++]=jsdecay;
00840
00841
00842
00843 }
00844
00845
00846 void EvtPhokhara::PhokharaInit(int dummy){
00847 static int first=1;
00848 if (first){
00849
00850 first=0;
00851
00852
00853 }
00854
00855 }
00856
00857
00858
00859 void EvtPhokhara::init_evt(EvtParticle* p){
00860 m_pion=getArg(0);
00861
00862
00863
00864
00865 if(m_pion<0 || m_pion>9){std::cout<<"mode index for phokhar 0~9, but you give "<<m_pion<<std::endl;abort();}
00866 EvtId myvpho=EvtPDL::getId("vpho");
00867 m_E = p->mass();
00869 m_tagged = 0;
00870 m_nm = 50000 ;
00871 m_nlo = 1;
00872 m_w = 0.0001;
00873 m_fsr = 1;
00874 m_fsrnlo = 1 ;
00875 m_NarrowRes = 1 ;
00876 m_FF_Kaon = 1 ;
00877 m_ivac = 0;
00878 m_FF_Pion = 0 ;
00879 m_f0_model = 0 ;
00880 m_q2min = 0.0;
00881 m_q2_min_c = 0.0447 ;
00882 m_q2_max_c = m_E*m_E;
00883 m_gmin = 0.001;
00884 m_phot1cut = 0.0;
00885 m_phot2cut = 180.0;
00886 m_pi1cut = 0.0 ;
00887 m_pi2cut = 180.0;
00888
00889 if(!(m_pion==0 || m_pion==1 || m_pion==6)){m_fsr = 0; m_fsrnlo = 0 ;}
00890 if( m_pion==9 ){m_nlo = 0 ;}
00891
00892 maxima_.iprint = 0;
00893 flags_.nlo = m_nlo;
00894 flags_.pion = m_pion;
00895 flags_.fsr = m_fsr;
00896 flags_.fsrnlo = m_fsrnlo;
00897 flags_.ivac = m_ivac;
00898 flags_.FF_pion = m_FF_Pion;
00899 flags_.f0_model = m_f0_model;
00900 flags_.FF_kaon = m_FF_Kaon;
00901 flags_.narr_res = m_NarrowRes;
00902
00903 ctes_.Sp = m_E*m_E; ;
00904
00905 cuts_.w = m_w;
00906 cuts_.q2min = m_q2min;
00907 cuts_.q2_min_c = m_q2_min_c;
00908 cuts_.q2_max_c = m_q2_max_c;
00909 cuts_.gmin = m_gmin;
00910 cuts_.phot1cut = m_phot1cut;
00911 cuts_.phot2cut = m_phot2cut;
00912 cuts_.pi1cut = m_pi1cut;
00913 cuts_.pi2cut = m_pi2cut;
00914
00915 INPUT();
00916
00917 cos1min = cos(cuts_.phot2cut*ctes_.pi/180.0);
00918 cos1max = cos(cuts_.phot1cut*ctes_.pi/180.0);
00919 cos2min = -1.0;
00920 cos2max = 1.0;
00921 cos3min = -1.0;
00922 cos3max = 1.0;
00923 if (flags_.pion == 0)
00924 qqmin = 4.0*ctes_.mmu*ctes_.mmu;
00925 else if (flags_.pion == 1)
00926 qqmin = 4.0*ctes_.mpi*ctes_.mpi;
00927 else if (flags_.pion == 2)
00928 qqmin = 4.0*(ctes_.mpi+ctes_.mpi0)*(ctes_.mpi+ctes_.mpi0);
00929 else if (flags_.pion == 3)
00930 qqmin = 16.0*ctes_.mpi*ctes_.mpi;
00931 else if (flags_.pion == 4)
00932 qqmin = 4.0*ctes_.mp*ctes_.mp;
00933 else if (flags_.pion == 5)
00934 qqmin = 4.0*ctes_.mnt*ctes_.mnt;
00935 else if (flags_.pion == 6)
00936 qqmin = 4.0*ctes_.mKp*ctes_.mKp;
00937 else if (flags_.pion == 7)
00938 qqmin = 4.0*ctes_.mKn*ctes_.mKn;
00939 else if (flags_.pion == 8)
00940 qqmin = (2.0*ctes_.mpi+ctes_.mpi0)*(2.0*ctes_.mpi+ctes_.mpi0);
00941 else if (flags_.pion == 9)
00942 qqmin = 4.0*ctes_.mlamb*ctes_.mlamb;
00943 qqmax = ctes_.Sp-2.0*sqrt(ctes_.Sp)*cuts_.gmin;
00944 if (cuts_.q2_max_c < qqmax)
00945 qqmax=cuts_.q2_max_c;
00946
00947
00948 if ( (cuts_.q2_min_c > qqmin) && (cuts_.q2_min_c < (ctes_.Sp*(1.0-2.0*(cuts_.gmin/sqrt(ctes_.Sp)+cuts_.w))) ) )
00949 qqmin = cuts_.q2_min_c;
00950 else {
00951 }
00952
00953
00954
00955
00956 for(int i = 0; i<2; i++)
00957 {
00958 maxima_.Mmax[i] = 1.0;
00959 maxima_.gross[i] = 0.0;
00960 maxima_.klein[i] = 0.0;
00961 }
00962
00963 if (flags_.nlo == 0)
00964 maxima_.Mmax[1]=0.0;
00965
00966 maxima_.tr[0] = 0.0;
00967 maxima_.tr[1] = 0.0;
00968 maxima_.count[0] = 0.0;
00969 maxima_.count[1] = 0.0;
00970
00971
00972
00973 maxima_.Mmax[0] = theMmax[m_pion][0];
00974 maxima_.Mmax[1] = theMmax[m_pion][1];
00975 if((flags_.pion == 1) && (flags_.fsrnlo == 1))
00976 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
00977 if((flags_.pion == 0) && (flags_.fsrnlo == 1))
00978 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
00979
00980 if((flags_.pion == 0) && (flags_.fsr == 1) && (flags_.fsrnlo == 0))
00981 maxima_.Mmax[1]=maxima_.Mmax[1]*1.2;
00982
00983 if((flags_.pion == 2) || (flags_.pion == 3)){
00984 maxima_.Mmax[0]=maxima_.Mmax[0]*1.1;
00985 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
00986 }
00987
00988 if(flags_.pion == 8){
00989 maxima_.Mmax[0]=maxima_.Mmax[0]*1.08;
00990 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
00991 }
00992
00993
00994 maxima_.tr[0] = 0.0;
00995 maxima_.tr[1] = 0.0;
00996 maxima_.count[0] = 0.0;
00997 maxima_.count[1] = 0.0;
00998 }
00999