File indexing completed on 2024-09-27 07:03:42
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <cstdio>
0021 #include <iostream>
0022 #include "starlightpythia.h"
0023 #include "pythiaInterface.h"
0024 #include "spectrumprotonnucleus.h"
0025 #include "eXevent.h"
0026 #include <cmath>
0027 #include <sstream>
0028
0029 starlightPythia::starlightPythia(const inputParameters& inputParametersInstance, beamBeamSystem& bbsystem) : eventChannel(inputParametersInstance, bbsystem)
0030 ,_spectrum(0)
0031 ,_doDoubleEvent(false)
0032 ,_minGammaEnergy(inputParametersInstance.minGammaEnergy())
0033 ,_maxGammaEnergy(inputParametersInstance.maxGammaEnergy())
0034 ,_fullEventRecord(false)
0035 {
0036 }
0037
0038 starlightPythia::~starlightPythia()
0039 {
0040
0041 }
0042
0043 int starlightPythia::init(std::string pythiaParams, bool fullEventRecord)
0044 {
0045 _fullEventRecord = fullEventRecord;
0046 _spectrum = new spectrumProtonNucleus(_randy,&_bbs);
0047
0048 _spectrum->setMinGammaEnergy(_minGammaEnergy);
0049 _spectrum->setMaxGammaEnergy(_maxGammaEnergy);
0050
0051 if(!_doDoubleEvent)
0052 {
0053 _spectrum->generateKsingle();
0054 }
0055 else
0056 {
0057 _spectrum->generateKdouble();
0058 }
0059
0060
0061 pythiaInterface::pygive("mstp(171)=1");
0062 pythiaInterface::pygive("mstp(172)=1");
0063
0064 std::stringstream ss(pythiaParams);
0065 std::string p;
0066 while(std::getline(ss, p, ';'))
0067 {
0068 if(p.size()>1)
0069 {
0070 pythiaInterface::pygive(p.c_str());
0071 }
0072 }
0073
0074 pythiaInterface::pyinit("FIXT", "gamma", "p", _maxGammaEnergy);
0075
0076 return 0;
0077 }
0078
0079 eXEvent starlightPythia::produceEvent()
0080 {
0081 eXEvent event;
0082
0083 if (!_doDoubleEvent)
0084 {
0085 double gammaE = 0;
0086 do
0087 {
0088 gammaE = _spectrum->drawKsingle();
0089 } while(isnan(gammaE));
0090 event.addGamma(gammaE);
0091
0092 char opt[32];
0093 std::sprintf(opt, "parp(171)=%f", gammaE/_maxGammaEnergy);
0094 pythiaInterface::pygive(opt);
0095 pythiaInterface::pyevnt();
0096
0097 int zdirection = (_bbs.beam1().Z()==1 ? 1 : -1);
0098 double boost = acosh(_bbs.beamLorentzGamma())*zdirection;
0099
0100 vector3 boostVector(0, 0, tanh(boost));
0101
0102 for(int idx = 0; idx < pyjets_.n; idx++)
0103 {
0104 if(pyjets_.k[0][idx] > 10 && _fullEventRecord==false) continue;
0105 int pdgCode = pyjets_.k[1][idx];
0106 int charge = 0;
0107 if( pdgCode == 12 || pdgCode == 14 || pdgCode == 16 || pdgCode == 22 || pdgCode == 111 || pdgCode == 130 || pdgCode == 321 || pdgCode == 2112)
0108 {
0109 charge = 0;
0110 }
0111 else
0112 {
0113 charge = (pdgCode > 0) - (pdgCode < 0);
0114 }
0115
0116 starlightParticle particle(pyjets_.p[0][idx], pyjets_.p[1][idx], -zdirection*pyjets_.p[2][idx], pyjets_.p[3][idx], pyjets_.p[4][idx], pyjets_.k[1][idx], charge);
0117 if(_fullEventRecord)
0118 {
0119 particle.setFirstDaughter(pyjets_.k[3][idx]);
0120 particle.setLastDaughter(pyjets_.k[4][idx]);
0121 particle.setStatus(pyjets_.k[0][idx]);
0122 }
0123 particle.Boost(boostVector);
0124 event.addParticle(particle);
0125 }
0126 }
0127 return event;
0128 }