CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TauolappInterface.cc
Go to the documentation of this file.
1 /* this is the code for the new Tauola++ */
2 
3 #include <iostream>
4 
6 
7 #include "Tauola/Tauola.h"
8 #include "Tauola/TauolaHepMCEvent.h"
9 #include "Tauola/Log.h"
10 #include "Tauola/TauolaHepMCParticle.h"
11 #include "Tauola/TauolaParticle.h"
15 
16 #include "CLHEP/Random/RandomEngine.h"
17 
18 #include "HepMC/GenEvent.h"
19 #include "HepMC/IO_HEPEVT.h"
20 #include "HepMC/HEPEVT_Wrapper.h"
21 
23 
24 // LHE Run
27 
28 // LHE Event
31 
32 using namespace gen;
33 using namespace edm;
34 using namespace std;
35 
36 CLHEP::HepRandomEngine* TauolappInterface::fRandomEngine = nullptr;
37 
38 extern "C" {
39 
40  void gen::ranmar_( float *rvec, int *lenv )
41  {
42  for(int i = 0; i < *lenv; i++)
43  *rvec++ = TauolappInterface::flat();
44  return;
45  }
46 
47  void gen::rmarin_( int*, int*, int* )
48  {
49  return;
50  }
51 
52 }
53 
55  fPolarization(false),
56  fPSet(0),
57  fIsInitialized(false),
58  fMDTAU(-1),
59  fSelectDecayByEvent(false),
60  lhe(NULL),
61  dmMatch(0.5),
62  dolhe(false),
63  dolheBosonCorr(false),
64  ntries(10)
65 {
66  if ( fPSet != 0 ) throw cms::Exception("TauolappInterfaceError") << "Attempt to override Tauola an existing ParameterSet\n" << std::endl;
67  fPSet = new ParameterSet(pset);
68 }
69 
71 
73  if ( fIsInitialized ) return; // do init only once
74  if ( fPSet == 0 ) throw cms::Exception("TauolappInterfaceError") << "Attempt to initialize Tauola with an empty ParameterSet\n" << std::endl;
75 
76  fIsInitialized = true;
77 
78  es.getData( fPDGTable ) ;
79 
80  Tauolapp::Tauola::setDecayingParticle(15);
81 
82  // LHE Information
83  dmMatch=fPSet->getUntrackedParameter<double>("dmMatch",0.5);
84  dolhe=fPSet->getUntrackedParameter<bool>("dolhe",false);
85  dolheBosonCorr=fPSet->getUntrackedParameter<bool>("dolheBosonCorr",true);
86  ntries=fPSet->getUntrackedParameter<int>("ntries",10);
87 
88  // polarization switch
89  // fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization") ? 1 : 0 ;
90  fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization");
91 
92  // read tau decay mode switches
93  //
94  ParameterSet cards = fPSet->getParameter< ParameterSet >("InputCards");
95 
96  fMDTAU = cards.getParameter< int >( "mdtau" );
97 
98  if ( fMDTAU == 0 || fMDTAU == 1 )
99  {
100  Tauolapp::Tauola::setSameParticleDecayMode( cards.getParameter< int >( "pjak1" ) ) ;
101  Tauolapp::Tauola::setOppositeParticleDecayMode( cards.getParameter< int >( "pjak2" ) ) ;
102  }
103 
104  Tauolapp::Tauola::spin_correlation.setAll(fPolarization);
105 
106  // some more options, copied over from an example
107  // Default values
108  //Tauola::setEtaK0sPi(0,0,0); // switches to decay eta K0_S and pi0 1/0 on/off.
109 
110  const HepPDT::ParticleData* PData = fPDGTable->particle(HepPDT::ParticleID( abs(Tauolapp::Tauola::getDecayingParticle()) )) ;
111  double lifetime = PData->lifetime().value();
112  Tauolapp::Tauola::setTauLifetime( lifetime );
113 
114  fPDGs.push_back( Tauolapp::Tauola::getDecayingParticle() );
115 
116  Tauolapp::Tauola::setRandomGenerator(gen::TauolappInterface::flat);
117 
118  if (fPSet->exists("parameterSets")){
119  std::vector<std::string> par = fPSet->getParameter< std::vector<std::string> >("parameterSets");
120  for (unsigned int ip=0; ip<par.size(); ++ip ){
121  std::string curSet = par[ip];
122  if(curSet=="setNewCurrents") Tauolapp::Tauola::setNewCurrents(fPSet->getParameter<int>(curSet));
123  }
124  }
125 
127 
128  Tauolapp::Tauola::spin_correlation.setAll(fPolarization);// Tauola switches this on during Tauola::initialise(); so we add this here to keep it on/off
129 
130  if (fPSet->exists("parameterSets")){
131  std::vector<std::string> par = fPSet->getParameter< std::vector<std::string> >("parameterSets");
132  for (unsigned int ip=0; ip<par.size(); ++ip ){
133  std::string curSet = par[ip];
134  if(curSet=="spinCorrelationSetAll") Tauolapp::Tauola::spin_correlation.setAll(fPSet->getParameter<bool>(curSet));
135  if(curSet=="spinCorrelationGAMMA") Tauolapp::Tauola::spin_correlation.GAMMA=fPSet->getParameter<bool>(curSet);
136  if(curSet=="spinCorrelationZ0") Tauolapp::Tauola::spin_correlation.Z0=fPSet->getParameter<bool>(curSet);
137  if(curSet=="spinCorrelationHIGGS") Tauolapp::Tauola::spin_correlation.HIGGS=fPSet->getParameter<bool>(curSet);
138  if(curSet=="spinCorrelationHIGGSH") Tauolapp::Tauola::spin_correlation.HIGGS_H=fPSet->getParameter<bool>(curSet);
139  if(curSet=="spinCorrelationHIGGSA") Tauolapp::Tauola::spin_correlation.HIGGS_A=fPSet->getParameter<bool>(curSet);
140  if(curSet=="spinCorrelationHIGGSPLUS") Tauolapp::Tauola::spin_correlation.HIGGS_PLUS=fPSet->getParameter<bool>(curSet);
141  if(curSet=="spinCorrelationHIGGSMINUS") Tauolapp::Tauola::spin_correlation.HIGGS_MINUS=fPSet->getParameter<bool>(curSet);
142  if(curSet=="spinCorrelationWPLUS") Tauolapp::Tauola::spin_correlation.W_PLUS=fPSet->getParameter<bool>(curSet);
143  if(curSet=="spinCorrelationWMINUS") Tauolapp::Tauola::spin_correlation.W_MINUS=fPSet->getParameter<bool>(curSet);
144 
145  if(curSet=="setHiggsScalarPseudoscalarPDG") Tauolapp::Tauola::setHiggsScalarPseudoscalarPDG(fPSet->getParameter<int>(curSet));
146  if(curSet=="setHiggsScalarPseudoscalarMixingAngle") Tauolapp::Tauola::setHiggsScalarPseudoscalarMixingAngle(fPSet->getParameter<double>(curSet));
147 
148  if(curSet=="setRadiation") Tauolapp::Tauola::setRadiation(fPSet->getParameter<bool>(curSet));
149  if(curSet=="setRadiationCutOff") Tauolapp::Tauola::setRadiationCutOff(fPSet->getParameter<double>(curSet));
150 
151  if(curSet=="setEtaK0sPi"){
152  std::vector<int> vpar = fPSet->getParameter<std::vector<int> >(curSet);
153  if(vpar.size()==3) Tauolapp::Tauola::setEtaK0sPi(vpar[0],vpar[1],vpar[2]);
154  else {std::cout << "WARNING invalid size for setEtaK0sPi: " << vpar.size() << " Require 3 elements " << std::endl;}
155  }
156 
157  if(curSet=="setTaukle"){
158  std::vector<double> vpar = fPSet->getParameter<std::vector<double> >(curSet);
159  if(vpar.size()==4) Tauolapp::Tauola::setTaukle(vpar[0], vpar[1], vpar[2], vpar[3]);
160  else {std::cout << "WARNING invalid size for setTaukle: " << vpar.size() << " Require 4 elements " << std::endl;}
161  }
162 
163  if(curSet=="setTauBr"){
165  std::vector<int> vJAK = cfg.getParameter<std::vector<int> >("JAK");
166  std::vector<double> vBR = cfg.getParameter<std::vector<double> >("BR");
167  if(vJAK.size() == vBR.size()){
168  for(unsigned int i=0;i<vJAK.size();i++) Tauolapp::Tauola::setTauBr(vJAK[i],vBR[i]);
169  }
170  else {std::cout << "WARNING invalid size for setTauBr - JAK: " << vJAK.size() << " BR: " << vBR.size() << std::endl;}
171  }
172  }
173  }
174 
175  // override decay modes if needs be
176  //
177  // we have to do it AFTER init because otherwises branching ratios are NOT filled in
178  //
179  if ( fMDTAU != 0 && fMDTAU != 1 )
180  {
181  decodeMDTAU( fMDTAU );
182  }
183 
184  Tauolapp::Log::LogWarning(false);
185 
186  return;
187 }
188 
190 {
191  if ( !fRandomEngine ) {
192  throw cms::Exception("LogicError")
193  << "TauolaInterface::flat: Attempt to generate random number when engine pointer is null\n"
194  << "This might mean that the code was modified to generate a random number outside the\n"
195  << "event and beginLuminosityBlock methods, which is not allowed.\n";
196  }
197  return fRandomEngine->flat();
198 }
199 
200 HepMC::GenEvent* TauolappInterface::decay( HepMC::GenEvent* evt ){
201  if ( !fIsInitialized ) return evt;
202  Tauolapp::Tauola::setRandomGenerator(gen::TauolappInterface::flat); // rest tauola++ random number incase other modules use tauola++
203  int NPartBefore = evt->particles_size();
204  int NVtxBefore = evt->vertices_size();
205 
206  // what do we do if Hep::GenEvent size is larger than 10K ???
207  // Tauola (& Photos, BTW) can only handle up to 10K via HEPEVT,
208  // and in case of CMS, it's only up to 4K !!!
209  // override decay mode if needs be
210 
213  }
214  if(dolhe && lhe!=NULL){
215  std::vector<HepMC::GenParticle> particles;
216  std::vector<int> m_idx;
217  std::vector<double> spinup=lhe->getHEPEUP()->SPINUP;
218  std::vector<int> pdg =lhe->getHEPEUP()->IDUP;
219  for(unsigned int i=0;i<spinup.size();i++){
220  particles.push_back(HepMC::GenParticle(HepMC::FourVector(lhe->getHEPEUP()->PUP.at(i)[0],lhe->getHEPEUP()->PUP.at(i)[1],lhe->getHEPEUP()->PUP.at(i)[2],lhe->getHEPEUP()->PUP.at(i)[3]),lhe->getHEPEUP()->IDUP.at(i)));
221  int status = lhe->getHEPEUP()->ISTUP.at(i);
222  particles.at(particles.size()-1).set_generated_mass(lhe->getHEPEUP()->PUP.at(i)[4]);
223  particles.at(particles.size()-1).set_status(status > 0 ? (status == 2 ? 3 : status) : 3);
224  m_idx.push_back(lhe->getHEPEUP()->MOTHUP.at(i).first-1);// correct for fortran index offset
225  }
226  // match to taus in hepmc and identify mother of taus
227  bool hastaus(false);
228  std::vector<HepMC::GenParticle*> match;
229  for(HepMC::GenEvent::particle_const_iterator iter = evt->particles_begin(); iter != evt->particles_end(); iter++) {
230  if(abs((*iter)->pdg_id())==15){
231  hastaus=true;
232  int mother_pid(0);
233  // check imediate parent to avoid parent tau ie tau->taugamma
234  for(HepMC::GenVertex::particle_iterator mother=(*iter)->production_vertex()->particles_begin(HepMC::parents); mother!=(*iter)->production_vertex()->particles_end(HepMC::parents); mother++) {
235  mother_pid = (*mother)->pdg_id();
236  if(mother_pid!=(*iter)->pdg_id()){
237  // match against lhe record
238  if(abs(mother_pid) == 24 || // W
239  abs(mother_pid) == 37 || // H+/-
240  abs(mother_pid) == 23 || // Z
241  abs(mother_pid) == 22 || // gamma
242  abs(mother_pid) == 25 || // H0 SM
243  abs(mother_pid) == 35 || // H0
244  abs(mother_pid) == 36 // A0
245  ){
246  bool isfound=false;
247  for(unsigned int k=0;k<match.size();k++){
248  if((*mother)==match.at(k))isfound=true;
249  }
250  if(!isfound) match.push_back(*mother);
251  }
252  }
253  }
254  }
255  }
256  if(hastaus){
257  // if is single gauge boson decay and match helicities
258  if(match.size()==1 && dolheBosonCorr){
259  for(int i=0; i<ntries;i++){
260  // re-decay taus then check if helicities match
261  auto * t_event = new Tauolapp::TauolaHepMCEvent(evt);
262  t_event->undecayTaus();
263  t_event->decayTaus();
264  bool ismatch=true;
265  for(unsigned int j=0;j<spinup.size();j++){
266  if(abs(pdg.at(j))==15){
267  double diffhelminus=(-1.0*(double)Tauolapp::Tauola::getHelMinus()-spinup.at(j)); // -1.0 to correct for tauola feature
268  double diffhelplus=((double)Tauolapp::Tauola::getHelPlus()-spinup.at(j));
269  if(pdg.at(j)==15 && diffhelminus>0.5) ismatch=false;
270  if(pdg.at(j)==-15 && diffhelplus>0.5) ismatch=false;
271  }
272  }
273  delete t_event;
274  if(ismatch) break;
275  }
276  }
277  else{
278  // If the event does not contain a single gauge boson the code will be run with
279  // remove all tau decays
280  auto * t_event = new Tauolapp::TauolaHepMCEvent(evt);
281  t_event->undecayTaus();
282  delete t_event;
283  // decay all taus manually based on the helicity
284  for(HepMC::GenEvent::particle_const_iterator iter = evt->particles_begin(); iter != evt->particles_end(); iter++) {
285  if(abs((*iter)->pdg_id())==15 && isLastTauInChain(*iter)){
286  TLorentzVector ltau((*iter)->momentum().px(),(*iter)->momentum().py(),(*iter)->momentum().pz(),(*iter)->momentum().e());
288  TLorentzVector mother(m->momentum().px(),m->momentum().py(),m->momentum().pz(),m->momentum().e());
289  TVector3 boost=-1.0*mother.BoostVector(); // boost into mother's CM frame
290  TLorentzVector ltau_lab=ltau;
291  ltau.Boost(boost);
292  mother.Boost(boost);
293  HepMC::GenEvent *tauevt=make_simple_tau_event(ltau,(*iter)->pdg_id(),(*iter)->status());
294  HepMC::GenParticle *p=(*(tauevt->particles_begin()));
295  Tauolapp::TauolaParticle *tp = new Tauolapp::TauolaHepMCParticle(p);
296  double helicity=MatchedLHESpinUp(*iter,particles,spinup,m_idx); // get helicity from lhe
297  if((*iter)->pdg_id()==15) helicity*=-1.0;
298  tp->undecay();
299  // use |S_{tau}|=0.999999 to avoid issues with numerical roundoff
300  Tauolapp::Tauola::decayOne(tp,true,0,0,((double)helicity)*0.999999);
301  boost*=-1.0; // boost back to lab frame
302  mother.Boost(boost);
303  update_particles((*iter),evt,p,boost);
304  //correct tau liftetime for boost (change rest frame from mothers to taus)
305  BoostProdToLabLifeTimeInDecays((*iter),ltau_lab,ltau);
306  delete tauevt;
307  }
308  }
309  }
310  }
311  }
312  else{
313  //construct tmp TAUOLA event
314  auto * t_event = new Tauolapp::TauolaHepMCEvent(evt);
315  //t_event->undecayTaus();
316  t_event->decayTaus();
317  delete t_event;
318  }
319 
320  for ( int iv=NVtxBefore+1; iv<=evt->vertices_size(); iv++ ){
321  HepMC::GenVertex* GenVtx = evt->barcode_to_vertex(-iv);
322  //
323  // now find decay products with funky barcode, weed out and replace with clones of sensible barcode
324  // we can NOT change the barcode while iterating, because iterators do depend on the barcoding
325  // thus we have to take a 2-step procedure
326  //
327  std::vector<int> BCodes;
328  BCodes.clear();
329  for (HepMC::GenVertex::particle_iterator pitr= GenVtx->particles_begin(HepMC::children);
330  pitr != GenVtx->particles_end(HepMC::children); ++pitr){
331  if ( (*pitr)->barcode() > 10000 ){
332  BCodes.push_back( (*pitr)->barcode() );
333  }
334  }
335  if ( BCodes.size() > 0 ){
336  for ( size_t ibc=0; ibc<BCodes.size(); ibc++ ){
337  HepMC::GenParticle* p1 = evt->barcode_to_particle( BCodes[ibc] );
338  int nbc = p1->barcode() - 10000 + NPartBefore;
339  p1->suggest_barcode( nbc );
340  }
341  }
342  }
343  return evt;
344 }
345 
347 {
348  return;
349 }
350 
352 {
353 
354  // Note-1:
355  // I have to hack the common block directly because set<...>DecayMode(...)
356  // only changes it in the Tauola++ instance but does NOT passes it over
357  // to the Fortran core - this it does only one, via initialize() stuff...
358  //
359  // So I'll do both ways of settings, just for consistency...
360  // but I probably need to communicate it to the Tauola(++) team...
361  //
362 
363  // Note-2:
364  // originally, the 1xx settings are meant for tau's from hard event,
365  // and the 2xx settings are for any tau in the event record;
366  //
367  // later one, we'll have to take this into account...
368  // but first I'll have to sort out what happens in the 1xx case
369  // to tau's coming outside of hard event (if any in the record)
370  //
371 
372  if ( mdtau == 101 || mdtau == 201 )
373  {
374  // override with electron mode for both tau's
375  //
376  Tauolapp::jaki_.jak1 = 1;
377  Tauolapp::jaki_.jak2 = 1;
378  Tauolapp::Tauola::setSameParticleDecayMode( 1 ) ;
379  Tauolapp::Tauola::setOppositeParticleDecayMode( 1 ) ;
380  return;
381  }
382 
383  if ( mdtau == 102 || mdtau == 202 )
384  {
385  // override with muon mode for both tau's
386  //
387  Tauolapp::jaki_.jak1 = 2;
388  Tauolapp::jaki_.jak2 = 2;
389  Tauolapp::Tauola::setSameParticleDecayMode( 2 ) ;
390  Tauolapp::Tauola::setOppositeParticleDecayMode( 2 ) ;
391  return;
392  }
393 
394  if ( mdtau == 111 || mdtau == 211 )
395  {
396  // override with electron mode for 1st tau
397  // and any mode for 2nd tau
398  //
399  Tauolapp::jaki_.jak1 = 1;
400  Tauolapp::jaki_.jak2 = 0;
401  Tauolapp::Tauola::setSameParticleDecayMode( 1 ) ;
402  Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
403  return;
404  }
405 
406  if ( mdtau == 112 || mdtau == 212 )
407  {
408  // override with muon mode for the 1st tau
409  // and any mode for the 2nd tau
410  //
411  Tauolapp::jaki_.jak1 = 2;
412  Tauolapp::jaki_.jak2 = 0;
413  Tauolapp::Tauola::setSameParticleDecayMode( 2 ) ;
414  Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
415  return;
416  }
417 
418  if ( mdtau == 121 || mdtau == 221 )
419  {
420  // override with any mode for the 1st tau
421  // and electron mode for the 2nd tau
422  //
423  Tauolapp::jaki_.jak1 = 0;
424  Tauolapp::jaki_.jak2 = 1;
425  Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
426  Tauolapp::Tauola::setOppositeParticleDecayMode( 1 ) ;
427  return;
428  }
429 
430  if ( mdtau == 122 || mdtau == 222 )
431  {
432  // override with any mode for the 1st tau
433  // and muon mode for the 2nd tau
434  //
435  Tauolapp::jaki_.jak1 = 0;
436  Tauolapp::jaki_.jak2 = 2;
437  Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
438  Tauolapp::Tauola::setOppositeParticleDecayMode( 2 ) ;
439  return;
440  }
441 
442  if ( mdtau == 140 || mdtau == 240 )
443  {
444  // override with pi+/- nutau mode for both tau's
445  //
446  Tauolapp::jaki_.jak1 = 3;
447  Tauolapp::jaki_.jak2 = 3;
448  Tauolapp::Tauola::setSameParticleDecayMode( 3 ) ;
449  Tauolapp::Tauola::setOppositeParticleDecayMode( 3 ) ;
450  return;
451  }
452 
453  if ( mdtau == 141 || mdtau == 241 )
454  {
455  // override with pi+/- nutau mode for the 1st tau
456  // and any mode for the 2nd tau
457  //
458  Tauolapp::jaki_.jak1 = 3;
459  Tauolapp::jaki_.jak2 = 0;
460  Tauolapp::Tauola::setSameParticleDecayMode( 3 ) ;
461  Tauolapp::Tauola::setOppositeParticleDecayMode( 0 ) ;
462  return;
463  }
464 
465  if ( mdtau == 142 || mdtau == 242 )
466  {
467  // override with any mode for the 1st tau
468  // and pi+/- nutau mode for 2nd tau
469  //
470  Tauolapp::jaki_.jak1 = 0;
471  Tauolapp::jaki_.jak2 = 3;
472  Tauolapp::Tauola::setSameParticleDecayMode( 0 ) ;
473  Tauolapp::Tauola::setOppositeParticleDecayMode( 3 ) ;
474  return;
475  }
476 
477  // OK, we come here for semi-inclusive modes
478  //
479 
480  // First of all, leptons and hadron modes sums
481  //
482  // re-scale branching ratios, just in case...
483  //
484  double sumBra = 0;
485 
486  // the number of decay modes is hardcoded at 22 because that's what it is right now in Tauola
487  // in the future, perhaps an asscess method would be useful - communicate to Tauola team...
488  //
489 
490  for ( int i=0; i<22; i++ )
491  {
492  sumBra += Tauolapp::taubra_.gamprt[i];
493  }
494  if ( sumBra == 0. ) return ; // perhaps need to throw ?
495  for ( int i=0; i<22; i++ )
496  {
497  double newBra = Tauolapp::taubra_.gamprt[i] / sumBra;
498  Tauolapp::Tauola::setTauBr( i+1, newBra );
499  }
500  sumBra = 1.0;
501 
502  double sumLeptonBra = Tauolapp::taubra_.gamprt[0] + Tauolapp::taubra_.gamprt[1];
503  double sumHadronBra = sumBra - sumLeptonBra;
504 
505  for ( int i=0; i<2; i++ )
506  {
507  fLeptonModes.push_back( i+1 );
508  fScaledLeptonBrRatios.push_back( (Tauolapp::taubra_.gamprt[i]/sumLeptonBra) );
509  }
510  for ( int i=2; i<22; i++ )
511  {
512  fHadronModes.push_back( i+1 );
513  fScaledHadronBrRatios.push_back( (Tauolapp::taubra_.gamprt[i]/sumHadronBra) );
514  }
515 
516  fSelectDecayByEvent = true;
517  return;
518 
519 }
520 
522 {
523 
524 
525  if ( fMDTAU == 100 || fMDTAU == 200 )
526  {
527  int mode = selectLeptonic();
528  Tauolapp::jaki_.jak1 = mode;
529  Tauolapp::Tauola::setSameParticleDecayMode( mode );
530  mode = selectLeptonic();
531  Tauolapp::jaki_.jak2 = mode;
532  Tauolapp::Tauola::setOppositeParticleDecayMode( mode );
533  return ;
534  }
535 
536  int modeL = selectLeptonic();
537  int modeH = selectHadronic();
538 
539  if ( fMDTAU == 110 || fMDTAU == 210 )
540  {
541  Tauolapp::jaki_.jak1 = modeL;
542  Tauolapp::jaki_.jak2 = 0;
543  Tauolapp::Tauola::setSameParticleDecayMode( modeL );
544  Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
545  return ;
546  }
547 
548  if ( fMDTAU == 120 || fMDTAU == 22 )
549  {
550  Tauolapp::jaki_.jak1 = 0;
551  Tauolapp::jaki_.jak2 = modeL;
552  Tauolapp::Tauola::setSameParticleDecayMode( 0 );
553  Tauolapp::Tauola::setOppositeParticleDecayMode( modeL );
554  return;
555  }
556 
557  if ( fMDTAU == 114 || fMDTAU == 214 )
558  {
559  Tauolapp::jaki_.jak1 = modeL;
560  Tauolapp::jaki_.jak2 = modeH;
561  Tauolapp::Tauola::setSameParticleDecayMode( modeL );
562  Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
563  return;
564  }
565 
566  if ( fMDTAU == 124 || fMDTAU == 224 )
567  {
568  Tauolapp::jaki_.jak1 = modeH;
569  Tauolapp::jaki_.jak2 = modeL;
570  Tauolapp::Tauola::setSameParticleDecayMode( modeH );
571  Tauolapp::Tauola::setOppositeParticleDecayMode( modeL );
572  return;
573  }
574 
575  if ( fMDTAU == 115 || fMDTAU == 215 )
576  {
577  Tauolapp::jaki_.jak1 = 1;
578  Tauolapp::jaki_.jak2 = modeH;
579  Tauolapp::Tauola::setSameParticleDecayMode( 1 );
580  Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
581  return;
582  }
583 
584  if ( fMDTAU == 125 || fMDTAU == 225 )
585  {
586  Tauolapp::jaki_.jak1 = modeH;
587  Tauolapp::jaki_.jak2 = 1;
588  Tauolapp::Tauola::setSameParticleDecayMode( modeH );
589  Tauolapp::Tauola::setOppositeParticleDecayMode( 1 );
590  return;
591  }
592 
593  if ( fMDTAU == 116 || fMDTAU == 216 )
594  {
595  Tauolapp::jaki_.jak1 = 2;
596  Tauolapp::jaki_.jak2 = modeH;
597  Tauolapp::Tauola::setSameParticleDecayMode( 2 );
598  Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
599  return;
600  }
601 
602  if ( fMDTAU == 126 || fMDTAU == 226 )
603  {
604  Tauolapp::jaki_.jak1 = modeH;
605  Tauolapp::jaki_.jak2 = 2;
606  Tauolapp::Tauola::setSameParticleDecayMode( modeH );
607  Tauolapp::Tauola::setOppositeParticleDecayMode( 2 );
608  return;
609  }
610 
611  if ( fMDTAU == 130 || fMDTAU == 230 )
612  {
613  Tauolapp::jaki_.jak1 = modeH;
614  Tauolapp::jaki_.jak2 = selectHadronic();
615  Tauolapp::Tauola::setSameParticleDecayMode( modeH );
616  Tauolapp::Tauola::setOppositeParticleDecayMode( Tauolapp::jaki_.jak2 );
617  return;
618  }
619 
620  if ( fMDTAU == 131 || fMDTAU == 231 )
621  {
622  Tauolapp::jaki_.jak1 = modeH;
623  Tauolapp::jaki_.jak2 = 0;
624  Tauolapp::Tauola::setSameParticleDecayMode( modeH );
625  Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
626  return;
627  }
628 
629  if ( fMDTAU == 132 || fMDTAU == 232 )
630  {
631  Tauolapp::jaki_.jak1 = 0;
632  Tauolapp::jaki_.jak2 = modeH;
633  Tauolapp::Tauola::setSameParticleDecayMode( 0 );
634  Tauolapp::Tauola::setOppositeParticleDecayMode( modeH );
635  return;
636  }
637 
638  // unlikely that we get here on unknown mdtau
639  // - there's a protection earlier
640  // but if we do, just set defaults
641  // probably need to spit a warning...
642  //
643  Tauolapp::Tauola::setSameParticleDecayMode( 0 );
644  Tauolapp::Tauola::setOppositeParticleDecayMode( 0 );
645 
646  return;
647 
648 
649 }
650 
652 {
653 
654  float prob = flat();
655 
656  if ( prob > 0. && prob <= fScaledLeptonBrRatios[0] )
657  {
658  return 1;
659  }
660  else if ( prob > fScaledLeptonBrRatios[1] && prob <=1. )
661  {
662  return 2;
663  }
664 
665  return 0;
666 }
667 
669 {
670 
671  float prob = 0.;
672  int len = 1;
673  ranmar_(&prob,&len);
674 
675  double sumBra = fScaledHadronBrRatios[0];
676  if ( prob > 0. && prob <= sumBra )
677  {
678  return fHadronModes[0];
679  }
680  else
681  {
682  int NN = fScaledHadronBrRatios.size();
683  for ( int i=1; i<NN; i++ )
684  {
685  if ( prob > sumBra && prob <= (sumBra+fScaledHadronBrRatios[i]) )
686  {
687  return fHadronModes[i];
688  }
689  sumBra += fScaledHadronBrRatios[i];
690  }
691  }
692 
693  return 0;
694 
695 }
696 
697 HepMC::GenEvent* TauolappInterface::make_simple_tau_event(const TLorentzVector &l,int pdgid,int status){
698  HepMC::GenEvent *event = new HepMC::GenEvent();
699  // make tau's four vector
700  HepMC::FourVector momentum_tau1(l.Px(),l.Py(),l.Pz(),l.E());
701  // make particles
702  HepMC::GenParticle *tau1 = new HepMC::GenParticle(momentum_tau1,pdgid,status);
703  // make the vertex
704  HepMC::GenVertex *vertex = new HepMC::GenVertex();
705  vertex->add_particle_out(tau1);
706  event->add_vertex(vertex);
707  return event;
708 }
709 
710 void TauolappInterface::update_particles(HepMC::GenParticle* partHep,HepMC::GenEvent* theEvent,HepMC::GenParticle* p,TVector3 &boost){
711  partHep->set_status(p->status());
712  if(p->end_vertex()){
713  if(!partHep->end_vertex()){
714  HepMC::GenVertex* vtx = new HepMC::GenVertex(p->end_vertex()->position());
715  theEvent->add_vertex(vtx);
716  vtx->add_particle_in(partHep);
717  }
718  if(p->end_vertex()->particles_out_size()!=0){
719  for(HepMC::GenVertex::particles_out_const_iterator d=p->end_vertex()->particles_out_const_begin(); d!=p->end_vertex()->particles_out_const_end();d++){
720  // Create daughter and add to event
721  TLorentzVector l((*d)->momentum().px(),(*d)->momentum().py(),(*d)->momentum().pz(),(*d)->momentum().e());
722  l.Boost(boost);
723  HepMC::FourVector momentum(l.Px(),l.Py(),l.Pz(),l.E());
724  HepMC::GenParticle *daughter = new HepMC::GenParticle(momentum,(*d)->pdg_id(),(*d)->status());
725  daughter->suggest_barcode(theEvent->particles_size()+1);
726  partHep->end_vertex()->add_particle_out(daughter);
727  if((*d)->end_vertex()) update_particles(daughter,theEvent,(*d),boost);
728  }
729  }
730  }
731 }
732 
734  if ( tau->end_vertex() ) {
735  HepMC::GenVertex::particle_iterator dau;
736  for (dau = tau->end_vertex()->particles_begin(HepMC::children); dau!= tau->end_vertex()->particles_end(HepMC::children); dau++ ) {
737  int dau_pid = (*dau)->pdg_id();
738  if(dau_pid == tau->pdg_id()) return false;
739  }
740  }
741  return true;
742 }
743 
744 double TauolappInterface::MatchedLHESpinUp(HepMC::GenParticle* tau, std::vector<HepMC::GenParticle> &p, std::vector<double> &spinup,
745  std::vector<int> &m_idx){
747  HepMC::GenParticle* mother=GetMother(Tau);
748  TLorentzVector t(tau->momentum().px(),tau->momentum().py(),tau->momentum().pz(),tau->momentum().e());
749  TLorentzVector m(mother->momentum().px(),mother->momentum().py(),mother->momentum().pz(),mother->momentum().e());
750  for(unsigned int i=0;i<p.size();i++){
751  if(tau->pdg_id()==p.at(i).pdg_id()){
752  if(mother->pdg_id()==p.at(m_idx.at(i)).pdg_id()){
753  TLorentzVector pm(p.at(m_idx.at(i)).momentum().px(),p.at(m_idx.at(i)).momentum().py(),p.at(m_idx.at(i)).momentum().pz(),p.at(m_idx.at(i)).momentum().e());
754  if(fabs(m.M()-pm.M())<dmMatch)return spinup.at(i);
755  }
756  }
757  }
758  return 0;
759 }
760 
762  if ( tau->production_vertex() ) {
763  HepMC::GenVertex::particle_iterator mother;
764  for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
765  if((*mother)->pdg_id()==tau->pdg_id()) return FirstTauInChain(*mother); // recursive call to get mother with different pdgid
766  }
767  }
768  return tau;
769 }
770 
772  if ( tau->production_vertex() ) {
773  HepMC::GenVertex::particle_iterator mother;
774  for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
775  if((*mother)->pdg_id() == tau->pdg_id()) return GetMother(*mother); // recursive call to get mother with different pdgid
776  return (*mother);
777  }
778  }
779  return tau;
780 }
781 
783  if(p->end_vertex() && p->production_vertex()){
784  HepMC::GenVertex* PGenVtx=p->production_vertex();
785  HepMC::GenVertex* EGenVtx=p->end_vertex();
786  double VxDec = PGenVtx->position().x()+lab.Px()/prod.Px()*(EGenVtx->position().x()-PGenVtx->position().x());
787  double VyDec = PGenVtx->position().y()+lab.Py()/prod.Py()*(EGenVtx->position().y()-PGenVtx->position().y());
788  double VzDec = PGenVtx->position().z()+lab.Pz()/prod.Pz()*(EGenVtx->position().z()-PGenVtx->position().z());
789  double VtDec = PGenVtx->position().t()+lab.Pt()/prod.Pt()*(EGenVtx->position().t()-PGenVtx->position().t());
790  EGenVtx->set_position(HepMC::FourVector(VxDec,VyDec,VzDec,VtDec));
791  for(HepMC::GenVertex::particle_iterator dau=p->end_vertex()->particles_begin(HepMC::children); dau!=p->end_vertex()->particles_end(HepMC::children); dau++){
792  BoostProdToLabLifeTimeInDecays((*dau),lab,prod); //recursively modify everything in the decay chain
793  }
794  }
795 }
T getParameter(std::string const &) const
static AlgebraicMatrix initialize()
TauolappInterface(const edm::ParameterSet &)
T getUntrackedParameter(std::string const &, T const &) const
int i
Definition: DBlmapReader.cc:9
edm::ParameterSet * fPSet
TPRegexp parents
Definition: eve_filter.cc:24
tuple cfg
Definition: looper.py:237
std::vector< double > fScaledHadronBrRatios
bool exists(std::string const &parameterName) const
checks if a parameter exists
#define NULL
Definition: scimark2.h:8
HepMC::GenEvent * make_simple_tau_event(const TLorentzVector &l, int pdgid, int status)
void rmarin_(int *, int *, int *)
const HEPEUP * getHEPEUP() const
Definition: LHEEvent.h:43
void getData(T &iHolder) const
Definition: EventSetup.h:78
float gamprt[30]
Definition: TauolaWrapper.h:78
return((rh^lh)&mask)
int jak2
Definition: TauolaWrapper.h:67
double p[5][pyjets_maxn]
std::vector< std::pair< int, int > > MOTHUP
Definition: LesHouches.h:236
std::vector< int > fPDGs
void BoostProdToLabLifeTimeInDecays(HepMC::GenParticle *p, TLorentzVector &lab, TLorentzVector &prod)
std::vector< int > fLeptonModes
std::vector< FiveVector > PUP
Definition: LesHouches.h:248
struct @582 taubra_
std::vector< double > SPINUP
Definition: LesHouches.h:261
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
int j
Definition: DBlmapReader.cc:9
std::vector< int > ISTUP
Definition: LesHouches.h:230
void ranmar_(float *rvec, int *lenv)
How EventSelector::AcceptEvent() decides whether to accept an event for output otherwise it is excluding the probing of A single or multiple positive and the trigger will pass if any such matching triggers are PASS or EXCEPTION[A criterion thatmatches no triggers at all is detected and causes a throw.] A single negative with an expectation of appropriate bit checking in the decision and the trigger will pass if any such matching triggers are FAIL or EXCEPTION A wildcarded negative criterion that matches more than one trigger in the trigger but the state exists so we define the behavior If all triggers are the negative crieriion will lead to accepting the event(this again matches the behavior of"!*"before the partial wildcard feature was incorporated).The per-event"cost"of each negative criterion with multiple relevant triggers is about the same as!*was in the past
HepPDT::ParticleData ParticleData
std::vector< int > IDUP
Definition: LesHouches.h:225
int k[5][pyjets_maxn]
void update_particles(HepMC::GenParticle *partHep, HepMC::GenEvent *theEvent, HepMC::GenParticle *p, TVector3 &boost)
void init(const edm::EventSetup &)
HepMC::GenParticle * FirstTauInChain(HepMC::GenParticle *tau)
HepMC::GenEvent * decay(HepMC::GenEvent *)
bool isLastTauInChain(const HepMC::GenParticle *tau)
std::vector< double > fScaledLeptonBrRatios
static CLHEP::HepRandomEngine * fRandomEngine
std::vector< int > fHadronModes
double p1[4]
Definition: TauolaWrapper.h:89
HepMC::GenParticle * GetMother(HepMC::GenParticle *tau)
tuple cout
Definition: gather_cfg.py:121
std::pair< typename Association::data_type::first_type, double > match(Reference key, Association association, bool bestMatchByMaxValue)
Generic matching function.
Definition: Utils.h:6
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
volatile std::atomic< bool > shutdown_flag false
tuple status
Definition: ntuplemaker.py:245
double MatchedLHESpinUp(HepMC::GenParticle *tau, std::vector< HepMC::GenParticle > &p, std::vector< double > &spinup, std::vector< int > &m_idx)
int mdtau
Definition: TauolaWrapper.h:59