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