Coding conventions
[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  }
104
105 //------------------------------------------------------------------------------
106  THerwig6::~THerwig6()
107  {
108 // Destroys the object, deletes and disposes all TParticles currently on list.
109
110     if (fParticles) {
111       fParticles->Delete();
112       delete fParticles;
113       fParticles = 0;
114    }
115  }
116
117 //______________________________________________________________________________
118 void THerwig6::GenerateEvent() 
119 {
120   
121   //initialize event
122   hwuine_();
123   // generate hard subprocess
124   hwepro_();
125   // generate parton cascades
126   hwbgen_();
127   // do heavy objects decay
128   hwdhob_();
129   // do cluster formation
130   hwcfor_();
131   // do cluster decays
132   hwcdec_();
133   // do unstable particle decays
134   hwdhad_();
135   // do heavy flavor hadrons decay
136   hwdhvy_();
137   // add soft underlying event
138   hwmevt_();
139   // finish event
140   hwufne_();
141 }
142
143 //______________________________________________________________________________
144 void THerwig6::OpenFortranFile(int lun, char* name) {
145   herwig6_open_fortran_file_(&lun, name, strlen(name));
146 }
147
148 //______________________________________________________________________________
149 void THerwig6::CloseFortranFile(int lun) {
150   herwig6_close_fortran_file_(&lun);
151 }
152
153 void THerwig6::Initialize(const char *beam, const char *target, double pbeam1, double pbeam2, int iproc)
154
155 {
156   // perform the initialization for Herwig6
157   // sets correct title. 
158   // after calling this method all parameters are set to their default
159   // values. If you want to modify any parameter you have to set the new
160   // value after calling Initialize and before PrepareRun.
161
162    char  cbeam[8];
163    strncpy(cbeam,beam,8);
164    char  ctarget[8];
165    strncpy(ctarget,target,8);
166    printf("\n Initializing Herwig !! \n");
167    if ( (!strncmp(beam, "E+"    ,2)) &&
168         (!strncmp(beam, "E-"    ,2)) &&
169         (!strncmp(beam, "MU+"   ,3)) &&
170         (!strncmp(beam, "MU-"   ,3)) &&
171         (!strncmp(beam, "NUE"   ,3)) &&
172         (!strncmp(beam, "NUEB"  ,4)) &&
173         (!strncmp(beam, "NUMU"  ,4)) &&
174         (!strncmp(beam, "NMUB"  ,4)) &&
175         (!strncmp(beam, "NTAU"  ,4)) &&
176         (!strncmp(beam, "NTAB"  ,4)) &&
177         (!strncmp(beam, "GAMA"  ,4)) &&
178         (!strncmp(beam, "P       ",8)) &&
179         (!strncmp(beam, "PBAR    ",8)) &&
180         (!strncmp(beam, "N"     ,1)) &&
181         (!strncmp(beam, "NBAR"  ,4)) &&
182         (!strncmp(beam, "PI+"   ,3)) &&
183         (!strncmp(beam, "PI-"   ,3)) ) {
184       printf("WARNING! In THerwig6:Initialize():\n");
185       printf(" specified beam=%s is unrecognized .\n",beam);
186       printf(" resetting to \"P\" .");
187       sprintf(cbeam,"P");
188    }
189
190    if ( (!strncmp(target, "E+"    ,2)) &&
191         (!strncmp(target, "E-"    ,2)) &&
192         (!strncmp(target, "MU+"   ,3)) &&
193         (!strncmp(target, "MU-"   ,3)) &&
194         (!strncmp(target, "NUE"   ,3)) &&
195         (!strncmp(target, "NUEB"  ,4)) &&
196         (!strncmp(target, "NUMU"  ,4)) &&
197         (!strncmp(target, "NMUB"  ,4)) &&
198         (!strncmp(target, "NTAU"  ,4)) &&
199         (!strncmp(target, "NTAB"  ,4)) &&
200         (!strncmp(target, "GAMA"  ,4)) &&
201         (!strncmp(target, "P       ",8)) &&
202         (!strncmp(target, "PBAR    ",8)) &&
203         (!strncmp(target, "N"     ,1)) &&
204         (!strncmp(target, "NBAR"  ,4)) &&
205         (!strncmp(target, "PI+"   ,3)) &&
206         (!strncmp(target, "PI-"   ,3)) ) {
207       printf("WARNING! In THerwig6:Initialize():\n");
208       printf(" specified target=%s is unrecognized .\n",target);
209       printf(" resetting to \"P\" .");
210       sprintf(ctarget,"P");
211    }
212
213    // initialization:
214    // type of beams
215    strncpy(fHwbmch->PART1,beam,8);
216    strncpy(fHwbmch->PART2,target,8);
217    // momentum of beams
218    fHwproc->PBEAM1=pbeam1;
219    fHwproc->PBEAM2=pbeam2;
220    // process to generate
221    fHwproc->IPROC=iproc;
222    // not used in the class definition 
223    fHwproc->MAXEV=1;
224
225    // reset all parameters
226    hwigin_();
227    
228    // set correct title
229    char atitle[132];
230    double win=pbeam1+pbeam2;
231    printf("\n %s - %s at %g GeV",beam,target,win);
232    sprintf(atitle,"%s-%s at %g GeV",cbeam,ctarget,win);
233    SetTitle(atitle); 
234 }
235
236 void THerwig6::PrepareRun()
237 {
238   // compute parameter dependent constants
239   hwuinc_();
240   // initialize elementary processes
241   hweini_();
242 }
243
244 //______________________________________________________________________________
245 TObjArray* THerwig6::ImportParticles(Option_t *option)
246 {
247 //
248 //  Default primary creation method. It reads the /HEPEVT/ common block which
249 //  has been filled by the GenerateEvent method. If the event generator does
250 //  not use the HEPEVT common block, This routine has to be overloaded by
251 //  the subclasses.
252 //  The default action is to store only the stable particles (ISTHEP = 1)
253 //  This can be demanded explicitly by setting the option = "Final"
254 //  If the option = "All", all the particles are stored.
255 //
256   fParticles->Clear();
257   Int_t numpart = fHepevt->NHEP;
258   TClonesArray &a = *((TClonesArray*)fParticles);
259   if (!strcmp(option,"") || !strcmp(option,"Final")) {
260     for (Int_t i = 0; i<=numpart; i++) {
261       if (fHepevt->ISTHEP[i] == 1) {
262 //
263 //  Use the common block values for the TParticle constructor
264 //
265         new(a[i]) TParticle(
266                                    fHepevt->IDHEP[i],
267                                    fHepevt->ISTHEP[i],
268                                    fHepevt->JMOHEP[i][0]-1,
269                                    fHepevt->JMOHEP[i][1]-1,
270                                    fHepevt->JDAHEP[i][0]-1,
271                                    fHepevt->JDAHEP[i][1]-1,
272
273                                    fHepevt->PHEP[i][0],
274                                    fHepevt->PHEP[i][1],
275                                    fHepevt->PHEP[i][2],
276                                    fHepevt->PHEP[i][3],
277                                    fHepevt->VHEP[i][0],
278                                    fHepevt->VHEP[i][1],
279                                    fHepevt->VHEP[i][2],
280                                    fHepevt->VHEP[i][3]);
281         }
282      }
283   }
284   else if (!strcmp(option,"All")) {
285     for (Int_t i = 0; i<=numpart; i++) {
286       new(a[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   return fParticles;
305 }
306
307 //______________________________________________________________________________
308 Int_t THerwig6::ImportParticles(TClonesArray *particles, Option_t *option)
309 {
310 //
311 //  Default primary creation method. It reads the /HEPEVT/ common block which
312 //  has been filled by the GenerateEvent method. If the event generator does
313 //  not use the HEPEVT common block, This routine has to be overloaded by
314 //  the subclasses.
315 //  The function loops on the generated particles and store them in
316 //  the TClonesArray pointed by the argument particles.
317 //  The default action is to store only the stable particles (ISTHEP = 1)
318 //  This can be demanded explicitly by setting the option = "Final"
319 //  If the option = "All", all the particles are stored.
320 //
321   if (particles == 0) return 0;
322   TClonesArray &refParticles = *particles;
323   refParticles.Clear();
324   Int_t numpart = fHepevt->NHEP;
325   if (!strcmp(option,"") || !strcmp(option,"Final")) {
326     for (Int_t i = 0; i< numpart; i++) {
327       if (fHepevt->ISTHEP[i] == 1) {
328 //
329 //  Use the common block values for the TParticle constructor
330 //
331         new(refParticles[i]) TParticle(
332                                    fHepevt->IDHEP[i],
333                                    fHepevt->ISTHEP[i],
334                                    fHepevt->JMOHEP[i][0]-1,
335                                    fHepevt->JMOHEP[i][1]-1,
336                                    fHepevt->JDAHEP[i][0]-1,
337                                    fHepevt->JDAHEP[i][1]-1,
338
339                                    fHepevt->PHEP[i][0],
340                                    fHepevt->PHEP[i][1],
341                                    fHepevt->PHEP[i][2],
342                                    fHepevt->PHEP[i][3],
343                                    fHepevt->VHEP[i][0],
344                                    fHepevt->VHEP[i][1],
345                                    fHepevt->VHEP[i][2],
346                                    fHepevt->VHEP[i][3]);
347         }
348      }
349   }
350   else if (!strcmp(option,"All")) {
351     for (Int_t i = 0; i< numpart; i++) {
352       new(refParticles[i]) TParticle(
353                                    fHepevt->IDHEP[i],
354                                    fHepevt->ISTHEP[i],
355                                    fHepevt->JMOHEP[i][0]-1,
356                                    fHepevt->JMOHEP[i][1]-1,
357                                    fHepevt->JDAHEP[i][0]-1,
358                                    fHepevt->JDAHEP[i][1]-1,
359
360                                    fHepevt->PHEP[i][0],
361                                    fHepevt->PHEP[i][1],
362                                    fHepevt->PHEP[i][2],
363                                    fHepevt->PHEP[i][3],
364                                    fHepevt->VHEP[i][0],
365                                    fHepevt->VHEP[i][1],
366                                    fHepevt->VHEP[i][2],
367                                    fHepevt->VHEP[i][3]);
368     }
369   }
370   return numpart;
371 }
372
373 void THerwig6::Hwigin()
374 {
375   hwigin_();
376 }
377
378 void THerwig6::Hwuinc()
379 {
380   hwuinc_();
381 }
382
383 void THerwig6::Hwusta(char* name)
384
385 {
386   hwusta_(name,8);
387 }
388
389 void THerwig6::Hweini()
390
391 {
392   hweini_();
393 }
394
395 void THerwig6::Hwuine()
396
397 {
398   hwuine_();
399 }
400
401 void THerwig6::Hwepro()
402
403 {
404   hwepro_();
405 }
406
407 void THerwig6::Hwbgen()
408
409 {
410   hwbgen_();
411 }
412
413 void THerwig6::Hwdhob()
414
415 {
416   hwdhob_();
417 }
418
419 void THerwig6::Hwcfor()
420
421 {
422   hwcfor_();
423 }
424
425 void THerwig6::Hwcdec()
426
427 {
428   hwcdec_();
429 }
430
431 void THerwig6::Hwdhad()
432
433 {
434   hwdhad_();
435 }
436
437 void THerwig6::Hwdhvy()
438
439 {
440   hwdhvy_();
441 }
442
443 void THerwig6::Hwmevt()
444
445 {
446   hwmevt_();
447 }
448
449 void THerwig6::Hwufne()
450
451 {
452   hwufne_();
453 }
454
455 void THerwig6::Hwefin()
456
457 {
458   hwefin_();
459 }
460
461 void THerwig6::SetupTest()
462 {
463   // exampe of running herwig and generating one event
464   // after changing some options
465   Initialize("P","PBAR",900.,900.,1500);
466   // here you can set some parameters
467   SetPTMIN(15.); // Min pt in hadronic jet production
468   SetYJMIN(-4.); // Min jet rapidity
469   SetYJMAX(4.);  // Max jet rapidity
470   // after you set your wished parameters
471   // herwig can do its work
472   PrepareRun();
473   int nEvToGenerate=1;
474   for (int i=0;i<nEvToGenerate;i++)
475     {
476       GenerateEvent();
477       // do your stuff. For ex:
478       int nOfPar=GetNumberOfParticles(); // from TGenerator
479       for (int j=0; j<nOfPar; j++)
480         {
481           TParticle* p=GetParticle(j);
482           // here you do whatever you want with the particle
483           p->Print();
484         };
485     };
486 }
487
488
489
490
491