EvtFDC.cc

Go to the documentation of this file.

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of models developed at BES collaboration
//      based on the EvtGen framework.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/BesCopyright
//      Copyright (A) 2006      Ping Rong-Gang @IHEP
//
// Module:  EvtFDC.cc
//
// Description: Model provided by user, see the mannual     
//
// Modification history:
//
//    Ping R.-G.       December, 2006       Module created
//
//------------------------------------------------------------------------
 
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtDiracParticle.hh"
#include "EvtGenBase/EvtScalarParticle.hh"
#include "EvtGenBase/EvtVectorParticle.hh"
#include "EvtGenBase/EvtTensorParticle.hh"
#include "EvtGenBase/EvtPhotonParticle.hh" 
#include "EvtGenBase/EvtNeutrinoParticle.hh"
#include "EvtGenBase/EvtStringParticle.hh"
#include "EvtGenBase/EvtRaritaSchwingerParticle.hh"
#include "EvtGenBase/EvtHighSpinParticle.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenModels/EvtFDC.hh"
#include "EvtGenBase/EvtRandom.hh"
#include <string>
#include <cstdlib>
using namespace std;  //::endl;
 
extern "C"{
 extern void setamp_(int *);
}
 
extern "C" {
 extern void init_fit_();
}
 
extern "C"{
 extern void initevtgen_();
}
 
extern "C" {
 extern double fdcevtgen_( double*,double*,double*, double*); 
}
 
EvtFDC::~EvtFDC() {}
 
void EvtFDC::getName(std::string& model_name){
 
  model_name="FDC";     
 
}
 
EvtDecayBase* EvtFDC::clone(){
 
  return new EvtFDC;
 
}
 
 
void EvtFDC::init(){
   init_fit_();
   initevtgen_();
   int ndiag = -1; //Feyman diagram number, started from 1, if -1, generate all intermediate states including interference effects
  if (getNArg()==1) {
    ndiag = getArg(0);
  }
   setamp_( &ndiag );
 
  EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
}
void EvtFDC::initProbMax(){
  noProbMax();
 
}
 
 
static int nrun=1;
static double max_amps=0.0;
 
void EvtFDC::setEvtMomentum(EvtVector4R *p4){
  for (int i=0; i<_nd; i++){
     e[i]  = _p4[i].get(0); 
     px[i] = _p4[i].get(1);
     py[i] = _p4[i].get(2);
     pz[i] = _p4[i].get(3);
  }
}
 
void EvtFDC::decay( EvtParticle *p ){
double amps,SamAmps,rd1;
// calculated the max amplitude square in 20000 events, is it enough larger?
  if(nrun==1){
     int ir,nd;
     std::cout<<"Wait for a moment for calculation of maxiumal amplitude squared"<<std::endl;
     for(ir=0;ir<=200000;ir++){
  loop0:
       p->initializePhaseSpace(getNDaug(),getDaugs());
       int nd=p->getNDaug();
       _nd=nd;
       for(int i=0;i<=nd-1;i++){
          _p4[i]=p->getDaug(i)->getP4Lab();
       }
       for (int i=0; i<_nd; i++){
          e[i]  = _p4[i].get(0);
          px[i] = _p4[i].get(1);
          py[i] = _p4[i].get(2);
          pz[i] = _p4[i].get(3);   
        }
       //for debuging
       amps = fdcevtgen_(e,px,py,pz); 
       //std::cout<<"nrun "<<nrun<<std::endl;
       if(amps<0) goto loop0;
       if(amps>max_amps) max_amps=amps*1.01;
       nrun++;
     }
  }
if(max_amps==0.0) {report(ERROR,"EvtGen")<<"The decay amplitude square should be positive number, but it is "<<max_amps<<endl;abort();}
//cout<<"max_amp="<<max_amps<<endl;
 
 
loop:
  p->initializePhaseSpace(getNDaug(),getDaugs());
  int i;
  for(i=0;i<=p->getNDaug()-1;i++){
       _p4[i]=p->getDaug(i)->getP4Lab();
      }
   setEvtMomentum(_p4);
//  Put phase space results into the daughters.
   amps = fdcevtgen_(e,px,py,pz);
   if(amps < 0) goto loop;
   SamAmps=amps/max_amps;
   rd1=EvtRandom::Flat(0.0, 1.0);
  if(rd1>=SamAmps) goto loop;
  if(amps>max_amps) max_amps=amps*1.01;
   nrun++; 
}