Updated instructions according to the suggestions/findings from Federico
[u/mrichter/AliRoot.git] / TFluka / TFluka.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 Revision 1.4  2002/11/04 16:00:46  iglez2
19 The conversion between ID and PDG now uses Fluka routines and arrays which is more consistent.
20
21 Revision 1.3  2002/10/22 15:12:14  alibrary
22 Introducing Riostream.h
23
24 Revision 1.2  2002/10/14 14:57:40  hristov
25 Merging the VirtualMC branch to the main development branch (HEAD)
26
27 Revision 1.1.2.8  2002/10/08 16:33:17  iglez2
28 LSOUIT is set to true before the second call to flukam.
29
30 Revision 1.1.2.7  2002/10/08 09:30:37  iglez2
31 Solved stupid missing ;
32
33 Revision 1.1.2.6  2002/10/07 13:40:22  iglez2
34 First implementations of the PDG <--> Fluka Id conversion routines
35
36 Revision 1.1.2.5  2002/09/26 16:26:03  iglez2
37 Added verbosity
38 Call to gAlice->Generator()->Generate()
39
40 Revision 1.1.2.4  2002/09/26 13:22:23  iglez2
41 Naive implementation of ProcessRun and ProcessEvent
42 Opening/Closing of input file (fInputFileName) with FORTRAN unit 5 before/after the first call to flukam inside Init()
43
44 Revision 1.1.2.3  2002/09/20 15:35:51  iglez2
45 Modification of LFDRTR. Value is passed to FLUKA !!!
46
47 Revision 1.1.2.2  2002/09/18 14:34:44  iglez2
48 Revised version with all pure virtual methods implemented
49
50 Revision 1.1.2.1  2002/07/24 08:49:41  alibrary
51 Adding TFluka to VirtualMC
52
53 Revision 1.1  2002/07/05 13:10:07  morsch
54 First commit of Fluka interface.
55
56 */
57
58 #include <Riostream.h>
59
60 #include "TFluka.h"
61 #include "TCallf77.h"      //For the fortran calls
62 #include "Fdblprc.h"       //(DBLPRC) fluka common
63 #include "Fiounit.h"       //(IOUNIT) fluka common
64 #include "Fepisor.h"       //(EPISOR) fluka common
65 #include "Fpart.h"         //(PART)   fluka common
66 #include "TVirtualMC.h"
67
68 #include "TG4GeometryManager.h" //For the geometry management
69 #include "TG4DetConstruction.h" //For the detector construction
70
71 #include "FGeometryInit.hh"
72
73 // Fluka methods that may be needed.
74 #ifndef WIN32 
75 # define flukam  flukam_
76 # define fluka_openinp fluka_openinp_
77 # define fluka_closeinp fluka_closeinp_
78 # define mcihad mcihad_
79 # define mpdgha mpdgha_
80 #else 
81 # define flukam  FLUKAM
82 # define fluka_openinp FLUKA_OPENINP
83 # define fluka_closeinp FLUKA_CLOSEINP
84 # define mcihad MCIHAD
85 # define mpdgha MPDGHA
86 #endif
87
88 extern "C" 
89 {
90   //
91   // Prototypes for FLUKA functions
92   //
93   void type_of_call flukam(const int&);
94   void type_of_call fluka_openinp(const int&, DEFCHARA);
95   void type_of_call fluka_closeinp(const int&);
96   int  type_of_call mcihad(const int&);
97   int  type_of_call mpdgha(const int&);
98 }
99
100 //
101 // Class implementation for ROOT
102 //
103 ClassImp(TFluka)
104
105 //
106 //----------------------------------------------------------------------------
107 // TFluka constructors and destructors.
108 //____________________________________________________________________________ 
109 TFluka::TFluka()
110   :TVirtualMC(),
111    fVerbosityLevel(0),
112    fInputFileName(""),
113    fDetector(0)
114
115   //
116   // Default constructor
117   //
118
119  
120 TFluka::TFluka(const char *title, Int_t verbosity)
121   :TVirtualMC("TFluka",title),
122    fVerbosityLevel(verbosity),
123    fInputFileName(""),
124    fDetector(0)
125 {
126   if (fVerbosityLevel >=3)
127     cout << "==> TFluka::TFluka(" << title << ") constructor called." << endl;
128
129   
130   // create geometry manager
131   if (fVerbosityLevel >=2)
132     cout << "\t* Creating G4 Geometry manager..." << endl;
133   fGeometryManager = new TG4GeometryManager();
134   if (fVerbosityLevel >=2)
135     cout << "\t* Creating G4 Detector..." << endl;
136   fDetector = new TG4DetConstruction();
137   FGeometryInit* geominit = FGeometryInit::GetInstance();
138   if (geominit)
139     geominit->setDetConstruction(fDetector);
140   else {
141     cerr << "ERROR: Could not create FGeometryInit!" << endl;
142     cerr << "       Exiting!!!" << endl;
143     abort();
144   }
145
146   if (fVerbosityLevel >=3)
147     cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
148 }
149
150 TFluka::~TFluka() {
151   if (fVerbosityLevel >=3)
152     cout << "==> TFluka::~TFluka() destructor called." << endl;
153
154   delete fGeometryManager;
155
156   if (fVerbosityLevel >=3)
157     cout << "<== TFluka::~TFluka() destructor called." << endl;
158 }
159
160 //
161 //_____________________________________________________________________________
162 // TFluka control methods
163 //____________________________________________________________________________ 
164 void TFluka::Init() {
165   if (fVerbosityLevel >=3)
166     cout << "==> TFluka::Init() called." << endl;
167
168   if (fVerbosityLevel >=2)
169     cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
170          << ") in fluka..." << endl;
171   GLOBAL.lfdrtr = true;
172
173   if (fVerbosityLevel >=2)
174     cout << "\t* Opening file " << fInputFileName << endl;
175   const char* fname = fInputFileName;
176   fluka_openinp(lunin, PASSCHARA(fname));
177
178   if (fVerbosityLevel >=2)
179     cout << "\t* Calling flukam..." << endl;
180   flukam(1);
181
182   if (fVerbosityLevel >=2)
183     cout << "\t* Closing file " << fInputFileName << endl;
184   fluka_closeinp(lunin);
185
186   if (fVerbosityLevel >=3)
187     cout << "<== TFluka::Init() called." << endl;
188 }
189
190 void TFluka::FinishGeometry() {
191   if (fVerbosityLevel >=3)
192     cout << "==> TFluka::FinishGeometry() called." << endl;
193
194   fGeometryManager->Ggclos();
195
196   if (fVerbosityLevel >=3)
197     cout << "<== TFluka::FinishGeometry() called." << endl;
198
199
200 void TFluka::BuildPhysics() {
201   if (fVerbosityLevel >=3)
202     cout << "==> TFluka::BuildPhysics() called." << endl;
203
204
205   if (fVerbosityLevel >=3)
206     cout << "<== TFluka::BuildPhysics() called." << endl;
207 }  
208
209 void TFluka::ProcessEvent() {
210   if (fVerbosityLevel >=3)
211     cout << "==> TFluka::ProcessEvent() called." << endl;
212
213   if (fVerbosityLevel >=3)
214     cout << "<== TFluka::ProcessEvent() called." << endl;
215 }
216
217
218 void TFluka::ProcessRun(Int_t nevent) {
219   if (fVerbosityLevel >=3)
220     cout << "==> TFluka::ProcessRun(" << nevent << ") called." 
221          << endl;
222
223   if (fVerbosityLevel >=2) {
224     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
225     cout << "\t* Calling flukam again..." << endl;
226   }
227   fApplication->GeneratePrimaries();
228   EPISOR.lsouit = true;
229   flukam(1);
230
231   if (fVerbosityLevel >=3)
232     cout << "<== TFluka::ProcessRun(" << nevent << ") called." 
233          << endl;
234 }
235
236 //_____________________________________________________________________________
237 // methods for building/management of geometry
238 //____________________________________________________________________________ 
239 // functions from GCONS 
240 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,  
241                     Float_t &dens, Float_t &radl, Float_t &absl,
242                     Float_t* ubuf, Int_t& nbuf) {
243 //
244   fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
245
246
247 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,  
248                     Double_t &dens, Double_t &radl, Double_t &absl,
249                     Double_t* ubuf, Int_t& nbuf) {
250 //
251   fGeometryManager->Gfmate(imat, name, a, z, dens, radl, absl, ubuf, nbuf);
252
253
254 // detector composition
255 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
256                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
257                       Float_t* buf, Int_t nwbuf) {
258 //
259   fGeometryManager
260     ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf); 
261
262 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
263                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
264                       Double_t* buf, Int_t nwbuf) {
265 //
266   fGeometryManager
267     ->Material(kmat, name, a, z, dens, radl, absl, buf, nwbuf); 
268
269
270 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a, 
271                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
272 //
273    fGeometryManager
274      ->Mixture(kmat, name, a, z, dens, nlmat, wmat); 
275
276 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a, 
277                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
278 //
279    fGeometryManager
280      ->Mixture(kmat, name, a, z, dens, nlmat, wmat); 
281
282
283 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
284                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
285                     Double_t stemax, Double_t deemax, Double_t epsil, 
286                     Double_t stmin, Float_t* ubuf, Int_t nbuf) { 
287   //
288   fGeometryManager
289     ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
290              epsil, stmin, ubuf, nbuf);
291
292 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
293                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
294                     Double_t stemax, Double_t deemax, Double_t epsil, 
295                     Double_t stmin, Double_t* ubuf, Int_t nbuf) { 
296   //
297   fGeometryManager
298     ->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
299              epsil, stmin, ubuf, nbuf);
300
301
302 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX, 
303                     Double_t thetaY, Double_t phiY, Double_t thetaZ, 
304                     Double_t phiZ) {
305 //                   
306   fGeometryManager
307     ->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ); 
308
309
310 void TFluka::Gstpar(Int_t itmed, const char *param, Double_t parval) {
311 //
312   fGeometryManager->Gstpar(itmed, param, parval); 
313 }    
314
315 // functions from GGEOM 
316 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
317                      Float_t *upar, Int_t np)  {
318 //
319   return fGeometryManager->Gsvolu(name, shape, nmed, upar, np); 
320 }
321 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
322                      Double_t *upar, Int_t np)  {
323 //
324   return fGeometryManager->Gsvolu(name, shape, nmed, upar, np); 
325 }
326  
327 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv, 
328                    Int_t iaxis) {
329 //
330   fGeometryManager->Gsdvn(name, mother, ndiv, iaxis); 
331
332
333 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv, 
334                     Int_t iaxis, Double_t c0i, Int_t numed) {
335 //
336   fGeometryManager->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed); 
337
338
339 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step, 
340                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
341 //                      
342   fGeometryManager->Gsdvt(name, mother, step, iaxis, numed, ndvmx); 
343
344
345 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step, 
346                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) { 
347 //
348   fGeometryManager->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx); 
349
350
351 void TFluka::Gsord(const char *name, Int_t iax) {
352 //
353   fGeometryManager->Gsord(name, iax); 
354
355
356 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,  
357                    Double_t x, Double_t y, Double_t z, Int_t irot, 
358                    const char *konly) {
359 //
360   fGeometryManager->Gspos(name, nr, mother, x, y, z, irot, konly); 
361
362
363 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
364                     Double_t x, Double_t y, Double_t z, Int_t irot,
365                     const char *konly, Float_t *upar, Int_t np)  {
366   //
367   fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
368
369 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
370                     Double_t x, Double_t y, Double_t z, Int_t irot,
371                     const char *konly, Double_t *upar, Int_t np)  {
372   //
373   fGeometryManager->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
374
375
376 void TFluka::Gsbool(const char* onlyVolName, const char* manyVolName) {
377 //
378   fGeometryManager->Gsbool(onlyVolName, manyVolName);
379 }
380
381 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t *ppckov,
382                          Float_t *absco, Float_t *effic, Float_t *rindex) {
383 //
384   fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
385 }  
386 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
387                          Double_t *absco, Double_t *effic, Double_t *rindex) {
388 //
389   fGeometryManager->SetCerenkov(itmed, npckov, ppckov, absco, effic, rindex);
390 }  
391     
392 // Euclid
393 void TFluka::WriteEuclid(const char* fileName, const char* topVol, 
394                           Int_t number, Int_t nlevel) {
395 //
396   fGeometryManager->WriteEuclid(fileName, topVol, number, nlevel); 
397
398
399
400
401
402
403
404 //____________________________________________________________________________ 
405 // ID <--> PDG transformations
406 //_____________________________________________________________________________
407 Int_t TFluka::IdFromPDG(Int_t pdg) const 
408 {
409   //
410   // Return Fluka code from PDG and pseudo ENDF code
411
412   // MCIHAD() goes from pdg to fluka internal.
413   Int_t intfluka = mcihad(pdg);
414   // KPTOIP array goes from internal to official
415   return GetFlukaKPTOIP(intfluka);
416 }
417
418 Int_t TFluka::PDGFromId(Int_t id) const 
419 {
420   //
421   // Return PDG code and pseudo ENDF code from Fluka code
422
423   //IPTOKP array goes from official to internal
424   Int_t intfluka = GetFlukaIPTOKP(id);
425   //MPKDHA() goes from internal to PDG
426   return mpdgha(intfluka);
427   
428 }