CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TauolaInterface.cc
Go to the documentation of this file.
1 /* for old tauola27
2 #include <iostream>
3 
4 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h"
5 
6 #include "GeneratorInterface/ExternalDecays/interface/TauolaInterface.h"
7 #include "GeneratorInterface/ExternalDecays/interface/TauolaWrapper.h"
8 
9 #include "FWCore/ServiceRegistry/interface/Service.h"
10 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
11 
12 #include "HepMC/GenEvent.h"
13 #include "HepMC/IO_HEPEVT.h"
14 #include "HepMC/HEPEVT_Wrapper.h"
15 
16 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
17 
18 #include "GeneratorInterface/ExternalDecays/interface/DecayRandomEngine.h"
19 
20 using namespace gen;
21 using namespace edm;
22 using namespace std;
23 
24 extern "C" {
25 
26  void ranmar_( float *rvec, int *lenv )
27  {
28 
29  for(int i = 0; i < *lenv; i++)
30  *rvec++ = decayRandomEngine->flat();
31 
32  return;
33 
34  }
35 
36  void rmarin_( int*, int*, int* )
37  {
38 
39  return;
40 
41  }
42 
43 }
44 
45 //
46 // General Note: While there're no explicit calls or otherwise "links" to Pythia6 anywhere,
47 // we're using Pythia6Service here because we run pretauola rather than "core" tauola;
48 // pretauola is an extension on top of tauola, which is tied to Pythia6 via several routines;
49 // most heavily use one is PYR - we can't avoid it (other Pythia6-tied routines we avoid)
50 //
51 
52 TauolaInterface::TauolaInterface( const ParameterSet& pset )
53  : fIsInitialized(false)
54 {
55  fPy6Service = new Pythia6Service;
56 
57  fPolarization = pset.getParameter<bool>("UseTauolaPolarization") ? 1 : 0 ;
58 
59  // set Tauola defaults
60  //
61  ki_taumod_.pjak1 = -1;
62  ki_taumod_.pjak2 = -1;
63  ki_taumod_.mdtau = -1;
64 
65  // read tau decay mode switches
66  //
67  ParameterSet cards = pset.getParameter< ParameterSet >("InputCards");
68  ki_taumod_.pjak1 = cards.getParameter< int >( "pjak1" );
69  ki_taumod_.pjak2 = cards.getParameter< int >( "pjak2" );
70  ki_taumod_.mdtau = cards.getParameter< int >( "mdtau" );
71 
72 }
73 
74 TauolaInterface::~TauolaInterface()
75 {
76  delete fPy6Service;
77 }
78 
79 void TauolaInterface::init( const edm::EventSetup& es )
80 {
81 
82  if ( fIsInitialized ) return; // do init only once
83 
84  if ( ki_taumod_.mdtau <= -1 ) // actually, need to throw exception !
85  return ;
86 
87  fPDGs.push_back( 15 ) ;
88  es.getData( fPDGTable ) ;
89 
90  cout << "----------------------------------------------" << endl;
91  cout << "Initializing Tauola" << endl;
92  if ( fPolarization == 0 )
93  {
94  cout << "Tauola: Polarization disabled" << endl;
95  }
96  else if ( fPolarization == 1 )
97  {
98  cout << "Tauola: Polarization enabled" << endl;
99  }
100 
101 // FIXME !!!
102 // This is a temporary hack - we're re-using master generator's seed to init RANMAR
103 // FIXME !!!
104 // This is now off because ranmar has been overriden (see code above) to use pyr_(...)
105 // - this way we're using guaranteed initialized rndm generator... BUT !!! in the long
106 // run we may want a separate random stream for tauola...
107 
108 // Service<RandomNumberGenerator> rng;
109 // int seed = rng->mySeed() ;
110 // int ntot=0, ntot2=0;
111 // rmarin_( &seed, &ntot, &ntot2 );
112 
113  int mode = -2;
114  taurep_( &mode ) ;
115  mode = -1;
116  // tauola_( &mode, &fPolarization );
117  // tauola_srs_( &mode, &fPolarization );
118  //
119  // We're using the call(...) method here because it'll make sure that Py6
120  // is initialized, and that's done only once, and will grab exatly that instance
121  //
122  fPy6Service->call( tauola_srs_, &mode, &fPolarization );
123 
124  fIsInitialized = true;
125 
126  return;
127 }
128 
129 HepMC::GenEvent* TauolaInterface::decay( HepMC::GenEvent* evt )
130 {
131 
132  // event record convertor
133  //
134  HepMC::IO_HEPEVT conv;
135 
136  if ( !fIsInitialized ) return conv.read_next_event();
137 
138  // We are using random numbers, we are fetched through Pythia6Service
139  // (through ranmar_ below) -> so grab the instance during decay()
140 
141  Pythia6Service::InstanceWrapper pythia6InstanceGuard( fPy6Service );
142 
143  // fill up HEPEVT common block
144  //
145  // IDEALLY, this should be the way to go
146  // BUT !!! this utility fills it up in the "reshuffled" order,
147  // and later on Tauola chocks on it
148  //
149  // Needs to be sorted out, eith in HepMC, or in Tauola, or both !!!
150  //
151  // At present, this thing blindly relies on the assumption that
152  // HEPEVT is always there - which wont be the case with Py8 or Hwg++
153  //
154  //HepMC::IO_HEPEVT conv;
155  //conv.write_event( evt ) ;
156 
157  int numPartBeforeTauola = HepMC::HEPEVT_Wrapper::number_entries();
158  // HepMC::HEPEVT_Wrapper::print_hepevt();
159 
160  int mode = 0;
161  // tauola_( &mode, &fPolarization );
162  fPy6Service->call( tauola_srs_, &mode, &fPolarization );
163 
164  int numPartAfterTauola = HepMC::HEPEVT_Wrapper::number_entries();
165  // HepMC::HEPEVT_Wrapper::print_hepevt();
166 
167  // before we do the conversion, we need to deal with decay vertexes
168  // since Tauola knows nothing about lifetimes, all decay vertexes are set to 0.
169  // nees to set them properly, knowing lifetime !
170  // here we do it on HEPEVT record, also for consistency, although it's probably
171  // even easier to deal with HepMC::GenEvent record
172 
173  // find 1st "non-doc" tau
174  //
175  bool foundTau = false;
176  for ( int ip=1; ip<=numPartAfterTauola; ip++ )
177  {
178  if ( std::abs( HepMC::HEPEVT_Wrapper::id( ip ) ) == 15
179  && HepMC::HEPEVT_Wrapper::status( ip ) != 3 )
180  {
181  foundTau = true;
182  break;
183  }
184  }
185 
186  if ( !foundTau )
187  {
188  // no tau found
189  // just give up here
190  //
191  return conv.read_next_event();
192  }
193 
194  std::vector<int> PrntIndx;
195 
196  for ( int ip=numPartAfterTauola; ip>numPartBeforeTauola; ip-- ) // Fortran indexing !
197  {
198 
199  // first of all, find out how many generations in decay chain
200  //
201  PrntIndx.clear();
202  int Prnt = HepMC::HEPEVT_Wrapper::first_parent(ip);
203  ip -= (HepMC::HEPEVT_Wrapper::number_children(Prnt)-1); // such that we don't go the same part again
204  PrntIndx.push_back( Prnt );
205  while ( abs( HepMC::HEPEVT_Wrapper::id(Prnt) ) != 15 ) // shortcut; need to loop over fPDGs...
206  {
207  int Prnt1 = HepMC::HEPEVT_Wrapper::first_parent( Prnt );
208  Prnt = Prnt1;
209  // such that the tau always appear at the start of the list
210  PrntIndx.insert( PrntIndx.begin(), Prnt );
211  ip -= HepMC::HEPEVT_Wrapper::number_children(Prnt); // such that we don't go the same part again
212  }
213  for ( size_t iprt=0; iprt<PrntIndx.size(); iprt++ )
214  {
215  int Indx = PrntIndx[iprt];
216  int PartID = HepMC::HEPEVT_Wrapper::id( Indx );
217  const HepPDT::ParticleData*
218  PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID))) ;
219  //
220  // prob = exp(-t/lifetime) ==> t = -lifetime * log(prob)
221  //
222  float prob = 0.;
223  int length=1;
224  ranmar_(&prob,&length);
225  double lifetime = PData->lifetime().value();
226  //
227  // in case of Py6, this would be copied into V(5,i)
228  // for HEPEVT, need to check...
229  //
230  double ct = -lifetime * std::log(prob);
231  //
232  double ee = HepMC::HEPEVT_Wrapper::e( Indx );
233  double px = HepMC::HEPEVT_Wrapper::px( Indx );
234  double py = HepMC::HEPEVT_Wrapper::py( Indx );
235  double pz = HepMC::HEPEVT_Wrapper::pz( Indx );
236  // double pp = std::sqrt( px*px + py*py + pz*pz );
237  double mass = HepMC::HEPEVT_Wrapper::m( Indx );
238  //
239  // this is in py6 terms:
240  // VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5)
241  //
242  double VxDec = HepMC::HEPEVT_Wrapper::x( Indx );
243  VxDec += ct * (px/mass);
244  double VyDec = HepMC::HEPEVT_Wrapper::y( Indx );
245  VyDec += ct * (py/mass);
246  double VzDec = HepMC::HEPEVT_Wrapper::z( Indx );
247  VzDec += ct * (pz/mass);
248  double VtDec = HepMC::HEPEVT_Wrapper::t( Indx );
249  VtDec += ct * (ee/mass);
250  for ( int idau=HepMC::HEPEVT_Wrapper::first_child( Indx );
251  idau<=HepMC::HEPEVT_Wrapper::last_child( Indx ); idau++ )
252  {
253  HepMC::HEPEVT_Wrapper::set_position( idau, VxDec, VyDec, VzDec, VtDec );
254  }
255  }
256  }
257 
258  return conv.read_next_event();
259 
260 }
261 
262 void TauolaInterface::statistics()
263 {
264  int mode = 1;
265  // tauola_( &mode, &fPolarization ) ;
266  // tauola_srs_( &mode, &fPolarization ) ;
267  fPy6Service->call( tauola_srs_, &mode, &fPolarization );
268  return;
269 }
270 
271 */
272 
273 /* this is the code for the new Tauola++ */
274 
275 #include <iostream>
276 
278 
279 #include "Tauola.h"
280 #include "TauolaHepMCEvent.h"
281 #include "Log.h"
282 
286 
287 #include "CLHEP/Random/RandomEngine.h"
288 
289 #include "HepMC/GenEvent.h"
290 #include "HepMC/IO_HEPEVT.h"
291 #include "HepMC/HEPEVT_Wrapper.h"
292 
294 
295 // #include "GeneratorInterface/ExternalDecays/interface/DecayRandomEngine.h"
296 
297 
298 extern "C" {
299 
300  void gen::ranmar_( float *rvec, int *lenv )
301  {
302  TauolaInterface* instance = TauolaInterface::getInstance();
303  for(int i = 0; i < *lenv; i++)
304  // *rvec++ = decayRandomEngine->flat();
305  *rvec++ = instance->flat();
306  return;
307  }
308 
309  void gen::rmarin_( int*, int*, int* )
310  {
311  return;
312  }
313 
314 }
315 
316 using namespace gen;
317 using namespace edm;
318 using namespace std;
319 
321 
322 
324  : fPolarization(false), fPSet(0), fIsInitialized(false), fMDTAU(-1), fSelectDecayByEvent(false)
325 {
326 
328  if(!rng.isAvailable()) {
329  throw cms::Exception("Configuration")
330  << "The RandomNumberProducer module requires the RandomNumberGeneratorService\n"
331  "which appears to be absent. Please add that service to your configuration\n"
332  "or remove the modules that require it." << std::endl;
333  }
334 
335  fRandomEngine = &rng->getEngine();
336 
337 }
338 
339 
340 //TauolaInterface::TauolaInterface( const ParameterSet& pset )
341 // : fIsInitialized(false)
342 //{
343 //
344 // Tauola::setDecayingParticle(15);
345 // // --> ??? Tauola::setRadiation(false);
346 //
347 // // polarization switch
348 // //
349 // // fPolarization = pset.getParameter<bool>("UseTauolaPolarization") ? 1 : 0 ;
350 // fPolarization = pset.getParameter<bool>("UseTauolaPolarization");
351 //
352 // // read tau decay mode switches
353 // //
354 // ParameterSet cards = pset.getParameter< ParameterSet >("InputCards");
355 // Tauola::setSameParticleDecayMode( cards.getParameter< int >( "pjak1" ) ) ;
356 // Tauola::setOppositeParticleDecayMode( cards.getParameter< int >( "pjak2" ) ) ;
357 //
358 // Tauola::setTauLifetime(0.0);
359 // Tauola::spin_correlation.setAll(fPolarization);
360 //
361 // // some more options, copied over from an example
362 // // - maybe will use later...
363 // //
364 // //Tauola::setEtaK0sPi(0,0,0); // switches to decay eta K0_S and pi0 1/0 on/off.
365 // //
366 //
367 //}
368 
369 
371 {
372 
373  if ( fInstance == 0 ) fInstance = new TauolaInterface() ;
374  return fInstance;
375 
376 }
377 
378 
380 {
381 
382  if ( fPSet != 0 ) delete fPSet;
383  if ( fInstance == this ) fInstance = 0;
384 
385 }
386 
388 {
389 
390  if ( fPSet != 0 )
391  {
392  throw cms::Exception("TauolaInterfaceError")
393  << "Attempt to override Tauola an existing ParameterSet\n"
394  << std::endl;
395  }
396 
397  fPSet = new ParameterSet(pset);
398 
399  return;
400 
401 }
402 
404 {
405 
406  if ( fIsInitialized ) return; // do init only once
407 
408  if ( fPSet == 0 )
409  {
410 
411  throw cms::Exception("TauolaInterfaceError")
412  << "Attempt to initialize Tauola with an empty ParameterSet\n"
413  << std::endl;
414  }
415 
416  fIsInitialized = true;
417 
418  es.getData( fPDGTable ) ;
419 
420  Tauola::setDecayingParticle(15);
421  // --> ??? Tauola::setRadiation(false);
422 
423  // polarization switch
424  //
425  // fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization") ? 1 : 0 ;
426  fPolarization = fPSet->getParameter<bool>("UseTauolaPolarization");
427 
428  // read tau decay mode switches
429  //
430  ParameterSet cards = fPSet->getParameter< ParameterSet >("InputCards");
431 
432  fMDTAU = cards.getParameter< int >( "mdtau" );
433 
434  if ( fMDTAU == 0 || fMDTAU == 1 )
435  {
436  Tauola::setSameParticleDecayMode( cards.getParameter< int >( "pjak1" ) ) ;
437  Tauola::setOppositeParticleDecayMode( cards.getParameter< int >( "pjak2" ) ) ;
438  }
439 
440  Tauola::setTauLifetime(0.0);
441  Tauola::spin_correlation.setAll(fPolarization);
442 
443  // some more options, copied over from an example
444  // - maybe will use later...
445  //
446  //Tauola::setEtaK0sPi(0,0,0); // switches to decay eta K0_S and pi0 1/0 on/off.
447  //
448 
449 //
450 // const HepPDT::ParticleData*
451 // PData = fPDGTable->particle(HepPDT::ParticleID( abs(Tauola::getDecayingParticle()) )) ;
452 // double lifetime = PData->lifetime().value();
453 // Tauola::setTauLifetime( lifetime );
454 
455  fPDGs.push_back( Tauola::getDecayingParticle() );
456 
457  Tauola::setRandomGenerator(&gen::TauolappInterface_RandGetter);
458  Tauola::initialize();
459 
460  Tauola::spin_correlation.setAll(fPolarization);// Tauola switches this on during Tauola::initialise(); so we add this here to keep it on/off
461 
462  // override decay modes if needs be
463  //
464  // we have to do it AFTER init because otherwises branching ratios are NOT filled in
465  //
466  if ( fMDTAU != 0 && fMDTAU != 1 )
467  {
468  decodeMDTAU( fMDTAU );
469  }
470 
471  Log::LogWarning(false);
472 
473  return;
474 }
475 
477 {
478 
479  if ( !fPSet )
480  {
481  // throw
482  throw cms::Exception("TauolaInterfaceError")
483  << "Attempt to run random number generator of un-initialized Tauola\n"
484  << std::endl;
485  }
486 
487  if ( !fIsInitialized )
488  {
489  // throw
490  throw cms::Exception("TauolaInterfaceError")
491  << "Attempt to run random number generator of un-initialized Tauola\n"
492  << std::endl;
493  }
494 
495  return fRandomEngine->flat();
496 
497 }
498 
499 HepMC::GenEvent* TauolaInterface::decay( HepMC::GenEvent* evt )
500 {
501 
502  if ( !fIsInitialized ) return evt;
503 
504  int NPartBefore = evt->particles_size();
505  int NVtxBefore = evt->vertices_size();
506 
507  // what do we do if Hep::GenEvent size is larger than 10K ???
508  // Tauola (& Photos, BTW) can only handle up to 10K via HEPEVT,
509  // and in case of CMS, it's only up to 4K !!!
510  //
511  // if ( NPartBefore > 10000 ) return evt;
512  //
513 
514  // override decay mode if needs be
515  if ( fSelectDecayByEvent )
516  {
518  }
519 
520  //construct tmp TAUOLA event
521  //
522  TauolaHepMCEvent * t_event = new TauolaHepMCEvent(evt);
523 
524  // another option: if one lets Pythia or another master gen to decay taus,
525  // we have to undecay them first
526  // t_event->undecayTaus();
527 
528  // run Tauola on the tmp event - HepMC::GenEvernt will be MODIFIED !!!
529  //
530  t_event->decayTaus();
531 
532  // delet tmp Tauola event
533  //
534  delete t_event;
535 
536  // do we also need to apply the lifetime and vtx position shift ???
537  // (see TauolaInterface, for example)
538  //
539  // NOTE: the procedure ASSYMES that vertex barcoding is COUNTIUOUS/SEQUENTIAL,
540  // and that the abs(barcode) corresponds to vertex "plain indexing"
541  //
542  for ( int iv=NVtxBefore+1; iv<=evt->vertices_size(); iv++ )
543  {
544  HepMC::GenVertex* GenVtx = evt->barcode_to_vertex(-iv);
545  HepMC::GenParticle* GenPart = *(GenVtx->particles_in_const_begin());
546  HepMC::GenVertex* ProdVtx = GenPart->production_vertex();
547  HepMC::FourVector PMom = GenPart->momentum();
548  double mass = GenPart->generated_mass();
549  const HepPDT::ParticleData*
550  PData = fPDGTable->particle(HepPDT::ParticleID(abs(GenPart->pdg_id()))) ;
551  double lifetime = PData->lifetime().value();
552  float prob = 0.;
553  int length=1;
554  ranmar_(&prob,&length);
555  double ct = -lifetime * std::log(prob);
556  double VxDec = GenVtx->position().x();
557  VxDec += ct * (PMom.px()/mass);
558  VxDec += ProdVtx->position().x();
559  double VyDec = GenVtx->position().y();
560  VyDec += ct * (PMom.py()/mass);
561  VyDec += ProdVtx->position().y();
562  double VzDec = GenVtx->position().z();
563  VzDec += ct * (PMom.pz()/mass);
564  VzDec += ProdVtx->position().z();
565  double VtDec = GenVtx->position().t();
566  VtDec += ct * (PMom.e()/mass);
567  VtDec += ProdVtx->position().t();
568  GenVtx->set_position( HepMC::FourVector(VxDec,VyDec,VzDec,VtDec) );
569  //
570  // now find decay products with funky barcode, weed out and replace with clones of sensible barcode
571  // we can NOT change the barcode while iterating, because iterators do depend on the barcoding
572  // thus we have to take a 2-step procedure
573  //
574  std::vector<int> BCodes;
575  BCodes.clear();
576  for (HepMC::GenVertex::particle_iterator pitr= GenVtx->particles_begin(HepMC::children);
577  pitr != GenVtx->particles_end(HepMC::children); ++pitr)
578  {
579  if ( (*pitr)->barcode() > 10000 )
580  {
581  BCodes.push_back( (*pitr)->barcode() );
582  }
583  }
584  if ( BCodes.size() > 0 )
585  {
586  for ( size_t ibc=0; ibc<BCodes.size(); ibc++ )
587  {
588  HepMC::GenParticle* p1 = evt->barcode_to_particle( BCodes[ibc] );
589  int nbc = p1->barcode() - 10000 + NPartBefore;
590  p1->suggest_barcode( nbc );
591  }
592  }
593  }
594 
595  return evt;
596 
597 }
598 
600 {
601  return;
602 }
603 
605 {
606 
607  // Note-1:
608  // I have to hack the common block directly because set<...>DecayMode(...)
609  // only changes it in the Tauola++ instance but does NOT passes it over
610  // to the Fortran core - this it does only one, via initialize() stuff...
611  //
612  // So I'll do both ways of settings, just for consistency...
613  // but I probably need to communicate it to the Tauola(++) team...
614  //
615 
616  // Note-2:
617  // originally, the 1xx settings are meant for tau's from hard event,
618  // and the 2xx settings are for any tau in the event record;
619  //
620  // later one, we'll have to take this into account...
621  // but first I'll have to sort out what happens in the 1xx case
622  // to tau's coming outside of hard event (if any in the record)
623  //
624 
625  if ( mdtau == 101 || mdtau == 201 )
626  {
627  // override with electron mode for both tau's
628  //
629  jaki_.jak1 = 1;
630  jaki_.jak2 = 1;
631  Tauola::setSameParticleDecayMode( 1 ) ;
632  Tauola::setOppositeParticleDecayMode( 1 ) ;
633  return;
634  }
635 
636  if ( mdtau == 102 || mdtau == 202 )
637  {
638  // override with muon mode for both tau's
639  //
640  jaki_.jak1 = 2;
641  jaki_.jak2 = 2;
642  Tauola::setSameParticleDecayMode( 2 ) ;
643  Tauola::setOppositeParticleDecayMode( 2 ) ;
644  return;
645  }
646 
647  if ( mdtau == 111 || mdtau == 211 )
648  {
649  // override with electron mode for 1st tau
650  // and any mode for 2nd tau
651  //
652  jaki_.jak1 = 1;
653  jaki_.jak2 = 0;
654  Tauola::setSameParticleDecayMode( 1 ) ;
655  Tauola::setOppositeParticleDecayMode( 0 ) ;
656  return;
657  }
658 
659  if ( mdtau == 112 || mdtau == 212 )
660  {
661  // override with muon mode for the 1st tau
662  // and any mode for the 2nd tau
663  //
664  jaki_.jak1 = 2;
665  jaki_.jak2 = 0;
666  Tauola::setSameParticleDecayMode( 2 ) ;
667  Tauola::setOppositeParticleDecayMode( 0 ) ;
668  return;
669  }
670 
671  if ( mdtau == 121 || mdtau == 221 )
672  {
673  // override with any mode for the 1st tau
674  // and electron mode for the 2nd tau
675  //
676  jaki_.jak1 = 0;
677  jaki_.jak2 = 1;
678  Tauola::setSameParticleDecayMode( 0 ) ;
679  Tauola::setOppositeParticleDecayMode( 1 ) ;
680  return;
681  }
682 
683  if ( mdtau == 122 || mdtau == 222 )
684  {
685  // override with any mode for the 1st tau
686  // and muon mode for the 2nd tau
687  //
688  jaki_.jak1 = 0;
689  jaki_.jak2 = 2;
690  Tauola::setSameParticleDecayMode( 0 ) ;
691  Tauola::setOppositeParticleDecayMode( 2 ) ;
692  return;
693  }
694 
695  if ( mdtau == 140 || mdtau == 240 )
696  {
697  // override with pi+/- nutau mode for both tau's
698  //
699  jaki_.jak1 = 3;
700  jaki_.jak2 = 3;
701  Tauola::setSameParticleDecayMode( 3 ) ;
702  Tauola::setOppositeParticleDecayMode( 3 ) ;
703  return;
704  }
705 
706  if ( mdtau == 141 || mdtau == 241 )
707  {
708  // override with pi+/- nutau mode for the 1st tau
709  // and any mode for the 2nd tau
710  //
711  jaki_.jak1 = 3;
712  jaki_.jak2 = 0;
713  Tauola::setSameParticleDecayMode( 3 ) ;
714  Tauola::setOppositeParticleDecayMode( 0 ) ;
715  return;
716  }
717 
718  if ( mdtau == 142 || mdtau == 242 )
719  {
720  // override with any mode for the 1st tau
721  // and pi+/- nutau mode for 2nd tau
722  //
723  jaki_.jak1 = 0;
724  jaki_.jak2 = 3;
725  Tauola::setSameParticleDecayMode( 0 ) ;
726  Tauola::setOppositeParticleDecayMode( 3 ) ;
727  return;
728  }
729 
730  // OK, we come here for semi-inclusive modes
731  //
732 
733  // First of all, leptons and hadron modes sums
734  //
735  // re-scale branching ratios, just in case...
736  //
737  double sumBra = 0;
738 
739  // the number of decay modes is hardcoded at 22 because that's what it is right now in Tauola
740  // in the future, perhaps an asscess method would be useful - communicate to Tauola team...
741  //
742 
743  for ( int i=0; i<22; i++ )
744  {
745  sumBra += taubra_.gamprt[i];
746  }
747  if ( sumBra == 0. ) return ; // perhaps need to throw ?
748  for ( int i=0; i<22; i++ )
749  {
750  double newBra = taubra_.gamprt[i] / sumBra;
751  Tauola::setTauBr( i+1, newBra );
752  }
753  sumBra = 1.0;
754 
755  double sumLeptonBra = taubra_.gamprt[0] + taubra_.gamprt[1];
756  double sumHadronBra = sumBra - sumLeptonBra;
757 
758  for ( int i=0; i<2; i++ )
759  {
760  fLeptonModes.push_back( i+1 );
761  fScaledLeptonBrRatios.push_back( (taubra_.gamprt[i]/sumLeptonBra) );
762  }
763  for ( int i=2; i<22; i++ )
764  {
765  fHadronModes.push_back( i+1 );
766  fScaledHadronBrRatios.push_back( (taubra_.gamprt[i]/sumHadronBra) );
767  }
768 
769  fSelectDecayByEvent = true;
770  return;
771 
772 }
773 
775 {
776 
777 
778  if ( fMDTAU == 100 || fMDTAU == 200 )
779  {
780  int mode = selectLeptonic();
781  jaki_.jak1 = mode;
782  Tauola::setSameParticleDecayMode( mode );
783  mode = selectLeptonic();
784  jaki_.jak2 = mode;
785  Tauola::setOppositeParticleDecayMode( mode );
786  return ;
787  }
788 
789  int modeL = selectLeptonic();
790  int modeH = selectHadronic();
791 
792  if ( fMDTAU == 110 || fMDTAU == 210 )
793  {
794  jaki_.jak1 = modeL;
795  jaki_.jak2 = 0;
796  Tauola::setSameParticleDecayMode( modeL );
797  Tauola::setOppositeParticleDecayMode( 0 );
798  return ;
799  }
800 
801  if ( fMDTAU == 120 || fMDTAU == 22 )
802  {
803  jaki_.jak1 = 0;
804  jaki_.jak2 = modeL;
805  Tauola::setSameParticleDecayMode( 0 );
806  Tauola::setOppositeParticleDecayMode( modeL );
807  return;
808  }
809 
810  if ( fMDTAU == 114 || fMDTAU == 214 )
811  {
812  jaki_.jak1 = modeL;
813  jaki_.jak2 = modeH;
814  Tauola::setSameParticleDecayMode( modeL );
815  Tauola::setOppositeParticleDecayMode( modeH );
816  return;
817  }
818 
819  if ( fMDTAU == 124 || fMDTAU == 224 )
820  {
821  jaki_.jak1 = modeH;
822  jaki_.jak2 = modeL;
823  Tauola::setSameParticleDecayMode( modeH );
824  Tauola::setOppositeParticleDecayMode( modeL );
825  return;
826  }
827 
828  if ( fMDTAU == 115 || fMDTAU == 215 )
829  {
830  jaki_.jak1 = 1;
831  jaki_.jak2 = modeH;
832  Tauola::setSameParticleDecayMode( 1 );
833  Tauola::setOppositeParticleDecayMode( modeH );
834  return;
835  }
836 
837  if ( fMDTAU == 125 || fMDTAU == 225 )
838  {
839  jaki_.jak1 = modeH;
840  jaki_.jak2 = 1;
841  Tauola::setSameParticleDecayMode( modeH );
842  Tauola::setOppositeParticleDecayMode( 1 );
843  return;
844  }
845 
846  if ( fMDTAU == 116 || fMDTAU == 216 )
847  {
848  jaki_.jak1 = 2;
849  jaki_.jak2 = modeH;
850  Tauola::setSameParticleDecayMode( 2 );
851  Tauola::setOppositeParticleDecayMode( modeH );
852  return;
853  }
854 
855  if ( fMDTAU == 126 || fMDTAU == 226 )
856  {
857  jaki_.jak1 = modeH;
858  jaki_.jak2 = 2;
859  Tauola::setSameParticleDecayMode( modeH );
860  Tauola::setOppositeParticleDecayMode( 2 );
861  return;
862  }
863 
864  if ( fMDTAU == 130 || fMDTAU == 230 )
865  {
866  jaki_.jak1 = modeH;
867  jaki_.jak2 = selectHadronic();
868  Tauola::setSameParticleDecayMode( modeH );
869  Tauola::setOppositeParticleDecayMode( jaki_.jak2 );
870  return;
871  }
872 
873  if ( fMDTAU == 131 || fMDTAU == 231 )
874  {
875  jaki_.jak1 = modeH;
876  jaki_.jak2 = 0;
877  Tauola::setSameParticleDecayMode( modeH );
878  Tauola::setOppositeParticleDecayMode( 0 );
879  return;
880  }
881 
882  if ( fMDTAU == 132 || fMDTAU == 232 )
883  {
884  jaki_.jak1 = 0;
885  jaki_.jak2 = modeH;
886  Tauola::setSameParticleDecayMode( 0 );
887  Tauola::setOppositeParticleDecayMode( modeH );
888  return;
889  }
890 
891  // unlikely that we get here on unknown mdtau
892  // - there's a protection earlier
893  // but if we do, just set defaults
894  // probably need to spit a warning...
895  //
896  Tauola::setSameParticleDecayMode( 0 );
897  Tauola::setOppositeParticleDecayMode( 0 );
898 
899  return;
900 
901 
902 }
903 
905 {
906 
907  float prob = flat();
908 
909  if ( prob > 0. && prob <= fScaledLeptonBrRatios[0] )
910  {
911  return 1;
912  }
913  else if ( prob > fScaledLeptonBrRatios[1] && prob <=1. )
914  {
915  return 2;
916  }
917 
918  return 0;
919 }
920 
922 {
923 
924  float prob = 0.;
925  int len = 1;
926  ranmar_(&prob,&len);
927 
928  double sumBra = fScaledHadronBrRatios[0];
929  if ( prob > 0. && prob <= sumBra )
930  {
931  return fHadronModes[0];
932  }
933  else
934  {
935  int NN = fScaledHadronBrRatios.size();
936  for ( int i=1; i<NN; i++ )
937  {
938  if ( prob > sumBra && prob <= (sumBra+fScaledHadronBrRatios[i]) )
939  {
940  return fHadronModes[i];
941  }
942  sumBra += fScaledHadronBrRatios[i];
943  }
944  }
945 
946  return 0;
947 
948 }
949 
952  return (double)instance->flat();
953 }
954 
955 /* */
T getParameter(std::string const &) const
int i
Definition: DBlmapReader.cc:9
HepMC::GenEvent * decay(HepMC::GenEvent *)
static PFTauRenderPlugin instance
void setPSet(const edm::ParameterSet &)
#define abs(x)
Definition: mlp_lapack.h:159
void rmarin_(int *, int *, int *)
std::vector< double > fScaledLeptonBrRatios
void getData(T &iHolder) const
Definition: EventSetup.h:67
std::vector< int > fLeptonModes
static TauolaInterface * fInstance
double TauolappInterface_RandGetter()
bool isAvailable() const
Definition: Service.h:47
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
virtual CLHEP::HepRandomEngine & getEngine() const =0
Use this to get the random number engine, this is the only function most users should call...
void ranmar_(float *rvec, int *lenv)
HepPDT::ParticleData ParticleData
CLHEP::HepRandomEngine * fRandomEngine
edm::ParameterSet * fPSet
std::vector< int > fPDGs
static TauolaInterface * getInstance()
std::vector< double > fScaledHadronBrRatios
double p1[4]
Definition: TauolaWrapper.h:89
void init(const edm::EventSetup &)
std::vector< int > fHadronModes
struct @351 taubra_
int mdtau
Definition: TauolaWrapper.h:59