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