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