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