Line data Source code
1 : #ifdef EVTGEN_TAUOLA
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: EvtTauolaEngine
12 : //
13 : // Description: Interface to the TAUOLA external generator, which
14 : // decays tau particles
15 : //
16 : // Modification history:
17 : //
18 : // John Back May 2011 Module created
19 : //
20 : //------------------------------------------------------------------------
21 :
22 : #include "EvtGenExternal/EvtTauolaEngine.hh"
23 :
24 : #include "EvtGenBase/EvtPDL.hh"
25 : #include "EvtGenBase/EvtVector4R.hh"
26 : #include "EvtGenBase/EvtDecayTable.hh"
27 : #include "EvtGenBase/EvtRandom.hh"
28 : #include "EvtGenBase/EvtReport.hh"
29 :
30 : #include "Tauola/Tauola.h"
31 : #include "Tauola/TauolaHepMCEvent.h"
32 : #include "Tauola/TauolaHepMCParticle.h"
33 : #include "Tauola/TauolaParticle.h"
34 :
35 : #include "HepMC/GenVertex.h"
36 : #include "HepMC/SimpleVector.h"
37 : #include "HepMC/Units.h"
38 :
39 : #include <iostream>
40 : #include <sstream>
41 : #include <string>
42 : #include <cmath>
43 :
44 : using std::endl;
45 :
46 0 : EvtTauolaEngine::EvtTauolaEngine(bool useEvtGenRandom) {
47 :
48 : // PDG standard code integer ID for tau particle
49 0 : _tauPDG = 15;
50 : // Number of possible decay modes in Tauola
51 0 : _nTauolaModes = 22;
52 :
53 0 : report(Severity::Info,"EvtGen")<<"Setting up TAUOLA."<<endl;
54 :
55 : // These three lines are not really necessary since they are the default.
56 : // But they are here so that we know what the initial conditions are.
57 0 : Tauolapp::Tauola::setDecayingParticle(_tauPDG); // tau PDG code
58 0 : Tauolapp::Tauola::setSameParticleDecayMode(Tauolapp::Tauola::All); // all modes allowed
59 0 : Tauolapp::Tauola::setOppositeParticleDecayMode(Tauolapp::Tauola::All); // all modes allowed
60 :
61 : // Initial the Tauola external generator
62 0 : if (useEvtGenRandom == true) {
63 :
64 0 : report(Severity::Info,"EvtGen")<<"Using EvtGen random number engine also for Tauola++"<<endl;
65 :
66 0 : Tauolapp::Tauola::setRandomGenerator(EvtRandom::Flat);
67 :
68 : }
69 :
70 0 : Tauolapp::Tauola::initialize();
71 :
72 : // Set-up possible decay modes when we have read the (user) decay file
73 0 : _initialised = false;
74 :
75 0 : }
76 :
77 0 : EvtTauolaEngine::~EvtTauolaEngine() {
78 :
79 0 : }
80 :
81 : void EvtTauolaEngine::initialise() {
82 :
83 : // Set up all possible tau decay modes.
84 : // This should be done just before the first doDecay() call,
85 : // since we want to make sure that any decay.dec files are processed
86 : // first to get lists of particle modes and their alias definitions
87 : // (for creating EvtParticles with the right history information).
88 :
89 0 : if (_initialised == false) {
90 :
91 0 : this->setUpPossibleTauModes();
92 :
93 0 : _initialised = true;
94 :
95 0 : }
96 :
97 0 : }
98 :
99 : void EvtTauolaEngine::setUpPossibleTauModes() {
100 :
101 : // Get the decay table list defined by the decay.dec files.
102 : // Only look for the first tau particle decay mode definitions with the Tauola name,
103 : // since that generator only allows the same BFs for both tau+ and tau- decays.
104 : // We can not choose a specific tau decay event-by-event, since this is
105 : // only possible before we call Tauola::initialize().
106 : // Otherwise, we could have selected a random mode ourselves for tau- and tau+
107 : // separately (via selecting a random number and comparing it to be less than
108 : // the cumulative BF) for each event.
109 :
110 0 : int nPDL = EvtPDL::entries();
111 : int iPDL(0);
112 :
113 : bool gotAnyTauolaModes(false);
114 :
115 0 : for (iPDL = 0; iPDL < nPDL; iPDL++) {
116 :
117 0 : EvtId particleId = EvtPDL::getEntry(iPDL);
118 0 : int PDGId = EvtPDL::getStdHep(particleId);
119 :
120 0 : if (abs(PDGId) == _tauPDG && gotAnyTauolaModes == false) {
121 :
122 0 : int aliasInt = particleId.getAlias();
123 :
124 : // Get the list of decay modes for this tau particle (alias)
125 0 : int nModes = EvtDecayTable::getInstance()->getNModes(aliasInt);
126 : int iMode(0), iTauMode(0);
127 :
128 : // Vector to store tau mode branching fractions.
129 : // The size of this vector equals the total number of possible
130 : // Tauola decay modes. Initialise all BFs to zero.
131 0 : std::vector<double> tauolaModeBFs(_nTauolaModes);
132 :
133 0 : for (iTauMode = 0; iTauMode < _nTauolaModes; iTauMode++) {
134 0 : tauolaModeBFs[iTauMode] = 0.0;
135 : }
136 :
137 : double totalTauModeBF(0.0);
138 :
139 : int nNonTauolaModes(0);
140 :
141 : // Loop through each decay mode
142 0 : for (iMode = 0; iMode < nModes; iMode++) {
143 :
144 0 : EvtDecayBase* decayModel = EvtDecayTable::getInstance()->findDecayModel(aliasInt, iMode);
145 0 : if (decayModel != 0) {
146 :
147 : // Check that the decay model name matches TAUOLA
148 0 : std::string modelName = decayModel->getName();
149 0 : if (modelName == "TAUOLA") {
150 :
151 0 : if (gotAnyTauolaModes == false) {gotAnyTauolaModes = true;}
152 :
153 : // Extract the decay mode integer type and branching fraction
154 0 : double BF = decayModel->getBranchingFraction();
155 0 : int modeArrayInt = this->getModeInt(decayModel) - 1;
156 :
157 0 : if (modeArrayInt >= 0 && modeArrayInt < _nTauolaModes) {
158 0 : tauolaModeBFs[modeArrayInt] = BF;
159 0 : totalTauModeBF += BF;
160 0 : }
161 :
162 0 : } else {
163 :
164 0 : nNonTauolaModes++;
165 :
166 : }
167 :
168 0 : } // Decay mode exists
169 :
170 : } // Loop over decay models
171 :
172 0 : if (gotAnyTauolaModes == true && nNonTauolaModes > 0) {
173 :
174 0 : report(Severity::Error, "EvtGen") << "Please remove all non-TAUOLA decay modes for particle "
175 0 : <<EvtPDL::name(particleId)<<endl;
176 0 : ::abort();
177 :
178 : }
179 :
180 : // Normalise all (non-zero) tau mode BFs to sum up to 1.0, and
181 : // let Tauola know about these normalised branching fractions
182 0 : if (totalTauModeBF > 0.0) {
183 :
184 0 : report(Severity::Info,"EvtGen")<<"Setting TAUOLA BF modes using the definitions for the particle "
185 0 : <<EvtPDL::name(particleId)<<endl;
186 :
187 0 : for (iTauMode = 0; iTauMode < _nTauolaModes; iTauMode++) {
188 :
189 0 : tauolaModeBFs[iTauMode] /= totalTauModeBF;
190 0 : double modeBF = tauolaModeBFs[iTauMode];
191 0 : report(Severity::Info,"EvtGen")<<"Setting TAUOLA BF for mode "<<iTauMode+1<<" = "<<modeBF<<endl;
192 0 : Tauolapp::Tauola::setTauBr(iTauMode+1, modeBF);
193 :
194 : }
195 :
196 0 : report(Severity::Info,"EvtGen")<<"Any other TAUOLA BF modes for other tau particle decay mode definitions will be ignored!"
197 0 : <<endl;
198 :
199 : }
200 :
201 0 : } // Got tau particle and have yet to get a TAUOLA mode
202 :
203 0 : } // Loop through PDL entries
204 :
205 0 : }
206 :
207 : int EvtTauolaEngine::getModeInt(EvtDecayBase* decayModel) {
208 :
209 : int modeInt(0);
210 :
211 0 : if (decayModel != 0) {
212 :
213 0 : int nVars = decayModel->getNArg();
214 :
215 0 : if (nVars > 0) {
216 0 : modeInt = static_cast<int>(decayModel->getArg(0));
217 0 : }
218 :
219 0 : }
220 :
221 0 : return modeInt;
222 :
223 : }
224 :
225 : bool EvtTauolaEngine::doDecay(EvtParticle* tauParticle) {
226 :
227 0 : if (_initialised == false) {this->initialise();}
228 :
229 0 : if (tauParticle == 0) {return false;}
230 :
231 : // Check that we have a tau particle.
232 0 : EvtId partId = tauParticle->getId();
233 0 : if (abs(EvtPDL::getStdHep(partId)) != _tauPDG) {return false;}
234 :
235 0 : int nTauDaug = tauParticle->getNDaug();
236 :
237 : // If the number of tau daughters is not zero, then we have already decayed
238 : // it using Tauola/another decay algorithm.
239 0 : if (nTauDaug > 0) {return true;}
240 :
241 0 : this->decayTauEvent(tauParticle);
242 :
243 0 : return true;
244 :
245 0 : }
246 :
247 : void EvtTauolaEngine::decayTauEvent(EvtParticle* tauParticle) {
248 :
249 : // Either we have a tau particle within a decay chain, or a single particle.
250 : // Create a dummy HepMC event & vertex for the parent particle, containing the tau as
251 : // one of the outgoing particles. If we have a decay chain, the parent will be the
252 : // incoming particle, while the daughters, including the tau, are outgoing particles.
253 : // For the single particle case, the incoming particle is null, while the single tau
254 : // is the only outgoing particle.
255 : // We can then pass this event to Tauola which should then decay the tau particle.
256 : // We also consider all other tau particles from the parent decay in the logic below.
257 :
258 : // Create the dummy event.
259 0 : HepMC::GenEvent* theEvent = new HepMC::GenEvent(HepMC::Units::GEV, HepMC::Units::MM);
260 :
261 : // Create the decay "vertex".
262 0 : HepMC::GenVertex* theVertex = new HepMC::GenVertex();
263 0 : theEvent->add_vertex(theVertex);
264 :
265 : // Get the parent of this tau particle
266 0 : EvtParticle* theParent = tauParticle->getParent();
267 :
268 : // Assign the parent particle as the incoming particle to the vertex.
269 0 : if (theParent != 0) {
270 0 : HepMC::GenParticle* hepMCParent = this->createGenParticle(theParent);
271 0 : theVertex->add_particle_in(hepMCParent);
272 0 : } else {
273 : // The tau particle has no parent. Set "itself" as the incoming particle for the first vertex.
274 : // This is needed, otherwise Tauola warns of momentum non-conservation for this (1st) vertex.
275 0 : HepMC::GenParticle* tauGenInit = this->createGenParticle(tauParticle);
276 0 : theVertex->add_particle_in(tauGenInit);
277 : }
278 :
279 : // Find all daughter particles and assign them as outgoing particles to the vertex.
280 : // This will include the tau particle we are currently processing.
281 : // If the parent decay has more than one tau particle, we need to include them as well.
282 : // This is important since Tauola needs the correct physics correlations: we do not
283 : // want Tauola to decay each particle separately if they are from tau pair combinations.
284 : // Tauola will process the event, and we will create EvtParticles from all tau decay
285 : // products, i.e. the tau particle we currently have and any other tau particles.
286 : // EvtGen will try to decay the other tau particle(s) by calling EvtTauola and therefore
287 : // this function. However, we check to see if the tau candidate has any daughters already.
288 : // If it does, then we have already set the tau decay products from Tauola.
289 :
290 : // Map to store (HepMC,EvtParticle) pairs for each tau candidate from the parent
291 : // decay. This is needed to find out what EvtParticle corresponds to a given tau HepMC
292 : // candidate: we do not want to recreate existing EvtParticle pointers.
293 0 : std::map<HepMC::GenParticle*, EvtParticle*> tauMap;
294 :
295 0 : if (theParent != 0) {
296 :
297 : // Find all tau particles in the decay tree and store them in the map
298 0 : int nDaug(theParent->getNDaug());
299 : int iDaug(0);
300 :
301 0 : for (iDaug = 0; iDaug < nDaug; iDaug++) {
302 :
303 0 : EvtParticle* theDaughter = theParent->getDaug(iDaug);
304 :
305 0 : if (theDaughter != 0) {
306 :
307 0 : HepMC::GenParticle* hepMCDaughter = this->createGenParticle(theDaughter);
308 0 : theVertex->add_particle_out(hepMCDaughter);
309 :
310 0 : EvtId theId = theDaughter->getId();
311 0 : int PDGInt = EvtPDL::getStdHep(theId);
312 :
313 0 : if (abs(PDGInt) == _tauPDG) {
314 : // Delete any siblings for the tau particle
315 0 : if (theDaughter->getNDaug() > 0) {theDaughter->deleteDaughters(false);}
316 0 : tauMap[hepMCDaughter] = theDaughter;
317 0 : } else {
318 : // Treat all other particles as "stable"
319 0 : hepMCDaughter->set_status(Tauolapp::TauolaParticle::STABLE);
320 : }
321 :
322 0 : } // theDaughter != 0
323 : } // Loop over daughters
324 :
325 0 : } else {
326 :
327 : // We only have the one tau particle. Store only this in the map.
328 0 : HepMC::GenParticle* singleTau = this->createGenParticle(tauParticle);
329 0 : theVertex->add_particle_out(singleTau);
330 0 : tauMap[singleTau] = tauParticle;
331 :
332 0 : }
333 :
334 : // Now pass the event to Tauola for processing
335 : // Create a Tauola event object
336 0 : Tauolapp::TauolaHepMCEvent tauolaEvent(theEvent);
337 :
338 : // Run the Tauola algorithm
339 0 : tauolaEvent.decayTaus();
340 :
341 : // Loop over all tau particles in the HepMC event and create their EvtParticle daughters.
342 : // Store all final "stable" descendent particles as the tau daughters, i.e.
343 : // let Tauola decay any resonances such as a_1 or rho.
344 : // If there is more than one tau particle in the event, then also create the
345 : // corresponding EvtParticles for their daughters as well. They will not be
346 : // re-decayed since we check at the start of this function if the tau particle has
347 : // any daughters before running Tauola decayTaus().
348 :
349 0 : HepMC::GenEvent::particle_iterator eventIter;
350 0 : for (eventIter = theEvent->particles_begin(); eventIter != theEvent->particles_end();
351 0 : ++eventIter) {
352 :
353 : // Check to see if we have a tau particle
354 0 : HepMC::GenParticle* aParticle = (*eventIter);
355 :
356 0 : if (aParticle != 0 && abs(aParticle->pdg_id()) == _tauPDG) {
357 :
358 : // Find out what EvtParticle corresponds to the HepMC particle.
359 : // We need this to create and attach EvtParticle daughters.
360 0 : EvtParticle* tauEvtParticle = tauMap[aParticle];
361 :
362 0 : if (tauEvtParticle != 0) {
363 :
364 : // Get the tau 4-momentum in the lab (first mother) frame. We need to boost
365 : // all the tau daughters to this frame, such that daug.getP4() is in the tau restframe.
366 0 : EvtVector4R tauP4CM = tauEvtParticle->getP4Lab();
367 0 : tauP4CM.set(tauP4CM.get(0), -tauP4CM.get(1), -tauP4CM.get(2), -tauP4CM.get(3));
368 :
369 : // Get the decay vertex for the tau particle
370 0 : HepMC::GenVertex* endVertex = aParticle->end_vertex();
371 0 : HepMC::GenVertex::particle_iterator tauIter;
372 :
373 0 : std::vector<EvtId> daugIdVect;
374 0 : std::vector<EvtVector4R> daugP4Vect;
375 :
376 : // Loop through all descendants
377 0 : for (tauIter = endVertex->particles_begin(HepMC::descendants);
378 0 : tauIter != endVertex->particles_end(HepMC::descendants); ++tauIter) {
379 :
380 0 : HepMC::GenParticle* tauDaug = (*tauIter);
381 :
382 : // Check to see if this descendant has its own decay vertex, e.g. rho resonance.
383 : // If so, skip this daughter and continue looping through the descendant list
384 : // until we reach the final "stable" products (e.g. pi pi from rho -> pi pi).
385 0 : HepMC::GenVertex* daugDecayVtx = tauDaug->end_vertex();
386 0 : if (daugDecayVtx != 0) {continue;}
387 :
388 : // Store the particle id and 4-momentum
389 0 : int tauDaugPDG = tauDaug->pdg_id();
390 0 : EvtId daugId = EvtPDL::evtIdFromStdHep(tauDaugPDG);
391 0 : daugIdVect.push_back(daugId);
392 :
393 0 : HepMC::FourVector tauDaugP4 = tauDaug->momentum();
394 0 : double tauDaug_px = tauDaugP4.px();
395 0 : double tauDaug_py = tauDaugP4.py();
396 0 : double tauDaug_pz = tauDaugP4.pz();
397 0 : double tauDaug_E = tauDaugP4.e();
398 :
399 0 : EvtVector4R daugP4(tauDaug_E, tauDaug_px, tauDaug_py, tauDaug_pz);
400 0 : daugP4Vect.push_back(daugP4);
401 :
402 0 : } // Loop over HepMC tau daughters
403 :
404 : // Create the tau EvtParticle daughters and assign their ids and 4-mtm
405 0 : int nDaug = daugIdVect.size();
406 :
407 0 : tauEvtParticle->makeDaughters(nDaug, daugIdVect);
408 :
409 : int iDaug(0);
410 0 : for (iDaug = 0; iDaug < nDaug; iDaug++) {
411 :
412 0 : EvtParticle* theDaugPart = tauEvtParticle->getDaug(iDaug);
413 :
414 0 : if (theDaugPart != 0) {
415 0 : EvtId theDaugId = daugIdVect[iDaug];
416 0 : EvtVector4R theDaugP4 = daugP4Vect[iDaug];
417 0 : theDaugP4.applyBoostTo(tauP4CM); // Boost the 4-mtm to the tau rest frame
418 0 : theDaugPart->init(theDaugId, theDaugP4);
419 0 : }
420 :
421 : } // Loop over tau daughters
422 :
423 0 : }
424 :
425 0 : } // We have a tau HepMC particle in the event
426 :
427 0 : }
428 :
429 0 : theEvent->clear();
430 0 : delete theEvent;
431 :
432 0 : }
433 :
434 : HepMC::GenParticle* EvtTauolaEngine::createGenParticle(EvtParticle* theParticle) {
435 :
436 : // Method to create an HepMC::GenParticle version of the given EvtParticle.
437 0 : if (theParticle == 0) {return 0;}
438 :
439 : // Get the 4-momentum (E, px, py, pz) for the EvtParticle
440 0 : EvtVector4R p4 = theParticle->getP4Lab();
441 :
442 : // Convert this to the HepMC 4-momentum
443 0 : double E = p4.get(0);
444 0 : double px = p4.get(1);
445 0 : double py = p4.get(2);
446 0 : double pz = p4.get(3);
447 :
448 0 : HepMC::FourVector hepMC_p4(px, py, pz, E);
449 :
450 0 : int PDGInt = EvtPDL::getStdHep(theParticle->getId());
451 :
452 : // Set the status flag for the particle.
453 : int status = Tauolapp::TauolaParticle::HISTORY;
454 :
455 0 : HepMC::GenParticle* genParticle = new HepMC::GenParticle(hepMC_p4, PDGInt, status);
456 :
457 : return genParticle;
458 :
459 0 : }
460 :
461 : #endif
|
