estarlight/src/starlightpythia.cpp

0001 /*
0002     <one line to give the library's name and an idea of what it does.>
0003     Copyright (C) 2011  Oystein Djuvsland <oystein.djuvsland@gmail.com>
0004
0005     This library is free software; you can redistribute it and/or
0006     modify it under the terms of the GNU Lesser General Public
0007     License as published by the Free Software Foundation; either
0008     version 2.1 of the License, or (at your option) any later version.
0009
0010     This library is distributed in the hope that it will be useful,
0011     but WITHOUT ANY WARRANTY; without even the implied warranty of
0012     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
0013     Lesser General Public License for more details.
0014
0015     You should have received a copy of the GNU Lesser General Public
0016     License along with this library; if not, write to the Free Software
0017     Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
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     // Set to run with varying energies
0061     pythiaInterface::pygive("mstp(171)=1"); // Varying energies
0062     pythiaInterface::pygive("mstp(172)=1"); // Set the energy before generation
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); // Fixed target, beam, target, beam momentum (GeV/c)
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); // Set the energy of the photon beam (gammaE/1000 * 1000.0);
0095       pythiaInterface::pyevnt(); // Generate event
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 }