Updates to HERWIG 6.507 and use with MC@NLO (R. Guernane)
[u/mrichter/AliRoot.git] / THerwig / THerwig6.cxx
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
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
51
52 #include "THerwig6.h"
53 #include "TClonesArray.h"
54 #include "TParticle.h"
55 #include "TObjArray.h"
56
57
58 ClassImp(THerwig6)
59
60 extern "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
67 THerwig6::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);
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);
120  }
121
122 //------------------------------------------------------------------------------
123  THerwig6::~THerwig6()
124  {
125    // Destructor. The data members of TGenerator are delete by itself
126  }
127
128 //______________________________________________________________________________
129 void 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 //______________________________________________________________________________
155 void THerwig6::OpenFortranFile(int lun, char* name) {
156   herwig6_open_fortran_file_(&lun, name, strlen(name));
157 }
158
159 //______________________________________________________________________________
160 void THerwig6::CloseFortranFile(int lun) {
161   herwig6_close_fortran_file_(&lun);
162 }
163
164 void 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
247 void THerwig6::PrepareRun()
248 {
249   // compute parameter dependent constants
250   hwuinc_();
251   // initialize elementary processes
252   hweini_();
253 }
254
255 //______________________________________________________________________________
256 TObjArray* 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 //______________________________________________________________________________
319 Int_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;
333   TClonesArray &refParticles = *particles;
334   refParticles.Clear();
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 //
342         new(refParticles[i]) TParticle(
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++) {
363       new(refParticles[i]) TParticle(
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
384 void THerwig6::Hwigin()
385 {
386   hwigin_();
387 }
388
389 void THerwig6::Hwuinc()
390 {
391   hwuinc_();
392 }
393
394 void THerwig6::Hwusta(char* name)
395
396 {
397   hwusta_(name,8);
398 }
399
400 void THerwig6::Hweini()
401
402 {
403   hweini_();
404 }
405
406 void THerwig6::Hwuine()
407
408 {
409   hwuine_();
410 }
411
412 void THerwig6::Hwepro()
413
414 {
415   hwepro_();
416 }
417
418 void THerwig6::Hwbgen()
419
420 {
421   hwbgen_();
422 }
423
424 void THerwig6::Hwdhob()
425
426 {
427   hwdhob_();
428 }
429
430 void THerwig6::Hwcfor()
431
432 {
433   hwcfor_();
434 }
435
436 void THerwig6::Hwcdec()
437
438 {
439   hwcdec_();
440 }
441
442 void THerwig6::Hwdhad()
443
444 {
445   hwdhad_();
446 }
447
448 void THerwig6::Hwdhvy()
449
450 {
451   hwdhvy_();
452 }
453
454 void THerwig6::Hwmevt()
455
456 {
457   hwmevt_();
458 }
459
460 void THerwig6::Hwufne()
461
462 {
463   hwufne_();
464 }
465
466 void THerwig6::Hwefin()
467
468 {
469   hwefin_();
470 }
471
472 void THerwig6::Hwiodk(int iopt)
473
474 {
475   hwiodk_(iopt);
476 }
477
478 void 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();
490   int nEvToGenerate=1;
491   for (int i=0;i<nEvToGenerate;i++)
492     {
493       GenerateEvent();
494       // do your stuff. For ex:
495       int nOfPar=GetNumberOfParticles(); // from TGenerator
496       for (int j=0; j<nOfPar; j++)
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