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 #include <TList.h>
  34
  35 #include "TFluka.h"
  36 #include "TFlukaCodes.h"
  37 #include "TCallf77.h"      //For the fortran calls
  38 #include "Fdblprc.h"       //(DBLPRC) fluka common
  39 #include "Fsourcm.h"       //(SOURCM) fluka common
  40 #include "Fgenstk.h"       //(GENSTK)  fluka common
  41 #include "Fiounit.h"       //(IOUNIT) fluka common
  42 #include "Fpaprop.h"       //(PAPROP) fluka common
  43 #include "Fpart.h"         //(PART)   fluka common
  44 #include "Ftrackr.h"       //(TRACKR) fluka common
  45 #include "Fpaprop.h"       //(PAPROP) fluka common
  46 #include "Ffheavy.h"       //(FHEAVY) fluka common
  47 #include "Fopphst.h"       //(OPPHST) fluka common
  48 #include "Fflkstk.h"       //(FLKSTK) fluka common
  49 #include "Fstepsz.h"       //(STEPSZ) fluka common
  50 #include "Fopphst.h"       //(OPPHST) fluka common
  51 #include "Fltclcm.h"       //(LTCLCM) fluka common
  52 #include "Falldlt.h"       //(ALLDLT) fluka common
  53
  54 #include "TVirtualMC.h"
  55 #include "TMCProcess.h"
  56 #include "TGeoManager.h"
  57 #include "TGeoMaterial.h"
  58 #include "TGeoMedium.h"
  59 #include "TFlukaMCGeometry.h"
  60 #include "TGeoMCGeometry.h"
  61 #include "TFlukaCerenkov.h"
  62 #include "TFlukaConfigOption.h"
  63 #include "TFlukaScoringOption.h"
  64 #include "TLorentzVector.h"
  65 #include "TArrayI.h"
  66 #include "TArrayD.h"
  67 #include "TDatabasePDG.h"
  68 #include "TStopwatch.h"
  69
  70
  71 // Fluka methods that may be needed.
  72 #ifndef WIN32
  73 # define flukam  flukam_
  74 # define fluka_openinp fluka_openinp_
  75 # define fluka_openout fluka_openout_
  76 # define fluka_closeinp fluka_closeinp_
  77 # define mcihad mcihad_
  78 # define mpdgha mpdgha_
  79 # define newplo newplo_
  80 # define genout genout_
  81 # define flkend flkend_
  82 #else
  83 # define flukam  FLUKAM
  84 # define fluka_openinp FLUKA_OPENINP
  85 # define fluka_openout FLUKA_OPENOUT
  86 # define fluka_closeinp FLUKA_CLOSEINP
  87 # define mcihad MCIHAD
  88 # define mpdgha MPDGHA
  89 # define newplo NEWPLO
  90 # define genout GENOUT
  91 # define flkend FLKEND
  92 #endif
  93
  94 extern "C"
  95 {
  96   //
  97   // Prototypes for FLUKA functions
  98   //
  99   void type_of_call flukam(const int&);
 100   void type_of_call newplo();
 101   void type_of_call genout();
 102   void type_of_call flkend();
 103   void type_of_call fluka_openinp(const int&, DEFCHARA);
 104   void type_of_call fluka_openout(const int&, DEFCHARA);
 105   void type_of_call fluka_closeinp(const int&);
 106   int  type_of_call mcihad(const int&);
 107   int  type_of_call mpdgha(const int&);
 108 }
 109
 110 //
 111 // Class implementation for ROOT
 112 //
 113 ClassImp(TFluka)
 114
 115 //
 116 //----------------------------------------------------------------------------
 117 // TFluka constructors and destructors.
 118 //______________________________________________________________________________
 119 TFluka::TFluka()
 120   :TVirtualMC(),
 121    fVerbosityLevel(0),
 122    fNEvent(0),
 123    fInputFileName(""),
 124    fCoreInputFileName(""),
 125    fCaller(kNoCaller),
 126    fIcode(kNoProcess),
 127    fNewReg(-1),
 128    fRull(0),
 129    fXsco(0),
 130    fYsco(0),
 131    fZsco(0),
 132    fPItime(0),
 133    fPIlength(0),
 134    fNPI(0),
 135    fTrackIsEntering(kFALSE),
 136    fTrackIsExiting(kFALSE),
 137    fTrackIsNew(kFALSE),
 138    fFieldFlag(kTRUE),
 139    fDummyBoundary(kFALSE),
 140    fStopped(kFALSE),
 141    fStopEvent(kFALSE),
 142    fStopRun(kFALSE),
 143    fPrimaryElectronIndex(-1),
 144    fLowEnergyNeutronTransport(kFALSE),
 145    fMaterials(0),
 146    fNVolumes(0),
 147    fCurrentFlukaRegion(-1),
 148    fNCerenkov(0),
 149    fGeom(0),
 150    fMCGeo(0),
 151    fUserConfig(0),
 152    fUserScore(0)
 153 {
 154   //
 155   // Default constructor
 156   //
 157     for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
 158 }
 159
 160 //______________________________________________________________________________
 161 TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
 162   :TVirtualMC("TFluka",title, isRootGeometrySupported),
 163    fVerbosityLevel(verbosity),
 164    fNEvent(0),
 165    fInputFileName(""),
 166    fCoreInputFileName(""),
 167    fCaller(kNoCaller),
 168    fIcode(kNoProcess),
 169    fNewReg(-1),
 170    fRull(0),
 171    fXsco(0),
 172    fYsco(0),
 173    fZsco(0),
 174    fPItime(0),
 175    fPIlength(0),
 176    fNPI(0),
 177    fTrackIsEntering(kFALSE),
 178    fTrackIsExiting(kFALSE),
 179    fTrackIsNew(kFALSE),
 180    fFieldFlag(kTRUE),
 181    fDummyBoundary(kFALSE),
 182    fStopped(kFALSE),
 183    fStopEvent(kFALSE),
 184    fStopRun(kFALSE),
 185    fPrimaryElectronIndex(-1),
 186    fLowEnergyNeutronTransport(kFALSE),
 187    fMaterials(0),
 188    fNVolumes(0),
 189    fCurrentFlukaRegion(-1),
 190    fNCerenkov(0),
 191    fGeom(0),
 192    fMCGeo(0),
 193    fUserConfig(new TObjArray(100)),
 194    fUserScore(new TObjArray(100))
 195 {
 196   // create geometry interface
 197     for (Int_t i = 0; i < 4; i++) fPint[i] = 0.;
 198
 199    if (fVerbosityLevel >=3)
 200        cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
 201    SetCoreInputFileName();
 202    SetInputFileName();
 203    fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kFALSE);
 204    fGeom  = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
 205    if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
 206    PrintHeader();
 207 }
 208
 209 //______________________________________________________________________________
 210 TFluka::~TFluka()
 211 {
 212     // Destructor
 213     if (fVerbosityLevel >=3)
 214         cout << "<== TFluka::~TFluka() destructor called." << endl;
 215     if (fMaterials) delete [] fMaterials;
 216
 217 //    delete fGeom;
 218 //    delete fMCGeo;
 219
 220     if (fUserConfig) {
 221         fUserConfig->Delete();
 222         delete fUserConfig;
 223     }
 224
 225     if (fUserScore) {
 226         fUserScore->Delete();
 227         delete fUserScore;
 228     }
 229 }
 230
 231 //
 232 //______________________________________________________________________________
 233 // TFluka control methods
 234 //______________________________________________________________________________
 235 void TFluka::Init() {
 236 //
 237 //  Geometry initialisation
 238 //
 239     if (fVerbosityLevel >=3) cout << "==> TFluka::Init() called." << endl;
 240
 241     if (!gGeoManager) new TGeoManager("geom", "FLUKA geometry");
 242     fApplication->ConstructGeometry();
 243     if (!gGeoManager->IsClosed()) {
 244        TGeoVolume *top = (TGeoVolume*)gGeoManager->GetListOfVolumes()->First();
 245        gGeoManager->SetTopVolume(top);
 246        gGeoManager->CloseGeometry("di");
 247     } else {
 248        TGeoNodeCache *cache = gGeoManager->GetCache();
 249        if (!cache->HasIdArray()) {
 250           Warning("Init", "Node ID tracking must be enabled with TFluka: enabling...\n");
 251           cache->BuildIdArray();
 252        }
 253     }
 254     fNVolumes = fGeom->NofVolumes();
 255     fGeom->CreateFlukaMatFile("flukaMat.inp");
 256     if (fVerbosityLevel >=3) {
 257        printf("== Number of volumes: %i\n ==", fNVolumes);
 258        cout << "\t* InitPhysics() - Prepare input file to be called" << endl;
 259     }
 260
 261     fApplication->InitGeometry();
 262     fApplication->ConstructOpGeometry();
 263     //
 264     // Add ions to PDG Data base
 265     //
 266      AddParticlesToPdgDataBase();
 267     //
 268 }
 269
 270
 271 //______________________________________________________________________________
 272 void TFluka::FinishGeometry() {
 273 //
 274 // Build-up table with region to medium correspondance
 275 //
 276   if (fVerbosityLevel >=3) {
 277     cout << "==> TFluka::FinishGeometry() called." << endl;
 278     printf("----FinishGeometry - applying misalignment if any\n");
 279     cout << "<== TFluka::FinishGeometry() called." << endl;
 280   }
 281   TVirtualMCApplication::Instance()->MisalignGeometry();
 282 }
 283
 284 //______________________________________________________________________________
 285 void TFluka::BuildPhysics() {
 286 //
 287 //  Prepare FLUKA input files and call FLUKA physics initialisation
 288 //
 289
 290     if (fVerbosityLevel >=3)
 291         cout << "==> TFluka::BuildPhysics() called." << endl;
 292
 293
 294     if (fVerbosityLevel >=3) {
 295         TList *medlist = gGeoManager->GetListOfMedia();
 296         TIter next(medlist);
 297         TGeoMedium*   med = 0x0;
 298         TGeoMaterial* mat = 0x0;
 299         Int_t ic = 0;
 300
 301         while((med = (TGeoMedium*)next()))
 302         {
 303             mat = med->GetMaterial();
 304             printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex());
 305             ic++;
 306         }
 307     }
 308
 309
 310     // Prepare input file with the current physics settings
 311
 312     InitPhysics();
 313 //  Open fortran files
 314     const char* fname = fInputFileName;
 315     fluka_openinp(lunin, PASSCHARA(fname));
 316     fluka_openout(11, PASSCHARA("fluka.out"));
 317 //  Read input cards
 318     cout << "==> TFluka::BuildPhysics() Read input cards." << endl;
 319     TStopwatch timer;
 320     timer.Start();
 321     GLOBAL.lfdrtr = true;
 322     flukam(1);
 323     cout << "<== TFluka::BuildPhysics() Read input cards End"
 324          << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
 325 //  Close input file
 326     fluka_closeinp(lunin);
 327 //  Finish geometry
 328     FinishGeometry();
 329 }
 330
 331 //______________________________________________________________________________
 332 void TFluka::ProcessEvent() {
 333 //
 334 // Process one event
 335 //
 336     if (fStopRun) {
 337         Warning("ProcessEvent", "User Run Abortion: No more events handled !\n");
 338         fNEvent += 1;
 339         return;
 340     }
 341
 342     if (fVerbosityLevel >=3)
 343         cout << "==> TFluka::ProcessEvent() called." << endl;
 344     fApplication->GeneratePrimaries();
 345     SOURCM.lsouit = true;
 346     flukam(1);
 347     if (fVerbosityLevel >=3)
 348         cout << "<== TFluka::ProcessEvent() called." << endl;
 349     //
 350     // Increase event number
 351     //
 352     fNEvent += 1;
 353 }
 354
 355 //______________________________________________________________________________
 356 Bool_t TFluka::ProcessRun(Int_t nevent) {
 357 //
 358 // Run steering
 359 //
 360
 361   if (fVerbosityLevel >=3)
 362     cout << "==> TFluka::ProcessRun(" << nevent << ") called."
 363          << endl;
 364
 365   if (fVerbosityLevel >=2) {
 366     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
 367     cout << "\t* Calling flukam again..." << endl;
 368   }
 369
 370   Int_t todo = TMath::Abs(nevent);
 371   for (Int_t ev = 0; ev < todo; ev++) {
 372       TStopwatch timer;
 373       timer.Start();
 374       fApplication->BeginEvent();
 375       ProcessEvent();
 376       fApplication->FinishEvent();
 377       cout << "Event: "<< ev
 378            << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
 379   }
 380
 381   if (fVerbosityLevel >=3)
 382     cout << "<== TFluka::ProcessRun(" << nevent << ") called."
 383          << endl;
 384
 385   // Write fluka specific scoring output
 386   genout();
 387   newplo();
 388   flkend();
 389
 390   return kTRUE;
 391 }
 392
 393 //_____________________________________________________________________________
 394 // methods for building/management of geometry
 395
 396 // functions from GCONS
 397 //____________________________________________________________________________
 398 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
 399                     Float_t &dens, Float_t &radl, Float_t &absl,
 400                     Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
 401 //
 402    TGeoMaterial *mat;
 403    TIter next (gGeoManager->GetListOfMaterials());
 404    while ((mat = (TGeoMaterial*)next())) {
 405      if (mat->GetUniqueID() == (UInt_t)imat) break;
 406    }
 407    if (!mat) {
 408       Error("Gfmate", "no material with index %i found", imat);
 409       return;
 410    }
 411    sprintf(name, "%s", mat->GetName());
 412    a = mat->GetA();
 413    z = mat->GetZ();
 414    dens = mat->GetDensity();
 415    radl = mat->GetRadLen();
 416    absl = mat->GetIntLen();
 417 }
 418
 419 //______________________________________________________________________________
 420 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
 421                     Double_t &dens, Double_t &radl, Double_t &absl,
 422                     Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
 423 //
 424    TGeoMaterial *mat;
 425    TIter next (gGeoManager->GetListOfMaterials());
 426    while ((mat = (TGeoMaterial*)next())) {
 427      if (mat->GetUniqueID() == (UInt_t)imat) break;
 428    }
 429    if (!mat) {
 430       Error("Gfmate", "no material with index %i found", imat);
 431       return;
 432    }
 433    sprintf(name, "%s", mat->GetName());
 434    a = mat->GetA();
 435    z = mat->GetZ();
 436    dens = mat->GetDensity();
 437    radl = mat->GetRadLen();
 438    absl = mat->GetIntLen();
 439 }
 440
 441 // detector composition
 442 //______________________________________________________________________________
 443 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 444                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 445                       Float_t* buf, Int_t nwbuf) {
 446 //
 447    Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);
 448    Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
 449    delete [] dbuf;
 450 }
 451
 452 //______________________________________________________________________________
 453 void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 454                       Double_t z, Double_t dens, Double_t radl, Double_t absl,
 455                       Double_t* /*buf*/, Int_t /*nwbuf*/) {
 456 //
 457 // Define a material
 458   TGeoMaterial *mat;
 459   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 460   if ((z-Int_t(z)) > 1E-3) {
 461      mat = fGeom->GetMakeWrongMaterial(z);
 462      if (mat) {
 463         mat->SetRadLen(radl,absl);
 464         mat->SetUniqueID(kmat);
 465         return;
 466      }
 467   }
 468   gGeoManager->Material(name, a, z, dens, kmat, radl, absl);
 469 }
 470
 471 //______________________________________________________________________________
 472 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
 473                      Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
 474 //
 475 // Define a material mixture
 476 //
 477   Double_t* da = fGeom->CreateDoubleArray(a, TMath::Abs(nlmat));
 478   Double_t* dz = fGeom->CreateDoubleArray(z, TMath::Abs(nlmat));
 479   Double_t* dwmat = fGeom->CreateDoubleArray(wmat, TMath::Abs(nlmat));
 480
 481   Mixture(kmat, name, da, dz, dens, nlmat, dwmat);
 482   for (Int_t i=0; i<nlmat; i++) {
 483     a[i] = da[i]; z[i] = dz[i]; wmat[i] = dwmat[i];
 484   }
 485
 486   delete [] da;
 487   delete [] dz;
 488   delete [] dwmat;
 489 }
 490
 491 //______________________________________________________________________________
 492 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
 493                      Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
 494 //
 495   // Defines mixture OR COMPOUND IMAT as composed by
 496   // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
 497   //
 498   // If NLMAT > 0 then wmat contains the proportion by
 499   // weights of each basic material in the mixture.
 500   //
 501   // If nlmat < 0 then WMAT contains the number of atoms
 502   // of a given kind into the molecule of the COMPOUND
 503   // In this case, WMAT in output is changed to relative
 504   // weigths.
 505   //
 506   Int_t i,j;
 507   if (nlmat < 0) {
 508      nlmat = - nlmat;
 509      Double_t amol = 0;
 510      for (i=0;i<nlmat;i++) {
 511         amol += a[i]*wmat[i];
 512      }
 513      for (i=0;i<nlmat;i++) {
 514         wmat[i] *= a[i]/amol;
 515      }
 516   }
 517   kmat = gGeoManager->GetListOfMaterials()->GetSize();
 518   // Check if we have elements with fractional Z
 519   TGeoMaterial *mat = 0;
 520   TGeoMixture *mix = 0;
 521   Bool_t mixnew = kFALSE;
 522   for (i=0; i<nlmat; i++) {
 523      if (z[i]-Int_t(z[i]) < 1E-3) continue;
 524      // We have found an element with fractional Z -> loop mixtures to look for it
 525      for (j=0; j<kmat; j++) {
 526         mat = (TGeoMaterial*)gGeoManager->GetListOfMaterials()->At(j);
 527         if (!mat) break;
 528         if (!mat->IsMixture()) continue;
 529         mix = (TGeoMixture*)mat;
 530         if (TMath::Abs(z[i]-mix->GetZ()) >1E-3) continue;
 531         mixnew = kTRUE;
 532         break;
 533      }
 534      if (!mixnew) Warning("Mixture","%s : cannot find component %i with fractional Z=%f\n", name, i, z[i]);
 535      break;
 536   }
 537   if (mixnew) {
 538      Int_t nlmatnew = nlmat+mix->GetNelements()-1;
 539      Double_t *anew = new Double_t[nlmatnew];
 540      Double_t *znew = new Double_t[nlmatnew];
 541      Double_t *wmatnew = new Double_t[nlmatnew];
 542      Int_t ind=0;
 543      for (j=0; j<nlmat; j++) {
 544         if (j==i) continue;
 545         anew[ind] = a[j];
 546         znew[ind] = z[j];
 547         wmatnew[ind] = wmat[j];
 548         ind++;
 549      }
 550      for (j=0; j<mix->GetNelements(); j++) {
 551         anew[ind] = mix->GetAmixt()[j];
 552         znew[ind] = mix->GetZmixt()[j];
 553         wmatnew[ind] = wmat[i]*mix->GetWmixt()[j];
 554         ind++;
 555      }
 556      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
 557      delete [] anew;
 558      delete [] znew;
 559      delete [] wmatnew;
 560      return;
 561   }
 562   // Now we need to compact identical elements within the mixture
 563   // First check if this happens
 564   mixnew = kFALSE;
 565   for (i=0; i<nlmat-1; i++) {
 566      for (j=i+1; j<nlmat; j++) {
 567         if (z[i] == z[j]) {
 568            mixnew = kTRUE;
 569            break;
 570         }
 571      }
 572      if (mixnew) break;
 573   }
 574   if (mixnew) {
 575      Int_t nlmatnew = 0;
 576      Double_t *anew = new Double_t[nlmat];
 577      Double_t *znew = new Double_t[nlmat];
 578      memset(znew, 0, nlmat*sizeof(Double_t));
 579      Double_t *wmatnew = new Double_t[nlmat];
 580      Bool_t skipi;
 581      for (i=0; i<nlmat; i++) {
 582         skipi = kFALSE;
 583         for (j=0; j<nlmatnew; j++) {
 584            if (z[i] == z[j]) {
 585               wmatnew[j] += wmat[i];
 586               skipi = kTRUE;
 587               break;
 588            }
 589         }
 590         if (skipi) continue;
 591         anew[nlmatnew] = a[i];
 592         znew[nlmatnew] = z[i];
 593         wmatnew[nlmatnew] = wmat[i];
 594         nlmatnew++;
 595      }
 596      Mixture(kmat, name, anew, znew, dens, nlmatnew, wmatnew);
 597      delete [] anew;
 598      delete [] znew;
 599      delete [] wmatnew;
 600      return;
 601    }
 602   gGeoManager->Mixture(name, a, z, dens, nlmat, wmat, kmat);
 603 }
 604
 605 //______________________________________________________________________________
 606 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 607                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 608                     Double_t stemax, Double_t deemax, Double_t epsil,
 609                     Double_t stmin, Float_t* ubuf, Int_t nbuf) {
 610   // Define a medium
 611   //
 612   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
 613   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 614              epsil, stmin, ubuf, nbuf);
 615 }
 616
 617 //______________________________________________________________________________
 618 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat,
 619                     Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
 620                     Double_t stemax, Double_t deemax, Double_t epsil,
 621                     Double_t stmin, Double_t* ubuf, Int_t nbuf) {
 622   // Define a medium
 623   //
 624   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
 625   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax,
 626              epsil, stmin, ubuf, nbuf);
 627 }
 628
 629 //______________________________________________________________________________
 630 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX,
 631                     Double_t thetaY, Double_t phiY, Double_t thetaZ,
 632                     Double_t phiZ) {
 633 //
 634   krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
 635   fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ);
 636 }
 637
 638 //______________________________________________________________________________
 639 void TFluka::Gstpar(Int_t itmed, const char* param, Double_t parval) {
 640 //
 641 //
 642 //
 643    Bool_t process = kFALSE;
 644    Bool_t modelp  = kFALSE;
 645
 646    if (strncmp(param, "DCAY",  4) == 0 ||
 647        strncmp(param, "PAIR",  4) == 0 ||
 648        strncmp(param, "COMP",  4) == 0 ||
 649        strncmp(param, "PHOT",  4) == 0 ||
 650        strncmp(param, "PFIS",  4) == 0 ||
 651        strncmp(param, "DRAY",  4) == 0 ||
 652        strncmp(param, "ANNI",  4) == 0 ||
 653        strncmp(param, "BREM",  4) == 0 ||
 654        strncmp(param, "MUNU",  4) == 0 ||
 655        strncmp(param, "CKOV",  4) == 0 ||
 656        strncmp(param, "HADR",  4) == 0 ||
 657        strncmp(param, "LOSS",  4) == 0 ||
 658        strncmp(param, "MULS",  4) == 0 ||
 659        strncmp(param, "RAYL",  4) == 0 ||
 660        strncmp(param, "STRA",  4) == 0)
 661    {
 662        process = kTRUE;
 663    }
 664
 665    if (strncmp(param, "PRIMIO_N",  8) == 0 ||
 666        strncmp(param, "PRIMIO_E",  8) == 0)
 667    {
 668        modelp = kTRUE;
 669    }
 670
 671    if (process) {
 672        // Process switch
 673        SetProcess(param, Int_t (parval), itmed);
 674    } else if (modelp) {
 675        // Model parameters
 676        SetModelParameter(param, parval, itmed);
 677    } else {
 678        // Cuts
 679        SetCut(param, parval, itmed);
 680    }
 681
 682
 683 }
 684
 685 // functions from GGEOM
 686 //_____________________________________________________________________________
 687 void TFluka::Gsatt(const char *name, const char *att, Int_t val)
 688 {
 689   // Set visualisation attributes for one volume
 690   char vname[5];
 691   fGeom->Vname(name,vname);
 692   char vatt[5];
 693   fGeom->Vname(att,vatt);
 694   gGeoManager->SetVolumeAttribute(vname, vatt, val);
 695 }
 696
 697 //______________________________________________________________________________
 698 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 699                      Float_t *upar, Int_t np)  {
 700 //
 701     return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
 702 }
 703
 704 //______________________________________________________________________________
 705 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,
 706                      Double_t *upar, Int_t np)  {
 707 //
 708     return fMCGeo->Gsvolu(name, shape, nmed, upar, np);
 709 }
 710
 711 //______________________________________________________________________________
 712 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv,
 713                    Int_t iaxis) {
 714 //
 715     fMCGeo->Gsdvn(name, mother, ndiv, iaxis);
 716 }
 717
 718 //______________________________________________________________________________
 719 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv,
 720                     Int_t iaxis, Double_t c0i, Int_t numed) {
 721 //
 722     fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed);
 723 }
 724
 725 //______________________________________________________________________________
 726 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step,
 727                    Int_t iaxis, Int_t numed, Int_t ndvmx) {
 728 //
 729     fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx);
 730 }
 731
 732 //______________________________________________________________________________
 733 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step,
 734                     Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
 735 //
 736     fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx);
 737 }
 738
 739 //______________________________________________________________________________
 740 void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
 741 //
 742 // Nothing to do with TGeo
 743 }
 744
 745 //______________________________________________________________________________
 746 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,
 747                    Double_t x, Double_t y, Double_t z, Int_t irot,
 748                    const char *konly) {
 749 //
 750   fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly);
 751 }
 752
 753 //______________________________________________________________________________
 754 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 755                     Double_t x, Double_t y, Double_t z, Int_t irot,
 756                     const char *konly, Float_t *upar, Int_t np)  {
 757   //
 758   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 759 }
 760
 761 //______________________________________________________________________________
 762 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,
 763                     Double_t x, Double_t y, Double_t z, Int_t irot,
 764                     const char *konly, Double_t *upar, Int_t np)  {
 765   //
 766   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np);
 767 }
 768
 769 //______________________________________________________________________________
 770 void TFluka::Gsbool(const char* /*onlyVolName*/, const char* /*manyVolName*/) {
 771 //
 772 // Nothing to do with TGeo
 773 }
 774
 775 //______________________________________________________________________
 776 Bool_t TFluka::GetTransformation(const TString &volumePath,TGeoHMatrix &mat)
 777 {
 778     // Returns the Transformation matrix between the volume specified
 779     // by the path volumePath and the Top or mater volume. The format
 780     // of the path volumePath is as follows (assuming ALIC is the Top volume)
 781     // "/ALIC_1/DDIP_1/S05I_2/S05H_1/S05G_3". Here ALIC is the top most
 782     // or master volume which has only 1 instance of. Of all of the daughter
 783     // volumes of ALICE, DDIP volume copy #1 is indicated. Similarly for
 784     // the daughter volume of DDIP is S05I copy #2 and so on.
 785     // Inputs:
 786     //   TString& volumePath  The volume path to the specific volume
 787     //                        for which you want the matrix. Volume name
 788     //                        hierarchy is separated by "/" while the
 789     //                        copy number is appended using a "_".
 790     // Outputs:
 791     //  TGeoHMatrix &mat      A matrix with its values set to those
 792     //                        appropriate to the Local to Master transformation
 793     // Return:
 794     //   A logical value if kFALSE then an error occurred and no change to
 795     //   mat was made.
 796
 797    // We have to preserve the modeler state
 798    return fMCGeo->GetTransformation(volumePath, mat);
 799 }
 800
 801 //______________________________________________________________________
 802 Bool_t TFluka::GetShape(const TString &volumePath,TString &shapeType,
 803                         TArrayD &par)
 804 {
 805     // Returns the shape and its parameters for the volume specified
 806     // by volumeName.
 807     // Inputs:
 808     //   TString& volumeName  The volume name
 809     // Outputs:
 810     //   TString &shapeType   Shape type
 811     //   TArrayD &par         A TArrayD of parameters with all of the
 812     //                        parameters of the specified shape.
 813     // Return:
 814     //   A logical indicating whether there was an error in getting this
 815     //   information
 816    return fMCGeo->GetShape(volumePath, shapeType, par);
 817 }
 818
 819 //______________________________________________________________________
 820 Bool_t TFluka::GetMaterial(const TString &volumeName,
 821                             TString &name,Int_t &imat,
 822                             Double_t &a,Double_t &z,Double_t &dens,
 823                             Double_t &radl,Double_t &inter,TArrayD &par)
 824 {
 825     // Returns the Material and its parameters for the volume specified
 826     // by volumeName.
 827     // Note, Geant3 stores and uses mixtures as an element with an effective
 828     // Z and A. Consequently, if the parameter Z is not integer, then
 829     // this material represents some sort of mixture.
 830     // Inputs:
 831     //   TString& volumeName  The volume name
 832     // Outputs:
 833     //   TSrting   &name       Material name
 834     //   Int_t     &imat       Material index number
 835     //   Double_t  &a          Average Atomic mass of material
 836     //   Double_t  &z          Average Atomic number of material
 837     //   Double_t  &dens       Density of material [g/cm^3]
 838     //   Double_t  &radl       Average radiation length of material [cm]
 839     //   Double_t  &inter      Average interaction length of material [cm]
 840     //   TArrayD   &par        A TArrayD of user defined parameters.
 841     // Return:
 842     //   kTRUE if no errors
 843    return fMCGeo->GetMaterial(volumeName,name,imat,a,z,dens,radl,inter,par);
 844 }
 845
 846 //______________________________________________________________________
 847 Bool_t TFluka::GetMedium(const TString &volumeName,TString &name,
 848                          Int_t &imed,Int_t &nmat,Int_t &isvol,Int_t &ifield,
 849                          Double_t &fieldm,Double_t &tmaxfd,Double_t &stemax,
 850                          Double_t &deemax,Double_t &epsil, Double_t &stmin,
 851                          TArrayD &par)
 852 {
 853     // Returns the Medium and its parameters for the volume specified
 854     // by volumeName.
 855     // Inputs:
 856     //   TString& volumeName  The volume name.
 857     // Outputs:
 858     //   TString  &name       Medium name
 859     //   Int_t    &nmat       Material number defined for this medium
 860     //   Int_t    &imed       The medium index number
 861     //   Int_t    &isvol      volume number defined for this medium
 862     //   Int_t    &iflield    Magnetic field flag
 863     //   Double_t &fieldm     Magnetic field strength
 864     //   Double_t &tmaxfd     Maximum angle of deflection per step
 865     //   Double_t &stemax     Maximum step size
 866     //   Double_t &deemax     Maximum fraction of energy allowed to be lost
 867     //                        to continuous process.
 868     //   Double_t &epsil      Boundary crossing precision
 869     //   Double_t &stmin      Minimum step size allowed
 870     //   TArrayD  &par        A TArrayD of user parameters with all of the
 871     //                        parameters of the specified medium.
 872     // Return:
 873     //   kTRUE if there where no errors
 874    return fMCGeo->GetMedium(volumeName,name,imed,nmat,isvol,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin,par);
 875 }
 876
 877 //______________________________________________________________________________
 878 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 879                          Float_t* absco, Float_t* effic, Float_t* rindex) {
 880 //
 881 // Set Cerenkov properties for medium itmed
 882 //
 883 // npckov: number of sampling points
 884 // ppckov: energy values
 885 // absco:  absorption length
 886 // effic:  quantum efficiency
 887 // rindex: refraction index
 888 //
 889 //
 890 //
 891 //  Create object holding Cerenkov properties
 892 //
 893
 894     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex);
 895 //
 896 //  Pass object to medium
 897     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
 898     medium->SetCerenkovProperties(cerenkovProperties);
 899 }
 900
 901 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 902                          Float_t* absco, Float_t* effic, Float_t* rindex, Float_t* rfl) {
 903 //
 904 // Set Cerenkov properties for medium itmed
 905 //
 906 // npckov: number of sampling points
 907 // ppckov: energy values
 908 // absco:  absorption length
 909 // effic:  quantum efficiency
 910 // rindex: refraction index
 911 // rfl:    reflectivity for boundary to medium itmed
 912 //
 913 //
 914 //  Create object holding Cerenkov properties
 915 //
 916     TFlukaCerenkov* cerenkovProperties = new TFlukaCerenkov(npckov, ppckov, absco, effic, rindex, rfl);
 917 //
 918 //  Pass object to medium
 919     TGeoMedium* medium = gGeoManager->GetMedium(itmed);
 920     medium->SetCerenkovProperties(cerenkovProperties);
 921 }
 922
 923
 924 //______________________________________________________________________________
 925 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t *ppckov,
 926                          Double_t *absco, Double_t *effic, Double_t *rindex) {
 927 //
 928 // Set Cerenkov properties for medium itmed
 929 //
 930 // npckov: number of sampling points
 931 // ppckov: energy values
 932 // absco:  absorption length
 933 // effic:  quantum efficiency
 934 // rindex: refraction index
 935 //
 936
 937 //
 938 //  Double_t version
 939   Float_t* fppckov = CreateFloatArray(ppckov, npckov);
 940   Float_t* fabsco  = CreateFloatArray(absco,  npckov);
 941   Float_t* feffic  = CreateFloatArray(effic,  npckov);
 942   Float_t* frindex = CreateFloatArray(rindex, npckov);
 943
 944   SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex);
 945
 946   delete [] fppckov;
 947   delete [] fabsco;
 948   delete [] feffic;
 949   delete [] frindex;
 950 }
 951
 952 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Double_t* ppckov,
 953                          Double_t* absco, Double_t* effic, Double_t* rindex, Double_t* rfl) {
 954 //
 955 // Set Cerenkov properties for medium itmed
 956 //
 957 // npckov: number of sampling points
 958 // ppckov: energy values
 959 // absco:  absorption length
 960 // effic:  quantum efficiency
 961 // rindex: refraction index
 962 // rfl:    reflectivity for boundary to medium itmed
 963 //
 964
 965 //
 966 // //  Double_t version
 967   Float_t* fppckov = CreateFloatArray(ppckov, npckov);
 968   Float_t* fabsco  = CreateFloatArray(absco,  npckov);
 969   Float_t* feffic  = CreateFloatArray(effic,  npckov);
 970   Float_t* frindex = CreateFloatArray(rindex, npckov);
 971   Float_t* frfl    = CreateFloatArray(rfl,    npckov);
 972
 973   SetCerenkov(itmed, npckov, fppckov, fabsco, feffic, frindex, frfl);
 974
 975   delete [] fppckov;
 976   delete [] fabsco;
 977   delete [] feffic;
 978   delete [] frindex;
 979   delete [] frfl;
 980 }
 981
 982 // Euclid
 983 //______________________________________________________________________________
 984 void TFluka::WriteEuclid(const char* /*fileName*/, const char* /*topVol*/,
 985                           Int_t /*number*/, Int_t /*nlevel*/) {
 986 //
 987 // Not with TGeo
 988    Warning("WriteEuclid", "Not implemented !");
 989 }
 990
 991
 992
 993 //_____________________________________________________________________________
 994 // methods needed by the stepping
 995 //____________________________________________________________________________
 996
 997 Int_t TFluka::GetMedium() const {
 998 //
 999 //  Get the medium number for the current fluka region
1000 //
1001     if (gGeoManager->IsOutside()) {
1002         return (-1);
1003     } else {
1004         return (fGeom->GetMedium()); // this I need to check due to remapping !!!
1005     }
1006 }
1007
1008 //____________________________________________________________________________
1009 Int_t TFluka::GetDummyRegion() const
1010 {
1011 // Returns index of the dummy region.
1012    return fGeom->GetDummyRegion();
1013 }
1014
1015 //____________________________________________________________________________
1016 Int_t TFluka::GetDummyLattice() const
1017 {
1018 // Returns index of the dummy lattice.
1019    return fGeom->GetDummyLattice();
1020 }
1021
1022 //____________________________________________________________________________
1023 // particle table usage
1024 // ID <--> PDG transformations
1025 //_____________________________________________________________________________
1026 Int_t TFluka::IdFromPDG(Int_t pdg) const
1027 {
1028     //
1029     // Return Fluka code from PDG and pseudo ENDF code
1030
1031     // Catch the feedback photons
1032     if (pdg == 50000051) return (kFLUKAoptical);
1033     // MCIHAD() goes from pdg to fluka internal.
1034     Int_t intfluka = mcihad(pdg);
1035     // KPTOIP array goes from internal to official
1036     return GetFlukaKPTOIP(intfluka);
1037 }
1038
1039 //______________________________________________________________________________
1040 Int_t TFluka::PDGFromId(Int_t id) const
1041 {
1042   //
1043   // Return PDG code and pseudo ENDF code from Fluka code
1044   //                      Alpha     He3       Triton    Deuteron  gen. ion  opt. photon
1045     Int_t idSpecial[6] = {GetIonPdg(2,4), GetIonPdg(2, 3), GetIonPdg(1,3), GetIonPdg(1,2), GetIonPdg(0,0), 50000050};
1046   // IPTOKP array goes from official to internal
1047
1048     if (id == kFLUKAoptical) {
1049 // Cerenkov photon
1050 //        if (fVerbosityLevel >= 3)
1051 //            printf("\n PDGFromId: Cerenkov Photon \n");
1052         return  50000050;
1053     }
1054 // Error id
1055     if (id == 0 || id < kFLUKAcodemin || id > kFLUKAcodemax) {
1056         if (fVerbosityLevel >= 3)
1057             printf("PDGFromId: Error id = 0 %5d %5d\n", id, fCaller);
1058         return -1;
1059     }
1060 // Good id
1061     if (id > 0) {
1062         Int_t intfluka = GetFlukaIPTOKP(id);
1063         if (intfluka == 0) {
1064             if (fVerbosityLevel >= 3)
1065                 printf("PDGFromId: Error intfluka = 0: %d\n", id);
1066             return -1;
1067         } else if (intfluka < 0) {
1068             if (fVerbosityLevel >= 3)
1069                 printf("PDGFromId: Error intfluka < 0: %d\n", id);
1070             return -1;
1071         }
1072 //        if (fVerbosityLevel >= 3)
1073 //            printf("mpdgha called with %d %d \n", id, intfluka);
1074         return mpdgha(intfluka);
1075     } else {
1076         // ions and optical photons
1077         return idSpecial[id - kFLUKAcodemin];
1078     }
1079 }
1080
1081 void TFluka::StopTrack()
1082 {
1083     // Set stopping conditions
1084     // Works for photons and charged particles
1085     fStopped = kTRUE;
1086 }
1087
1088 //_____________________________________________________________________________
1089 // methods for physics management
1090 //____________________________________________________________________________
1091 //
1092 // set methods
1093 //
1094
1095 void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
1096 {
1097 //  Set process user flag for material imat
1098 //
1099 //
1100 //  Update if already in the list
1101 //
1102     TIter next(fUserConfig);
1103     TFlukaConfigOption* proc;
1104     while((proc = (TFlukaConfigOption*)next()))
1105     {
1106         if (proc->Medium() == imed) {
1107             proc->SetProcess(flagName, flagValue);
1108             return;
1109         }
1110     }
1111     proc = new TFlukaConfigOption(imed);
1112     proc->SetProcess(flagName, flagValue);
1113     fUserConfig->Add(proc);
1114 }
1115
1116 //______________________________________________________________________________
1117 Bool_t TFluka::SetProcess(const char* flagName, Int_t flagValue)
1118 {
1119 //  Set process user flag
1120 //
1121 //
1122     SetProcess(flagName, flagValue, -1);
1123     return kTRUE;
1124 }
1125
1126 //______________________________________________________________________________
1127 void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
1128 {
1129 // Set user cut value for material imed
1130 //
1131     TIter next(fUserConfig);
1132     TFlukaConfigOption* proc;
1133     while((proc = (TFlukaConfigOption*)next()))
1134     {
1135         if (proc->Medium() == imed) {
1136             proc->SetCut(cutName, cutValue);
1137             return;
1138         }
1139     }
1140
1141     proc = new TFlukaConfigOption(imed);
1142     proc->SetCut(cutName, cutValue);
1143     fUserConfig->Add(proc);
1144 }
1145
1146
1147 //______________________________________________________________________________
1148 void TFluka::SetModelParameter(const char* parName, Double_t parValue, Int_t imed)
1149 {
1150 // Set model parameter for material imed
1151 //
1152     TIter next(fUserConfig);
1153     TFlukaConfigOption* proc;
1154     while((proc = (TFlukaConfigOption*)next()))
1155     {
1156         if (proc->Medium() == imed) {
1157             proc->SetModelParameter(parName, parValue);
1158             return;
1159         }
1160     }
1161
1162     proc = new TFlukaConfigOption(imed);
1163     proc->SetModelParameter(parName, parValue);
1164     fUserConfig->Add(proc);
1165 }
1166
1167 //______________________________________________________________________________
1168 Bool_t TFluka::SetCut(const char* cutName, Double_t cutValue)
1169 {
1170 // Set user cut value
1171 //
1172 //
1173     SetCut(cutName, cutValue, -1);
1174     return kTRUE;
1175 }
1176
1177
1178 void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what)
1179 {
1180 //
1181 // Adds a user scoring option to the list
1182 //
1183     TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr,outfile,what);
1184     fUserScore->Add(opt);
1185 }
1186 //______________________________________________________________________________
1187 void TFluka::SetUserScoring(const char* option, const char* sdum, Int_t npr, char* outfile, Float_t* what,
1188                             const char* det1, const char* det2, const char* det3)
1189 {
1190 //
1191 // Adds a user scoring option to the list
1192 //
1193     TFlukaScoringOption* opt = new TFlukaScoringOption(option, sdum, npr, outfile, what, det1, det2, det3);
1194     fUserScore->Add(opt);
1195 }
1196
1197 //______________________________________________________________________________
1198 Double_t TFluka::Xsec(char*, Double_t, Int_t, Int_t)
1199 {
1200   Warning("Xsec", "Not yet implemented.!\n"); return -1.;
1201 }
1202
1203
1204 //______________________________________________________________________________
1205 void TFluka::InitPhysics()
1206 {
1207 //
1208 // Physics initialisation with preparation of FLUKA input cards
1209 //
1210 // Construct file names
1211     FILE *pFlukaVmcCoreInp, *pFlukaVmcFlukaMat, *pFlukaVmcInp;
1212     TString sFlukaVmcTmp = "flukaMat.inp";
1213     TString sFlukaVmcInp = GetInputFileName();
1214     TString sFlukaVmcCoreInp = GetCoreInputFileName();
1215
1216 // Open files
1217     if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
1218         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
1219         exit(1);
1220     }
1221     if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
1222         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcTmp.Data());
1223         exit(1);
1224     }
1225     if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
1226         Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcInp.Data());
1227         exit(1);
1228     }
1229
1230 // Copy core input file
1231     Char_t sLine[255];
1232     Float_t fEventsPerRun;
1233
1234     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
1235         if (strncmp(sLine,"GEOEND",6) != 0)
1236             fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
1237         else {
1238             fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
1239             goto flukamat;
1240         }
1241     } // end of while until GEOEND card
1242
1243
1244  flukamat:
1245     while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
1246         fprintf(pFlukaVmcInp,"%s\n",sLine);
1247     }
1248
1249     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
1250         if (strncmp(sLine,"START",5) != 0)
1251             fprintf(pFlukaVmcInp,"%s\n",sLine);
1252         else {
1253             sscanf(sLine+10,"%10f",&fEventsPerRun);
1254             goto fin;
1255         }
1256     } //end of while until START card
1257
1258  fin:
1259
1260
1261 // Pass information to configuration objects
1262
1263     Float_t fLastMaterial = fGeom->GetLastMaterialIndex();
1264     TFlukaConfigOption::SetStaticInfo(pFlukaVmcInp, 3, fLastMaterial, fGeom);
1265
1266     TIter next(fUserConfig);
1267     TFlukaConfigOption* proc;
1268     while((proc = dynamic_cast<TFlukaConfigOption*> (next()))) proc->WriteFlukaInputCards();
1269 //
1270 // Process Fluka specific scoring options
1271 //
1272     TFlukaScoringOption::SetStaticInfo(pFlukaVmcInp, fGeom);
1273     Float_t loginp        = -49.0;
1274     Int_t inp             = 0;
1275     Int_t nscore          = fUserScore->GetEntries();
1276
1277     TFlukaScoringOption *mopo = 0;
1278     TFlukaScoringOption *mopi = 0;
1279
1280     for (Int_t isc = 0; isc < nscore; isc++)
1281     {
1282         mopo = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
1283         char*    fileName = mopo->GetFileName();
1284         Int_t    size     = strlen(fileName);
1285         Float_t  lun      = -1.;
1286 //
1287 // Check if new output file has to be opened
1288         for (Int_t isci = 0; isci < isc; isci++) {
1289
1290
1291             mopi = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isci));
1292             if(strncmp(mopi->GetFileName(), fileName, size)==0) {
1293                 //
1294                 // No, the file already exists
1295                 lun = mopi->GetLun();
1296                 mopo->SetLun(lun);
1297                 break;
1298             }
1299         } // inner loop
1300
1301         if (lun == -1.) {
1302             // Open new output file
1303             inp++;
1304             mopo->SetLun(loginp + inp);
1305             mopo->WriteOpenFlukaFile();
1306         }
1307         mopo->WriteFlukaInputCards();
1308     }
1309
1310 // Add RANDOMIZ card
1311     fprintf(pFlukaVmcInp,"RANDOMIZ  %10.1f%10.0f\n", 1., Float_t(gRandom->GetSeed()));
1312 // Add START and STOP card
1313     fprintf(pFlukaVmcInp,"START     %10.1f\n",fEventsPerRun);
1314     fprintf(pFlukaVmcInp,"STOP      \n");
1315
1316
1317 // Close files
1318    fclose(pFlukaVmcCoreInp);
1319    fclose(pFlukaVmcFlukaMat);
1320    fclose(pFlukaVmcInp);
1321
1322
1323 //
1324 // Initialisation needed for Cerenkov photon production and transport
1325     TObjArray *matList = GetFlukaMaterials();
1326     Int_t nmaterial =  matList->GetEntriesFast();
1327     fMaterials = new Int_t[nmaterial+25];
1328
1329     for (Int_t im = 0; im < nmaterial; im++)
1330     {
1331         TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
1332         Int_t idmat = material->GetIndex();
1333         fMaterials[idmat] = im;
1334     }
1335 } // end of InitPhysics
1336
1337
1338 //______________________________________________________________________________
1339 void TFluka::SetMaxStep(Double_t step)
1340 {
1341 // Set the maximum step size
1342 //    if (step > 1.e4) return;
1343
1344 //    Int_t mreg=0, latt=0;
1345 //    fGeom->GetCurrentRegion(mreg, latt);
1346
1347
1348     Int_t mreg = fGeom->GetCurrentRegion();
1349     STEPSZ.stepmx[mreg - 1] = step;
1350 }
1351
1352
1353 Double_t TFluka::MaxStep() const
1354 {
1355 // Return the maximum for current medium
1356     Int_t mreg, latt;
1357     fGeom->GetCurrentRegion(mreg, latt);
1358     return (STEPSZ.stepmx[mreg - 1]);
1359 }
1360
1361 //______________________________________________________________________________
1362 void TFluka::SetMaxNStep(Int_t)
1363 {
1364 // SetMaxNStep is dummy procedure in TFluka !
1365   if (fVerbosityLevel >=3)
1366   cout << "SetMaxNStep is dummy procedure in TFluka !" << endl;
1367 }
1368
1369 //______________________________________________________________________________
1370 void TFluka::SetUserDecay(Int_t)
1371 {
1372 // SetUserDecay is dummy procedure in TFluka !
1373   if (fVerbosityLevel >=3)
1374   cout << "SetUserDecay is dummy procedure in TFluka !" << endl;
1375 }
1376
1377 //
1378 // dynamic properties
1379 //
1380 //______________________________________________________________________________
1381 void TFluka::TrackPosition(TLorentzVector& position) const
1382 {
1383 // Return the current position in the master reference frame of the
1384 // track being transported
1385 // TRACKR.atrack = age of the particle
1386 // TRACKR.xtrack = x-position of the last point
1387 // TRACKR.ytrack = y-position of the last point
1388 // TRACKR.ztrack = z-position of the last point
1389   FlukaCallerCode_t caller = GetCaller();
1390   if (caller == kENDRAW    || caller == kUSDRAW ||
1391       caller == kBXExiting || caller == kBXEntering ||
1392       caller == kUSTCKV) {
1393       position.SetX(GetXsco());
1394       position.SetY(GetYsco());
1395       position.SetZ(GetZsco());
1396       position.SetT(TRACKR.atrack);
1397   }
1398   else if (caller == kMGDRAW) {
1399       Int_t i = -1;
1400       if ((i = fPrimaryElectronIndex) > -1) {
1401           // Primary Electron Ionisation
1402           Double_t x, y, z, t;
1403           GetPrimaryElectronPosition(i, x, y, z, t);
1404           position.SetX(x);
1405           position.SetY(y);
1406           position.SetZ(z);
1407           position.SetT(t);
1408       } else {
1409           position.SetX(TRACKR.xtrack[TRACKR.ntrack]);
1410           position.SetY(TRACKR.ytrack[TRACKR.ntrack]);
1411           position.SetZ(TRACKR.ztrack[TRACKR.ntrack]);
1412           position.SetT(TRACKR.atrack);
1413       }
1414   }
1415   else if (caller == kSODRAW) {
1416       Int_t ist = FLKSTK.npflka;
1417       position.SetX(FLKSTK.xflk[ist]);
1418       position.SetY(FLKSTK.yflk[ist]);
1419       position.SetZ(FLKSTK.zflk[ist]);
1420       position.SetT(FLKSTK.agestk[ist]);
1421   } else if (caller == kMGResumedTrack) {
1422       position.SetX(TRACKR.spausr[0]);
1423       position.SetY(TRACKR.spausr[1]);
1424       position.SetZ(TRACKR.spausr[2]);
1425       position.SetT(TRACKR.spausr[3]);
1426   }
1427   else
1428       Warning("TrackPosition","position not available");
1429 }
1430
1431 //______________________________________________________________________________
1432 void TFluka::TrackPosition(Double_t& x, Double_t& y, Double_t& z) const
1433 {
1434 // Return the current position in the master reference frame of the
1435 // track being transported
1436 // TRACKR.atrack = age of the particle
1437 // TRACKR.xtrack = x-position of the last point
1438 // TRACKR.ytrack = y-position of the last point
1439 // TRACKR.ztrack = z-position of the last point
1440   FlukaCallerCode_t caller = GetCaller();
1441   if (caller == kENDRAW    || caller == kUSDRAW ||
1442       caller == kBXExiting || caller == kBXEntering ||
1443       caller == kUSTCKV) {
1444       x = GetXsco();
1445       y = GetYsco();
1446       z = GetZsco();
1447   }
1448   else if (caller == kMGDRAW) {
1449       Int_t i = -1;
1450       if ((i = fPrimaryElectronIndex) > -1) {
1451           Double_t t;
1452           GetPrimaryElectronPosition(i, x, y, z, t);
1453       } else {
1454           x = TRACKR.xtrack[TRACKR.ntrack];
1455           y = TRACKR.ytrack[TRACKR.ntrack];
1456           z = TRACKR.ztrack[TRACKR.ntrack];
1457       }
1458   }
1459   else if (caller == kSODRAW) {
1460       Int_t ist = FLKSTK.npflka;
1461       x = FLKSTK.xflk[ist];
1462       y = FLKSTK.yflk[ist];
1463       z = FLKSTK.zflk[ist];
1464   }
1465   else if (caller == kMGResumedTrack) {
1466       x = TRACKR.spausr[0];
1467       y = TRACKR.spausr[1];
1468       z = TRACKR.spausr[2];
1469   }
1470   else
1471       Warning("TrackPosition","position not available");
1472 }
1473
1474 //______________________________________________________________________________
1475 void TFluka::TrackMomentum(TLorentzVector& momentum) const
1476 {
1477 // Return the direction and the momentum (GeV/c) of the track
1478 // currently being transported
1479 // TRACKR.ptrack = momentum of the particle (not always defined, if
1480 //               < 0 must be obtained from etrack)
1481 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1482 // TRACKR.etrack = total energy of the particle
1483 // TRACKR.jtrack = identity number of the particle
1484 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1485   FlukaCallerCode_t  caller = GetCaller();
1486   FlukaProcessCode_t icode  = GetIcode();
1487
1488   if (caller  != kEEDRAW         &&
1489       caller  != kMGResumedTrack &&
1490       caller  != kSODRAW         &&
1491       caller  != kUSDRAW         &&
1492       (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
1493       if (TRACKR.ptrack >= 0) {
1494           momentum.SetPx(TRACKR.ptrack*TRACKR.cxtrck);
1495           momentum.SetPy(TRACKR.ptrack*TRACKR.cytrck);
1496           momentum.SetPz(TRACKR.ptrack*TRACKR.cztrck);
1497           momentum.SetE(TRACKR.etrack);
1498           return;
1499       }
1500       else {
1501           Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) * ParticleMassFPC(TRACKR.jtrack));
1502           momentum.SetPx(p*TRACKR.cxtrck);
1503           momentum.SetPy(p*TRACKR.cytrck);
1504           momentum.SetPz(p*TRACKR.cztrck);
1505           momentum.SetE(TRACKR.etrack);
1506           return;
1507       }
1508   } else if  (caller == kMGResumedTrack) {
1509       momentum.SetPx(TRACKR.spausr[4]);
1510       momentum.SetPy(TRACKR.spausr[5]);
1511       momentum.SetPz(TRACKR.spausr[6]);
1512       momentum.SetE (TRACKR.spausr[7]);
1513       return;
1514   } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1515       momentum.SetPx(0.);
1516       momentum.SetPy(0.);
1517       momentum.SetPz(0.);
1518       momentum.SetE(TrackMass());
1519
1520   } else if (caller == kSODRAW) {
1521       Int_t ist  = FLKSTK.npflka;
1522       Double_t p = FLKSTK.pmoflk[ist];
1523       Int_t ifl  = FLKSTK.iloflk[ist];
1524       Double_t m = PAPROP.am[ifl + 6];
1525       Double_t e = TMath::Sqrt(p * p + m * m);
1526       momentum.SetPx(p * FLKSTK.txflk[ist]);
1527       momentum.SetPy(p * FLKSTK.tyflk[ist]);
1528       momentum.SetPz(p * FLKSTK.tzflk[ist]);
1529       momentum.SetE(e);
1530   } else if (caller == kUSDRAW) {
1531       if (icode == kEMFSCObrems  ||
1532           icode == kEMFSCOmoller ||
1533           icode == kEMFSCObhabha ||
1534           icode == kEMFSCOcompton )
1535       {
1536           momentum.SetPx(fPint[0]);
1537           momentum.SetPy(fPint[1]);
1538           momentum.SetPz(fPint[2]);
1539           momentum.SetE(fPint[3]);
1540       } else if (icode == kKASKADdray  ||
1541                  icode == kKASKADbrems ||
1542                  icode == kKASKADpair) {
1543           momentum.SetPx(GENSTK.plr[0] * GENSTK.cxr[0]);
1544           momentum.SetPy(GENSTK.plr[0] * GENSTK.cyr[0]);
1545           momentum.SetPz(GENSTK.plr[0] * GENSTK.czr[0]);
1546           momentum.SetE (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
1547       } else {
1548           Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack
1549                             - ParticleMassFPC(TRACKR.jtrack)
1550                             * ParticleMassFPC(TRACKR.jtrack));
1551           momentum.SetPx(p*TRACKR.cxtrck);
1552           momentum.SetPy(p*TRACKR.cytrck);
1553           momentum.SetPz(p*TRACKR.cztrck);
1554           momentum.SetE(TRACKR.etrack);
1555       }
1556   }
1557   else
1558     Warning("TrackMomentum","momentum not available");
1559 }
1560
1561 //______________________________________________________________________________
1562 void TFluka::TrackMomentum(Double_t& px, Double_t& py, Double_t& pz, Double_t& e) const
1563 {
1564 // Return the direction and the momentum (GeV/c) of the track
1565 // currently being transported
1566 // TRACKR.ptrack = momentum of the particle (not always defined, if
1567 //               < 0 must be obtained from etrack)
1568 // TRACKR.cx,y,ztrck = direction cosines of the current particle
1569 // TRACKR.etrack = total energy of the particle
1570 // TRACKR.jtrack = identity number of the particle
1571 // PAPROP.am[TRACKR.jtrack] = particle mass in gev
1572   FlukaCallerCode_t   caller = GetCaller();
1573   FlukaProcessCode_t  icode  = GetIcode();
1574   if (caller != kEEDRAW         &&
1575       caller != kMGResumedTrack &&
1576       caller != kSODRAW         &&
1577       caller != kUSDRAW         &&
1578       (caller != kENDRAW || (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2))) {
1579     if (TRACKR.ptrack >= 0) {
1580       px = TRACKR.ptrack*TRACKR.cxtrck;
1581       py = TRACKR.ptrack*TRACKR.cytrck;
1582       pz = TRACKR.ptrack*TRACKR.cztrck;
1583       e  = TRACKR.etrack;
1584       return;
1585     }
1586     else {
1587       Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) *  ParticleMassFPC(TRACKR.jtrack));
1588       px = p*TRACKR.cxtrck;
1589       py = p*TRACKR.cytrck;
1590       pz = p*TRACKR.cztrck;
1591       e  = TRACKR.etrack;
1592       return;
1593     }
1594   } else if (caller == kMGResumedTrack) {
1595       px = TRACKR.spausr[4];
1596       py = TRACKR.spausr[5];
1597       pz = TRACKR.spausr[6];
1598       e  = TRACKR.spausr[7];
1599       return;
1600   } else if (caller == kENDRAW && (icode == kEMFSCOstopping1 || icode == kEMFSCOstopping2)) {
1601       px = 0.;
1602       py = 0.;
1603       pz = 0.;
1604       e  = TrackMass();
1605   } else if (caller == kSODRAW) {
1606       Int_t ist  = FLKSTK.npflka;
1607       Double_t p = FLKSTK.pmoflk[ist];
1608       Int_t ifl  = FLKSTK.iloflk[ist];
1609       Double_t m = PAPROP.am[ifl + 6];
1610                e = TMath::Sqrt(p * p + m * m);
1611       px = p * FLKSTK.txflk[ist];
1612       py = p * FLKSTK.tyflk[ist];
1613       pz = p * FLKSTK.tzflk[ist];
1614   } else if (caller == kUSDRAW) {
1615       if (icode == kEMFSCObrems  ||
1616           icode == kEMFSCOmoller ||
1617           icode == kEMFSCObhabha ||
1618           icode == kEMFSCOcompton )
1619       {
1620           px = fPint[0];
1621           py = fPint[1];
1622           pz = fPint[2];
1623           e  = fPint[3];
1624       } else if (icode == kKASKADdray  ||
1625                  icode == kKASKADbrems ||
1626                  icode == kKASKADpair) {
1627           px = GENSTK.plr[0] * GENSTK.cxr[0];
1628           py = GENSTK.plr[0] * GENSTK.cyr[0];
1629           pz = GENSTK.plr[0] * GENSTK.czr[0];
1630           e  = GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6];
1631       } else {
1632           Double_t p = sqrt(TRACKR.etrack * TRACKR.etrack - ParticleMassFPC(TRACKR.jtrack) *  ParticleMassFPC(TRACKR.jtrack));
1633           px = p*TRACKR.cxtrck;
1634           py = p*TRACKR.cytrck;
1635           pz = p*TRACKR.cztrck;
1636           e  = TRACKR.etrack;
1637       }
1638   }
1639   else
1640       Warning("TrackMomentum","momentum not available");
1641 }
1642
1643 //______________________________________________________________________________
1644 Double_t TFluka::TrackStep() const
1645 {
1646 // Return the length in centimeters of the current step
1647 // TRACKR.ctrack = total curved path
1648     FlukaCallerCode_t caller = GetCaller();
1649     if (caller == kMGDRAW) {
1650         Int_t i;
1651         if ((i = fPrimaryElectronIndex) > -1) {
1652             if (i > 0) {
1653                 return (fPIlength[i] - fPIlength[i-1]);
1654             } else {
1655                 Double_t s (TRACKR.ctrack - (fPIlength[fNPI - 1] - fPIlength[0]));
1656                 return s;
1657             }
1658         } else {
1659             return TRACKR.ctrack;
1660         }
1661     } else if (caller == kBXEntering || caller == kBXExiting ||
1662                caller == kENDRAW     || caller == kUSDRAW ||
1663                caller == kUSTCKV     || caller == kMGResumedTrack ||
1664                caller == kSODRAW)
1665     {
1666         return 0.0;
1667     } else {
1668         Warning("TrackStep", "track step not available");
1669         return 0.0;
1670     }
1671 }
1672
1673 //______________________________________________________________________________
1674 Double_t TFluka::TrackLength() const
1675 {
1676 // TRACKR.cmtrck = cumulative curved path since particle birth
1677   FlukaCallerCode_t caller = GetCaller();
1678   if (caller == kMGDRAW) {
1679       Int_t i;
1680       if ((i = fPrimaryElectronIndex) > -1) {
1681           return fPIlength[i];
1682       } else {
1683           return TRACKR.cmtrck;
1684       }
1685
1686   } else if (caller == kBXEntering || caller == kBXExiting ||
1687              caller == kENDRAW || caller == kUSDRAW || caller == kUSTCKV)
1688       return TRACKR.cmtrck;
1689   else if (caller == kMGResumedTrack)
1690       return TRACKR.spausr[8];
1691   else if (caller == kSODRAW)
1692       return 0.0;
1693   else {
1694       Warning("TrackLength", "track length not available for caller %5d \n", caller);
1695       return 0.0;
1696   }
1697 }
1698
1699
1700 //______________________________________________________________________________
1701 Double_t TFluka::TrackTime() const
1702 {
1703 // Return the current time of flight of the track being transported
1704 // TRACKR.atrack = age of the particle
1705   FlukaCallerCode_t caller = GetCaller();
1706   if (caller == kMGDRAW) {
1707       Int_t i;
1708       if ((i = fPrimaryElectronIndex) > -1) {
1709           Double_t t = fPItime[i];
1710           return t;
1711       } else {
1712           return TRACKR.atrack;
1713       }
1714   } else if (caller == kBXEntering || caller == kBXExiting ||
1715              caller == kENDRAW     || caller == kUSDRAW    ||
1716              caller == kUSTCKV)
1717     return TRACKR.atrack;
1718   else if (caller == kMGResumedTrack)
1719     return TRACKR.spausr[3];
1720   else if (caller == kSODRAW) {
1721       return (FLKSTK.agestk[FLKSTK.npflka]);
1722   }
1723   else {
1724     Warning("TrackTime", "track time not available");
1725     return 0.0;
1726   }
1727 }
1728
1729 //______________________________________________________________________________
1730 Double_t TFluka::Edep() const
1731 {
1732 // Energy deposition
1733 // if TRACKR.ntrack = 0, TRACKR.mtrack = 0:
1734 // -->local energy deposition (the value and the point are not recorded in TRACKR)
1735 //    but in the variable "rull" of the procedure "endraw.cxx"
1736 // if TRACKR.ntrack > 0, TRACKR.mtrack = 0:
1737 // -->no energy loss along the track
1738 // if TRACKR.ntrack > 0, TRACKR.mtrack > 0:
1739 // -->energy loss distributed along the track
1740 // TRACKR.dtrack = energy deposition of the jth deposition event
1741
1742   // If coming from bxdraw we have 2 steps of 0 length and 0 edep
1743   // If coming from usdraw we just signal particle production - no edep
1744   // If just first time after resuming, no edep for the primary
1745   FlukaCallerCode_t caller = GetCaller();
1746
1747   if (caller == kBXExiting || caller == kBXEntering ||
1748       caller == kUSDRAW    || caller == kMGResumedTrack ||
1749       caller == kSODRAW)
1750       return 0.0;
1751   Double_t sum = 0;
1752   Int_t i = -1;
1753
1754   // Material with primary ionisation activated but number of primary electrons nprim = 0
1755   if (fPrimaryElectronIndex == -2) return 0.0;
1756   // nprim > 0
1757   if ((i = fPrimaryElectronIndex) > -1) {
1758       // Primary ionisation
1759       sum = GetPrimaryElectronKineticEnergy(i);
1760       if (sum > 100.) {
1761           printf("edep > 100. %d %d %f \n", i, ALLDLT.nalldl, sum);
1762       }
1763       return sum;
1764   } else {
1765       // Normal ionisation
1766       if (TRACKR.mtrack > 1) printf("Edep: %6d\n", TRACKR.mtrack);
1767
1768       for ( Int_t j=0;j<TRACKR.mtrack;j++) {
1769           sum +=TRACKR.dtrack[j];
1770       }
1771       if (TRACKR.ntrack == 0 && TRACKR.mtrack == 0)
1772           return fRull + sum;
1773       else {
1774           return sum;
1775       }
1776   }
1777 }
1778
1779 //______________________________________________________________________________
1780 Int_t TFluka::CorrectFlukaId() const
1781 {
1782    // since we don't put photons and e- created bellow transport cut on the vmc stack
1783    // and there is a call to endraw for energy deposition for each of them
1784    // and they have the track number of their parent, but different identity (pdg)
1785    // so we want to assign also their parent identity.
1786
1787    if( (IsTrackStop())
1788         && TRACKR.ispusr[mkbmx2 - 4] == TRACKR.ispusr[mkbmx2 - 1]
1789         && TRACKR.jtrack != TRACKR.ispusr[mkbmx2 - 3] ) {
1790       if (fVerbosityLevel >=3)
1791          cout << "CorrectFlukaId() for icode=" << GetIcode()
1792                << " track=" << TRACKR.ispusr[mkbmx2 - 1]
1793                << " current PDG=" << PDGFromId(TRACKR.jtrack)
1794                << " assign parent PDG=" << PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
1795       return TRACKR.ispusr[mkbmx2 - 3]; // assign parent identity
1796    }
1797    if (TRACKR.jtrack <= 64){
1798        return TRACKR.jtrack;
1799    } else {
1800        return TRACKR.j0trck;
1801    }
1802 }
1803
1804
1805 //______________________________________________________________________________
1806 Int_t TFluka::TrackPid() const
1807 {
1808 // Return the id of the particle transported
1809 // TRACKR.jtrack = identity number of the particle
1810   FlukaCallerCode_t caller = GetCaller();
1811   if (caller != kEEDRAW && caller != kSODRAW) {
1812      return PDGFromId( CorrectFlukaId() );
1813   }
1814   else if (caller == kSODRAW) {
1815       return PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
1816   }
1817   else
1818     return -1000;
1819 }
1820
1821 //______________________________________________________________________________
1822 Double_t TFluka::TrackCharge() const
1823 {
1824 // Return charge of the track currently transported
1825 // PAPROP.ichrge = electric charge of the particle
1826 // TRACKR.jtrack = identity number of the particle
1827
1828   FlukaCallerCode_t caller = GetCaller();
1829   if (caller != kEEDRAW && caller != kSODRAW)
1830      return PAPROP.ichrge[CorrectFlukaId() + 6];
1831   else if (caller == kSODRAW) {
1832       Int_t ifl =  PDGFromId(FLKSTK.iloflk[FLKSTK.npflka]);
1833       return PAPROP.ichrge[ifl + 6];
1834   }
1835   else
1836     return -1000.0;
1837 }
1838
1839 //______________________________________________________________________________
1840 Double_t TFluka::TrackMass() const
1841 {
1842 // PAPROP.am = particle mass in GeV
1843 // TRACKR.jtrack = identity number of the particle
1844   FlukaCallerCode_t caller = GetCaller();
1845   if (caller != kEEDRAW && caller != kSODRAW)
1846      return PAPROP.am[CorrectFlukaId()+6];
1847   else if (caller == kSODRAW) {
1848       Int_t ifl =  FLKSTK.iloflk[FLKSTK.npflka];
1849       return PAPROP.am[ifl + 6];
1850   }
1851   else
1852     return -1000.0;
1853 }
1854
1855 //______________________________________________________________________________
1856 Double_t TFluka::Etot() const
1857 {
1858 // TRACKR.etrack = total energy of the particle
1859   FlukaCallerCode_t  caller = GetCaller();
1860   FlukaProcessCode_t icode  = GetIcode();
1861   if (caller != kEEDRAW && caller != kSODRAW && caller != kUSDRAW)
1862   {
1863       return TRACKR.etrack;
1864   } else if (caller == kUSDRAW) {
1865       if (icode == kEMFSCObrems  ||
1866           icode == kEMFSCOmoller ||
1867           icode == kEMFSCObhabha ||
1868           icode == kEMFSCOcompton ) {
1869           return  fPint[3];
1870       }
1871       else if (icode == kKASKADdray  ||
1872                icode == kKASKADbrems ||
1873                icode == kKASKADpair) {
1874           return (GENSTK.tki[0] + PAPROP.am[GENSTK.kpart[0]+6]);
1875       } else {
1876           return TRACKR.etrack;
1877       }
1878
1879   }
1880   else if (caller == kSODRAW) {
1881       Int_t ist  = FLKSTK.npflka;
1882       Double_t p = FLKSTK.pmoflk[ist];
1883       Int_t ifl  = FLKSTK.iloflk[ist];
1884       Double_t m = PAPROP.am[ifl + 6];
1885       Double_t e = TMath::Sqrt(p * p + m * m);
1886       return e;
1887   }
1888   printf("Etot %5d %5d \n", caller, icode);
1889
1890   return -1000.0;
1891 }
1892
1893 //
1894 // track status
1895 //
1896 //______________________________________________________________________________
1897 Bool_t   TFluka::IsNewTrack() const
1898 {
1899 // Return true for the first call of Stepping()
1900    return fTrackIsNew;
1901 }
1902
1903 void     TFluka::SetTrackIsNew(Bool_t flag)
1904 {
1905 // Return true for the first call of Stepping()
1906    fTrackIsNew = flag;
1907
1908 }
1909
1910
1911 //______________________________________________________________________________
1912 Bool_t   TFluka::IsTrackInside() const
1913 {
1914 // True if the track is not at the boundary of the current volume
1915 // In Fluka a step is always inside one kind of material
1916 // If the step would go behind the region of one material,
1917 // it will be shortened to reach only the boundary.
1918 // Therefore IsTrackInside() is always true.
1919   FlukaCallerCode_t caller = GetCaller();
1920   if (caller == kBXEntering || caller == kBXExiting)
1921     return 0;
1922   else
1923     return 1;
1924 }
1925
1926 //______________________________________________________________________________
1927 Bool_t   TFluka::IsTrackEntering() const
1928 {
1929 // True if this is the first step of the track in the current volume
1930
1931   FlukaCallerCode_t caller = GetCaller();
1932   if (caller == kBXEntering)
1933     return 1;
1934   else return 0;
1935 }
1936
1937 //______________________________________________________________________________
1938 Bool_t   TFluka::IsTrackExiting() const
1939 {
1940 // True if track is exiting volume
1941 //
1942   FlukaCallerCode_t caller = GetCaller();
1943   if (caller == kBXExiting)
1944     return 1;
1945   else return 0;
1946 }
1947
1948 //______________________________________________________________________________
1949 Bool_t   TFluka::IsTrackOut() const
1950 {
1951 // True if the track is out of the setup
1952 // means escape
1953   FlukaProcessCode_t icode = GetIcode();
1954
1955   if (icode == kKASKADescape ||
1956       icode == kEMFSCOescape ||
1957       icode == kKASNEUescape ||
1958       icode == kKASHEAescape ||
1959       icode == kKASOPHescape)
1960        return 1;
1961   else return 0;
1962 }
1963
1964 //______________________________________________________________________________
1965 Bool_t   TFluka::IsTrackDisappeared() const
1966 {
1967 // All inelastic interactions and decays
1968 // fIcode from usdraw
1969   FlukaProcessCode_t icode = GetIcode();
1970   if (icode == kKASKADinelint    || // inelastic interaction
1971       icode == kKASKADdecay      || // particle decay
1972       icode == kKASKADdray       || // delta ray generation by hadron
1973       icode == kKASKADpair       || // direct pair production
1974       icode == kKASKADbrems      || // bremsstrahlung (muon)
1975       icode == kEMFSCObrems      || // bremsstrahlung (electron)
1976       icode == kEMFSCOmoller     || // Moller scattering
1977       icode == kEMFSCObhabha     || // Bhaba scattering
1978       icode == kEMFSCOanniflight || // in-flight annihilation
1979       icode == kEMFSCOannirest   || // annihilation at rest
1980       icode == kEMFSCOpair       || // pair production
1981       icode == kEMFSCOcompton    || // Compton scattering
1982       icode == kEMFSCOphotoel    || // Photoelectric effect
1983       icode == kKASNEUhadronic   || // hadronic interaction
1984       icode == kKASHEAdray          // delta-ray
1985       ) return 1;
1986   else return 0;
1987 }
1988
1989 //______________________________________________________________________________
1990 Bool_t   TFluka::IsTrackStop() const
1991 {
1992 // True if the track energy has fallen below the threshold
1993 // means stopped by signal or below energy threshold
1994   FlukaProcessCode_t icode = GetIcode();
1995   if (icode == kKASKADstopping  || // stopping particle
1996       icode == kKASKADtimekill  || // time kill
1997       icode == kEMFSCOstopping1 || // below user-defined cut-off
1998       icode == kEMFSCOstopping2 || // below user cut-off
1999       icode == kEMFSCOtimekill  || // time kill
2000       icode == kKASNEUstopping  || // neutron below threshold
2001       icode == kKASNEUtimekill  || // time kill
2002       icode == kKASHEAtimekill  || // time kill
2003       icode == kKASOPHtimekill) return 1; // time kill
2004   else return 0;
2005 }
2006
2007 //______________________________________________________________________________
2008 Bool_t   TFluka::IsTrackAlive() const
2009 {
2010 // Means not disappeared or not out
2011     FlukaProcessCode_t icode = GetIcode();
2012
2013     if (IsTrackOut()               ||
2014         IsTrackStop()              ||
2015         icode == kKASKADinelint    || // inelastic interaction
2016         icode == kKASKADdecay      || // particle decay
2017         icode == kEMFSCOanniflight || // in-flight annihilation
2018         icode == kEMFSCOannirest   || // annihilation at rest
2019         icode == kEMFSCOpair       || // pair production
2020         icode == kEMFSCOphotoel    || // Photoelectric effect
2021         icode == kKASNEUhadronic      // hadronic interaction
2022         )
2023     {
2024         // Exclude the cases for which the particle has disappeared (paused) but will reappear later (= alive).
2025         return 0;
2026     } else {
2027         return 1;
2028     }
2029 }
2030
2031 //
2032 // secondaries
2033 //
2034
2035 //______________________________________________________________________________
2036 Int_t TFluka::NSecondaries() const
2037
2038 {
2039 // Number of secondary particles generated in the current step
2040 // GENSTK.np = number of secondaries except light and heavy ions
2041 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
2042     FlukaCallerCode_t caller = GetCaller();
2043     if (caller == kUSDRAW)  // valid only after usdraw
2044         return GENSTK.np + FHEAVY.npheav;
2045     else if (caller == kUSTCKV) {
2046         // Cerenkov Photon production
2047         return fNCerenkov;
2048     }
2049     return 0;
2050 } // end of NSecondaries
2051
2052 //______________________________________________________________________________
2053 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
2054                 TLorentzVector& position, TLorentzVector& momentum)
2055 {
2056 // Copy particles from secondary stack to vmc stack
2057 //
2058
2059     FlukaCallerCode_t caller = GetCaller();
2060     if (caller == kUSDRAW) {  // valid only after usdraw
2061         if (GENSTK.np > 0) {
2062             // Hadronic interaction
2063             if (isec >= 0 && isec < GENSTK.np) {
2064                 particleId = PDGFromId(GENSTK.kpart[isec]);
2065                 position.SetX(fXsco);
2066                 position.SetY(fYsco);
2067                 position.SetZ(fZsco);
2068                 position.SetT(TRACKR.atrack);
2069                 momentum.SetPx(GENSTK.plr[isec]*GENSTK.cxr[isec]);
2070                 momentum.SetPy(GENSTK.plr[isec]*GENSTK.cyr[isec]);
2071                 momentum.SetPz(GENSTK.plr[isec]*GENSTK.czr[isec]);
2072                 momentum.SetE(GENSTK.tki[isec] + PAPROP.am[GENSTK.kpart[isec]+6]);
2073             }
2074             else if (isec >= GENSTK.np && isec < GENSTK.np + FHEAVY.npheav) {
2075                 Int_t jsec = isec - GENSTK.np;
2076                 particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
2077                 position.SetX(fXsco);
2078                 position.SetY(fYsco);
2079                 position.SetZ(fZsco);
2080                 position.SetT(TRACKR.atrack);
2081                 momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
2082                 momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
2083                 momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
2084                 if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
2085                     momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
2086                 else if (FHEAVY.tkheav[jsec] > 6)
2087                     momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
2088             }
2089             else
2090                 Warning("GetSecondary","isec out of range");
2091         }
2092     } else if (caller == kUSTCKV) {
2093         Int_t index = OPPHST.lstopp - isec;
2094         position.SetX(OPPHST.xoptph[index]);
2095         position.SetY(OPPHST.yoptph[index]);
2096         position.SetZ(OPPHST.zoptph[index]);
2097         position.SetT(OPPHST.agopph[index]);
2098         Double_t p = OPPHST.poptph[index];
2099
2100         momentum.SetPx(p * OPPHST.txopph[index]);
2101         momentum.SetPy(p * OPPHST.tyopph[index]);
2102         momentum.SetPz(p * OPPHST.tzopph[index]);
2103         momentum.SetE(p);
2104     }
2105     else
2106         Warning("GetSecondary","no secondaries available");
2107
2108 } // end of GetSecondary
2109
2110
2111 //______________________________________________________________________________
2112 TMCProcess TFluka::ProdProcess(Int_t) const
2113
2114 {
2115 // Name of the process that has produced the secondary particles
2116 // in the current step
2117
2118     Int_t mugamma = (TRACKR.jtrack == kFLUKAphoton ||
2119                      TRACKR.jtrack == kFLUKAmuplus ||
2120                      TRACKR.jtrack == kFLUKAmuminus);
2121     FlukaProcessCode_t icode = GetIcode();
2122
2123     if      (icode == kKASKADdecay)                                   return kPDecay;
2124     else if (icode == kKASKADpair || icode == kEMFSCOpair)            return kPPair;
2125     else if (icode == kEMFSCOcompton)                                 return kPCompton;
2126     else if (icode == kEMFSCOphotoel)                                 return kPPhotoelectric;
2127     else if (icode == kKASKADbrems      || icode == kEMFSCObrems)     return kPBrem;
2128     else if (icode == kKASKADdray       || icode == kKASHEAdray)      return kPDeltaRay;
2129     else if (icode == kEMFSCOmoller     || icode == kEMFSCObhabha)    return kPDeltaRay;
2130     else if (icode == kEMFSCOanniflight || icode == kEMFSCOannirest)  return kPAnnihilation;
2131     else if (icode == kKASKADinelint) {
2132         if (!mugamma)                                                 return kPHadronic;
2133         else if (TRACKR.jtrack == kFLUKAphoton)                       return kPPhotoFission;
2134         else                                                          return kPMuonNuclear;
2135     }
2136     else if (icode == kEMFSCOrayleigh)                                return kPRayleigh;
2137 // Fluka codes 100, 300 and 400 still to be investigasted
2138     else                                                              return kPNoProcess;
2139 }
2140
2141
2142 Int_t TFluka::StepProcesses(TArrayI &proc) const
2143 {
2144   //
2145   // Return processes active in the current step
2146   //
2147     FlukaProcessCode_t icode   = GetIcode();
2148     FlukaCallerCode_t  caller  = GetCaller();
2149     proc.Set(1);
2150     TMCProcess iproc;
2151     if (caller == kBXEntering || caller == kBXExiting || caller == kEEDRAW || caller == kSODRAW) {
2152         iproc = kPTransportation;
2153     }
2154     else if (caller == kUSTCKV) {
2155         iproc = kPCerenkov;
2156     } else {
2157         switch (icode) {
2158         case kEMFSCO:
2159             if (Edep() > 0.) {
2160                 iproc = kPEnergyLoss;
2161             } else {
2162                 iproc = kPTransportation;
2163             }
2164             break;
2165         case kKASKAD:
2166             if (Edep() > 0.) {
2167                 iproc = kPEnergyLoss;
2168             } else {
2169                 iproc = kPTransportation;
2170             }
2171             break;
2172         case kKASHEA:
2173         case kKASNEU:
2174         case kKASOPH:
2175         case kKASKADescape:
2176         case kEMFSCOescape:
2177         case kKASNEUescape:
2178         case kKASHEAescape:
2179         case kKASOPHescape:
2180             iproc = kPTransportation;
2181             break;
2182         case kKASKADtimekill:
2183         case kEMFSCOtimekill:
2184         case kKASNEUtimekill:
2185         case kKASHEAtimekill:
2186         case kKASOPHtimekill:
2187             iproc =  kPTOFlimit;
2188             break;
2189         case kKASKADstopping:
2190         case kEMFSCOstopping1:
2191         case kEMFSCOstopping2:
2192         case kKASNEUstopping:
2193             iproc = kPStop;
2194             break;
2195         case kKASKADinelint:
2196         case kKASNEUhadronic:
2197             iproc = kPHadronic;
2198             break;
2199         case kKASKADinelarecoil:
2200             iproc = kPHadronic;
2201             break;
2202         case kKASKADnelint:
2203             iproc = kPHElastic;
2204             break;
2205         case kKASOPHabsorption:
2206             iproc = kPLightAbsorption;
2207             break;
2208         case kKASOPHrefraction:
2209             iproc = kPLightRefraction;
2210             break;
2211         case kEMFSCOlocaldep :
2212             iproc = kPPhotoelectric;
2213             break;
2214         default:
2215             iproc = ProdProcess(0);
2216         }
2217     }
2218
2219     proc[0] = iproc;
2220     return 1;
2221 }
2222 //______________________________________________________________________________
2223 Int_t TFluka::VolId2Mate(Int_t id) const
2224 {
2225 //
2226 // Returns the material number for a given volume ID
2227 //
2228    return fMCGeo->VolId2Mate(id);
2229 }
2230
2231 //______________________________________________________________________________
2232 const char* TFluka::VolName(Int_t id) const
2233 {
2234 //
2235 // Returns the volume name for a given volume ID
2236 //
2237    return fMCGeo->VolName(id);
2238 }
2239
2240 Int_t TFluka::MediumId(const Text_t* mediumName) const
2241 {
2242     //
2243     // Return the unique medium id for medium with name mediumName
2244     TList *medlist = gGeoManager->GetListOfMedia();
2245     TGeoMedium* med = (TGeoMedium*) medlist->FindObject(mediumName);
2246     if (med) {
2247         return (med->GetId());
2248     } else {
2249         return (-1);
2250     }
2251 }
2252
2253 //______________________________________________________________________________
2254 Int_t TFluka::VolId(const Text_t* volName) const
2255 {
2256 //
2257 // Converts from volume name to volume ID.
2258 // Time consuming. (Only used during set-up)
2259 // Could be replaced by hash-table
2260 //
2261     char sname[20];
2262     Int_t len;
2263     strncpy(sname, volName, len = strlen(volName));
2264     sname[len] = 0;
2265     while (sname[len - 1] == ' ') sname[--len] = 0;
2266     return fMCGeo->VolId(sname);
2267 }
2268
2269 //______________________________________________________________________________
2270 Int_t TFluka::CurrentVolID(Int_t& copyNo) const
2271 {
2272 //
2273 // Return the logical id and copy number corresponding to the current fluka region
2274 //
2275   if (gGeoManager->IsOutside()) return 0;
2276   TGeoNode *node = gGeoManager->GetCurrentNode();
2277   copyNo = node->GetNumber();
2278   Int_t id = node->GetVolume()->GetNumber();
2279   return id;
2280 }
2281
2282 //______________________________________________________________________________
2283 Int_t TFluka::CurrentVolOffID(Int_t off, Int_t& copyNo) const
2284 {
2285 //
2286 // Return the logical id and copy number of off'th mother
2287 // corresponding to the current fluka region
2288 //
2289   if (off<0 || off>gGeoManager->GetLevel()) return 0;
2290   if (off==0) return CurrentVolID(copyNo);
2291   TGeoNode *node = gGeoManager->GetMother(off);
2292   if (!node) return 0;
2293   copyNo = node->GetNumber();
2294   return node->GetVolume()->GetNumber();
2295 }
2296
2297 //______________________________________________________________________________
2298 const char* TFluka::CurrentVolName() const
2299 {
2300 //
2301 // Return the current volume name
2302 //
2303   if (gGeoManager->IsOutside()) return "OutOfWorld";
2304   return gGeoManager->GetCurrentVolume()->GetName();
2305 }
2306
2307 //______________________________________________________________________________
2308 const char* TFluka::CurrentVolOffName(Int_t off) const
2309 {
2310 //
2311 // Return the volume name of the off'th mother of the current volume
2312 //
2313   if (off<0 || off>gGeoManager->GetLevel()) return 0;
2314   if (off==0) return CurrentVolName();
2315   TGeoNode *node = gGeoManager->GetMother(off);
2316   if (!node) return 0;
2317   return node->GetVolume()->GetName();
2318 }
2319
2320 const char* TFluka::CurrentVolPath() {
2321   // Return the current volume path
2322   return gGeoManager->GetPath();
2323 }
2324 //______________________________________________________________________________
2325 Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z,
2326                       Float_t & dens, Float_t & radl, Float_t & absl) const
2327 {
2328 //
2329 //  Return the current medium number and material properties
2330 //
2331   Int_t copy;
2332   Int_t id  =  TFluka::CurrentVolID(copy);
2333   Int_t med =  TFluka::VolId2Mate(id);
2334   TGeoVolume*     vol = gGeoManager->GetCurrentVolume();
2335   TGeoMaterial*   mat = vol->GetMaterial();
2336   a    = mat->GetA();
2337   z    = mat->GetZ();
2338   dens = mat->GetDensity();
2339   radl = mat->GetRadLen();
2340   absl = mat->GetIntLen();
2341
2342   return med;
2343 }
2344
2345 //______________________________________________________________________________
2346 void TFluka::Gmtod(Float_t* xm, Float_t* xd, Int_t iflag)
2347 {
2348 // Transforms a position from the world reference frame
2349 // to the current volume reference frame.
2350 //
2351 //  Geant3 desription:
2352 //  ==================
2353 //       Computes coordinates XD (in DRS)
2354 //       from known coordinates XM in MRS
2355 //       The local reference system can be initialized by
2356 //         - the tracking routines and GMTOD used in GUSTEP
2357 //         - a call to GMEDIA(XM,NUMED)
2358 //         - a call to GLVOLU(NLEVEL,NAMES,NUMBER,IER)
2359 //             (inverse routine is GDTOM)
2360 //
2361 //        If IFLAG=1  convert coordinates
2362 //           IFLAG=2  convert direction cosinus
2363 //
2364 // ---
2365    Double_t xmL[3], xdL[3];
2366    Int_t i;
2367    for (i=0;i<3;i++) xmL[i]=xm[i];
2368    if (iflag == 1) gGeoManager->MasterToLocal(xmL,xdL);
2369    else            gGeoManager->MasterToLocalVect(xmL,xdL);
2370    for (i=0;i<3;i++) xd[i] = xdL[i];
2371 }
2372
2373 //______________________________________________________________________________
2374 void TFluka::Gmtod(Double_t* xm, Double_t* xd, Int_t iflag)
2375 {
2376 //
2377 // See Gmtod(Float_t*, Float_t*, Int_t)
2378 //
2379    if (iflag == 1) gGeoManager->MasterToLocal(xm,xd);
2380    else            gGeoManager->MasterToLocalVect(xm,xd);
2381 }
2382
2383 //______________________________________________________________________________
2384 void TFluka::Gdtom(Float_t* xd, Float_t* xm, Int_t iflag)
2385 {
2386 // Transforms a position from the current volume reference frame
2387 // to the world reference frame.
2388 //
2389 //  Geant3 desription:
2390 //  ==================
2391 //  Computes coordinates XM (Master Reference System
2392 //  knowing the coordinates XD (Detector Ref System)
2393 //  The local reference system can be initialized by
2394 //    - the tracking routines and GDTOM used in GUSTEP
2395 //    - a call to GSCMED(NLEVEL,NAMES,NUMBER)
2396 //        (inverse routine is GMTOD)
2397 //
2398 //   If IFLAG=1  convert coordinates
2399 //      IFLAG=2  convert direction cosinus
2400 //
2401 // ---
2402    Double_t xmL[3], xdL[3];
2403    Int_t i;
2404    for (i=0;i<3;i++) xdL[i] = xd[i];
2405    if (iflag == 1) gGeoManager->LocalToMaster(xdL,xmL);
2406    else            gGeoManager->LocalToMasterVect(xdL,xmL);
2407    for (i=0;i<3;i++) xm[i]=xmL[i];
2408 }
2409
2410 //______________________________________________________________________________
2411 void TFluka::Gdtom(Double_t* xd, Double_t* xm, Int_t iflag)
2412 {
2413 //
2414 // See Gdtom(Float_t*, Float_t*, Int_t)
2415 //
2416    if (iflag == 1) gGeoManager->LocalToMaster(xd,xm);
2417    else            gGeoManager->LocalToMasterVect(xd,xm);
2418 }
2419
2420 //______________________________________________________________________________
2421 TObjArray *TFluka::GetFlukaMaterials()
2422 {
2423 //
2424 // Get array of Fluka materials
2425    return fGeom->GetMatList();
2426 }
2427
2428 //______________________________________________________________________________
2429 void TFluka::SetMreg(Int_t l, Int_t lttc)
2430 {
2431 // Set current fluka region
2432    fCurrentFlukaRegion = l;
2433    fGeom->SetMreg(l,lttc);
2434 }
2435
2436
2437
2438
2439 //______________________________________________________________________________
2440 TString TFluka::ParticleName(Int_t pdg) const
2441 {
2442     // Return particle name for particle with pdg code pdg.
2443     Int_t ifluka = IdFromPDG(pdg);
2444     return TString((CHPPRP.btype[ifluka - kFLUKAcodemin]), 8);
2445 }
2446
2447
2448 //______________________________________________________________________________
2449 Double_t TFluka::ParticleMass(Int_t pdg) const
2450 {
2451     // Return particle mass for particle with pdg code pdg.
2452     Int_t ifluka = IdFromPDG(pdg);
2453     return (PAPROP.am[ifluka - kFLUKAcodemin]);
2454 }
2455
2456 //______________________________________________________________________________
2457 Double_t TFluka::ParticleMassFPC(Int_t fpc) const
2458 {
2459     // Return particle mass for particle with Fluka particle code fpc
2460     return (PAPROP.am[fpc - kFLUKAcodemin]);
2461 }
2462
2463 //______________________________________________________________________________
2464 Double_t TFluka::ParticleCharge(Int_t pdg) const
2465 {
2466     // Return particle charge for particle with pdg code pdg.
2467     Int_t ifluka = IdFromPDG(pdg);
2468     return Double_t(PAPROP.ichrge[ifluka - kFLUKAcodemin]);
2469 }
2470
2471 //______________________________________________________________________________
2472 Double_t TFluka::ParticleLifeTime(Int_t pdg) const
2473 {
2474     // Return particle lifetime for particle with pdg code pdg.
2475     Int_t ifluka = IdFromPDG(pdg);
2476     return (PAPROP.tmnlf[ifluka - kFLUKAcodemin]);
2477 }
2478
2479 //______________________________________________________________________________
2480 void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
2481 {
2482     // Retrieve particle properties for particle with pdg code pdg.
2483
2484     strcpy(name, ParticleName(pdg).Data());
2485     type   = ParticleMCType(pdg);
2486     mass   = ParticleMass(pdg);
2487     charge = ParticleCharge(pdg);
2488     tlife  = ParticleLifeTime(pdg);
2489 }
2490
2491 //______________________________________________________________________________
2492 void TFluka::PrintHeader()
2493 {
2494     //
2495     // Print a header
2496     printf("\n");
2497     printf("\n");
2498     printf("------------------------------------------------------------------------------\n");
2499     printf("- You are using the TFluka Virtual Monte Carlo Interface to FLUKA.           -\n");
2500     printf("- Please see the file fluka.out for FLUKA output and licensing information.  -\n");
2501     printf("------------------------------------------------------------------------------\n");
2502     printf("\n");
2503     printf("\n");
2504 }
2505
2506
2507 #define pshckp pshckp_
2508 #define ustckv ustckv_
2509
2510
2511 extern "C" {
2512   void pshckp(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
2513               Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
2514               Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
2515   {
2516     //
2517     // Pushes one cerenkov photon to the stack
2518     //
2519
2520     TFluka* fluka =  (TFluka*) gMC;
2521     TVirtualMCStack* cppstack = fluka->GetStack();
2522     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
2523     cppstack->PushTrack(0, parent, 50000050,
2524                         px, py, pz, e,
2525                         vx, vy, vz, tof,
2526                         polx, poly, polz,
2527                         kPCerenkov, ntr, wgt, 0);
2528     if (fluka->GetVerbosityLevel() >= 3)
2529             printf("pshckp: track=%d parent=%d lattc=%d %s\n", ntr, parent, TRACKR.lt1trk, fluka->CurrentVolName());
2530   }
2531
2532     void ustckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
2533     {
2534         //
2535         // Calls stepping in order to signal cerenkov production
2536         //
2537         TFluka *fluka = (TFluka*)gMC;
2538         fluka->SetMreg(mreg, TRACKR.lt1trk); //LTCLCM.mlatm1);
2539         fluka->SetXsco(x);
2540         fluka->SetYsco(y);
2541         fluka->SetZsco(z);
2542         fluka->SetNCerenkov(nphot);
2543         fluka->SetCaller(kUSTCKV);
2544         if (fluka->GetVerbosityLevel() >= 3)
2545             printf("ustckv: %10d mreg=%d lattc=%d  newlat=%d (%f, %f, %f) edep=%f vol=%s\n",
2546                     nphot, mreg, TRACKR.lt1trk, LTCLCM.newlat, x, y, z, fluka->Edep(), fluka->CurrentVolName());
2547
2548     // check region lattice consistency (debug Ernesto)
2549     // *****************************************************
2550    Int_t nodeId;
2551    Int_t volId = fluka->CurrentVolID(nodeId);
2552    Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
2553
2554    if( mreg != volId  && !gGeoManager->IsOutside() ) {
2555        cout << "  ustckv:   track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
2556             << " icode=" << fluka->GetIcode() << " gNstep=" << fluka->GetNstep() << endl
2557             << "               fluka   mreg=" << mreg << " mlttc=" << TRACKR.lt1trk << endl
2558             << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
2559             << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
2560             << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
2561             << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
2562             << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
2563         if( TRACKR.lt1trk == crtlttc ) cout << "   *************************************************************" << endl;
2564     }
2565     // *****************************************************
2566
2567
2568
2569         (TVirtualMCApplication::Instance())->Stepping();
2570     }
2571 }
2572
2573 //______________________________________________________________________________
2574 void TFluka::AddParticlesToPdgDataBase() const
2575 {
2576
2577 //
2578 // Add particles to the PDG data base
2579
2580     TDatabasePDG *pdgDB = TDatabasePDG::Instance();
2581
2582     const Double_t kAu2Gev   = 0.9314943228;
2583     const Double_t khSlash   = 1.0545726663e-27;
2584     const Double_t kErg2Gev  = 1/1.6021773349e-3;
2585     const Double_t khShGev   = khSlash*kErg2Gev;
2586     const Double_t kYear2Sec = 3600*24*365.25;
2587 //
2588 // Ions
2589 //
2590   pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE,
2591                      0,3,"Ion",GetIonPdg(1,2));
2592   pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE,
2593                      khShGev/(12.33*kYear2Sec),3,"Ion",GetIonPdg(1,3));
2594   pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE,
2595                      khShGev/(12.33*kYear2Sec),6,"Ion",GetIonPdg(2,4));
2596   pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE,
2597                      0,6,"Ion",GetIonPdg(2,3));
2598 //
2599 //
2600 //
2601 // Special particles
2602 //
2603   pdgDB->AddParticle("Cherenkov","Cherenkov",0,kFALSE,
2604                      0,0,"Special",GetSpecialPdg(50));
2605   pdgDB->AddParticle("FeedbackPhoton","FeedbackPhoton",0,kFALSE,
2606                      0,0,"Special",GetSpecialPdg(51));
2607 }
2608
2609 void TFluka::AddIon(Int_t a, Int_t z) const
2610 {
2611
2612     // Add a new ion
2613     TDatabasePDG *pdgDB = TDatabasePDG::Instance();
2614     const Double_t kAu2Gev   = 0.9314943228;
2615     Int_t pdg =  GetIonPdg(z, a);
2616     if (pdgDB->GetParticle(pdg)) return;
2617
2618     pdgDB->AddParticle(Form("Iion A  = %5d Z = %5d", a, z),"Ion", Float_t(a) * kAu2Gev + 8.071e-3, kTRUE,
2619                        0, 3 * z, "Ion", pdg);
2620 }
2621
2622 //
2623 // Info about primary ionization electrons
2624 //
2625
2626 //______________________________________________________________________________
2627 Int_t TFluka::GetNPrimaryElectrons()
2628 {
2629     // Get number of primary electrons
2630     return ALLDLT.nalldl;
2631 }
2632
2633 //______________________________________________________________________________
2634 Double_t TFluka::GetPrimaryElectronKineticEnergy(Int_t i) const
2635 {
2636     // Returns kinetic energy of primary electron i
2637
2638     Double_t ekin = -1.;
2639
2640     if (i >= 0 && i < ALLDLT.nalldl) {
2641         ekin =  ALLDLT.talldl[i];
2642     } else {
2643         Warning("GetPrimaryElectronKineticEnergy",
2644                 "Primary electron index out of range %d %d \n",
2645                 i, ALLDLT.nalldl);
2646     }
2647     return ekin;
2648 }
2649
2650 void TFluka::GetPrimaryElectronPosition(Int_t i, Double_t& x, Double_t& y, Double_t& z, Double_t& t) const
2651 {
2652     // Returns position  of primary electron i
2653         if (i >= 0 && i < ALLDLT.nalldl) {
2654             x = ALLDLT.xalldl[i];
2655             y = ALLDLT.yalldl[i];
2656             z = ALLDLT.zalldl[i];
2657             t = ALLDLT.talldl[i];
2658             return;
2659         } else {
2660             Warning("GetPrimaryElectronPosition",
2661                     "Primary electron index out of range %d %d \n",
2662                     i, ALLDLT.nalldl);
2663             return;
2664         }
2665         return;
2666 }
2667
2668 Int_t TFluka::GetIonPdg(Int_t z, Int_t a, Int_t i) const
2669 {
2670 // Acording to
2671 // http://cepa.fnal.gov/psm/stdhep/pdg/montecarlorpp-2006.pdf
2672
2673   return 1000000000 + 10*1000*z + 10*a + i;
2674 }
2675
2676 //__________________________________________________________________
2677 Int_t TFluka::GetSpecialPdg(Int_t number) const
2678 {
2679 // Numbering for special particles
2680
2681   return 50000000 + number;
2682 }
2683
2684
2685 void  TFluka::PrimaryIonisationStepping(Int_t nprim)
2686 {
2687 // Call Stepping for primary ionisation electrons
2688 // Protection against nprim > mxalld
2689 // Multiple steps for nprim > 0
2690     Int_t i;
2691     fNPI = nprim;
2692     if (nprim > 0) {
2693         CalcPrimaryIonisationTime();
2694         for (i = 0; i < nprim; i++) {
2695             SetCurrentPrimaryElectronIndex(i);
2696             (TVirtualMCApplication::Instance())->Stepping();
2697             if (i == 0) SetTrackIsNew(kFALSE);
2698         }
2699     } else {
2700         // No primary electron ionisation
2701         // Call Stepping anyway but flag nprim = 0 as index = -2
2702         SetCurrentPrimaryElectronIndex(-2);
2703         (TVirtualMCApplication::Instance())->Stepping();
2704     }
2705     // Reset the index
2706     SetCurrentPrimaryElectronIndex(-1);
2707 }
2708
2709 //______________________________________________________________________
2710 Float_t* TFluka::CreateFloatArray(Double_t* array, Int_t size) const
2711 {
2712 // Converts Double_t* array to Float_t*,
2713 // !! The new array has to be deleted by user.
2714 // ---
2715
2716   Float_t* floatArray;
2717   if (size>0) {
2718     floatArray = new Float_t[size];
2719     for (Int_t i=0; i<size; i++)
2720       if (array[i] >= FLT_MAX )
2721         floatArray[i] = FLT_MAX/100.;
2722       else
2723         floatArray[i] = array[i];
2724   }
2725   else {
2726     //floatArray = 0;
2727     floatArray = new Float_t[1];
2728   }
2729   return floatArray;
2730 }
2731
2732 void TFluka::CalcPrimaryIonisationTime()
2733 {
2734     // Calculates the primary ionisation time
2735     if (fPItime) delete [] fPItime;
2736     fPItime = new Double_t[fNPI];
2737     if (fPIlength) delete [] fPIlength;
2738     fPIlength = new Double_t[fNPI];
2739     //
2740     Double_t px, py, pz, e, t;
2741     TrackMomentum(px, py, pz, e);
2742     Double_t p    = TMath::Sqrt(px * px + py * py + pz * pz);
2743     Double_t beta = p / e;
2744     Double_t x0, y0, z0;
2745     fPItime[fNPI -1]   = TRACKR.atrack;
2746     fPIlength[fNPI -1] = TRACKR.cmtrck;
2747     GetPrimaryElectronPosition(fNPI - 1, x0, y0, z0, t);
2748     if (fNPI > 1) {
2749         for (Int_t i = fNPI - 2; i > -1; i--) {
2750             Double_t x, y, z, t;
2751             GetPrimaryElectronPosition(i, x, y, z, t);
2752             Double_t ds = TMath::Sqrt((x-x0) * (x-x0) + (y-y0) * (y-y0) + (z-z0) * (z-z0));
2753             fPItime[i]   = fPItime[i+1]   - ds / (beta * 2.99792458e10);
2754             fPIlength[i] = fPIlength[i+1] - ds;
2755             x0 = x; y0 = y; z0 = z;
2756         }
2757     }
2758
2759 }