1 // definition of c++ Class THerwig6 to be used in ROOT
2 // this is a c++ interface to the F77 Herwig6 program
3 // author: j. g. contreras jgcn@moni.mda.cinvestav.mx
4 // date: december 22, 2000
5 // Class THerwig6 is an interface to the Herwig program
7 // C-----------------------------------------------------------------------
10 // C a Monte Carlo event generator for simulating
11 // C +---------------------------------------------------+
12 // C | Hadron Emission Reactions With Interfering Gluons |
13 // C +---------------------------------------------------+
14 // C I.G. Knowles(*), G. Marchesini(+), M.H. Seymour($) and B.R. Webber(#)
15 // C-----------------------------------------------------------------------
16 // C with Minimal Supersymmetric Standard Model Matrix Elements by
17 // C S. Moretti($) and K. Odagiri($)
18 // C-----------------------------------------------------------------------
19 // C R parity violating Supersymmetric Decays and Matrix Elements by
21 // C-----------------------------------------------------------------------
22 // C matrix element corrections to top decay and Drell-Yan type processes
23 // C by G. Corcella(+)
24 // C-----------------------------------------------------------------------
25 // C Deep Inelastic Scattering and Heavy Flavour Electroproduction by
26 // C G. Abbiendi(@) and L. Stanco(%)
27 // C-----------------------------------------------------------------------
28 // C and Jet Photoproduction in Lepton-Hadron Collisions by J. Chyla(~)
29 // C-----------------------------------------------------------------------
30 // C(*) Department of Physics & Astronomy, University of Edinburgh
31 // C(+) Dipartimento di Fisica, Universita di Milano
32 // C($) Rutherford Appleton Laboratory
33 // C(#) Cavendish Laboratory, Cambridge
34 // C(&) Department of Physics, University of Oxford
35 // C(@) Dipartimento di Fisica, Universita di Bologna
36 // C(%) Dipartimento di Fisica, Universita di Padova
37 // C(~) Institute of Physics, Prague
38 // C-----------------------------------------------------------------------
39 // C Version 6.100 - 16th December 1999
40 // C-----------------------------------------------------------------------
42 // C G.Marchesini, B.R.Webber, G.Abbiendi, I.G.Knowles, M.H.Seymour,
43 // C and L.Stanco, Computer Physics Communications 67 (1992) 465.
44 // C-----------------------------------------------------------------------
45 // C Please send e-mail about this program to one of the authors at the
46 // C following Internet addresses:
47 // C I.Knowles@ed.ac.uk Giuseppe.Marchesini@mi.infn.it
48 // C M.Seymour@rl.ac.uk webber@hep.phy.cam.ac.uk
49 // C-----------------------------------------------------------------------
54 #include "TClonesArray.h"
55 #include "TParticle.h"
56 #include "TObjArray.h"
57 #include "Riostream.h"
63 void herwig6_open_fortran_file_ (int* lun, char* name, int);
64 void herwig6_close_fortran_file_(int* lun);
67 THerwig6 *THerwig6::fgInstance = 0;
69 THerwig6::THerwig6() : TGenerator("Herwig6","Herwig6")
72 // THerwig6 constructor: creates a TClonesArray in which it will store all
73 // particles. Note that there may be only one functional THerwig6 object
74 // at a time, so it's not use to create more than one instance of it.
76 delete fParticles; // was allocated as TObjArray in TGenerator
77 fParticles = new TClonesArray("TParticle",50);
79 // initialize common-blocks
82 cout << "WARNING: creating second instance of THerwig6" << endl;
86 THerwig6::THerwig6(const THerwig6 & source): TGenerator(source)
88 Fatal("THerwig6","Copy constructor not implemented yet");
90 //------------------------------------------------------------------------------
93 // Destructor. The data members of TGenerator are delete by itself
97 //------------------------------------------------------------------------------
98 THerwig6 *THerwig6::Instance()
100 return fgInstance ? fgInstance : new THerwig6;
102 //______________________________________________________________________________
103 void THerwig6::GenerateEvent()
108 // generate hard subprocess
110 // generate parton cascades
112 // do heavy objects decay
114 // do cluster formation
118 // do unstable particle decays
120 // do heavy flavor hadrons decay
122 // add soft underlying event
127 //______________________________________________________________________________
128 void THerwig6::OpenFortranFile(int lun, char* name) {
129 herwig6_open_fortran_file_(&lun, name, strlen(name));
132 //______________________________________________________________________________
133 void THerwig6::CloseFortranFile(int lun) {
134 herwig6_close_fortran_file_(&lun);
137 void THerwig6::Initialize(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc)
140 // perform the initialization for Herwig6
141 // sets correct title.
142 // after calling this method all parameters are set to their default
143 // values. If you want to modify any parameter you have to set the new
144 // value after calling Initialize and before PrepareRun.
147 strncpy(cbeam,beam, 8);
149 strncpy(ctarget,target, 8);
150 printf("\n Initializing Herwig !! \n");
151 if ( (!strncmp(beam, "E+" ,2)) &&
152 (!strncmp(beam, "E-" ,2)) &&
153 (!strncmp(beam, "MU+" ,3)) &&
154 (!strncmp(beam, "MU-" ,3)) &&
155 (!strncmp(beam, "NUE" ,3)) &&
156 (!strncmp(beam, "NUEB" ,4)) &&
157 (!strncmp(beam, "NUMU" ,4)) &&
158 (!strncmp(beam, "NMUB" ,4)) &&
159 (!strncmp(beam, "NTAU" ,4)) &&
160 (!strncmp(beam, "NTAB" ,4)) &&
161 (!strncmp(beam, "GAMA" ,4)) &&
162 (!strncmp(beam, "P ",8)) &&
163 (!strncmp(beam, "PBAR ",8)) &&
164 (!strncmp(beam, "N" ,1)) &&
165 (!strncmp(beam, "NBAR" ,4)) &&
166 (!strncmp(beam, "PI+" ,3)) &&
167 (!strncmp(beam, "PI-" ,3)) ) {
168 printf("WARNING! In THerwig6:Initialize():\n");
169 printf(" specified beam=%s is unrecognized .\n",beam);
170 printf(" resetting to \"P\" .");
171 snprintf(cbeam, 8, "P");
174 if ( (!strncmp(target, "E+" ,2)) &&
175 (!strncmp(target, "E-" ,2)) &&
176 (!strncmp(target, "MU+" ,3)) &&
177 (!strncmp(target, "MU-" ,3)) &&
178 (!strncmp(target, "NUE" ,3)) &&
179 (!strncmp(target, "NUEB" ,4)) &&
180 (!strncmp(target, "NUMU" ,4)) &&
181 (!strncmp(target, "NMUB" ,4)) &&
182 (!strncmp(target, "NTAU" ,4)) &&
183 (!strncmp(target, "NTAB" ,4)) &&
184 (!strncmp(target, "GAMA" ,4)) &&
185 (!strncmp(target, "P ",8)) &&
186 (!strncmp(target, "PBAR ",8)) &&
187 (!strncmp(target, "N" ,1)) &&
188 (!strncmp(target, "NBAR" ,4)) &&
189 (!strncmp(target, "PI+" ,3)) &&
190 (!strncmp(target, "PI-" ,3)) ) {
191 printf("WARNING! In THerwig6:Initialize():\n");
192 printf(" specified target=%s is unrecognized .\n",target);
193 printf(" resetting to \"P\" .");
194 snprintf(ctarget,8, "P");
199 strncpy(HWBMCH.PART1,beam, 8);
200 strncpy(HWBMCH.PART2,target, 8);
202 HWPROC.PBEAM1=pbeam1;
203 HWPROC.PBEAM2=pbeam2;
204 // process to generate
206 // not used in the class definition
209 // reset all parameters
214 double win=pbeam1+pbeam2;
215 printf("\n %s - %s at %g GeV \n",beam,target,win);
216 //sprintf(atitle,"%s-%s at %g GeV",cbeam,ctarget,win);
220 void THerwig6::InitializeJimmy(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc)
223 // perform the initialization for Herwig6
224 // sets correct title.
225 // after calling this method all parameters are set to their default
226 // values. If you want to modify any parameter you have to set the new
227 // value after calling Initialize and before PrepareRun.
230 strncpy(cbeam,beam,8);
232 strncpy(ctarget,target,8);
233 printf("\n Initializing Herwig !! \n");
234 if ( (!strncmp(beam, "E+" ,2)) &&
235 (!strncmp(beam, "E-" ,2)) &&
236 (!strncmp(beam, "MU+" ,3)) &&
237 (!strncmp(beam, "MU-" ,3)) &&
238 (!strncmp(beam, "NUE" ,3)) &&
239 (!strncmp(beam, "NUEB" ,4)) &&
240 (!strncmp(beam, "NUMU" ,4)) &&
241 (!strncmp(beam, "NMUB" ,4)) &&
242 (!strncmp(beam, "NTAU" ,4)) &&
243 (!strncmp(beam, "NTAB" ,4)) &&
244 (!strncmp(beam, "GAMA" ,4)) &&
245 (!strncmp(beam, "P ",8)) &&
246 (!strncmp(beam, "PBAR ",8)) &&
247 (!strncmp(beam, "N" ,1)) &&
248 (!strncmp(beam, "NBAR" ,4)) &&
249 (!strncmp(beam, "PI+" ,3)) &&
250 (!strncmp(beam, "PI-" ,3)) ) {
251 printf("WARNING! In THerwig6:Initialize():\n");
252 printf(" specified beam=%s is unrecognized .\n",beam);
253 printf(" resetting to \"P\" .");
254 snprintf(cbeam, 8, "P");
257 if ( (!strncmp(target, "E+" ,2)) &&
258 (!strncmp(target, "E-" ,2)) &&
259 (!strncmp(target, "MU+" ,3)) &&
260 (!strncmp(target, "MU-" ,3)) &&
261 (!strncmp(target, "NUE" ,3)) &&
262 (!strncmp(target, "NUEB" ,4)) &&
263 (!strncmp(target, "NUMU" ,4)) &&
264 (!strncmp(target, "NMUB" ,4)) &&
265 (!strncmp(target, "NTAU" ,4)) &&
266 (!strncmp(target, "NTAB" ,4)) &&
267 (!strncmp(target, "GAMA" ,4)) &&
268 (!strncmp(target, "P ",8)) &&
269 (!strncmp(target, "PBAR ",8)) &&
270 (!strncmp(target, "N" ,1)) &&
271 (!strncmp(target, "NBAR" ,4)) &&
272 (!strncmp(target, "PI+" ,3)) &&
273 (!strncmp(target, "PI-" ,3)) ) {
274 printf("WARNING! In THerwig6:Initialize():\n");
275 printf(" specified target=%s is unrecognized .\n",target);
276 printf(" resetting to \"P\" .");
277 snprintf(ctarget, 8, "P");
282 strncpy(HWBMCH.PART1,beam,8);
283 strncpy(HWBMCH.PART2,target,8);
285 HWPROC.PBEAM1=pbeam1;
286 HWPROC.PBEAM2=pbeam2;
287 // process to generate
289 // not used in the class definition
292 // reset all parameters
294 // JIMMY initialization
299 double win=pbeam1+pbeam2;
300 printf("\n %s - %s at %g GeV",beam,target,win);
301 // sprintf(atitle,"%s-%s at %g GeV",cbeam,ctarget,win);
305 void THerwig6::PrepareRun()
307 // compute parameter dependent constants
309 // initialize elementary processes
313 void THerwig6::PrepareRunJimmy()
315 // compute parameter dependent constants
317 // initialize elementary processes
319 // more initializations for JIMMY
322 //______________________________________________________________________________
323 TObjArray* THerwig6::ImportParticles(Option_t *option)
326 // Default primary creation method. It reads the /HEPEVT/ common block which
327 // has been filled by the GenerateEvent method. If the event generator does
328 // not use the HEPEVT common block, This routine has to be overloaded by
330 // The default action is to store only the stable particles (ISTHEP = 1)
331 // This can be demanded explicitly by setting the option = "Final"
332 // If the option = "All", all the particles are stored.
335 Int_t numpart = HEPEVT.NHEP;
336 TClonesArray &a = *((TClonesArray*)fParticles);
337 if (!strcmp(option,"") || !strcmp(option,"Final")) {
338 for (Int_t i = 0; i < numpart; i++) {
339 if (HEPEVT.ISTHEP[i] == 1) {
341 // Use the common block values for the TParticle constructor
346 HEPEVT.JMOHEP[i][0]-1,
347 HEPEVT.JMOHEP[i][1]-1,
348 HEPEVT.JDAHEP[i][0]-1,
349 HEPEVT.JDAHEP[i][1]-1,
362 else if (!strcmp(option,"All")) {
363 for (Int_t i = 0; i < numpart; i++) {
367 HEPEVT.JMOHEP[i][0]-1,
368 HEPEVT.JMOHEP[i][1]-1,
369 HEPEVT.JDAHEP[i][0]-1,
370 HEPEVT.JDAHEP[i][1]-1,
385 //______________________________________________________________________________
386 Int_t THerwig6::ImportParticles(TClonesArray *particles, Option_t *option)
389 // Default primary creation method. It reads the /HEPEVT/ common block which
390 // has been filled by the GenerateEvent method. If the event generator does
391 // not use the HEPEVT common block, This routine has to be overloaded by
393 // The function loops on the generated particles and store them in
394 // the TClonesArray pointed by the argument particles.
395 // The default action is to store only the stable particles (ISTHEP = 1)
396 // This can be demanded explicitly by setting the option = "Final"
397 // If the option = "All", all the particles are stored.
399 if (particles == 0) return 0;
400 TClonesArray &refParticles = *particles;
401 refParticles.Clear();
402 Int_t numpart = HEPEVT.NHEP;
403 if (!strcmp(option,"") || !strcmp(option,"Final")) {
404 for (Int_t i = 0; i < numpart; i++) {
405 if (HEPEVT.ISTHEP[i] == 1) {
407 // Use the common block values for the TParticle constructor
409 new(refParticles[i]) TParticle(
412 HEPEVT.JMOHEP[i][0]-1,
413 HEPEVT.JMOHEP[i][1]-1,
414 HEPEVT.JDAHEP[i][0]-1,
415 HEPEVT.JDAHEP[i][1]-1,
428 else if (!strcmp(option,"All")) {
429 for (Int_t i = 0; i< numpart; i++) {
430 new(refParticles[i]) TParticle(
433 HEPEVT.JMOHEP[i][0]-1,
434 HEPEVT.JMOHEP[i][1]-1,
435 HEPEVT.JDAHEP[i][0]-1,
436 HEPEVT.JDAHEP[i][1]-1,
445 HEPEVT.VHEP[i][3]); //
451 void THerwig6::Hwigin()
456 void THerwig6::Hwuinc()
461 void THerwig6::Hwusta(const char* name)
467 void THerwig6::Hweini()
473 void THerwig6::Hwuine()
479 void THerwig6::Hwepro()
485 void THerwig6::Hwbgen()
491 void THerwig6::Hwdhob()
497 void THerwig6::Hwcfor()
503 void THerwig6::Hwcdec()
509 void THerwig6::Hwdhad()
515 void THerwig6::Hwdhvy()
521 void THerwig6::Hwmevt()
527 void THerwig6::Hwufne()
533 void THerwig6::Hwefin()
539 void THerwig6::Hwiodk(int iopt)
545 void THerwig6::SetupTest()
547 // exampe of running herwig and generating one event
548 // after changing some options
549 Initialize("P","PBAR",900.,900.,1500);
550 // here you can set some parameters
551 SetPTMIN(15.); // Min pt in hadronic jet production
552 SetYJMIN(-4.); // Min jet rapidity
553 SetYJMAX(4.); // Max jet rapidity
554 // after you set your wished parameters
555 // herwig can do its work
558 for (int i=0;i<nEvToGenerate;i++)
561 // do your stuff. For ex:
562 int nOfPar=GetNumberOfParticles(); // from TGenerator
563 for (int j=0; j<nOfPar; j++)
565 TParticle* p=GetParticle(j);
566 // here you do whatever you want with the particle
574 void THerwig6::Jminit()
579 void THerwig6::Jimmin()
584 void THerwig6::Jmefin()
589 void THerwig6::PrintEvt()
595 // acces to hep common block
596 int THerwig6::GetNEVHEP () const { return HEPEVT.NEVHEP; }
597 int THerwig6::GetNhep () const { return HEPEVT.NHEP; }
598 int THerwig6::GetISTHEP (int i)const { return HEPEVT.ISTHEP[i-1]; }
599 int THerwig6::GetIDHEP (int i)const { return HEPEVT.IDHEP[i-1]; }
600 int THerwig6::GetJMOHEP (int i, int j) const
601 { return HEPEVT.JMOHEP[i-1][j-1]; }
602 int THerwig6::GetJDAHEP (int i, int j) const
603 { return HEPEVT.JDAHEP[i-1][j-1]; }
604 double THerwig6::GetPHEP (int i, int j) const
605 { return HEPEVT.PHEP[i-1][j-1]; }
606 double THerwig6::GetVHEP (int i, int j) const
607 { return HEPEVT.VHEP[i-1][j-1]; }
609 // access to Herwig6 common-blocks
610 // WARNING: Some arrays start in 1, others in 0. Look up the manual!
614 int THerwig6::GetIPART1 () const { return HWBEAM.IPART1; }
615 int THerwig6::GetIPART2 () const { return HWBEAM.IPART2; }
618 char* THerwig6::GetPART1 () const { return HWBMCH.PART1; }
619 char* THerwig6::GetPART2 () const { return HWBMCH.PART2; }
623 double THerwig6::GetEBEAM1 () const { return HWPROC.EBEAM1; }
624 double THerwig6::GetEBEAM2 () const { return HWPROC.EBEAM2; }
625 double THerwig6::GetPBEAM1 () const { return HWPROC.PBEAM1; }
626 double THerwig6::GetPBEAM2 () const { return HWPROC.PBEAM2; }
627 int THerwig6::GetIPROC () const { return HWPROC.IPROC; }
628 int THerwig6::GetMAXEV () const { return HWPROC.MAXEV; }
631 double THerwig6::GetQCDLAM () const { return HWPRAM.QCDLAM; }
632 void THerwig6::SetQCDLAM (double q) const { HWPRAM.QCDLAM = q; }
633 double THerwig6::GetVQCUT () const { return HWPRAM.VQCUT; }
634 void THerwig6::SetVQCUT (double v) const { HWPRAM.VQCUT = v; }
635 double THerwig6::GetVGCUT () const { return HWPRAM.VGCUT; }
636 void THerwig6::SetVGCUT (double v) const { HWPRAM.VGCUT = v; }
637 double THerwig6::GetVPCUT () const { return HWPRAM.VPCUT; }
638 void THerwig6::SetVPCUT (double v) const { HWPRAM.VPCUT = v; }
639 double THerwig6::GetCLMAX () const { return HWPRAM.CLMAX; }
640 void THerwig6::SetCLMAX (double c) const { HWPRAM.CLMAX = c; }
641 double THerwig6::GetCLPOW () const { return HWPRAM.CLPOW; }
642 void THerwig6::SetCLPOW (double c) const { HWPRAM.CLPOW = c; }
643 double THerwig6::GetPSPLT (int i) const { return HWPRAM.PSPLT[i-1];}
644 void THerwig6::SetPSPLT (int i, double p) const { HWPRAM.PSPLT[i-1] = p;}
645 double THerwig6::GetQDIQK () const { return HWPRAM.QDIQK; }
646 void THerwig6::SetQDIQK (double q) const { HWPRAM.QDIQK = q; }
647 double THerwig6::GetPDIQK () const { return HWPRAM.PDIQK; }
648 void THerwig6::SetPDIQK (double p) const { HWPRAM.PDIQK = p; }
649 double THerwig6::GetQSPAC () const { return HWPRAM.QSPAC; }
650 void THerwig6::SetQSPAC (double q) const { HWPRAM.QSPAC = q; }
651 double THerwig6::GetPTRMS () const { return HWPRAM.PTRMS; }
652 void THerwig6::SetPTRMS (double p) const { HWPRAM.PTRMS = p; }
653 double THerwig6::GetENSOF () const { return HWPRAM.ENSOF; }
654 void THerwig6::SetENSOF (double e) const { HWPRAM.ENSOF = e; }
655 int THerwig6::GetIPRINT () const { return HWPRAM.IPRINT; }
656 void THerwig6::SetIPRINT (int i) const { HWPRAM.IPRINT = i; }
657 int THerwig6::GetMODPDF (int i) const { return HWPRAM.MODPDF[i-1];}
658 void THerwig6::SetMODPDF (int i, int j) const { HWPRAM.MODPDF[i-1] = j; }
659 int THerwig6::GetNSTRU () const { return HWPRAM.NSTRU; }
660 void THerwig6::SetNSTRU (int i) const { HWPRAM.NSTRU = i; }
663 char* THerwig6::GetAUTPDF (int i) const { return HWPRCH.AUTPDF[i-1]; }
664 void THerwig6::SetAUTPDF(int i,const char* s)const { strncpy(HWPRCH.AUTPDF[i-1], s, 19);}
665 char* THerwig6::GetBDECAY () const { return HWPRCH.BDECAY; }
668 double THerwig6::GetAVWGT () const { return HWEVNT.AVWGT; }
669 int THerwig6::GetMAXPR () const { return HWEVNT.MAXPR; }
670 void THerwig6::SetMAXPR (int i) const { HWEVNT.MAXPR = i; }
671 int THerwig6::GetMAXER () const { return HWEVNT.MAXER; }
672 void THerwig6::SetMAXER (int i) const { HWEVNT.MAXER = i; }
673 int THerwig6::GetNRN (int i) const { return HWEVNT.NRN[i-1]; }
674 void THerwig6::SetNRN (int i, int j) const { HWEVNT.NRN[i-1] = j; }
675 double THerwig6::GetEVWGT () const { return HWEVNT.EVWGT; }
677 int THerwig6::GetIDHW (int i) const { return HWEVNT.IDHW[i]; }
679 int THerwig6::GetIERROR () const { return HWEVNT.IERROR; }
682 double THerwig6::GetPTMIN () const { return HWHARD.PTMIN; }
683 void THerwig6::SetPTMIN (double d) const { HWHARD.PTMIN = d; }
684 double THerwig6::GetPTMAX () const { return HWHARD.PTMAX; }
685 void THerwig6::SetPTMAX (double d) const { HWHARD.PTMAX = d; }
686 double THerwig6::GetPTPOW () const { return HWHARD.PTPOW; }
687 void THerwig6::SetPTPOW (double d) const { HWHARD.PTPOW = d; }
688 double THerwig6::GetYJMIN () const { return HWHARD.YJMIN; }
689 void THerwig6::SetYJMIN (double d) const { HWHARD.YJMIN = d; }
690 double THerwig6::GetYJMAX () const { return HWHARD.YJMAX; }
691 void THerwig6::SetYJMAX (double d) const { HWHARD.YJMAX = d; }
692 double THerwig6::GetQ2MIN () const { return HWHARD.Q2MIN; }
693 void THerwig6::SetQ2MIN (double d) const { HWHARD.Q2MIN = d; }
694 double THerwig6::GetQ2MAX () const { return HWHARD.Q2MAX; }
695 void THerwig6::SetQ2MAX (double d) const { HWHARD.Q2MAX = d; }
696 double THerwig6::GetYBMIN () const { return HWHARD.YBMIN; }
697 void THerwig6::SetYBMIN (double d) const { HWHARD.YBMIN = d; }
698 double THerwig6::GetYBMAX () const { return HWHARD.YBMAX; }
699 void THerwig6::SetYBMAX (double d) const { HWHARD.YBMAX = d; }
700 double THerwig6::GetZJMAX () const { return HWHARD.ZJMAX; }
701 void THerwig6::SetZJMAX (double d) const { HWHARD.ZJMAX = d; }
702 int THerwig6::GetIHPRO () const { return HWHARD.IHPRO; }
704 double THerwig6::GetRMASS (int i) const { return HWPROP.RMASS[i]; }
705 void THerwig6::SetRMASS (int i, double r) const { HWPROP.RMASS[i] = r; }
708 void THerwig6::GetRNAME (int i, char a[9]) const { for (int j=0;j<8;j++) a[j] = HWUNAM.RNAME[i][j]; a[8] = '\0';}