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 }