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_pi0pi0pipi.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_pi0pi0pipi::nevtgen=0;
00058 int EvtPhokhara_pi0pi0pipi::nphokharadecays=0;
00059 EvtDecayBasePtr* EvtPhokhara_pi0pi0pipi::phokharadecays=0;
00060 int EvtPhokhara_pi0pi0pipi::ntable=0;
00061
00062 int EvtPhokhara_pi0pi0pipi::ncommand=0;
00063 int EvtPhokhara_pi0pi0pipi::lcommand=0;
00064 std::string* EvtPhokhara_pi0pi0pipi::commands=0;
00065 int EvtPhokhara_pi0pi0pipi::nevt=0;
00066
00067 EvtPhokhara_pi0pi0pipi::EvtPhokhara_pi0pi0pipi(){}
00068 EvtPhokhara_pi0pi0pipi::~EvtPhokhara_pi0pi0pipi(){
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_pi0pi0pipi::getName(std::string& model_name){
00096
00097 model_name="PHOKHARA_pi0pi0pipi";
00098
00099 }
00100
00101 EvtDecayBase* EvtPhokhara_pi0pi0pipi::clone(){
00102
00103 return new EvtPhokhara_pi0pi0pipi;
00104
00105 }
00106
00107
00108 void EvtPhokhara_pi0pi0pipi::initProbMax(){
00109
00110 noProbMax();
00111
00112 }
00113
00114
00115 void EvtPhokhara_pi0pi0pipi::init_mode(EvtParticle* p){
00116 m_pion=2;
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 = 12345678;
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_pi0pi0pipi::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_pi0pi0pipi::commandName(){
00524
00525 return std::string("PhokharaPar");
00526
00527 }
00528
00529 void EvtPhokhara_pi0pi0pipi::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_pi0pi0pipi::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
00563 std::cout<<"PHOKHARA : pi0pi0pi+pi- mode "<<std::endl;
00564 int istdheppar=EvtPDL::getStdHep(p->getId());
00565 int ntrials = 0;
00566 int tr_old[2];
00567 tr_old[0] = (int)maxima_.tr[0];
00568 tr_old[1] = (int)maxima_.tr[1];
00569
00570 while( ntrials < 10000)
00571 {
00572 ievent++;
00573 RANLXDF(Ar_r,1);
00574 Ar[1] = Ar_r[0];
00575
00576 if (Ar[1] <= (maxima_.Mmax[0]/(maxima_.Mmax[0]+maxima_.Mmax[1]))) {
00577 maxima_.count[0] = maxima_.count[0]+1.0;
00578 GEN_1PH(2,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
00579 }
00580 else {
00581 maxima_.count[1] = maxima_.count[1]+1.0;
00582 GEN_2PH(2,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
00583 }
00584
00585 if( ((int)maxima_.tr[0]+(int)maxima_.tr[1]) > (tr_old[0]+tr_old[1]) )
00586 {
00587 goto storedEvents;
00588 }
00589 ntrials ++;
00590 }
00591
00592 std::cout <<"FATAL: Could not satisfy cuts after " << ntrials << "trials. Terminate." <<std::endl;
00593
00594 storedEvents:
00595 int more=0;
00596 int numstable=0;
00597 int numparton=0;
00598 EvtId evtnumstable[100];
00599 EvtVector4R p4[20];
00600
00601
00602 if (flags_.pion == 0) {
00603
00604 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-13);
00605 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00606 numstable++;
00607
00608 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(13);
00609 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00610 numstable++;
00611 }
00612 if (flags_.pion == 1) {
00613
00614 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
00615 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00616 numstable++;
00617
00618 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
00619 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00620 numstable++;
00621 }
00622 if (flags_.pion == 2) {
00623
00624 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
00625 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00626 numstable++;
00627
00628 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
00629 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00630 numstable++;
00631
00632 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
00633 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
00634 numstable++;
00635
00636 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
00637 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
00638 numstable++;
00639 }
00640 if (flags_.pion == 3) {
00641
00642 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
00643 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00644 numstable++;
00645
00646 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
00647 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00648 numstable++;
00649
00650 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
00651 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
00652 numstable++;
00653
00654 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
00655 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
00656 numstable++;
00657 }
00658 if (flags_.pion == 4) {
00659
00660 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2212);
00661 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00662 numstable++;
00663
00664 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2212);
00665 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00666 numstable++;
00667 }
00668 if (flags_.pion == 5) {
00669
00670 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2112);
00671 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00672 numstable++;
00673
00674 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2112);
00675 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00676 numstable++;
00677 }
00678 if (flags_.pion == 6) {
00679
00680 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(321);
00681 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00682 numstable++;
00683
00684 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-321);
00685 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00686 numstable++;
00687 }
00688 if (flags_.pion == 7) {
00689
00690 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(311);
00691 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
00692 numstable++;
00693
00694 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-311);
00695 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
00696 numstable++;
00697 }
00698 if (flags_.pion == 8) {
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(111);
00709 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
00710 numstable++;
00711 }
00712 if (flags_.pion == 9) {
00713
00714 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
00715 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
00716 numstable++;
00717
00718 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2212);
00719 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
00720 numstable++;
00721
00722 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
00723 p4[numstable].set(ctes_.momenta[0][9],ctes_.momenta[1][9], ctes_.momenta[2][9], ctes_.momenta[3][9]);
00724 numstable++;
00725
00726 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2212);
00727 p4[numstable].set(ctes_.momenta[0][10],ctes_.momenta[1][10], ctes_.momenta[2][10], ctes_.momenta[3][10]);
00728 numstable++;
00729 }
00730
00731
00732 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
00733 p4[numstable].set(ctes_.momenta[0][2],ctes_.momenta[1][2], ctes_.momenta[2][2], ctes_.momenta[3][2]);
00734 numstable++;
00735 if( ctes_.momenta[0][3] != 0 )
00736 {
00737 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
00738 p4[numstable].set(ctes_.momenta[0][3],ctes_.momenta[1][3], ctes_.momenta[2][3], ctes_.momenta[3][3]);
00739 numstable++;
00740 }
00741
00742 int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
00743 more=(channel!=-1);
00744 if(more) {std::cout<<"Existence of mode "<<channel<<" in exclusive decay list has the same final state as this one"<<std::endl;abort(); }
00745
00746 p->makeDaughters(numstable,evtnumstable);
00747
00748 int ndaugFound=0;
00749 EvtVector4R SUMP4(0,0,0,0);
00750 for(int i=0;i<numstable;i++){
00751 p->getDaug(i)->init(evtnumstable[i],p4[i]);
00752 ndaugFound++;
00753 }
00754 if ( ndaugFound == 0 ) {
00755 report(ERROR,"EvtGen") << "Phokhara has failed to do a decay ";
00756 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass()<<endl;
00757 assert(0);
00758 }
00759
00760 nevtgen++;
00761 return ;
00762
00763 }
00764
00765
00766
00767 void EvtPhokhara_pi0pi0pipi::store(EvtDecayBase* jsdecay){
00768
00769 if (nphokharadecays==ntable){
00770
00771 EvtDecayBasePtr* newphokharadecays=new EvtDecayBasePtr[2*ntable+10];
00772 int i;
00773 for(i=0;i<ntable;i++){
00774 newphokharadecays[i]=phokharadecays[i];
00775 }
00776 ntable=2*ntable+10;
00777 delete [] phokharadecays;
00778 phokharadecays=newphokharadecays;
00779 }
00780
00781 phokharadecays[nphokharadecays++]=jsdecay;
00782
00783
00784
00785 }
00786
00787
00788 void EvtPhokhara_pi0pi0pipi::PhokharaInit(int dummy){
00789 static int first=1;
00790 if (first){
00791
00792 first=0;
00793
00794
00795 }
00796
00797 }
00798
00799
00800
00801 void EvtPhokhara_pi0pi0pipi::init_evt(EvtParticle* p){
00802 m_pion=2;
00803
00804
00805
00806
00807 if(m_pion<0 || m_pion>9){std::cout<<"mode index for phokhar 0~9, but you give "<<m_pion<<std::endl;abort();}
00808 EvtId myvpho=EvtPDL::getId("vpho");
00809 m_E = p->mass();
00811 m_tagged = 0;
00812 m_nm = 50000 ;
00813 m_nlo = 1;
00814 m_w = 0.0001;
00815 m_fsr = 1;
00816 m_fsrnlo = 1 ;
00817 m_NarrowRes = 0 ;
00818 m_FF_Kaon = 1 ;
00819 m_ivac = 0;
00820 m_FF_Pion = 0 ;
00821 m_f0_model = 0 ;
00822 m_q2min = 0.0;
00823 m_q2_min_c = 0.0447 ;
00824 m_q2_max_c = m_E*m_E;
00825 m_gmin = 0.001;
00826 m_phot1cut = 0.0;
00827 m_phot2cut = 180.0;
00828 m_pi1cut = 0.0 ;
00829 m_pi2cut = 180.0;
00830
00831 if(!(m_pion==0 || m_pion==1 || m_pion==6)){m_fsr = 0; m_fsrnlo = 0 ;}
00832 if( m_pion==9 ){m_nlo = 0 ;}
00833
00834 maxima_.iprint = 0;
00835 flags_.nlo = m_nlo;
00836 flags_.pion = m_pion;
00837 flags_.fsr = m_fsr;
00838 flags_.fsrnlo = m_fsrnlo;
00839 flags_.ivac = m_ivac;
00840 flags_.FF_pion = m_FF_Pion;
00841 flags_.f0_model = m_f0_model;
00842 flags_.FF_kaon = m_FF_Kaon;
00843 flags_.narr_res = m_NarrowRes;
00844
00845 ctes_.Sp = m_E*m_E; ;
00846
00847 cuts_.w = m_w;
00848 cuts_.q2min = m_q2min;
00849 cuts_.q2_min_c = m_q2_min_c;
00850 cuts_.q2_max_c = m_q2_max_c;
00851 cuts_.gmin = m_gmin;
00852 cuts_.phot1cut = m_phot1cut;
00853 cuts_.phot2cut = m_phot2cut;
00854 cuts_.pi1cut = m_pi1cut;
00855 cuts_.pi2cut = m_pi2cut;
00856
00857 INPUT();
00858
00859 cos1min = cos(cuts_.phot2cut*ctes_.pi/180.0);
00860 cos1max = cos(cuts_.phot1cut*ctes_.pi/180.0);
00861 cos2min = -1.0;
00862 cos2max = 1.0;
00863 cos3min = -1.0;
00864 cos3max = 1.0;
00865 if (flags_.pion == 0)
00866 qqmin = 4.0*ctes_.mmu*ctes_.mmu;
00867 else if (flags_.pion == 1)
00868 qqmin = 4.0*ctes_.mpi*ctes_.mpi;
00869 else if (flags_.pion == 2)
00870 qqmin = 4.0*(ctes_.mpi+ctes_.mpi0)*(ctes_.mpi+ctes_.mpi0);
00871 else if (flags_.pion == 3)
00872 qqmin = 16.0*ctes_.mpi*ctes_.mpi;
00873 else if (flags_.pion == 4)
00874 qqmin = 4.0*ctes_.mp*ctes_.mp;
00875 else if (flags_.pion == 5)
00876 qqmin = 4.0*ctes_.mnt*ctes_.mnt;
00877 else if (flags_.pion == 6)
00878 qqmin = 4.0*ctes_.mKp*ctes_.mKp;
00879 else if (flags_.pion == 7)
00880 qqmin = 4.0*ctes_.mKn*ctes_.mKn;
00881 else if (flags_.pion == 8)
00882 qqmin = (2.0*ctes_.mpi+ctes_.mpi0)*(2.0*ctes_.mpi+ctes_.mpi0);
00883 else if (flags_.pion == 9)
00884 qqmin = 4.0*ctes_.mlamb*ctes_.mlamb;
00885 qqmax = ctes_.Sp-2.0*sqrt(ctes_.Sp)*cuts_.gmin;
00886 if (cuts_.q2_max_c < qqmax)
00887 qqmax=cuts_.q2_max_c;
00888
00889
00890 if ( (cuts_.q2_min_c > qqmin) && (cuts_.q2_min_c < (ctes_.Sp*(1.0-2.0*(cuts_.gmin/sqrt(ctes_.Sp)+cuts_.w))) ) )
00891 qqmin = cuts_.q2_min_c;
00892 else {
00893 }
00894
00895
00896
00897
00898
00899 for(int i = 0; i<2; i++)
00900 {
00901 maxima_.Mmax[i] = 1.0;
00902 maxima_.gross[i] = 0.0;
00903 maxima_.klein[i] = 0.0;
00904 }
00905
00906 if (flags_.nlo == 0)
00907 maxima_.Mmax[1]=0.0;
00908
00909 maxima_.tr[0] = 0.0;
00910 maxima_.tr[1] = 0.0;
00911 maxima_.count[0] = 0.0;
00912 maxima_.count[1] = 0.0;
00913
00914
00915
00916 maxima_.Mmax[0] = theMmax0;
00917 maxima_.Mmax[1] = theMmax1;
00918 if((flags_.pion == 1) && (flags_.fsrnlo == 1))
00919 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
00920 if((flags_.pion == 0) && (flags_.fsrnlo == 1))
00921 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
00922
00923 if((flags_.pion == 0) && (flags_.fsr == 1) && (flags_.fsrnlo == 0))
00924 maxima_.Mmax[1]=maxima_.Mmax[1]*1.2;
00925
00926 if((flags_.pion == 2) || (flags_.pion == 3)){
00927 maxima_.Mmax[0]=maxima_.Mmax[0]*1.1;
00928 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
00929 }
00930
00931 if(flags_.pion == 8){
00932 maxima_.Mmax[0]=maxima_.Mmax[0]*1.08;
00933 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
00934 }
00935
00936
00937 maxima_.tr[0] = 0.0;
00938 maxima_.tr[1] = 0.0;
00939 maxima_.count[0] = 0.0;
00940 maxima_.count[1] = 0.0;
00941 }
00942