Interface to HERWIG
[u/mrichter/AliRoot.git] / THerwig / THerwig6.cxx
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