001
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017 #include "Pythia_i/Pythia.h"
018 #include "GeneratorModules/GeneratorName.h"
019
020
021 #include "GaudiKernel/MsgStream.h"
022
023
024 #include "StoreGate/StoreGateSvc.h"
025
026 #include "HepMC/IO_HEPEVT.h"
027 #include "HepMC/IO_Ascii.h"
028 #include "HepMC/HEPEVT_Wrapper.h"
029
030
031
032
033 #include "GeneratorUtils/StringParse.h"
034 #include <stdlib.h>
035
036
037 #include "CLHEP/Random/RandFlat.h"
038 #include "AthenaKernel/IAtRndmGenSvc.h"
039
040
041 IAtRndmGenSvc* Pythia::p_AtRndmGenSvc = 0;
042 std::string Pythia::pythia_stream = "PYTHIA_INIT";
043 extern "C" double atl_pyr_( int* )
044 {
045 CLHEP::HepRandomEngine* engine = Pythia::p_AtRndmGenSvc->GetEngine(Pythia::pythia_stream);
046 return CLHEP::RandFlat::shoot(engine);
047 }
048
049 extern "C" {
050 void initpyblock_(int*, const char*);
051 int openrandom_(int*, int*, const char*);
052 void pyinit_(const char*,
053 const char*,
054 const char*,
055 double* ,
056 int,
057 int,
058 int
059 );
060 void pyevnt_();
061 void pyevnw_();
062 void pystat_(int*);
063 void pylist_(int*);
064 void pyupda_(int*, int*);
065 void pyrget_(int*, int*);
066 void pyrset_(int*, int*);
067 int pycomp_(int*);
068
069 void pyhepc_(int*);
070 int opdcay_(const char*, int*, char*, int, int);
071 void cldcay_(int*);
072 void rinpar_();
073 void subnum_();
074 void rheader_();
075 void extproc_(int*);
076 void opensusyfile_(const char* ,int * , int);
077 }
078
079
080
081 using HepMC::IO_HEPEVT;
082 using HepMC::IO_Ascii;
083
084
085
086
087 Atlas_HEPEVT* Pythia::atlas_HEPEVT = new Atlas_HEPEVT();
088
089
090 Pythia::Pythia(const std::string& name,
091 ISvcLocator* pSvcLocator): GenModule(name,pSvcLocator)
092 {
093
094 pythia_stream = "PYTHIA_INIT";
095 declareProperty("SetAtlasDefaults", m_AtlasDefaults = true );
096 declareProperty("SusyInputFile", m_read_Filesusy = " ");
097 declareProperty("PythiaCommand", m_pythiaCommandVector);
098
099 m_firstlistevent = -1;
100 m_lastlistevent = -1;
101 }
102
103 Pythia::~Pythia(){
104
105 }
106
107
108
109
110 Pydat1& Pythia::pydat1() {
111 return m_pydat1;
112 }
113 Pydat2& Pythia::pydat2() {
114 return m_pydat2;
115 }
116 Pydat3& Pythia::pydat3() {
117 return m_pydat3;
118 }
119 Pysubs& Pythia::pysubs() {
120 return m_pysubs;
121 }
122 Pypars& Pythia::pypars() {
123 return m_pypars;
124 }
125 Pydatr& Pythia::pydatr() {
126 return m_pydatr;
127 }
128 Pymssm& Pythia::pymssm() {
129 return m_pymssm;
130 }
131 Pypued& Pythia::pypued() {
132 return m_pypued;
133 }
134 Pymsrv& Pythia::pymsrv() {
135 return m_pymsrv;
136 }
137 Pyint1& Pythia::pyint1() {
138 return m_pyint1;
139 }
140 Pyint2& Pythia::pyint2() {
141 return m_pyint2;
142 }
143 Pyint5& Pythia::pyint5() {
144 return m_pyint5;
145 }
146 Pytcsm& Pythia::pytcsm() {
147 return m_pytcsm;
148 }
149 AcerMC_acset& Pythia::acermc_acset() {
150 return m_acermc_acset;
151 }
152
153
154 StatusCode Pythia::genInitialize() {
155
156
157
158 MsgStream log(messageService(), name());
159 log << MSG::INFO << " PYTHIA INITIALISING. \n" << endreq;
160
161 static const bool CREATEIFNOTTHERE(true);
162 StatusCode RndmStatus = service("AtRndmGenSvc",
163 Pythia::p_AtRndmGenSvc,
164 CREATEIFNOTTHERE);
165 if (!RndmStatus.isSuccess() || 0 == Pythia::p_AtRndmGenSvc)
166 {
167 log << MSG::ERROR << " Could not initialize Random Number Service" << endreq;
168 return RndmStatus;
169 }
170 Pythia::pythia_stream = "PYTHIA_INIT";
171
172
173 m_frame = "CMS ";
174 m_beam = "P ";
175 m_target = "P ";
176 m_win=14000.;
177 m_initlistlevel=11;
178 m_pystatlistlevel.push_back(1);
179 m_fortout=0;
180
181 m_RndmFileName=" ";
182 m_RndmFileNumber=0;
183 m_RndmSwitch=0;
184 m_RndmFirstEvent=1;
185 m_RndmSkipEvents=0;
186 m_RndmMOVE=0;
187 m_RndmFileLength=0;
188 m_ExternalProcess = 0;
189
190
191
192 m_randomseed=19780503;
193 this->pydatr().mrpy(1) = m_randomseed;
194
195
196 m_msel=6;
197 this->pysubs().msel() = m_msel;
198
199
200 if (m_AtlasDefaults) {
201
202 this->pydat2().pmas(6,1) =175.;
203 this->pydat2().pmas(24,1) =80.42;
204 this->pydat2().pmas(24,2) =2.124;
205 this->pydat2().pmas(23,1) =91.19;
206 this->pydat2().pmas(23,2) =2.495;
207 this->pypars().mstp(128) =1;
208 this->pydat1().mstu(21) =1;
209
210 this->pypars().mstp(81) =21;
211
212
213
214
215 for (unsigned int i = 0; i < m_pythiaCommandVector.size(); i++) {
216 StringParse mystring(m_pythiaCommandVector[i]);
217 string myblock=mystring.piece(1);
218 string myentry=mystring.piece(2);
219 int myint1=mystring.intpiece(3);
220 int myint2=mystring.intpiece(4);
221 if (myblock == "pypars"){
222 if(myentry == "mstp" && myint1 == 81){
223 this->pypars().mstp(myint1)=myint2;
224 }
225 }
226 }
227
228
229
230 if (this->pypars().mstp(81) >= 20) {
231
232 this->pypars().mstp(70)=2;
233 this->pypars().mstp(72)=0;
234 this->pypars().mstp(82)=4;
235 this->pypars().mstp(84)=1;
236 this->pypars().mstp(85)=1;
237 this->pypars().mstp(86)=2;
238 this->pypars().mstp(87)=4;
239 this->pypars().mstp(88)=0;
240 this->pypars().mstp(89)=1;
241 this->pypars().mstp(90)=1;
242 this->pypars().mstp(95)=1;
243 this->pypars().parp(78)=0.2;
244 this->pypars().parp(80)=0.01;
245 this->pypars().parp(82)=1.9;
246 this->pypars().parp(83)=0.3;
247 this->pypars().parp(84)=0.5;
248 this->pypars().parp(89)=1800.;
249 this->pypars().parp(90)=0.22;
250 this->pydat1().mstj(11)=3;
251 this->pydat1().mstj(22)=2;
252 this->pydat1().parj(54)=-0.07;
253 this->pydat1().parj(55)=-0.006;
254 this->pydat1().parj(81)=0.14;
255 this->pypars().mstp(52)=2;
256 this->pypars().mstp(54)=2;
257 this->pypars().mstp(56)=2;
258 this->pypars().mstp(51)=10042;
259 this->pypars().mstp(53)=10042;
260 this->pypars().mstp(55)=10042;
261 } else {
262 this->pypars().mstp(81) = 1;
263 this->pypars().mstp(82) = 4;
264 this->pypars().mstp(86) = 2;
265 this->pypars().parp(67) = 3.;
266 this->pypars().parp(82) = 2.2;
267 this->pypars().parp(83) = 0.5;
268 this->pypars().parp(84) = 0.3;
269 this->pypars().parp(85) = 0.33;
270 this->pypars().parp(86) = 0.66;
271 this->pypars().parp(89) = 1800.;
272 this->pypars().parp(90) = 0.24;
273 this->pypars().mstp(52) = 2;
274 this->pypars().mstp(54) = 2;
275 this->pypars().mstp(56) = 2;
276 this->pypars().mstp(51) = 10042;
277 this->pypars().mstp(53) = 10042;
278 this->pypars().mstp(55) = 10042;
279 this->pydat1().mstj(11) = 3;
280 this->pydat1().mstj(22) = 2;
281 this->pydat1().parj(54) = -0.07;
282 this->pydat1().parj(55) = -0.006;
283 }
284 }
285
286
287
288
289
290
291 for (unsigned int i = 0; i < m_pythiaCommandVector.size(); i++) {
292 log << MSG::INFO << " Command is: " << m_pythiaCommandVector[i] << endreq;
293 StringParse mystring(m_pythiaCommandVector[i]);
294 string myblock=mystring.piece(1);
295 string myentry=mystring.piece(2);
296 string mystr=mystring.piece(3);
297 int myint1=mystring.intpiece(3);
298 int myint2=mystring.intpiece(4);
299 int myint3=mystring.intpiece(5);
300 int myint4=mystring.intpiece(6);
301 int myint5=mystring.intpiece(7);
302 double myfl1=mystring.numpiece(4);
303 double myfl2=mystring.numpiece(5);
304 double myfl0=mystring.numpiece(3);
305
306
307 log << MSG:: INFO << myblock << " block " << myentry << " item " << myint1 << " value " << myfl1 <<endreq;
308 if (myblock=="pyinit") {
309 if(myentry=="user"){
310 m_frame="USER ";
311 m_external=mystr;
312 m_ExternalProcess = generator_int(mystr);
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346 extproc_(&m_ExternalProcess);
347 }
348 else if(myentry=="FIXT"){
349 m_frame="FIXT";
350 }
351 else if(myentry=="pbar"){
352 m_beam="P~- ";
353 }
354 else if(myentry=="win"){
355 m_win=myfl0;
356 }
357 else if(myentry=="pylisti"){
358 m_initlistlevel=myint1;
359 }
360 else if(myentry=="pylistf"){
361 m_eventlistlevel=myint1;
362 }
363 else if(myentry=="dumpr"){
364 m_firstlistevent=myint1;
365 m_lastlistevent=myint2;
366 }
367 else if(myentry=="output"){
368 m_envval1=mystr;
369 m_fortout=49;
370 this->pydat1().mstu(11)=m_fortout;
371 }
372 else if(myentry=="rndm_IO"){
373 m_RndmFileName=mystr;
374 m_RndmFileNumber=36;
375 m_RndmSwitch=myint2;
376 m_RndmFirstEvent=myint3;
377 m_RndmSkipEvents=myint4;
378 m_RndmMOVE=myint5;
379 if (m_RndmSwitch == 1) ++m_RndmSkipEvents;
380
381 if (m_RndmSwitch > 0)
382 {
383 if (m_RndmFirstEvent < 1 || m_RndmSkipEvents < 0)
384 {
385 log << MSG:: ERROR
386 << " INCOSISTENT SET OF rndm_IO PARAMETERS : FirstEvent < 1 || SkipEvents < 0. DUMPING RNDM SEEDS SWITCHED OFF"
387 << myentry <<endreq ;
388 m_RndmSwitch = 0;
389 }
390 }
391
392 log << MSG::INFO << " !!!!!!!!!! WARNING ON PYTHIA RANDOM NUMBERS !!!!!!!! " << endreq;
393 log << MSG::INFO << " THE ATHENA/PYTHIA DOES NOT USE ANY MORE THE STANDARD " << endreq;
394 log << MSG::INFO << " PYTHIA RANDOM NUMBER SERVICES. SINCE RELEASE 5.2.0 THE" << endreq;
395 log << MSG::INFO << " ATHENA SERVICE AtRndmGenSvc IS USED. PLEASE FOR MORE " << endreq;
396 log << MSG::INFO << " DETAILS LOOK IN " << endreq;
397 log << MSG::INFO << " http://atlassw1.phy.bnl.gov/lxrsource/current/atlas/Generators/GeneratorModules/doc/Pythia.pdf "
398 << endreq;
399 log << MSG::INFO << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endreq;
400
401 m_RndmSwitch = 0;
402 }
403 else {
404 log << MSG:: ERROR << "PYTHIA ERROR, entry PYINIT has USER PBAR PYLISTI PYLISTF PYSTAT OUTPUT DUMPR WIN AND RNDM_IO: YOU HAVE SPECIFIED "<< myentry <<endreq ;
405 }
406 }
407 else if (myblock == "pysubs") {
408 if(myentry == "msel"){
409 this->pysubs().msel()=myint1;
410 }
411 else if (myentry == "msub"){
412 this->pysubs().msub(myint1)=myint2;
413 }
414 else if (myentry == "ckin"){
415 this->pysubs().ckin(myint1)=myfl1;
416 }
417 else if (myentry == "kfin"){
418 this->pysubs().kfin(myint1,myint2)=myint3;
419 }
420 else {
421 log << MSG:: ERROR << "PYTHIA ERROR, block PYSUBS has MSEL, MSUB, KFIN AND CKIN: YOU HAVE SPECIFIED "<< myentry <<endreq ;
422 }
423 }
424 else if (myblock == "pypars"){
425 if(myentry == "mstp"){
426 this->pypars().mstp(myint1)=myint2;
427 }
428 else if(myentry == "msti"){
429 this->pypars().msti(myint1)=myint2;
430 }
431 else if(myentry == "parp"){
432 this->pypars().parp(myint1)=myfl1;
433 }
434 else if(myentry == "pari"){
435 this->pypars().pari(myint1)=myfl1;
436 }
437 else {
438 log << MSG:: ERROR << "PYTHIA ERROR, block PYPARS has MSTP,PARP, MSTI AND PARI: YOU HAVE SPECIFIED "<< myentry <<endreq ;
439 }
440 }
441 else if (myblock == "pydat1"){
442 if(myentry == "mstu"){
443 this->pydat1().mstu(myint1)=myint2;
444 }
445 else if(myentry == "mstj"){
446 this->pydat1().mstj(myint1)=myint2;
447 }
448 else if(myentry == "paru"){
449 this->pydat1().paru(myint1)=myfl1;
450 }
451 else if(myentry == "parj"){
452 this->pydat1().parj(myint1)=myfl1;
453 }
454 else {
455 log << MSG:: ERROR << "PYTHIA ERROR, block PYDAT1 HAS MSTU, MSTJ, PARU AND PARJ: YOU HAVE SPECIFIED "<< myentry <<endreq ;
456 }
457 }
458 else if (myblock == "pydat2"){
459 if(myentry == "kchg"){
460 this->pydat2().kchg(myint1,myint2)=myint3;
461 }
462 else if(myentry == "pmas"){
463 this->pydat2().pmas(pycomp_(&myint1),myint2)=myfl2;
464 }
465 else if(myentry == "parf"){
466 this->pydat2().parf(myint1)=myfl1;
467 }
468 else if(myentry == "vckm"){
469 this->pydat2().vckm(myint1,myint2)=myfl2;
470 }
471 else {
472 log << MSG:: ERROR << "PYTHIA ERROR, block PYDAT2 HAS KCHG, PMAS, PARF AND VCKM: YOU HAVE SPECIFIED "<< myentry <<endreq ;
473 }
474 }
475 else if (myblock == "pydat3"){
476 if(myentry == "mdcy"){
477 this->pydat3().mdcy(myint1,myint2)=myint3;
478 }
479 else if(myentry == "mdme"){
480 this->pydat3().mdme(myint1,myint2)=myint3;
481 }
482 else if(myentry == "brat"){
483 this->pydat3().brat(myint1)=myfl1;
484 }
485 else if(myentry == "kfdp"){
486 this->pydat3().kfdp(myint1,myint2)=myint3;
487 }
488 else {
489 log << MSG:: ERROR << "PYTHIA ERROR, block PYDAT3 HAS KFDP, MDCY, BRAT AND MDME : YOU HAVE SPECIFIED "<< myentry <<endreq ;
490 }
491 }
492 else if (myblock == "pydatr"){
493 if(myentry == "mrpy"){
494 log << MSG::INFO << " !!!!!!!!!! WARNING ON PYTHIA RANDOM NUMBERS !!!!!!!! " << endreq;
495 log << MSG::INFO << " THE ATHENA/PYTHIA DOES NOT USE ANY MORE THE STANDARD " << endreq;
496 log << MSG::INFO << " PYTHIA RANDOM NUMBER SERVICES. SINCE RELEASE 5.2.0 THE" << endreq;
497 log << MSG::INFO << " ATHENA SERVICE AtRndmGenSvc IS USED. PLEASE FOR MORE " << endreq;
498 log << MSG::INFO << " DETAILS LOOK IN " << endreq;
499 log << MSG::INFO << " http://atlassw1.phy.bnl.gov/lxrsource/current/atlas/Generators/GeneratorModules/doc/Pythia.pdf "
500 << endreq;
501 log << MSG::INFO << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endreq;
502
503 this->pydatr().mrpy(myint1)=myint2;
504 }
505 else if(myentry == "rrpy"){
506 this->pydatr().rrpy(myint1)=myfl1;
507 }
508 else {
509 log << MSG:: ERROR << "PYTHIA ERROR, block PYDATR HAS MRPY AND RRPY : YOU HAVE SPECIFIED "<< myentry <<endreq ;
510 }
511 }
512 else if (myblock == "pymssm"){
513 if (myentry=="imss"){
514 if (myint1 == 21 || myint1 == 22) {
515 log << MSG::WARNING << "The seting of imss(21) and imss(22) is not allowed. When imss(1)=11 is chosen imss(21) and imss(22) are set to 66 by default" <<endreq;
516 } else {
517 this->pymssm().imss(myint1)=myint2;
518 }
519 }
520 else if (myentry=="rmss"){
521 this->pymssm().rmss(myint1)=myfl1;
522 }
523 else {
524 log << MSG:: ERROR << "PYTHIA ERROR, block PYMSSM has IMSS AND RMSS: YOU HAVE SPECIFIED "<< myentry <<endreq ;
525 }
526 }
527 else if (myblock == "pypued"){
528 if (myentry=="iued")
529 {
530 this->pypued().iued(myint1)=myint2;
531 }
532 else if (myentry=="rued"){
533 this->pypued().rued(myint1)=myfl1;
534 }
535 else {
536 log << MSG:: ERROR << "PYTHIA ERROR, block PYPUED has IUED AND RUED: YOU HAVE SPECIFIED "<< myentry <<endreq ;
537 }
538 }
539
540 else if(myblock == "pymsrv") {
541
542 if (myentry=="rvlam"){
543 this->pymsrv().rvlam(myint1)=myfl1;
544 log << MSG::INFO << "setting Yukawa coupling lambda to "<< myfl1 <<endreq ;
545 }
546
547 else if (myentry=="rvlamp"){
548 this->pymsrv().rvlamp(myint1)=myfl1;
549 log << MSG::INFO << "setting Yukawa coupling lambda^prime to "<< myfl1 <<endreq ;
550 }
551
552 else if (myentry=="rvlamb"){
553 this->pymsrv().rvlamb(myint1)=myfl1;
554 log << MSG::INFO << "setting Yukawa coupling lambda^doubelprime to "<< myfl1 <<endreq ;
555 }
556 else {
557 log << MSG:: ERROR << "PYTHIA ERROR, block PYMSRV has RVLAM, RVLAMP AND RVLAMB: YOU HAVE SPECIFIED "<< myentry <<endreq ;
558 }
559 }
560 else if (myblock == "pyint2"){
561 if (myentry == "iset"){
562 this->pyint2().iset(myint1)=myint2;
563 }
564 else if (myentry == "kfpr"){
565 this->pyint2().kfpr(myint1,myint2)=myint3;
566 }
567 else if (myentry == "coef"){
568 this->pyint2().coef(myint1,myint2)=myfl2;
569 }
570 else if (myentry == "icol"){
571 this->pyint2().icol(myint1,myint2,myint3)=myint4;
572 }
573 else {
574 log << MSG:: ERROR << "PYTHIA ERROR, block PYINT2 has ISET KFPR COEF AND ICOL: YOU HAVE SPECIFIED "<< myentry <<endreq ;
575 }
576 }
577 else if (myblock == "pystat"){
578 m_pystatlistlevel.clear();
579 for (unsigned i = 2; i <= mystring.string_size(); ++i)
580 m_pystatlistlevel.push_back(mystring.intpiece(i));
581 }
582 else if (myblock == "pytcsm"){
583 if (myentry == "itcm"){
584 this->pytcsm().itcm(myint1)=myint2;
585 }
586 else if (myentry=="rtcm"){
587 this->pytcsm().rtcm(myint1)=myfl1;
588 }
589 else {
590 log << MSG:: ERROR << "PYTHIA ERROR, block PYTCM has ITCM AND RTCM: YOU HAVE SPECIFIED "<< myentry <<endreq ;
591 }
592
593 }
594
595
596 else if (myblock == "acermc" && myentry=="acset12") {
597 if (m_ExternalProcess == ACERMC) {
598 this->acermc_acset().acset12()=myint1;
599 }
600 }
601
602 else {
603 log << MSG:: ERROR << " ERROR in PYTHIA PARAMETERS " << myblock << " is and invalid common block name !" << endreq;
604 }
605 }
606
607
608
609
610 if (this->pymssm().imss(1) == 11) {
611 int syunit = 66;
612 this->pymssm().imss(21) = syunit;
613 this->pymssm().imss(22) = syunit;
614 const std::string& fileName = m_read_Filesusy;
615 opensusyfile_(fileName.c_str(),&syunit,fileName.size());
616 }
617
618
619 if (this->pypars().mstp(81) == 0 && this->pypars().mstp(82) > 2)
620 {
621 this->pypars().mstp(82) = 1;
622 log << MSG::INFO << " !!!!!!!!!! WARNING ON PYTHIA LOOPING !!!!!!!! " << endreq;
623 log << MSG::INFO << " YOU HAVE SWITCHED OFF MULTIPLE INTERACTIONS, mstp(81) = 0 " << endreq;
624 log << MSG::INFO << " THE DEFAULT ATLAS MULTIPLE INTERACTIONS SCENARIO, mstp(82) = 4 " << endreq;
625 log << MSG::INFO << " CHANGED TO mstp(82) = 1, BECAUSE PYTHIA IS LOOPING WHEN " << endreq;
626 log << MSG::INFO << " mstp(81) = 0 and mstp(82) > 2 " << endreq;
627 log << MSG::INFO << " !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! " << endreq;
628 }
629
630
631
632
633
634
635
636
637
638
639
640 const char* envval = m_envval1.c_str();
641 initpyblock_(&m_fortout,envval);
642
643 if (m_RndmSwitch > 0)
644 {
645 const char* RndmFileName = m_RndmFileName.c_str();
646 m_RndmFileLength = openrandom_(&m_RndmSwitch, &m_RndmFileNumber, RndmFileName);
647 if (m_RndmSwitch ==2 ) log << MSG:: INFO << " THEY ARE " << m_RndmFileLength
648 << " EVENTS STORED IN THE PYTHIA RANDOM NUMBER FILE \n" << endreq;
649 }
650
651 const char* frame = m_frame.c_str();
652 const char* beam = m_beam.c_str();
653 double winval = m_win;
654 const char* target = m_target.c_str();
655 int minlist = m_initlistlevel;
656
657 if(m_external=="comphep"){
658 rinpar_();
659 subnum_();
660 rheader_();
661 }
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688 CLHEP::HepRandomEngine* engine = Pythia::p_AtRndmGenSvc->GetEngine(Pythia::pythia_stream);
689 const long* sip = engine->getSeeds();
690 long int si1 = sip[0];
691 long int si2 = sip[1];
692
693 pyinit_(frame,beam,target,&winval,
694 strlen(frame),strlen(beam),strlen(target)
695 );
696
697
698 Pythia::p_AtRndmGenSvc->CreateStream(si1, si2, Pythia::pythia_stream);
699
700 Pythia::pythia_stream = "PYTHIA";
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725 pylist_(&minlist);
726 m_events = 0;
727
728
729
730 const char *p_envar1="PYDAT";
731 char *p_envval1=0;
732 p_envval1=getenv(p_envar1);
733
734
735
736
737
738
739
740
741
742 HepMC::HEPEVT_Wrapper::set_sizeof_int(sizeof(int));
743 HepMC::HEPEVT_Wrapper::set_sizeof_real(8);
744 HepMC::HEPEVT_Wrapper::set_max_number_entries(10000);
745
746 return StatusCode::SUCCESS;
747 }
748
749
750
751 StatusCode Pythia::callGenerator() {
752
753 MsgStream log(messageService(), name());
754 log << MSG::DEBUG << " PYTHIA generating. \n" << endreq;
755
756
757 if (m_RndmSwitch > 0) RandomNumberIO();
758
759
760 CLHEP::HepRandomEngine* engine = Pythia::p_AtRndmGenSvc->GetEngine(Pythia::pythia_stream);
761 const long* s = engine->getSeeds();
762 m_seeds.clear();
763 m_seeds.push_back(s[0]);
764 m_seeds.push_back(s[1]);
765
766
767 pyevnt_();
768
769
770 ++m_events;
771 int mevlist=m_eventlistlevel;
772
773 if ( m_events >= m_firstlistevent &&
774 m_events <= m_lastlistevent) {
775 log<< MSG:: INFO << "PYEVNT event no. " << m_events << " will be listed" << endreq;
776 if (m_frame == "USER ")
777 {
778 int uspr = 7;
779 pylist_(&uspr);
780 }
781 pylist_(&mevlist);
782 }
783
784 int mconv=1;
785
786
787 pyhepc_(&mconv);
788 if ( HepMC::HEPEVT_Wrapper::number_entries() <= 0 ) {
789 std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
790 std::cout << "!!! " << generator_name(m_ExternalProcess) << " TERMINATES NORMALY: NO MORE EVENTS IN FILE !!!" << std::endl;
791 std::cout << "!!! PLEASE IGNORE ANY ATHENA ERROR MESSAGES LIKE !!!" << std::endl;
792 std::cout << "!!! AthenaEventLoopMgr ERROR Terminating event processing loop due to errors !!!" << std::endl;
793 std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl;
794 return StatusCode::FAILURE;
795 }
796 HepMC::HEPEVT_Wrapper::set_event_number(m_events);
797
798 return StatusCode::SUCCESS;
799 }
800
801
802 StatusCode Pythia::genFinalize() {
803
804 MsgStream log(messageService(), name());
805 log << MSG::INFO << " PYTHIA Ending. \n" << endreq;
806 for (std::vector<int>::iterator i = m_pystatlistlevel.begin(); i != m_pystatlistlevel.end(); ++i)
807 {
808 log << MSG::INFO <<"Call PYSTAT at endRun with level " << *i << endreq;
809 pystat_(&(*i));
810 }
811 std::cout << "MetaData: cross-section (nb)= " << 1000000. * this->pyint5().xsec(0,3) << std::endl;
812 return StatusCode::SUCCESS;
813 }
814
815 StatusCode Pythia::fillEvt(GenEvent* evt) {
816
817 MsgStream log(messageService(), name());
818
819 log << MSG::DEBUG << " PYTHIA Atlas_HEPEVT Filling. \n" << endreq;
820 store_Atlas_HEPEVT();
821
822 log << MSG::DEBUG << " PYTHIA Filling. \n" << endreq;
823 HepMC::IO_HEPEVT hepio;
824 hepio.set_print_inconsistency_errors(0);
825
826
827
828
829
830
831
832
833 hepio.fill_next_event(evt);
834 int pr_id = PYTHIA + m_ExternalProcess + this->pyint1().mint(1);
835 if (m_ExternalProcess > 0) pr_id = PYTHIA + m_ExternalProcess;
836 evt->set_signal_process_id(pr_id);
837 evt->set_random_states(m_seeds);
838 evt->weights().push_back(this->pypars().pari(7));
839 evt->weights().push_back(this->pypars().pari(9));
840 evt->weights().push_back(this->pypars().pari(10));
841
842
843
844
845
846 int id1=this->pypars().msti(15);
847 int id2=this->pypars().msti(16);
848 double x1=this->pypars().pari(33);
849 double x2=this->pypars().pari(34);
850 double q=this->pypars().pari(23);
851 double pdf1=(double)this->pypars().mstp(51);
852 double pdf2=0.;
853 HepMC::PdfInfo tmp_pdi(id1,id2,x1,x2,q,pdf1,pdf2);
854 evt->set_pdf_info(tmp_pdi);
855
856
857
858
859
860
861
862
863
864
865
866
867
868 GeVToMeV(evt);
869
870 return StatusCode::SUCCESS;
871 }
872
873 void
874 Pythia::RandomNumberIO()
875 {
876 MsgStream log(messageService(), name());
877
878
879
880 if (m_RndmSwitch == 1)
881 {
882 int WriteEvent = m_events - m_RndmFirstEvent + 1;
883 if (WriteEvent >= 0)
884 {
885 WriteEvent = WriteEvent%m_RndmSkipEvents;
886 if (WriteEvent == 0) pyrget_(&m_RndmFileNumber, &m_RndmMOVE);
887 }
888 }
889
890
891 else if (m_RndmSwitch == 2)
892 {
893 int MOVE = m_RndmMOVE;
894 if (MOVE == 0)
895 {
896 if (m_RndmFirstEvent > 0)
897 {
898 if (m_events == 0)
899 {
900 MOVE = m_RndmFirstEvent - 1;
901 }
902 else
903 {
904 MOVE = m_RndmSkipEvents;
905 }
906 }
907 int c_pos = this->pydatr().mrpy(6)+MOVE;
908 if (c_pos < m_RndmFileLength)
909 {
910 log<< MSG:: INFO << "EVENT " << m_events+1
911 << " will be generated with the random seeds read from record "
912 << c_pos+1 << " in file " << m_RndmFileName << endreq;
913 pyrset_(&m_RndmFileNumber, &MOVE);
914 }
915 }
916 else
917 {
918 pyrset_(&m_RndmFileNumber, &MOVE);
919 }
920
921 }
922 }
923
924 void
925 Pythia::store_Atlas_HEPEVT(void)
926 {
927
928
929
930
931
932
933
934 atlas_HEPEVT->fill();
935
936 Atlas_HEPEVT* Ahep = new Atlas_HEPEVT();
937 *(Ahep)=*(atlas_HEPEVT);
938 static const std::string keyid = "Pythia";
939 StatusCode sc = m_sgSvc->record(Ahep, keyid);
940 if (!sc.isSuccess()) {
941 MsgStream msg(messageService(), name());
942 msg << MSG::WARNING << " Could not record Atlas_HEPEVT" << endreq;
943 }
944
945 }
Due to the LXR bug, the updates fail sometimes to remove references to deleted files. The Saturday's full rebuilds fix these problems |
This page was automatically generated by the
LXR engine.
|
|