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