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