Line data Source code
1 : #include "HepMC/GenEvent.h"
2 : #include "PhotosHepMCParticle.h"
3 : #include "Log.h"
4 : #include "Photos.h"
5 :
6 : namespace Photospp
7 : {
8 :
9 0 : PhotosHepMCParticle::PhotosHepMCParticle(){
10 0 : m_particle = new HepMC::GenParticle();
11 0 : }
12 :
13 0 : PhotosHepMCParticle::PhotosHepMCParticle(int pdg_id, int status, double mass){
14 0 : m_particle = new HepMC::GenParticle();
15 0 : m_particle->set_pdg_id(pdg_id);
16 0 : m_particle->set_status(status);
17 0 : m_particle->set_generated_mass(mass);
18 0 : }
19 :
20 0 : PhotosHepMCParticle::PhotosHepMCParticle(HepMC::GenParticle * particle){
21 0 : m_particle = particle;
22 0 : }
23 :
24 0 : PhotosHepMCParticle::~PhotosHepMCParticle(){
25 0 : clear(m_mothers);
26 0 : clear(m_daughters);
27 : // clear(m_created_particles);
28 0 : }
29 :
30 :
31 : //delete the TauolaHepMCParticle objects
32 : void PhotosHepMCParticle::clear(std::vector<PhotosParticle*> v){
33 0 : while(v.size()!=0){
34 0 : PhotosParticle * temp = v.back();
35 0 : v.pop_back();
36 0 : delete temp;
37 : }
38 0 : }
39 :
40 : HepMC::GenParticle * PhotosHepMCParticle::getHepMC(){
41 0 : return m_particle;
42 : }
43 :
44 : void PhotosHepMCParticle::setMothers(vector<PhotosParticle*> mothers){
45 :
46 : /******** Deal with mothers ***********/
47 :
48 0 : clear(m_mothers);
49 :
50 : //If there are mothers
51 0 : if(mothers.size()>0){
52 :
53 : HepMC::GenParticle * part;
54 0 : part=dynamic_cast<PhotosHepMCParticle*>(mothers.at(0))->getHepMC();
55 :
56 : //Use end vertex of first mother as production vertex for particle
57 0 : HepMC::GenVertex * production_vertex = part->end_vertex();
58 : HepMC::GenVertex * orig_production_vertex = production_vertex;
59 :
60 0 : if(!production_vertex){ //if it does not exist create it
61 0 : production_vertex = new HepMC::GenVertex();
62 0 : part->parent_event()->add_vertex(production_vertex);
63 0 : }
64 :
65 : //Loop over all mothers to check that the end points to the right place
66 0 : vector<PhotosParticle*>::iterator mother_itr;
67 0 : for(mother_itr = mothers.begin(); mother_itr != mothers.end();
68 0 : mother_itr++){
69 :
70 : HepMC::GenParticle * moth;
71 0 : moth = dynamic_cast<PhotosHepMCParticle*>(*mother_itr)->getHepMC();
72 :
73 0 : if(moth->end_vertex()!=orig_production_vertex)
74 0 : Log::Fatal("PhotosHepMCParticle::setMothers(): Mother production_vertices point to difference places. Can not override. Please delete vertices first.",1);
75 : else
76 0 : production_vertex->add_particle_in(moth);
77 :
78 : //update status info
79 0 : if(moth->status()==PhotosParticle::STABLE)
80 0 : moth->set_status(PhotosParticle::DECAYED);
81 : }
82 0 : production_vertex->add_particle_out(m_particle);
83 0 : }
84 0 : }
85 :
86 :
87 :
88 : void PhotosHepMCParticle::addDaughter(PhotosParticle* daughter){
89 :
90 : //add to this classes internal list as well.
91 0 : m_daughters.push_back(daughter);
92 :
93 : //this assumes there is already an end vertex for the particle
94 :
95 0 : if(!m_particle->end_vertex())
96 0 : Log::Fatal("PhotosHepMCParticle::addDaughter(): This method assumes an end_vertex exists. Maybe you really want to use setDaughters.",2);
97 :
98 0 : HepMC::GenParticle * daugh = (dynamic_cast<PhotosHepMCParticle*>(daughter))->getHepMC();
99 0 : m_particle->end_vertex()->add_particle_out(daugh);
100 :
101 0 : }
102 :
103 : void PhotosHepMCParticle::setDaughters(vector<PhotosParticle*> daughters){
104 :
105 0 : if(!m_particle->parent_event())
106 0 : Log::Fatal("PhotosHepMCParticle::setDaughters(): New particle needs the event set before it's daughters can be added",3);
107 :
108 0 : clear(m_daughters);
109 :
110 : //If there are daughters
111 0 : if(daughters.size()>0){
112 :
113 : //Use production vertex of first daughter as end vertex for particle
114 : HepMC::GenParticle * first_daughter;
115 0 : first_daughter = (dynamic_cast<PhotosHepMCParticle*>(daughters.at(0)))->getHepMC();
116 :
117 : HepMC::GenVertex * end_vertex;
118 0 : end_vertex=first_daughter->production_vertex();
119 : HepMC::GenVertex * orig_end_vertex = end_vertex;
120 :
121 0 : if(!end_vertex){ //if it does not exist create it
122 0 : end_vertex = new HepMC::GenVertex();
123 0 : m_particle->parent_event()->add_vertex(end_vertex);
124 0 : }
125 :
126 : //Loop over all daughters to check that the end points to the right place
127 0 : vector<PhotosParticle*>::iterator daughter_itr;
128 0 : for(daughter_itr = daughters.begin(); daughter_itr != daughters.end();
129 0 : daughter_itr++){
130 :
131 : HepMC::GenParticle * daug;
132 0 : daug = dynamic_cast<PhotosHepMCParticle*>(*daughter_itr)->getHepMC();
133 :
134 :
135 0 : if(daug->production_vertex()!=orig_end_vertex)
136 0 : Log::Fatal("PhotosHepMCParticle::setDaughters(): Daughter production_vertices point to difference places. Can not override. Please delete vertices first.",4);
137 : else
138 0 : end_vertex->add_particle_out(daug);
139 : }
140 0 : end_vertex->add_particle_in(m_particle);
141 0 : }
142 :
143 0 : }
144 :
145 : std::vector<PhotosParticle*> PhotosHepMCParticle::getMothers(){
146 :
147 0 : if(m_mothers.size()==0&&m_particle->production_vertex()){
148 :
149 0 : HepMC::GenVertex::particles_in_const_iterator pcle_itr;
150 0 : pcle_itr=m_particle->production_vertex()->particles_in_const_begin();
151 :
152 0 : HepMC::GenVertex::particles_in_const_iterator pcle_itr_end;
153 0 : pcle_itr_end=m_particle->production_vertex()->particles_in_const_end();
154 :
155 0 : for(;pcle_itr != pcle_itr_end; pcle_itr++){
156 0 : m_mothers.push_back(new PhotosHepMCParticle(*pcle_itr));
157 : }
158 0 : }
159 0 : return m_mothers;
160 : }
161 :
162 : std::vector<PhotosParticle*> PhotosHepMCParticle::getDaughters(){
163 :
164 0 : if(m_daughters.size()==0&&m_particle->end_vertex()){
165 0 : HepMC::GenVertex::particles_out_const_iterator pcle_itr;
166 0 : pcle_itr=m_particle->end_vertex()->particles_out_const_begin();
167 :
168 0 : HepMC::GenVertex::particles_out_const_iterator pcle_itr_end;
169 0 : pcle_itr_end=m_particle->end_vertex()->particles_out_const_end();
170 :
171 0 : for(;pcle_itr != pcle_itr_end; pcle_itr++){
172 :
173 : // ommit particles if their status code is ignored by Photos
174 0 : if( Photos::isStatusCodeIgnored( (*pcle_itr)->status() ) ) continue;
175 :
176 0 : m_daughters.push_back(new PhotosHepMCParticle(*pcle_itr));
177 0 : }
178 0 : }
179 0 : return m_daughters;
180 :
181 : }
182 :
183 : std::vector<PhotosParticle*> PhotosHepMCParticle::getAllDecayProducts(){
184 :
185 0 : m_decay_products.clear();
186 :
187 0 : if(!hasDaughters()) // if this particle has no daughters
188 0 : return m_decay_products;
189 :
190 0 : std::vector<PhotosParticle*> daughters = getDaughters();
191 :
192 : // copy daughters to list of all decay products
193 0 : m_decay_products.insert(m_decay_products.end(),daughters.begin(),daughters.end());
194 :
195 : // Now, get all daughters recursively, without duplicates.
196 : // That is, for each daughter:
197 : // 1) get list of her daughters
198 : // 2) for each particle on this list:
199 : // a) check if it is already on the list
200 : // b) if it's not, add her to the end of the list
201 0 : for(unsigned int i=0;i<m_decay_products.size();i++)
202 : {
203 0 : std::vector<PhotosParticle*> daughters2 = m_decay_products[i]->getDaughters();
204 :
205 0 : if(!m_decay_products[i]->hasDaughters()) continue;
206 0 : for(unsigned int j=0;j<daughters2.size();j++)
207 : {
208 : bool add=true;
209 0 : for(unsigned int k=0;k<m_decay_products.size();k++)
210 0 : if( daughters2[j]->getBarcode() == m_decay_products[k]->getBarcode() )
211 : {
212 : add=false;
213 0 : break;
214 : }
215 :
216 0 : if(add) m_decay_products.push_back(daughters2[j]);
217 : }
218 0 : }
219 :
220 0 : return m_decay_products;
221 0 : }
222 :
223 : bool PhotosHepMCParticle::checkMomentumConservation(){
224 :
225 0 : if(!m_particle->end_vertex()) return true;
226 :
227 : // HepMC version of check_momentum_conservation
228 : // Ommitting history entries (status == 3)
229 :
230 : double sumpx = 0, sumpy = 0, sumpz = 0, sume = 0;
231 0 : for( HepMC::GenVertex::particles_in_const_iterator part1 = m_particle->end_vertex()->particles_in_const_begin();
232 0 : part1 != m_particle->end_vertex()->particles_in_const_end(); part1++ ){
233 :
234 0 : if( Photos::isStatusCodeIgnored((*part1)->status()) ) continue;
235 :
236 0 : sumpx += (*part1)->momentum().px();
237 0 : sumpy += (*part1)->momentum().py();
238 0 : sumpz += (*part1)->momentum().pz();
239 0 : sume += (*part1)->momentum().e();
240 0 : }
241 :
242 0 : for( HepMC::GenVertex::particles_out_const_iterator part2 = m_particle->end_vertex()->particles_out_const_begin();
243 0 : part2 != m_particle->end_vertex()->particles_out_const_end(); part2++ ){
244 :
245 0 : if( Photos::isStatusCodeIgnored((*part2)->status()) ) continue;
246 :
247 0 : sumpx -= (*part2)->momentum().px();
248 0 : sumpy -= (*part2)->momentum().py();
249 0 : sumpz -= (*part2)->momentum().pz();
250 0 : sume -= (*part2)->momentum().e();
251 0 : }
252 :
253 0 : if( sqrt( sumpx*sumpx + sumpy*sumpy + sumpz*sumpz + sume*sume) > Photos::momentum_conservation_threshold ) {
254 0 : Log::Warning()<<"Momentum not conserved in the vertex:"<<endl;
255 0 : Log::RedirectOutput(Log::Warning(false));
256 0 : m_particle->end_vertex()->print();
257 0 : Log::RevertOutput();
258 0 : return false;
259 : }
260 :
261 0 : return true;
262 0 : }
263 :
264 : void PhotosHepMCParticle::setPdgID(int pdg_id){
265 0 : m_particle->set_pdg_id(pdg_id);
266 0 : }
267 :
268 : void PhotosHepMCParticle::setMass(double mass){
269 0 : m_particle->set_generated_mass(mass);
270 0 : }
271 :
272 : void PhotosHepMCParticle::setStatus(int status){
273 0 : m_particle->set_status(status);
274 0 : }
275 :
276 : int PhotosHepMCParticle::getPdgID(){
277 0 : return m_particle->pdg_id();
278 : }
279 :
280 : int PhotosHepMCParticle::getStatus(){
281 0 : return m_particle->status();
282 : }
283 :
284 : int PhotosHepMCParticle::getBarcode(){
285 0 : return m_particle->barcode();
286 : }
287 :
288 :
289 : PhotosHepMCParticle * PhotosHepMCParticle::createNewParticle(
290 : int pdg_id, int status, double mass,
291 : double px, double py, double pz, double e){
292 :
293 0 : PhotosHepMCParticle * new_particle = new PhotosHepMCParticle();
294 0 : new_particle->getHepMC()->set_pdg_id(pdg_id);
295 0 : new_particle->getHepMC()->set_status(status);
296 0 : new_particle->getHepMC()->set_generated_mass(mass);
297 :
298 0 : HepMC::FourVector momentum(px,py,pz,e);
299 0 : new_particle->getHepMC()->set_momentum(momentum);
300 :
301 0 : m_created_particles.push_back(new_particle);
302 0 : return new_particle;
303 0 : }
304 :
305 : void PhotosHepMCParticle::createHistoryEntry(){
306 :
307 0 : if(!m_particle->production_vertex())
308 : {
309 0 : Log::Warning()<<"PhotosHepMCParticle::createHistoryEntry(): particle without production vertex."<<endl;
310 0 : return;
311 : }
312 :
313 0 : HepMC::GenParticle *part = new HepMC::GenParticle(*m_particle);
314 0 : part->set_status(Photos::historyEntriesStatus);
315 0 : m_particle->production_vertex()->add_particle_out(part);
316 0 : }
317 :
318 : void PhotosHepMCParticle::createSelfDecayVertex(PhotosParticle *out)
319 : {
320 0 : if(m_particle->end_vertex())
321 : {
322 0 : Log::Error()<<"PhotosHepMCParticle::createSelfDecayVertex: particle already has end vertex!"<<endl;
323 0 : return;
324 : }
325 :
326 0 : if(getHepMC()->parent_event()==NULL)
327 : {
328 0 : Log::Error()<<"PhotosHepMCParticle::createSelfDecayVertex: particle not in the HepMC event!"<<endl;
329 0 : return;
330 : }
331 :
332 : // Add new vertex and new particle to HepMC
333 0 : HepMC::GenParticle *outgoing = new HepMC::GenParticle( *(dynamic_cast<PhotosHepMCParticle*>(out)->m_particle) );
334 0 : HepMC::GenVertex *v = new HepMC::GenVertex();
335 :
336 : // Copy vertex position from parent vertex
337 0 : v->set_position( m_particle->production_vertex()->position() );
338 :
339 0 : v->add_particle_in (m_particle);
340 0 : v->add_particle_out(outgoing);
341 :
342 0 : getHepMC()->parent_event()->add_vertex(v);
343 :
344 : // If this particle was stable, set its status to 2
345 0 : if(getStatus()==1) setStatus(2);
346 0 : }
347 :
348 : void PhotosHepMCParticle::print(){
349 0 : m_particle->print();
350 0 : }
351 :
352 :
353 : /******** Getter and Setter methods: ***********************/
354 :
355 : inline double PhotosHepMCParticle::getPx(){
356 0 : return m_particle->momentum().px();
357 : }
358 :
359 : inline double PhotosHepMCParticle::getPy(){
360 0 : return m_particle->momentum().py();
361 : }
362 :
363 : double PhotosHepMCParticle::getPz(){
364 0 : return m_particle->momentum().pz();
365 : }
366 :
367 : double PhotosHepMCParticle::getE(){
368 0 : return m_particle->momentum().e();
369 : }
370 :
371 : void PhotosHepMCParticle::setPx(double px){
372 : //make new momentum as something is wrong with
373 : //the HepMC momentum setters
374 :
375 0 : HepMC::FourVector momentum(m_particle->momentum());
376 0 : momentum.setPx(px);
377 0 : m_particle->set_momentum(momentum);
378 0 : }
379 :
380 : void PhotosHepMCParticle::setPy(double py){
381 0 : HepMC::FourVector momentum(m_particle->momentum());
382 0 : momentum.setPy(py);
383 0 : m_particle->set_momentum(momentum);
384 0 : }
385 :
386 :
387 : void PhotosHepMCParticle::setPz(double pz){
388 0 : HepMC::FourVector momentum(m_particle->momentum());
389 0 : momentum.setPz(pz);
390 0 : m_particle->set_momentum(momentum);
391 0 : }
392 :
393 : void PhotosHepMCParticle::setE(double e){
394 0 : HepMC::FourVector momentum(m_particle->momentum());
395 0 : momentum.setE(e);
396 0 : m_particle->set_momentum(momentum);
397 0 : }
398 :
399 : double PhotosHepMCParticle::getMass()
400 : {
401 0 : return m_particle->generated_mass();
402 : }
403 :
404 : } // namespace Photospp
|
