File indexing completed on 2021-10-19 10:28:30
0001
0002
0003 #include <iostream>
0004
0005 #include "GeneratorInterface/TauolaInterface/interface/TauolappInterface.h"
0006
0007 #include "Tauola/Tauola.h"
0008 #include "Tauola/TauolaHepMCEvent.h"
0009 #include "Tauola/Log.h"
0010 #include "Tauola/TauolaHepMCParticle.h"
0011 #include "Tauola/TauolaParticle.h"
0012 #include "FWCore/ServiceRegistry/interface/Service.h"
0013 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0014 #include "FWCore/Utilities/interface/Exception.h"
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include "CLHEP/Random/RandomEngine.h"
0018
0019 #include "HepMC/GenEvent.h"
0020 #include "HepMC/IO_HEPEVT.h"
0021 #include "HepMC/HEPEVT_Wrapper.h"
0022
0023 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
0024
0025
0026 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
0027 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0028
0029
0030 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0031 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0032
0033 using namespace gen;
0034 using namespace edm;
0035 using namespace std;
0036
0037 CLHEP::HepRandomEngine* TauolappInterface::fRandomEngine = nullptr;
0038
0039 extern "C" {
0040
0041 void gen::ranmar_(float* rvec, int* lenv) {
0042 for (int i = 0; i < *lenv; i++)
0043 *rvec++ = TauolappInterface::flat();
0044 return;
0045 }
0046
0047 void gen::rmarin_(int*, int*, int*) { return; }
0048 }
0049
0050 TauolappInterface::TauolappInterface(const edm::ParameterSet& pset, edm::ConsumesCollector iCollector)
0051 : fPolarization(false),
0052 fPDGTableToken(iCollector.esConsumes<edm::Transition::BeginLuminosityBlock>()),
0053 fPSet(nullptr),
0054 fIsInitialized(false),
0055 fMDTAU(-1),
0056 fSelectDecayByEvent(false),
0057 lhe(nullptr),
0058 dmMatch(0.5),
0059 dolhe(false),
0060 dolheBosonCorr(false),
0061 ntries(10) {
0062 fPSet = new ParameterSet(pset);
0063 }
0064
0065 TauolappInterface::~TauolappInterface() {
0066 if (fPSet != nullptr)
0067 delete fPSet;
0068 }
0069
0070 void TauolappInterface::init(const edm::EventSetup& es) {
0071 if (fIsInitialized)
0072 return;
0073 if (fPSet == nullptr)
0074 throw cms::Exception("TauolappInterfaceError") << "Attempt to initialize Tauola with an empty ParameterSet\n"
0075 << std::endl;
0076
0077 fIsInitialized = true;
0078
0079 fPDGTable = es.getHandle(fPDGTableToken);
0080
0081 Tauolapp::Tauola::setDecayingParticle(15);
0082
0083
0084 dmMatch = fPSet->getUntrackedParameter<double>("dmMatch", 0.5);
0085 dolhe = fPSet->getUntrackedParameter<bool>("dolhe", false);
0086 dolheBosonCorr = fPSet->getUntrackedParameter<bool>("dolheBosonCorr", true);
0087 ntries = fPSet->getUntrackedParameter<int>("ntries", 10);
0088
0089
0090
0091 fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization");
0092
0093
0094
0095 ParameterSet cards = fPSet->getParameter<ParameterSet>("InputCards");
0096
0097 fMDTAU = cards.getParameter<int>("mdtau");
0098
0099 if (fMDTAU == 0 || fMDTAU == 1) {
0100 Tauolapp::Tauola::setSameParticleDecayMode(cards.getParameter<int>("pjak1"));
0101 Tauolapp::Tauola::setOppositeParticleDecayMode(cards.getParameter<int>("pjak2"));
0102 }
0103
0104 Tauolapp::Tauola::spin_correlation.setAll(fPolarization);
0105
0106
0107
0108
0109
0110 const HepPDT::ParticleData* PData =
0111 fPDGTable->particle(HepPDT::ParticleID(abs(Tauolapp::Tauola::getDecayingParticle())));
0112 double lifetime = PData->lifetime().value();
0113 Tauolapp::Tauola::setTauLifetime(lifetime);
0114
0115 fPDGs.push_back(Tauolapp::Tauola::getDecayingParticle());
0116
0117 Tauolapp::Tauola::setRandomGenerator(gen::TauolappInterface::flat);
0118
0119 if (fPSet->exists("parameterSets")) {
0120 std::vector<std::string> par = fPSet->getParameter<std::vector<std::string> >("parameterSets");
0121 for (unsigned int ip = 0; ip < par.size(); ++ip) {
0122 std::string curSet = par[ip];
0123 if (curSet == "setNewCurrents")
0124 Tauolapp::Tauola::setNewCurrents(fPSet->getParameter<int>(curSet));
0125 }
0126 }
0127
0128 Tauolapp::Tauola::initialize();
0129
0130 Tauolapp::Tauola::spin_correlation.setAll(
0131 fPolarization);
0132
0133 if (fPSet->exists("parameterSets")) {
0134 std::vector<std::string> par = fPSet->getParameter<std::vector<std::string> >("parameterSets");
0135 for (unsigned int ip = 0; ip < par.size(); ++ip) {
0136 std::string curSet = par[ip];
0137 if (curSet == "spinCorrelationSetAll")
0138 Tauolapp::Tauola::spin_correlation.setAll(fPSet->getParameter<bool>(curSet));
0139 if (curSet == "spinCorrelationGAMMA")
0140 Tauolapp::Tauola::spin_correlation.GAMMA = fPSet->getParameter<bool>(curSet);
0141 if (curSet == "spinCorrelationZ0")
0142 Tauolapp::Tauola::spin_correlation.Z0 = fPSet->getParameter<bool>(curSet);
0143 if (curSet == "spinCorrelationHIGGS")
0144 Tauolapp::Tauola::spin_correlation.HIGGS = fPSet->getParameter<bool>(curSet);
0145 if (curSet == "spinCorrelationHIGGSH")
0146 Tauolapp::Tauola::spin_correlation.HIGGS_H = fPSet->getParameter<bool>(curSet);
0147 if (curSet == "spinCorrelationHIGGSA")
0148 Tauolapp::Tauola::spin_correlation.HIGGS_A = fPSet->getParameter<bool>(curSet);
0149 if (curSet == "spinCorrelationHIGGSPLUS")
0150 Tauolapp::Tauola::spin_correlation.HIGGS_PLUS = fPSet->getParameter<bool>(curSet);
0151 if (curSet == "spinCorrelationHIGGSMINUS")
0152 Tauolapp::Tauola::spin_correlation.HIGGS_MINUS = fPSet->getParameter<bool>(curSet);
0153 if (curSet == "spinCorrelationWPLUS")
0154 Tauolapp::Tauola::spin_correlation.W_PLUS = fPSet->getParameter<bool>(curSet);
0155 if (curSet == "spinCorrelationWMINUS")
0156 Tauolapp::Tauola::spin_correlation.W_MINUS = fPSet->getParameter<bool>(curSet);
0157
0158 if (curSet == "setHiggsScalarPseudoscalarPDG")
0159 Tauolapp::Tauola::setHiggsScalarPseudoscalarPDG(fPSet->getParameter<int>(curSet));
0160 if (curSet == "setHiggsScalarPseudoscalarMixingAngle")
0161 Tauolapp::Tauola::setHiggsScalarPseudoscalarMixingAngle(fPSet->getParameter<double>(curSet));
0162
0163 if (curSet == "setRadiation")
0164 Tauolapp::Tauola::setRadiation(fPSet->getParameter<bool>(curSet));
0165 if (curSet == "setRadiationCutOff")
0166 Tauolapp::Tauola::setRadiationCutOff(fPSet->getParameter<double>(curSet));
0167
0168 if (curSet == "setEtaK0sPi") {
0169 std::vector<int> vpar = fPSet->getParameter<std::vector<int> >(curSet);
0170 if (vpar.size() == 3)
0171 Tauolapp::Tauola::setEtaK0sPi(vpar[0], vpar[1], vpar[2]);
0172 else {
0173 std::cout << "WARNING invalid size for setEtaK0sPi: " << vpar.size() << " Require 3 elements " << std::endl;
0174 }
0175 }
0176
0177 if (curSet == "setTaukle") {
0178 std::vector<double> vpar = fPSet->getParameter<std::vector<double> >(curSet);
0179 if (vpar.size() == 4)
0180 Tauolapp::Tauola::setTaukle(vpar[0], vpar[1], vpar[2], vpar[3]);
0181 else {
0182 std::cout << "WARNING invalid size for setTaukle: " << vpar.size() << " Require 4 elements " << std::endl;
0183 }
0184 }
0185
0186 if (curSet == "setTauBr") {
0187 edm::ParameterSet cfg = fPSet->getParameter<edm::ParameterSet>(curSet);
0188 std::vector<int> vJAK = cfg.getParameter<std::vector<int> >("JAK");
0189 std::vector<double> vBR = cfg.getParameter<std::vector<double> >("BR");
0190 if (vJAK.size() == vBR.size()) {
0191 for (unsigned int i = 0; i < vJAK.size(); i++)
0192 Tauolapp::Tauola::setTauBr(vJAK[i], vBR[i]);
0193 } else {
0194 std::cout << "WARNING invalid size for setTauBr - JAK: " << vJAK.size() << " BR: " << vBR.size() << std::endl;
0195 }
0196 }
0197 }
0198 }
0199
0200
0201
0202
0203
0204 if (fMDTAU != 0 && fMDTAU != 1) {
0205 decodeMDTAU(fMDTAU);
0206 }
0207
0208 Tauolapp::Log::LogWarning(false);
0209
0210 return;
0211 }
0212
0213 double TauolappInterface::flat() {
0214 if (!fRandomEngine) {
0215 throw cms::Exception("LogicError")
0216 << "TauolaInterface::flat: Attempt to generate random number when engine pointer is null\n"
0217 << "This might mean that the code was modified to generate a random number outside the\n"
0218 << "event and beginLuminosityBlock methods, which is not allowed.\n";
0219 }
0220 return fRandomEngine->flat();
0221 }
0222
0223 HepMC::GenEvent* TauolappInterface::decay(HepMC::GenEvent* evt) {
0224 if (!fIsInitialized)
0225 return evt;
0226 Tauolapp::Tauola::setRandomGenerator(
0227 gen::TauolappInterface::flat);
0228 int NPartBefore = evt->particles_size();
0229 int NVtxBefore = evt->vertices_size();
0230
0231
0232
0233
0234
0235
0236 if (fSelectDecayByEvent) {
0237 selectDecayByMDTAU();
0238 }
0239 if (dolhe && lhe != nullptr) {
0240 std::vector<HepMC::GenParticle> particles;
0241 std::vector<int> m_idx;
0242 std::vector<double> spinup = lhe->getHEPEUP()->SPINUP;
0243 std::vector<int> pdg = lhe->getHEPEUP()->IDUP;
0244 for (unsigned int i = 0; i < spinup.size(); i++) {
0245 particles.push_back(HepMC::GenParticle(HepMC::FourVector(lhe->getHEPEUP()->PUP.at(i)[0],
0246 lhe->getHEPEUP()->PUP.at(i)[1],
0247 lhe->getHEPEUP()->PUP.at(i)[2],
0248 lhe->getHEPEUP()->PUP.at(i)[3]),
0249 lhe->getHEPEUP()->IDUP.at(i)));
0250 int status = lhe->getHEPEUP()->ISTUP.at(i);
0251 particles.at(particles.size() - 1).set_generated_mass(lhe->getHEPEUP()->PUP.at(i)[4]);
0252 particles.at(particles.size() - 1).set_status(status > 0 ? (status == 2 ? 3 : status) : 3);
0253 m_idx.push_back(lhe->getHEPEUP()->MOTHUP.at(i).first - 1);
0254 }
0255
0256 bool hastaus(false);
0257 std::vector<HepMC::GenParticle*> match;
0258 for (HepMC::GenEvent::particle_const_iterator iter = evt->particles_begin(); iter != evt->particles_end(); iter++) {
0259 if (abs((*iter)->pdg_id()) == 15) {
0260 hastaus = true;
0261 int mother_pid(0);
0262
0263 for (HepMC::GenVertex::particle_iterator mother = (*iter)->production_vertex()->particles_begin(HepMC::parents);
0264 mother != (*iter)->production_vertex()->particles_end(HepMC::parents);
0265 mother++) {
0266 mother_pid = (*mother)->pdg_id();
0267 if (mother_pid != (*iter)->pdg_id()) {
0268
0269 if (abs(mother_pid) == 24 ||
0270 abs(mother_pid) == 37 ||
0271 abs(mother_pid) == 23 ||
0272 abs(mother_pid) == 22 ||
0273 abs(mother_pid) == 25 ||
0274 abs(mother_pid) == 35 ||
0275 abs(mother_pid) == 36
0276 ) {
0277 bool isfound = false;
0278 for (unsigned int k = 0; k < match.size(); k++) {
0279 if ((*mother) == match.at(k))
0280 isfound = true;
0281 }
0282 if (!isfound)
0283 match.push_back(*mother);
0284 }
0285 }
0286 }
0287 }
0288 }
0289 if (hastaus) {
0290
0291 if (match.size() == 1 && dolheBosonCorr) {
0292 for (int i = 0; i < ntries; i++) {
0293
0294 auto* t_event = new Tauolapp::TauolaHepMCEvent(evt);
0295 t_event->undecayTaus();
0296 t_event->decayTaus();
0297 bool ismatch = true;
0298 for (unsigned int j = 0; j < spinup.size(); j++) {
0299 if (abs(pdg.at(j)) == 15) {
0300 double diffhelminus = (-1.0 * (double)Tauolapp::Tauola::getHelMinus() -
0301 spinup.at(j));
0302 double diffhelplus = ((double)Tauolapp::Tauola::getHelPlus() - spinup.at(j));
0303 if (pdg.at(j) == 15 && diffhelminus > 0.5)
0304 ismatch = false;
0305 if (pdg.at(j) == -15 && diffhelplus > 0.5)
0306 ismatch = false;
0307 }
0308 }
0309 delete t_event;
0310 if (ismatch)
0311 break;
0312 }
0313 } else {
0314
0315
0316 auto* t_event = new Tauolapp::TauolaHepMCEvent(evt);
0317 t_event->undecayTaus();
0318 delete t_event;
0319
0320 for (HepMC::GenEvent::particle_const_iterator iter = evt->particles_begin(); iter != evt->particles_end();
0321 iter++) {
0322 if (abs((*iter)->pdg_id()) == 15 && isLastTauInChain(*iter)) {
0323 TLorentzVector ltau(
0324 (*iter)->momentum().px(), (*iter)->momentum().py(), (*iter)->momentum().pz(), (*iter)->momentum().e());
0325 HepMC::GenParticle* m = GetMother(*iter);
0326 TLorentzVector mother(m->momentum().px(), m->momentum().py(), m->momentum().pz(), m->momentum().e());
0327 TVector3 boost = -1.0 * mother.BoostVector();
0328 TLorentzVector ltau_lab = ltau;
0329 ltau.Boost(boost);
0330 mother.Boost(boost);
0331 HepMC::GenEvent* tauevt = make_simple_tau_event(ltau, (*iter)->pdg_id(), (*iter)->status());
0332 HepMC::GenParticle* p = (*(tauevt->particles_begin()));
0333 Tauolapp::TauolaParticle* tp = new Tauolapp::TauolaHepMCParticle(p);
0334 double helicity = MatchedLHESpinUp(*iter, particles, spinup, m_idx);
0335 if ((*iter)->pdg_id() == 15)
0336 helicity *= -1.0;
0337 tp->undecay();
0338
0339 Tauolapp::Tauola::decayOne(tp, true, 0, 0, ((double)helicity) * 0.999999);
0340 boost *= -1.0;
0341 mother.Boost(boost);
0342 update_particles((*iter), evt, p, boost);
0343
0344 BoostProdToLabLifeTimeInDecays((*iter), ltau_lab, ltau);
0345 delete tauevt;
0346 }
0347 }
0348 }
0349 }
0350 } else {
0351
0352 auto* t_event = new Tauolapp::TauolaHepMCEvent(evt);
0353
0354 t_event->decayTaus();
0355 delete t_event;
0356 }
0357
0358 for (int iv = NVtxBefore + 1; iv <= evt->vertices_size(); iv++) {
0359 HepMC::GenVertex* GenVtx = evt->barcode_to_vertex(-iv);
0360
0361
0362
0363
0364
0365 std::vector<int> BCodes;
0366 BCodes.clear();
0367 for (HepMC::GenVertex::particle_iterator pitr = GenVtx->particles_begin(HepMC::children);
0368 pitr != GenVtx->particles_end(HepMC::children);
0369 ++pitr) {
0370 if ((*pitr)->barcode() > 10000) {
0371 BCodes.push_back((*pitr)->barcode());
0372 }
0373 }
0374 if (!BCodes.empty()) {
0375 for (size_t ibc = 0; ibc < BCodes.size(); ibc++) {
0376 HepMC::GenParticle* p1 = evt->barcode_to_particle(BCodes[ibc]);
0377 int nbc = p1->barcode() - 10000 + NPartBefore;
0378 p1->suggest_barcode(nbc);
0379 }
0380 }
0381 }
0382
0383 for (HepMC::GenEvent::particle_const_iterator p = evt->particles_begin(); p != evt->particles_end(); ++p) {
0384 if ((*p)->end_vertex() && (*p)->status() == 1)
0385 (*p)->set_status(2);
0386 if ((*p)->end_vertex() && (*p)->end_vertex()->particles_out_size() == 0)
0387 edm::LogWarning("TauolappInterface::decay error: empty end vertex!");
0388 }
0389
0390 return evt;
0391 }
0392
0393 void TauolappInterface::statistics() { return; }
0394
0395 void TauolappInterface::decodeMDTAU(int mdtau) {
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410
0411
0412
0413
0414 if (mdtau == 101 || mdtau == 201) {
0415
0416
0417 Tauolapp::jaki_.jak1 = 1;
0418 Tauolapp::jaki_.jak2 = 1;
0419 Tauolapp::Tauola::setSameParticleDecayMode(1);
0420 Tauolapp::Tauola::setOppositeParticleDecayMode(1);
0421 return;
0422 }
0423
0424 if (mdtau == 102 || mdtau == 202) {
0425
0426
0427 Tauolapp::jaki_.jak1 = 2;
0428 Tauolapp::jaki_.jak2 = 2;
0429 Tauolapp::Tauola::setSameParticleDecayMode(2);
0430 Tauolapp::Tauola::setOppositeParticleDecayMode(2);
0431 return;
0432 }
0433
0434 if (mdtau == 111 || mdtau == 211) {
0435
0436
0437
0438 Tauolapp::jaki_.jak1 = 1;
0439 Tauolapp::jaki_.jak2 = 0;
0440 Tauolapp::Tauola::setSameParticleDecayMode(1);
0441 Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0442 return;
0443 }
0444
0445 if (mdtau == 112 || mdtau == 212) {
0446
0447
0448
0449 Tauolapp::jaki_.jak1 = 2;
0450 Tauolapp::jaki_.jak2 = 0;
0451 Tauolapp::Tauola::setSameParticleDecayMode(2);
0452 Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0453 return;
0454 }
0455
0456 if (mdtau == 121 || mdtau == 221) {
0457
0458
0459
0460 Tauolapp::jaki_.jak1 = 0;
0461 Tauolapp::jaki_.jak2 = 1;
0462 Tauolapp::Tauola::setSameParticleDecayMode(0);
0463 Tauolapp::Tauola::setOppositeParticleDecayMode(1);
0464 return;
0465 }
0466
0467 if (mdtau == 122 || mdtau == 222) {
0468
0469
0470
0471 Tauolapp::jaki_.jak1 = 0;
0472 Tauolapp::jaki_.jak2 = 2;
0473 Tauolapp::Tauola::setSameParticleDecayMode(0);
0474 Tauolapp::Tauola::setOppositeParticleDecayMode(2);
0475 return;
0476 }
0477
0478 if (mdtau == 140 || mdtau == 240) {
0479
0480
0481 Tauolapp::jaki_.jak1 = 3;
0482 Tauolapp::jaki_.jak2 = 3;
0483 Tauolapp::Tauola::setSameParticleDecayMode(3);
0484 Tauolapp::Tauola::setOppositeParticleDecayMode(3);
0485 return;
0486 }
0487
0488 if (mdtau == 141 || mdtau == 241) {
0489
0490
0491
0492 Tauolapp::jaki_.jak1 = 3;
0493 Tauolapp::jaki_.jak2 = 0;
0494 Tauolapp::Tauola::setSameParticleDecayMode(3);
0495 Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0496 return;
0497 }
0498
0499 if (mdtau == 142 || mdtau == 242) {
0500
0501
0502
0503 Tauolapp::jaki_.jak1 = 0;
0504 Tauolapp::jaki_.jak2 = 3;
0505 Tauolapp::Tauola::setSameParticleDecayMode(0);
0506 Tauolapp::Tauola::setOppositeParticleDecayMode(3);
0507 return;
0508 }
0509
0510
0511
0512
0513
0514
0515
0516
0517 double sumBra = 0;
0518
0519
0520
0521
0522
0523 for (int i = 0; i < 22; i++) {
0524 sumBra += Tauolapp::taubra_.gamprt[i];
0525 }
0526 if (sumBra == 0.)
0527 return;
0528 for (int i = 0; i < 22; i++) {
0529 double newBra = Tauolapp::taubra_.gamprt[i] / sumBra;
0530 Tauolapp::Tauola::setTauBr(i + 1, newBra);
0531 }
0532 sumBra = 1.0;
0533
0534 double sumLeptonBra = Tauolapp::taubra_.gamprt[0] + Tauolapp::taubra_.gamprt[1];
0535 double sumHadronBra = sumBra - sumLeptonBra;
0536
0537 for (int i = 0; i < 2; i++) {
0538 fLeptonModes.push_back(i + 1);
0539 fScaledLeptonBrRatios.push_back((Tauolapp::taubra_.gamprt[i] / sumLeptonBra));
0540 }
0541 for (int i = 2; i < 22; i++) {
0542 fHadronModes.push_back(i + 1);
0543 fScaledHadronBrRatios.push_back((Tauolapp::taubra_.gamprt[i] / sumHadronBra));
0544 }
0545
0546 fSelectDecayByEvent = true;
0547 return;
0548 }
0549
0550 void TauolappInterface::selectDecayByMDTAU() {
0551 if (fMDTAU == 100 || fMDTAU == 200) {
0552 int mode = selectLeptonic();
0553 Tauolapp::jaki_.jak1 = mode;
0554 Tauolapp::Tauola::setSameParticleDecayMode(mode);
0555 mode = selectLeptonic();
0556 Tauolapp::jaki_.jak2 = mode;
0557 Tauolapp::Tauola::setOppositeParticleDecayMode(mode);
0558 return;
0559 }
0560
0561 int modeL = selectLeptonic();
0562 int modeH = selectHadronic();
0563
0564 if (fMDTAU == 110 || fMDTAU == 210) {
0565 Tauolapp::jaki_.jak1 = modeL;
0566 Tauolapp::jaki_.jak2 = 0;
0567 Tauolapp::Tauola::setSameParticleDecayMode(modeL);
0568 Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0569 return;
0570 }
0571
0572 if (fMDTAU == 120 || fMDTAU == 22) {
0573 Tauolapp::jaki_.jak1 = 0;
0574 Tauolapp::jaki_.jak2 = modeL;
0575 Tauolapp::Tauola::setSameParticleDecayMode(0);
0576 Tauolapp::Tauola::setOppositeParticleDecayMode(modeL);
0577 return;
0578 }
0579
0580 if (fMDTAU == 114 || fMDTAU == 214) {
0581 Tauolapp::jaki_.jak1 = modeL;
0582 Tauolapp::jaki_.jak2 = modeH;
0583 Tauolapp::Tauola::setSameParticleDecayMode(modeL);
0584 Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0585 return;
0586 }
0587
0588 if (fMDTAU == 124 || fMDTAU == 224) {
0589 Tauolapp::jaki_.jak1 = modeH;
0590 Tauolapp::jaki_.jak2 = modeL;
0591 Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0592 Tauolapp::Tauola::setOppositeParticleDecayMode(modeL);
0593 return;
0594 }
0595
0596 if (fMDTAU == 115 || fMDTAU == 215) {
0597 Tauolapp::jaki_.jak1 = 1;
0598 Tauolapp::jaki_.jak2 = modeH;
0599 Tauolapp::Tauola::setSameParticleDecayMode(1);
0600 Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0601 return;
0602 }
0603
0604 if (fMDTAU == 125 || fMDTAU == 225) {
0605 Tauolapp::jaki_.jak1 = modeH;
0606 Tauolapp::jaki_.jak2 = 1;
0607 Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0608 Tauolapp::Tauola::setOppositeParticleDecayMode(1);
0609 return;
0610 }
0611
0612 if (fMDTAU == 116 || fMDTAU == 216) {
0613 Tauolapp::jaki_.jak1 = 2;
0614 Tauolapp::jaki_.jak2 = modeH;
0615 Tauolapp::Tauola::setSameParticleDecayMode(2);
0616 Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0617 return;
0618 }
0619
0620 if (fMDTAU == 126 || fMDTAU == 226) {
0621 Tauolapp::jaki_.jak1 = modeH;
0622 Tauolapp::jaki_.jak2 = 2;
0623 Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0624 Tauolapp::Tauola::setOppositeParticleDecayMode(2);
0625 return;
0626 }
0627
0628 if (fMDTAU == 130 || fMDTAU == 230) {
0629 Tauolapp::jaki_.jak1 = modeH;
0630 Tauolapp::jaki_.jak2 = selectHadronic();
0631 Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0632 Tauolapp::Tauola::setOppositeParticleDecayMode(Tauolapp::jaki_.jak2);
0633 return;
0634 }
0635
0636 if (fMDTAU == 131 || fMDTAU == 231) {
0637 Tauolapp::jaki_.jak1 = modeH;
0638 Tauolapp::jaki_.jak2 = 0;
0639 Tauolapp::Tauola::setSameParticleDecayMode(modeH);
0640 Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0641 return;
0642 }
0643
0644 if (fMDTAU == 132 || fMDTAU == 232) {
0645 Tauolapp::jaki_.jak1 = 0;
0646 Tauolapp::jaki_.jak2 = modeH;
0647 Tauolapp::Tauola::setSameParticleDecayMode(0);
0648 Tauolapp::Tauola::setOppositeParticleDecayMode(modeH);
0649 return;
0650 }
0651
0652
0653
0654
0655
0656
0657 Tauolapp::Tauola::setSameParticleDecayMode(0);
0658 Tauolapp::Tauola::setOppositeParticleDecayMode(0);
0659
0660 return;
0661 }
0662
0663 int TauolappInterface::selectLeptonic() {
0664 float prob = flat();
0665
0666 if (prob > 0. && prob <= fScaledLeptonBrRatios[0]) {
0667 return 1;
0668 } else if (prob > fScaledLeptonBrRatios[1] && prob <= 1.) {
0669 return 2;
0670 }
0671
0672 return 0;
0673 }
0674
0675 int TauolappInterface::selectHadronic() {
0676 float prob = 0.;
0677 int len = 1;
0678 ranmar_(&prob, &len);
0679
0680 double sumBra = fScaledHadronBrRatios[0];
0681 if (prob > 0. && prob <= sumBra) {
0682 return fHadronModes[0];
0683 } else {
0684 int NN = fScaledHadronBrRatios.size();
0685 for (int i = 1; i < NN; i++) {
0686 if (prob > sumBra && prob <= (sumBra + fScaledHadronBrRatios[i])) {
0687 return fHadronModes[i];
0688 }
0689 sumBra += fScaledHadronBrRatios[i];
0690 }
0691 }
0692
0693 return 0;
0694 }
0695
0696 HepMC::GenEvent* TauolappInterface::make_simple_tau_event(const TLorentzVector& l, int pdgid, int status) {
0697 HepMC::GenEvent* event = new HepMC::GenEvent();
0698
0699 HepMC::FourVector momentum_tau1(l.Px(), l.Py(), l.Pz(), l.E());
0700
0701 HepMC::GenParticle* tau1 = new HepMC::GenParticle(momentum_tau1, pdgid, status);
0702
0703 HepMC::GenVertex* vertex = new HepMC::GenVertex();
0704 vertex->add_particle_out(tau1);
0705 event->add_vertex(vertex);
0706 return event;
0707 }
0708
0709 void TauolappInterface::update_particles(HepMC::GenParticle* partHep,
0710 HepMC::GenEvent* theEvent,
0711 HepMC::GenParticle* p,
0712 TVector3& boost) {
0713 partHep->set_status(p->status());
0714 if (p->end_vertex()) {
0715 if (!partHep->end_vertex()) {
0716 HepMC::GenVertex* vtx = new HepMC::GenVertex(p->end_vertex()->position());
0717 theEvent->add_vertex(vtx);
0718 vtx->add_particle_in(partHep);
0719 }
0720 if (p->end_vertex()->particles_out_size() != 0) {
0721 for (HepMC::GenVertex::particles_out_const_iterator d = p->end_vertex()->particles_out_const_begin();
0722 d != p->end_vertex()->particles_out_const_end();
0723 d++) {
0724
0725 TLorentzVector l((*d)->momentum().px(), (*d)->momentum().py(), (*d)->momentum().pz(), (*d)->momentum().e());
0726 l.Boost(boost);
0727 HepMC::FourVector momentum(l.Px(), l.Py(), l.Pz(), l.E());
0728 HepMC::GenParticle* daughter = new HepMC::GenParticle(momentum, (*d)->pdg_id(), (*d)->status());
0729 daughter->suggest_barcode(theEvent->particles_size() + 1);
0730 partHep->end_vertex()->add_particle_out(daughter);
0731 if ((*d)->end_vertex())
0732 update_particles(daughter, theEvent, (*d), boost);
0733 }
0734 }
0735 }
0736 }
0737
0738 bool TauolappInterface::isLastTauInChain(const HepMC::GenParticle* tau) {
0739 if (tau->end_vertex()) {
0740 HepMC::GenVertex::particle_iterator dau;
0741 for (dau = tau->end_vertex()->particles_begin(HepMC::children);
0742 dau != tau->end_vertex()->particles_end(HepMC::children);
0743 dau++) {
0744 int dau_pid = (*dau)->pdg_id();
0745 if (dau_pid == tau->pdg_id())
0746 return false;
0747 }
0748 }
0749 return true;
0750 }
0751
0752 double TauolappInterface::MatchedLHESpinUp(HepMC::GenParticle* tau,
0753 std::vector<HepMC::GenParticle>& p,
0754 std::vector<double>& spinup,
0755 std::vector<int>& m_idx) {
0756 HepMC::GenParticle* Tau = FirstTauInChain(tau);
0757 HepMC::GenParticle* mother = GetMother(Tau);
0758 TLorentzVector t(tau->momentum().px(), tau->momentum().py(), tau->momentum().pz(), tau->momentum().e());
0759 TLorentzVector m(mother->momentum().px(), mother->momentum().py(), mother->momentum().pz(), mother->momentum().e());
0760 for (unsigned int i = 0; i < p.size(); i++) {
0761 if (tau->pdg_id() == p.at(i).pdg_id()) {
0762 if (mother->pdg_id() == p.at(m_idx.at(i)).pdg_id()) {
0763 TLorentzVector pm(p.at(m_idx.at(i)).momentum().px(),
0764 p.at(m_idx.at(i)).momentum().py(),
0765 p.at(m_idx.at(i)).momentum().pz(),
0766 p.at(m_idx.at(i)).momentum().e());
0767 if (fabs(m.M() - pm.M()) < dmMatch)
0768 return spinup.at(i);
0769 }
0770 }
0771 }
0772 return 0;
0773 }
0774
0775 HepMC::GenParticle* TauolappInterface::FirstTauInChain(HepMC::GenParticle* tau) {
0776 if (tau->production_vertex()) {
0777 HepMC::GenVertex::particle_iterator mother;
0778 for (mother = tau->production_vertex()->particles_begin(HepMC::parents);
0779 mother != tau->production_vertex()->particles_end(HepMC::parents);
0780 mother++) {
0781 if ((*mother)->pdg_id() == tau->pdg_id())
0782 return FirstTauInChain(*mother);
0783 }
0784 }
0785 return tau;
0786 }
0787
0788 HepMC::GenParticle* TauolappInterface::GetMother(HepMC::GenParticle* tau) {
0789 if (tau->production_vertex()) {
0790 HepMC::GenVertex::particle_iterator mother;
0791 for (mother = tau->production_vertex()->particles_begin(HepMC::parents);
0792 mother != tau->production_vertex()->particles_end(HepMC::parents);
0793 mother++) {
0794 if ((*mother)->pdg_id() == tau->pdg_id())
0795 return GetMother(*mother);
0796 return (*mother);
0797 }
0798 }
0799 return tau;
0800 }
0801
0802 void TauolappInterface::BoostProdToLabLifeTimeInDecays(HepMC::GenParticle* p,
0803 TLorentzVector& lab,
0804 TLorentzVector& prod) {
0805 if (p->end_vertex() && p->production_vertex()) {
0806 HepMC::GenVertex* PGenVtx = p->production_vertex();
0807 HepMC::GenVertex* EGenVtx = p->end_vertex();
0808 double VxDec = PGenVtx->position().x() + lab.Px() / prod.Px() * (EGenVtx->position().x() - PGenVtx->position().x());
0809 double VyDec = PGenVtx->position().y() + lab.Py() / prod.Py() * (EGenVtx->position().y() - PGenVtx->position().y());
0810 double VzDec = PGenVtx->position().z() + lab.Pz() / prod.Pz() * (EGenVtx->position().z() - PGenVtx->position().z());
0811 double VtDec = PGenVtx->position().t() + lab.Pt() / prod.Pt() * (EGenVtx->position().t() - PGenVtx->position().t());
0812 EGenVtx->set_position(HepMC::FourVector(VxDec, VyDec, VzDec, VtDec));
0813 for (HepMC::GenVertex::particle_iterator dau = p->end_vertex()->particles_begin(HepMC::children);
0814 dau != p->end_vertex()->particles_end(HepMC::children);
0815 dau++) {
0816 BoostProdToLabLifeTimeInDecays((*dau), lab, prod);
0817 }
0818 }
0819 }