Commit a2468719 authored by Leif Lonnblad's avatar Leif Lonnblad
Browse files

This commit was manufactured by cvs2svn to create branch 'CLHEP_1_900'.

parent 5e8697e7
// -*- C++ -*-
// $Id: PuncturedSmearedExp.hh,v 1.1 2004/02/05 15:35:25 boudreau Exp $
// ------------------------------------------------------------------------------//
// This function-object makes an exponential with acceptance holes ("punctures") //
// smeared by a resolution function. //
// //
// //
// Joe Boudreau. //
// //
// ------------------------------------------------------------------------------//
#ifndef _PuncturedSmearedExp_h_
#define _PuncturedSmearedExp_h_
#include "CLHEP/GenericFunctions/AbsFunction.hh"
#include "CLHEP/GenericFunctions/Parameter.hh"
namespace Genfun {
/**
* @author
* @ingroup genfun
*/
class PuncturedSmearedExp: public AbsFunction {
FUNCTION_OBJECT_DEF(PuncturedSmearedExp)
public:
// Constructor
PuncturedSmearedExp();
// Copy constructor
PuncturedSmearedExp(const PuncturedSmearedExp &right);
// Destructor:
virtual ~PuncturedSmearedExp();
// Retreive function value
virtual double operator ()(double argument) const;
virtual double operator ()(const Argument & a) const {return operator() (a[0]);}
// Lifetime of exponential:
Parameter & lifetime();
const Parameter & lifetime() const;
// Width of the gaussian:
Parameter & sigma();
const Parameter & sigma() const;
// Puncture this thing:
void puncture(double min, double max);
// Get the puncture parameters:
Parameter & min(unsigned int i);
Parameter & max(unsigned int i);
const Parameter & min(unsigned int i) const;
const Parameter & max(unsigned int i) const;
private:
// These are for calculating mixing terms.
double pow(double x, int n) const ;
double erfc(double x) const ;
// It is illegal to assign an adjustable constant
const PuncturedSmearedExp & operator=(const PuncturedSmearedExp &right);
Parameter _lifetime;
Parameter _sigma;
std::vector<Parameter> _punctures;
};
} // namespace Genfun
#endif
// -*- C++ -*-
// $Id: PuncturedSmearedExp.cc,v 1.2 2004/04/20 15:03:52 pfeiffer Exp $
#include "CLHEP/GenericFunctions/PuncturedSmearedExp.hh"
#include <sstream>
#include <cmath>
namespace Genfun {
FUNCTION_OBJECT_IMP(PuncturedSmearedExp)
PuncturedSmearedExp::PuncturedSmearedExp() :
_lifetime ("Lifetime", 1.0, 0.0), // Bounded from below by zero, by default
_sigma ("Sigma", 1.0, 0.0) // Bounded from below by zero, by default
{
}
PuncturedSmearedExp::PuncturedSmearedExp(const PuncturedSmearedExp & right):
_lifetime (right._lifetime),
_sigma (right._sigma),
_punctures (right._punctures)
{
}
void PuncturedSmearedExp::puncture(double min, double max) {
std::ostringstream mn, mx;
mn << "Min_" << _punctures.size()/2;
mx << "Max_" << _punctures.size()/2;
_punctures.push_back(Parameter(mn.str(), min, 0.0, 10.0));
_punctures.push_back(Parameter(mx.str(), max, 0.0, 10.0));
}
Parameter & PuncturedSmearedExp::min(unsigned int i) {
return _punctures[2*i];
}
const Parameter & PuncturedSmearedExp::min(unsigned int i) const {
return _punctures[2*i];
}
Parameter & PuncturedSmearedExp::max(unsigned int i) {
return _punctures[2*i+1];
}
const Parameter & PuncturedSmearedExp::max(unsigned int i) const {
return _punctures[2*i+1];
}
PuncturedSmearedExp::~PuncturedSmearedExp() {
}
Parameter & PuncturedSmearedExp::lifetime() {
return _lifetime;
}
const Parameter & PuncturedSmearedExp::lifetime() const {
return _lifetime;
}
Parameter & PuncturedSmearedExp::sigma() {
return _sigma;
}
const Parameter & PuncturedSmearedExp::sigma() const {
return _sigma;
}
double PuncturedSmearedExp::operator() (double argument) const {
// Fetch the paramters. This operator does not convolve numerically.
static const double sqrtTwo = sqrt(2.0);
double sigma = _sigma.getValue();
double tau = _lifetime.getValue();
double x = argument;
// Copy:
std::vector<double> punctures(_punctures.size());
for (size_t i=0;i<_punctures.size();i++) punctures[i]=_punctures[i].getValue();
// Overlap elimination:
bool overlap=true;
while (overlap) {
overlap=false;
for (size_t i=0;i<punctures.size()/2;i++) {
std::sort(punctures.begin()+2*i, punctures.begin()+2*i+2);
double min1=punctures[2*i];
double max1=punctures[2*i+1];
for (size_t j=i+1;j<punctures.size()/2;j++) {
std::sort(punctures.begin()+2*j, punctures.begin()+2*j+2);
double min2=punctures[2*j];
double max2=punctures[2*j+1];
if ((min2>min1 && max1>min2) || (min1>min2 && max2<min1)) {
punctures[2*i] =std::min(min1,min2);
punctures[2*i+1]=std::max(max1,max2);
std::vector<double>::iterator t0=punctures.begin()+2*j, t1=t0+2;
punctures.erase(t0,t1);
overlap=true;
break;
}
}
if (overlap) break;
}
}
double expG=0,norm=0;
for (size_t i=0;i<punctures.size()/2;i++) {
double a = punctures[2*i];
double b = punctures[2*i+1];
double alpha = (a/sigma + sigma/tau)/sqrtTwo;
double beta = (b/sigma + sigma/tau)/sqrtTwo;
double delta = 1/sqrtTwo/sigma;
norm += 2*tau*exp(1/(4*delta*delta*tau*tau))*(exp(-alpha/(delta*tau)) - exp(-beta/(delta*tau)));
expG += ((erfc(alpha-delta*x) - erfc(beta-delta*x))*exp(-x/tau) );
}
// protection:
if (norm==0) return norm;
return expG/norm;
}
double PuncturedSmearedExp::erfc(double x) const {
// This is accurate to 7 places.
// See Numerical Recipes P 221
double t, z, ans;
z = (x < 0) ? -x : x;
t = 1.0/(1.0+.5*z);
ans = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
t*(-0.82215223+t*0.17087277 ))) ))) )));
if ( x < 0 ) ans = 2.0 - ans;
return ans;
}
double PuncturedSmearedExp::pow(double x,int n) const {
double val=1;
for(int i=0;i<n;i++){
val=val*x;
}
return val;
}
} // namespace Genfun
#ifndef PARTICLENAME_HH
#define PARTICLENAME_HH
// $Id: ParticleName.hh,v 1.1 2004/04/14 23:56:27 garren Exp $
// ----------------------------------------------------------------------
//
// ParticleName.hh
// Author: Lynn Garren and Walter Brown
//
// Create a map that gives a standard name for each pre-defined
// particle ID number. This map is initialized if and only if
// the public functions are called. Because the map is static,
// the initialization happens only once.
//
//
// ----------------------------------------------------------------------
#include <string>
#include <map>
#include <iostream>
namespace HepPDT {
// get a known HepPDT Particle name
std::string particleName( const int );
// list all known names
void listHepPDTParticleNames( std::ostream & os );
// verify that this number has a valid name
bool validParticleName( const int );
// access the ParticleNameMap for other purposes
typedef std::map< int, std::string > ParticleNameMap;
ParticleNameMap const & getParticleNameMap();
} // namespace HepPDT
#endif // PARTICLENAME_HH
// $Id: ParticleTranslation.hh,v 1.1 2004/04/14 23:56:27 garren Exp $
// ----------------------------------------------------------------------
//
// ParticleTranslation.hh
// Author: Lynn Garren
//
// Use this class to maintain translation info
// This class is used when testing translations
// ----------------------------------------------------------------------
#ifndef PARTICLETRANSLATION_HH
#define PARTICLETRANSLATION_HH
#include <string>
#include <iostream>
#include "CLHEP/HepPDT/ParticleID.hh"
namespace HepPDT {
/**
* @author
* @ingroup heppdt
*/
class ParticleTranslation {
public:
// --- birth/death:
//
inline ParticleTranslation( ParticleID, int, std::string &, std::string & mc );
inline ParticleTranslation( );
inline ~ParticleTranslation();
// --- copying:
//
inline void swap ( ParticleTranslation & rhs );
inline ParticleTranslation( const ParticleTranslation & orig );
inline ParticleTranslation & operator = ( const ParticleTranslation & rhs );
// --- accessors:
//
const std::string name() const { return itsPID.PDTname(); }
const std::string & originalName() const { return itsOriginalName; }
const std::string & monteCarloName() const { return itsMonteCarloName; }
ParticleID ID() const { return itsPID; }
int pid() const { return itsPID.pid(); }
int oid() const { return itsOriginalID; }
inline void write( std::ostream & os ) const;
private:
// --- class-specific data:
//
ParticleID itsPID;
int itsOriginalID;
std::string itsOriginalName;
std::string itsMonteCarloName;
}; // ParticleTranslation
inline void swap( ParticleTranslation & first, ParticleTranslation & second ) {
first.swap( second );
}
} // HepPDT
#include "CLHEP/HepPDT/ParticleTranslation.icc"
#endif // PARTICLETRANSLATION_HH
// $Id: ParticleTranslation.icc,v 1.1 2004/04/14 23:56:27 garren Exp $
// ----------------------------------------------------------------------
//
// ParticleTranslation.icc
// Author: Lynn Garren
//
// ----------------------------------------------------------------------
#include <algorithm> // swap()
namespace HepPDT {
ParticleTranslation::ParticleTranslation(ParticleID pid, int oid,
std::string & onm, std::string & mc )
: itsPID ( pid ),
itsOriginalID ( oid ),
itsOriginalName ( onm ),
itsMonteCarloName ( mc )
{ ; }
ParticleTranslation::ParticleTranslation( )
: itsPID ( ParticleID(0) ),
itsOriginalID ( 0 ),
itsOriginalName ( "" ),
itsMonteCarloName ( "" )
{ ; }
ParticleTranslation::~ParticleTranslation()
{ ; }
void ParticleTranslation::swap( ParticleTranslation & other )
{
itsPID.swap( other.itsPID );
std::swap(itsOriginalID, other.itsOriginalID);
std::swap(itsOriginalName, other.itsOriginalName);
std::swap(itsMonteCarloName , other.itsMonteCarloName);
}
ParticleTranslation::ParticleTranslation( const ParticleTranslation & orig )
: itsPID ( orig.itsPID ),
itsOriginalID ( orig.itsOriginalID ),
itsOriginalName ( orig.itsOriginalName ),
itsMonteCarloName ( orig.itsMonteCarloName )
{ ; }
ParticleTranslation & ParticleTranslation::operator=( const ParticleTranslation & rhs )
{
ParticleTranslation temp( rhs );
swap( temp );
return *this;
}
void ParticleTranslation::write( std::ostream & os ) const
{
os << " " << itsMonteCarloName << ": ";
os.width(20);
os << itsOriginalName;
os.width(12);
os << itsOriginalID;
os << " HepPDT: ";
os.width(20);
os << name();
os.width(12);
os << pid() << std::endl;
return;
}
} // HepPDT
#ifndef TRANSLATIONLIST_HH
#define TRANSLATIONLIST_HH
// $Id: TranslationList.hh,v 1.1 2004/04/14 23:56:27 garren Exp $
// ----------------------------------------------------------------------
//
// TranslationList.hh
// Author: Lynn Garren
//
// list translations for various MonteCarlo input types
// DO NOT mix these functions with the addXXXParticles functions
// These functions will read a table file and write a translation list
//
//
// ----------------------------------------------------------------------
#include <iostream>
namespace HepPDT {
/**
* @author
* @ingroup heppdt
*/
// --- free functions
//
bool listPDGTranslation ( std::istream &, std::ostream & );
bool listPythiaTranslation ( std::istream &, std::ostream & );
} // namespace HepPDT
#endif // TRANSLATIONLIST_HH
// $Id: examListHerwig.cc,v 1.1 2004/04/14 23:56:27 garren Exp $
// -------------------------------------------------------------------
//
// List the herwig translation
// This requires access to the herwig library
//
// main program must be in C++
//
// -------------------------------------------------------------------
#include <fstream>
#include <iostream>
#include "CLHEP/HepPDT/ParticleName.hh"
#include "CLHEP/HepPDT/ParticleID.hh"
extern "C" {
void list_herwig_init_ ( int * nevt );
void list_herwig_end_ ( );
void get_list_size_ ( int * );
void get_herwig_name_( int * ihwg, int * id, char *name );
}
int main()
{
int nevt=20;
int idir=1;
int i, j, iend, isize;
int hid, id;
char cname[10];
std::string hname;
std::string pn;
static char outfile[] = "examListHerwig.out";
std::string title = "HepPDT listing of Herwig translations";
// initialize herwig
list_herwig_init_ ( & nevt );
// open the output stream
std::ofstream os( outfile );
if( !os ) {
std::cout << "error opening output file" << std::endl;
exit(1);
}
get_list_size_( & isize );
os << " " << title << std::endl;
os << " number of Herwig particles: " << isize << std::endl;
for( i=1, iend=isize+1; i<isize; ++i ) {
// get info from herwig
for( j=0; j<10; ++j) { cname[j] = '\0'; }
get_herwig_name_( & i, & hid, cname );
hname = std::string( cname );
id = HepPDT::translateHerwigtoPDT( hid );
pn = HepPDT::particleName( id );
os << "Herwig: ";
os.width(7);
os << i ;
os.width(12);
os << hid << " " << hname;
os << " HepPDT: " ;
os.width(12);
os << id << " " << pn << std::endl;
}
list_herwig_end_();
return 0;
}
This diff is collapsed.
subroutine list_herwig_init(nevt)
c
c initialization for the herwig C++ listing
c
#include "herwig65.inc"
integer lnhwrt,lnhrd,lnhout,lnhdcy
common/heplun/lnhwrt,lnhrd,lnhout,lnhdcy
external hwudat
integer n
integer istr,nevt
C
C initialize HEP logical units
lnhwrt=0
lnhrd=0
lnhdcy=0
lnhout=22
lhwout=lnhout
C open(unit=lnhout,file='examHerwigToStdHep.lpt',status='new')
C
c call hptrlsth
C
return
end
subroutine get_list_size( isize )
c return the maximum size of herwig's particle list
#include "herwig65.inc"
integer isize
isize = NRES
return
end
subroutine get_herwig_name( ihwg, id, name )
c ihwg is the index into herwig's short list
#include "herwig65.inc"
integer id, ihwg
character*8 name
id = 0
call HWUIDT(2,id,ihwg,name)
return
end
subroutine list_herwig_end
integer lnhwrt,lnhrd,lnhout,lnhdcy
common/heplun/lnhwrt,lnhrd,lnhout,lnhdcy
C---terminate elementary process
c call hwefin
C close(unit=lnhout)
return
end
C----------------------------------------------------------------------
subroutine hwabeg
C... user's routine for initialization
end
subroutine hwaend
C... user's routine for terminal calculations, histogram output, etc
end
subroutine hwanal
C... user's routine to analyse data from event
end
C----------------------------------------------------------------------
// $Id: examListIsajet.cc,v 1.1 2004/04/14 23:56:27 garren Exp $
// -------------------------------------------------------------------
//
// List the isajet translation
// This requires access to the isajet library
//
// main program must be in C++
//
// -------------------------------------------------------------------
#include <fstream>
#include <iostream>
#include "CLHEP/HepPDT/ParticleName.hh"
#include "CLHEP/HepPDT/ParticleID.hh"
extern "C" {
void list_isajet_init_ ( );
void get_isajet_name_( int * i, int * id, int * aid, char *name, char *aname );