Interface to HERWIG
[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
7#include <stream.h>
8#include "THerwig6.h"
9#include "TClonesArray.h"
10#include "TParticle.h"
11
12
13ClassImp(THerwig6)
14
15extern "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
22THerwig6::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//______________________________________________________________________________
73void 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//______________________________________________________________________________
99void THerwig6::OpenFortranFile(int lun, char* name) {
100 herwig6_open_fortran_file_(&lun, name, strlen(name));
101}
102
103//______________________________________________________________________________
104void THerwig6::CloseFortranFile(int lun) {
105 herwig6_close_fortran_file_(&lun);
106}
107
108void 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
191void THerwig6::PrepareRun()
192{
193 // compute parameter dependent constants
194 hwuinc_();
195 // initialize elementary processes
196 hweini_();
197}
198
199//______________________________________________________________________________
200TObjArray* 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//______________________________________________________________________________
263Int_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
328void THerwig6::Hwigin()
329{
330 hwigin_();
331}
332
333void THerwig6::Hwuinc()
334{
335 hwuinc_();
336}
337
338void THerwig6::Hwusta(char* name)
339
340{
341 hwusta_(name,8);
342}
343
344void THerwig6::Hweini()
345
346{
347 hweini_();
348}
349
350void THerwig6::Hwuine()
351
352{
353 hwuine_();
354}
355
356void THerwig6::Hwepro()
357
358{
359 hwepro_();
360}
361
362void THerwig6::Hwbgen()
363
364{
365 hwbgen_();
366}
367
368void THerwig6::Hwdhob()
369
370{
371 hwdhob_();
372}
373
374void THerwig6::Hwcfor()
375
376{
377 hwcfor_();
378}
379
380void THerwig6::Hwcdec()
381
382{
383 hwcdec_();
384}
385
386void THerwig6::Hwdhad()
387
388{
389 hwdhad_();
390}
391
392void THerwig6::Hwdhvy()
393
394{
395 hwdhvy_();
396}
397
398void THerwig6::Hwmevt()
399
400{
401 hwmevt_();
402}
403
404void THerwig6::Hwufne()
405
406{
407 hwufne_();
408}
409
410void THerwig6::Hwefin()
411
412{
413 hwefin_();
414}
415
416void 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