Line data Source code
1 : #include "TauolaHepMCParticle.h"
2 : #include "TauolaLog.h"
3 :
4 : namespace Tauolapp
5 : {
6 :
7 0 : TauolaHepMCParticle::TauolaHepMCParticle(){
8 0 : m_particle = new HepMC::GenParticle();
9 0 : }
10 :
11 0 : TauolaHepMCParticle::~TauolaHepMCParticle(){
12 :
13 : //delete the mother and daughter pointers
14 0 : while(m_mothers.size()!=0){
15 0 : TauolaParticle * temp = m_mothers.back();
16 0 : m_mothers.pop_back();
17 0 : delete temp;
18 : }
19 0 : while(m_daughters.size()!=0){
20 0 : TauolaParticle * temp = m_daughters.back();
21 0 : m_daughters.pop_back();
22 0 : delete temp;
23 : }
24 :
25 0 : while(m_created_particles.size()!=0){
26 0 : TauolaHepMCParticle * temp = (TauolaHepMCParticle*) m_created_particles.back();
27 0 : m_created_particles.pop_back();
28 0 : if(temp->getHepMC()->barcode()==0) delete temp->getHepMC();
29 0 : delete temp;
30 : }
31 :
32 0 : }
33 :
34 : // NOTE: Not executed by release examples
35 0 : TauolaHepMCParticle::TauolaHepMCParticle(int pdg_id, int status, double mass){
36 0 : m_particle = new HepMC::GenParticle();
37 0 : m_particle->set_pdg_id(pdg_id);
38 0 : m_particle->set_status(status);
39 0 : m_particle->set_generated_mass(mass);
40 0 : }
41 :
42 0 : TauolaHepMCParticle::TauolaHepMCParticle(HepMC::GenParticle * particle){
43 0 : m_particle = particle;
44 0 : }
45 :
46 : HepMC::GenParticle * TauolaHepMCParticle::getHepMC(){
47 0 : return m_particle;
48 : }
49 :
50 : void TauolaHepMCParticle::undecay(){
51 0 : std::vector<TauolaParticle*> daughters = getDaughters();
52 0 : std::vector<TauolaParticle*>::iterator dIter = daughters.begin();
53 :
54 0 : for(; dIter != daughters.end(); dIter++)
55 0 : (*dIter)->undecay();
56 :
57 0 : if(m_particle->end_vertex())
58 : {
59 0 : while(m_particle->end_vertex()->particles_out_size())
60 : {
61 0 : HepMC::GenParticle *p = m_particle->end_vertex()->remove_particle(*(m_particle->end_vertex()->particles_out_const_begin()));
62 0 : delete p;
63 : }
64 0 : delete m_particle->end_vertex();
65 : }
66 :
67 0 : m_daughters.clear();
68 0 : m_particle->set_status(TauolaParticle::STABLE);
69 :
70 0 : for(unsigned int i=0;i<daughters.size();i++)
71 0 : delete daughters[i];
72 0 : }
73 :
74 : void TauolaHepMCParticle::setMothers(vector<TauolaParticle*> mothers){
75 :
76 : /******** Deal with mothers ***********/
77 :
78 : //If there are mothers
79 0 : if(mothers.size()>0){
80 :
81 : HepMC::GenParticle * part;
82 0 : part=dynamic_cast<TauolaHepMCParticle*>(mothers.at(0))->getHepMC();
83 :
84 : //Use end vertex of first mother as production vertex for particle
85 0 : HepMC::GenVertex * production_vertex = part->end_vertex();
86 : HepMC::GenVertex * orig_production_vertex = production_vertex;
87 :
88 : //If production_vertex does not exist - create it
89 : //If it's tau decay - set the time and position including the tau lifetime correction
90 : //otherwise - copy the time and position of decaying particle
91 0 : if(!production_vertex){
92 0 : production_vertex = new HepMC::GenVertex();
93 0 : HepMC::FourVector point = part->production_vertex()->position();
94 0 : production_vertex->set_position(point);
95 0 : part->parent_event()->add_vertex(production_vertex);
96 0 : }
97 :
98 : //Loop over all mothers to check that the end points to the right place
99 0 : vector<TauolaParticle*>::iterator mother_itr;
100 0 : for(mother_itr = mothers.begin(); mother_itr != mothers.end();
101 0 : mother_itr++){
102 :
103 : HepMC::GenParticle * moth;
104 0 : moth = dynamic_cast<TauolaHepMCParticle*>(*mother_itr)->getHepMC();
105 :
106 0 : if(moth->end_vertex()!=orig_production_vertex)
107 0 : Log::Fatal("Mother production_vertices point to difference places. Can not override. Please delete vertices first.",1);
108 : else
109 0 : production_vertex->add_particle_in(moth);
110 :
111 : //update status info
112 0 : if(moth->status()==TauolaParticle::STABLE)
113 0 : moth->set_status(TauolaParticle::DECAYED);
114 : }
115 0 : production_vertex->add_particle_out(m_particle);
116 0 : }
117 0 : }
118 :
119 : void TauolaHepMCParticle::setDaughters(vector<TauolaParticle*> daughters){
120 :
121 0 : if(!m_particle->parent_event())
122 0 : Log::Fatal("New particle needs the event set before it's daughters can be added",2);
123 :
124 : //If there are daughters
125 0 : if(daughters.size()>0){
126 : // NOTE: Not executed by release examples
127 : // because daughters.size() is always 0
128 :
129 : //Use production vertex of first daughter as end vertex for particle
130 : HepMC::GenParticle * first_daughter;
131 0 : first_daughter = (dynamic_cast<TauolaHepMCParticle*>(daughters.at(0)))->getHepMC();
132 :
133 : HepMC::GenVertex * end_vertex;
134 0 : end_vertex=first_daughter->production_vertex();
135 : HepMC::GenVertex * orig_end_vertex = end_vertex;
136 :
137 0 : if(!end_vertex){ //if it does not exist create it
138 0 : end_vertex = new HepMC::GenVertex();
139 0 : m_particle->parent_event()->add_vertex(end_vertex);
140 0 : }
141 :
142 : //Loop over all daughters to check that the end points to the right place
143 0 : vector<TauolaParticle*>::iterator daughter_itr;
144 0 : for(daughter_itr = daughters.begin(); daughter_itr != daughters.end();
145 0 : daughter_itr++){
146 :
147 : HepMC::GenParticle * daug;
148 0 : daug = dynamic_cast<TauolaHepMCParticle*>(*daughter_itr)->getHepMC();
149 :
150 :
151 0 : if(daug->production_vertex()!=orig_end_vertex)
152 0 : Log::Fatal("Daughter production_vertices point to difference places. Can not override. Please delete vertices first.",3);
153 : else
154 0 : end_vertex->add_particle_out(daug);
155 : }
156 0 : end_vertex->add_particle_in(m_particle);
157 0 : }
158 0 : }
159 :
160 : std::vector<TauolaParticle*> TauolaHepMCParticle::getMothers(){
161 :
162 0 : if(m_mothers.size()==0&&m_particle->production_vertex()){
163 0 : HepMC::GenVertex::particles_in_const_iterator pcle_itr;
164 0 : pcle_itr=m_particle->production_vertex()->particles_in_const_begin();
165 :
166 0 : HepMC::GenVertex::particles_in_const_iterator pcle_itr_end;
167 0 : pcle_itr_end=m_particle->production_vertex()->particles_in_const_end();
168 :
169 0 : for(;pcle_itr != pcle_itr_end; pcle_itr++){
170 0 : m_mothers.push_back(new TauolaHepMCParticle(*pcle_itr));
171 : }
172 0 : }
173 0 : return m_mothers;
174 : }
175 :
176 : std::vector<TauolaParticle*> TauolaHepMCParticle::getDaughters(){
177 :
178 0 : if(m_daughters.size()==0&&m_particle->end_vertex()){
179 0 : HepMC::GenVertex::particles_out_const_iterator pcle_itr;
180 0 : pcle_itr=m_particle->end_vertex()->particles_out_const_begin();
181 :
182 0 : HepMC::GenVertex::particles_out_const_iterator pcle_itr_end;
183 0 : pcle_itr_end=m_particle->end_vertex()->particles_out_const_end();
184 :
185 0 : for(;pcle_itr != pcle_itr_end; pcle_itr++){
186 0 : m_daughters.push_back(new TauolaHepMCParticle(*pcle_itr));
187 : }
188 0 : }
189 0 : return m_daughters;
190 : }
191 :
192 : void TauolaHepMCParticle::checkMomentumConservation(){
193 :
194 0 : if(!m_particle->end_vertex()) return;
195 :
196 : // HepMC version of check_momentum_conservation
197 : // with added energy check
198 :
199 : double sumpx = 0, sumpy = 0, sumpz = 0, sume = 0;
200 0 : for( HepMC::GenVertex::particles_in_const_iterator part1 = m_particle->end_vertex()->particles_in_const_begin();
201 0 : part1 != m_particle->end_vertex()->particles_in_const_end(); part1++ ){
202 :
203 0 : sumpx += (*part1)->momentum().px();
204 0 : sumpy += (*part1)->momentum().py();
205 0 : sumpz += (*part1)->momentum().pz();
206 0 : sume += (*part1)->momentum().e();
207 : }
208 :
209 0 : for( HepMC::GenVertex::particles_out_const_iterator part2 = m_particle->end_vertex()->particles_out_const_begin();
210 0 : part2 != m_particle->end_vertex()->particles_out_const_end(); part2++ ){
211 :
212 0 : sumpx -= (*part2)->momentum().px();
213 0 : sumpy -= (*part2)->momentum().py();
214 0 : sumpz -= (*part2)->momentum().pz();
215 0 : sume -= (*part2)->momentum().e();
216 : }
217 :
218 0 : if( sqrt( sumpx*sumpx + sumpy*sumpy + sumpz*sumpz + sume*sume) > Tauola::momentum_conservation_threshold ) {
219 0 : Log::Warning()<<"Momentum not conserved in the vertex:"<<endl;
220 0 : Log::RedirectOutput(Log::Warning(false));
221 0 : m_particle->end_vertex()->print();
222 0 : Log::RevertOutput();
223 0 : return;
224 : }
225 :
226 0 : return;
227 0 : }
228 :
229 : // NOTE: Not executed by release examples
230 : void TauolaHepMCParticle::setPdgID(int pdg_id){
231 0 : m_particle->set_pdg_id(pdg_id);
232 0 : }
233 :
234 : void TauolaHepMCParticle::setMass(double mass){
235 0 : m_particle->set_generated_mass(mass);
236 0 : }
237 :
238 : // NOTE: Not executed by release examples
239 : void TauolaHepMCParticle::setStatus(int status){
240 0 : m_particle->set_status(status);
241 0 : }
242 :
243 : int TauolaHepMCParticle::getPdgID(){
244 0 : return m_particle->pdg_id();
245 : }
246 :
247 : int TauolaHepMCParticle::getStatus(){
248 0 : return m_particle->status();
249 : }
250 :
251 : int TauolaHepMCParticle::getBarcode(){
252 0 : return m_particle->barcode();
253 : }
254 :
255 : // Set (X,T) Position of tau decay trees
256 : void TauolaHepMCParticle::decayEndgame(){
257 :
258 0 : double lifetime = Tauola::tau_lifetime * (-log( Tauola::randomDouble() ));
259 0 : HepMC::FourVector tau_momentum = m_particle->momentum();
260 :
261 0 : double mass = sqrt(abs( tau_momentum.e()*tau_momentum.e()
262 0 : - tau_momentum.px()*tau_momentum.px()
263 0 : - tau_momentum.py()*tau_momentum.py()
264 0 : - tau_momentum.pz()*tau_momentum.pz()
265 : ) );
266 :
267 : // Get previous position
268 0 : HepMC::FourVector previous_position = m_particle->production_vertex()->position();
269 :
270 : // Calculate new position
271 0 : HepMC::FourVector new_position(previous_position.x()+tau_momentum.px()/mass*lifetime,
272 0 : previous_position.y()+tau_momentum.py()/mass*lifetime,
273 0 : previous_position.z()+tau_momentum.pz()/mass*lifetime,
274 0 : previous_position.t()+tau_momentum.e() /mass*lifetime);
275 :
276 : // Set new position
277 0 : m_particle->end_vertex()->set_position(new_position);
278 0 : recursiveSetPosition(m_particle,new_position);
279 0 : }
280 :
281 : void TauolaHepMCParticle::recursiveSetPosition(HepMC::GenParticle *p, HepMC::FourVector pos){
282 :
283 0 : if(!p->end_vertex()) return;
284 :
285 : // Iterate over all outgoing particles
286 0 : for(HepMC::GenVertex::particles_out_const_iterator pp = p->end_vertex()->particles_out_const_begin();
287 0 : pp != p->end_vertex()->particles_out_const_end();
288 0 : ++pp){
289 0 : if( !(*pp)->end_vertex() ) continue;
290 :
291 : // Set position
292 0 : (*pp)->end_vertex()->set_position(pos);
293 0 : recursiveSetPosition(*pp,pos);
294 0 : }
295 0 : }
296 :
297 : TauolaHepMCParticle * TauolaHepMCParticle::createNewParticle(
298 : int pdg_id, int status, double mass,
299 : double px, double py, double pz, double e){
300 :
301 0 : TauolaHepMCParticle * new_particle = new TauolaHepMCParticle();
302 0 : new_particle->getHepMC()->set_pdg_id(pdg_id);
303 0 : new_particle->getHepMC()->set_status(status);
304 0 : new_particle->getHepMC()->set_generated_mass(mass);
305 :
306 0 : HepMC::FourVector momentum(px,py,pz,e);
307 0 : new_particle->getHepMC()->set_momentum(momentum);
308 :
309 0 : m_created_particles.push_back(new_particle);
310 :
311 0 : return new_particle;
312 0 : }
313 :
314 : void TauolaHepMCParticle::print(){
315 0 : m_particle->print();
316 0 : }
317 :
318 :
319 : /******** Getter and Setter methods: ***********************/
320 :
321 : inline double TauolaHepMCParticle::getPx(){
322 0 : return m_particle->momentum().px();
323 : }
324 :
325 : inline double TauolaHepMCParticle::getPy(){
326 0 : return m_particle->momentum().py();
327 : }
328 :
329 : double TauolaHepMCParticle::getPz(){
330 0 : return m_particle->momentum().pz();
331 : }
332 :
333 : double TauolaHepMCParticle::getE(){
334 0 : return m_particle->momentum().e();
335 : }
336 :
337 : void TauolaHepMCParticle::setPx(double px){
338 : //make new momentum as something is wrong with
339 : //the HepMC momentum setters
340 :
341 0 : HepMC::FourVector momentum(m_particle->momentum());
342 0 : momentum.setPx(px);
343 0 : m_particle->set_momentum(momentum);
344 0 : }
345 :
346 : void TauolaHepMCParticle::setPy(double py){
347 0 : HepMC::FourVector momentum(m_particle->momentum());
348 0 : momentum.setPy(py);
349 0 : m_particle->set_momentum(momentum);
350 0 : }
351 :
352 :
353 : void TauolaHepMCParticle::setPz(double pz){
354 0 : HepMC::FourVector momentum(m_particle->momentum());
355 0 : momentum.setPz(pz);
356 0 : m_particle->set_momentum(momentum);
357 0 : }
358 :
359 : void TauolaHepMCParticle::setE(double e){
360 0 : HepMC::FourVector momentum(m_particle->momentum());
361 0 : momentum.setE(e);
362 0 : m_particle->set_momentum(momentum);
363 0 : }
364 :
365 : } // namespace Tauolapp
|
