STARLIGHT/starlight/src/starlightdpmjet.cpp

   1 /*
   2     <one line to give the program's name and a brief idea of what it does.>
   3     Copyright (C) <year>  <name of author>
   4
   5     This program is free software; you can redistribute it and/or modify
   6     it under the terms of the GNU General Public License as published by
   7     the Free Software Foundation; either version 2 of the License, or
   8     (at your option) any later version.
   9
  10     This program is distributed in the hope that it will be useful,
  11     but WITHOUT ANY WARRANTY; without even the implied warranty of
  12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13     GNU General Public License for more details.
  14
  15     You should have received a copy of the GNU General Public License along
  16     with this program; if not, write to the Free Software Foundation, Inc.,
  17     51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
  18
  19 */
  20
  21 #include "starlightdpmjet.h"
  22 #include "spectrum.h"
  23 #include <iostream>
  24 #include <spectrumprotonnucleus.h>
  25 #include "starlightconfig.h"
  26
  27
  28 extern "C"
  29 {
  30     extern struct
  31         {
  32
  33             double slpx[200000];
  34             double slpy[200000];
  35             double slpz[200000];
  36             double sle[200000];
  37             double slm[200000];
  38             int slpid[200000];
  39             int slcharge[200000];
  40
  41         } dpmjetparticle_;
  42
  43     void dt_produceevent_(float* gammaE, int* nparticles);
  44     void dt_getparticle_(int *ipart, int *res);
  45     void dt_initialise_();
  46 }
  47
  48 starlightDpmJet::starlightDpmJet(beamBeamSystem& beamsystem ) : eventChannel(beamsystem)
  49         ,_spectrum(0)
  50         ,_doDoubleEvent(true)
  51         ,_minGammaEnergy(6.0)
  52         ,_maxGammaEnergy(600000.0)
  53         ,_protonMode(false)
  54 {
  55
  56 }
  57
  58 int starlightDpmJet::init()
  59 {
  60    if(_protonMode)
  61    {
  62     _spectrum = new spectrumProtonNucleus(&_bbs);
  63    }
  64    else
  65    {
  66     _spectrum = new spectrum(&_bbs);
  67    }
  68
  69    _spectrum->setMinGammaEnergy(_minGammaEnergy);
  70    _spectrum->setMaxGammaEnergy(_maxGammaEnergy);
  71
  72    if(!_doDoubleEvent)
  73    {
  74     _spectrum->generateKsingle();
  75    }
  76    else
  77    {
  78     _spectrum->generateKdouble();
  79    }
  80
  81    return 0;
  82
  83 }
  84
  85
  86 upcEvent starlightDpmJet::produceEvent()
  87 {
  88
  89     upcEvent event;
  90
  91     if (!_doDoubleEvent)
  92     {
  93       //int zdirection = (Randy.Rndom()) < 0.5 ? -1 : 1;
  94       int zdirection = 1;
  95         float gammaE = _spectrum->drawKsingle();
  96         event = produceSingleEvent(zdirection, gammaE);
  97         //        std::cout << "Gamma energy: " << gammaE << std::endl;
  98     }
  99     else
 100     {
 101         event = produceDoubleEvent();
 102     }
 103
 104     return event;
 105 }
 106
 107 upcEvent starlightDpmJet::produceSingleEvent(int zdirection, float gammaE)
 108 {
 109
 110     upcEvent event;
 111     event.addGamma(gammaE);
 112
 113     int nParticles = 0;
 114
 115     dt_produceevent_(&gammaE, &nParticles);
 116
 117
 118     //In which direction do we go?
 119     double rapidity = _bbs.beam1().rapidity()*zdirection;
 120
 121     for (int i = 0; i < nParticles; i++)
 122     {
 123         starlightParticle particle(dpmjetparticle_.slpx[i], dpmjetparticle_.slpy[i], zdirection*dpmjetparticle_.slpz[i], dpmjetparticle_.sle[i], dpmjetparticle_.slm[i], dpmjetparticle_.slpid[i], dpmjetparticle_.slcharge[i]);
 124         vector3 boostVector(0, 0, tanh(-rapidity));
 125         particle.Boost(boostVector);
 126         event.addParticle(particle);
 127     }
 128     return event;
 129 }
 130
 131 upcEvent starlightDpmJet::produceDoubleEvent()
 132 {
 133     upcEvent event;
 134
 135     float gammaE1 = 0.0;
 136     float gammaE2 = 0.0;
 137
 138     _spectrum->drawKdouble(gammaE1, gammaE2);
 139
 140 //    std::cout << "Gamma1 energy: " << gammaE1 << std::endl;
 141     //std::cout << "Gamma2 energy: " << gammaE2 << std::endl;
 142
 143     //In which direction do we go?
 144     int zdirection = (randyInstance.Rndom()) < 0.5 ? -1 : 1;
 145
 146     event = produceSingleEvent(zdirection, gammaE1);
 147
 148     zdirection = zdirection *-1;
 149
 150     event = event + produceSingleEvent(zdirection, gammaE2);
 151
 152     return event;
 153 }
 154
 155