17 #include "CLHEP/Random/Random.h"
22 #include "HepMC/GenEvent.h"
40 std::cout <<
" EvtGenProducer starting ... " << std::endl;
50 <<
"The EvtGenProducer module requires the RandomNumberGeneratorService\n"
51 "which is not present in the configuration file. You must add the service\n"
52 "in the configuration file if you want to run EvtGenProducer";
55 CLHEP::HepRandomEngine& m_engine = rngen->
getEngine();
56 m_flat =
new CLHEP::RandFlat(m_engine, 0., 1.);
67 std::vector<double>());
68 if (polarize_ids.size() != polarize_pol.size()) {
70 <<
"EvtGenProducer requires that the particles_to_polarize and particle_polarization\n"
71 "vectors be the same size. Please fix this in your configuration.";
73 for (
unsigned int ndx = 0; ndx < polarize_ids.size(); ndx++) {
74 if (polarize_pol[ndx] < -1. || polarize_pol[ndx] > 1.) {
76 <<
"EvtGenProducer error: particle polarizations must be in the range -1 < P < 1";
78 polarizations.insert(std::pair<int, float>(polarize_ids[ndx], polarize_pol[ndx]));
90 std::vector<std::string> forced_names = pset.
getParameter< std::vector<std::string> >(
"list_forced_decays");
92 m_EvtGen =
new EvtGen (decay_table_s.c_str(),pdt_s.c_str(),the_engine);
95 if (!useDefault) m_EvtGen->readUDecay( user_decay_s.c_str() );
97 std::vector<std::string>::const_iterator
i;
100 for (i=forced_names.begin(); i!=forced_names.end(); ++
i)
106 if (found.getId()==-1)
109 <<
"name in part list for forced decays not found: " << *
i;
111 if (found.getId()==found.getAlias())
114 <<
"name in part list for forced decays is not an alias: " << *
i;
116 forced_Evt.push_back(found);
117 forced_Hep.push_back(EvtPDL::getStdHep(found));
126 m_PDGs.push_back( 300553 ) ;
127 m_PDGs.push_back( 511 ) ;
128 m_PDGs.push_back( 521 ) ;
129 m_PDGs.push_back( 523 ) ;
130 m_PDGs.push_back( 513 ) ;
131 m_PDGs.push_back( 533 ) ;
132 m_PDGs.push_back( 531 ) ;
134 m_PDGs.push_back( 15 ) ;
136 m_PDGs.push_back( 413 ) ;
137 m_PDGs.push_back( 423 ) ;
138 m_PDGs.push_back( 433 ) ;
139 m_PDGs.push_back( 411 ) ;
140 m_PDGs.push_back( 421 ) ;
141 m_PDGs.push_back( 431 ) ;
142 m_PDGs.push_back( 10411 );
143 m_PDGs.push_back( 10421 );
144 m_PDGs.push_back( 10413 );
145 m_PDGs.push_back( 10423 );
146 m_PDGs.push_back( 20413 );
147 m_PDGs.push_back( 20423 );
149 m_PDGs.push_back( 415 );
150 m_PDGs.push_back( 425 );
151 m_PDGs.push_back( 10431 );
152 m_PDGs.push_back( 20433 );
153 m_PDGs.push_back( 10433 );
154 m_PDGs.push_back( 435 );
156 m_PDGs.push_back( 310 );
157 m_PDGs.push_back( 311 );
158 m_PDGs.push_back( 313 );
159 m_PDGs.push_back( 323 );
160 m_PDGs.push_back( 10321 );
161 m_PDGs.push_back( 10311 );
162 m_PDGs.push_back( 10313 );
163 m_PDGs.push_back( 10323 );
164 m_PDGs.push_back( 20323 );
165 m_PDGs.push_back( 20313 );
166 m_PDGs.push_back( 325 );
167 m_PDGs.push_back( 315 );
169 m_PDGs.push_back( 100313 );
170 m_PDGs.push_back( 100323 );
171 m_PDGs.push_back( 30313 );
172 m_PDGs.push_back( 30323 );
173 m_PDGs.push_back( 30343 );
174 m_PDGs.push_back( 30353 );
175 m_PDGs.push_back( 30363 );
177 m_PDGs.push_back( 111 );
178 m_PDGs.push_back( 221 );
179 m_PDGs.push_back( 113 );
180 m_PDGs.push_back( 213 );
181 m_PDGs.push_back( 223 );
182 m_PDGs.push_back( 331 );
183 m_PDGs.push_back( 333 );
184 m_PDGs.push_back( 20213 );
185 m_PDGs.push_back( 20113 );
186 m_PDGs.push_back( 215 );
187 m_PDGs.push_back( 115 );
188 m_PDGs.push_back( 10213 );
189 m_PDGs.push_back( 10113 );
190 m_PDGs.push_back( 9000111 );
191 m_PDGs.push_back( 9000211 );
192 m_PDGs.push_back( 9010221 );
193 m_PDGs.push_back( 10221 );
194 m_PDGs.push_back( 20223 );
195 m_PDGs.push_back( 20333 );
196 m_PDGs.push_back( 225 );
197 m_PDGs.push_back( 9020221 );
198 m_PDGs.push_back( 335 );
199 m_PDGs.push_back( 10223 );
200 m_PDGs.push_back( 10333 );
201 m_PDGs.push_back( 100213 );
202 m_PDGs.push_back( 100113 );
204 m_PDGs.push_back( 441 );
205 m_PDGs.push_back( 100441 );
206 m_PDGs.push_back( 443 );
207 m_PDGs.push_back( 100443 );
208 m_PDGs.push_back( 9000443 );
209 m_PDGs.push_back( 9010443 );
210 m_PDGs.push_back( 9020443 );
211 m_PDGs.push_back( 10441 );
212 m_PDGs.push_back( 20443 );
213 m_PDGs.push_back( 445 );
215 m_PDGs.push_back( 30443 );
216 m_PDGs.push_back( 551 );
217 m_PDGs.push_back( 553 );
218 m_PDGs.push_back( 100553 );
219 m_PDGs.push_back( 200553 );
220 m_PDGs.push_back( 10551 );
221 m_PDGs.push_back( 20553 );
222 m_PDGs.push_back( 555 );
223 m_PDGs.push_back( 10553 );
225 m_PDGs.push_back( 110551 );
226 m_PDGs.push_back( 120553 );
227 m_PDGs.push_back( 100555 );
228 m_PDGs.push_back( 210551 );
229 m_PDGs.push_back( 220553 );
230 m_PDGs.push_back( 200555 );
231 m_PDGs.push_back( 30553 );
232 m_PDGs.push_back( 20555 );
234 m_PDGs.push_back( 557 );
235 m_PDGs.push_back( 130553 );
236 m_PDGs.push_back( 120555 );
237 m_PDGs.push_back( 100557 );
238 m_PDGs.push_back( 110553 );
239 m_PDGs.push_back( 210553 );
240 m_PDGs.push_back( 10555 );
241 m_PDGs.push_back( 110555 );
243 m_PDGs.push_back( 4122 );
244 m_PDGs.push_back( 4132 );
246 m_PDGs.push_back( 4112 );
247 m_PDGs.push_back( 4212 );
248 m_PDGs.push_back( 4232 );
249 m_PDGs.push_back( 4222 );
250 m_PDGs.push_back( 4322 );
251 m_PDGs.push_back( 4312 );
253 m_PDGs.push_back( 13122 );
254 m_PDGs.push_back( 13124 );
255 m_PDGs.push_back( 23122 );
256 m_PDGs.push_back( 33122 );
257 m_PDGs.push_back( 43122 );
258 m_PDGs.push_back( 53122 );
259 m_PDGs.push_back( 13126 );
260 m_PDGs.push_back( 13212 );
261 m_PDGs.push_back( 13241 );
263 m_PDGs.push_back( 3126 );
264 m_PDGs.push_back( 3124 );
265 m_PDGs.push_back( 3122 );
266 m_PDGs.push_back( 3222 );
267 m_PDGs.push_back( 2214 );
268 m_PDGs.push_back( 2224 );
269 m_PDGs.push_back( 3324 );
270 m_PDGs.push_back( 2114 );
271 m_PDGs.push_back( 1114 );
272 m_PDGs.push_back( 3112 );
273 m_PDGs.push_back( 3212 );
274 m_PDGs.push_back( 3114 );
275 m_PDGs.push_back( 3224 );
276 m_PDGs.push_back( 3214 );
277 m_PDGs.push_back( 3216 );
278 m_PDGs.push_back( 3322 );
279 m_PDGs.push_back( 3312 );
280 m_PDGs.push_back( 3314 );
281 m_PDGs.push_back( 3334 );
283 m_PDGs.push_back( 4114 );
284 m_PDGs.push_back( 4214 );
285 m_PDGs.push_back( 4224 );
286 m_PDGs.push_back( 4314 );
287 m_PDGs.push_back( 4324 );
288 m_PDGs.push_back( 4332 );
289 m_PDGs.push_back( 4334 );
292 m_PDGs.push_back( 10443 );
294 m_PDGs.push_back( 5122 );
295 m_PDGs.push_back( 5132 );
296 m_PDGs.push_back( 5232 );
297 m_PDGs.push_back( 5332 );
298 m_PDGs.push_back( 5222 );
299 m_PDGs.push_back( 5112 );
300 m_PDGs.push_back( 5212 );
301 m_PDGs.push_back( 541 );
302 m_PDGs.push_back( 14122 );
303 m_PDGs.push_back( 14124 );
304 m_PDGs.push_back( 5312 );
305 m_PDGs.push_back( 5322 );
306 m_PDGs.push_back( 10521 );
307 m_PDGs.push_back( 20523 );
308 m_PDGs.push_back( 10523 );
310 m_PDGs.push_back( 525 );
311 m_PDGs.push_back( 10511 );
312 m_PDGs.push_back( 20513 );
313 m_PDGs.push_back( 10513 );
314 m_PDGs.push_back( 515 );
315 m_PDGs.push_back( 10531 );
316 m_PDGs.push_back( 20533 );
317 m_PDGs.push_back( 10533 );
318 m_PDGs.push_back( 535 );
319 m_PDGs.push_back( 543 );
320 m_PDGs.push_back( 545 );
321 m_PDGs.push_back( 5114 );
322 m_PDGs.push_back( 5224 );
323 m_PDGs.push_back( 5214 );
324 m_PDGs.push_back( 5314 );
325 m_PDGs.push_back( 5324 );
326 m_PDGs.push_back( 5334 );
327 m_PDGs.push_back( 10541 );
328 m_PDGs.push_back( 10543 );
329 m_PDGs.push_back( 20543 );
331 m_PDGs.push_back( 4424 );
332 m_PDGs.push_back( 4422 );
333 m_PDGs.push_back( 4414 );
334 m_PDGs.push_back( 4412 );
335 m_PDGs.push_back( 4432 );
336 m_PDGs.push_back( 4434 );
338 m_PDGs.push_back( 130 );
341 if ( pset.
exists(
"operates_on_particles") )
343 std::vector<int> tmpPIDs = pset.
getParameter< std::vector<int> >(
"operates_on_particles");
344 if ( tmpPIDs.size() > 0 )
346 if ( tmpPIDs[0] > 0 )
359 std::cout <<
" EvtGenProducer terminating ... " << std::endl;
369 if (usePythia) EvtPythia::pythiaInit(0);
393 nPythia = evt->particles_size();
402 for (HepMC::GenEvent::particle_const_iterator
p= evt->particles_begin();
p != evt->particles_end(); ++
p)
404 status = (*p)->status();
409 idHep = (*p)->pdg_id();
411 for(
int i=0;
i<nforced;
i++)
413 if(idHep == forced_Hep[
i])
415 update_candlist(i,*
p);
421 idEvt = EvtPDL::evtIdFromStdHep(idHep);
422 ipart = idEvt.getId();
423 if (ipart==-1)
continue;
424 if (EvtDecayTable::getNMode(ipart)==0)
continue;
425 addToHepMC(*
p,idEvt,evt,
true);
433 int which = (int)(nlist*m_flat->fire());
434 if (which == nlist) which = nlist-1;
436 for(
int k=0;
k < nlist;
k++)
438 if(
k == which || !usePythia) {
439 addToHepMC(listp[
k],forced_Evt[
index[k]],evt,
false);
443 int id_non_alias = forced_Evt[
index[
k]].getId();
444 EvtId non_alias(id_non_alias,id_non_alias);
445 addToHepMC(listp[
k],non_alias,evt,
false);
457 EvtSpinType::spintype stype = EvtPDL::getSpinType(idEvt);
458 EvtParticle* partEvt;
460 case EvtSpinType::SCALAR:
461 partEvt =
new EvtScalarParticle();
464 partEvt =
new EvtStringParticle();
466 case EvtSpinType::DIRAC:
467 partEvt =
new EvtDiracParticle();
469 case EvtSpinType::VECTOR:
470 partEvt =
new EvtVectorParticle();
472 case EvtSpinType::RARITASCHWINGER:
473 partEvt =
new EvtRaritaSchwingerParticle();
475 case EvtSpinType::TENSOR:
476 partEvt =
new EvtTensorParticle();
478 case EvtSpinType::SPIN5HALF:
case EvtSpinType::SPIN3:
case EvtSpinType::SPIN7HALF:
case EvtSpinType::SPIN4:
479 partEvt =
new EvtHighSpinParticle();
482 std::cout <<
"Unknown spintype in EvtSpinType!" << std::endl;
488 HepMC::FourVector momHep = partHep->momentum();
489 momEvt.set(momHep.t(),momHep.x(),momHep.y(),momHep.z());
491 HepMC::GenVertex* initVert = partHep->production_vertex();
492 HepMC::FourVector posHep = initVert->position();
493 posEvt.set(posHep.t(),posHep.x(),posHep.y(),posHep.z());
494 partEvt->init(idEvt,momEvt);
495 if (stype == EvtSpinType::DIRAC
496 && polarizations.find(partHep->pdg_id()) != polarizations.end()) {
501 float pol = polarizations.find(partHep->pdg_id())->
second;
510 EvtSpinDensity theSpinDensity;
511 theSpinDensity.SetDim(2);
512 theSpinDensity.Set(0, 0, EvtComplex(1./2. + polVec.
z()/2., 0.));
513 theSpinDensity.Set(0, 1, EvtComplex(polVec.
x()/2., -polVec.
y()/2.));
514 theSpinDensity.Set(1, 0, EvtComplex(polVec.
x()/2., polVec.
y()/2.));
515 theSpinDensity.Set(1, 1, EvtComplex(1./2. - polVec.
z()/2., 0.));
517 partEvt->setSpinDensityForwardHelicityBasis(theSpinDensity);
520 partEvt->setDiagonalSpinDensity();
526 if (del_daug) go_through_daughters(partEvt);
529 static EvtStdHep evtstdhep;
532 partEvt->makeStdHep(evtstdhep);
535 partEvt->deleteTree();
540 HepMC::GenVertex* theVerts[200];
541 for (
int ivert = 0; ivert < 200; ivert++) {
545 for (
int ipart = 0; ipart < evtstdhep.getNPart(); ipart++) {
546 int theMum = evtstdhep.getFirstMother(ipart);
547 if (theMum != -1 && !theVerts[theMum]) {
548 EvtVector4R theVpos = evtstdhep.getX4(ipart) + posEvt;
550 new HepMC::GenVertex(HepMC::FourVector(theVpos.get(1),
558 partHep->set_status(2);
559 if (theVerts[0]) theVerts[0]->add_particle_in( partHep );
561 for (
int ipart2 = 1; ipart2 < evtstdhep.getNPart(); ipart2++) {
562 int idHep = evtstdhep.getStdHepID(ipart2);
565 evtstdhep.getP4(ipart2).get(2),
566 evtstdhep.getP4(ipart2).get(3),
567 evtstdhep.getP4(ipart2).get(0)),
569 evtstdhep.getIStat(ipart2));
571 thePart->suggest_barcode(npartial + nPythia);
572 int theMum2 = evtstdhep.getFirstMother(ipart2);
573 if (theMum2 != -1 && theVerts[theMum2]) theVerts[theMum2]->add_particle_out( thePart );
574 if (theVerts[ipart2]) theVerts[ipart2]->add_particle_in( thePart );
578 for (
int ipart3 = 0; ipart3 < evtstdhep.getNPart(); ipart3++) {
579 if (theVerts[ipart3]) theEvent->add_vertex( theVerts[ipart3] );
597 int NDaug=part->getNDaug();
600 EvtParticle* Daughter;
601 for (
int i=0;
i<NDaug;
i++)
603 Daughter=part->getDaug(
i);
604 int idHep = EvtPDL::getStdHep(Daughter->getId());
606 for(
int k=0;
k<nforced;
k++)
608 if(idHep == forced_Hep[
k])
611 Daughter->deleteDaughters();
614 if (!found) go_through_daughters(Daughter);
624 bool isThere =
false;
627 if (listp[
check]->barcode() == thePart->barcode()) isThere =
true;
631 listp[nlist] = thePart;
632 index[nlist++] = theIndex;
638 <<
"more than 10 candidates to be forced ";
T getParameter(std::string const &) const
T getUntrackedParameter(std::string const &, T const &) const
HepMC::GenEvent * decay(HepMC::GenEvent *)
void addToHepMC(HepMC::GenParticle *partHep, EvtId idEvt, HepMC::GenEvent *theEvent, bool del_daug)
void update_candlist(int theIndex, HepMC::GenParticle *thePart)
static unsigned int getId(void)
bool exists(std::string const ¶meterName) const
checks if a parameter exists
void go_through_daughters(EvtParticle *part)
U second(std::pair< T, U > const &p)
bool check(const DataFrame &df, bool capcheck, bool dvercheck)
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
virtual CLHEP::HepRandomEngine & getEngine() const =0
Use this to get the random number engine, this is the only function most users should call...
Vector3DBase unit() const
EvtGenInterface(const edm::ParameterSet &)
std::string fullPath() const