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"
66 void herwig6_open_fortran_file_ (int* lun, char* name, int);
67 void herwig6_close_fortran_file_(int* lun);
70 THerwig6 *THerwig6::fgInstance = 0;
72 THerwig6::THerwig6() : TGenerator("Herwig6","Herwig6")
75 // THerwig6 constructor: creates a TClonesArray in which it will store all
76 // particles. Note that there may be only one functional THerwig6 object
77 // at a time, so it's not use to create more than one instance of it.
79 delete fParticles; // was allocated as TObjArray in TGenerator
80 fParticles = new TClonesArray("TParticle",50);
82 // initialize common-blocks
85 cout << "WARNING: creating second instance of THerwig6" << endl;
89 THerwig6::THerwig6(const THerwig6 & source): TGenerator(source)
91 Fatal("THerwig6","Copy constructor not implemented yet");
93 //------------------------------------------------------------------------------
96 // Destructor. The data members of TGenerator are delete by itself
100 //------------------------------------------------------------------------------
101 THerwig6 *THerwig6::Instance()
103 return fgInstance ? fgInstance : new THerwig6;
105 //______________________________________________________________________________
106 void THerwig6::GenerateEvent()
111 // generate hard subprocess
113 // generate parton cascades
115 // do heavy objects decay
117 // do cluster formation
121 // do unstable particle decays
123 // do heavy flavor hadrons decay
125 // add soft underlying event
130 //______________________________________________________________________________
131 void THerwig6::OpenFortranFile(int lun, char* name) {
132 herwig6_open_fortran_file_(&lun, name, strlen(name));
135 //______________________________________________________________________________
136 void THerwig6::CloseFortranFile(int lun) {
137 herwig6_close_fortran_file_(&lun);
140 void THerwig6::Initialize(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc)
143 // perform the initialization for Herwig6
144 // sets correct title.
145 // after calling this method all parameters are set to their default
146 // values. If you want to modify any parameter you have to set the new
147 // value after calling Initialize and before PrepareRun.
150 strncpy(cbeam,beam, 8);
152 strncpy(ctarget,target, 8);
153 printf("\n Initializing Herwig !! \n");
154 if ( (!strncmp(beam, "E+" ,2)) &&
155 (!strncmp(beam, "E-" ,2)) &&
156 (!strncmp(beam, "MU+" ,3)) &&
157 (!strncmp(beam, "MU-" ,3)) &&
158 (!strncmp(beam, "NUE" ,3)) &&
159 (!strncmp(beam, "NUEB" ,4)) &&
160 (!strncmp(beam, "NUMU" ,4)) &&
161 (!strncmp(beam, "NMUB" ,4)) &&
162 (!strncmp(beam, "NTAU" ,4)) &&
163 (!strncmp(beam, "NTAB" ,4)) &&
164 (!strncmp(beam, "GAMA" ,4)) &&
165 (!strncmp(beam, "P ",8)) &&
166 (!strncmp(beam, "PBAR ",8)) &&
167 (!strncmp(beam, "N" ,1)) &&
168 (!strncmp(beam, "NBAR" ,4)) &&
169 (!strncmp(beam, "PI+" ,3)) &&
170 (!strncmp(beam, "PI-" ,3)) ) {
171 printf("WARNING! In THerwig6:Initialize():\n");
172 printf(" specified beam=%s is unrecognized .\n",beam);
173 printf(" resetting to \"P\" .");
174 snprintf(cbeam, 8, "P");
177 if ( (!strncmp(target, "E+" ,2)) &&
178 (!strncmp(target, "E-" ,2)) &&
179 (!strncmp(target, "MU+" ,3)) &&
180 (!strncmp(target, "MU-" ,3)) &&
181 (!strncmp(target, "NUE" ,3)) &&
182 (!strncmp(target, "NUEB" ,4)) &&
183 (!strncmp(target, "NUMU" ,4)) &&
184 (!strncmp(target, "NMUB" ,4)) &&
185 (!strncmp(target, "NTAU" ,4)) &&
186 (!strncmp(target, "NTAB" ,4)) &&
187 (!strncmp(target, "GAMA" ,4)) &&
188 (!strncmp(target, "P ",8)) &&
189 (!strncmp(target, "PBAR ",8)) &&
190 (!strncmp(target, "N" ,1)) &&
191 (!strncmp(target, "NBAR" ,4)) &&
192 (!strncmp(target, "PI+" ,3)) &&
193 (!strncmp(target, "PI-" ,3)) ) {
194 printf("WARNING! In THerwig6:Initialize():\n");
195 printf(" specified target=%s is unrecognized .\n",target);
196 printf(" resetting to \"P\" .");
197 snprintf(ctarget,8, "P");
202 memcpy(HWBMCH.PART1,beam, 8);
203 memcpy(HWBMCH.PART2,target, 8);
205 HWPROC.PBEAM1=pbeam1;
206 HWPROC.PBEAM2=pbeam2;
207 // process to generate
209 // not used in the class definition
212 // reset all parameters
217 double win=pbeam1+pbeam2;
218 printf("\n %s - %s at %g GeV \n",beam,target,win);
219 //sprintf(atitle,"%s-%s at %g GeV",cbeam,ctarget,win);
223 void THerwig6::InitializeJimmy(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc)
226 // perform the initialization for Herwig6
227 // sets correct title.
228 // after calling this method all parameters are set to their default
229 // values. If you want to modify any parameter you have to set the new
230 // value after calling Initialize and before PrepareRun.
233 strncpy(cbeam,beam,8);
235 strncpy(ctarget,target,8);
236 printf("\n Initializing Herwig !! \n");
237 if ( (!strncmp(beam, "E+" ,2)) &&
238 (!strncmp(beam, "E-" ,2)) &&
239 (!strncmp(beam, "MU+" ,3)) &&
240 (!strncmp(beam, "MU-" ,3)) &&
241 (!strncmp(beam, "NUE" ,3)) &&
242 (!strncmp(beam, "NUEB" ,4)) &&
243 (!strncmp(beam, "NUMU" ,4)) &&
244 (!strncmp(beam, "NMUB" ,4)) &&
245 (!strncmp(beam, "NTAU" ,4)) &&
246 (!strncmp(beam, "NTAB" ,4)) &&
247 (!strncmp(beam, "GAMA" ,4)) &&
248 (!strncmp(beam, "P ",8)) &&
249 (!strncmp(beam, "PBAR ",8)) &&
250 (!strncmp(beam, "N" ,1)) &&
251 (!strncmp(beam, "NBAR" ,4)) &&
252 (!strncmp(beam, "PI+" ,3)) &&
253 (!strncmp(beam, "PI-" ,3)) ) {
254 printf("WARNING! In THerwig6:Initialize():\n");
255 printf(" specified beam=%s is unrecognized .\n",beam);
256 printf(" resetting to \"P\" .");
257 snprintf(cbeam, 8, "P");
260 if ( (!strncmp(target, "E+" ,2)) &&
261 (!strncmp(target, "E-" ,2)) &&
262 (!strncmp(target, "MU+" ,3)) &&
263 (!strncmp(target, "MU-" ,3)) &&
264 (!strncmp(target, "NUE" ,3)) &&
265 (!strncmp(target, "NUEB" ,4)) &&
266 (!strncmp(target, "NUMU" ,4)) &&
267 (!strncmp(target, "NMUB" ,4)) &&
268 (!strncmp(target, "NTAU" ,4)) &&
269 (!strncmp(target, "NTAB" ,4)) &&
270 (!strncmp(target, "GAMA" ,4)) &&
271 (!strncmp(target, "P ",8)) &&
272 (!strncmp(target, "PBAR ",8)) &&
273 (!strncmp(target, "N" ,1)) &&
274 (!strncmp(target, "NBAR" ,4)) &&
275 (!strncmp(target, "PI+" ,3)) &&
276 (!strncmp(target, "PI-" ,3)) ) {
277 printf("WARNING! In THerwig6:Initialize():\n");
278 printf(" specified target=%s is unrecognized .\n",target);
279 printf(" resetting to \"P\" .");
280 snprintf(ctarget, 8, "P");
285 memcpy(HWBMCH.PART1,beam,8);
286 memcpy(HWBMCH.PART2,target,8);
288 HWPROC.PBEAM1=pbeam1;
289 HWPROC.PBEAM2=pbeam2;
290 // process to generate
292 // not used in the class definition
295 // reset all parameters
297 // JIMMY initialization
302 double win=pbeam1+pbeam2;
303 printf("\n %s - %s at %g GeV",beam,target,win);
304 // sprintf(atitle,"%s-%s at %g GeV",cbeam,ctarget,win);
308 void THerwig6::PrepareRun()
310 // compute parameter dependent constants
312 // initialize elementary processes
316 void THerwig6::PrepareRunJimmy()
318 // compute parameter dependent constants
320 // initialize elementary processes
322 // more initializations for JIMMY
325 //______________________________________________________________________________
326 TObjArray* THerwig6::ImportParticles(Option_t *option)
329 // Default primary creation method. It reads the /HEPEVT/ common block which
330 // has been filled by the GenerateEvent method. If the event generator does
331 // not use the HEPEVT common block, This routine has to be overloaded by
333 // The default action is to store only the stable particles (ISTHEP = 1)
334 // This can be demanded explicitly by setting the option = "Final"
335 // If the option = "All", all the particles are stored.
338 Int_t numpart = HEPEVT.NHEP;
339 TClonesArray &a = *((TClonesArray*)fParticles);
340 if (!strcmp(option,"") || !strcmp(option,"Final")) {
341 for (Int_t i = 0; i < numpart; i++) {
342 if (HEPEVT.ISTHEP[i] == 1) {
344 // Use the common block values for the TParticle constructor
349 HEPEVT.JMOHEP[i][0]-1,
350 HEPEVT.JMOHEP[i][1]-1,
351 HEPEVT.JDAHEP[i][0]-1,
352 HEPEVT.JDAHEP[i][1]-1,
365 else if (!strcmp(option,"All")) {
366 for (Int_t i = 0; i < numpart; i++) {
370 HEPEVT.JMOHEP[i][0]-1,
371 HEPEVT.JMOHEP[i][1]-1,
372 HEPEVT.JDAHEP[i][0]-1,
373 HEPEVT.JDAHEP[i][1]-1,
388 //______________________________________________________________________________
389 Int_t THerwig6::ImportParticles(TClonesArray *particles, Option_t *option)
392 // Default primary creation method. It reads the /HEPEVT/ common block which
393 // has been filled by the GenerateEvent method. If the event generator does
394 // not use the HEPEVT common block, This routine has to be overloaded by
396 // The function loops on the generated particles and store them in
397 // the TClonesArray pointed by the argument particles.
398 // The default action is to store only the stable particles (ISTHEP = 1)
399 // This can be demanded explicitly by setting the option = "Final"
400 // If the option = "All", all the particles are stored.
402 if (particles == 0) return 0;
403 TClonesArray &refParticles = *particles;
404 refParticles.Clear();
405 Int_t numpart = HEPEVT.NHEP;
406 if (!strcmp(option,"") || !strcmp(option,"Final")) {
407 for (Int_t i = 0; i < numpart; i++) {
408 if (HEPEVT.ISTHEP[i] == 1) {
410 // Use the common block values for the TParticle constructor
412 new(refParticles[i]) TParticle(
415 HEPEVT.JMOHEP[i][0]-1,
416 HEPEVT.JMOHEP[i][1]-1,
417 HEPEVT.JDAHEP[i][0]-1,
418 HEPEVT.JDAHEP[i][1]-1,
431 else if (!strcmp(option,"All")) {
432 for (Int_t i = 0; i< numpart; i++) {
433 new(refParticles[i]) TParticle(
436 HEPEVT.JMOHEP[i][0]-1,
437 HEPEVT.JMOHEP[i][1]-1,
438 HEPEVT.JDAHEP[i][0]-1,
439 HEPEVT.JDAHEP[i][1]-1,
448 HEPEVT.VHEP[i][3]); //
454 void THerwig6::Hwigin()
459 void THerwig6::Hwuinc()
464 void THerwig6::Hwusta(const char* name)
470 void THerwig6::Hweini()
476 void THerwig6::Hwuine()
482 void THerwig6::Hwepro()
488 void THerwig6::Hwbgen()
494 void THerwig6::Hwdhob()
500 void THerwig6::Hwcfor()
506 void THerwig6::Hwcdec()
512 void THerwig6::Hwdhad()
518 void THerwig6::Hwdhvy()
524 void THerwig6::Hwmevt()
530 void THerwig6::Hwufne()
536 void THerwig6::Hwefin()
542 void THerwig6::Hwiodk(int iopt)
548 void THerwig6::SetupTest()
550 // exampe of running herwig and generating one event
551 // after changing some options
552 Initialize("P","PBAR",900.,900.,1500);
553 // here you can set some parameters
554 SetPTMIN(15.); // Min pt in hadronic jet production
555 SetYJMIN(-4.); // Min jet rapidity
556 SetYJMAX(4.); // Max jet rapidity
557 // after you set your wished parameters
558 // herwig can do its work
561 for (int i=0;i<nEvToGenerate;i++)
564 // do your stuff. For ex:
565 int nOfPar=GetNumberOfParticles(); // from TGenerator
566 for (int j=0; j<nOfPar; j++)
568 TParticle* p=GetParticle(j);
569 // here you do whatever you want with the particle
577 void THerwig6::Jminit()
582 void THerwig6::Jimmin()
587 void THerwig6::Jmefin()
592 void THerwig6::PrintEvt()
598 // acces to hep common block
599 int THerwig6::GetNEVHEP () const { return HEPEVT.NEVHEP; }
600 int THerwig6::GetNhep () const { return HEPEVT.NHEP; }
601 int THerwig6::GetISTHEP (int i)const { return HEPEVT.ISTHEP[i-1]; }
602 int THerwig6::GetIDHEP (int i)const { return HEPEVT.IDHEP[i-1]; }
603 int THerwig6::GetJMOHEP (int i, int j) const
604 { return HEPEVT.JMOHEP[i-1][j-1]; }
605 int THerwig6::GetJDAHEP (int i, int j) const
606 { return HEPEVT.JDAHEP[i-1][j-1]; }
607 double THerwig6::GetPHEP (int i, int j) const
608 { return HEPEVT.PHEP[i-1][j-1]; }
609 double THerwig6::GetVHEP (int i, int j) const
610 { return HEPEVT.VHEP[i-1][j-1]; }
612 // access to Herwig6 common-blocks
613 // WARNING: Some arrays start in 1, others in 0. Look up the manual!
617 int THerwig6::GetIPART1 () const { return HWBEAM.IPART1; }
618 int THerwig6::GetIPART2 () const { return HWBEAM.IPART2; }
621 char* THerwig6::GetPART1 () const { return HWBMCH.PART1; }
622 char* THerwig6::GetPART2 () const { return HWBMCH.PART2; }
626 double THerwig6::GetEBEAM1 () const { return HWPROC.EBEAM1; }
627 double THerwig6::GetEBEAM2 () const { return HWPROC.EBEAM2; }
628 double THerwig6::GetPBEAM1 () const { return HWPROC.PBEAM1; }
629 double THerwig6::GetPBEAM2 () const { return HWPROC.PBEAM2; }
630 int THerwig6::GetIPROC () const { return HWPROC.IPROC; }
631 int THerwig6::GetMAXEV () const { return HWPROC.MAXEV; }
634 double THerwig6::GetQCDLAM () const { return HWPRAM.QCDLAM; }
635 void THerwig6::SetQCDLAM (double q) const { HWPRAM.QCDLAM = q; }
636 double THerwig6::GetVQCUT () const { return HWPRAM.VQCUT; }
637 void THerwig6::SetVQCUT (double v) const { HWPRAM.VQCUT = v; }
638 double THerwig6::GetVGCUT () const { return HWPRAM.VGCUT; }
639 void THerwig6::SetVGCUT (double v) const { HWPRAM.VGCUT = v; }
640 double THerwig6::GetVPCUT () const { return HWPRAM.VPCUT; }
641 void THerwig6::SetVPCUT (double v) const { HWPRAM.VPCUT = v; }
642 double THerwig6::GetCLMAX () const { return HWPRAM.CLMAX; }
643 void THerwig6::SetCLMAX (double c) const { HWPRAM.CLMAX = c; }
644 double THerwig6::GetCLPOW () const { return HWPRAM.CLPOW; }
645 void THerwig6::SetCLPOW (double c) const { HWPRAM.CLPOW = c; }
646 double THerwig6::GetPSPLT (int i) const { return HWPRAM.PSPLT[i-1];}
647 void THerwig6::SetPSPLT (int i, double p) const { HWPRAM.PSPLT[i-1] = p;}
648 double THerwig6::GetQDIQK () const { return HWPRAM.QDIQK; }
649 void THerwig6::SetQDIQK (double q) const { HWPRAM.QDIQK = q; }
650 double THerwig6::GetPDIQK () const { return HWPRAM.PDIQK; }
651 void THerwig6::SetPDIQK (double p) const { HWPRAM.PDIQK = p; }
652 double THerwig6::GetQSPAC () const { return HWPRAM.QSPAC; }
653 void THerwig6::SetQSPAC (double q) const { HWPRAM.QSPAC = q; }
654 double THerwig6::GetPTRMS () const { return HWPRAM.PTRMS; }
655 void THerwig6::SetPTRMS (double p) const { HWPRAM.PTRMS = p; }
656 double THerwig6::GetENSOF () const { return HWPRAM.ENSOF; }
657 void THerwig6::SetENSOF (double e) const { HWPRAM.ENSOF = e; }
658 int THerwig6::GetIPRINT () const { return HWPRAM.IPRINT; }
659 void THerwig6::SetIPRINT (int i) const { HWPRAM.IPRINT = i; }
660 int THerwig6::GetMODPDF (int i) const { return HWPRAM.MODPDF[i-1];}
661 void THerwig6::SetMODPDF (int i, int j) const { HWPRAM.MODPDF[i-1] = j; }
662 int THerwig6::GetNSTRU () const { return HWPRAM.NSTRU; }
663 void THerwig6::SetNSTRU (int i) const { HWPRAM.NSTRU = i; }
666 char* THerwig6::GetAUTPDF (int i) const { return HWPRCH.AUTPDF[i-1]; }
667 void THerwig6::SetAUTPDF(int i,const char* s)const { strncpy(HWPRCH.AUTPDF[i-1], s, 19);}
668 char* THerwig6::GetBDECAY () const { return HWPRCH.BDECAY; }
671 double THerwig6::GetAVWGT () const { return HWEVNT.AVWGT; }
672 int THerwig6::GetMAXPR () const { return HWEVNT.MAXPR; }
673 void THerwig6::SetMAXPR (int i) const { HWEVNT.MAXPR = i; }
674 int THerwig6::GetMAXER () const { return HWEVNT.MAXER; }
675 void THerwig6::SetMAXER (int i) const { HWEVNT.MAXER = i; }
676 int THerwig6::GetNRN (int i) const { return HWEVNT.NRN[i-1]; }
677 void THerwig6::SetNRN (int i, int j) const { HWEVNT.NRN[i-1] = j; }
678 double THerwig6::GetEVWGT () const { return HWEVNT.EVWGT; }
680 int THerwig6::GetIDHW (int i) const { return HWEVNT.IDHW[i]; }
682 int THerwig6::GetIERROR () const { return HWEVNT.IERROR; }
685 double THerwig6::GetPTMIN () const { return HWHARD.PTMIN; }
686 void THerwig6::SetPTMIN (double d) const { HWHARD.PTMIN = d; }
687 double THerwig6::GetPTMAX () const { return HWHARD.PTMAX; }
688 void THerwig6::SetPTMAX (double d) const { HWHARD.PTMAX = d; }
689 double THerwig6::GetPTPOW () const { return HWHARD.PTPOW; }
690 void THerwig6::SetPTPOW (double d) const { HWHARD.PTPOW = d; }
691 double THerwig6::GetYJMIN () const { return HWHARD.YJMIN; }
692 void THerwig6::SetYJMIN (double d) const { HWHARD.YJMIN = d; }
693 double THerwig6::GetYJMAX () const { return HWHARD.YJMAX; }
694 void THerwig6::SetYJMAX (double d) const { HWHARD.YJMAX = d; }
695 double THerwig6::GetQ2MIN () const { return HWHARD.Q2MIN; }
696 void THerwig6::SetQ2MIN (double d) const { HWHARD.Q2MIN = d; }
697 double THerwig6::GetQ2MAX () const { return HWHARD.Q2MAX; }
698 void THerwig6::SetQ2MAX (double d) const { HWHARD.Q2MAX = d; }
699 double THerwig6::GetYBMIN () const { return HWHARD.YBMIN; }
700 void THerwig6::SetYBMIN (double d) const { HWHARD.YBMIN = d; }
701 double THerwig6::GetYBMAX () const { return HWHARD.YBMAX; }
702 void THerwig6::SetYBMAX (double d) const { HWHARD.YBMAX = d; }
703 double THerwig6::GetZJMAX () const { return HWHARD.ZJMAX; }
704 void THerwig6::SetZJMAX (double d) const { HWHARD.ZJMAX = d; }
705 int THerwig6::GetIHPRO () const { return HWHARD.IHPRO; }
707 double THerwig6::GetRMASS (int i) const { return HWPROP.RMASS[i]; }
708 void THerwig6::SetRMASS (int i, double r) const { HWPROP.RMASS[i] = r; }
711 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';}