New TPC parameters
[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);
e2054d85 103 fHwgrav = (Hwgrav_t*) herwig6_common_block_address_((char*)"HWGRAV",6);
104 fHw6202 = (Hw6202_t*) herwig6_common_block_address_((char*)"HW6202",6);
105 fHw6203 = (Hw6203_t*) herwig6_common_block_address_((char*)"HW6203",6);
106 fHw6300 = (Hw6300_t*) herwig6_common_block_address_((char*)"HW6300",6);
107 fHwpmrs = (Hwpmrs_t*) herwig6_common_block_address_((char*)"HWPMRS",6);
108 fHwcirc = (Hwcirc_t*) herwig6_common_block_address_((char*)"HWCIRC",6);
109 fHwdspb = (Hwdspb_t*) herwig6_common_block_address_((char*)"HWDSPB",6);
110 fHwdsp2 = (Hwdsp2_t*) herwig6_common_block_address_((char*)"HWDSP2",6);
111 fHwdsp3 = (Hwdsp3_t*) herwig6_common_block_address_((char*)"HWDSP3",6);
112 fHwdsp4 = (Hwdsp4_t*) herwig6_common_block_address_((char*)"HWDSP4",6);
113 fHwdspn = (Hwdspn_t*) herwig6_common_block_address_((char*)"HWDSPN",6);
114 fHwspin = (Hwspin_t*) herwig6_common_block_address_((char*)"HWSPIN",6);
115 fHwstau = (Hwstau_t*) herwig6_common_block_address_((char*)"HWSTAU",6);
116 fHwgupr = (Hwgupr_t*) herwig6_common_block_address_((char*)"HWGUPR",6);
117 fHw6500 = (Hw6500_t*) herwig6_common_block_address_((char*)"HW6500",6);
118 fHw6504 = (Hw6504_t*) herwig6_common_block_address_((char*)"HW6504",6);
119 fHw6506 = (Hw6506_t*) herwig6_common_block_address_((char*)"HW6506",6);
2d8d4543 120 }
121
122//------------------------------------------------------------------------------
123 THerwig6::~THerwig6()
124 {
bbc19b09 125 // Destructor. The data members of TGenerator are delete by itself
2d8d4543 126 }
127
128//______________________________________________________________________________
129void THerwig6::GenerateEvent()
130{
131
132 //initialize event
133 hwuine_();
134 // generate hard subprocess
135 hwepro_();
136 // generate parton cascades
137 hwbgen_();
138 // do heavy objects decay
139 hwdhob_();
140 // do cluster formation
141 hwcfor_();
142 // do cluster decays
143 hwcdec_();
144 // do unstable particle decays
145 hwdhad_();
146 // do heavy flavor hadrons decay
147 hwdhvy_();
148 // add soft underlying event
149 hwmevt_();
150 // finish event
151 hwufne_();
152}
153
154//______________________________________________________________________________
155void THerwig6::OpenFortranFile(int lun, char* name) {
156 herwig6_open_fortran_file_(&lun, name, strlen(name));
157}
158
159//______________________________________________________________________________
160void THerwig6::CloseFortranFile(int lun) {
161 herwig6_close_fortran_file_(&lun);
162}
163
164void THerwig6::Initialize(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc)
165
166{
167 // perform the initialization for Herwig6
168 // sets correct title.
169 // after calling this method all parameters are set to their default
170 // values. If you want to modify any parameter you have to set the new
171 // value after calling Initialize and before PrepareRun.
172
173 char cbeam[8];
174 strncpy(cbeam,beam,8);
175 char ctarget[8];
176 strncpy(ctarget,target,8);
177 printf("\n Initializing Herwig !! \n");
178 if ( (!strncmp(beam, "E+" ,2)) &&
179 (!strncmp(beam, "E-" ,2)) &&
180 (!strncmp(beam, "MU+" ,3)) &&
181 (!strncmp(beam, "MU-" ,3)) &&
182 (!strncmp(beam, "NUE" ,3)) &&
183 (!strncmp(beam, "NUEB" ,4)) &&
184 (!strncmp(beam, "NUMU" ,4)) &&
185 (!strncmp(beam, "NMUB" ,4)) &&
186 (!strncmp(beam, "NTAU" ,4)) &&
187 (!strncmp(beam, "NTAB" ,4)) &&
188 (!strncmp(beam, "GAMA" ,4)) &&
189 (!strncmp(beam, "P ",8)) &&
190 (!strncmp(beam, "PBAR ",8)) &&
191 (!strncmp(beam, "N" ,1)) &&
192 (!strncmp(beam, "NBAR" ,4)) &&
193 (!strncmp(beam, "PI+" ,3)) &&
194 (!strncmp(beam, "PI-" ,3)) ) {
195 printf("WARNING! In THerwig6:Initialize():\n");
196 printf(" specified beam=%s is unrecognized .\n",beam);
197 printf(" resetting to \"P\" .");
198 sprintf(cbeam,"P");
199 }
200
201 if ( (!strncmp(target, "E+" ,2)) &&
202 (!strncmp(target, "E-" ,2)) &&
203 (!strncmp(target, "MU+" ,3)) &&
204 (!strncmp(target, "MU-" ,3)) &&
205 (!strncmp(target, "NUE" ,3)) &&
206 (!strncmp(target, "NUEB" ,4)) &&
207 (!strncmp(target, "NUMU" ,4)) &&
208 (!strncmp(target, "NMUB" ,4)) &&
209 (!strncmp(target, "NTAU" ,4)) &&
210 (!strncmp(target, "NTAB" ,4)) &&
211 (!strncmp(target, "GAMA" ,4)) &&
212 (!strncmp(target, "P ",8)) &&
213 (!strncmp(target, "PBAR ",8)) &&
214 (!strncmp(target, "N" ,1)) &&
215 (!strncmp(target, "NBAR" ,4)) &&
216 (!strncmp(target, "PI+" ,3)) &&
217 (!strncmp(target, "PI-" ,3)) ) {
218 printf("WARNING! In THerwig6:Initialize():\n");
219 printf(" specified target=%s is unrecognized .\n",target);
220 printf(" resetting to \"P\" .");
221 sprintf(ctarget,"P");
222 }
223
224 // initialization:
225 // type of beams
226 strncpy(fHwbmch->PART1,beam,8);
227 strncpy(fHwbmch->PART2,target,8);
228 // momentum of beams
229 fHwproc->PBEAM1=pbeam1;
230 fHwproc->PBEAM2=pbeam2;
231 // process to generate
232 fHwproc->IPROC=iproc;
233 // not used in the class definition
234 fHwproc->MAXEV=1;
235
236 // reset all parameters
237 hwigin_();
238
239 // set correct title
240 char atitle[132];
241 double win=pbeam1+pbeam2;
242 printf("\n %s - %s at %g GeV",beam,target,win);
243 sprintf(atitle,"%s-%s at %g GeV",cbeam,ctarget,win);
244 SetTitle(atitle);
245}
246
247void THerwig6::PrepareRun()
248{
249 // compute parameter dependent constants
250 hwuinc_();
251 // initialize elementary processes
252 hweini_();
253}
254
255//______________________________________________________________________________
256TObjArray* THerwig6::ImportParticles(Option_t *option)
257{
258//
259// Default primary creation method. It reads the /HEPEVT/ common block which
260// has been filled by the GenerateEvent method. If the event generator does
261// not use the HEPEVT common block, This routine has to be overloaded by
262// the subclasses.
263// The default action is to store only the stable particles (ISTHEP = 1)
264// This can be demanded explicitly by setting the option = "Final"
265// If the option = "All", all the particles are stored.
266//
267 fParticles->Clear();
268 Int_t numpart = fHepevt->NHEP;
269 TClonesArray &a = *((TClonesArray*)fParticles);
270 if (!strcmp(option,"") || !strcmp(option,"Final")) {
271 for (Int_t i = 0; i<=numpart; i++) {
272 if (fHepevt->ISTHEP[i] == 1) {
273//
274// Use the common block values for the TParticle constructor
275//
276 new(a[i]) TParticle(
277 fHepevt->IDHEP[i],
278 fHepevt->ISTHEP[i],
279 fHepevt->JMOHEP[i][0]-1,
280 fHepevt->JMOHEP[i][1]-1,
281 fHepevt->JDAHEP[i][0]-1,
282 fHepevt->JDAHEP[i][1]-1,
283
284 fHepevt->PHEP[i][0],
285 fHepevt->PHEP[i][1],
286 fHepevt->PHEP[i][2],
287 fHepevt->PHEP[i][3],
288 fHepevt->VHEP[i][0],
289 fHepevt->VHEP[i][1],
290 fHepevt->VHEP[i][2],
291 fHepevt->VHEP[i][3]);
292 }
293 }
294 }
295 else if (!strcmp(option,"All")) {
296 for (Int_t i = 0; i<=numpart; i++) {
297 new(a[i]) TParticle(
298 fHepevt->IDHEP[i],
299 fHepevt->ISTHEP[i],
300 fHepevt->JMOHEP[i][0]-1,
301 fHepevt->JMOHEP[i][1]-1,
302 fHepevt->JDAHEP[i][0]-1,
303 fHepevt->JDAHEP[i][1]-1,
304
305 fHepevt->PHEP[i][0],
306 fHepevt->PHEP[i][1],
307 fHepevt->PHEP[i][2],
308 fHepevt->PHEP[i][3],
309 fHepevt->VHEP[i][0],
310 fHepevt->VHEP[i][1],
311 fHepevt->VHEP[i][2],
312 fHepevt->VHEP[i][3]);
313 }
314 }
315 return fParticles;
316}
317
318//______________________________________________________________________________
319Int_t THerwig6::ImportParticles(TClonesArray *particles, Option_t *option)
320{
321//
322// Default primary creation method. It reads the /HEPEVT/ common block which
323// has been filled by the GenerateEvent method. If the event generator does
324// not use the HEPEVT common block, This routine has to be overloaded by
325// the subclasses.
326// The function loops on the generated particles and store them in
327// the TClonesArray pointed by the argument particles.
328// The default action is to store only the stable particles (ISTHEP = 1)
329// This can be demanded explicitly by setting the option = "Final"
330// If the option = "All", all the particles are stored.
331//
332 if (particles == 0) return 0;
884ccf6a 333 TClonesArray &refParticles = *particles;
334 refParticles.Clear();
2d8d4543 335 Int_t numpart = fHepevt->NHEP;
336 if (!strcmp(option,"") || !strcmp(option,"Final")) {
337 for (Int_t i = 0; i< numpart; i++) {
338 if (fHepevt->ISTHEP[i] == 1) {
339//
340// Use the common block values for the TParticle constructor
341//
884ccf6a 342 new(refParticles[i]) TParticle(
2d8d4543 343 fHepevt->IDHEP[i],
344 fHepevt->ISTHEP[i],
345 fHepevt->JMOHEP[i][0]-1,
346 fHepevt->JMOHEP[i][1]-1,
347 fHepevt->JDAHEP[i][0]-1,
348 fHepevt->JDAHEP[i][1]-1,
349
350 fHepevt->PHEP[i][0],
351 fHepevt->PHEP[i][1],
352 fHepevt->PHEP[i][2],
353 fHepevt->PHEP[i][3],
354 fHepevt->VHEP[i][0],
355 fHepevt->VHEP[i][1],
356 fHepevt->VHEP[i][2],
357 fHepevt->VHEP[i][3]);
358 }
359 }
360 }
361 else if (!strcmp(option,"All")) {
362 for (Int_t i = 0; i< numpart; i++) {
884ccf6a 363 new(refParticles[i]) TParticle(
2d8d4543 364 fHepevt->IDHEP[i],
365 fHepevt->ISTHEP[i],
366 fHepevt->JMOHEP[i][0]-1,
367 fHepevt->JMOHEP[i][1]-1,
368 fHepevt->JDAHEP[i][0]-1,
369 fHepevt->JDAHEP[i][1]-1,
370
371 fHepevt->PHEP[i][0],
372 fHepevt->PHEP[i][1],
373 fHepevt->PHEP[i][2],
374 fHepevt->PHEP[i][3],
375 fHepevt->VHEP[i][0],
376 fHepevt->VHEP[i][1],
377 fHepevt->VHEP[i][2],
378 fHepevt->VHEP[i][3]);
379 }
380 }
381 return numpart;
382}
383
384void THerwig6::Hwigin()
385{
386 hwigin_();
387}
388
389void THerwig6::Hwuinc()
390{
391 hwuinc_();
392}
393
394void THerwig6::Hwusta(char* name)
395
396{
397 hwusta_(name,8);
398}
399
400void THerwig6::Hweini()
401
402{
403 hweini_();
404}
405
406void THerwig6::Hwuine()
407
408{
409 hwuine_();
410}
411
412void THerwig6::Hwepro()
413
414{
415 hwepro_();
416}
417
418void THerwig6::Hwbgen()
419
420{
421 hwbgen_();
422}
423
424void THerwig6::Hwdhob()
425
426{
427 hwdhob_();
428}
429
430void THerwig6::Hwcfor()
431
432{
433 hwcfor_();
434}
435
436void THerwig6::Hwcdec()
437
438{
439 hwcdec_();
440}
441
442void THerwig6::Hwdhad()
443
444{
445 hwdhad_();
446}
447
448void THerwig6::Hwdhvy()
449
450{
451 hwdhvy_();
452}
453
454void THerwig6::Hwmevt()
455
456{
457 hwmevt_();
458}
459
460void THerwig6::Hwufne()
461
462{
463 hwufne_();
464}
465
466void THerwig6::Hwefin()
467
468{
469 hwefin_();
470}
471
e2054d85 472void THerwig6::Hwiodk(int iopt)
473
474{
475 hwiodk_(iopt);
476}
477
2d8d4543 478void THerwig6::SetupTest()
479{
480 // exampe of running herwig and generating one event
481 // after changing some options
482 Initialize("P","PBAR",900.,900.,1500);
483 // here you can set some parameters
484 SetPTMIN(15.); // Min pt in hadronic jet production
485 SetYJMIN(-4.); // Min jet rapidity
486 SetYJMAX(4.); // Max jet rapidity
487 // after you set your wished parameters
488 // herwig can do its work
489 PrepareRun();
884ccf6a 490 int nEvToGenerate=1;
491 for (int i=0;i<nEvToGenerate;i++)
2d8d4543 492 {
493 GenerateEvent();
494 // do your stuff. For ex:
884ccf6a 495 int nOfPar=GetNumberOfParticles(); // from TGenerator
496 for (int j=0; j<nOfPar; j++)
2d8d4543 497 {
498 TParticle* p=GetParticle(j);
499 // here you do whatever you want with the particle
500 p->Print();
501 };
502 };
503}
504
505
506
507
508