CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
Pythia6Gun.cc
Go to the documentation of this file.
1 /*
2  * $Date: 2013/03/01 22:23:10 $
3  * $Revision: 1.20 $
4  * \author Julia Yarba
5  */
6 
7 #include <iostream>
8 
9 #include "Pythia6Gun.h"
10 
11 //#include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
14 
18 
20 
21 using namespace edm;
22 using namespace gen;
23 
24 
25 Pythia6Gun::Pythia6Gun( const ParameterSet& pset ) :
26  fPy6Service( new Pythia6Service(pset) ),
27  fEvt(0)
28  // fPDGTable( new DefaultConfig::ParticleDataTable("PDG Table") )
29 {
30 
31  ParameterSet pgun_params =
32  pset.getParameter<ParameterSet>("PGunParameters");
33 
34  // although there's the method ParameterSet::empty(),
35  // it looks like it's NOT even necessary to check if it is,
36  // before trying to extract parameters - if it is empty,
37  // the default values seem to be taken
38  //
39  fMinPhi = pgun_params.getParameter<double>("MinPhi"); // ,-3.14159265358979323846);
40  fMaxPhi = pgun_params.getParameter<double>("MaxPhi"); // , 3.14159265358979323846);
41 
42  fHepMCVerbosity = pset.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false ) ;
43  fPylistVerbosity = pset.getUntrackedParameter<int>( "pythiaPylistVerbosity", 0 ) ;
44  fMaxEventsToPrint = pset.getUntrackedParameter<int>( "maxEventsToPrint", 0 );
45 
46 // Turn off banner printout
47  if (!call_pygive("MSTU(12)=12345"))
48  {
49  throw edm::Exception(edm::errors::Configuration,"PythiaError")
50  <<" pythia did not accept MSTU(12)=12345";
51  }
52 
53  produces<HepMCProduct>();
54 
55 }
56 
58 {
59  if ( fPy6Service ) delete fPy6Service;
60  //
61  // note that GenEvent or any undelaying (GenVertex, GenParticle) do NOT
62  // need to be cleaned, as it'll be done automatically by HepMCProduct
63  //
64 }
65 
66 
68 {
69  // es.getData( fPDGTable ) ;
70  return ;
71 
72 }
73 
74 void Pythia6Gun::beginRun( Run const&, EventSetup const& es )
75 {
76  assert ( fPy6Service ) ;
77 
78  Pythia6Service::InstanceWrapper guard(fPy6Service); // grab Py6 instance
79 
83 
84  call_pygive("MSTU(10)=1");
85 
86  call_pyinit("NONE", "", "", 0.0);
87 
88  std::cout << " FYI: MSTU(10)=1 is ENFORCED in Py6-PGuns, for technical reasons"
89  << std::endl;
90  return;
91 }
92 
93 void Pythia6Gun::endRun( Run const&, EventSetup const& es )
94 {
95 
96  // here put in GenRunInfoProduct
97 
98  call_pystat(1);
99 
100  return;
101 }
102 
104 {
105 
106  for ( int iprt=fPartIDs.size(); iprt<pyjets.n; iprt++ ) // the pointer is shifted by -1, c++ style
107  {
108  int parent = pyjets.k[2][iprt];
109  if ( parent != 0 )
110  {
111  // pull up parent particle
112  //
113  HepMC::GenParticle* parentPart = fEvt->barcode_to_particle( parent );
114  parentPart->set_status( 2 ); // reset status, to mark that it's decayed
115 
116  HepMC::GenVertex* DecVtx = new HepMC::GenVertex(HepMC::FourVector(pyjets.v[0][iprt],
117  pyjets.v[1][iprt],
118  pyjets.v[2][iprt],
119  pyjets.v[3][iprt]));
120  DecVtx->add_particle_in( parentPart ); // this will cleanup end_vertex if exists,
121  // and replace with the new one
122  // I presume barcode will be given automatically
123 
124  HepMC::FourVector pmom(pyjets.p[0][iprt],pyjets.p[1][iprt],
125  pyjets.p[2][iprt],pyjets.p[3][iprt] );
126 
127  int dstatus = 0;
128  if ( pyjets.k[0][iprt] >= 1 && pyjets.k[0][iprt] <= 10 )
129  {
130  dstatus = 1;
131  }
132  else if ( pyjets.k[0][iprt] >= 11 && pyjets.k[0][iprt] <= 20 )
133  {
134  dstatus = 2;
135  }
136  else if ( pyjets.k[0][iprt] >= 21 && pyjets.k[0][iprt] <= 30 )
137  {
138  dstatus = 3;
139  }
140  else if ( pyjets.k[0][iprt] >= 31 && pyjets.k[0][iprt] <= 100 )
141  {
142  dstatus = pyjets.k[0][iprt];
143  }
144  HepMC::GenParticle* daughter =
145  new HepMC::GenParticle(pmom,
146  HepPID::translatePythiatoPDT( pyjets.k[1][iprt] ),
147  dstatus);
148  daughter->suggest_barcode( iprt+1 );
149  DecVtx->add_particle_out( daughter );
150  // give particle barcode as well !
151 
152  int iprt1;
153  for ( iprt1=iprt+1; iprt1<pyjets.n; iprt1++ ) // the pointer is shifted by -1, c++ style
154  {
155  if ( pyjets.k[2][iprt1] != parent ) break; // another parent particle, break the loop
156 
157  HepMC::FourVector pmomN(pyjets.p[0][iprt1],pyjets.p[1][iprt1],
158  pyjets.p[2][iprt1],pyjets.p[3][iprt1] );
159 
160  dstatus = 0;
161  if ( pyjets.k[0][iprt1] >= 1 && pyjets.k[0][iprt1] <= 10 )
162  {
163  dstatus = 1;
164  }
165  else if ( pyjets.k[0][iprt1] >= 11 && pyjets.k[0][iprt1] <= 20 )
166  {
167  dstatus = 2;
168  }
169  else if ( pyjets.k[0][iprt1] >= 21 && pyjets.k[0][iprt1] <= 30 )
170  {
171  dstatus = 3;
172  }
173  else if ( pyjets.k[0][iprt1] >= 31 && pyjets.k[0][iprt1] <= 100 )
174  {
175  dstatus = pyjets.k[0][iprt1];
176  }
177  HepMC::GenParticle* daughterN =
178  new HepMC::GenParticle(pmomN,
179  HepPID::translatePythiatoPDT( pyjets.k[1][iprt1] ),
180  dstatus);
181  daughterN->suggest_barcode( iprt1+1 );
182  DecVtx->add_particle_out( daughterN );
183  }
184 
185  iprt = iprt1-1; // reset counter such that it doesn't go over the same child more than once
186  // don't forget to offset back into c++ counting, as it's already +1 forward
187 
188  fEvt->add_vertex( DecVtx );
189 
190  }
191  }
192 
193  return;
194 
195 }
196 
198 {
199 
200  generateEvent() ;
201 
202  fEvt->set_beam_particles(0,0);
203  fEvt->set_event_number(evt.id().event()) ;
204  fEvt->set_signal_process_id(pypars.msti[0]) ;
205 
207 
208  int evtN = evt.id().event();
209  if ( evtN <= fMaxEventsToPrint )
210  {
211  if ( fPylistVerbosity )
212  {
214  }
215  if ( fHepMCVerbosity )
216  {
217  if ( fEvt ) fEvt->print();
218  }
219  }
220 
221  loadEvent( evt );
222 }
223 
225 {
226 
227  std::auto_ptr<HepMCProduct> bare_product(new HepMCProduct());
228 
229  if(fEvt) bare_product->addHepMCData( fEvt );
230 
231  evt.put(bare_product);
232 
233 
234  return;
235 
236 }
237 
239  double& ee, double& eta, double& phi )
240 {
241 
242  if ( ip < 2 ) return 0;
243 
244 // translate PDG to Py6
245  int py6PID = HepPID::translatePDTtoPythia( particleID );
246 // Check if particle is its own anti-particle.
247  int pythiaCode = pycomp_(py6PID); // this is py6 internal validity check, it takes Pythia6 pid
248  // so actually I'll need to convert
249  int has_antipart = pydat2.kchg[3-1][pythiaCode-1];
250  int particleID2 = has_antipart ? -1 * particleID : particleID; // this is PDG, for HepMC::GenEvent
251  int py6PID2 = has_antipart ? -1 * py6PID : py6PID; // this py6 id, for py1ent
252  double the = 2.*atan(exp(eta));
253  phi = phi + M_PI;
254  if (phi > 2.* M_PI) {phi = phi - 2.* M_PI;}
255 
256  // copy over mass of the previous one, because then py6 will pick it up
257  pyjets.p[4][ip-1] = pyjets.p[4][ip-2];
258 
259  py1ent_(ip, py6PID2, ee, the, phi);
260 
261  double px = pyjets.p[0][ip-1]; // pt*cos(phi) ;
262  double py = pyjets.p[1][ip-1]; // pt*sin(phi) ;
263  double pz = pyjets.p[2][ip-1]; // mom*cos(the) ;
264  HepMC::FourVector ap(px,py,pz,ee) ;
265  HepMC::GenParticle* APart =
266  new HepMC::GenParticle(ap,particleID2,1);
267  APart->suggest_barcode( ip ) ;
268 
269  return APart;
270 
271 }
T getParameter(std::string const &) const
EventNumber_t event() const
Definition: EventID.h:44
T getUntrackedParameter(std::string const &, T const &) const
void beginJob() override
Definition: Pythia6Gun.cc:67
void endRun(edm::Run const &, edm::EventSetup const &) override
Definition: Pythia6Gun.cc:93
list parent
Definition: dbtoconf.py:74
void produce(edm::Event &, const edm::EventSetup &) override
Definition: Pythia6Gun.cc:197
void beginRun(edm::Run const &, edm::EventSetup const &) override
Definition: Pythia6Gun.cc:74
bool call_pygive(const std::string &line)
HepMC::GenParticle * addAntiParticle(int &, int &, double &, double &, double &)
Definition: Pythia6Gun.cc:238
HepMC::GenEvent * fEvt
Definition: Pythia6Gun.h:61
bool fHepMCVerbosity
Definition: Pythia6Gun.h:68
int fPylistVerbosity
Definition: Pythia6Gun.h:69
T eta() const
virtual void generateEvent()=0
void call_pylist(int mode)
virtual ~Pythia6Gun()
Definition: Pythia6Gun.cc:57
void attachPy6DecaysToGenEvent()
Definition: Pythia6Gun.cc:103
OrphanHandle< PROD > put(std::auto_ptr< PROD > product)
Put a new product.
Definition: Event.h:94
double fMinPhi
Definition: Pythia6Gun.h:56
int fMaxEventsToPrint
Definition: Pythia6Gun.h:70
void loadEvent(edm::Event &)
Definition: Pythia6Gun.cc:224
int pycomp_(int &)
#define pypars
double fMaxPhi
Definition: Pythia6Gun.h:57
#define M_PI
Definition: BFit3D.cc:3
void py1ent_(int &ip, int &kf, double &pe, double &the, double &phi)
edm::EventID id() const
Definition: EventBase.h:56
Pythia6Service * fPy6Service
Definition: Pythia6Gun.h:51
tuple cout
Definition: gather_cfg.py:121
std::vector< int > fPartIDs
Definition: Pythia6Gun.h:55
Definition: Run.h:36
Definition: DDAxes.h:10