#include #include #include #include "PhotosRandom.h" #include "PhotosEvent.h" #include "Photos.h" #include "Log.h" using std::vector; using std::cout; using std::endl; using std::ios_base; namespace Photospp { Photos Photos::_instance; vector* > *Photos::supBremList = 0; vector* > *Photos::forceBremList = 0; vector* > *Photos::forceMassList = 0; vector *Photos::ignoreStatusCodeList = 0; bool Photos::isSuppressed=false; bool Photos::massFrom4Vector=true; double Photos::momentum_conservation_threshold = 0.1; bool Photos::meCorrectionWtForZ=false; bool Photos::meCorrectionWtForW=false; bool Photos::meCorrectionWtForScalar=false; bool Photos::isCreateHistoryEntries=false; int Photos::historyEntriesStatus = 3; double (*Photos::randomDouble)() = PhotosRandom::randomReal; Photos::Photos() { setAlphaQED (0.00729735039); setInfraredCutOff (0.01); setInterference (true); setDoubleBrem (true); setQuatroBrem (false); setTopProcessRadiation(true); setCorrectionWtForW (true); // setExponentiation(true) moved to initialize() due to communication // problems with Fortran under MacOS. phokey_.iexp = 1; } void Photos::initialize() { // Should return if already initialized? /******************************************************************************* All the following parameter setters can be called after PHOINI. Initialization of kinematic correction against rounding errors. The set values will be used later if called with zero. Default parameter is 1 (no correction) optionally 2, 3, 4 In case of exponentiation new version 5 is needed in most cases. Definition given here will be thus overwritten in such a case below in routine PHOCIN *******************************************************************************/ if(!phokey_.iexp) initializeKinematicCorrections(1); else setExponentiation(true); int buf=1; iphqrk_(&buf); // Blocks emission from quarks if buf=1 (default); enables if buf=2 // Physical treatment will be 3, option 2 is not realistic and for tests only buf=2; iphekl_(&buf); // Blocks emission in pi0 to gamma e+ e- if parameter is >1 (enables otherwise) // Initialize status counter for warning messages for(int i=0;i<10;i++) phosta_.status[i]=0; // elementary security level, should remain 1 but we may want to have a method to change. phosta_.ifstop=1; pholun_.phlun=6; // Logical output unit for printing error messages // Define pi and 2*pi phpico_.pi =3.14159265358979324; phpico_.twopi=6.28318530717958648; // Further initialization done automatically // see places with - VARIANT A - VARIANT B - all over to switch between options //----------- SLOWER VARIANT A, but stable ------------------ //--- it is limiting choice for small XPHCUT in fixed orer //--- modes of operation // Best choice is if FINT=2**N where N+1 is maximal number // of charged daughters // see report on overweihted events if(phokey_.interf) maxWtInterference(2.0); else maxWtInterference(1.0); /* ----------- FASTER VARIANT B ------------------ -- it is good for tests of fixed order and small XPHCUT -- but is less promising for more complex cases of interference -- sometimes fails because of that if(phokey_.interf) maxWtInterference(1.8); else maxWtInterference(0.0); ------------END VARIANTS A B ----------------------- */ //------------------------------------------------------------------------------ // Print PHOTOS header //------------------------------------------------------------------------------ int coutPrec = cout.precision(6); ios_base::fmtflags flags = cout.setf(ios_base::floatfield); cout< particles = p->getDecayTree(); vector branches = PhotosBranch::createBranches(particles); for(int i=0;i<(int)branches.size();i++) branches.at(i)->process(); } void Photos::suppressBremForDecay(int count, int motherID, ... ) { va_list arg; va_start(arg, motherID); vector *v = new vector(); v->push_back(motherID); for(int i = 0;ipush_back(va_arg(arg,int)); } va_end(arg); v->push_back(0); if(!supBremList) supBremList = new vector< vector* >(); supBremList->push_back(v); } void Photos::suppressBremForBranch(int count, int motherID, ... ) { va_list arg; va_start(arg, motherID); vector *v = new vector(); v->push_back(motherID); for(int i = 0;ipush_back(va_arg(arg,int)); } va_end(arg); v->push_back(1); if(!supBremList) supBremList = new vector< vector* >(); supBremList->push_back(v); } void Photos::forceBremForDecay(int count, int motherID, ... ) { va_list arg; va_start(arg, motherID); vector *v = new vector(); v->push_back(motherID); for(int i = 0;ipush_back(va_arg(arg,int)); } va_end(arg); v->push_back(0); if(!forceBremList) forceBremList = new vector< vector* >(); forceBremList->push_back(v); } void Photos::forceBremForBranch(int count, int motherID, ... ) { va_list arg; va_start(arg, motherID); vector *v = new vector(); v->push_back(motherID); for(int i = 0;ipush_back(va_arg(arg,int)); } va_end(arg); v->push_back(1); if(!forceBremList) forceBremList = new vector< vector* >(); forceBremList->push_back(v); } void Photos::createHistoryEntries(bool flag, int status) { if(status<3) { Log::Warning()<<"Photos::createHistoryEntries: status must be >=3"<=3"<(); // Do not add duplicate entries to the list for(unsigned int i=0;isize();i++) if( status==ignoreStatusCodeList->at(i) ) return; ignoreStatusCodeList->push_back(status); } void Photos::deIgnoreParticlesOfStatus(int status) { if(!ignoreStatusCodeList) return; for(unsigned int i=0;isize();i++) { if( status==ignoreStatusCodeList->at(i) ) { ignoreStatusCodeList->erase(ignoreStatusCodeList->begin()+i); return; } } } bool Photos::isStatusCodeIgnored(int status) { if(!ignoreStatusCodeList) return false; for(unsigned int i=0;isize();i++) if( status==ignoreStatusCodeList->at(i) ) return true; return false; } void Photos::setRandomGenerator( double (*gen)() ) { if(gen==NULL) randomDouble = PhotosRandom::randomReal; else randomDouble = gen; } void Photos::setExponentiation(bool expo) { phokey_.iexp = (int) expo; if(expo) { setDoubleBrem(false); setQuatroBrem(false); setInfraredCutOff(0.0000001); initializeKinematicCorrections(5); phokey_.expeps=0.0001; } } void Photos::setMeCorrectionWtForW(bool corr) { meCorrectionWtForW=corr; } void Photos::setMeCorrectionWtForZ(bool corr) { meCorrectionWtForZ=corr; } void Photos::setMeCorrectionWtForScalar(bool corr) { meCorrectionWtForScalar=corr; } void Photos::setStopAtCriticalError(bool stop) { phosta_.ifstop=(int)stop; if(!stop) { Log::Info()<<"PHOTOS production mode. Elementary test of data flow from event record disabled. "<* >(); forceMassList->push_back( new pair(pdgid, -1.0) ); } void Photos::forceMass(int pdgid, double mass) { if(mass<0.0) { Log::Warning()<<"Photos::forceMass: Mass must be > 0.0"<* >(); forceMassList->push_back( new pair(pdgid, mass) ); } } // namespace Photospp