get tables from the aliroot directory if they are not in the current one
[u/mrichter/AliRoot.git] / THerwig / THerwig6.cxx
CommitLineData
2d8d4543 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
884ccf6a 5// Class THerwig6 is an interface to the Herwig program
6//
7// C-----------------------------------------------------------------------
8// C H E R W I G
9// 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
20// C P. Richardson(&)
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-----------------------------------------------------------------------
41// C Main reference:
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-----------------------------------------------------------------------
50
2d8d4543 51
2d8d4543 52#include "THerwig6.h"
53#include "TClonesArray.h"
54#include "TParticle.h"
884ccf6a 55#include "TObjArray.h"
2d8d4543 56
57
58ClassImp(THerwig6)
59
60extern "C" {
61 void* herwig6_common_block_address_(char*, int len);
62 void herwig6_open_fortran_file_ (int* lun, char* name, int);
63 void herwig6_close_fortran_file_(int* lun);
64}
65
66
67THerwig6::THerwig6() : TGenerator("Herwig6","Herwig6") {
68
69// THerwig6 constructor: creates a TClonesArray in which it will store all
70// particles. Note that there may be only one functional THerwig6 object
71// at a time, so it's not use to create more than one instance of it.
72
73 delete fParticles; // was allocated as TObjArray in TGenerator
74
75 fParticles = new TClonesArray("TParticle",50);
76
77 // initialize common-blocks
78 fHepevt = (Hepevt_t*) herwig6_common_block_address_((char*)"HEPEVT",6);
79 fHwbeam = (Hwbeam_t*) herwig6_common_block_address_((char*)"HWBEAM",6);
80 fHwbmch = (Hwbmch_t*) herwig6_common_block_address_((char*)"HWBMCH",6);
81 fHwproc = (Hwproc_t*) herwig6_common_block_address_((char*)"HWPROC",6);
82 fHwpram = (Hwpram_t*) herwig6_common_block_address_((char*)"HWPRAM",6);
83 fHwprch = (Hwprch_t*) herwig6_common_block_address_((char*)"HWPRCH",6);
84 fHwpart = (Hwpart_t*) herwig6_common_block_address_((char*)"HWPART",6);
85 fHwparp = (Hwparp_t*) herwig6_common_block_address_((char*)"HWPARP",6);
86 fHwbosc = (Hwbosc_t*) herwig6_common_block_address_((char*)"HWBOSC",6);
87 fHwparc = (Hwparc_t*) herwig6_common_block_address_((char*)"HWPARC",6);
88 fHwbrch = (Hwbrch_t*) herwig6_common_block_address_((char*)"HWBRCH",6);
89 fHwevnt = (Hwevnt_t*) herwig6_common_block_address_((char*)"HWEVNT",6);
90 fHwhard = (Hwhard_t*) herwig6_common_block_address_((char*)"HWHARD",6);
91 fHwprop = (Hwprop_t*) herwig6_common_block_address_((char*)"HWPROP",6);
92 fHwunam = (Hwunam_t*) herwig6_common_block_address_((char*)"HWUNAM",6);
93 fHwupdt = (Hwupdt_t*) herwig6_common_block_address_((char*)"HWUPDT",6);
94 fHwuwts = (Hwuwts_t*) herwig6_common_block_address_((char*)"HWUWTS",6);
95 fHwuclu = (Hwuclu_t*) herwig6_common_block_address_((char*)"HWUCLU",6);
96 fHwdist = (Hwdist_t*) herwig6_common_block_address_((char*)"HWDIST",6);
97 fHwqdks = (Hwqdks_t*) herwig6_common_block_address_((char*)"HWQDKS",6);
98 fHwusud = (Hwusud_t*) herwig6_common_block_address_((char*)"HWUSUD",6);
99 fHwsusy = (Hwsusy_t*) herwig6_common_block_address_((char*)"HWSUSY",6);
100 fHwrpar = (Hwrpar_t*) herwig6_common_block_address_((char*)"HWRPAR",6);
101 fHwminb = (Hwminb_t*) herwig6_common_block_address_((char*)"HWMINB",6);
102 fHwclus = (Hwclus_t*) herwig6_common_block_address_((char*)"HWCLUS",6);
103 }
104
105//------------------------------------------------------------------------------
106 THerwig6::~THerwig6()
107 {
108// Destroys the object, deletes and disposes all TParticles currently on list.
109
110 if (fParticles) {
111 fParticles->Delete();
112 delete fParticles;
113 fParticles = 0;
114 }
115 }
116
117//______________________________________________________________________________
118void THerwig6::GenerateEvent()
119{
120
121 //initialize event
122 hwuine_();
123 // generate hard subprocess
124 hwepro_();
125 // generate parton cascades
126 hwbgen_();
127 // do heavy objects decay
128 hwdhob_();
129 // do cluster formation
130 hwcfor_();
131 // do cluster decays
132 hwcdec_();
133 // do unstable particle decays
134 hwdhad_();
135 // do heavy flavor hadrons decay
136 hwdhvy_();
137 // add soft underlying event
138 hwmevt_();
139 // finish event
140 hwufne_();
141}
142
143//______________________________________________________________________________
144void THerwig6::OpenFortranFile(int lun, char* name) {
145 herwig6_open_fortran_file_(&lun, name, strlen(name));
146}
147
148//______________________________________________________________________________
149void THerwig6::CloseFortranFile(int lun) {
150 herwig6_close_fortran_file_(&lun);
151}
152
153void THerwig6::Initialize(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc)
154
155{
156 // perform the initialization for Herwig6
157 // sets correct title.
158 // after calling this method all parameters are set to their default
159 // values. If you want to modify any parameter you have to set the new
160 // value after calling Initialize and before PrepareRun.
161
162 char cbeam[8];
163 strncpy(cbeam,beam,8);
164 char ctarget[8];
165 strncpy(ctarget,target,8);
166 printf("\n Initializing Herwig !! \n");
167 if ( (!strncmp(beam, "E+" ,2)) &&
168 (!strncmp(beam, "E-" ,2)) &&
169 (!strncmp(beam, "MU+" ,3)) &&
170 (!strncmp(beam, "MU-" ,3)) &&
171 (!strncmp(beam, "NUE" ,3)) &&
172 (!strncmp(beam, "NUEB" ,4)) &&
173 (!strncmp(beam, "NUMU" ,4)) &&
174 (!strncmp(beam, "NMUB" ,4)) &&
175 (!strncmp(beam, "NTAU" ,4)) &&
176 (!strncmp(beam, "NTAB" ,4)) &&
177 (!strncmp(beam, "GAMA" ,4)) &&
178 (!strncmp(beam, "P ",8)) &&
179 (!strncmp(beam, "PBAR ",8)) &&
180 (!strncmp(beam, "N" ,1)) &&
181 (!strncmp(beam, "NBAR" ,4)) &&
182 (!strncmp(beam, "PI+" ,3)) &&
183 (!strncmp(beam, "PI-" ,3)) ) {
184 printf("WARNING! In THerwig6:Initialize():\n");
185 printf(" specified beam=%s is unrecognized .\n",beam);
186 printf(" resetting to \"P\" .");
187 sprintf(cbeam,"P");
188 }
189
190 if ( (!strncmp(target, "E+" ,2)) &&
191 (!strncmp(target, "E-" ,2)) &&
192 (!strncmp(target, "MU+" ,3)) &&
193 (!strncmp(target, "MU-" ,3)) &&
194 (!strncmp(target, "NUE" ,3)) &&
195 (!strncmp(target, "NUEB" ,4)) &&
196 (!strncmp(target, "NUMU" ,4)) &&
197 (!strncmp(target, "NMUB" ,4)) &&
198 (!strncmp(target, "NTAU" ,4)) &&
199 (!strncmp(target, "NTAB" ,4)) &&
200 (!strncmp(target, "GAMA" ,4)) &&
201 (!strncmp(target, "P ",8)) &&
202 (!strncmp(target, "PBAR ",8)) &&
203 (!strncmp(target, "N" ,1)) &&
204 (!strncmp(target, "NBAR" ,4)) &&
205 (!strncmp(target, "PI+" ,3)) &&
206 (!strncmp(target, "PI-" ,3)) ) {
207 printf("WARNING! In THerwig6:Initialize():\n");
208 printf(" specified target=%s is unrecognized .\n",target);
209 printf(" resetting to \"P\" .");
210 sprintf(ctarget,"P");
211 }
212
213 // initialization:
214 // type of beams
215 strncpy(fHwbmch->PART1,beam,8);
216 strncpy(fHwbmch->PART2,target,8);
217 // momentum of beams
218 fHwproc->PBEAM1=pbeam1;
219 fHwproc->PBEAM2=pbeam2;
220 // process to generate
221 fHwproc->IPROC=iproc;
222 // not used in the class definition
223 fHwproc->MAXEV=1;
224
225 // reset all parameters
226 hwigin_();
227
228 // set correct title
229 char atitle[132];
230 double win=pbeam1+pbeam2;
231 printf("\n %s - %s at %g GeV",beam,target,win);
232 sprintf(atitle,"%s-%s at %g GeV",cbeam,ctarget,win);
233 SetTitle(atitle);
234}
235
236void THerwig6::PrepareRun()
237{
238 // compute parameter dependent constants
239 hwuinc_();
240 // initialize elementary processes
241 hweini_();
242}
243
244//______________________________________________________________________________
245TObjArray* THerwig6::ImportParticles(Option_t *option)
246{
247//
248// Default primary creation method. It reads the /HEPEVT/ common block which
249// has been filled by the GenerateEvent method. If the event generator does
250// not use the HEPEVT common block, This routine has to be overloaded by
251// the subclasses.
252// The default action is to store only the stable particles (ISTHEP = 1)
253// This can be demanded explicitly by setting the option = "Final"
254// If the option = "All", all the particles are stored.
255//
256 fParticles->Clear();
257 Int_t numpart = fHepevt->NHEP;
258 TClonesArray &a = *((TClonesArray*)fParticles);
259 if (!strcmp(option,"") || !strcmp(option,"Final")) {
260 for (Int_t i = 0; i<=numpart; i++) {
261 if (fHepevt->ISTHEP[i] == 1) {
262//
263// Use the common block values for the TParticle constructor
264//
265 new(a[i]) TParticle(
266 fHepevt->IDHEP[i],
267 fHepevt->ISTHEP[i],
268 fHepevt->JMOHEP[i][0]-1,
269 fHepevt->JMOHEP[i][1]-1,
270 fHepevt->JDAHEP[i][0]-1,
271 fHepevt->JDAHEP[i][1]-1,
272
273 fHepevt->PHEP[i][0],
274 fHepevt->PHEP[i][1],
275 fHepevt->PHEP[i][2],
276 fHepevt->PHEP[i][3],
277 fHepevt->VHEP[i][0],
278 fHepevt->VHEP[i][1],
279 fHepevt->VHEP[i][2],
280 fHepevt->VHEP[i][3]);
281 }
282 }
283 }
284 else if (!strcmp(option,"All")) {
285 for (Int_t i = 0; i<=numpart; i++) {
286 new(a[i]) TParticle(
287 fHepevt->IDHEP[i],
288 fHepevt->ISTHEP[i],
289 fHepevt->JMOHEP[i][0]-1,
290 fHepevt->JMOHEP[i][1]-1,
291 fHepevt->JDAHEP[i][0]-1,
292 fHepevt->JDAHEP[i][1]-1,
293
294 fHepevt->PHEP[i][0],
295 fHepevt->PHEP[i][1],
296 fHepevt->PHEP[i][2],
297 fHepevt->PHEP[i][3],
298 fHepevt->VHEP[i][0],
299 fHepevt->VHEP[i][1],
300 fHepevt->VHEP[i][2],
301 fHepevt->VHEP[i][3]);
302 }
303 }
304 return fParticles;
305}
306
307//______________________________________________________________________________
308Int_t THerwig6::ImportParticles(TClonesArray *particles, Option_t *option)
309{
310//
311// Default primary creation method. It reads the /HEPEVT/ common block which
312// has been filled by the GenerateEvent method. If the event generator does
313// not use the HEPEVT common block, This routine has to be overloaded by
314// the subclasses.
315// The function loops on the generated particles and store them in
316// the TClonesArray pointed by the argument particles.
317// The default action is to store only the stable particles (ISTHEP = 1)
318// This can be demanded explicitly by setting the option = "Final"
319// If the option = "All", all the particles are stored.
320//
321 if (particles == 0) return 0;
884ccf6a 322 TClonesArray &refParticles = *particles;
323 refParticles.Clear();
2d8d4543 324 Int_t numpart = fHepevt->NHEP;
325 if (!strcmp(option,"") || !strcmp(option,"Final")) {
326 for (Int_t i = 0; i< numpart; i++) {
327 if (fHepevt->ISTHEP[i] == 1) {
328//
329// Use the common block values for the TParticle constructor
330//
884ccf6a 331 new(refParticles[i]) TParticle(
2d8d4543 332 fHepevt->IDHEP[i],
333 fHepevt->ISTHEP[i],
334 fHepevt->JMOHEP[i][0]-1,
335 fHepevt->JMOHEP[i][1]-1,
336 fHepevt->JDAHEP[i][0]-1,
337 fHepevt->JDAHEP[i][1]-1,
338
339 fHepevt->PHEP[i][0],
340 fHepevt->PHEP[i][1],
341 fHepevt->PHEP[i][2],
342 fHepevt->PHEP[i][3],
343 fHepevt->VHEP[i][0],
344 fHepevt->VHEP[i][1],
345 fHepevt->VHEP[i][2],
346 fHepevt->VHEP[i][3]);
347 }
348 }
349 }
350 else if (!strcmp(option,"All")) {
351 for (Int_t i = 0; i< numpart; i++) {
884ccf6a 352 new(refParticles[i]) TParticle(
2d8d4543 353 fHepevt->IDHEP[i],
354 fHepevt->ISTHEP[i],
355 fHepevt->JMOHEP[i][0]-1,
356 fHepevt->JMOHEP[i][1]-1,
357 fHepevt->JDAHEP[i][0]-1,
358 fHepevt->JDAHEP[i][1]-1,
359
360 fHepevt->PHEP[i][0],
361 fHepevt->PHEP[i][1],
362 fHepevt->PHEP[i][2],
363 fHepevt->PHEP[i][3],
364 fHepevt->VHEP[i][0],
365 fHepevt->VHEP[i][1],
366 fHepevt->VHEP[i][2],
367 fHepevt->VHEP[i][3]);
368 }
369 }
370 return numpart;
371}
372
373void THerwig6::Hwigin()
374{
375 hwigin_();
376}
377
378void THerwig6::Hwuinc()
379{
380 hwuinc_();
381}
382
383void THerwig6::Hwusta(char* name)
384
385{
386 hwusta_(name,8);
387}
388
389void THerwig6::Hweini()
390
391{
392 hweini_();
393}
394
395void THerwig6::Hwuine()
396
397{
398 hwuine_();
399}
400
401void THerwig6::Hwepro()
402
403{
404 hwepro_();
405}
406
407void THerwig6::Hwbgen()
408
409{
410 hwbgen_();
411}
412
413void THerwig6::Hwdhob()
414
415{
416 hwdhob_();
417}
418
419void THerwig6::Hwcfor()
420
421{
422 hwcfor_();
423}
424
425void THerwig6::Hwcdec()
426
427{
428 hwcdec_();
429}
430
431void THerwig6::Hwdhad()
432
433{
434 hwdhad_();
435}
436
437void THerwig6::Hwdhvy()
438
439{
440 hwdhvy_();
441}
442
443void THerwig6::Hwmevt()
444
445{
446 hwmevt_();
447}
448
449void THerwig6::Hwufne()
450
451{
452 hwufne_();
453}
454
455void THerwig6::Hwefin()
456
457{
458 hwefin_();
459}
460
461void THerwig6::SetupTest()
462{
463 // exampe of running herwig and generating one event
464 // after changing some options
465 Initialize("P","PBAR",900.,900.,1500);
466 // here you can set some parameters
467 SetPTMIN(15.); // Min pt in hadronic jet production
468 SetYJMIN(-4.); // Min jet rapidity
469 SetYJMAX(4.); // Max jet rapidity
470 // after you set your wished parameters
471 // herwig can do its work
472 PrepareRun();
884ccf6a 473 int nEvToGenerate=1;
474 for (int i=0;i<nEvToGenerate;i++)
2d8d4543 475 {
476 GenerateEvent();
477 // do your stuff. For ex:
884ccf6a 478 int nOfPar=GetNumberOfParticles(); // from TGenerator
479 for (int j=0; j<nOfPar; j++)
2d8d4543 480 {
481 TParticle* p=GetParticle(j);
482 // here you do whatever you want with the particle
483 p->Print();
484 };
485 };
486}
487
488
489
490
491