TFluka/TFluka.cxx

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