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