Skip to content

Commit

Permalink
---
Browse files Browse the repository at this point in the history 
yaml
---
r: 62302
b: refs/heads/l1tmuon-upgrade-dev
c: 7f3a270
h: refs/heads/l1tmuon-upgrade-dev
  • Loading branch information
Julia Yarba committed Mar 18, 2009
1 parent 5b668ff commit 5c70168
Show file tree
Hide file tree
Showing 4 changed files with 127 additions and 9 deletions.
2 changes: 1 addition & 1 deletion [refs]
View file
Original file line number Diff line number Diff line change
@@ -1,2 +1,2 @@
---
refs/heads/l1tmuon-upgrade-dev: 5391ca82dce0423f6ed012af12d45020d1d71263
refs/heads/l1tmuon-upgrade-dev: 7f3a270d87d504697e3461ecb1701d7e6fdea7f2
13 changes: 12 additions & 1 deletion trunk/GeneratorInterface/ExternalDecays/interface/TauolaInterface.h
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,14 @@
#ifndef gen_ExternalDecays_TauolaInterface_h
#define gen_ExternalDecays_TauolaInterface_h

// #include "HepPDT/defs.h"
// #include "HepPDT/TableBuilder.hh"
#include "HepPDT/ParticleDataTable.hh"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/Framework/interface/EventSetup.h"


namespace HepMC
{
Expand All @@ -19,7 +26,7 @@ namespace gen {
TauolaInterface( const edm::ParameterSet& );
~TauolaInterface();

void init();
void init( const edm::EventSetup& );
const std::vector<int>& operatesOnParticles() { return fPDGs; }
HepMC::GenEvent* decay( const HepMC::GenEvent* );
void statistics() ;
Expand All @@ -29,6 +36,10 @@ namespace gen {
//
std::vector<int> fPDGs;
int fPolarization;

edm::ESHandle<HepPDT::ParticleDataTable> fPDGTable ;
//CLHEP::HepRandomEngine* fRandomEngine ;
//CLHEP::RandFlat* fRandomGenerator;

};

Expand Down
5 changes: 4 additions & 1 deletion trunk/GeneratorInterface/ExternalDecays/interface/TauolaWrapper.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,10 +26,13 @@
//-------------------------------------------------------------------------------
//

// main function of TAUOLA tau decay library
// main function(s) of TAUOLA/pretauola tau decay library

extern "C" {
void tauola_(int*, int*);
void tauola_srs_(int*,int*);
void taurep_(int*);
void ranmar_(float*,int*);
}
#define tauola tauola_

Expand Down
116 changes: 110 additions & 6 deletions trunk/GeneratorInterface/ExternalDecays/src/TauolaInterface.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@
#include "HepMC/IO_HEPEVT.h"
#include "HepMC/HEPEVT_Wrapper.h"

#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"

using namespace gen;
using namespace edm;
Expand Down Expand Up @@ -37,14 +38,20 @@ TauolaInterface::~TauolaInterface()
{
}

void TauolaInterface::init()
void TauolaInterface::init( const edm::EventSetup& es )
{

if ( ki_taumod_.mdtau <= -1 ) // actually, need to throw exception !
return ;

fPDGs.push_back( 15 ) ;

es.getData( fPDGTable ) ;
/*
Service<RandomNumberGenerator> rng;
long seed = (long)(rng->mySeed()) ;
fRandomEngine = new HepJamesRandom(seed) ;
fRandomGenerator = new RandFlat(fRandomEngine) ;
*/
cout << "----------------------------------------------" << endl;
cout << "Initializing Tauola" << endl;
if ( fPolarization == 0 )
Expand All @@ -55,8 +62,11 @@ void TauolaInterface::init()
{
cout << "Tauola: Polarization enabled" << endl;
}
int mode = -1;
tauola_( &mode, &fPolarization );
int mode = -2;
taurep_( &mode ) ;
mode = -1;
// tauola_( &mode, &fPolarization );
tauola_srs_( &mode, &fPolarization );

return;
}
Expand All @@ -82,18 +92,112 @@ HepMC::GenEvent* TauolaInterface::decay( const HepMC::GenEvent* evt )
// HepMC::HEPEVT_Wrapper::print_hepevt();

int mode = 0;
tauola_( &mode, &fPolarization );
// tauola_( &mode, &fPolarization );
tauola_srs_( &mode, &fPolarization );

int numPartAfterTauola = HepMC::HEPEVT_Wrapper::number_entries();
// HepMC::HEPEVT_Wrapper::print_hepevt();

HepMC::IO_HEPEVT conv;

// before we do the conversion, we need to deal with decay vertexes
// since Tauola knows nothing about lifetimes, all decay vertexes are set to 0.
// nees to set them properly, knowing lifetime !
// here we do it on HEPEVT record, also for consistency, although it's probably
// even easier to deal with HepMC::GenEvent record

// find 1st "non-doc" tau
//
bool foundTau = false;
for ( int ip=1; ip<=numPartAfterTauola; ip++ )
{
if ( std::abs( HepMC::HEPEVT_Wrapper::id( ip ) ) == 15
&& HepMC::HEPEVT_Wrapper::status( ip ) != 3 )
{
foundTau = true;
break;
}
}

if ( !foundTau )
{
// no tau found
// just give up here
//
return conv.read_next_event();
}

std::vector<int> PrntIndx;

for ( int ip=numPartAfterTauola; ip>numPartBeforeTauola; ip-- ) // Fortran indexing !
{

// first of all, find out how many generations in decay chain
//
PrntIndx.clear();
int Prnt = HepMC::HEPEVT_Wrapper::first_parent(ip);
ip -= (HepMC::HEPEVT_Wrapper::number_children(Prnt)-1); // such that we don't go the same part again
PrntIndx.push_back( Prnt );
while ( abs( HepMC::HEPEVT_Wrapper::id(Prnt) ) != 15 ) // shortcut; need to loop over fPDGs...
{
int Prnt1 = HepMC::HEPEVT_Wrapper::first_parent( Prnt );
Prnt = Prnt1;
// such that the tau always appear at the start of the list
PrntIndx.insert( PrntIndx.begin(), Prnt );
ip -= HepMC::HEPEVT_Wrapper::number_children(Prnt); // such that we don't go the same part again
}
for ( int iprt=0; iprt<PrntIndx.size(); iprt++ )
{
int Indx = PrntIndx[iprt];
int PartID = HepMC::HEPEVT_Wrapper::id( Indx );
const HepPDT::ParticleData*
PData = fPDGTable->particle(HepPDT::ParticleID(abs(PartID))) ;
//
// prob = exp(-t/lifetime) ==> t = -lifetime * log(prob)
//
float prob = 0.;
int length=1;
ranmar_(&prob,&length);
double lifetime = PData->lifetime().value();
//
// in case of Py6, this would be copied into V(5,i)
// for HEPEVT, need to check...
//
double ct = -lifetime * std::log(prob);
//
double ee = HepMC::HEPEVT_Wrapper::e( Indx );
double px = HepMC::HEPEVT_Wrapper::px( Indx );
double py = HepMC::HEPEVT_Wrapper::py( Indx );
double pz = HepMC::HEPEVT_Wrapper::pz( Indx );
double pp = std::sqrt( px*px + py*py + pz*pz );
double mass = HepMC::HEPEVT_Wrapper::m( Indx );
//
// this is in py6 terms:
// VDCY(J)=V(IP,J)+V(IP,5)*P(IP,J)/P(IP,5)
//
double VxDec = HepMC::HEPEVT_Wrapper::x( Indx );
VxDec += ct * (px/mass);
double VyDec = HepMC::HEPEVT_Wrapper::y( Indx );
VyDec += ct * (py/mass);
double VzDec = HepMC::HEPEVT_Wrapper::z( Indx );
VzDec += ct * (pz/mass);
double VtDec = HepMC::HEPEVT_Wrapper::t( Indx );
VtDec += ct * (ee/mass);
for ( int idau=HepMC::HEPEVT_Wrapper::first_child( Indx );
idau<=HepMC::HEPEVT_Wrapper::last_child( Indx ); idau++ )
{
HepMC::HEPEVT_Wrapper::set_position( idau, VxDec, VyDec, VzDec, VtDec );
}
}
}

return conv.read_next_event();

}

void TauolaInterface::statistics()
{
int mode = 1;
tauola_( &mode, &fPolarization ) ;
// tauola_( &mode, &fPolarization ) ;
tauola_srs_( &mode, &fPolarization ) ;
}

0 comments on commit 5c70168

Please sign in to comment.