BOSS 7.0.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtPhokhara_4pi Class Reference

#include <EvtPhokhara_4pi.hh>

+ Inheritance diagram for EvtPhokhara_4pi:

Public Member Functions

 EvtPhokhara_4pi ()
 
virtual ~EvtPhokhara_4pi ()
 
void getName (std::string &name)
 
EvtDecayBaseclone ()
 
void decay (EvtParticle *p)
 
std::string commandName ()
 
void command (std::string cmd)
 
void init ()
 
void init_mode (EvtParticle *p)
 
void init_evt (EvtParticle *p)
 
void initProbMax ()
 
int getTotalEvt ()
 
void PhokharaInit (int dummy)
 
void ExclusiveDecay (EvtParticle *p)
 
void init_pars ()
 
- Public Member Functions inherited from EvtDecayIncoherent
void makeDecay (EvtParticle *p)
 
virtual ~EvtDecayIncoherent ()
 
void setDaughterSpinDensity (int daughter)
 
int isDaughterSpinDensitySet (int daughter)
 
- Public Member Functions inherited from EvtDecayBase
virtual void getName (std::string &name)=0
 
virtual void decay (EvtParticle *p)=0
 
virtual void makeDecay (EvtParticle *p)=0
 
virtual EvtDecayBaseclone ()=0
 
virtual void init ()
 
virtual void initProbMax ()
 
virtual std::string commandName ()
 
virtual void command (std::string cmd)
 
double getProbMax (double prob)
 
double resetProbMax (double prob)
 
 EvtDecayBase ()
 
virtual ~EvtDecayBase ()
 
virtual bool matchingDecay (const EvtDecayBase &other) const
 
EvtId getParentId ()
 
double getBranchingFraction ()
 
void disableCheckQ ()
 
void checkQ ()
 
int getNDaug ()
 
EvtIdgetDaugs ()
 
EvtId getDaug (int i)
 
int getNArg ()
 
int getPHOTOS ()
 
void setPHOTOS ()
 
void setVerbose ()
 
void setSummary ()
 
double * getArgs ()
 
std::string * getArgsStr ()
 
double getArg (int j)
 
std::string getArgStr (int j)
 
std::string getModelName ()
 
int getDSum ()
 
int summary ()
 
int verbose ()
 
void saveDecayInfo (EvtId ipar, int ndaug, EvtId *daug, int narg, std::vector< std::string > &args, std::string name, double brfr)
 
void printSummary ()
 
void setProbMax (double prbmx)
 
void noProbMax ()
 
void checkNArg (int a1, int a2=-1, int a3=-1, int a4=-1)
 
void checkNDaug (int d1, int d2=-1)
 
void checkSpinParent (EvtSpinType::spintype sp)
 
void checkSpinDaughter (int d1, EvtSpinType::spintype sp)
 
virtual int nRealDaughters ()
 

Additional Inherited Members

- Static Public Member Functions inherited from EvtDecayBase
static void findMasses (EvtParticle *p, int ndaugs, EvtId daugs[10], double masses[10])
 
static void findMass (EvtParticle *p)
 
static double findMaxMass (EvtParticle *p)
 
- Protected Member Functions inherited from EvtDecayBase
bool daugsDecayedByParentModel ()
 
- Protected Attributes inherited from EvtDecayBase
bool _daugsDecayedByParentModel
 

Detailed Description

Definition at line 33 of file EvtPhokhara_4pi.hh.

Constructor & Destructor Documentation

◆ EvtPhokhara_4pi()

EvtPhokhara_4pi::EvtPhokhara_4pi ( )

Definition at line 67 of file EvtPhokhara_4pi.cc.

67{}

Referenced by clone().

◆ ~EvtPhokhara_4pi()

EvtPhokhara_4pi::~EvtPhokhara_4pi ( )
virtual

Definition at line 68 of file EvtPhokhara_4pi.cc.

68 {
69 int i;
70 //the deletion of commands is really uggly!
71
72 if (nphokharadecays==0) {
73 delete [] commands;
74 commands=0;
75 return;
76 }
77
78 for(i=0;i<nphokharadecays;i++){
79 if (phokharadecays[i]==this){
80 phokharadecays[i]=phokharadecays[nphokharadecays-1];
81 nphokharadecays--;
82 if (nphokharadecays==0) {
83 delete [] commands;
84 commands=0;
85 }
86 return;
87 }
88 }
89
90 report(ERROR,"EvtGen") << "Error in destroying Phokhara model!"<<endl;
91
92}
ostream & report(Severity severity, const char *facility)
Definition: EvtReport.cc:36
@ ERROR
Definition: EvtReport.hh:49

Member Function Documentation

◆ clone()

EvtDecayBase * EvtPhokhara_4pi::clone ( )
virtual

Implements EvtDecayBase.

Definition at line 101 of file EvtPhokhara_4pi.cc.

101 {
102
103 return new EvtPhokhara_4pi;
104
105}

◆ command()

void EvtPhokhara_4pi::command ( std::string  cmd)
virtual

Reimplemented from EvtDecayBase.

Definition at line 508 of file EvtPhokhara_4pi.cc.

508 {
509
510 if (ncommand==lcommand){
511
512 lcommand=10+2*lcommand;
513
514 std::string* newcommands=new std::string[lcommand];
515
516 int i;
517
518 for(i=0;i<ncommand;i++){
519 newcommands[i]=commands[i];
520 }
521
522 delete [] commands;
523
524 commands=newcommands;
525
526 }
527
528 commands[ncommand]=cmd;
529
530 ncommand++;
531
532}

◆ commandName()

std::string EvtPhokhara_4pi::commandName ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 502 of file EvtPhokhara_4pi.cc.

502 {
503
504 return std::string("PhokharaPar");
505
506}

◆ decay()

void EvtPhokhara_4pi::decay ( EvtParticle p)
virtual

Implements EvtDecayBase.

Definition at line 536 of file EvtPhokhara_4pi.cc.

536 {
537 EvtId myvpho=EvtPDL::getId("vpho");
538 if(p->getId()!=myvpho) {std::cout<<"Parent particle is required to be vpho for Phokhara model"<<std::endl;abort();}
539 if(nevtgen==0) {init_mode(p);}
540 else{init_evt(p);}
541 std::cout<<"PHOKHARA : 2(pi+pi-) mode "<<std::endl;
542
543
544 int istdheppar=EvtPDL::getStdHep(p->getId());
545 int ntrials = 0;
546 int tr_old[2];
547 tr_old[0] = (int)maxima_.tr[0];
548 tr_old[1] = (int)maxima_.tr[1];
549
550 while( ntrials < 10000)
551 {
552 ievent++;
553 RANLXDF(Ar_r,1);
554 Ar[1] = Ar_r[0];
555
556 if (Ar[1] <= (maxima_.Mmax[0]/(maxima_.Mmax[0]+maxima_.Mmax[1]))) {
557 maxima_.count[0] = maxima_.count[0]+1.0;
558 GEN_1PH(2,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
559 }
560 else {
561 maxima_.count[1] = maxima_.count[1]+1.0;
562 GEN_2PH(2,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
563 }
564
565 if( ((int)maxima_.tr[0]+(int)maxima_.tr[1]) > (tr_old[0]+tr_old[1]) ) // event accepted after cuts
566 {
567 goto storedEvents;
568 }
569 ntrials ++;
570 }
571
572 std::cout <<"FATAL: Could not satisfy cuts after " << ntrials << "trials. Terminate." <<std::endl;
573 //----
574 storedEvents:
575 int more=0;
576 int numstable=0;
577 int numparton=0;
578 EvtId evtnumstable[100];//
579 EvtVector4R p4[20];
580
581 // except ISR photos, products depending on channel
582 if (flags_.pion == 0) { // mu+ mu-
583 // mu+
584 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-13);
585 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
586 numstable++;
587 // mu -
588 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(13);
589 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
590 numstable++;
591 }
592 if (flags_.pion == 1) { // pi+ pi-
593 // pi+
594 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
595 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
596 numstable++;
597 // pi -
598 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
599 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
600 numstable++;
601 }
602 if (flags_.pion == 2) { // pi+ pi-2pi0
603 // pi0
604 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
605 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
606 numstable++;
607 // pi0
608 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
609 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
610 numstable++;
611 // pi-
612 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
613 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
614 numstable++;
615 // pi +
616 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
617 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
618 numstable++;
619 }
620 if (flags_.pion == 3) { // 2(pi+ pi-)
621 // pi+
622 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
623 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
624 numstable++;
625 // pi-
626 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
627 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
628 numstable++;
629 // pi+
630 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
631 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
632 numstable++;
633 // pi -
634 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
635 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
636 numstable++;
637 }
638 if (flags_.pion == 4) { // ppbar
639 // pbar
640 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2212);
641 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
642 numstable++;
643 // p
644 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2212);
645 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
646 numstable++;
647 }
648 if (flags_.pion == 5) { // nnbar
649 // pbar
650 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2112);
651 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
652 numstable++;
653 // p
654 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2112);
655 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
656 numstable++;
657 }
658 if (flags_.pion == 6) { // K+ K-
659 // K+
660 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(321);
661 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
662 numstable++;
663 // K -
664 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-321);
665 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
666 numstable++;
667 }
668 if (flags_.pion == 7) { // K0K0bar
669 // Kbar
670 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(311);
671 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
672 numstable++;
673 // K0
674 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-311);
675 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
676 numstable++;
677 }
678 if (flags_.pion == 8) { // pi+ pi-pi0
679 // pi+
680 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
681 p4[numstable].set(ctes_.momenta[0][5],ctes_.momenta[1][5], ctes_.momenta[2][5], ctes_.momenta[3][5]);
682 numstable++;
683 // pi-
684 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
685 p4[numstable].set(ctes_.momenta[0][6],ctes_.momenta[1][6], ctes_.momenta[2][6], ctes_.momenta[3][6]);
686 numstable++;
687 // pi0
688 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(111);
689 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
690 numstable++;
691 }
692 if (flags_.pion == 9) { //Lambda Lambdabar-> pi+ pi- ppbar
693 // pi+
694 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(211);
695 p4[numstable].set(ctes_.momenta[0][7],ctes_.momenta[1][7], ctes_.momenta[2][7], ctes_.momenta[3][7]);
696 numstable++;
697 // pbar
698 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-2212);
699 p4[numstable].set(ctes_.momenta[0][8],ctes_.momenta[1][8], ctes_.momenta[2][8], ctes_.momenta[3][8]);
700 numstable++;
701 // pi-
702 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(-211);
703 p4[numstable].set(ctes_.momenta[0][9],ctes_.momenta[1][9], ctes_.momenta[2][9], ctes_.momenta[3][9]);
704 numstable++;
705 // p
706 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(2212);
707 p4[numstable].set(ctes_.momenta[0][10],ctes_.momenta[1][10], ctes_.momenta[2][10], ctes_.momenta[3][10]);
708 numstable++;
709 }
710
711 // ISR gamma
712 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
713 p4[numstable].set(ctes_.momenta[0][2],ctes_.momenta[1][2], ctes_.momenta[2][2], ctes_.momenta[3][2]);
714 numstable++;
715 if( ctes_.momenta[0][3] != 0 ) // second photon exists
716 {
717 evtnumstable[numstable]=EvtPDL::evtIdFromStdHep(22);
718 p4[numstable].set(ctes_.momenta[0][3],ctes_.momenta[1][3], ctes_.momenta[2][3], ctes_.momenta[3][3]);
719 numstable++;
720 }
721
722 int channel=EvtDecayTable::inChannelList(p->getId(),numstable,evtnumstable);
723 more=(channel!=-1);
724 if(more) {std::cout<<"Existence of mode "<<channel<<" in exclusive decay list has the same final state as this one"<<std::endl;abort(); }
725
726 p->makeDaughters(numstable,evtnumstable);
727 //double weight = p->initializePhaseSpace(getNDaug(),getDaugs());
728
729 int ndaugFound=0;
730 EvtVector4R SUMP4(0,0,0,0);
731 for(int i=0;i<numstable;i++){
732 p->getDaug(i)->init(evtnumstable[i],p4[i]);
733 ndaugFound++;
734 }
735 if ( ndaugFound == 0 ) {
736 report(ERROR,"EvtGen") << "Phokhara has failed to do a decay ";
737 report(ERROR,"EvtGen") << EvtPDL::name(p->getId()).c_str() << " " << p->mass()<<endl;
738 assert(0);
739 }
740
741 nevtgen++;
742 return ;
743
744}
struct @17 flags_
#define GEN_1PH(I, QQMIN, QQMAX, COS1MIN, COS1MAX, COS3MIN, COS3MAX)
#define RANLXDF(AR, VAL)
struct @16 maxima_
struct @10 ctes_
#define GEN_2PH(I, QQMIN, COS1MIN, COS1MAX, COS2MIN, COS2MAX, COS3MIN, COS3MAX)
static int inChannelList(EvtId parent, int ndaug, EvtId *daugs)
Definition: EvtId.hh:27
static int getStdHep(EvtId id)
Definition: EvtPDL.hh:56
static EvtId evtIdFromStdHep(int stdhep)
Definition: EvtPDL.cc:244
static std::string name(EvtId i)
Definition: EvtPDL.hh:64
static EvtId getId(const std::string &name)
Definition: EvtPDL.cc:287
void makeDaughters(int ndaug, EvtId *id)
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtId getId() const
Definition: EvtParticle.cc:113
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:85
double mass() const
Definition: EvtParticle.cc:127
void init_evt(EvtParticle *p)
void init_mode(EvtParticle *p)
void set(int i, double d)
Definition: EvtVector4R.hh:183

◆ ExclusiveDecay()

void EvtPhokhara_4pi::ExclusiveDecay ( EvtParticle p)

◆ getName()

void EvtPhokhara_4pi::getName ( std::string &  name)
virtual

Implements EvtDecayBase.

Definition at line 95 of file EvtPhokhara_4pi.cc.

95 {
96
97 model_name="PHOKHARA_4pi";
98
99}

◆ getTotalEvt()

int EvtPhokhara_4pi::getTotalEvt ( )
inline

Definition at line 52 of file EvtPhokhara_4pi.hh.

52{return nevt;}

◆ init()

void EvtPhokhara_4pi::init ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 477 of file EvtPhokhara_4pi.cc.

477 {
478 checkNArg(0);
479
480 std::string locvp=getenv("BESEVTGENROOT");
481 system("cat $BESEVTGENROOT/share/phokhara.param>phokhara.param");
482
483
484 if (getParentId().isAlias()){
485
486 report(ERROR,"EvtGen") << "EvtPhokhara finds that you are decaying the"<<endl
487 << " aliased particle "
488 << EvtPDL::name(getParentId()).c_str()
489 << " with the Phokhara model"<<endl
490 << " this does not work, please modify decay table."
491 << endl;
492 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
493 ::abort();
494
495 }
496
497 store(this);
498
499}
EvtId getParentId()
Definition: EvtDecayBase.hh:60
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)

◆ init_evt()

void EvtPhokhara_4pi::init_evt ( EvtParticle p)

======== list parameters to be initialized

Definition at line 782 of file EvtPhokhara_4pi.cc.

782 {
783 m_pion=3;
784 // mu+mu-(0),pi+pi-(1),2pi0pi+pi-(2),
785 // 2pi+2pi-(3),ppbar(4),nnbar(5),
786 // K+K-(6),K0K0bar(7),pi+pi-pi0(8),
787 // Lamb Lambbar->pi-pi+ppbar(9)
788 if(m_pion<0 || m_pion>9){std::cout<<"mode index for phokhar 0~9, but you give "<<m_pion<<std::endl;abort();}
789 EvtId myvpho=EvtPDL::getId("vpho");
790 m_E = p->mass();//EvtPDL::getMeanMass(myvpho);
791///======== list parameters to be initialized
792 init_pars();
793// --- input parameter initialization -----------
794 maxima_.iprint = 0;
795 flags_.nlo = m_nlo;
796 flags_.pion = m_pion;
797 flags_.fsr = m_fsr;
798 flags_.fsrnlo = m_fsrnlo;
799 flags_.ivac = m_ivac;
800 flags_.FF_pion = m_FF_Pion;
801 flags_.f0_model = m_f0_model;
802 flags_.FF_kaon = m_FF_Kaon;
803 flags_.narr_res = m_NarrowRes;
804
805 ctes_.Sp = m_E*m_E; ;
806
807 cuts_.w = m_w;
808 cuts_.q2min = m_q2min;
809 cuts_.q2_min_c = m_q2_min_c;
810 cuts_.q2_max_c = m_q2_max_c;
811 cuts_.gmin = m_gmin;
812 cuts_.phot1cut = m_phot1cut;
813 cuts_.phot2cut = m_phot2cut;
814 cuts_.pi1cut = m_pi1cut;
815 cuts_.pi2cut = m_pi2cut;
816
817 INPUT();
818
819 cos1min = cos(cuts_.phot2cut*ctes_.pi/180.0); // photon1 angle cuts in the
820 cos1max = cos(cuts_.phot1cut*ctes_.pi/180.0); // LAB rest frame
821 cos2min = -1.0; // photon2 angle limits
822 cos2max = 1.0; //
823 cos3min = -1.0; // hadrons/muons angle limits
824 cos3max = 1.0; // in their rest frame
825 if (flags_.pion == 0) // virtual photon energy cut
826 qqmin = 4.0*ctes_.mmu*ctes_.mmu;
827 else if (flags_.pion == 1)
828 qqmin = 4.0*ctes_.mpi*ctes_.mpi;
829 else if (flags_.pion == 2)
830 qqmin = 4.0*(ctes_.mpi+ctes_.mpi0)*(ctes_.mpi+ctes_.mpi0);
831 else if (flags_.pion == 3)
832 qqmin = 16.0*ctes_.mpi*ctes_.mpi;
833 else if (flags_.pion == 4)
834 qqmin = 4.0*ctes_.mp*ctes_.mp;
835 else if (flags_.pion == 5)
836 qqmin = 4.0*ctes_.mnt*ctes_.mnt;
837 else if (flags_.pion == 6)
838 qqmin = 4.0*ctes_.mKp*ctes_.mKp;
839 else if (flags_.pion == 7)
840 qqmin = 4.0*ctes_.mKn*ctes_.mKn;
841 else if (flags_.pion == 8)
842 qqmin = (2.0*ctes_.mpi+ctes_.mpi0)*(2.0*ctes_.mpi+ctes_.mpi0);
843 else if (flags_.pion == 9)
844 qqmin = 4.0*ctes_.mlamb*ctes_.mlamb;
845 qqmax = ctes_.Sp-2.0*sqrt(ctes_.Sp)*cuts_.gmin; // if only one photon
846 if (cuts_.q2_max_c < qqmax)
847 qqmax=cuts_.q2_max_c; // external cuts
848
849 // -------------------
850 if ( (cuts_.q2_min_c > qqmin) && (cuts_.q2_min_c < (ctes_.Sp*(1.0-2.0*(cuts_.gmin/sqrt(ctes_.Sp)+cuts_.w))) ) )
851 qqmin = cuts_.q2_min_c;
852 else {
853 }
854
855
856 // =================================================
857// --- finding the maximum -------------------------
858 for(int i = 0; i<2; i++)
859 {
860 maxima_.Mmax[i] = 1.0;
861 maxima_.gross[i] = 0.0;
862 maxima_.klein[i] = 0.0;
863 }
864 if (flags_.nlo == 0)
865 maxima_.Mmax[1]=0.0; // only 1 photon events generated
866
867 maxima_.tr[0] = 0.0;
868 maxima_.tr[1] = 0.0;
869 maxima_.count[0] = 0.0;
870 maxima_.count[1] = 0.0;
871
872 // =================================================
873 // --- for the second run ---
874 maxima_.Mmax[0] = theMmax0;
875 maxima_.Mmax[1] = theMmax1;
876 if((flags_.pion == 1) && (flags_.fsrnlo == 1))
877 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
878 if((flags_.pion == 0) && (flags_.fsrnlo == 1))
879 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
880
881 if((flags_.pion == 0) && (flags_.fsr == 1) && (flags_.fsrnlo == 0))
882 maxima_.Mmax[1]=maxima_.Mmax[1]*1.2;
883
884 if((flags_.pion == 2) || (flags_.pion == 3)){
885 maxima_.Mmax[0]=maxima_.Mmax[0]*1.1;
886 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
887 }
888
889 if(flags_.pion == 8){
890 maxima_.Mmax[0]=maxima_.Mmax[0]*1.08;
891 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
892 }
893// --- end of the second run -----------------------
894
895 maxima_.tr[0] = 0.0;
896 maxima_.tr[1] = 0.0;
897 maxima_.count[0] = 0.0;
898 maxima_.count[1] = 0.0;
899}
#define INPUT
Definition: BesBdkRc.cxx:35
struct @11 cuts_
double cos(const BesAngle a)

Referenced by decay().

◆ init_mode()

void EvtPhokhara_4pi::init_mode ( EvtParticle p)

======== list parameters to be initialized

Definition at line 115 of file EvtPhokhara_4pi.cc.

115 {
116 m_pion=3;
117 // mu+mu-(0),pi+pi-(1),2pi0pi+pi-(2),
118 // 2pi+2pi-(3),ppbar(4),nnbar(5),
119 // K+K-(6),K0K0bar(7),pi+pi-pi0(8),
120 // Lamb Lambbar->pi-pi+ppbar(9)
121 if(m_pion<0 || m_pion>9){std::cout<<"mode index for phokhar 0~9, but you give "<<m_pion<<std::endl;abort();}
122 EvtId myvpho=EvtPDL::getId("vpho");
123 m_E = p->mass();//EvtPDL::getMeanMass(myvpho);
124///======== list parameters to be initialized
125 init_pars();
126 // --- input parameter initialization -----------
127 m_initSeed = 123456789;
128 RLXDINIT(1, m_initSeed);
129 maxima_.iprint = 0;
130 flags_.nlo = m_nlo;
131 flags_.pion = m_pion;
132 flags_.fsr = m_fsr;
133 flags_.fsrnlo = m_fsrnlo;
134 flags_.ivac = m_ivac;
135 flags_.FF_pion = m_FF_Pion;
136 flags_.f0_model = m_f0_model;
137 flags_.FF_kaon = m_FF_Kaon;
138 flags_.narr_res = m_NarrowRes;
139
140 ctes_.Sp = m_E*m_E; ;
141
142 cuts_.w = m_w;
143 cuts_.q2min = m_q2min;
144 cuts_.q2_min_c = m_q2_min_c;
145 cuts_.q2_max_c = m_q2_max_c;
146 cuts_.gmin = m_gmin;
147 cuts_.phot1cut = m_phot1cut;
148 cuts_.phot2cut = m_phot2cut;
149 cuts_.pi1cut = m_pi1cut;
150 cuts_.pi2cut = m_pi2cut;
151
152 INPUT();
153
154 // --- print run data ------------------------------
155 cout << "-------------------------------------------------------------" << endl;
156 if (flags_.pion == 0)
157 cout << " PHOKHARA 7.0 : e^+ e^- -> mu^+ mu^- gamma" << endl;
158 else if (flags_.pion == 1)
159 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- gamma" << endl;
160 else if (flags_.pion == 2)
161 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- 2pi^0 gamma" << endl;
162 else if (flags_.pion == 3)
163 cout << " PHOKHARA 7.0: e^+ e^- -> 2pi^+ 2pi^- gamma" << endl;
164 else if (flags_.pion == 4)
165 cout << " PHOKHARA 7.0: e^+ e^- -> p pbar gamma" << endl;
166 else if (flags_.pion == 5)
167 cout << " PHOKHARA 7.0: e^+ e^- -> n nbar gamma" << endl;
168 else if (flags_.pion == 6)
169 cout << " PHOKHARA 7.0: e^+ e^- -> K^+ K^- gamma" << endl;
170 else if (flags_.pion == 7)
171 cout << " PHOKHARA 7.0: e^+ e^- -> K_0 K_0bar gamma" << endl;
172 else if (flags_.pion == 8)
173 cout << " PHOKHARA 7.0: e^+ e^- -> pi^+ pi^- pi^0 gamma" << endl;
174 else if (flags_.pion == 9) {
175 cout << "PHOKHARA 7.0 : e^+ e^- ->" << endl;
176 cout << " Lambda (-> pi^- p) Lambda bar (-> pi^+ pbar) gamma" << endl;
177 } else
178 cout << " PHOKHARA 7.0: not yet implemented" << endl;
179
180 // --------------------------------
181 cout << "--------------------------------------------------------------" << endl;
182 printf(" %s %f %s\n","cms total energy = ",sqrt(ctes_.Sp)," GeV ");
183
184 //if((cuts_.gmin/2.0/ctes_.ebeam) < 0.000098){
185 if(cuts_.gmin<0.001){
186 cerr << " minimal missing energy set too small, abort" << endl;
187 abort();
188 }
189 printf(" %s %f %s\n","minimal tagged photon energy = ",cuts_.gmin," GeV ");
190 printf(" %s %f,%f\n","angular cuts on tagged photon = ",cuts_.phot1cut,cuts_.phot2cut);
191
192 // --------------------------------
193 if (flags_.pion == 0)
194 printf(" %s %f,%f\n","angular cuts on muons = ",cuts_.pi1cut,cuts_.pi2cut);
195 else if (flags_.pion == 4)
196 printf(" %s %f,%f\n","angular cuts on protons = ",cuts_.pi1cut,cuts_.pi2cut);
197 else if (flags_.pion == 5)
198 printf(" %s %f,%f\n","angular cuts on neutrons = ", cuts_.pi1cut,cuts_.pi2cut);
199 else if ((flags_.pion == 6)||(flags_.pion == 7))
200 printf(" %s %f,%f\n","angular cuts on kaons = ", cuts_.pi1cut,cuts_.pi2cut);
201 else if (flags_.pion == 9)
202 printf(" %s %f,%f\n","angular cuts on pions and protons = ", cuts_.pi1cut,cuts_.pi2cut);
203 else
204 printf(" %s %f,%f\n","angular cuts on pions = ", cuts_.pi1cut,cuts_.pi2cut);
205 if (flags_.pion == 0)
206 printf(" %s %f %s\n","min. muons-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2");
207 else if (flags_.pion == 4)
208 printf(" %s %f %s\n","min. protons-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2" );
209 else if (flags_.pion == 5)
210 printf(" %s %f %s\n","min. neutrons-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2" );
211 else if ((flags_.pion == 6)||(flags_.pion == 7))
212 printf(" %s %f %s\n","min. kaons-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2" );
213 else if (flags_.pion == 9)
214 printf(" %s %f %s\n"," min. lambdas-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2" );
215 else
216 printf(" %s %f %s\n","min. pions-tagged photon inv.mass^2 = ", cuts_.q2min," GeV^2" );
217 cos1min = cos(cuts_.phot2cut*ctes_.pi/180.0); // photon1 angle cuts in the
218 cos1max = cos(cuts_.phot1cut*ctes_.pi/180.0); // LAB rest frame
219 cos2min = -1.0; // photon2 angle limits
220 cos2max = 1.0; //
221 cos3min = -1.0; // hadrons/muons angle limits
222 cos3max = 1.0; // in their rest frame
223 if (flags_.pion == 0) // virtual photon energy cut
224 qqmin = 4.0*ctes_.mmu*ctes_.mmu;
225 else if (flags_.pion == 1)
226 qqmin = 4.0*ctes_.mpi*ctes_.mpi;
227 else if (flags_.pion == 2)
228 qqmin = 4.0*(ctes_.mpi+ctes_.mpi0)*(ctes_.mpi+ctes_.mpi0);
229 else if (flags_.pion == 3)
230 qqmin = 16.0*ctes_.mpi*ctes_.mpi;
231 else if (flags_.pion == 4)
232 qqmin = 4.0*ctes_.mp*ctes_.mp;
233 else if (flags_.pion == 5)
234 qqmin = 4.0*ctes_.mnt*ctes_.mnt;
235 else if (flags_.pion == 6)
236 qqmin = 4.0*ctes_.mKp*ctes_.mKp;
237 else if (flags_.pion == 7)
238 qqmin = 4.0*ctes_.mKn*ctes_.mKn;
239 else if (flags_.pion == 8)
240 qqmin = (2.0*ctes_.mpi+ctes_.mpi0)*(2.0*ctes_.mpi+ctes_.mpi0);
241 else if (flags_.pion == 9)
242 qqmin = 4.0*ctes_.mlamb*ctes_.mlamb;
243 qqmax = ctes_.Sp-2.0*sqrt(ctes_.Sp)*cuts_.gmin; // if only one photon
244 if (cuts_.q2_max_c < qqmax)
245 qqmax=cuts_.q2_max_c; // external cuts
246
247 // -------------------
248 if ( (cuts_.q2_min_c > qqmin) && (cuts_.q2_min_c < (ctes_.Sp*(1.0-2.0*(cuts_.gmin/sqrt(ctes_.Sp)+cuts_.w))) ) )
249 qqmin = cuts_.q2_min_c;
250 else {
251 cerr << "------------------------------" << endl;
252 cerr << " Q^2_min TOO SMALL" << endl;
253 cerr << " Q^2_min CHANGED BY PHOKHARA = " << qqmin << " GeV^2" << endl;
254 cerr << "------------------------------" << endl;
255 }
256 // -------------------
257 if(qqmax <= qqmin){
258 cerr << " Q^2_max to small " << endl;
259 cerr << " Q^2_max = " << qqmax << endl;
260 cerr << " Q^2_min = " << qqmin << endl;
261 abort();
262 }
263
264 // -------------------
265 if (flags_.pion == 0) {
266 printf(" %s %f %s\n", "minimal muon-pair invariant mass^2 = ", qqmin," GeV^2");
267 printf(" %s %f %s\n", "maximal muon-pair invariant mass^2 = ", qqmax," GeV^2");
268 } else if (flags_.pion == 1) {
269 printf(" %s %f %s\n", "minimal pion-pair invariant mass^2 = ", qqmin," GeV^2");
270 printf(" %s %f %s\n", "maximal pion-pair invariant mass^2 = ", qqmax," GeV^2");
271 } else if (flags_.pion == 4) {
272 printf(" %s %f %s\n", "minimal proton-pair invariant mass^2 = ", qqmin," GeV^2");
273 printf(" %s %f %s\n", "maximal proton-pair invariant mass^2 = ", qqmax," GeV^2");
274 } else if (flags_.pion == 5) {
275 printf(" %s %f %s\n", "minimal neutron-pair invariant mass^2 = ", qqmin," GeV^2");
276 printf(" %s %f %s\n", "maximal neutron-pair invariant mass^2 = ", qqmax," GeV^2");
277 } else if ((flags_.pion == 6) || (flags_.pion == 7)) {
278 printf(" %s %f %s\n", "minimal kaon-pair invariant mass^2 = ", qqmin," GeV^2");
279 printf(" %s %f %s\n", "maximal kaon-pair invariant mass^2 = ", qqmax," GeV^2");
280 } else if(flags_.pion == 8){
281 printf(" %s %f %s\n", "minimal three-pion invariant mass^2 = ", qqmin," GeV^2");
282 printf(" %s %f %s\n", "maximal three-pion invariant mass^2 = ", qqmax," GeV^2");
283 } else if(flags_.pion == 9){
284 printf(" %s %f %s\n", "minimal lambda-pair invariant mass^2 = ", qqmin," GeV^2");
285 printf(" %s %f %s\n", "maximal lambda-pair invariant mass^2 = ", qqmax," GeV^2");
286 } else {
287 printf(" %s %f %s\n", "minimal four-pion invariant mass^2 = ", qqmin," GeV^2" );
288 printf(" %s %f %s\n", "maximal four-pion invariant mass^2 = ", qqmax," GeV^2");
289 }
290 // -------------------
291 if (flags_.nlo == 0) {
292 cout << "Born" << endl;
293 if(flags_.fsrnlo != 0){
294 cerr << "WRONG FSRNLO flag - only fsrnlo=0 allowed for Born" << endl;
295 abort();
296 }
297 }
298 // -------------------
299 if((flags_.pion == 9) && (flags_.nlo != 0)) {
300 cerr << "WRONG NLO flag - only Born allowed for Lambdas"<< endl;
301 cerr << "If you feel that you need NLO"<< endl;
302 cerr << "please contact the authors"<< endl;
303 abort();
304 }
305 // -------------------
306 if (flags_.nlo == 1) printf(" %s %f\n", "NLO: soft photon cutoff w = ",cuts_.w);
307 if ((flags_.pion <= 1) || (flags_.pion == 6))
308 {
309
310 if( ! ((flags_.fsr == 1) || (flags_.fsr == 2) || (flags_.fsrnlo == 0)
311 || (flags_.fsr == 1) || (flags_.fsrnlo == 1)
312 || (flags_.fsr == 0) || (flags_.fsrnlo == 0))) {
313 cerr << "WRONG combination of FSR, FSRNLO flags" <<endl;
314 abort();
315 }
316 // ------------------
317 if (flags_.fsr == 0)
318 cout << "ISR only" << endl;
319 else if (flags_.fsr == 1)
320 cout << "ISR+FSR" << endl;
321 else if (flags_.fsr == 2) {
322 if (flags_.nlo == 0)
323 cout << "ISR+INT+FSR" << endl;
324 else {
325 cerr << "WRONG FSR flag: interference is included only for nlo=0" << endl;
326 abort();
327 }
328 }
329 else {
330 cerr << "WRONG FSR flag" << flags_.fsr << endl;
331 abort();
332 }
333 if(flags_.fsrnlo == 1)
334 cout << "IFSNLO included" << endl;
335 }
336 else
337 {
338 if((flags_.fsr == 0) && (flags_.fsrnlo == 0))
339 cout << "ISR only" << endl;
340 else {
341 cerr << "FSR is implemented only for pi+pi-, mu+mu- and K+K- modes" << endl;
342 abort();
343 }
344 }
345 // ------------------
346 if(flags_.ivac == 0){
347 cout << "Vacuum polarization is NOT included" << endl;}
348 else if(flags_.ivac == 1){
349 cout << "Vacuum polarization by Fred Jegerlehner (http://www-com.physik.hu-berlin.de/fjeger/alphaQEDn.uu)" << endl;}
350 else if(flags_.ivac == 2){
351 cout << "Vacuum polarization (VP_HLMNT_v1_3nonr) by Daisuke Nomura and Thomas Teubner" << endl;}
352 else {
353 cout << "WRONG vacuum polarization switch" << endl;
354 abort();
355 }
356
357// -----------------
358 if(flags_.pion == 1){
359 if(flags_.FF_pion == 0)
360 cout << "Kuhn-Santamaria PionFormFactor" << endl;
361 else if(flags_.FF_pion == 1)
362 cout << "Gounaris-Sakurai PionFormFactor old" << endl;
363 else if(flags_.FF_pion == 2)
364 cout << "Gounaris-Sakurai PionFormFactor new" << endl;
365 else {
366 cout << "WRONG PionFormFactor switch" << endl;
367 abort();
368 }
369 if(flags_.fsr != 0){
370 if(flags_.f0_model == 0)
371 cout << "f0+f0(600): K+K- model" << endl;
372 else if(flags_.f0_model == 1)
373 cout << "f0+f0(600): \"no structure\" model" << endl;
374 else if(flags_.f0_model == 2)
375 cout << "NO f0+f0(600)" << endl;
376 else if(flags_.f0_model == 3)
377 cout << "only f0, KLOE: Cesare Bini-private communication" << endl;
378 else {
379 cout << "WRONG f0+f0(600) switch" << endl;
380 abort();
381 }
382 }
383 }
384//-------
385 if((flags_.pion == 6) || (flags_.pion==7)){
386 if(flags_.FF_kaon == 0)
387 cout << "constrained KaonFormFactor" << endl;
388 else if(flags_.FF_kaon == 1) {
389 cout << "unconstrained KaonFormFactor" << endl;}
390 else if(flags_.FF_kaon == 2) {
391 cout << "Kuhn-Khodjamirian-Bruch KaonFormFactor" << endl;}
392 else{
393 cout << "WRONG KaonFormFactor switch" << endl;
394 abort();
395 }
396 }
397// --------------------
398 if((flags_.pion == 0) || (flags_.pion==1) || (flags_.pion == 6) || (flags_.pion == 7)){
399 if(flags_.narr_res == 1){
400 cout << "THE NARROW RESONANCE J/PSI INCLUDED" << endl;}
401 else if(flags_.narr_res == 2){
402 cout << "THE NARROW RESONANCE PSI(2S) INCLUDED" << endl;}
403 else if(flags_.narr_res != 0){
404 cout << "wrong flag narr_res: only 0, 1 or 2 allowed" << endl;
405 abort();
406 }
407 }
408// ------
409 cout << "--------------------------------------------------------------" << endl;
410//
411 // =================================================
412// --- finding the maximum -------------------------
413 for(int i = 0; i<2; i++)
414 {
415 maxima_.Mmax[i] = 1.0;
416 maxima_.gross[i] = 0.0;
417 maxima_.klein[i] = 0.0;
418 }
419
420 if (flags_.nlo == 0)
421 maxima_.Mmax[1]=0.0; // only 1 photon events generated
422
423 maxima_.tr[0] = 0.0;
424 maxima_.tr[1] = 0.0;
425 maxima_.count[0] = 0.0;
426 maxima_.count[1] = 0.0;
427
428 // =================================================
429 // --- beginning the MC loop event generation ------
430 for(int j = 1; j <= m_nm; j++)
431 {
432 RANLXDF(Ar_r,1);
433 Ar[1] = Ar_r[0];
434 if (Ar[1] <= (maxima_.Mmax[0]/(maxima_.Mmax[0]+maxima_.Mmax[1]))) {
435 maxima_.count[0] = maxima_.count[0]+1.0;
436 GEN_1PH(1,qqmin,qqmax,cos1min,cos1max,cos3min,cos3max);
437 }
438 else {
439 maxima_.count[1] = maxima_.count[1]+1.0;
440 GEN_2PH(1,qqmin,cos1min,cos1max,cos2min,cos2max,cos3min,cos3max);
441 }
442 }
443 // --- end of the MC loop --------------------------
444 // =================================================
445 // --- for the second run ---
446 maxima_.Mmax[0] = maxima_.gross[0]+.05*sqrt(maxima_.gross[0]*maxima_.gross[0]);
447 maxima_.Mmax[1] = maxima_.gross[1]+(.03+.02*ctes_.Sp)*sqrt(maxima_.gross[1]*maxima_.gross[1]);
448 theMmax0=maxima_.Mmax[0];
449 theMmax1=maxima_.Mmax[1];
450 if((flags_.pion == 1) && (flags_.fsrnlo == 1))
451 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
452 if((flags_.pion == 0) && (flags_.fsrnlo == 1))
453 maxima_.Mmax[1]=maxima_.Mmax[1]*1.5;
454
455 if((flags_.pion == 0) && (flags_.fsr == 1) && (flags_.fsrnlo == 0))
456 maxima_.Mmax[1]=maxima_.Mmax[1]*1.2;
457
458 if((flags_.pion == 2) || (flags_.pion == 3)){
459 maxima_.Mmax[0]=maxima_.Mmax[0]*1.1;
460 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
461 }
462
463 if(flags_.pion == 8){
464 maxima_.Mmax[0]=maxima_.Mmax[0]*1.08;
465 maxima_.Mmax[1]=maxima_.Mmax[1]*1.1;
466 }
467// --- end of the second run -----------------------
468
469 maxima_.tr[0] = 0.0;
470 maxima_.tr[1] = 0.0;
471 maxima_.count[0] = 0.0;
472 maxima_.count[1] = 0.0;
473}
#define RLXDINIT(LUXURY, SEED)

Referenced by decay().

◆ init_pars()

void EvtPhokhara_4pi::init_pars ( )

Definition at line 902 of file EvtPhokhara_4pi.cc.

902 {
903 m_tagged = 0;
904 m_nm = 50000 ; // # of events to determine the maximum
905 m_nlo = 1; // Born(0), NLO(1)
906 m_w = 0.00001; // soft photon cutoff
907 m_fsr = 0; // ISR only(0), ISR+FSR(1), ISR+INT+FSR(2)
908 m_fsrnlo = 0 ; // yes(1), no(0)
909 m_NarrowRes = 0 ;// none(0) jpsi(1) psip(2)
910 m_FF_Kaon = 1 ; // constrained (0), unconstrained (1), Kuhn-Khodjamirian-Bruch (2)
911 m_ivac = 0; // yes(1), no(0)
912 m_FF_Pion = 0 ; // KS Pionformfactor(0), GS Pionformfactor old(1) and new(2)
913 m_f0_model = 0 ; // f0+f0(600): KK model(0), no structure(1), no f0+f0(600)(2), f0 KLOE(3)
914 m_q2min = (4*0.135)*(4*0.135); // minimal hadrons(muons)-gamma-inv mass squared (GeV^2)
915 m_q2_min_c = m_q2min ; // minimal inv. mass squared of the hadrons(muons)(GeV^2)
916 m_q2_max_c = m_E*m_E; // maximal inv. mass squared of the hadrons(muons)(GeV^2)
917 m_gmin = 0.001; // minimal photon energy/missing energy (GeV) //total Ecm will reduce m_gmin
918 m_phot1cut = 0.0; // maximal photon angle/missing momentum angle (grad)
919 m_phot2cut = 180.0; // maximal photon angle/missing momentum angle (grad)
920 m_pi1cut = 0.0 ; // minimal hadrons(muons) angle (grad)
921 m_pi2cut = 180.0; // maximal hadrons(muons) angle (grad)
922 if(!(m_pion==0 || m_pion==1 || m_pion==6)){m_fsr = 0; m_fsrnlo = 0 ;}
923 if( m_pion==9 ){m_nlo = 0 ;}
924}

Referenced by init_evt(), and init_mode().

◆ initProbMax()

void EvtPhokhara_4pi::initProbMax ( )
virtual

Reimplemented from EvtDecayBase.

Definition at line 108 of file EvtPhokhara_4pi.cc.

108 {
109
110 noProbMax();
111
112}
void noProbMax()

◆ PhokharaInit()

void EvtPhokhara_4pi::PhokharaInit ( int  dummy)

Definition at line 769 of file EvtPhokhara_4pi.cc.

769 {
770 static int first=1;
771 if (first){
772
773 first=0;
774 //for(int i=0;i<ncommand;i++)
775 // lugive0_(commands[i].c_str(),strlen(commands[i].c_str()));
776 }
777
778}
Index first(Pair i)
Definition: EvtCyclic3.cc:195

The documentation for this class was generated from the following files: