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 /* $Id$ */
  17
  18 //
  19 // Realisation of the TVirtualMC interface for the FLUKA code
  20 // (See official web side http://www.fluka.org/).
  21 //
  22 // This implementation makes use of the TGeo geometry modeller.
  23 // User configuration is via automatic generation of FLUKA input cards.
  24 //
  25 // Authors:
  26 // A. Fasso
  27 // E. Futo
  28 // A. Gheata
  29 // A. Morsch
  30 //
  31
  32 #include <Riostream.h>
  33
  34 #include "TFluka.h"
  35 #include "TCallf77.h"      //For the fortran calls
  36 #include "Fdblprc.h"       //(DBLPRC) fluka common
  37 #include "Fepisor.h"       //(EPISOR) fluka common
  38 #include "Ffinuc.h"        //(FINUC)  fluka common
  39 #include "Fiounit.h"       //(IOUNIT) fluka common
  40 #include "Fpaprop.h"       //(PAPROP) fluka common
  41 #include "Fpart.h"         //(PART)   fluka common
  42 #include "Ftrackr.h"       //(TRACKR) fluka common
  43 #include "Fpaprop.h"       //(PAPROP) fluka common
  44 #include "Ffheavy.h"       //(FHEAVY) fluka common
  45 #include "Fopphst.h"       //(OPPHST) fluka common
  46 #include "Fstack.h"        //(STACK)  fluka common
  47 #include "Fstepsz.h"       //(STEPSZ) fluka common
  48 #include "Fopphst.h"       //(OPPHST) fluka common
  49
  50 #include "TVirtualMC.h"
  51 #include "TMCProcess.h"
  52 #include "TGeoManager.h"
  53 #include "TGeoMaterial.h"
  54 #include "TGeoMedium.h"
  55 #include "TFlukaMCGeometry.h"
  56 #include "TGeoMCGeometry.h"
  57 #include "TFlukaCerenkov.h"
  58 #include "TFlukaConfigOption.h"
  59 #include "TFlukaScoringOption.h"
  60 #include "TLorentzVector.h"
  61 #include "TArrayI.h"
  62
  63 // Fluka methods that may be needed.
  64 #ifndef WIN32
  65 # define flukam  flukam_
  66 # define fluka_openinp fluka_openinp_
  67 # define fluka_closeinp fluka_closeinp_
  68 # define mcihad mcihad_
  69 # define mpdgha mpdgha_
  70 # define newplo newplo_
  71 #else
  72 # define flukam  FLUKAM
  73 # define fluka_openinp FLUKA_OPENINP
  74 # define fluka_closeinp FLUKA_CLOSEINP
  75 # define mcihad MCIHAD
  76 # define mpdgha MPDGHA
  77 # define newplo NEWPLO
  78 #endif
  79
  80 extern "C"
  81 {
  82   //
  83   // Prototypes for FLUKA functions
  84   //
  85   void type_of_call flukam(const int&);
  86   void type_of_call newplo();
  87   void type_of_call fluka_openinp(const int&, DEFCHARA);
  88   void type_of_call fluka_closeinp(const int&);
  89   int  type_of_call mcihad(const int&);
  90   int  type_of_call mpdgha(const int&);
  91 }
  92
  93 //
  94 // Class implementation for ROOT
  95 //
  96 ClassImp(TFluka)
  97
  98 //
  99 //----------------------------------------------------------------------------
 100 // TFluka constructors and destructors.
 101 //______________________________________________________________________________
 102 TFluka::TFluka()
 103   :TVirtualMC(),
 104    fVerbosityLevel(0),
 105    fInputFileName(""),
 106    fUserConfig(0),
 107    fUserScore(0)
 108 {
 109   //
 110   // Default constructor
 111   //
 112    fGeneratePemf = kFALSE;
 113    fNVolumes = 0;
 114    fCurrentFlukaRegion = -1;
 115    fGeom = 0;
 116    fMCGeo = 0;
 117    fMaterials = 0;
 118    fDummyBoundary = 0;
 119    fFieldFlag = 1;
 120    fStopped   = 0;
 121    fStopEvent = 0;
 122    fStopRun   = 0;
 123    fNEvent    = 0;
 124 }
 125
 126 //______________________________________________________________________________
 127 TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
 128   :TVirtualMC("TFluka",title, isRootGeometrySupported),
 129    fVerbosityLevel(verbosity),
 130    fInputFileName(""),
 131    fTrackIsEntering(0),
 132    fTrackIsExiting(0),
 133    fTrackIsNew(0),
 134    fUserConfig(new TObjArray(100)),
 135    fUserScore(new TObjArray(100))
 136 {
 137   // create geometry interface
 138    if (fVerbosityLevel >=3)
 139        cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
 140    SetCoreInputFileName();
 141    SetInputFileName();
 142    SetGeneratePemf(kFALSE);
 143    fNVolumes      = 0;
 144    fCurrentFlukaRegion = -1;
 145    fDummyBoundary = 0;
 146    fFieldFlag = 1;
 147    fGeneratePemf = kFALSE;
 148    fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kTRUE);
 149    fGeom  = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
 150    if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
 151    fMaterials = 0;
 152    fStopped   = 0;
 153    fStopEvent = 0;
 154    fStopRun   = 0;
 155    fNEvent    = 0;
 156 }
 157
 158 //______________________________________________________________________________
 159 TFluka::~TFluka() {
 160 // Destructor
 161     if (fVerbosityLevel >=3)
 162         cout << "<== TFluka::~TFluka() destructor called." << endl;
 163
 164     delete fGeom;
 165     delete fMCGeo;
 166
 167     if (fUserConfig) {
 168         fUserConfig->Delete();
 169         delete fUserConfig;
 170     }
 171
 172
 173 }
 174
 175 //
 176 //______________________________________________________________________________
 177 // TFluka control methods
 178 //______________________________________________________________________________
 179 void TFluka::Init() {
 180 //
 181 //  Geometry initialisation
 182 //
 183     if (fVerbosityLevel >=3) cout << "==> TFluka::Init() called." << endl;
 184
 185     if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
 186     fApplication->ConstructGeometry();
 187     TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
 188     gGeoManager->SetTopVolume(top);
 189     gGeoManager->CloseGeometry("di");
 190     gGeoManager->DefaultColors();  // to be removed
 191     fNVolumes = fGeom->NofVolumes();
 192     fGeom->CreateFlukaMatFile("flukaMat.inp");
 193     if (fVerbosityLevel >=3) {
 194        printf("== Number of volumes: %i\n ==", fNVolumes);
 195        cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
 196     }
 197     // now we have TGeo geometry created and we have to patch FlukaVmc.inp
 198     // with the material mapping file FlukaMat.inp
 199 }
 200
 201
 202 //______________________________________________________________________________
 203 void TFluka::FinishGeometry() {
 204 //
 205 // Build-up table with region to medium correspondance
 206 //
 207   if (fVerbosityLevel >=3) {
 208     cout << "==> TFluka::FinishGeometry() called." << endl;
 209     printf("----FinishGeometry - nothing to do with TGeo\n");
 210     cout << "<== TFluka::FinishGeometry() called." << endl;
 211   }
 212 }
 213
 214 //______________________________________________________________________________
 215 void TFluka::BuildPhysics() {
 216 //
 217 //  Prepare FLUKA input files and call FLUKA physics initialisation
 218 //
 219
 220     if (fVerbosityLevel >=3)
 221         cout << "==> TFluka::BuildPhysics() called." << endl;
 222 // Prepare input file with the current physics settings
 223     InitPhysics();
 224     cout << "\t* InitPhysics() - Prepare input file was called" << endl;
 225
 226     if (fVerbosityLevel >=2)
 227         cout << "\t* Changing lfdrtr = (" << (GLOBAL.lfdrtr?'T':'F')
 228              << ") in fluka..." << endl;
 229     GLOBAL.lfdrtr = true;
 230
 231     if (fVerbosityLevel >=2)
 232         cout << "\t* Opening file " << fInputFileName << endl;
 233     const char* fname = fInputFileName;
 234     fluka_openinp(lunin, PASSCHARA(fname));
 235
 236     if (fVerbosityLevel >=2)
 237         cout << "\t* Calling flukam..." << endl;
 238     flukam(1);
 239
 240     if (fVerbosityLevel >=2)
 241         cout << "\t* Closing file " << fInputFileName << endl;
 242     fluka_closeinp(lunin);
 243
 244     FinishGeometry();
 245
 246     if (fVerbosityLevel >=3)
 247         cout << "<== TFluka::Init() called." << endl;
 248
 249
 250     if (fVerbosityLevel >=3)
 251         cout << "<== TFluka::BuildPhysics() called." << endl;
 252 }
 253
 254 //______________________________________________________________________________
 255 void TFluka::ProcessEvent() {
 256 //
 257 // Process one event
 258 //
 259     if (fStopRun) {
 260         printf("User Run Abortion: No more events handled !\n");
 261         fNEvent += 1;
 262         return;
 263     }
 264
 265     if (fVerbosityLevel >=3)
 266         cout << "==> TFluka::ProcessEvent() called." << endl;
 267     fApplication->GeneratePrimaries();
 268     EPISOR.lsouit = true;
 269     flukam(1);
 270     if (fVerbosityLevel >=3)
 271         cout << "<== TFluka::ProcessEvent() called." << endl;
 272     //
 273     // Increase event number
 274     //
 275     fNEvent += 1;
 276 }
 277
 278 //______________________________________________________________________________
 279 Bool_t TFluka::ProcessRun(Int_t nevent) {
 280 //
 281 // Run steering
 282 //
 283
 284   if (fVerbosityLevel >=3)
 285     cout << "==> TFluka::ProcessRun(" << nevent << ") called."
 286          << endl;
 287
 288   if (fVerbosityLevel >=2) {
 289     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
 290     cout << "\t* Calling flukam again..." << endl;
 291   }
 292
 293   fApplication->InitGeometry();
 294   Int_t todo = TMath::Abs(nevent);
 295   for (Int_t ev = 0; ev < todo; ev++) {
 296       fApplication->BeginEvent();
 297       ProcessEvent();
 298       fApplication->FinishEvent();
 299   }
 300
 301   if (fVerbosityLevel >=3)
 302     cout << "<== TFluka::ProcessRun(" << nevent << ") called."
 303          << endl;
 304   // Write fluka specific scoring output
 305   newplo();
 306
 307   return kTRUE;
 308 }
 309
 310 //_____________________________________________________________________________
 311 // methods for building/management of geometry
 312
 313 // functions from GCONS
 314 //____________________________________________________________________________
 315 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
 316                     Float_t &dens, Float_t &radl, Float_t &absl,
 317                     Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
 318 //
 319    TGeoMaterial *mat;
 320    TIter next (gGeoManager->GetListOfMaterials());
 321    while ((mat = (TGeoMaterial*)next())) {
 322      if (mat->GetUniqueID() == (UInt_t)imat) break;
 323    }
 324    if (!mat) {
 325       Error("Gfmate", "no material with index %i found", imat);
 326       return;
 327    }
 328    sprintf(name, "%s", mat->GetName());
 329    a = mat->GetA();
 330    z = mat->GetZ();
 331    dens = mat->GetDensity();
 332    radl = mat->GetRadLen();
 333    absl = mat->GetIntLen();
 334 }
 335
 336 //______________________________________________________________________________
 337 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
 338                     Double_t &dens, Double_t &radl, Double_t &absl,
 339                     Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
 340 //
 341    TGeoMaterial *mat;
 342    TIter next (gGeoManager->GetListOfMaterials());
 343    while ((mat = (TGeoMaterial*)next())) {
 344      if (mat->GetUniqueID() == (UInt_t)imat) break;
 345    }
 346    if (!mat) {
 347       Error("Gfmate", "no material with index %i found", imat);
 348       return;
 349    }
 350    sprintf(name, "%s", mat->GetName());
 351    a = mat->GetA();
 352    z = mat->GetZ();
 353    dens = mat->GetDensity();
 354    radl = mat->GetRadLen();
 355    absl = mat->GetIntLen();
 356 }
 357
 358 // detector composition
 359 //______________________________________________________________________________
 360 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 361                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 362                       Float_t* buf, Int_t nwbuf) {
 363 //
 364    Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);
 365    Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
 366    delete [] dbuf;
 367 }
 368
 369 //______________________________________________________________________________
 370 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 371                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 372                       Double_t* /*buf*/, Int_t /*nwbuf*/) {
 373 //
 374 // Define a material
 375   TGeoMaterial *mat;
 376   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 377   if ((z-Int_t(z)) > 1E-3) {
 378      mat = fGeom->GetMakeWrongMaterial(z);
 379      if (mat) {
 380         mat->SetRadLen(radl,absl);
 381         mat->SetUniqueID(kmat);
 382         return;
 383      }
 384   }
 385   gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
 386 }
 387
 388 //______________________________________________________________________________
 389 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
 390                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
 391 //
 392 // Define a material mixture
 393 //
 394   Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));
 395   Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));
 396   Double_t* dwmat = fGeom->CreateDoubleArray(wmat, TMath::Abs(nlmat));
 397
 398   Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
 399   for (Int_t i=0; i<nlmat; i++) {
 400     a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
 401   }
 402
 403   delete [] da;
 404   delete [] dz;
 405   delete [] dwmat;
 406 }
 407
 408 //______________________________________________________________________________
 409 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
 410                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
 411 //
 412   // Defines mixture OR COMPOUND IMAT as composed by
 413   // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
 414   //
 415   // If NLMAT > 0 then wmat contains the proportion by
 416   // weights of each basic material in the mixture.
 417   //
 418   // If nlmat < 0 then WMAT contains the number of atoms
 419   // of a given kind into the molecule of the COMPOUND
 420   // In this case, WMAT in output is changed to relative
 421   // weigths.
 422   //
 423   Int_t i,j;
 424   if (nlmat < 0) {
 425      nlmat = - nlmat;
 426      Double_t amol = 0;
 427      for (i=0;i<nlmat;i++) {
 428         amol += a[i]*wmat[i];
 429      }
 430      for (i=0;i<nlmat;i++) {
 431         wmat[i] *= a[i]/amol;
 432      }
 433   }
 434   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 435   // Check if we have elements with fractional Z
 436   TGeoMaterial *mat = 0;
 437   TGeoMixture *mix = 0;
 438   Bool_t mixnew = kFALSE;
 439   for (i=0; i<nlmat; i++) {
 440      if (z[i]-Int_t(z[i]) < 1E-3) continue;
 441      // We have found an element with fractional Z -> loop mixtures to look for it
 442      for (j=0; j<kmat; j++) {
 443         mat = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(j);
 444         if (!mat) break;
 445         if (!mat->IsMixture()) continue;
 446         mix = (TGeoMixture*)mat;
 447         if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
 448 //        printf(" FOUND component %i as mixture %s\n", i, mat->GetName());
 449         mixnew = kTRUE;
 450         break;
 451      }
 452      if (!mixnew) Warning("Mixture","%s : cannot find component %i with fractional Z=%f\n", name, i, z[i]);
 453      break;
 454   }
 455   if (mixnew) {
 456      Int_t nlmatnew = nlmat+mix->GetNelements()-1;
 457      Double_t *anew = new Double_t[nlmatnew];
 458      Double_t *znew = new Double_t[nlmatnew];
 459      Double_t *wmatnew = new Double_t[nlmatnew];
 460      Int_t ind=0;
 461      for (j=0; j<nlmat; j++) {
 462         if (j==i) continue;
 463         anew[ind] = a[j];
 464         znew[ind] = z[j];
 465         wmatnew[ind] = wmat[j];
 466         ind++;
 467      }
 468      for (j=0; j<mix->GetNelements(); j++) {
 469         anew[ind] = mix->GetAmixt()[j];
 470         znew[ind] = mix->GetZmixt()[j];
 471         wmatnew[ind] = wmat[i]*mix->GetWmixt()[j];
 472         ind++;
 473      }
 474      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
 475      delete [] anew;
 476      delete [] znew;
 477      delete [] wmatnew;
 478      return;
 479   }
 480   // Now we need to compact identical elements within the mixture
 481   // First check if this happens
 482   mixnew = kFALSE;
 483   for (i=0; i<nlmat-1; i++) {
 484      for (j=i+1; j<nlmat; j++) {
 485         if (z[i] == z[j]) {
 486            mixnew = kTRUE;
 487            break;
 488         }
 489      }
 490      if (mixnew) break;
 491   }
 492   if (mixnew) {
 493      Int_t nlmatnew = 0;
 494      Double_t *anew = new Double_t[nlmat];
 495      Double_t *znew = new Double_t[nlmat];
 496      memset(znew, 0, nlmat*sizeof(Double_t));
 497      Double_t *wmatnew = new Double_t[nlmat];
 498      Bool_t skipi;
 499      for (i=0; i<nlmat; i++) {
 500         skipi = kFALSE;
 501         for (j=0; j<nlmatnew; j++) {
 502            if (z[i] == z[j]) {
 503               wmatnew[j] += wmat[i];
 504               skipi = kTRUE;
 505               break;
 506            }
 507         }
 508         if (skipi) continue;
 509         anew[nlmatnew] = a[i];
 510         znew[nlmatnew] = z[i];
 511         wmatnew[nlmatnew] = wmat[i];
 512         nlmatnew++;
 513      }
 514      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
 515      delete [] anew;
 516      delete [] znew;
 517      delete [] wmatnew;
 518      return;
 519    }
 520    gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
 521 }
 522
 523 //______________________________________________________________________________
 524 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 525                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 526                     Double_t stemax, Double_t deemax, Double_t epsil,
 527                     Double_t stmin, Float_t* ubuf, Int_t nbuf) {
 528   // Define a medium
 529   //
 530   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
 531   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 532              epsil, stmin, ubuf, nbuf);
 533 }
 534
 535 //______________________________________________________________________________
 536 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 537                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 538                     Double_t stemax, Double_t deemax, Double_t epsil,
 539                     Double_t stmin, Double_t* ubuf, Int_t nbuf) {
 540   // Define a medium
 541   //
 542   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
 543   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 544              epsil, stmin, ubuf, nbuf);
 545 }
 546
 547 //______________________________________________________________________________
 548 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
 549                     Double_t thetaY, Double_t phiY, Double_t thetaZ,
 550                     Double_t phiZ) {
 551 //
 552   krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
 553   fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
 554 }
 555
 556 //______________________________________________________________________________
 557 void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
 558 //
 559 //
 560 // Check if material is used
 561    if (fVerbosityLevel >=3)
 562        printf("Gstpar called with %6d %5s %12.4e %6d\n", itmed, param, parval, fGeom->GetFlukaMaterial(itmed));
 563    Int_t* reglist;
 564    Int_t nreg;
 565    reglist = fGeom->GetMaterialList(fGeom->GetFlukaMaterial(itmed), nreg);
 566    if (nreg == 0) {
 567        return;
 568    }
 569
 570 //
 571    Bool_t process = kFALSE;
 572    if (strncmp(param, "DCAY",  4) == 0 ||
 573        strncmp(param, "PAIR",  4) == 0 ||
 574        strncmp(param, "COMP",  4) == 0 ||
 575        strncmp(param, "PHOT",  4) == 0 ||
 576        strncmp(param, "PFIS",  4) == 0 ||
 577        strncmp(param, "DRAY",  4) == 0 ||
 578        strncmp(param, "ANNI",  4) == 0 ||
 579        strncmp(param, "BREM",  4) == 0 ||
 580        strncmp(param, "MUNU",  4) == 0 ||
 581        strncmp(param, "CKOV",  4) == 0 ||
 582        strncmp(param, "HADR",  4) == 0 ||
 583        strncmp(param, "LOSS",  4) == 0 ||
 584        strncmp(param, "MULS",  4) == 0 ||
 585        strncmp(param, "RAYL",  4) == 0)
 586    {
 587        process = kTRUE;
 588    }
 589    if (process) {
 590        SetProcess(param, Int_t (parval), fGeom->GetFlukaMaterial(itmed));
 591    } else {
 592        SetCut(param, parval, fGeom->GetFlukaMaterial(itmed));
 593    }
 594 }
 595
 596 // functions from GGEOM
 597 //_____________________________________________________________________________
 598 void TFluka::Gsatt(const char *name, const char *att, Int_t val)
 599 {
 600   // Set visualisation attributes for one volume
 601   char vname[5];
 602   fGeom->Vname(name,vname);
 603   char vatt[5];
 604   fGeom->Vname(att,vatt);
 605   gGeoManager->SetVolumeAttribute(vname, vatt, val);
 606 }
 607
 608 //______________________________________________________________________________
 609 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 610                      Float_t *upar, Int_t np)  {
 611 //
 612     return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
 613 }
 614
 615 //______________________________________________________________________________
 616 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 617                      Double_t *upar, Int_t np)  {
 618 //
 619     return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
 620 }
 621
 622 //______________________________________________________________________________
 623 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
 624                    Int_t iaxis) {
 625 //
 626     fMCGeo->Gsdvn(name, mother, ndiv, iaxis);
 627 }
 628
 629 //______________________________________________________________________________
 630 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
 631                     Int_t iaxis, Double_t c0i, Int_t numed) {
 632 //
 633     fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
 634 }
 635
 636 //______________________________________________________________________________
 637 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
 638                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
 639 //
 640     fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
 641 }
 642
 643 //______________________________________________________________________________
 644 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
 645                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
 646 //
 647     fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
 648 }
 649
 650 //______________________________________________________________________________
 651 void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
 652 //
 653 // Nothing to do with TGeo
 654 }
 655
 656 //______________________________________________________________________________
 657 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
 658                    Double_t x, Double_t y, Double_t z, Int_t irot,
 659                    const char *konly) {
 660 //
 661   fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly);
 662 }
 663
 664 //______________________________________________________________________________
 665 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 666                     Double_t x, Double_t y, Double_t z, Int_t irot,
 667                     const char *konly, Float_t *upar, Int_t np)  {
 668   //
 669   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 670 }
 671
 672 //______________________________________________________________________________
 673 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 674                     Double_t x, Double_t y, Double_t z, Int_t irot,
 675                     const char *konly, Double_t *upar, Int_t np)  {
 676   //
 677   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 678 }
 679
 680 //______________________________________________________________________________
 681 void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
 682 //
 683 // Nothing to do with TGeo
 684 }
 685
 686 //______________________________________________________________________________
 687 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 688                          Float_t* absco, Float_t* effic, Float_t* rindex) {
 689 //
 690 // Set Cerenkov properties for medium itmed
 691 //
 692 // npckov: number of sampling points
 693 // ppckov: energy values
 694 // absco:  absorption length
 695 // effic:  quantum efficiency
 696 // rindex: refraction index
 697 //
 698 //
 699 //
 700 //  Create object holding Cerenkov properties
 701 //
 702     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
 703 //
 704 //  Pass object to medium
 705     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
 706     medium->SetCerenkovProperties(cerenkovProperties);
 707 }
 708
 709 //______________________________________________________________________________
 710 void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
 711                          Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
 712 //
 713 // Not implemented with TGeo - what G4 did ? Any FLUKA card generated?
 714    Warning("SetCerenkov", "Not implemented with TGeo");
 715 }
 716
 717 // Euclid
 718 //______________________________________________________________________________
 719 void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/,
 720                           Int_t /*number*/, Int_t /*nlevel*/) {
 721 //
 722 // Not with TGeo
 723    Warning("WriteEuclid", "Not implemented with TGeo");
 724 }
 725
 726
 727
 728 //_____________________________________________________________________________
 729 // methods needed by the stepping
 730 //____________________________________________________________________________
 731
 732 Int_t TFluka::GetMedium() const {
 733 //
 734 //  Get the medium number for the current fluka region
 735 //
 736     return fGeom->GetMedium(); // this I need to check due to remapping !!!
 737 }
 738
 739
 740
 741 //____________________________________________________________________________
 742 // particle table usage
 743 // ID <--> PDG transformations
 744 //_____________________________________________________________________________
 745 Int_t TFluka::IdFromPDG(Int_t pdg) const
 746 {
 747     //
 748     // Return Fluka code from PDG and pseudo ENDF code
 749
 750     // Catch the feedback photons
 751     if (pdg == 50000051) return (-1);
 752     // MCIHAD() goes from pdg to fluka internal.
 753     Int_t intfluka = mcihad(pdg);
 754     // KPTOIP array goes from internal to official
 755     return GetFlukaKPTOIP(intfluka);
 756 }
 757
 758 //______________________________________________________________________________
 759 Int_t TFluka::PDGFromId(Int_t id) const
 760 {
 761   //
 762   // Return PDG code and pseudo ENDF code from Fluka code
 763   //                      Alpha     He3       Triton    Deuteron  gen. ion  opt. photon
 764     Int_t idSpecial[6] = {10020040, 10020030, 10010030, 10010020, 10000000, 50000050};
 765   // IPTOKP array goes from official to internal
 766
 767     if (id == -1) {
 768 // Cerenkov photon
 769         if (fVerbosityLevel >= 3)
 770             printf("\n PDGFromId: Cerenkov Photon \n");
 771         return  50000050;
 772     }
 773 // Error id
 774     if (id == 0 || id < -6 || id > 250) {
 775         if (fVerbosityLevel >= 3)
 776             printf("PDGFromId: Error id = 0\n");
 777         return -1;
 778     }
 779 // Good id
 780     if (id > 0) {
 781         Int_t intfluka = GetFlukaIPTOKP(id);
 782         if (intfluka == 0) {
 783             if (fVerbosityLevel >= 3)
 784                 printf("PDGFromId: Error intfluka = 0: %d\n", id);
 785             return -1;
 786         } else if (intfluka < 0) {
 787             if (fVerbosityLevel >= 3)
 788                 printf("PDGFromId: Error intfluka < 0: %d\n", id);
 789             return -1;
 790         }
 791         if (fVerbosityLevel >= 3)
 792             printf("mpdgha called with %d %d \n", id, intfluka);
 793         // MPDGHA() goes from fluka internal to pdg.
 794         return mpdgha(intfluka);
 795     } else {
 796         // ions and optical photons
 797         return idSpecial[id + 6];
 798     }
 799 }
 800
 801 void TFluka::StopTrack()
 802 {
 803     // Set stopping conditions
 804     // Works for photons and charged particles
 805     fStopped = kTRUE;
 806 }
 807
 808 //_____________________________________________________________________________
 809 // methods for physics management
 810 //____________________________________________________________________________
 811 //
 812 // set methods
 813 //
 814
 815 void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
 816 {
 817 //  Set process user flag for material imat
 818 //
 819 //
 820 //  Update if already in the list
 821 //
 822     TIter next(fUserConfig);
 823     TFlukaConfigOption* proc;
 824     while((proc = (TFlukaConfigOption*)next()))
 825     {
 826         if (proc->Medium() == imed) {
 827             proc->SetProcess(flagName, flagValue);
 828             return;
 829         }
 830     }
 831     proc = new TFlukaConfigOption(imed);
 832     proc->SetProcess(flagName, flagValue);
 833     fUserConfig->Add(proc);
 834 }
 835
 836 //______________________________________________________________________________
 837 Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
 838 {
 839 //  Set process user flag
 840 //
 841 //
 842     SetProcess(flagName, flagValue, -1);
 843     return kTRUE;
 844 }
 845
 846 //______________________________________________________________________________
 847 void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
 848 {
 849 // Set user cut value for material imed
 850 //
 851     TIter next(fUserConfig);
 852     TFlukaConfigOption* proc;
 853     while((proc = (TFlukaConfigOption*)next()))
 854     {
 855         if (proc->Medium() == imed) {
 856             proc->SetCut(cutName, cutValue);
 857             return;
 858         }
 859     }
 860
 861     proc = new TFlukaConfigOption(imed);
 862     proc->SetCut(cutName, cutValue);
 863     fUserConfig->Add(proc);
 864 }
 865
 866 //______________________________________________________________________________
 867 Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
 868 {
 869 // Set user cut value
 870 //
 871 //
 872     SetCut(cutName, cutValue, -1);
 873     return kTRUE;
 874 }
 875
 876
 877 void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what)
 878 {
 879 //
 880 // Adds a user scoring option to the list
 881 //
 882     TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr,outfile,what);
 883     fUserScore->Add(opt);
 884 }
 885 //______________________________________________________________________________
 886 void TFluka::SetUserScoring(const char* option, Int_t npr, char* outfile, Float_t* what, const char* det1, const char* det2, const char* det3)
 887 {
 888 //
 889 // Adds a user scoring option to the list
 890 //
 891     TFlukaScoringOption* opt = new TFlukaScoringOption(option, "User Scoring", npr, outfile, what, det1, det2, det3);
 892     fUserScore->Add(opt);
 893 }
 894
 895 //______________________________________________________________________________
 896 Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
 897 {
 898   printf("WARNING: Xsec not yet implemented !\n"); return -1.;
 899 }
 900
 901
 902 //______________________________________________________________________________
 903 void TFluka::InitPhysics()
 904 {
 905 //
 906 // Physics initialisation with preparation of FLUKA input cards
 907 //
 908     printf("=>InitPhysics\n");
 909
 910 // Construct file names
 911     FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
 912     TString sFlukaVmcCoreInp = getenv("ALICE_ROOT");
 913     sFlukaVmcCoreInp +="/TFluka/input/";
 914     TString sFlukaVmcTmp = "flukaMat.inp";
 915     TString sFlukaVmcInp = GetInputFileName();
 916     sFlukaVmcCoreInp += GetCoreInputFileName();
 917
 918 // Open files
 919     if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
 920         printf("\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
 921         exit(1);
 922     }
 923     if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
 924         printf("\nCannot open file %s\n",sFlukaVmcTmp.Data());
 925         exit(1);
 926     }
 927     if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
 928         printf("\nCannot open file %s\n",sFlukaVmcInp.Data());
 929         exit(1);
 930     }
 931
 932 // Copy core input file
 933     Char_t sLine[255];
 934     Float_t fEventsPerRun;
 935
 936     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
 937         if (strncmp(sLine,"GEOEND",6) != 0)
 938             fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
 939         else {
 940             fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
 941             goto flukamat;
 942         }
 943     } // end of while until GEOEND card
 944
 945
 946  flukamat:
 947     while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
 948         fprintf(pFlukaVmcInp,"%s\n",sLine);
 949     }
 950
 951     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
 952         if (strncmp(sLine,"START",5) != 0)
 953             fprintf(pFlukaVmcInp,"%s\n",sLine);
 954         else {
 955             sscanf(sLine+10,"%10f",&fEventsPerRun);
 956             goto fin;
 957         }
 958     } //end of while until START card
 959
 960  fin:
 961
 962
 963 // Pass information to configuration objects
 964
 965     Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
 966     TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom);
 967
 968     TIter next(fUserConfig);
 969     TFlukaConfigOption* proc;
 970     while((proc = dynamic_cast<TFlukaConfigOption*> (next()))) proc->WriteFlukaInputCards();
 971 //
 972 // Process Fluka specific scoring options
 973 //
 974     TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom);
 975     Float_t loginp        = 49.0;
 976     Int_t inp             = 0;
 977     Int_t nscore          = fUserScore->GetEntries();
 978
 979     TFlukaScoringOption *mopo = 0x0;
 980     TFlukaScoringOption *mopi = 0x0;
 981
 982     for (Int_t isc = 0; isc < nscore; isc++)
 983     {
 984         mopo = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
 985         char*    fileName = mopo->GetFileName();
 986         Int_t    size     = strlen(fileName);
 987         Float_t  lun      = -1.;
 988 //
 989 // Check if new output file has to be opened
 990         for (Int_t isci = 0; isci < isc; isci++) {
 991             mopi = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
 992             if(strncmp(mopi->GetFileName(), fileName, size)==0) {
 993                 //
 994                 // No, the file already exists
 995                 lun = mopi->GetLun();
 996                 mopo->SetLun(lun);
 997                 break;
 998             }
 999         } // inner loop
1000
1001         if (lun == -1.) {
1002             // Open new output file
1003             inp++;
1004             mopo->SetLun(loginp + inp);
1005             mopo->WriteOpenFlukaFile();
1006         }
1007         mopo->WriteFlukaInputCards();
1008     }
1009
1010 // Add START and STOP card
1011     fprintf(pFlukaVmcInp,"START     %10.1f\n",fEventsPerRun);
1012     fprintf(pFlukaVmcInp,"STOP      \n");
1013
1014
1015 // Close files
1016    fclose(pFlukaVmcCoreInp);
1017    fclose(pFlukaVmcFlukaMat);
1018    fclose(pFlukaVmcInp);
1019
1020
1021 //
1022 // Initialisation needed for Cerenkov photon production and transport
1023     TObjArray *matList = GetFlukaMaterials();
1024     Int_t nmaterial =  matList->GetEntriesFast();
1025     fMaterials = new Int_t[nmaterial+3];
1026
1027     for (Int_t im = 0; im < nmaterial; im++)
1028     {
1029         TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
1030         Int_t idmat = material->GetIndex();
1031         fMaterials[idmat] = im;
1032     }
1033 } // end of InitPhysics
1034
1035
1036 //______________________________________________________________________________
1037 void TFluka::SetMaxStep(Double_t step)
1038 {
1039 // Set the maximum step size
1040     if (step > 1.e4) return;
1041
1042     Int_t mreg, latt;
1043     fGeom->GetCurrentRegion(mreg, latt);
1044     STEPSZ.stepmx[mreg - 1] = step;
1045 }
1046
1047
1048 Double_t TFluka::MaxStep() const
1049 {
1050 // Return the maximum for current medium
1051     Int_t mreg, latt;
1052     fGeom->GetCurrentRegion(mreg, latt);
1053     return (STEPSZ.stepmx[mreg - 1]);
1054 }
1055
1056 //______________________________________________________________________________
1057 void TFluka::SetMaxNStep(Int_t)
1058 {
1059 // SetMaxNStep is dummy procedure in TFluka !
1060   if (fVerbosityLevel >=3)
1061   cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1062 }
1063
1064 //______________________________________________________________________________
1065 void TFluka::SetUserDecay(Int_t)
1066 {
1067 // SetUserDecay is dummy procedure in TFluka !
1068   if (fVerbosityLevel >=3)
1069   cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1070 }
1071
1072 //
1073 // dynamic properties
1074 //
1075 //______________________________________________________________________________
1076 void TFluka::TrackPosition(TLorentzVector& position) const
1077 {
1078 // Return the current position in the master reference frame of the
1079 // track being transported
1080 // TRACKR.atrack = age of the particle
1081 // TRACKR.xtrack = x-position of the last point
1082 // TRACKR.ytrack = y-position of the last point
1083 // TRACKR.ztrack = z-position of the last point
1084   Int_t caller = GetCaller();
1085   if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
1086     position.SetX(GetXsco());
1087     position.SetY(GetYsco());
1088     position.SetZ(GetZsco());
1089     position.SetT(TRACKR.atrack);
1090   }
1091   else if (caller == 4) { // mgdraw
1092     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1093     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1094     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1095     position.SetT(TRACKR.atrack);
1096   }
1097   else if (caller == 5) { // sodraw
1098     position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1099     position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1100     position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1101     position.SetT(0);
1102   }
1103   else
1104     Warning("TrackPosition","position not available");
1105 }
1106
1107 //______________________________________________________________________________
1108 void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1109 {
1110 // Return the current position in the master reference frame of the
1111 // track being transported
1112 // TRACKR.atrack = age of the particle
1113 // TRACKR.xtrack = x-position of the last point
1114 // TRACKR.ytrack = y-position of the last point
1115 // TRACKR.ztrack = z-position of the last point
1116   Int_t caller = GetCaller();
1117   if (caller == 3 || caller == 6 || caller == 11 || caller == 12) { //bxdraw,endraw,usdraw
1118     x = GetXsco();
1119     y = GetYsco();
1120     z = GetZsco();
1121   }
1122   else if (caller == 4 || caller == 5) { // mgdraw, sodraw
1123     x = TRACKR.xtrack[TRACKR.ntrack];
1124     y = TRACKR.ytrack[TRACKR.ntrack];
1125     z = TRACKR.ztrack[TRACKR.ntrack];
1126   }
1127   else
1128     Warning("TrackPosition","position not available");
1129 }
1130
1131 //______________________________________________________________________________
1132 void TFluka::TrackMomentum(TLorentzVector& momentum) const
1133 {
1134 // Return the direction and the momentum (GeV/c) of the track
1135 // currently being transported
1136 // TRACKR.ptrack = momentum of the particle (not always defined, if
1137 //               < 0 must be obtained from etrack)
1138 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1139 // TRACKR.etrack = total energy of the particle
1140 // TRACKR.jtrack = identity number of the particle
1141 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1142   Int_t caller = GetCaller();
1143   if (caller != 2) { // not eedraw
1144     if (TRACKR.ptrack >= 0) {
1145       momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1146       momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1147       momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1148       momentum.SetE(TRACKR.etrack);
1149       return;
1150     }
1151     else {
1152       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1153       momentum.SetPx(p*TRACKR.cxtrck);
1154       momentum.SetPy(p*TRACKR.cytrck);
1155       momentum.SetPz(p*TRACKR.cztrck);
1156       momentum.SetE(TRACKR.etrack);
1157       return;
1158     }
1159   }
1160   else
1161     Warning("TrackMomentum","momentum not available");
1162 }
1163
1164 //______________________________________________________________________________
1165 void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1166 {
1167 // Return the direction and the momentum (GeV/c) of the track
1168 // currently being transported
1169 // TRACKR.ptrack = momentum of the particle (not always defined, if
1170 //               < 0 must be obtained from etrack)
1171 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1172 // TRACKR.etrack = total energy of the particle
1173 // TRACKR.jtrack = identity number of the particle
1174 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1175   Int_t caller = GetCaller();
1176   if (caller != 2) { // not eedraw
1177     if (TRACKR.ptrack >= 0) {
1178       px = TRACKR.ptrack*TRACKR.cxtrck;
1179       py = TRACKR.ptrack*TRACKR.cytrck;
1180       pz = TRACKR.ptrack*TRACKR.cztrck;
1181       e = TRACKR.etrack;
1182       return;
1183     }
1184     else {
1185       Double_t p = sqrt(TRACKR.etrack*TRACKR.etrack - PAPROP.am[TRACKR.jtrack+6]*PAPROP.am[TRACKR.jtrack+6]);
1186       px = p*TRACKR.cxtrck;
1187       py = p*TRACKR.cytrck;
1188       pz = p*TRACKR.cztrck;
1189       e = TRACKR.etrack;
1190       return;
1191     }
1192   }
1193   else
1194     Warning("TrackMomentum","momentum not available");
1195 }
1196
1197 //______________________________________________________________________________
1198 Double_t TFluka::TrackStep() const
1199 {
1200 // Return the length in centimeters of the current step
1201 // TRACKR.ctrack = total curved path
1202   Int_t caller = GetCaller();
1203   if (caller == 11 || caller==12 || caller == 3 || caller == 6) //bxdraw,endraw,usdraw
1204     return 0.0;
1205   else if (caller == 4) //mgdraw
1206     return TRACKR.ctrack;
1207   else
1208     return -1.0;
1209 }
1210
1211 //______________________________________________________________________________
1212 Double_t TFluka::TrackLength() const
1213 {
1214 // TRACKR.cmtrck = cumulative curved path since particle birth
1215   Int_t caller = GetCaller();
1216   if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
1217     return TRACKR.cmtrck;
1218   else
1219     return -1.0;
1220 }
1221
1222 //______________________________________________________________________________
1223 Double_t TFluka::TrackTime() const
1224 {
1225 // Return the current time of flight of the track being transported
1226 // TRACKR.atrack = age of the particle
1227   Int_t caller = GetCaller();
1228   if (caller == 11 || caller==12 || caller == 3 || caller == 4 || caller == 6) //bxdraw,endraw,mgdraw,usdraw
1229     return TRACKR.atrack;
1230   else
1231     return -1;
1232 }
1233
1234 //______________________________________________________________________________
1235 Double_t TFluka::Edep() const
1236 {
1237 // Energy deposition
1238 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1239 // -->local energy deposition (the value and the point are not recorded in TRACKR)
1240 //    but in the variable "rull" of the procedure "endraw.cxx"
1241 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1242 // -->no energy loss along the track
1243 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1244 // -->energy loss distributed along the track
1245 // TRACKR.dtrack = energy deposition of the jth deposition event
1246
1247   // If coming from bxdraw we have 2 steps of 0 length and 0 edep
1248   Int_t caller = GetCaller();
1249   if (caller == 11 || caller==12) return 0.0;
1250   Double_t sum = 0;
1251   for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1252     sum +=TRACKR.dtrack[j];
1253   }
1254   if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1255     return fRull + sum;
1256   else {
1257     return sum;
1258   }
1259 }
1260
1261 //______________________________________________________________________________
1262 Int_t TFluka::TrackPid() const
1263 {
1264 // Return the id of the particle transported
1265 // TRACKR.jtrack = identity number of the particle
1266   Int_t caller = GetCaller();
1267   if (caller != 2) { // not eedraw
1268       return PDGFromId(TRACKR.jtrack);
1269   }
1270   else
1271     return -1000;
1272 }
1273
1274 //______________________________________________________________________________
1275 Double_t TFluka::TrackCharge() const
1276 {
1277 // Return charge of the track currently transported
1278 // PAPROP.ichrge = electric charge of the particle
1279 // TRACKR.jtrack = identity number of the particle
1280   Int_t caller = GetCaller();
1281   if (caller != 2)  // not eedraw
1282     return PAPROP.ichrge[TRACKR.jtrack+6];
1283   else
1284     return -1000.0;
1285 }
1286
1287 //______________________________________________________________________________
1288 Double_t TFluka::TrackMass() const
1289 {
1290 // PAPROP.am = particle mass in GeV
1291 // TRACKR.jtrack = identity number of the particle
1292   Int_t caller = GetCaller();
1293   if (caller != 2)  // not eedraw
1294     return PAPROP.am[TRACKR.jtrack+6];
1295   else
1296     return -1000.0;
1297 }
1298
1299 //______________________________________________________________________________
1300 Double_t TFluka::Etot() const
1301 {
1302 // TRACKR.etrack = total energy of the particle
1303   Int_t caller = GetCaller();
1304   if (caller != 2)  // not eedraw
1305     return TRACKR.etrack;
1306   else
1307     return -1000.0;
1308 }
1309
1310 //
1311 // track status
1312 //
1313 //______________________________________________________________________________
1314 Bool_t   TFluka::IsNewTrack() const
1315 {
1316 // Return true for the first call of Stepping()
1317    return fTrackIsNew;
1318 }
1319
1320 void     TFluka::SetTrackIsNew(Bool_t flag)
1321 {
1322 // Return true for the first call of Stepping()
1323    fTrackIsNew = flag;
1324
1325 }
1326
1327
1328 //______________________________________________________________________________
1329 Bool_t   TFluka::IsTrackInside() const
1330 {
1331 // True if the track is not at the boundary of the current volume
1332 // In Fluka a step is always inside one kind of material
1333 // If the step would go behind the region of one material,
1334 // it will be shortened to reach only the boundary.
1335 // Therefore IsTrackInside() is always true.
1336   Int_t caller = GetCaller();
1337   if (caller == 11 || caller==12)  // bxdraw
1338     return 0;
1339   else
1340     return 1;
1341 }
1342
1343 //______________________________________________________________________________
1344 Bool_t   TFluka::IsTrackEntering() const
1345 {
1346 // True if this is the first step of the track in the current volume
1347
1348   Int_t caller = GetCaller();
1349   if (caller == 11)  // bxdraw entering
1350     return 1;
1351   else return 0;
1352 }
1353
1354 //______________________________________________________________________________
1355 Bool_t   TFluka::IsTrackExiting() const
1356 {
1357 // True if track is exiting volume
1358 //
1359   Int_t caller = GetCaller();
1360   if (caller == 12)  // bxdraw exiting
1361     return 1;
1362   else return 0;
1363 }
1364
1365 //______________________________________________________________________________
1366 Bool_t   TFluka::IsTrackOut() const
1367 {
1368 // True if the track is out of the setup
1369 // means escape
1370 // Icode = 14: escape - call from Kaskad
1371 // Icode = 23: escape - call from Emfsco
1372 // Icode = 32: escape - call from Kasneu
1373 // Icode = 40: escape - call from Kashea
1374 // Icode = 51: escape - call from Kasoph
1375   if (fIcode == 14 ||
1376       fIcode == 23 ||
1377       fIcode == 32 ||
1378       fIcode == 40 ||
1379       fIcode == 51) return 1;
1380   else return 0;
1381 }
1382
1383 //______________________________________________________________________________
1384 Bool_t   TFluka::IsTrackDisappeared() const
1385 {
1386 // means all inelastic interactions and decays
1387 // fIcode from usdraw
1388   if (fIcode == 101 || // inelastic interaction
1389       fIcode == 102 || // particle decay
1390       fIcode == 103 || // delta ray generation by hadron
1391       fIcode == 104 || // direct pair production
1392       fIcode == 105 || // bremsstrahlung (muon)
1393       fIcode == 208 || // bremsstrahlung (electron)
1394       fIcode == 214 || // in-flight annihilation
1395       fIcode == 215 || // annihilation at rest
1396       fIcode == 217 || // pair production
1397       fIcode == 219 || // Compton scattering
1398       fIcode == 221 || // Photoelectric effect
1399       fIcode == 300 || // hadronic interaction
1400       fIcode == 400    // delta-ray
1401       ) return 1;
1402   else return 0;
1403 }
1404
1405 //______________________________________________________________________________
1406 Bool_t   TFluka::IsTrackStop() const
1407 {
1408 // True if the track energy has fallen below the threshold
1409 // means stopped by signal or below energy threshold
1410 // Icode = 12: stopping particle       - call from Kaskad
1411 // Icode = 15: time kill               - call from Kaskad
1412 // Icode = 21: below threshold, iarg=1 - call from Emfsco
1413 // Icode = 22: below threshold, iarg=2 - call from Emfsco
1414 // Icode = 24: time kill               - call from Emfsco
1415 // Icode = 31: below threshold         - call from Kasneu
1416 // Icode = 33: time kill               - call from Kasneu
1417 // Icode = 41: time kill               - call from Kashea
1418 // Icode = 52: time kill               - call from Kasoph
1419   if (fIcode == 12 ||
1420       fIcode == 15 ||
1421       fIcode == 21 ||
1422       fIcode == 22 ||
1423       fIcode == 24 ||
1424       fIcode == 31 ||
1425       fIcode == 33 ||
1426       fIcode == 41 ||
1427       fIcode == 52) return 1;
1428   else return 0;
1429 }
1430
1431 //______________________________________________________________________________
1432 Bool_t   TFluka::IsTrackAlive() const
1433 {
1434 // means not disappeared or not out
1435   if (IsTrackDisappeared() || IsTrackOut() ) return 0;
1436   else return 1;
1437 }
1438
1439 //
1440 // secondaries
1441 //
1442
1443 //______________________________________________________________________________
1444 Int_t TFluka::NSecondaries() const
1445
1446 {
1447 // Number of secondary particles generated in the current step
1448 // FINUC.np = number of secondaries except light and heavy ions
1449 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
1450     Int_t caller = GetCaller();
1451     if (caller == 6)  // valid only after usdraw
1452         return FINUC.np + FHEAVY.npheav;
1453     else if (caller == 50) {
1454         // Cerenkov Photon production
1455         return fNCerenkov;
1456     }
1457     return 0;
1458 } // end of NSecondaries
1459
1460 //______________________________________________________________________________
1461 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
1462                 TLorentzVector& position, TLorentzVector& momentum)
1463 {
1464 // Copy particles from secondary stack to vmc stack
1465 //
1466
1467     Int_t caller = GetCaller();
1468     if (caller == 6) {  // valid only after usdraw
1469         if (FINUC.np > 0) {
1470             // Hadronic interaction
1471             if (isec >= 0 && isec < FINUC.np) {
1472                 particleId = PDGFromId(FINUC.kpart[isec]);
1473                 position.SetX(fXsco);
1474                 position.SetY(fYsco);
1475                 position.SetZ(fZsco);
1476                 position.SetT(TRACKR.atrack);
1477                 momentum.SetPx(FINUC.plr[isec]*FINUC.cxr[isec]);
1478                 momentum.SetPy(FINUC.plr[isec]*FINUC.cyr[isec]);
1479                 momentum.SetPz(FINUC.plr[isec]*FINUC.czr[isec]);
1480                 momentum.SetE(FINUC.tki[isec] + PAPROP.am[FINUC.kpart[isec]+6]);
1481             }
1482             else if (isec >= FINUC.np && isec < FINUC.np + FHEAVY.npheav) {
1483                 Int_t jsec = isec - FINUC.np;
1484                 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
1485                 position.SetX(fXsco);
1486                 position.SetY(fYsco);
1487                 position.SetZ(fZsco);
1488                 position.SetT(TRACKR.atrack);
1489                 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
1490                 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
1491                 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
1492                 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
1493                     momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
1494                 else if (FHEAVY.tkheav[jsec] > 6)
1495                     momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
1496             }
1497             else
1498                 Warning("GetSecondary","isec out of range");
1499         }
1500     } else if (caller == 50) {
1501         Int_t index = OPPHST.lstopp - isec;
1502         position.SetX(OPPHST.xoptph[index]);
1503         position.SetY(OPPHST.yoptph[index]);
1504         position.SetZ(OPPHST.zoptph[index]);
1505         position.SetT(OPPHST.agopph[index]);
1506         Double_t p = OPPHST.poptph[index];
1507
1508         momentum.SetPx(p * OPPHST.txopph[index]);
1509         momentum.SetPy(p * OPPHST.tyopph[index]);
1510         momentum.SetPz(p * OPPHST.tzopph[index]);
1511         momentum.SetE(p);
1512     }
1513     else
1514         Warning("GetSecondary","no secondaries available");
1515
1516 } // end of GetSecondary
1517
1518
1519 //______________________________________________________________________________
1520 TMCProcess TFluka::ProdProcess(Int_t) const
1521
1522 {
1523 // Name of the process that has produced the secondary particles
1524 // in the current step
1525
1526     Int_t mugamma = (TRACKR.jtrack == 7 || TRACKR.jtrack == 10 || TRACKR.jtrack == 11);
1527
1528     if      (fIcode == 102)                  return kPDecay;
1529     else if (fIcode == 104 || fIcode == 217) return kPPair;
1530     else if (fIcode == 219)                  return kPCompton;
1531     else if (fIcode == 221)                  return kPPhotoelectric;
1532     else if (fIcode == 105 || fIcode == 208) return kPBrem;
1533     else if (fIcode == 103 || fIcode == 400) return kPDeltaRay;
1534     else if (fIcode == 210 || fIcode == 212) return kPDeltaRay;
1535     else if (fIcode == 214 || fIcode == 215) return kPAnnihilation;
1536     else if (fIcode == 101)                  return kPHadronic;
1537     else if (fIcode == 101) {
1538         if (!mugamma)                        return kPHadronic;
1539         else if (TRACKR.jtrack == 7)         return kPPhotoFission;
1540         else return kPMuonNuclear;
1541     }
1542     else if (fIcode == 225)                  return kPRayleigh;
1543 // Fluka codes 100, 300 and 400 still to be investigasted
1544     else                                     return kPNoProcess;
1545 }
1546
1547
1548 Int_t TFluka::StepProcesses(TArrayI &proc) const
1549 {
1550   //
1551   // Return processes active in the current step
1552   //
1553     proc.Set(1);
1554     TMCProcess iproc;
1555     switch (fIcode) {
1556     case 15:
1557     case 24:
1558     case 33:
1559     case 41:
1560     case 52:
1561         iproc =  kPTOFlimit;
1562         break;
1563     case 12:
1564     case 14:
1565     case 21:
1566     case 22:
1567     case 23:
1568     case 31:
1569     case 32:
1570     case 40:
1571     case 51:
1572         iproc = kPStop;
1573         break;
1574     case 50:
1575         iproc = kPLightAbsorption;
1576         break;
1577     case 59:
1578         iproc = kPLightRefraction;
1579     case 20:
1580         iproc = kPPhotoelectric;
1581         break;
1582     default:
1583         iproc = ProdProcess(0);
1584     }
1585     proc[0] = iproc;
1586     return 1;
1587 }
1588 //______________________________________________________________________________
1589 Int_t TFluka::VolId2Mate(Int_t id) const
1590 {
1591 //
1592 // Returns the material number for a given volume ID
1593 //
1594    return fMCGeo->VolId2Mate(id);
1595 }
1596
1597 //______________________________________________________________________________
1598 const char* TFluka::VolName(Int_t id) const
1599 {
1600 //
1601 // Returns the volume name for a given volume ID
1602 //
1603    return fMCGeo->VolName(id);
1604 }
1605
1606 //______________________________________________________________________________
1607 Int_t TFluka::VolId(const Text_t* volName) const
1608 {
1609 //
1610 // Converts from volume name to volume ID.
1611 // Time consuming. (Only used during set-up)
1612 // Could be replaced by hash-table
1613 //
1614     char sname[20];
1615     Int_t len;
1616     strncpy(sname, volName, len = strlen(volName));
1617     sname[len] = 0;
1618     while (sname[len - 1] == ' ') sname[--len] = 0;
1619     return fMCGeo->VolId(sname);
1620 }
1621
1622 //______________________________________________________________________________
1623 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
1624 {
1625 //
1626 // Return the logical id and copy number corresponding to the current fluka region
1627 //
1628   if (gGeoManager->IsOutside()) return 0;
1629   TGeoNode *node = gGeoManager->GetCurrentNode();
1630   copyNo = node->GetNumber();
1631   Int_t id = node->GetVolume()->GetNumber();
1632   return id;
1633 }
1634
1635 //______________________________________________________________________________
1636 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
1637 {
1638 //
1639 // Return the logical id and copy number of off'th mother
1640 // corresponding to the current fluka region
1641 //
1642   if (off<0 || off>gGeoManager->GetLevel()) return 0;
1643   if (off==0) return CurrentVolID(copyNo);
1644   TGeoNode *node = gGeoManager->GetMother(off);
1645   if (!node) return 0;
1646   copyNo = node->GetNumber();
1647   return node->GetVolume()->GetNumber();
1648 }
1649
1650 //______________________________________________________________________________
1651 const char* TFluka::CurrentVolName() const
1652 {
1653 //
1654 // Return the current volume name
1655 //
1656   if (gGeoManager->IsOutside()) return 0;
1657   return gGeoManager->GetCurrentVolume()->GetName();
1658 }
1659
1660 //______________________________________________________________________________
1661 const char* TFluka::CurrentVolOffName(Int_t off) const
1662 {
1663 //
1664 // Return the volume name of the off'th mother of the current volume
1665 //
1666   if (off<0 || off>gGeoManager->GetLevel()) return 0;
1667   if (off==0) return CurrentVolName();
1668   TGeoNode *node = gGeoManager->GetMother(off);
1669   if (!node) return 0;
1670   return node->GetVolume()->GetName();
1671 }
1672
1673 //______________________________________________________________________________
1674 Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z,
1675                       Float_t & dens, Float_t & radl, Float_t & absl) const
1676 {
1677 //
1678 //  Return the current medium number and material properties
1679 //
1680   Int_t copy;
1681   Int_t id  =  TFluka::CurrentVolID(copy);
1682   Int_t med =  TFluka::VolId2Mate(id);
1683   TGeoVolume*     vol = gGeoManager->GetCurrentVolume();
1684   TGeoMaterial*   mat = vol->GetMaterial();
1685   a    = mat->GetA();
1686   z    = mat->GetZ();
1687   dens = mat->GetDensity();
1688   radl = mat->GetRadLen();
1689   absl = mat->GetIntLen();
1690
1691   return med;
1692 }
1693
1694 //______________________________________________________________________________
1695 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
1696 {
1697 // Transforms a position from the world reference frame
1698 // to the current volume reference frame.
1699 //
1700 //  Geant3 desription:
1701 //  ==================
1702 //       Computes coordinates XD (in DRS)
1703 //       from known coordinates XM in MRS
1704 //       The local reference system can be initialized by
1705 //         - the tracking routines and GMTOD used in GUSTEP
1706 //         - a call to GMEDIA(XM,NUMED)
1707 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
1708 //             (inverse routine is GDTOM)
1709 //
1710 //        If IFLAG=1  convert coordinates
1711 //           IFLAG=2  convert direction cosinus
1712 //
1713 // ---
1714    Double_t xmL[3], xdL[3];
1715    Int_t i;
1716    for (i=0;i<3;i++) xmL[i]=xm[i];
1717    if (iflag == 1) gGeoManager->MasterToLocal(xmL,xdL);
1718    else            gGeoManager->MasterToLocalVect(xmL,xdL);
1719    for (i=0;i<3;i++) xd[i] = xdL[i];
1720 }
1721
1722 //______________________________________________________________________________
1723 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
1724 {
1725    if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
1726    else            gGeoManager->MasterToLocalVect(xm,xd);
1727 }
1728
1729 //______________________________________________________________________________
1730 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
1731 {
1732 // Transforms a position from the current volume reference frame
1733 // to the world reference frame.
1734 //
1735 //  Geant3 desription:
1736 //  ==================
1737 //  Computes coordinates XM (Master Reference System
1738 //  knowing the coordinates XD (Detector Ref System)
1739 //  The local reference system can be initialized by
1740 //    - the tracking routines and GDTOM used in GUSTEP
1741 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
1742 //        (inverse routine is GMTOD)
1743 //
1744 //   If IFLAG=1  convert coordinates
1745 //      IFLAG=2  convert direction cosinus
1746 //
1747 // ---
1748    Double_t xmL[3], xdL[3];
1749    Int_t i;
1750    for (i=0;i<3;i++) xdL[i] = xd[i];
1751    if (iflag == 1) gGeoManager->LocalToMaster(xdL,xmL);
1752    else            gGeoManager->LocalToMasterVect(xdL,xmL);
1753    for (i=0;i<3;i++) xm[i]=xmL[i];
1754 }
1755
1756 //______________________________________________________________________________
1757 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
1758 {
1759    if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
1760    else            gGeoManager->LocalToMasterVect(xd,xm);
1761 }
1762
1763 //______________________________________________________________________________
1764 TObjArray *TFluka::GetFlukaMaterials()
1765 {
1766    return fGeom->GetMatList();
1767 }
1768
1769 //______________________________________________________________________________
1770 void TFluka::SetMreg(Int_t l)
1771 {
1772 // Set current fluka region
1773    fCurrentFlukaRegion = l;
1774    fGeom->SetMreg(l);
1775 }
1776
1777
1778
1779
1780 TString TFluka::ParticleName(Int_t pdg) const
1781 {
1782     // Return particle name for particle with pdg code pdg.
1783     Int_t ifluka = IdFromPDG(pdg);
1784     return TString((CHPPRP.btype[ifluka+6]), 8);
1785 }
1786
1787
1788 Double_t TFluka::ParticleMass(Int_t pdg) const
1789 {
1790     // Return particle mass for particle with pdg code pdg.
1791     Int_t ifluka = IdFromPDG(pdg);
1792     return (PAPROP.am[ifluka+6]);
1793 }
1794
1795 Double_t TFluka::ParticleCharge(Int_t pdg) const
1796 {
1797     // Return particle charge for particle with pdg code pdg.
1798     Int_t ifluka = IdFromPDG(pdg);
1799     return Double_t(PAPROP.ichrge[ifluka+6]);
1800 }
1801
1802 Double_t TFluka::ParticleLifeTime(Int_t pdg) const
1803 {
1804     // Return particle lifetime for particle with pdg code pdg.
1805     Int_t ifluka = IdFromPDG(pdg);
1806     return (PAPROP.thalf[ifluka+6]);
1807 }
1808
1809 void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
1810 {
1811     // Retrieve particle properties for particle with pdg code pdg.
1812
1813     strcpy(name, ParticleName(pdg).Data());
1814     type   = ParticleMCType(pdg);
1815     mass   = ParticleMass(pdg);
1816     charge = ParticleCharge(pdg);
1817     tlife  = ParticleLifeTime(pdg);
1818 }
1819
1820
1821
1822 #define pushcerenkovphoton pushcerenkovphoton_
1823 #define usersteppingckv    usersteppingckv_
1824
1825
1826 extern "C" {
1827     void pushcerenkovphoton(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
1828                             Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
1829                             Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
1830     {
1831         //
1832         // Pushes one cerenkov photon to the stack
1833         //
1834
1835         TFluka* fluka =  (TFluka*) gMC;
1836         TVirtualMCStack* cppstack = fluka->GetStack();
1837         Int_t parent =  TRACKR.ispusr[mkbmx2-1];
1838         cppstack->PushTrack(0, parent, 50000050,
1839                             px, py, pz, e,
1840                             vx, vy, vz, tof,
1841                             polx, poly, polz,
1842                             kPCerenkov, ntr, wgt, 0);
1843     }
1844
1845     void usersteppingckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
1846     {
1847         //
1848         // Calls stepping in order to signal cerenkov production
1849         //
1850         TFluka *fluka = (TFluka*)gMC;
1851         fluka->SetMreg(mreg);
1852         fluka->SetXsco(x);
1853         fluka->SetYsco(y);
1854         fluka->SetZsco(z);
1855         fluka->SetNCerenkov(nphot);
1856         fluka->SetCaller(50);
1857         printf("userstepping ckv: %10d %10d %13.3f %13.3f %13.2f\n", nphot, mreg, x, y, z);
1858         (TVirtualMCApplication::Instance())->Stepping();
1859     }
1860 }
1861