Line data Source code
1 : #ifdef EVTGEN_PHOTOS
2 : //--------------------------------------------------------------------------
3 : //
4 : // Environment:
5 : // This software is part of the EvtGen package. If you use all or part
6 : // of it, please give an appropriate acknowledgement.
7 : //
8 : // Copyright Information: See EvtGen/COPYRIGHT
9 : // Copyright (C) 2011 University of Warwick, UK
10 : //
11 : // Module: EvtPhotosEngine
12 : //
13 : // Description: Interface to the PHOTOS external generator
14 : //
15 : // Modification history:
16 : //
17 : // John Back May 2011 Module created
18 : //
19 : //------------------------------------------------------------------------
20 :
21 : #include "EvtGenExternal/EvtPhotosEngine.hh"
22 :
23 : #include "EvtGenBase/EvtPDL.hh"
24 : #include "EvtGenBase/EvtVector4R.hh"
25 : #include "EvtGenBase/EvtPhotonParticle.hh"
26 : #include "EvtGenBase/EvtReport.hh"
27 : #include "EvtGenBase/EvtRandom.hh"
28 :
29 : #include "Photos.h"
30 : #include "PhotosHepMCEvent.h"
31 : #include "PhotosHepMCParticle.h"
32 : #include "PhotosParticle.h"
33 :
34 : #include "HepMC/GenVertex.h"
35 : #include "HepMC/SimpleVector.h"
36 : #include "HepMC/Units.h"
37 :
38 : #include <iostream>
39 : #include <sstream>
40 : #include <vector>
41 :
42 : using std::endl;
43 :
44 0 : EvtPhotosEngine::EvtPhotosEngine(std::string photonType, bool useEvtGenRandom) {
45 :
46 0 : _photonType = photonType;
47 0 : _gammaId = EvtId(-1,-1);
48 0 : _mPhoton = 0.0;
49 :
50 0 : report(Severity::Info,"EvtGen")<<"Setting up PHOTOS."<<endl;
51 :
52 0 : if (useEvtGenRandom == true) {
53 :
54 0 : report(Severity::Info,"EvtGen")<<"Using EvtGen random number engine also for Photos++"<<endl;
55 :
56 0 : Photospp::Photos::setRandomGenerator(EvtRandom::Flat);
57 :
58 : }
59 :
60 0 : Photospp::Photos::initialize();
61 :
62 : // Set minimum photon energy (0.1 keV at 1 GeV scale)
63 0 : Photospp::Photos::setInfraredCutOff(1.0e-7);
64 : // Increase the maximum possible value of the interference weight
65 0 : Photospp::Photos::maxWtInterference(64.0); // 2^n, where n = number of charges (+,-)
66 0 : Photospp::Photos::setInterference(true);
67 0 : Photospp::Photos::setExponentiation(true);
68 :
69 0 : _initialised = false;
70 :
71 0 : }
72 :
73 0 : EvtPhotosEngine::~EvtPhotosEngine() {
74 :
75 0 : }
76 :
77 : void EvtPhotosEngine::initialise() {
78 :
79 0 : if (_initialised == false) {
80 :
81 0 : _gammaId = EvtPDL::getId(_photonType);
82 :
83 0 : if (_gammaId == EvtId(-1,-1)) {
84 0 : report(Severity::Info,"EvtGen")<<"Error in EvtPhotosEngine. Do not recognise the photon type "
85 0 : <<_photonType<<". Setting this to \"gamma\". "<<endl;
86 0 : _gammaId = EvtPDL::getId("gamma");
87 0 : }
88 :
89 0 : _mPhoton = EvtPDL::getMeanMass(_gammaId);
90 :
91 0 : _initialised = true;
92 :
93 0 : }
94 :
95 0 : }
96 :
97 : bool EvtPhotosEngine::doDecay(EvtParticle* theMother) {
98 :
99 0 : if (_initialised == false) {this->initialise();}
100 :
101 0 : if (theMother == 0) {return false;}
102 :
103 : // Create a dummy HepMC GenEvent containing a single vertex, with the mother
104 : // assigned as the incoming particle and its daughters as outgoing particles.
105 : // We then pass this event to Photos for processing.
106 : // It will return a modified version of the event, updating the momentum of
107 : // the original particles and will contain any new photon particles.
108 : // We add these extra photons to the mother particle daughter list.
109 :
110 : // Skip running Photos if the particle has no daughters, since we can't add FSR.
111 : // Also skip Photos if the particle has too many daughters (>= 10) to avoid a problem
112 : // with a hard coded upper limit in the PHOENE subroutine.
113 0 : int nDaug(theMother->getNDaug());
114 0 : if (nDaug == 0 || nDaug >= 10) {return false;}
115 :
116 : // Create the dummy event.
117 0 : HepMC::GenEvent* theEvent = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);
118 :
119 : // Create the decay "vertex".
120 0 : HepMC::GenVertex* theVertex = new HepMC::GenVertex();
121 0 : theEvent->add_vertex(theVertex);
122 :
123 : // Add the mother particle as the incoming particle to the vertex.
124 0 : HepMC::GenParticle* hepMCMother = this->createGenParticle(theMother, true);
125 0 : theVertex->add_particle_in(hepMCMother);
126 :
127 : // Find all daughter particles and assign them as outgoing particles to the vertex.
128 : int iDaug(0);
129 0 : for (iDaug = 0; iDaug < nDaug; iDaug++) {
130 :
131 0 : EvtParticle* theDaughter = theMother->getDaug(iDaug);
132 0 : HepMC::GenParticle* hepMCDaughter = this->createGenParticle(theDaughter, false);
133 0 : theVertex->add_particle_out(hepMCDaughter);
134 :
135 : }
136 :
137 : // Now pass the event to Photos for processing
138 : // Create a Photos event object
139 0 : Photospp::PhotosHepMCEvent photosEvent(theEvent);
140 :
141 : // Run the Photos algorithm
142 0 : photosEvent.process();
143 :
144 : // See if Photos has created new particles. If not, do nothing extra.
145 0 : int nDecayPart = theVertex->particles_out_size();
146 : int iLoop(0);
147 :
148 0 : if (nDecayPart > nDaug) {
149 :
150 : // We have extra particles from Photos
151 :
152 : // Get the iterator of outgoing particles for this vertex
153 0 : HepMC::GenVertex::particles_out_const_iterator outIter;
154 0 : for (outIter = theVertex->particles_out_const_begin();
155 0 : outIter != theVertex->particles_out_const_end(); ++outIter) {
156 :
157 : // Get the next HepMC GenParticle
158 0 : HepMC::GenParticle *outParticle = *outIter;
159 :
160 : // Get the three-momentum Photos result for this particle
161 : double px(0.0), py(0.0), pz(0.0);
162 0 : if (outParticle != 0) {
163 0 : HepMC::FourVector HepMCP4 = outParticle->momentum();
164 0 : px = HepMCP4.px();
165 0 : py = HepMCP4.py();
166 0 : pz = HepMCP4.pz();
167 0 : }
168 :
169 : // Create an empty 4-momentum vector for the new/modified daughters
170 0 : EvtVector4R newP4;
171 :
172 0 : if (iLoop < nDaug) {
173 :
174 : // Original daughters
175 0 : EvtParticle* daugParticle = theMother->getDaug(iLoop);
176 0 : if (daugParticle != 0) {
177 :
178 : // Keep the original particle mass, but set the three-momentum
179 : // according to what Photos has modified. However, this will
180 : // violate energy conservation (from what Photos has provided).
181 0 : double mass = daugParticle->mass();
182 0 : double energy = sqrt(mass*mass + px*px + py*py + pz*pz);
183 0 : newP4.set(energy, px, py, pz);
184 : // Set the new four-momentum (FSR applied)
185 0 : daugParticle->setP4WithFSR(newP4);
186 :
187 0 : }
188 :
189 0 : } else {
190 :
191 : // Extra photon particle. Setup the four-momentum object.
192 0 : double energy = sqrt(_mPhoton*_mPhoton + px*px + py*py + pz*pz);
193 0 : newP4.set(energy, px, py, pz);
194 :
195 : // Create a new photon particle and add it to the list of daughters
196 0 : EvtPhotonParticle* gamma = new EvtPhotonParticle();
197 0 : gamma->init(_gammaId, newP4);
198 0 : gamma->setFSRP4toZero();
199 0 : gamma->addDaug(theMother); // Let the mother know about this new particle
200 :
201 : }
202 :
203 : // Increment the loop counter for detecting additional photon particles
204 0 : iLoop++;
205 :
206 0 : }
207 :
208 0 : }
209 :
210 : // Cleanup
211 0 : theEvent->clear();
212 0 : delete theEvent;
213 :
214 : return true;
215 :
216 0 : }
217 :
218 : HepMC::GenParticle* EvtPhotosEngine::createGenParticle(EvtParticle* theParticle, bool incoming) {
219 :
220 : // Method to create an HepMC::GenParticle version of the given EvtParticle.
221 0 : if (theParticle == 0) {return 0;}
222 :
223 : // Get the 4-momentum (E, px, py, pz) for the EvtParticle
224 0 : EvtVector4R p4(0.0, 0.0, 0.0, 0.0);
225 :
226 0 : if (incoming == true) {
227 0 : p4 = theParticle->getP4Restframe();
228 0 : } else {
229 0 : p4 = theParticle->getP4();
230 : }
231 :
232 : // Convert this to the HepMC 4-momentum
233 0 : double E = p4.get(0);
234 0 : double px = p4.get(1);
235 0 : double py = p4.get(2);
236 0 : double pz = p4.get(3);
237 :
238 0 : HepMC::FourVector hepMC_p4(px, py, pz, E);
239 :
240 0 : int PDGInt = EvtPDL::getStdHep(theParticle->getId());
241 :
242 : // Set the status flag for the particle. This is required, otherwise Photos++
243 : // will crash from out-of-bounds array index problems.
244 : int status = Photospp::PhotosParticle::HISTORY;
245 0 : if (incoming == false) {status = Photospp::PhotosParticle::STABLE;}
246 :
247 0 : HepMC::GenParticle* genParticle = new HepMC::GenParticle(hepMC_p4, PDGInt, status);
248 :
249 : return genParticle;
250 :
251 0 : }
252 :
253 : #endif
|
