- EffC++ warnings corrected
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 21 Sep 2006 09:32:57 +0000 (09:32 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 21 Sep 2006 09:32:57 +0000 (09:32 +0000)
- Improved synchronization with fluka
(E. Lopez Torres)

18 files changed:
TFluka/TFluka.cxx
TFluka/TFluka.h
TFluka/TFlukaCerenkov.cxx
TFluka/TFlukaCerenkov.h
TFluka/TFlukaCodes.h
TFluka/TFlukaConfigOption.cxx
TFluka/TFlukaConfigOption.h
TFluka/TFlukaMCGeometry.cxx
TFluka/TFlukaMCGeometry.h
TFluka/TFlukaScoringOption.cxx
TFluka/TFlukaScoringOption.h
TFluka/bxdraw.cxx
TFluka/endraw.cxx
TFluka/mgdraw.cxx
TFluka/source.cxx
TFluka/stupre.cxx
TFluka/stuprf.cxx
TFluka/usdraw.cxx

index f583d63..3c9d7c7 100644 (file)
@@ -64,6 +64,8 @@
 #include "TArrayI.h"
 #include "TArrayD.h"
 #include "TDatabasePDG.h"
+#include "TStopwatch.h"
+
 
 // Fluka methods that may be needed.
 #ifndef WIN32 
@@ -74,6 +76,8 @@
 # define mcihad mcihad_
 # define mpdgha mpdgha_
 # define newplo newplo_
+# define genout genout_
+# define flkend flkend_
 #else 
 # define flukam  FLUKAM
 # define fluka_openinp FLUKA_OPENINP
@@ -82,6 +86,8 @@
 # define mcihad MCIHAD
 # define mpdgha MPDGHA
 # define newplo NEWPLO
+# define genout GENOUT
+# define flkend FLKEND
 #endif
 
 extern "C" 
@@ -91,6 +97,8 @@ extern "C"
   //
   void type_of_call flukam(const int&);
   void type_of_call newplo();
+  void type_of_call genout();
+  void type_of_call flkend();
   void type_of_call fluka_openinp(const int&, DEFCHARA);
   void type_of_call fluka_openout(const int&, DEFCHARA);
   void type_of_call fluka_closeinp(const int&);
@@ -110,36 +118,68 @@ ClassImp(TFluka)
 TFluka::TFluka()
   :TVirtualMC(),
    fVerbosityLevel(0),
+   fNEvent(0),
    fInputFileName(""),
+   fCoreInputFileName(""),
+   fCaller(kNoCaller),
+   fIcode(kNoProcess),
+   fNewReg(-1),
+   fRull(0),
+   fXsco(0),
+   fYsco(0),
+   fZsco(0),
+   fTrackIsEntering(kFALSE),
+   fTrackIsExiting(kFALSE),
+   fTrackIsNew(kFALSE),
+   fFieldFlag(kTRUE),
+   fGeneratePemf(kFALSE),
+   fDummyBoundary(kFALSE),
+   fStopped(kFALSE),
+   fStopEvent(kFALSE),
+   fStopRun(kFALSE),
+   fMaterials(0),
+   fNVolumes(0),
+   fCurrentFlukaRegion(-1),
+   fNCerenkov(0),
+   fGeom(0),
+   fMCGeo(0),
    fUserConfig(0), 
    fUserScore(0)
 { 
   //
   // Default constructor
   //
-   fGeneratePemf = kFALSE;
-   fNVolumes = 0;
-   fCurrentFlukaRegion = -1;
-   fNewReg = -1;
-   fGeom = 0;
-   fMCGeo = 0;
-   fMaterials = 0;
-   fDummyBoundary = 0;
-   fFieldFlag = 1;
-   fStopped   = 0;
-   fStopEvent = 0;
-   fStopRun   = 0;
-   fNEvent    = 0;
 } 
  
 //______________________________________________________________________________ 
 TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupported)
   :TVirtualMC("TFluka",title, isRootGeometrySupported),
    fVerbosityLevel(verbosity),
+   fNEvent(0),
    fInputFileName(""),
-   fTrackIsEntering(0),
-   fTrackIsExiting(0),
-   fTrackIsNew(0),
+   fCoreInputFileName(""),
+   fCaller(kNoCaller),
+   fIcode(kNoProcess),
+   fNewReg(-1),
+   fRull(0),
+   fXsco(0),
+   fYsco(0),
+   fZsco(0),
+   fTrackIsEntering(kFALSE),
+   fTrackIsExiting(kFALSE),
+   fTrackIsNew(kFALSE),
+   fFieldFlag(kTRUE),
+   fGeneratePemf(kFALSE),
+   fDummyBoundary(kFALSE),
+   fStopped(kFALSE),
+   fStopEvent(kFALSE),
+   fStopRun(kFALSE),
+   fMaterials(0),
+   fNVolumes(0),
+   fCurrentFlukaRegion(-1),
+   fNCerenkov(0),
+   fGeom(0),
+   fMCGeo(0),
    fUserConfig(new TObjArray(100)),
    fUserScore(new TObjArray(100)) 
 {
@@ -148,41 +188,31 @@ TFluka::TFluka(const char *title, Int_t verbosity, Bool_t isRootGeometrySupporte
        cout << "<== TFluka::TFluka(" << title << ") constructor called." << endl;
    SetCoreInputFileName();
    SetInputFileName();
-   SetGeneratePemf(kFALSE);
-   fNVolumes      = 0;
-   fCurrentFlukaRegion = -1;
-   fNewReg = -1;
-   fDummyBoundary = 0;
-   fFieldFlag = 1;
-   fGeneratePemf = kFALSE;
    fMCGeo = new TGeoMCGeometry("MCGeo", "TGeo Implementation of VirtualMCGeometry", kFALSE);
    fGeom  = new TFlukaMCGeometry("geom", "FLUKA VMC Geometry");
    if (verbosity > 2) fGeom->SetDebugMode(kTRUE);
-   fMaterials = 0;
-   fStopped   = 0;
-   fStopEvent = 0;
-   fStopRun   = 0;
-   fNEvent    = 0;
    PrintHeader();
 }
 
 //______________________________________________________________________________ 
-TFluka::~TFluka() {
-// Destructor
+TFluka::~TFluka()
+{
+    // Destructor
     if (fVerbosityLevel >=3)
-       cout << "<== TFluka::~TFluka() destructor called." << endl;
+        cout << "<== TFluka::~TFluka() destructor called." << endl;
+    if (fMaterials) delete [] fMaterials;
     
     delete fGeom;
     delete fMCGeo;
     
     if (fUserConfig) {
-       fUserConfig->Delete();
-       delete fUserConfig;
+        fUserConfig->Delete();
+        delete fUserConfig;
     }
     
     if (fUserScore) {
-       fUserScore->Delete();
-       delete fUserScore;
+        fUserScore->Delete();
+        delete fUserScore;
     }
 }
 
@@ -244,22 +274,22 @@ void TFluka::BuildPhysics() {
 //
     
     if (fVerbosityLevel >=3)
-       cout << "==> TFluka::BuildPhysics() called." << endl;
+        cout << "==> TFluka::BuildPhysics() called." << endl;
 
     
     if (fVerbosityLevel >=3) {
-       TList *medlist = gGeoManager->GetListOfMedia();
-       TIter next(medlist);
-       TGeoMedium*   med = 0x0;
-       TGeoMaterial* mat = 0x0;
-       Int_t ic = 0;
-       
-       while((med = (TGeoMedium*)next()))
-       {
-           mat = med->GetMaterial();
-           printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex());
-           ic++;
-       }
+        TList *medlist = gGeoManager->GetListOfMedia();
+        TIter next(medlist);
+        TGeoMedium*   med = 0x0;
+        TGeoMaterial* mat = 0x0;
+        Int_t ic = 0;
+        
+        while((med = (TGeoMedium*)next()))
+        {
+            mat = med->GetMaterial();
+            printf("Medium %5d %12s %5d %5d\n", ic, (med->GetName()), med->GetId(), mat->GetIndex());
+            ic++;
+        }
     }
     
     //
@@ -277,8 +307,13 @@ void TFluka::BuildPhysics() {
     fluka_openinp(lunin, PASSCHARA(fname));
     fluka_openout(11, PASSCHARA("fluka.out"));
 //  Read input cards    
+    cout << "==> TFluka::BuildPhysics() Read input cards." << endl;
+    TStopwatch timer;
+    timer.Start();
     GLOBAL.lfdrtr = true;
     flukam(1);
+    cout << "<== TFluka::BuildPhysics() Read input cards End"
+         << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
 //  Close input file
     fluka_closeinp(lunin);
 //  Finish geometry    
@@ -291,18 +326,18 @@ void TFluka::ProcessEvent() {
 // Process one event
 //
     if (fStopRun) {
-       Warning("ProcessEvent", "User Run Abortion: No more events handled !\n");
-       fNEvent += 1;
-       return;
+        Warning("ProcessEvent", "User Run Abortion: No more events handled !\n");
+        fNEvent += 1;
+        return;
     }
 
     if (fVerbosityLevel >=3)
-       cout << "==> TFluka::ProcessEvent() called." << endl;
+        cout << "==> TFluka::ProcessEvent() called." << endl;
     fApplication->GeneratePrimaries();
     SOURCM.lsouit = true;
     flukam(1);
     if (fVerbosityLevel >=3)
-       cout << "<== TFluka::ProcessEvent() called." << endl;
+        cout << "<== TFluka::ProcessEvent() called." << endl;
     //
     // Increase event number
     //
@@ -317,7 +352,7 @@ Bool_t TFluka::ProcessRun(Int_t nevent) {
 
   if (fVerbosityLevel >=3)
     cout << "==> TFluka::ProcessRun(" << nevent << ") called." 
-        << endl;
+         << endl;
 
   if (fVerbosityLevel >=2) {
     cout << "\t* GLOBAL.fdrtr = " << (GLOBAL.lfdrtr?'T':'F') << endl;
@@ -326,16 +361,23 @@ Bool_t TFluka::ProcessRun(Int_t nevent) {
 
   Int_t todo = TMath::Abs(nevent);
   for (Int_t ev = 0; ev < todo; ev++) {
+      TStopwatch timer;
+      timer.Start();
       fApplication->BeginEvent();
       ProcessEvent();
       fApplication->FinishEvent();
+      cout << "Event: "<< ev
+           << Form(" R:%.2fs C:%.2fs", timer.RealTime(),timer.CpuTime()) << endl;
   }
 
   if (fVerbosityLevel >=3)
     cout << "<== TFluka::ProcessRun(" << nevent << ") called." 
-        << endl;
+         << endl;
+  
   // Write fluka specific scoring output
+  genout();
   newplo();
+  flkend();
   
   return kTRUE;
 }
@@ -346,8 +388,8 @@ Bool_t TFluka::ProcessRun(Int_t nevent) {
 // functions from GCONS 
 //____________________________________________________________________________ 
 void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,  
-                   Float_t &dens, Float_t &radl, Float_t &absl,
-                   Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
+                    Float_t &dens, Float_t &radl, Float_t &absl,
+                    Float_t* /*ubuf*/, Int_t& /*nbuf*/) {
 //
    TGeoMaterial *mat;
    TIter next (gGeoManager->GetListOfMaterials());
@@ -368,8 +410,8 @@ void TFluka::Gfmate(Int_t imat, char *name, Float_t &a, Float_t &z,
 
 //______________________________________________________________________________ 
 void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,  
-                   Double_t &dens, Double_t &radl, Double_t &absl,
-                   Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
+                    Double_t &dens, Double_t &radl, Double_t &absl,
+                    Double_t* /*ubuf*/, Int_t& /*nbuf*/) {
 //
    TGeoMaterial *mat;
    TIter next (gGeoManager->GetListOfMaterials());
@@ -391,8 +433,8 @@ void TFluka::Gfmate(Int_t imat, char *name, Double_t &a, Double_t &z,
 // detector composition
 //______________________________________________________________________________ 
 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
-                     Double_t z, Double_t dens, Double_t radl, Double_t absl,
-                     Float_t* buf, Int_t nwbuf) {
+                      Double_t z, Double_t dens, Double_t radl, Double_t absl,
+                      Float_t* buf, Int_t nwbuf) {
 //
    Double_t* dbuf = fGeom->CreateDoubleArray(buf, nwbuf);  
    Material(kmat, name, a, z, dens, radl, absl, dbuf, nwbuf);
@@ -401,8 +443,8 @@ void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 
 //______________________________________________________________________________ 
 void TFluka::Material(Int_t& kmat, const char* name, Double_t a, 
-                     Double_t z, Double_t dens, Double_t radl, Double_t absl,
-                     Double_t* /*buf*/, Int_t /*nwbuf*/) {
+                      Double_t z, Double_t dens, Double_t radl, Double_t absl,
+                      Double_t* /*buf*/, Int_t /*nwbuf*/) {
 //
 // Define a material
   TGeoMaterial *mat;
@@ -420,7 +462,7 @@ void TFluka::Material(Int_t& kmat, const char* name, Double_t a,
 
 //______________________________________________________________________________ 
 void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a, 
-                    Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
+                     Float_t *z, Double_t dens, Int_t nlmat, Float_t *wmat) {
 //
 // Define a material mixture
 //
@@ -440,7 +482,7 @@ void TFluka::Mixture(Int_t& kmat, const char *name, Float_t *a,
 
 //______________________________________________________________________________ 
 void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a, 
-                    Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
+                     Double_t *z, Double_t dens, Int_t nlmat, Double_t *wmat) {
 //
   // Defines mixture OR COMPOUND IMAT as composed by 
   // THE BASIC NLMAT materials defined by arrays A,Z and WMAT
@@ -554,33 +596,33 @@ void TFluka::Mixture(Int_t& kmat, const char *name, Double_t *a,
 
 //______________________________________________________________________________ 
 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
-                   Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
-                   Double_t stemax, Double_t deemax, Double_t epsil, 
-                   Double_t stmin, Float_t* ubuf, Int_t nbuf) { 
+                    Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
+                    Double_t stemax, Double_t deemax, Double_t epsil,
+                    Double_t stmin, Float_t* ubuf, Int_t nbuf) {
   // Define a medium
   // 
   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
-            epsil, stmin, ubuf, nbuf);
+             epsil, stmin, ubuf, nbuf);
 } 
 
 //______________________________________________________________________________ 
 void TFluka::Medium(Int_t& kmed, const char *name, Int_t nmat, 
-                   Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd, 
-                   Double_t stemax, Double_t deemax, Double_t epsil, 
-                   Double_t stmin, Double_t* ubuf, Int_t nbuf) { 
+                    Int_t isvol, Int_t ifield, Double_t fieldm, Double_t tmaxfd,
+                    Double_t stemax, Double_t deemax, Double_t epsil,
+                    Double_t stmin, Double_t* ubuf, Int_t nbuf) {
   // Define a medium
   // 
   kmed = gGeoManager->GetListOfMedia()->GetSize()+1;
   fMCGeo->Medium(kmed, name, nmat, isvol, ifield, fieldm, tmaxfd, stemax, deemax, 
-            epsil, stmin, ubuf, nbuf);
+             epsil, stmin, ubuf, nbuf);
 } 
 
 //______________________________________________________________________________ 
 void TFluka::Matrix(Int_t& krot, Double_t thetaX, Double_t phiX, 
-                   Double_t thetaY, Double_t phiY, Double_t thetaZ, 
-                   Double_t phiZ) {
-//                  
+                    Double_t thetaY, Double_t phiY, Double_t thetaZ,
+                    Double_t phiZ) {
+//        
   krot = gGeoManager->GetListOfMatrices()->GetEntriesFast();
   fMCGeo->Matrix(krot, thetaX, phiX, thetaY, phiY, thetaZ, phiZ); 
 } 
@@ -645,42 +687,42 @@ void TFluka::Gsatt(const char *name, const char *att, Int_t val)
 
 //______________________________________________________________________________ 
 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
-                    Float_t *upar, Int_t np)  {
+                     Float_t *upar, Int_t np)  {
 //
     return fMCGeo->Gsvolu(name, shape, nmed, upar, np); 
 }
 
 //______________________________________________________________________________ 
 Int_t TFluka::Gsvolu(const char *name, const char *shape, Int_t nmed,  
-                    Double_t *upar, Int_t np)  {
+                     Double_t *upar, Int_t np)  {
 //
     return fMCGeo->Gsvolu(name, shape, nmed, upar, np); 
 }
  
 //______________________________________________________________________________ 
 void TFluka::Gsdvn(const char *name, const char *mother, Int_t ndiv, 
-                  Int_t iaxis) {
+                   Int_t iaxis) {
 //
     fMCGeo->Gsdvn(name, mother, ndiv, iaxis); 
 } 
 
 //______________________________________________________________________________ 
 void TFluka::Gsdvn2(const char *name, const char *mother, Int_t ndiv, 
-                   Int_t iaxis, Double_t c0i, Int_t numed) {
+                    Int_t iaxis, Double_t c0i, Int_t numed) {
 //
     fMCGeo->Gsdvn2(name, mother, ndiv, iaxis, c0i, numed); 
 } 
 
 //______________________________________________________________________________ 
 void TFluka::Gsdvt(const char *name, const char *mother, Double_t step, 
-                  Int_t iaxis, Int_t numed, Int_t ndvmx) {
-//     
+                   Int_t iaxis, Int_t numed, Int_t ndvmx) {
+//        
     fMCGeo->Gsdvt(name, mother, step, iaxis, numed, ndvmx); 
 } 
 
 //______________________________________________________________________________ 
 void TFluka::Gsdvt2(const char *name, const char *mother, Double_t step, 
-                   Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) { 
+                    Int_t iaxis, Double_t c0, Int_t numed, Int_t ndvmx) {
 //
     fMCGeo->Gsdvt2(name, mother, step, iaxis, c0, numed, ndvmx); 
 } 
@@ -693,24 +735,24 @@ void TFluka::Gsord(const char * /*name*/, Int_t /*iax*/) {
 
 //______________________________________________________________________________ 
 void TFluka::Gspos(const char *name, Int_t nr, const char *mother,  
-                  Double_t x, Double_t y, Double_t z, Int_t irot, 
-                  const char *konly) {
+                   Double_t x, Double_t y, Double_t z, Int_t irot,
+                   const char *konly) {
 //
   fMCGeo->Gspos(name, nr, mother, x, y, z, irot, konly); 
 } 
 
 //______________________________________________________________________________ 
 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
-                   Double_t x, Double_t y, Double_t z, Int_t irot,
-                   const char *konly, Float_t *upar, Int_t np)  {
+                    Double_t x, Double_t y, Double_t z, Int_t irot,
+                    const char *konly, Float_t *upar, Int_t np)  {
   //
   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
 } 
 
 //______________________________________________________________________________ 
 void TFluka::Gsposp(const char *name, Int_t nr, const char *mother,  
-                   Double_t x, Double_t y, Double_t z, Int_t irot,
-                   const char *konly, Double_t *upar, Int_t np)  {
+                    Double_t x, Double_t y, Double_t z, Int_t irot,
+                    const char *konly, Double_t *upar, Int_t np)  {
   //
   fMCGeo->Gsposp(name, nr, mother, x, y, z, irot, konly, upar, np); 
 } 
@@ -825,7 +867,7 @@ Bool_t TFluka::GetMedium(const TString &volumeName,TString &name,
 
 //______________________________________________________________________________ 
 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
-                        Float_t* absco, Float_t* effic, Float_t* rindex) {
+                         Float_t* absco, Float_t* effic, Float_t* rindex) {
 //
 // Set Cerenkov properties for medium itmed
 //
@@ -847,7 +889,7 @@ void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 }  
 
 void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
-                        Float_t* absco, Float_t* effic, Float_t* rindex, Float_t* rfl) {
+                         Float_t* absco, Float_t* effic, Float_t* rindex, Float_t* rfl) {
 //
 // Set Cerenkov properties for medium itmed
 //
@@ -871,13 +913,13 @@ void TFluka::SetCerenkov(Int_t itmed, Int_t npckov, Float_t* ppckov,
 
 //______________________________________________________________________________ 
 void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t * /*ppckov*/,
-                        Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
+                         Double_t * /*absco*/, Double_t * /*effic*/, Double_t * /*rindex*/) {
 //
 //  Double_t version not implemented
 }  
 
 void TFluka::SetCerenkov(Int_t /*itmed*/, Int_t /*npckov*/, Double_t* /*ppckov*/,
-                        Double_t* /*absco*/, Double_t* /*effic*/, Double_t* /*rindex*/, Double_t* /*rfl*/) { 
+                         Double_t* /*absco*/, Double_t* /*effic*/, Double_t* /*rindex*/, Double_t* /*rfl*/) {
 //
 // //  Double_t version not implemented
 }
@@ -946,34 +988,34 @@ Int_t TFluka::PDGFromId(Int_t id) const
 
     if (id == kFLUKAoptical) {
 // Cerenkov photon
-       if (fVerbosityLevel >= 3)
-           printf("\n PDGFromId: Cerenkov Photon \n");
-       return  50000050;
+//        if (fVerbosityLevel >= 3)
+//            printf("\n PDGFromId: Cerenkov Photon \n");
+        return  50000050;
     }
 // Error id    
     if (id == 0 || id < kFLUKAcodemin || id > kFLUKAcodemax) {
-       if (fVerbosityLevel >= 3)
-           printf("PDGFromId: Error id = 0\n");
-       return -1;
+        if (fVerbosityLevel >= 3)
+            printf("PDGFromId: Error id = 0\n");
+        return -1;
     }
 // Good id    
     if (id > 0) {
-       Int_t intfluka = GetFlukaIPTOKP(id);
-       if (intfluka == 0) {
-           if (fVerbosityLevel >= 3)
-               printf("PDGFromId: Error intfluka = 0: %d\n", id);
-           return -1;
-       } else if (intfluka < 0) {
-           if (fVerbosityLevel >= 3)
-               printf("PDGFromId: Error intfluka < 0: %d\n", id);
-           return -1;
-       }
-//     if (fVerbosityLevel >= 3)
-//         printf("mpdgha called with %d %d \n", id, intfluka);
-       return mpdgha(intfluka);
+        Int_t intfluka = GetFlukaIPTOKP(id);
+        if (intfluka == 0) {
+            if (fVerbosityLevel >= 3)
+                printf("PDGFromId: Error intfluka = 0: %d\n", id);
+            return -1;
+        } else if (intfluka < 0) {
+            if (fVerbosityLevel >= 3)
+                printf("PDGFromId: Error intfluka < 0: %d\n", id);
+            return -1;
+        }
+//        if (fVerbosityLevel >= 3)
+//            printf("mpdgha called with %d %d \n", id, intfluka);
+        return mpdgha(intfluka);
     } else {
-       // ions and optical photons
-       return idSpecial[id - kFLUKAcodemin];
+        // ions and optical photons
+        return idSpecial[id - kFLUKAcodemin];
     }
 }
 
@@ -1002,10 +1044,10 @@ void TFluka::SetProcess(const char* flagName, Int_t flagValue, Int_t imed)
     TFlukaConfigOption* proc;
     while((proc = (TFlukaConfigOption*)next()))
     { 
-       if (proc->Medium() == imed) {
-           proc->SetProcess(flagName, flagValue);
-           return;
-       }
+        if (proc->Medium() == imed) {
+            proc->SetProcess(flagName, flagValue);
+            return;
+        }
     }
     proc = new TFlukaConfigOption(imed);
     proc->SetProcess(flagName, flagValue);
@@ -1031,10 +1073,10 @@ void TFluka::SetCut(const char* cutName, Double_t cutValue, Int_t imed)
     TFlukaConfigOption* proc;
     while((proc = (TFlukaConfigOption*)next()))
     { 
-       if (proc->Medium() == imed) {
-           proc->SetCut(cutName, cutValue);
-           return;
-       }
+        if (proc->Medium() == imed) {
+            proc->SetCut(cutName, cutValue);
+            return;
+        }
     }
 
     proc = new TFlukaConfigOption(imed);
@@ -1052,10 +1094,10 @@ void TFluka::SetModelParameter(const char* parName, Double_t parValue, Int_t ime
     TFlukaConfigOption* proc;
     while((proc = (TFlukaConfigOption*)next()))
     { 
-       if (proc->Medium() == imed) {
-           proc->SetModelParameter(parName, parValue);
-           return;
-       }
+        if (proc->Medium() == imed) {
+            proc->SetModelParameter(parName, parValue);
+            return;
+        }
     }
 
     proc = new TFlukaConfigOption(imed);
@@ -1115,16 +1157,16 @@ void TFluka::InitPhysics()
     
 // Open files 
     if ((pFlukaVmcCoreInp = fopen(sFlukaVmcCoreInp.Data(),"r")) == NULL) {
-       Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
-       exit(1);
+        Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcCoreInp.Data());
+        exit(1);
     }
     if ((pFlukaVmcFlukaMat = fopen(sFlukaVmcTmp.Data(),"r")) == NULL) {
-       Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcTmp.Data());
-       exit(1);
+        Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcTmp.Data());
+        exit(1);
     }
     if ((pFlukaVmcInp = fopen(sFlukaVmcInp.Data(),"w")) == NULL) {
-       Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcInp.Data());
-       exit(1);
+        Warning("InitPhysics", "\nCannot open file %s\n",sFlukaVmcInp.Data());
+        exit(1);
     }
 
 // Copy core input file 
@@ -1132,27 +1174,27 @@ void TFluka::InitPhysics()
     Float_t fEventsPerRun;
     
     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) {
-       if (strncmp(sLine,"GEOEND",6) != 0)
-           fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
-       else {
-           fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
-           goto flukamat;
-       }
+        if (strncmp(sLine,"GEOEND",6) != 0)
+            fprintf(pFlukaVmcInp,"%s",sLine); // copy until GEOEND card
+        else {
+            fprintf(pFlukaVmcInp,"GEOEND\n");   // add GEOEND card
+            goto flukamat;
+        }
     } // end of while until GEOEND card
     
 
  flukamat:
     while ((fgets(sLine,255,pFlukaVmcFlukaMat)) != NULL) { // copy flukaMat.inp file
-       fprintf(pFlukaVmcInp,"%s\n",sLine);
+        fprintf(pFlukaVmcInp,"%s\n",sLine);
     }
     
     while ((fgets(sLine,255,pFlukaVmcCoreInp)) != NULL) { 
-       if (strncmp(sLine,"START",5) != 0)
-           fprintf(pFlukaVmcInp,"%s\n",sLine);
-       else {
-           sscanf(sLine+10,"%10f",&fEventsPerRun);
-           goto fin;
-       }
+        if (strncmp(sLine,"START",5) != 0)
+            fprintf(pFlukaVmcInp,"%s\n",sLine);
+        else {
+            sscanf(sLine+10,"%10f",&fEventsPerRun);
+            goto fin;
+        }
     } //end of while until START card
     
  fin:
@@ -1179,32 +1221,32 @@ void TFluka::InitPhysics()
 
     for (Int_t isc = 0; isc < nscore; isc++) 
     {
-       mopo = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
-       char*    fileName = mopo->GetFileName();
-       Int_t    size     = strlen(fileName);
-       Float_t  lun      = -1.;
+        mopo = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isc));
+        char*    fileName = mopo->GetFileName();
+        Int_t    size     = strlen(fileName);
+        Float_t  lun      = -1.;
 //
 // Check if new output file has to be opened
-       for (Int_t isci = 0; isci < isc; isci++) {
-
-           
-           mopi = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isci));
-           if(strncmp(mopi->GetFileName(), fileName, size)==0) {
-               // 
-               // No, the file already exists
-               lun = mopi->GetLun();
-               mopo->SetLun(lun);
-               break;
-           }
-       } // inner loop
-
-       if (lun == -1.) {
-           // Open new output file
-           inp++;
-           mopo->SetLun(loginp + inp);
-           mopo->WriteOpenFlukaFile();
-       }
-       mopo->WriteFlukaInputCards();
+        for (Int_t isci = 0; isci < isc; isci++) {
+
+        
+            mopi = dynamic_cast<TFlukaScoringOption*> (fUserScore->At(isci));
+            if(strncmp(mopi->GetFileName(), fileName, size)==0) {
+                //
+                // No, the file already exists
+                lun = mopi->GetLun();
+                mopo->SetLun(lun);
+                break;
+            }
+        } // inner loop
+
+        if (lun == -1.) {
+            // Open new output file
+            inp++;
+            mopo->SetLun(loginp + inp);
+            mopo->WriteOpenFlukaFile();
+        }
+        mopo->WriteFlukaInputCards();
     }
 
 // Add RANDOMIZ card
@@ -1228,9 +1270,9 @@ void TFluka::InitPhysics()
     
     for (Int_t im = 0; im < nmaterial; im++)
     {
-       TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
-       Int_t idmat = material->GetIndex();
-       fMaterials[idmat] = im;
+        TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
+        Int_t idmat = material->GetIndex();
+        fMaterials[idmat] = im;
     }
 } // end of InitPhysics
 
@@ -1239,10 +1281,11 @@ void TFluka::InitPhysics()
 void TFluka::SetMaxStep(Double_t step)
 {
 // Set the maximum step size
-    if (step > 1.e4) return;
+//    if (step > 1.e4) return;
     
-    Int_t mreg, latt;
-    fGeom->GetCurrentRegion(mreg, latt);
+//    Int_t mreg=0, latt=0;
+//    fGeom->GetCurrentRegion(mreg, latt);
+    Int_t mreg = fGeom->GetCurrentRegion();
     STEPSZ.stepmx[mreg - 1] = step;
 }
 
@@ -1526,7 +1569,7 @@ Int_t TFluka::CorrectFlukaId() const
    // since we don't put photons and e- created bellow transport cut on the vmc stack
    // and there is a call to endraw for energy deposition for each of them
    // and they have the track number of their parent, but different identity (pdg)
-   // so we want to assign also their parent identity also.
+   // so we want to assign also their parent identity.
    if( (IsTrackStop() )
         && TRACKR.ispusr[mkbmx2 - 4] == TRACKR.ispusr[mkbmx2 - 1]
         && TRACKR.jtrack != TRACKR.ispusr[mkbmx2 - 3] ) {
@@ -1725,69 +1768,69 @@ Int_t TFluka::NSecondaries() const
 // FHEAVY.npheav = number of secondaries for light and heavy secondary ions
     FlukaCallerCode_t caller = GetCaller();
     if (caller == kUSDRAW)  // valid only after usdraw
-       return GENSTK.np + FHEAVY.npheav;
+        return GENSTK.np + FHEAVY.npheav;
     else if (caller == kUSTCKV) {
-       // Cerenkov Photon production
-       return fNCerenkov;
+        // Cerenkov Photon production
+        return fNCerenkov;
     }
     return 0;
 } // end of NSecondaries
 
 //______________________________________________________________________________ 
 void TFluka::GetSecondary(Int_t isec, Int_t& particleId,
-               TLorentzVector& position, TLorentzVector& momentum)
+                TLorentzVector& position, TLorentzVector& momentum)
 {
 // Copy particles from secondary stack to vmc stack
 //
 
     FlukaCallerCode_t caller = GetCaller();
     if (caller == kUSDRAW) {  // valid only after usdraw
-       if (GENSTK.np > 0) {
-           // Hadronic interaction
-           if (isec >= 0 && isec < GENSTK.np) {
-               particleId = PDGFromId(GENSTK.kpart[isec]);
-               position.SetX(fXsco);
-               position.SetY(fYsco);
-               position.SetZ(fZsco);
-               position.SetT(TRACKR.atrack);
-               momentum.SetPx(GENSTK.plr[isec]*GENSTK.cxr[isec]);
-               momentum.SetPy(GENSTK.plr[isec]*GENSTK.cyr[isec]);
-               momentum.SetPz(GENSTK.plr[isec]*GENSTK.czr[isec]);
-               momentum.SetE(GENSTK.tki[isec] + PAPROP.am[GENSTK.kpart[isec]+6]);
-           }
-           else if (isec >= GENSTK.np && isec < GENSTK.np + FHEAVY.npheav) {
-               Int_t jsec = isec - GENSTK.np;
-               particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
-               position.SetX(fXsco);
-               position.SetY(fYsco);
-               position.SetZ(fZsco);
-               position.SetT(TRACKR.atrack);
-               momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
-               momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
-               momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
-               if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6) 
-                   momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
-               else if (FHEAVY.tkheav[jsec] > 6)
-                   momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
-           }
-           else
-               Warning("GetSecondary","isec out of range");
-       } 
+        if (GENSTK.np > 0) {
+            // Hadronic interaction
+            if (isec >= 0 && isec < GENSTK.np) {
+                particleId = PDGFromId(GENSTK.kpart[isec]);
+                position.SetX(fXsco);
+                position.SetY(fYsco);
+                position.SetZ(fZsco);
+                position.SetT(TRACKR.atrack);
+                momentum.SetPx(GENSTK.plr[isec]*GENSTK.cxr[isec]);
+                momentum.SetPy(GENSTK.plr[isec]*GENSTK.cyr[isec]);
+                momentum.SetPz(GENSTK.plr[isec]*GENSTK.czr[isec]);
+                momentum.SetE(GENSTK.tki[isec] + PAPROP.am[GENSTK.kpart[isec]+6]);
+            }
+            else if (isec >= GENSTK.np && isec < GENSTK.np + FHEAVY.npheav) {
+                Int_t jsec = isec - GENSTK.np;
+                particleId = FHEAVY.kheavy[jsec]; // this is Fluka id !!!
+                position.SetX(fXsco);
+                position.SetY(fYsco);
+                position.SetZ(fZsco);
+                position.SetT(TRACKR.atrack);
+                momentum.SetPx(FHEAVY.pheavy[jsec]*FHEAVY.cxheav[jsec]);
+                momentum.SetPy(FHEAVY.pheavy[jsec]*FHEAVY.cyheav[jsec]);
+                momentum.SetPz(FHEAVY.pheavy[jsec]*FHEAVY.czheav[jsec]);
+                if (FHEAVY.tkheav[jsec] >= 3 && FHEAVY.tkheav[jsec] <= 6)
+                    momentum.SetE(FHEAVY.tkheav[jsec] + PAPROP.am[jsec+6]);
+                else if (FHEAVY.tkheav[jsec] > 6)
+                    momentum.SetE(FHEAVY.tkheav[jsec] + FHEAVY.amnhea[jsec]); // to be checked !!!
+            }
+            else
+                Warning("GetSecondary","isec out of range");
+        }
     } else if (caller == kUSTCKV) {
-       Int_t index = OPPHST.lstopp - isec;
-       position.SetX(OPPHST.xoptph[index]);
-       position.SetY(OPPHST.yoptph[index]);
-       position.SetZ(OPPHST.zoptph[index]);
-       position.SetT(OPPHST.agopph[index]);
-       Double_t p = OPPHST.poptph[index];
-       
-       momentum.SetPx(p * OPPHST.txopph[index]);
-       momentum.SetPy(p * OPPHST.tyopph[index]);
-       momentum.SetPz(p * OPPHST.tzopph[index]);
-       momentum.SetE(p);
+        Int_t index = OPPHST.lstopp - isec;
+        position.SetX(OPPHST.xoptph[index]);
+        position.SetY(OPPHST.yoptph[index]);
+        position.SetZ(OPPHST.zoptph[index]);
+        position.SetT(OPPHST.agopph[index]);
+        Double_t p = OPPHST.poptph[index];
+        
+        momentum.SetPx(p * OPPHST.txopph[index]);
+        momentum.SetPy(p * OPPHST.tyopph[index]);
+        momentum.SetPz(p * OPPHST.tzopph[index]);
+        momentum.SetE(p);
     }
     else
-       Warning("GetSecondary","no secondaries available");
+        Warning("GetSecondary","no secondaries available");
     
 } // end of GetSecondary
 
@@ -1800,8 +1843,8 @@ TMCProcess TFluka::ProdProcess(Int_t) const
 // in the current step
 
     Int_t mugamma = (TRACKR.jtrack == kFLUKAphoton || 
-                    TRACKR.jtrack == kFLUKAmuplus || 
-                    TRACKR.jtrack == kFLUKAmuminus);
+                     TRACKR.jtrack == kFLUKAmuplus ||
+                     TRACKR.jtrack == kFLUKAmuminus);
     FlukaProcessCode_t icode = GetIcode();
 
     if      (icode == kKASKADdecay)                                   return kPDecay;
@@ -1813,9 +1856,9 @@ TMCProcess TFluka::ProdProcess(Int_t) const
     else if (icode == kEMFSCOmoller     || icode == kEMFSCObhabha)    return kPDeltaRay;
     else if (icode == kEMFSCOanniflight || icode == kEMFSCOannirest)  return kPAnnihilation;
     else if (icode == kKASKADinelint) {
-       if (!mugamma)                                                 return kPHadronic;
-       else if (TRACKR.jtrack == kFLUKAphoton)                       return kPPhotoFission;
-       else                                                          return kPMuonNuclear;
+        if (!mugamma)                                                 return kPHadronic;
+        else if (TRACKR.jtrack == kFLUKAphoton)                       return kPPhotoFission;
+        else                                                          return kPMuonNuclear;
     }
     else if (icode == kEMFSCOrayleigh)                                return kPRayleigh;
 // Fluka codes 100, 300 and 400 still to be investigasted
@@ -1837,8 +1880,8 @@ Int_t TFluka::StepProcesses(TArrayI &proc) const
     case kKASNEUtimekill:
     case kKASHEAtimekill:
     case kKASOPHtimekill:
-       iproc =  kPTOFlimit;
-       break;
+        iproc =  kPTOFlimit;
+        break;
     case kKASKADstopping:
     case kKASKADescape:
     case kEMFSCOstopping1:
@@ -1848,18 +1891,18 @@ Int_t TFluka::StepProcesses(TArrayI &proc) const
     case kKASNEUescape:
     case kKASHEAescape:
     case kKASOPHescape:
-       iproc = kPStop;
-       break;
+        iproc = kPStop;
+        break;
     case kKASOPHabsorption:
-       iproc = kPLightAbsorption;
-       break;
+        iproc = kPLightAbsorption;
+        break;
     case kKASOPHrefraction:
-       iproc = kPLightRefraction;
+        iproc = kPLightRefraction;
     case kEMSCOlocaledep : 
-       iproc = kPPhotoelectric;
-       break;
+        iproc = kPPhotoelectric;
+        break;
     default:
-       iproc = ProdProcess(0);
+        iproc = ProdProcess(0);
     }
     proc[0] = iproc;
     return 1;
@@ -1955,7 +1998,7 @@ const char* TFluka::CurrentVolPath() {
 }
 //______________________________________________________________________________ 
 Int_t TFluka::CurrentMaterial(Float_t & a, Float_t & z, 
-                     Float_t & dens, Float_t & radl, Float_t & absl) const
+                      Float_t & dens, Float_t & radl, Float_t & absl) const
 {
 //
 //  Return the current medium number and material properties
@@ -2068,6 +2111,7 @@ void TFluka::SetMreg(Int_t l, Int_t lttc)
 
 
 
+//______________________________________________________________________________
 TString TFluka::ParticleName(Int_t pdg) const
 {
     // Return particle name for particle with pdg code pdg.
@@ -2076,6 +2120,7 @@ TString TFluka::ParticleName(Int_t pdg) const
 }
  
 
+//______________________________________________________________________________
 Double_t TFluka::ParticleMass(Int_t pdg) const
 {
     // Return particle mass for particle with pdg code pdg.
@@ -2083,12 +2128,14 @@ Double_t TFluka::ParticleMass(Int_t pdg) const
     return (PAPROP.am[ifluka - kFLUKAcodemin]);
 }
 
+//______________________________________________________________________________
 Double_t TFluka::ParticleMassFPC(Int_t fpc) const
 {
     // Return particle mass for particle with Fluka particle code fpc
     return (PAPROP.am[fpc - kFLUKAcodemin]);
 }
 
+//______________________________________________________________________________
 Double_t TFluka::ParticleCharge(Int_t pdg) const
 {
     // Return particle charge for particle with pdg code pdg.
@@ -2096,6 +2143,7 @@ Double_t TFluka::ParticleCharge(Int_t pdg) const
     return Double_t(PAPROP.ichrge[ifluka - kFLUKAcodemin]);
 }
 
+//______________________________________________________________________________
 Double_t TFluka::ParticleLifeTime(Int_t pdg) const
 {
     // Return particle lifetime for particle with pdg code pdg.
@@ -2103,6 +2151,7 @@ Double_t TFluka::ParticleLifeTime(Int_t pdg) const
     return (PAPROP.tmnlf[ifluka - kFLUKAcodemin]);
 }
 
+//______________________________________________________________________________
 void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t& charge, Float_t& tlife)
 {
     // Retrieve particle properties for particle with pdg code pdg.
@@ -2114,6 +2163,7 @@ void TFluka::Gfpart(Int_t pdg, char* name, Int_t& type, Float_t& mass, Float_t&
     tlife  = ParticleLifeTime(pdg);
 }
 
+//______________________________________________________________________________
 void TFluka::PrintHeader()
 {
     //
@@ -2135,8 +2185,8 @@ void TFluka::PrintHeader()
 
 extern "C" {
   void pshckp(Double_t & px, Double_t & py, Double_t & pz, Double_t & e,
-             Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
-             Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
+              Double_t & vx, Double_t & vy, Double_t & vz, Double_t & tof,
+              Double_t & polx, Double_t & poly, Double_t & polz, Double_t & wgt, Int_t& ntr)
   {
     //
     // Pushes one cerenkov photon to the stack
@@ -2146,30 +2196,56 @@ extern "C" {
     TVirtualMCStack* cppstack = fluka->GetStack();
     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
     cppstack->PushTrack(0, parent, 50000050,
-                       px, py, pz, e,
-                       vx, vy, vz, tof,
-                       polx, poly, polz,
-                       kPCerenkov, ntr, wgt, 0);
+                        px, py, pz, e,
+                        vx, vy, vz, tof,
+                        polx, poly, polz,
+                        kPCerenkov, ntr, wgt, 0);
+    if (fluka->GetVerbosityLevel() >= 3)
+            printf("pshckp: track=%d parent=%d lattc=%d %s\n", ntr, parent, TRACKR.lt1trk, fluka->CurrentVolName());
   }
     
     void ustckv(Int_t & nphot, Int_t & mreg, Double_t & x, Double_t & y, Double_t & z)
     {
-       //
-       // Calls stepping in order to signal cerenkov production
-       //
-       TFluka *fluka = (TFluka*)gMC;
-       fluka->SetMreg(mreg,LTCLCM.mlatm1);
-       fluka->SetXsco(x);
-       fluka->SetYsco(y);
-       fluka->SetZsco(z);
-       fluka->SetNCerenkov(nphot);
-       fluka->SetCaller(kUSTCKV);
-       if (fluka->GetVerbosityLevel() >= 3) 
-       (TVirtualMCApplication::Instance())->Stepping();
-       
+        //
+        // Calls stepping in order to signal cerenkov production
+        //
+        TFluka *fluka = (TFluka*)gMC;
+        fluka->SetMreg(mreg, TRACKR.lt1trk); //LTCLCM.mlatm1);
+        fluka->SetXsco(x);
+        fluka->SetYsco(y);
+        fluka->SetZsco(z);
+        fluka->SetNCerenkov(nphot);
+        fluka->SetCaller(kUSTCKV);
+        if (fluka->GetVerbosityLevel() >= 3)
+            printf("ustckv: %10d mreg=%d lattc=%d  newlat=%d (%f, %f, %f) edep=%f vol=%s\n",
+                    nphot, mreg, TRACKR.lt1trk, LTCLCM.newlat, x, y, z, fluka->Edep(), fluka->CurrentVolName());
+   
+    // check region lattice consistency (debug Ernesto)
+    // *****************************************************
+   Int_t nodeId;
+   Int_t volId = fluka->CurrentVolID(nodeId);
+   Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
+
+   if( mreg != volId  && !gGeoManager->IsOutside() ) {
+       cout << "  ustckv:   track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+            << " icode=" << fluka->GetIcode() << " gNstep=" << fluka->GetNstep() << endl
+            << "               fluka   mreg=" << mreg << " mlttc=" << TRACKR.lt1trk << endl
+            << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
+            << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
+            << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
+            << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
+            << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
+        if( TRACKR.lt1trk == crtlttc ) cout << "   *************************************************************" << endl;
+    }
+    // *****************************************************
+
+
+
+        (TVirtualMCApplication::Instance())->Stepping();
     }
 }
 
+//______________________________________________________________________________
 void TFluka::AddParticlesToPdgDataBase() const
 {
 
@@ -2199,24 +2275,28 @@ void TFluka::AddParticlesToPdgDataBase() const
                      0,6,"Ion",kion+20030);
 }
 
-  //
-  // Info about primary ionization electrons
-  Int_t TFluka::GetNPrimaryElectrons()
+//
+// Info about primary ionization electrons
+//
+
+//______________________________________________________________________________
+Int_t TFluka::GetNPrimaryElectrons()
 {
     // Get number of primary electrons
     return ALLDLT.nalldl;
 }
 
+//______________________________________________________________________________
 Double_t GetPrimaryElectronKineticEnergy(Int_t i)
 {
     Double_t ekin = -1.;
     // Returns kinetic energy of primary electron i
     if (i >= 0 && i < ALLDLT.nalldl) {
-       ekin =  ALLDLT.talldl[i];
+        ekin =  ALLDLT.talldl[i];
     } else {
-       Warning("GetPrimaryElectronKineticEnergy", 
-               "Primary electron index out of range %d %d \n", 
-               i, ALLDLT.nalldl);
+        Warning("GetPrimaryElectronKineticEnergy",
+                "Primary electron index out of range %d %d \n",
+                i, ALLDLT.nalldl);
     }
     return ekin;
 }
index 34da8ee..512b4ce 100644 (file)
 
 #include "TVirtualMC.h"
 #include "TFlukaCodes.h"
+#include "TFlukaMCGeometry.h"
+
 //Forward declaration
 class TGeoMCGeometry;
-class TFlukaMCGeometry;
+//class TFlukaMCGeometry;
 class TGeoMaterial;
 
 class TFluka : public TVirtualMC {
@@ -29,6 +31,8 @@ class TFluka : public TVirtualMC {
   virtual ~TFluka();
   virtual Bool_t IsRootGeometrySupported() const { return kTRUE;}
   
+  Int_t         GetNstep() { return fGeom->GetNstep(); } // to be removed
+  
   //
   // Methods for building/management of geometry
   // ------------------------------------------------
@@ -326,65 +330,61 @@ class TFluka : public TVirtualMC {
   Int_t GetVerbosityLevel() const {return fVerbosityLevel;}
   void SetVerbosityLevel(Int_t l) {fVerbosityLevel = l;}
 
-  // - Fluka Draw procedures identifiers
-  // bxdraw = 1  inside
-  // bxdraw = 11 entering
-  // bxdraw = 12 exiting
-  // eedraw = 2
-  // endraw = 3
-  // mgdraw = 4
-  // sodraw = 5
-  // usdraw = 6
+  //
+  // - Fluka Draw procedures identifiers, see TFlukaCodes.h
+  //
   FlukaCallerCode_t GetCaller() const {return fCaller;}
-  void SetCaller(FlukaCallerCode_t l) {fCaller = l;}
-  
   FlukaProcessCode_t GetIcode() const {return fIcode;}
-  void  SetIcode(FlukaProcessCode_t l) {fIcode = l;}
-
   Int_t GetMreg() const {return fCurrentFlukaRegion;}
-  void SetMreg(Int_t l, Int_t lttc);
-
   Int_t GetNewreg() const {return fNewReg;}
-  void SetNewreg(Int_t l, Int_t /*lttc*/) {fNewReg = l;}
-
   Double_t GetRull() const {return fRull;}
-  void SetRull(Double_t r) {fRull = r;}
-
   Double_t GetXsco() const {return fXsco;}
-  void SetXsco(Double_t x) {fXsco = x;}
-
   Double_t GetYsco() const {return fYsco;}
-  void SetYsco(Double_t y) {fYsco = y;}
-
   Double_t GetZsco() const {return fZsco;}
+  Int_t              GetCurrentFlukaRegion() const {return fCurrentFlukaRegion;}
+  // - Fluka Draw Setters
+  void  SetCurrentFlukaRegion(Int_t reg) {fCurrentFlukaRegion=reg;}
+  void  SetCaller(FlukaCallerCode_t l) {fCaller = l;}
+  void  SetIcode(FlukaProcessCode_t l) {fIcode = l;}
+  void  SetMreg(Int_t l, Int_t lttc);
+  void  SetNewreg(Int_t l, Int_t /*lttc*/) {fNewReg = l;}
+  void  SetRull(Double_t r) {fRull = r;}
+  void  SetXsco(Double_t x) {fXsco = x;}
+  void  SetYsco(Double_t y) {fYsco = y;}
   void SetZsco(Double_t z) {fZsco = z;}
 
-  void SetCurrentFlukaRegion(Int_t reg) {fCurrentFlukaRegion=reg;}
-  Int_t GetCurrentFlukaRegion() const {return fCurrentFlukaRegion;}
+  void  SetTrackIsEntering(){fTrackIsEntering = kTRUE; fTrackIsExiting = kFALSE;}
+  void  SetTrackIsExiting() {fTrackIsExiting  = kTRUE; fTrackIsEntering = kFALSE;}
+  void  SetTrackIsInside()  {fTrackIsExiting  = kFALSE; fTrackIsEntering = kFALSE;}
+  void  SetTrackIsNew(Bool_t flag = kTRUE);
 
   void   SetDummyBoundary(Int_t mode) {fDummyBoundary = mode;}
   Int_t  GetDummyBoundary() const {return fDummyBoundary;}
   Bool_t IsDummyBoundary() const {return (fDummyBoundary==0)?kFALSE:kTRUE;}
   
+
   void   SetGeneratePemf(Bool_t flag = kTRUE) {fGeneratePemf = flag;}
   Bool_t IsGeneratePemf() const {return fGeneratePemf;}
   
   void   EnableField(Bool_t flag=kTRUE) {fFieldFlag = flag;}
   Bool_t IsFieldEnabled() const {return fFieldFlag;}
-  void SetTrackIsEntering(){fTrackIsEntering = kTRUE; fTrackIsExiting = kFALSE;}
-  void SetTrackIsExiting() {fTrackIsExiting  = kTRUE; fTrackIsEntering = kFALSE;}
-  void SetTrackIsInside()  {fTrackIsExiting  = kFALSE; fTrackIsEntering = kFALSE;}
-  void SetTrackIsNew(Bool_t flag = kTRUE);
   
   Int_t GetMaterialIndex(Int_t idmat) const {return fMaterials[idmat];}
-  TObjArray *GetFlukaMaterials();
+
+  TObjArray *          GetFlukaMaterials();
   virtual void SetRootGeometry() {;} // Dummy
   virtual Int_t        NofVolDaughters(const char* volName) const;
   virtual const char*  VolDaughterName(const char* volName, Int_t i) const;
   virtual Int_t        VolDaughterCopyNo(const char* volName, Int_t i) const;
   virtual const char*  CurrentVolPath();
   virtual void         ForceDecayTime(Float_t){;}
+
   private:
+   
+  // Copy constructor and operator= declared but not implemented (-Weff++ flag)
+  TFluka(const TFluka &mc); //: TVirtualMC(mc) {;}
+  TFluka & operator=(const TFluka &); // {return (*this);}
   void PrintHeader();
   void AddParticlesToPdgDataBase() const;
   //
@@ -392,23 +392,19 @@ class TFluka : public TVirtualMC {
   Int_t    GetNPrimaryElectrons();
   Double_t GetPrimaryElectronKineticEnergy(Int_t i);
   //
-  TFluka(const TFluka &mc): TVirtualMC(mc) {;}
-  TFluka & operator=(const TFluka &) {return (*this);}
-
 
-  
   Int_t   fVerbosityLevel; //Verbosity level (0 lowest - 3 highest)
   Int_t   fNEvent;         //Current event number
   TString fInputFileName;     //Name of the real input file 
   TString fCoreInputFileName; //Name of the input file 
 
-  FlukaCallerCode_t     fCaller; //Parameter to indicate who is the caller of the Fluka Draw
-  FlukaProcessCode_t    fIcode;  //Fluka Draw procedures formal parameter 
-  Int_t    fNewReg; //Fluka Draw procedures formal parameter
-  Double_t fRull;   //Fluka Draw procedures formal parameter
-  Double_t fXsco;   //Fluka Draw procedures formal parameter
-  Double_t fYsco;   //Fluka Draw procedures formal parameter
-  Double_t fZsco;   //Fluka Draw procedures formal parameter
+  FlukaCallerCode_t     fCaller;           // Parameter to indicate who is the caller of the Fluka Draw
+  FlukaProcessCode_t    fIcode;            // Fluka Draw procedures formal parameter 
+  Int_t                 fNewReg;           // Fluka Draw procedures formal parameter
+  Double_t              fRull;             // Fluka Draw procedures formal parameter
+  Double_t              fXsco;             // Fluka Draw procedures formal parameter
+  Double_t              fYsco;             // Fluka Draw procedures formal parameter
+  Double_t              fZsco;             // Fluka Draw procedures formal parameter
   Bool_t   fTrackIsEntering;  // Flag for track entering
   Bool_t   fTrackIsExiting;   // Flag for track exiting  
   Bool_t   fTrackIsNew;       // Flag for new track
@@ -433,6 +429,7 @@ class TFluka : public TVirtualMC {
   TObjArray* fUserConfig;            // List of user physics configuration 
   TObjArray* fUserScore;             // List of user scoring options
   
+
   ClassDef(TFluka,1)  //C++ interface to Fluka montecarlo
 
 
index 8c17ce3..53d26b5 100644 (file)
@@ -51,79 +51,100 @@ TFlukaCerenkov::TFlukaCerenkov()
 
 
 TFlukaCerenkov::TFlukaCerenkov(Int_t npckov, Float_t *ppckov, Float_t *absco, Float_t *effic, Float_t *rindex)
+    : fSamples(npckov),
+      fIsMetal(kFALSE),
+      fIsSensitive(kFALSE),
+      fEnergy(new Float_t[fSamples]),
+      fWaveLength(new Float_t[fSamples]),
+      fAbsorptionCoefficient(new Float_t[fSamples]),
+      fQuantumEfficiency(new Float_t[fSamples]),
+      fRefractionIndex(new Float_t[fSamples]),
+      fReflectivity(new Float_t[fSamples]),
+      fMaximumEfficiency(0.)
 {
-// Constructor    
-    fSamples = npckov;
-    fEnergy                = new Float_t[fSamples];
-    fWaveLength            = new Float_t[fSamples];
-    fAbsorptionCoefficient = new Float_t[fSamples];
-    fRefractionIndex       = new Float_t[fSamples];
-    fQuantumEfficiency     = new Float_t[fSamples];
-    fReflectivity          = new Float_t[fSamples];    
-    
+// Constructor
+
+//    fSamples = npckov;
+//    fEnergy                = new Float_t[fSamples];
+//    fWaveLength            = new Float_t[fSamples];
+//    fAbsorptionCoefficient = new Float_t[fSamples];
+//    fRefractionIndex       = new Float_t[fSamples];
+//    fQuantumEfficiency     = new Float_t[fSamples];
+//    fReflectivity          = new Float_t[fSamples];    
+
     for (Int_t i = 0; i < fSamples; i++) {
-       fEnergy[i]             = ppckov[i];
-       fWaveLength[i]         = khc / ppckov[i];
-       if (absco[i] > 0.) {
-           fAbsorptionCoefficient[i]   = 1./ absco[i];
-       } else {
-           fAbsorptionCoefficient[i]   = 1.e15;
-       }
-       fRefractionIndex[i]    = rindex[i];
-       fQuantumEfficiency[i]  = effic[i];
-       //
-       // Find local maximum quantum efficiency
-       if (effic[i] > fMaximumEfficiency) fMaximumEfficiency = effic[i];
-       //
-       // Flag is sensitive if quantum efficiency 0 < eff < 1 for at least one value.
-       if (effic[i] < 1. && effic[i] > 0.) fIsSensitive = 1;
-       // G3 way to define metal
-       if (rindex[0] == 0.) {
-           fIsMetal = kTRUE;
-           fReflectivity[i] = absco[i];
-       }
+        fEnergy[i]             = ppckov[i];
+        fWaveLength[i]         = khc / ppckov[i];
+        if (absco[i] > 0.) {
+            fAbsorptionCoefficient[i]   = 1./ absco[i];
+        } else {
+            fAbsorptionCoefficient[i]   = 1.e15;
+        }
+        fRefractionIndex[i]    = rindex[i];
+        fQuantumEfficiency[i]  = effic[i];
+        //
+        // Find local maximum quantum efficiency
+        if (effic[i] > fMaximumEfficiency) fMaximumEfficiency = effic[i];
+        //
+        // Flag is sensitive if quantum efficiency 0 < eff < 1 for at least one value.
+        if (effic[i] < 1. && effic[i] > 0.) fIsSensitive = 1;
+        // G3 way to define metal
+        if (rindex[0] == 0.) {
+            fIsMetal = kTRUE;
+            fReflectivity[i] = absco[i];
+        }
     }
     // Find global  maximum quantum efficiency
     if (fMaximumEfficiency > GetGlobalMaximumEfficiency()) {
-       SetGlobalMaximumEfficiency(fMaximumEfficiency);
+        SetGlobalMaximumEfficiency(fMaximumEfficiency);
     }
 }
 
 TFlukaCerenkov::TFlukaCerenkov(Int_t npckov, Float_t *ppckov, Float_t *absco, Float_t *effic, Float_t *rindex, Float_t *refl)
+    : fSamples(npckov),
+      fIsMetal(kFALSE),
+      fIsSensitive(kFALSE),
+      fEnergy(new Float_t[fSamples]),
+      fWaveLength(new Float_t[fSamples]),
+      fAbsorptionCoefficient(new Float_t[fSamples]),
+      fQuantumEfficiency(new Float_t[fSamples]),
+      fRefractionIndex(new Float_t[fSamples]),
+      fReflectivity(new Float_t[fSamples]),
+      fMaximumEfficiency(0.)
 {
     // Constructor including reflectivity
-    fSamples = npckov;
-    fEnergy                = new Float_t[fSamples];
-    fWaveLength            = new Float_t[fSamples];
-    fAbsorptionCoefficient = new Float_t[fSamples];
-    fRefractionIndex       = new Float_t[fSamples];
-    fQuantumEfficiency     = new Float_t[fSamples];
+//    fSamples = npckov;
+//    fEnergy                = new Float_t[fSamples];
+//    fWaveLength            = new Float_t[fSamples];
+//    fAbsorptionCoefficient = new Float_t[fSamples];
+//    fRefractionIndex       = new Float_t[fSamples];
+//    fQuantumEfficiency     = new Float_t[fSamples];
     
     
     for (Int_t i = 0; i < fSamples; i++) {
-       fEnergy[i]             = ppckov[i];
-       fWaveLength[i]         = khc / ppckov[i];
-       if (absco[i] > 0.) {
-           fAbsorptionCoefficient[i]   = 1./ absco[i];
-       } else {
-           fAbsorptionCoefficient[i]   = 1.e15;
-       }
-       fRefractionIndex[i]    = rindex[i];
-       fQuantumEfficiency[i]  = effic[i];
-       //
-       // Find local maximum quantum efficiency
-       if (effic[i] > fMaximumEfficiency) fMaximumEfficiency = effic[i];
-       //
-       // Flag is sensitive if quantum efficiency 0 < eff < 1 for at least one value.
-       if (effic[i] < 1. && effic[i] > 0.) fIsSensitive = 1;
-       //
+        fEnergy[i]             = ppckov[i];
+        fWaveLength[i]         = khc / ppckov[i];
+        if (absco[i] > 0.) {
+            fAbsorptionCoefficient[i]   = 1./ absco[i];
+        } else {
+            fAbsorptionCoefficient[i]   = 1.e15;
+        }
+        fRefractionIndex[i]    = rindex[i];
+        fQuantumEfficiency[i]  = effic[i];
+        //
+        // Find local maximum quantum efficiency
+        if (effic[i] > fMaximumEfficiency) fMaximumEfficiency = effic[i];
+        //
+        // Flag is sensitive if quantum efficiency 0 < eff < 1 for at least one value.
+        if (effic[i] < 1. && effic[i] > 0.) fIsSensitive = 1;
+        //
 
     }
     // Find global  maximum quantum efficiency
     if (fMaximumEfficiency > GetGlobalMaximumEfficiency()) {
-       SetGlobalMaximumEfficiency(fMaximumEfficiency);
+        SetGlobalMaximumEfficiency(fMaximumEfficiency);
     }
-    fReflectivity     = new Float_t[fSamples];
+//    fReflectivity     = new Float_t[fSamples];
     for (Int_t i = 0; i < fSamples; i++) fReflectivity[i] = refl[i];
     fIsMetal = kTRUE;
 }
@@ -208,15 +229,15 @@ Float_t TFlukaCerenkov::Interpolate(Float_t value, Float_t* array1, Float_t* arr
 // Interpolate array values 
 //
     if (value < array1[0] && value >= array1[fSamples - 1]) {
-       Warning("Interpolate()", "Photon energy out of range. Returning 0.");
-       return (0.);
+        Warning("Interpolate()", "Photon energy out of range. Returning 0.");
+        return (0.);
     }
     
     Int_t i = TMath::BinarySearch(fSamples, array1, value);
     if (i == (fSamples-1)) {
-       return (array2[i]);
+        return (array2[i]);
     } else {
-       return (array2[i] + (array2[i+1] - array2[i]) / (array1[i+1] - array1[i]) * (value - array1[i]));
+        return (array2[i] + (array2[i+1] - array2[i]) / (array1[i+1] - array1[i]) * (value - array1[i]));
     }
 }
 
index 84c00a5..c1cb4e7 100644 (file)
@@ -48,9 +48,11 @@ public:
     virtual Double_t  GetMaximumEfficiency() const             {return fMaximumEfficiency;}
     static  Double_t  GetGlobalMaximumEfficiency()             {return fgGlobalMaximumEfficiency;}
     static  void      SetGlobalMaximumEfficiency(Double_t eff) {fgGlobalMaximumEfficiency = eff;}
+
  protected:
+
     virtual Float_t  Interpolate(Float_t energy, Float_t* array1, Float_t* array2);
-    
+
  protected:
     Int_t        fSamples;                  // Number of sampling points
     Bool_t       fIsMetal;                  // Flag for metals
@@ -64,7 +66,12 @@ public:
     Double_t     fMaximumEfficiency;        // Local maximum quantum efficiency
     // static 
     static Double_t fgGlobalMaximumEfficiency; // Global maximum quantum efficiency
-    
+
+ private:
+    // Copy constructor and operator= declared but not implemented (-Weff++ flag)
+    TFlukaCerenkov(const TFlukaCerenkov&);
+    TFlukaCerenkov& operator=(const TFlukaCerenkov&);
+
     ClassDef(TFlukaCerenkov, 1)          // CerenkovProperties
 };
        
index 749df16..51ba7bd 100644 (file)
@@ -53,6 +53,7 @@ typedef enum {
 FlukaProcessCode_t;      
 
 typedef enum {
+    kNoCaller       =  0,
     kEEDRAW         =  2,      // Stepping called from eedraw 
     kENDRAW         =  3,      // Stepping called from endraw      
     kMGDRAW         =  4,      // Stepping called from mgdraw 
index 1303794..83b20be 100644 (file)
@@ -42,11 +42,15 @@ ClassImp(TFlukaConfigOption)
 
 
 TFlukaConfigOption::TFlukaConfigOption()
+  :fMedium(-1),
+   fCMatMin(-1),
+   fCMatMax(-1),
+   fCMaterial(0)
 {
     // Default constructor
-    fMedium  = -1;
-    fCMatMin = -1.;
-    fCMatMax = -1.;    
+//    fMedium  = -1;
+//    fCMatMin = -1.;
+//    fCMatMax = -1.;    
     Int_t i;
     for (i = 0; i < 11; i++) fCutValue[i]    = -1.;
     for (i = 0; i < 15; i++) fProcessFlag[i] = -1;
@@ -54,11 +58,15 @@ TFlukaConfigOption::TFlukaConfigOption()
 
 
 TFlukaConfigOption::TFlukaConfigOption(Int_t medium)
+  :fMedium(medium),
+   fCMatMin(-1),
+   fCMatMax(-1),
+   fCMaterial(0)
 {
     // Constructor
-    fMedium = medium;
-    fCMatMin = -1.;
-    fCMatMax = -1.;    
+//    fMedium = medium;
+//    fCMatMin = -1.;
+//    fCMatMax = -1.;    
     Int_t i;
     for (i = 0; i < 11; i++) fCutValue[i]    = -1.;
     for (i = 0; i < 15; i++) fProcessFlag[i] = -1;
@@ -68,14 +76,14 @@ void TFlukaConfigOption::SetCut(const char* flagname, Double_t val)
 {
     // Set a cut value
     const TString cuts[11] = 
-       {"CUTGAM", "CUTELE", "CUTNEU", "CUTHAD", "CUTMUO", "BCUTE", "BCUTM", "DCUTE", "DCUTM", "PPCUTM", "TOFMAX"};
+       {"CUTGAM", "CUTELE", "CUTNEU", "CUTHAD", "CUTMUO", "BCUTE", "BCUTM", "DCUTE", "DCUTM", "PPCUTM", "TOFMAX"};
     Int_t i;
     for (i = 0; i < 11; i++) {
-       if (cuts[i].CompareTo(flagname) == 0) {
-           fCutValue[i] = val;
-           if (fMedium == -1) fgDCutValue[i] = val;
-           break;
-       }
+       if (cuts[i].CompareTo(flagname) == 0) {
+           fCutValue[i] = val;
+           if (fMedium == -1) fgDCutValue[i] = val;
+           break;
+       }
     }
 }
 
@@ -85,10 +93,10 @@ void TFlukaConfigOption::SetModelParameter(const char* flagname, Double_t val)
     const TString parms[2] = {"PRIMIO_N", "PRIMIO_E"};
     Int_t i;
     for (i = 0; i < 2; i++) {
-       if (parms[i].CompareTo(flagname) == 0) {
-           fModelParameter[i] = val;
-           break;
-       }
+       if (parms[i].CompareTo(flagname) == 0) {
+           fModelParameter[i] = val;
+           break;
+       }
     }
 }
 
@@ -97,16 +105,16 @@ void TFlukaConfigOption::SetProcess(const char* flagname, Int_t flag)
 {
     // Set a process flag
     const TString process[15] = 
-       {"DCAY", "PAIR", "COMP", "PHOT", "PFIS", "DRAY", "ANNI", "BREM", "MUNU", "CKOV", 
-        "HADR", "LOSS", "MULS", "RAYL", "STRA"};
+       {"DCAY", "PAIR", "COMP", "PHOT", "PFIS", "DRAY", "ANNI", "BREM", "MUNU", "CKOV",
+        "HADR", "LOSS", "MULS", "RAYL", "STRA"};
     
     Int_t i;
     for (i = 0; i < 15; i++) {
-       if (process[i].CompareTo(flagname) == 0) {
-           fProcessFlag[i] = flag;
-           if (fMedium == -1) fgDProcessFlag[i] = flag;
-           break;
-       }
+       if (process[i].CompareTo(flagname) == 0) {
+           fProcessFlag[i] = flag;
+           if (fMedium == -1) fgDProcessFlag[i] = flag;
+           break;
+       }
     }
 }
 
@@ -124,52 +132,52 @@ void TFlukaConfigOption::WriteFlukaInputCards()
     TGeoMaterial* mat    = 0x0;
 
     if (fMedium != -1) {
-       TFluka* fluka = (TFluka*) gMC;
-       fMedium = fgGeom->GetFlukaMaterial(fMedium);
-       //
-       // Check if material is actually used
-       Int_t* reglist;
-       Int_t nreg;     
-       reglist = fgGeom->GetMaterialList(fMedium, nreg);
-       if (nreg == 0) {
-           // Material not used -- return
-           return;
-       }
-       //
-       // Find material
-       TObjArray *matList = fluka->GetFlukaMaterials();
-       Int_t nmaterial =  matList->GetEntriesFast();
-       fCMaterial = 0;
-       for (Int_t im = 0; im < nmaterial; im++)
-       {
-           fCMaterial = dynamic_cast<TGeoMaterial*> (matList->At(im));
-           Int_t idmat = fCMaterial->GetIndex();
-           if (idmat == fMedium) break;            
-       } // materials
+       TFluka* fluka = (TFluka*) gMC;
+       fMedium = fgGeom->GetFlukaMaterial(fMedium);
+       //
+       // Check if material is actually used
+       Int_t* reglist;
+       Int_t nreg;
+       reglist = fgGeom->GetMaterialList(fMedium, nreg);
+       if (nreg == 0) {
+           // Material not used -- return
+           return;
+       }
+       //
+       // Find material
+       TObjArray *matList = fluka->GetFlukaMaterials();
+       Int_t nmaterial =  matList->GetEntriesFast();
+       fCMaterial = 0;
+       for (Int_t im = 0; im < nmaterial; im++)
+       {
+           fCMaterial = dynamic_cast<TGeoMaterial*> (matList->At(im));
+           Int_t idmat = fCMaterial->GetIndex();
+           if (idmat == fMedium) break;
+       } // materials
         //
-       // Find medium
-       TList *medlist = gGeoManager->GetListOfMedia();
-       TIter next(medlist);
-       while((med = (TGeoMedium*)next()))
-       {
-           mat = med->GetMaterial();
-           if (mat->GetIndex() == fMedium) {
-               medium = med;
-               break;
-           }
-       } // media
-       //
-       // Check if sensitive
-       if (medium->GetParam(0) != 0.) mediumIsSensitive = kTRUE;
-
-
-       fprintf(fgFile,"*\n*Material specific process and cut settings for #%8d %s\n", fMedium, fCMaterial->GetName());
-       fCMatMin = fMedium;
-       fCMatMax = fMedium;
+       // Find medium
+       TList *medlist = gGeoManager->GetListOfMedia();
+       TIter next(medlist);
+       while((med = (TGeoMedium*)next()))
+       {
+           mat = med->GetMaterial();
+           if (mat->GetIndex() == fMedium) {
+              medium = med;
+              break;
+           }
+       } // media
+       //
+       // Check if sensitive
+       if (medium->GetParam(0) != 0.) mediumIsSensitive = kTRUE;
+
+
+       fprintf(fgFile,"*\n*Material specific process and cut settings for #%8d %s\n", fMedium, fCMaterial->GetName());
+       fCMatMin = fMedium;
+       fCMatMax = fMedium;
     } else {
-       fprintf(fgFile,"*\n*Global process and cut settings \n");
-       fCMatMin = fgMatMin;
-       fCMatMax = fgMatMax;
+       fprintf(fgFile,"*\n*Global process and cut settings \n");
+       fCMatMin = fgMatMin;
+       fCMatMax = fgMatMax;
     }
 
 //
@@ -178,11 +186,11 @@ void TFlukaConfigOption::WriteFlukaInputCards()
 //
 //  First make sure that all cuts are taken into account
     if (DefaultProcessFlag(kPAIR) > 0 && fProcessFlag[kPAIR] == -1 && (fCutValue[kCUTELE] >= 0. || fCutValue[kPPCUTM] >= 0.)) 
-       fProcessFlag[kPAIR] = DefaultProcessFlag(kPAIR);
+       fProcessFlag[kPAIR] = DefaultProcessFlag(kPAIR);
     if (DefaultProcessFlag(kBREM) > 0 && fProcessFlag[kBREM] == -1 && (fCutValue[kBCUTE]  >= 0. || fCutValue[kBCUTM] >= 0.)) 
-       fProcessFlag[kBREM] = DefaultProcessFlag(kBREM);
+       fProcessFlag[kBREM] = DefaultProcessFlag(kBREM);
     if (DefaultProcessFlag(kDRAY) > 0 && fProcessFlag[kDRAY] == -1 && (fCutValue[kDCUTE]  >= 0. || fCutValue[kDCUTM] >= 0.)) 
-       fProcessFlag[kDRAY] = DefaultProcessFlag(kDRAY);
+       fProcessFlag[kDRAY] = DefaultProcessFlag(kDRAY);
 //    
 //
     if (fProcessFlag[kDCAY] != -1) ProcessDCAY();
@@ -222,9 +230,9 @@ void TFlukaConfigOption::ProcessDCAY()
     // Process DCAY option
     fprintf(fgFile,"*\n* --- DCAY --- Decays. Flag = %5d\n", fProcessFlag[kDCAY]);
     if (fProcessFlag[kDCAY] == 0) {
-       printf("Decays cannot be switched off \n");
+       printf("Decays cannot be switched off \n");
     } else {
-       fprintf(fgFile, "* Decays are on by default\n");
+       fprintf(fgFile, "* Decays are on by default\n");
     }
 }
 
@@ -233,14 +241,14 @@ void TFlukaConfigOption::ProcessPAIR()
 {
     // Process PAIR option
     fprintf(fgFile,"*\n* --- PAIR --- Pair production by gammas, muons and hadrons. Flag = %5d, PPCUTM = %13.4g, PPCUTE = %13.4g\n", 
-           fProcessFlag[kPAIR], fCutValue[kCUTELE], fCutValue[kPPCUTM]);
+           fProcessFlag[kPAIR], fCutValue[kCUTELE], fCutValue[kPPCUTM]);
     //
     // gamma -> e+ e-
     //
     if (fProcessFlag[kPAIR] > 0) {
-       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.0, fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.0, fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 1e10,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 1e10,  fCMatMin, fCMatMax, 1.);
     }
     
     //
@@ -251,7 +259,7 @@ void TFlukaConfigOption::ProcessPAIR()
     //
 
     if (fProcessFlag[kBREM] == -1 ) fProcessFlag[kBREM] = fgDProcessFlag[kBREM];
-    if (fCutValue[kBCUTM]   == -1.) fCutValue[kBCUTM]   = fgDCutValue[kBCUTM]; 
+    if (fCutValue[kBCUTM]   == -1.) fCutValue[kBCUTM]   = fgDCutValue[kBCUTM];       
 
 
     Float_t flag = -3.;    
@@ -266,7 +274,7 @@ void TFlukaConfigOption::ProcessPAIR()
     // Energy cut for pair prodution
     //
     Float_t cutP = fCutValue[kPPCUTM];
-    if (fCutValue[kPPCUTM]   == -1.) cutP = fgDCutValue[kPPCUTM];      
+    if (fCutValue[kPPCUTM]   == -1.) cutP = fgDCutValue[kPPCUTM];       
     // In G3 this is the cut on the total energy of the e+e- pair
     // In FLUKA the cut is on the kinetic energy of the electron and poistron
     cutP = cutP / 2. - 0.51099906e-3;
@@ -278,12 +286,12 @@ void TFlukaConfigOption::ProcessPAIR()
     //
     Float_t cutB = 0.;
     if (flag > 1.) {
-       fprintf(fgFile,"*\n* +++ BREM --- Bremsstrahlung by muons/hadrons. Flag = %5d, BCUTM = %13.4g \n",
-           fProcessFlag[kBREM], fCutValue[kBCUTM]);
+       fprintf(fgFile,"*\n* +++ BREM --- Bremsstrahlung by muons/hadrons. Flag = %5d, BCUTM = %13.4g \n",
+           fProcessFlag[kBREM], fCutValue[kBCUTM]);
 
-       cutB =  fCutValue[kBCUTM];
-       // No explicite production of gammas
-       if (fProcessFlag[kBREM] == 2) cutB = -1.;
+       cutB =  fCutValue[kBCUTM];
+       // No explicite production of gammas
+       if (fProcessFlag[kBREM] == 2) cutB = -1.;
     }
 
     fprintf(fgFile,"PAIRBREM  %10.1f%10.4g%10.4g%10.1f%10.1f\n",flag, cutP, cutB, fCMatMin, fCMatMax);
@@ -294,26 +302,26 @@ void TFlukaConfigOption::ProcessBREM()
 {
     // Process BREM option
     fprintf(fgFile,"*\n* --- BREM --- Bremsstrahlung by e+/- and muons/hadrons. Flag = %5d, BCUTE = %13.4g, BCUTM = %13.4g \n",
-           fProcessFlag[kBREM], fCutValue[kBCUTE], fCutValue[kBCUTM]);
+           fProcessFlag[kBREM], fCutValue[kBCUTE], fCutValue[kBCUTM]);
 
     //
     // e+/- -> e+/- gamma
     //
     Float_t cutB = fCutValue[kBCUTE];
-    if (fCutValue[kBCUTE]   == -1.) cutB = fgDCutValue[kBCUTE];        
+    if (fCutValue[kBCUTE]   == -1.) cutB = fgDCutValue[kBCUTE];       
     
     
     if (TMath::Abs(fProcessFlag[kBREM]) > 0) {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",cutB,  0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",cutB,  0., 0.,  fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fELPO-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
     }
 
     //
     // Bremsstrahlung by muons and hadrons
     //
     cutB = fCutValue[kBCUTM];
-    if (fCutValue[kBCUTM]   == -1.) cutB = fgDCutValue[kBCUTM];        
+    if (fCutValue[kBCUTM]   == -1.) cutB = fgDCutValue[kBCUTM];       
     if (fProcessFlag[kBREM] == 2) cutB = -1.;
     Float_t flag = 2.;
     if (fProcessFlag[kBREM] == 0) flag = -2.;
@@ -331,9 +339,9 @@ void TFlukaConfigOption::ProcessCOMP()
     //
 
     if (fProcessFlag[kCOMP] > 0) {
-       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0.   , 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0.   , 0., 0.,  fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -347,9 +355,9 @@ void TFlukaConfigOption::ProcessPHOT()
     //
 
     if (fProcessFlag[kPHOT] > 0) {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fPHOT-THR\n",0., 0., 0.,  fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"EMFCUT    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 1.e10, 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1fPHOT-THR\n",0., 1.e10, 0.,  fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -363,9 +371,9 @@ void TFlukaConfigOption::ProcessANNI()
     //
 
     if (fProcessFlag[kANNI] > 0) {
-       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",0.   , 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",0.   , 0., 0.,  fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.1f%10.1f%10.1f%10.1f%10.1fANNH-THR\n",1.e10, 0., 0.,  fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -380,9 +388,9 @@ void TFlukaConfigOption::ProcessPFIS()
     //
 
     if (fProcessFlag[kPFIS] > 0) {
-       fprintf(fgFile,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",(Float_t) fProcessFlag[kPFIS], 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",(Float_t) fProcessFlag[kPFIS], 0., 0.,  fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0., 0.,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"PHOTONUC  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0., 0.,  fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -395,9 +403,9 @@ void TFlukaConfigOption::ProcessMUNU()
     // Muon nuclear interactions
     //
     if (fProcessFlag[kMUNU] > 0) {
-       fprintf(fgFile,"MUPHOTON  %10.1f%10.3f%10.3f%10.1f%10.1f%10.1f\n",(Float_t )fProcessFlag[kMUNU], 0.25, 0.75,  fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"MUPHOTON  %10.1f%10.3f%10.3f%10.1f%10.1f%10.1f\n",(Float_t )fProcessFlag[kMUNU], 0.25, 0.75,  fCMatMin, fCMatMax, 1.);
     } else {
-       fprintf(fgFile,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0.,   0.,    fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"MUPHOTON  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",-1.                          , 0.,   0.,    fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -414,12 +422,12 @@ void TFlukaConfigOption::ProcessRAYL()
     //
     // Loop over regions of a given material
     for (Int_t k = 0; k < nreg; k++) {
-       Float_t ireg = reglist[k];
-       if (fProcessFlag[kRAYL] > 0) {
-           fprintf(fgFile,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n", 1., ireg, ireg, 1.);
-       } else {
-           fprintf(fgFile,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n", 3., ireg, ireg, 1.);
-       }
+       Float_t ireg = reglist[k];
+       if (fProcessFlag[kRAYL] > 0) {
+           fprintf(fgFile,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n", 1., ireg, ireg, 1.);
+       } else {
+           fprintf(fgFile,"EMFRAY    %10.1f%10.1f%10.1f%10.1f\n", 3., ireg, ireg, 1.);
+       }
     }
 }
 
@@ -437,54 +445,54 @@ void TFlukaConfigOption::ProcessCKOV()
     Int_t nmaterial =  matList->GetEntriesFast();
     for (Int_t im = 0; im < nmaterial; im++)
     {
-       TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
-       Int_t idmat = material->GetIndex();
+       TGeoMaterial* material = dynamic_cast<TGeoMaterial*> (matList->At(im));
+       Int_t idmat = material->GetIndex();
 //
 // Check if global option
-       if (fMedium != -1 && idmat != fMedium) continue;
-       
-       TFlukaCerenkov* cerenkovProp;
-       if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
-       //
-       // This medium has Cerenkov properties 
-       //
-       //
-       if (fMedium == -1 || (fMedium != -1 && fProcessFlag[kCKOV] > 0)) {
-           // Write OPT-PROD card for each medium 
-           Float_t  emin  = cerenkovProp->GetMinimumEnergy();
-           Float_t  emax  = cerenkovProp->GetMaximumEnergy();        
-           fprintf(fgFile, "OPT-PROD  %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0., 
-                   Float_t(idmat), Float_t(idmat), 0.); 
-           //
-           // Write OPT-PROP card for each medium 
-           // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx)
-           //
-           fprintf(fgFile, "OPT-PROP  %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n",  
-                   cerenkovProp->GetMinimumWavelength(), cerenkovProp->GetMaximumWavelength(), cerenkovProp->GetMaximumWavelength(), 
-                   Float_t(idmat), Float_t(idmat), 0.0);
-           
-
-           fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", -100., -100., -100., 
-                   Float_t(idmat), Float_t(idmat), 0.0);
-           
-           for (Int_t j = 0; j < 3; j++) {
-               fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n", -100., -100., -100., 
-                       Float_t(idmat), Float_t(idmat), 0.0);
-           }
-
-
-           // Photon detection efficiency user defined     
-           if (cerenkovProp->IsSensitive())
-               fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n", -100., -100., -100., 
-                       Float_t(idmat), Float_t(idmat), 0.0);
-           // Material has a reflective surface
-           if (cerenkovProp->IsMetal())
-               fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n", -100., -100., -100., 
-                       Float_t(idmat), Float_t(idmat), 0.0);
-
-       } else {
-           fprintf(fgFile,"OPT-PROD  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
-       }
+       if (fMedium != -1 && idmat != fMedium) continue;
+       
+       TFlukaCerenkov* cerenkovProp;
+       if (!(cerenkovProp = dynamic_cast<TFlukaCerenkov*>(material->GetCerenkovProperties()))) continue;
+       //
+       // This medium has Cerenkov properties
+       //
+       //
+       if (fMedium == -1 || (fMedium != -1 && fProcessFlag[kCKOV] > 0)) {
+           // Write OPT-PROD card for each medium
+           Float_t  emin  = cerenkovProp->GetMinimumEnergy();
+           Float_t  emax  = cerenkovProp->GetMaximumEnergy();
+           fprintf(fgFile, "OPT-PROD  %10.4g%10.4g%10.4g%10.4g%10.4g%10.4gCERENKOV\n", emin, emax, 0.,
+                  Float_t(idmat), Float_t(idmat), 0.);
+           //
+           // Write OPT-PROP card for each medium
+           // Forcing FLUKA to call user routines (queffc.cxx, rflctv.cxx, rfrndx.cxx)
+           //
+           fprintf(fgFile, "OPT-PROP  %10.4g%10.4g%10.4g%10.1f%10.1f%10.1fWV-LIMIT\n",
+                  cerenkovProp->GetMinimumWavelength(), cerenkovProp->GetMaximumWavelength(), cerenkovProp->GetMaximumWavelength(),
+                  Float_t(idmat), Float_t(idmat), 0.0);
+       
+
+           fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", -100., -100., -100.,
+                  Float_t(idmat), Float_t(idmat), 0.0);
+       
+           for (Int_t j = 0; j < 3; j++) {
+              fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f&\n", -100., -100., -100.,
+                     Float_t(idmat), Float_t(idmat), 0.0);
+           }
+
+
+           // Photon detection efficiency user defined
+           if (cerenkovProp->IsSensitive())
+              fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fSENSITIV\n", -100., -100., -100.,
+                     Float_t(idmat), Float_t(idmat), 0.0);
+           // Material has a reflective surface
+           if (cerenkovProp->IsMetal())
+              fprintf(fgFile, "OPT-PROP  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fMETAL\n", -100., -100., -100.,
+                     Float_t(idmat), Float_t(idmat), 0.0);
+
+       } else {
+           fprintf(fgFile,"OPT-PROD  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fCERE-OFF\n",0., 0., 0., fCMatMin, fCMatMax, 1.);
+       }
     }
 }
 
@@ -495,10 +503,10 @@ void TFlukaConfigOption::ProcessHADR()
     fprintf(fgFile,"*\n* --- HADR --- Hadronic interactions. Flag = %5d\n", fProcessFlag[kHADR]);
     
     if (fProcessFlag[kHADR] > 0) {
-       fprintf(fgFile,"*\n*Hadronic interaction is ON by default in FLUKA\n");
+       fprintf(fgFile,"*\n*Hadronic interaction is ON by default in FLUKA\n");
     } else {
-       if (fMedium != -1) printf("Hadronic interactions cannot be switched off material by material !\n");
-       fprintf(fgFile,"THRESHOL  %10.1f%10.1f%10.1f%10.1e%10.1f\n",0., 0., 0., 1.e10, 0.);
+       if (fMedium != -1) printf("Hadronic interactions cannot be switched off material by material !\n");
+       fprintf(fgFile,"THRESHOL  %10.1f%10.1f%10.1f%10.1e%10.1f\n",0., 0., 0., 1.e10, 0.);
     }
 }
 
@@ -512,9 +520,9 @@ void TFlukaConfigOption::ProcessMULS()
     // Multiple scattering
     //
     if (fProcessFlag[kMULS] > 0) {
-       fprintf(fgFile,"*\n*Multiple scattering is  ON by default in FLUKA\n");
+       fprintf(fgFile,"*\n*Multiple scattering is  ON by default in FLUKA\n");
     } else {
-       fprintf(fgFile,"MULSOPT   %10.1f%10.1f%10.1f%10.1f%10.1f\n",0., 3., 3., fCMatMin, fCMatMax);
+       fprintf(fgFile,"MULSOPT   %10.1f%10.1f%10.1f%10.1f%10.1f\n",0., 3., 3., fCMatMin, fCMatMax);
     }
 }
 
@@ -522,7 +530,7 @@ void TFlukaConfigOption::ProcessLOSS()
 {
     // Process LOSS option
     fprintf(fgFile,"*\n* --- LOSS --- Ionisation energy loss. Flags: LOSS = %5d, DRAY = %5d, STRA = %5d; Cuts: DCUTE = %13.4g, DCUTM = %13.4g \n",
-           fProcessFlag[kLOSS], fProcessFlag[kDRAY], fProcessFlag[kSTRA], fCutValue[kDCUTE], fCutValue[kDCUTM]);
+           fProcessFlag[kLOSS], fProcessFlag[kDRAY], fProcessFlag[kSTRA], fCutValue[kDCUTE], fCutValue[kDCUTM]);
     //
     // Ionisation energy loss
     //
@@ -531,27 +539,27 @@ void TFlukaConfigOption::ProcessLOSS()
     
     if (fProcessFlag[kLOSS] == 1 || fProcessFlag[kLOSS] == 3 || fProcessFlag[kLOSS] > 10) {
     // Restricted fluctuations
-       fProcessFlag[kDRAY] = 1;
+       fProcessFlag[kDRAY] = 1;
     } else if (fProcessFlag[kLOSS] == 2) {
     // Full fluctuations
-       fProcessFlag[kDRAY] = 0;
-       fCutValue[kDCUTE] = 1.e10;
-       fCutValue[kDCUTM] = 1.e10;      
+       fProcessFlag[kDRAY] = 0;
+       fCutValue[kDCUTE] = 1.e10;
+       fCutValue[kDCUTM] = 1.e10;
     } else {
-       if (fProcessFlag[kDRAY] == 1) {
-           fProcessFlag[kLOSS] = 1;
-       } else if (fProcessFlag[kDRAY] == 0) {
-           fProcessFlag[kLOSS] = 2;
-           fCutValue[kDCUTE] = 1.e10;
-           fCutValue[kDCUTM] = 1.e10;  
-       }
+       if (fProcessFlag[kDRAY] == 1) {
+           fProcessFlag[kLOSS] = 1;
+       } else if (fProcessFlag[kDRAY] == 0) {
+           fProcessFlag[kLOSS] = 2;
+           fCutValue[kDCUTE] = 1.e10;
+           fCutValue[kDCUTM] = 1.e10;
+       }
     }
     
     if (fCutValue[kDCUTE] == -1.) fCutValue[kDCUTE] = fgDCutValue[kDCUTE];
     if (fCutValue[kDCUTM] == -1.) fCutValue[kDCUTM] = fgDCutValue[kDCUTM];    
     
     Float_t cutM =  fCutValue[kDCUTM];    
-       
+       
 
     if (fProcessFlag[kSTRA] == -1) fProcessFlag[kSTRA] = fgDProcessFlag[kSTRA];
     if (fProcessFlag[kSTRA] ==  0) fProcessFlag[kSTRA] = 1;
@@ -562,34 +570,34 @@ void TFlukaConfigOption::ProcessLOSS()
 //
 // Restricted energy loss fluctuations 
 //
-       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax);
-       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax);
+       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
     } else if (fProcessFlag[kLOSS] > 10) {
 //
 // Primary ionisation electron generation
 //
-       // Ionisation model
-       Float_t ioModel = Float_t (fProcessFlag[kLOSS]-10);
-       //  Effective 1st ionisation potential
-       Float_t ePot    = ModelParameter(kPRIMIOE);
-       // Number of primary ionisations per cm for a mip 
-       Float_t nPrim   = ModelParameter(kPRIMION); 
-       
-       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax);
-       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPRIM-ION\n", ePot, nPrim, ioModel, fCMatMin, fCMatMax, 1.);
-       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
+       // Ionisation model
+       Float_t ioModel = Float_t (fProcessFlag[kLOSS]-10);
+       //  Effective 1st ionisation potential
+       Float_t ePot    = ModelParameter(kPRIMIOE);
+       // Number of primary ionisations per cm for a mip
+       Float_t nPrim   = ModelParameter(kPRIMION);
+       
+       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n", 1., 1., stra, fCMatMin, fCMatMax);
+       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1fPRIM-ION\n", ePot, nPrim, ioModel, fCMatMin, fCMatMax, 1.);
+       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", cutM, 0., 0., fCMatMin, fCMatMax, 1.);
     } else if (fProcessFlag[kLOSS] == 4) {
 //
 // No fluctuations
 //
-       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., -1., stra, fCMatMin, fCMatMax);        
-       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);      
+       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",-1., -1., stra, fCMatMin, fCMatMax);
+       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
     }  else {
 //
 // Full fluctuations
 //
-       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",1., 1., stra, fCMatMin, fCMatMax);  
-       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);      
+       fprintf(fgFile,"IONFLUCT  %10.1f%10.1f%10.1f%10.1f%10.1f\n",1., 1., stra, fCMatMin, fCMatMax);
+       fprintf(fgFile,"DELTARAY  %10.4g%10.1f%10.1f%10.1f%10.1f%10.1f\n", 1.e10, 0., 0., fCMatMin, fCMatMax, 1.);
     }
 }
 
@@ -600,18 +608,18 @@ void TFlukaConfigOption::ProcessCUTGAM()
 //
     fprintf(fgFile,"*\n*Cut for Gammas. CUTGAM = %13.4g\n", fCutValue[kCUTGAM]);
     if (fMedium == -1) {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 
-               0., fCutValue[kCUTGAM], 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
+              0., fCutValue[kCUTGAM], 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
 
     } else {
-       Int_t nreg, *reglist;
-       Float_t ireg;
-       reglist = fgGeom->GetMaterialList(fMedium, nreg);
-       // Loop over regions of a given material
-       for (Int_t k = 0; k < nreg; k++) {
-           ireg = reglist[k];
-           fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.);
-       }
+       Int_t nreg, *reglist;
+       Float_t ireg;
+       reglist = fgGeom->GetMaterialList(fMedium, nreg);
+       // Loop over regions of a given material
+       for (Int_t k = 0; k < nreg; k++) {
+           ireg = reglist[k];
+           fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 0.,fCutValue[kCUTGAM], 0., ireg, ireg, 1.);
+       }
     }
     
     // Transport production cut used for pemf
@@ -620,7 +628,7 @@ void TFlukaConfigOption::ProcessCUTGAM()
     //  For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be
     Float_t parFudgem = (fCutValue[kCUTGAM] > 1.e-4)? 1.0 : 0.0 ;
     fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 
-           0., fCutValue[kCUTGAM], parFudgem, fCMatMin, fCMatMax, 1.);
+           0., fCutValue[kCUTGAM], parFudgem, fCMatMin, fCMatMax, 1.);
 }
 
 void TFlukaConfigOption::ProcessCUTELE()
@@ -629,17 +637,17 @@ void TFlukaConfigOption::ProcessCUTELE()
 //
     fprintf(fgFile,"*\n*Cut for e+/e-. CUTELE = %13.4g\n", fCutValue[kCUTELE]);
     if (fMedium == -1) {
-       fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", 
-               -fCutValue[kCUTELE], 0., 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
+       fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
+              -fCutValue[kCUTELE], 0., 0., 0., Float_t(fgGeom->NofVolumes()), 1.);
     } else {
-       Int_t nreg, *reglist;
-       Float_t ireg;
-       reglist = fgGeom->GetMaterialList(fMedium, nreg);
-       // Loop over regions of a given material
-       for (Int_t k = 0; k < nreg; k++) {
-           ireg = reglist[k];
-           fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.);
-       }
+       Int_t nreg, *reglist;
+       Float_t ireg;
+       reglist = fgGeom->GetMaterialList(fMedium, nreg);
+       // Loop over regions of a given material
+       for (Int_t k = 0; k < nreg; k++) {
+           ireg = reglist[k];
+           fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n", -fCutValue[kCUTELE], 0., 0., ireg, ireg, 1.);
+       }
     }
     // Transport production cut used for pemf
     //
@@ -647,7 +655,7 @@ void TFlukaConfigOption::ProcessCUTELE()
     //  For production and transport cut-offs larger than 100 keV it must be set = 1.0, while in the keV region it must be
     Float_t parFudgem = (fCutValue[kCUTELE] > 1.e-4)? 1.0 : 0.0;
     fprintf(fgFile,"EMFCUT    %10.4g%10.4g%10.1f%10.1f%10.1f%10.1fPROD-CUT\n", 
-           -fCutValue[kCUTELE], 0., parFudgem, fCMatMin, fCMatMax, 1.);
+           -fCutValue[kCUTELE], 0., parFudgem, fCMatMin, fCMatMax, 1.);
 }
 
 void TFlukaConfigOption::ProcessCUTNEU()
@@ -695,49 +703,49 @@ void TFlukaConfigOption::ProcessCUTNEU()
         // 9.0 = Antineutron
         fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -neutronCut,  8.0,  9.0);
         fprintf(fgFile,"LOW-BIAS  %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
-               Float_t(groupCut), 73.0, 0.95, 2., Float_t(fgGeom->NofVolumes()), 1.);
+              Float_t(groupCut), 73.0, 0.95, 2., Float_t(fgGeom->NofVolumes()), 1.);
         //
-       //
-       // 12.0 = Kaon zero long
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 12.0, 12.0);
-       // 17.0 = Lambda, 18.0 = Antilambda
-       // 19.0 = Kaon zero short
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 17.0, 19.0);
-       // 22.0 = Sigma zero, Pion zero, Kaon zero
-       // 25.0 = Antikaon zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 22.0, 25.0);
-       // 32.0 = Antisigma zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 32.0, 32.0);
-       // 34.0 = Xi zero
-       // 35.0 = AntiXi zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 34.0, 35.0);
-       // 47.0 = D zero
-       // 48.0 = AntiD zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 47.0, 48.0);
-       // 53.0 = Xi_c zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 53.0, 53.0);
-       // 55.0 = Xi'_c zero
-       // 56.0 = Omega_c zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 55.0, 56.0);
-       // 59.0 = AntiXi_c zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 59.0, 59.0);
-       // 61.0 = AntiXi'_c zero
-       // 62.0 = AntiOmega_c zero
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 61.0, 62.0);
+       //
+       // 12.0 = Kaon zero long
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 12.0, 12.0);
+       // 17.0 = Lambda, 18.0 = Antilambda
+       // 19.0 = Kaon zero short
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 17.0, 19.0);
+       // 22.0 = Sigma zero, Pion zero, Kaon zero
+       // 25.0 = Antikaon zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 22.0, 25.0);
+       // 32.0 = Antisigma zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 32.0, 32.0);
+       // 34.0 = Xi zero
+       // 35.0 = AntiXi zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 34.0, 35.0);
+       // 47.0 = D zero
+       // 48.0 = AntiD zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 47.0, 48.0);
+       // 53.0 = Xi_c zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 53.0, 53.0);
+       // 55.0 = Xi'_c zero
+       // 56.0 = Omega_c zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 55.0, 56.0);
+       // 59.0 = AntiXi_c zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 59.0, 59.0);
+       // 61.0 = AntiXi'_c zero
+       // 62.0 = AntiOmega_c zero
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 61.0, 62.0);
     } else {
         Int_t nreg, *reglist;
         Float_t ireg;
         reglist = fgGeom->GetMaterialList(fMedium, nreg);
         // Loop over regions of a given material
         for (Int_t k = 0; k < nreg; k++) {
-         ireg = reglist[k];
-         fprintf(fgFile,"LOW-BIAS  %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
-                 Float_t(groupCut), 73.0, 0.95, ireg, ireg, 1.);
-       }
-
-       Warning("ProcessCUTNEU", 
-               "Material #%4d %s: Cut on neutral hadrons (Ekin > %9.3e) material by material only implemented for low-energy neutrons !\n", 
-               fMedium, fCMaterial->GetName(), cut);
+         ireg = reglist[k];
+         fprintf(fgFile,"LOW-BIAS  %10.4g%10.4g%10.1f%10.1f%10.1f%10.1f\n",
+                Float_t(groupCut), 73.0, 0.95, ireg, ireg, 1.);
+       }
+
+       Warning("ProcessCUTNEU",
+              "Material #%4d %s: Cut on neutral hadrons (Ekin > %9.3e) material by material only implemented for low-energy neutrons !\n",
+              fMedium, fCMaterial->GetName(), cut);
     }
 }
 
@@ -747,39 +755,39 @@ void TFlukaConfigOption::ProcessCUTHAD()
     fprintf(fgFile,"*\n*Cut for charge hadrons. CUTHAD = %13.4g\n", fCutValue[kCUTHAD]);
     Float_t cut = fCutValue[kCUTHAD];
     if (fMedium == -1) {
-       // 1.0 = Proton
-       // 2.0 = Antiproton
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut,  1.0,  2.0);
-       // 13.0 = Positive Pion, Negative Pion, Positive Kaon
-       // 16.0 = Negative Kaon
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 13.0, 16.0);
-       // 20.0 = Negative Sigma
-       // 21.0 = Positive Sigma
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 20.0, 21.0);
-       // 31.0 = Antisigma minus
-       // 33.0 = Antisigma plus
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 31.0, 31.0);
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 33.0, 33.0);
-       // 36.0 = Negative Xi, Positive Xi, Omega minus
-       // 39.0 = Antiomega
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 36.0, 39.0);
-       // 45.0 = D plus
-       // 46.0 = D minus
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 45.0, 46.0);
-       // 49.0 = D_s plus, D_s minus, Lambda_c plus
-       // 52.0 = Xi_c plus
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 49.0, 52.0);
-       // 54.0 = Xi'_c plus
-       // 60.0 = AntiXi'_c minus
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 54.0, 54.0);
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 60.0, 60.0);
-       // 57.0 = Antilambda_c minus
-       // 58.0 = AntiXi_c minus
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0);
+       // 1.0 = Proton
+       // 2.0 = Antiproton
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut,  1.0,  2.0);
+       // 13.0 = Positive Pion, Negative Pion, Positive Kaon
+       // 16.0 = Negative Kaon
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 13.0, 16.0);
+       // 20.0 = Negative Sigma
+       // 21.0 = Positive Sigma
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 20.0, 21.0);
+       // 31.0 = Antisigma minus
+       // 33.0 = Antisigma plus
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 31.0, 31.0);
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 33.0, 33.0);
+       // 36.0 = Negative Xi, Positive Xi, Omega minus
+       // 39.0 = Antiomega
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 36.0, 39.0);
+       // 45.0 = D plus
+       // 46.0 = D minus
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 45.0, 46.0);
+       // 49.0 = D_s plus, D_s minus, Lambda_c plus
+       // 52.0 = Xi_c plus
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 49.0, 52.0);
+       // 54.0 = Xi'_c plus
+       // 60.0 = AntiXi'_c minus
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 54.0, 54.0);
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 60.0, 60.0);
+       // 57.0 = Antilambda_c minus
+       // 58.0 = AntiXi_c minus
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n", -cut, 57.0, 58.0);
     } else {
       Warning("ProcessCUTHAD", 
-             "Material #%4d %s: Cut on charged hadrons (Ekin > 9.3e) material by material not yet implemented !\n", 
-             fMedium, fCMaterial->GetName(), cut); 
+              "Material #%4d %s: Cut on charged hadrons (Ekin > %9.3e) material by material not yet implemented !\n",
+             fMedium, fCMaterial->GetName(), cut);
     }
 }
 
@@ -789,10 +797,10 @@ void TFlukaConfigOption::ProcessCUTMUO()
     fprintf(fgFile,"*\n*Cut for muons. CUTMUO = %13.4g\n", fCutValue[kCUTMUO]);
     Float_t cut = fCutValue[kCUTMUO];
     if (fMedium == -1) {
-       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n",-cut, 10.0, 11.0);
+       fprintf(fgFile,"PART-THR  %10.4g%10.1f%10.1f\n",-cut, 10.0, 11.0);
     } else {
-       Warning("ProcessCUTMUO", "Material #%4d %s: Cut on muons (Ekin > %9.3e) material by material not yet implemented !\n", 
-               fMedium, fCMaterial->GetName(), cut);
+       Warning("ProcessCUTMUO", "Material #%4d %s: Cut on muons (Ekin > %9.3e) material by material not yet implemented !\n",
+              fMedium, fCMaterial->GetName(), cut);
     }
     
     
index 09db597..ab1b58b 100644 (file)
@@ -41,7 +41,7 @@ public:
     // Getters
     Double_t Cut(FlukaCutOption_t i)              const {return fCutValue[i];}
     Int_t    Flag(FlukaProcessOption_t i)         const {return fProcessFlag[i];}
-    Int_t    ModelParameter(FlukaModelOption_t i) const {return fModelParameter[i];}
+    Double_t  ModelParameter(FlukaModelOption_t i) const {return fModelParameter[i];}
     Int_t    Medium()                             const {return fMedium;}    
     // Setters
     void     SetCut(const char* flagname, Double_t val);
@@ -79,6 +79,8 @@ public:
        {fgFile = file; fgMatMin = matMin; fgMatMax = matMax; fgGeom = geom;}
     static Double_t DefaultCut(FlukaCutOption_t i) {return fgDCutValue[i];}
     static Int_t    DefaultProcessFlag(FlukaProcessOption_t i) {return fgDProcessFlag[i];}
+    
  protected:
     Double_t fCutValue[11];            // User cut
     Int_t    fProcessFlag[15];         // User flag assigned to processes
@@ -95,6 +97,12 @@ public:
     static Float_t   fgMatMax;            // Maximum meterial number
     static FILE*     fgFile;              // Output file
     static TFlukaMCGeometry* fgGeom;      // Pointer to geometry     
+
+ private:
+    // Copy constructor and operator= declared but not implemented (-Weff++ flag)
+    TFlukaConfigOption(const TFlukaConfigOption&);
+    TFlukaConfigOption& operator=(const TFlukaConfigOption&);
+    
     ClassDef(TFlukaConfigOption, 1)    // Fluka Configuration Option
 };
        
index d73dc9c..a170f77 100644 (file)
@@ -33,6 +33,7 @@
 #include "TObjString.h"
 #include "Fsourcm.h"
 #include "Ftrackr.h"
+#include "Fstepsz.h"       //(STEPSZ) fluka common
 
 #ifndef WIN32 
 # define idnrwr idnrwr_
@@ -79,11 +80,11 @@ extern "C"
    void  type_of_call   g1wr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/, 
                              Double_t * /*pV*/,  Int_t & /*oldReg*/ , const Int_t & /*oldLttc*/, Double_t & /*propStep*/,
                              Int_t & /*nascFlag*/, Double_t & /*retStep*/, Int_t & /*newReg*/,
-                                 Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/,
+                                  Double_t & /*saf*/, Int_t & /*newLttc*/, Int_t & /*LttcFlag*/,
                              Double_t *s /*Lt*/, Int_t * /*jrLt*/);
    
    void  type_of_call g1rtwr();
-   void  type_of_call conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/); 
+   void  type_of_call conhwr(Int_t & /*intHist*/, Int_t & /*incrCount*/);
    void  type_of_call inihwr(Int_t & /*intHist*/);
    void  type_of_call jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/,
                              Int_t & /*flukaReg*/);
@@ -95,17 +96,17 @@ extern "C"
                              Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
    void  type_of_call lkmgwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
-                                      Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
+                                       Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
    void  type_of_call   lkwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
                              Double_t * /*pV*/, const Int_t & /*oldReg*/, const Int_t & /*oldLttc*/,
-                                 Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
+                                  Int_t & /*flagErr*/, Int_t & /*newReg*/, Int_t & /*newLttc*/);
 //   void  type_of_call magfld(const Double_t & /*pX*/, const Double_t & /*pY*/, const Double_t & /*pZ*/,
 //                             Double_t & /*cosBx*/, Double_t & /*cosBy*/, Double_t & /*cosBz*/, 
-//                             Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/);       
+//                             Double_t & /*Bmag*/, Int_t & /*reg*/, Int_t & /*idiscflag*/);        
    void  type_of_call nrmlwr(Double_t & /*pSx*/, Double_t & /*pSy*/, Double_t & /*pSz*/,
                              Double_t & /*pVx*/, Double_t & /*pVy*/, Double_t & /*pVz*/,
-                                 Double_t * /*norml*/, const Int_t & /*oldReg*/, 
-                                 const Int_t & /*newReg*/, Int_t & /*flagErr*/);
+                                  Double_t * /*norml*/, const Int_t & /*oldReg*/,
+                                  const Int_t & /*newReg*/, Int_t & /*flagErr*/);
    void  type_of_call rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
                              Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/);
    Int_t type_of_call isvhwr(const Int_t & /*fCheck*/, const Int_t & /*intHist*/);
@@ -122,47 +123,49 @@ TFlukaMCGeometry* TFlukaMCGeometry::fgInstance= NULL;
 
 //_____________________________________________________________________________
 TFlukaMCGeometry::TFlukaMCGeometry(const char *name, const char *title) 
-                 : TNamed(name, title)
+   :TNamed(name, title), 
+   fDebug(kFALSE),
+   fLastMaterial(0),
+   fDummyRegion(0),
+   fCurrentRegion(0),
+   fCurrentLattice(0),
+   fNextRegion(0),
+   fNextLattice(0),
+   fRegionList(0),
+   fIndmat(0),
+   fMatList(new TObjArray(256)),
+   fMatNames(new TObjArray(256))
 {
   //
   // Standard constructor
   //
-  fDebug        = kFALSE;
-  fLastMaterial = 0;
-  fCurrentRegion   = 0;
-  fCurrentLattice  = 0;
-  fDummyRegion  = 0;
-  fNextRegion   = 0;
-  fNextLattice  = 0;
-  fRegionList   = 0;
-  fIndmat = 0;
   gFluka = (TFluka*)gMC;
   gMCGeom = this;
   gNstep = 0;
-  fMatList = new TObjArray(256);
-  fMatNames = new TObjArray(256);
 }
 
 //_____________________________________________________________________________
 TFlukaMCGeometry::TFlukaMCGeometry()
-{    
+  :TNamed(),
+   fDebug(kFALSE),
+   fLastMaterial(0),
+   fDummyRegion(0),
+   fCurrentRegion(0),
+   fCurrentLattice(0),
+   fNextRegion(0),
+   fNextLattice(0),
+   fRegionList(0),
+   fIndmat(0),
+   fMatList(0),
+   fMatNames(0)
+
+{
   //
   // Default constructor
   //
-  fDebug        = kFALSE;
-  fLastMaterial = 0;
-  fCurrentRegion   = 0;
-  fCurrentLattice  = 0;
-  fDummyRegion  = 0;
-  fNextRegion   = 0;
-  fNextLattice  = 0;
-  fRegionList   = 0;
-  fIndmat = 0;
   gFluka = (TFluka*)gMC;
   gMCGeom = this;
   gNstep = 0;
-  fMatList = 0;
-  fMatNames = 0;
 }
 
 //_____________________________________________________________________________
@@ -181,14 +184,6 @@ TFlukaMCGeometry::~TFlukaMCGeometry()
 //
 // private methods
 //
-//_____________________________________________________________________________
-TFlukaMCGeometry::TFlukaMCGeometry(const TFlukaMCGeometry &)
-  : TNamed()
-{    
-  //
-  // Copy constructor
-  //
-}
 
 //_____________________________________________________________________________
 Double_t* TFlukaMCGeometry::CreateDoubleArray(Float_t* array, Int_t size) const
@@ -253,12 +248,12 @@ Int_t *TFlukaMCGeometry::GetRegionList(Int_t imed, Int_t &nreg)
        if ((med = vol->GetMedium()) == 0) continue;
        imedium = med->GetId();
        if (imedium == imed) {
-          ireg = vol->GetNumber();
-          fRegionList[nreg++] = ireg;
+           ireg = vol->GetNumber();
+           fRegionList[nreg++] = ireg;
        }
    }
    return fRegionList;
-}         
+}
 
 //_____________________________________________________________________________
 Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg)
@@ -274,12 +269,12 @@ Int_t *TFlukaMCGeometry::GetMaterialList(Int_t imat, Int_t &nreg)
        if ((med = vol->GetMedium()) == 0) continue;
        imaterial = med->GetMaterial()->GetIndex();
        if (imaterial == imat) {
-          ireg = vol->GetNumber();
-          fRegionList[nreg++] = ireg;
+           ireg = vol->GetNumber();
+           fRegionList[nreg++] = ireg;
        }
    }
    return fRegionList;
-}         
+}
  
 //_____________________________________________________________________________
 Int_t TFlukaMCGeometry::NofVolumes() const 
@@ -827,32 +822,32 @@ void TFlukaMCGeometry::CreatePemfFile()
     TString sname;
     
     for (i = fIndmat; i < fLastMaterial - 2; i++) {
-       printf("Write Peg Files %d\n", i);
-       
-       mat = (TGeoMaterial*)fMatList->At(i);
-       if (!mat->IsUsed()) continue;
-       sname = "mat";
-       sprintf(number, "%d", i);
-       sname.Append(number);
-       cout << endl;
-       cout << endl;
-       cout << "******************************************************************************" << endl;
-       cout << "******************************************************************************" << endl;
-       cout << endl;
-       WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas);
-       sname.Prepend("$FLUPRO/pemf/rpemf peg/");
-       gSystem->Exec(sname.Data());
-       
-       // check if the pemf file was created
-       TString sname = Form("peg/mat%d.pemf", i);
-       ifstream in( sname.Data() );
-       if ( in ) {
-           countMatOK++;
-           in.close();
-       } else {
-           cout << "ERROR Fail to create the pemf file " << sname << endl;
-           countPemfError++; 
-       }
+        printf("Write Peg Files %d\n", i);
+        
+        mat = (TGeoMaterial*)fMatList->At(i);
+        if (!mat->IsUsed()) continue;
+        sname = "mat";
+        sprintf(number, "%d", i);
+        sname.Append(number);
+        cout << endl;
+        cout << endl;
+        cout << "******************************************************************************" << endl;
+        cout << "******************************************************************************" << endl;
+        cout << endl;
+        WritePegFile(i, &countNoStern, &countElemError, &countMixError, &countGas);
+        sname.Prepend("$FLUPRO/pemf/rpemf peg/");
+        gSystem->Exec(sname.Data());
+        
+        // check if the pemf file was created
+        TString sname = Form("peg/mat%d.pemf", i);
+        ifstream in( sname.Data() );
+        if ( in ) {
+            countMatOK++;
+            in.close();
+        } else {
+            cout << "ERROR Fail to create the pemf file " << sname << endl;
+            countPemfError++;
+        }
     }
     cout << "Materials (pemf created)   " << countMatOK         << endl;
     cout << "Not Sternheimer par. found  " << countNoStern   << endl;
@@ -1025,11 +1020,11 @@ void TFlukaMCGeometry::WritePegFile(Int_t imat, Int_t *NoStern, Int_t *ElemError
    while((proc = (TFlukaConfigOption*)next()))
    { 
        if (proc->Medium() == mat->GetIndex()) {
-          ap = proc->Cut(kCUTGAM);
-          ae = proc->Cut(kCUTELE);
-          if (ap == -1.) ap =  TFlukaConfigOption::DefaultCut(kCUTGAM);
-          if (ae == -1.) ae =  TFlukaConfigOption::DefaultCut(kCUTELE);
-          break;
+           ap = proc->Cut(kCUTGAM);
+           ae = proc->Cut(kCUTELE);
+           if (ap == -1.) ap =  TFlukaConfigOption::DefaultCut(kCUTGAM);
+           if (ae == -1.) ae =  TFlukaConfigOption::DefaultCut(kCUTELE);
+           break;
        }
    }
 
@@ -1073,6 +1068,10 @@ Double_t * TFlukaMCGeometry::GetISSB(Double_t rho, Int_t nElem, Double_t *zelem,
       Double_t    afact;              // AFACT
       Double_t    sk;                 // SK
       Double_t    delta0;             // DELTA0
+
+      sternheimerData():
+         longname(""), nelems(0), density(0), iev(0), cbar(0),
+         x0(0), x1(0), afact(0), sk(0), delta0(0) {}
    };
    
    TString     shortname;
@@ -1291,9 +1290,9 @@ void TFlukaMCGeometry::FlukaMatName(TString &str) const
     TObjArray * tokens = str.Tokenize("_");
     Int_t ntok = tokens->GetEntries();
     if (ntok > 0) {
-       TString head = ((TObjString*) tokens->At(0))->GetString();
-       Int_t nhead = head.Length();
-       str = str.Remove(0, nhead + 1);
+        TString head = ((TObjString*) tokens->At(0))->GetString();
+        Int_t nhead = head.Length();
+        str = str.Remove(0, nhead + 1);
     }
     tokens->Clear();
     delete tokens;
@@ -1326,7 +1325,7 @@ void TFlukaMCGeometry::FlukaMatName(TString &str) const
       }   
    }   
 }   
-         
+
 //______________________________________________________________________________
 void TFlukaMCGeometry::Vname(const char *name, char *vname) const
 {
@@ -1342,6 +1341,12 @@ void TFlukaMCGeometry::Vname(const char *name, char *vname) const
 }
 
 
+//______________________________________________________________________________
+Int_t TFlukaMCGeometry::GetNstep()
+{
+   // return gNstep for debug propose
+   return gNstep;
+}
 
 // FLUKA GEOMETRY WRAPPERS - to replace FLUGG wrappers
 
@@ -1363,12 +1368,19 @@ Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/)
 void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz, 
           Double_t *pV,  Int_t &oldReg , const Int_t &oldLttc, Double_t &propStep,
           Int_t &/*nascFlag*/, Double_t &retStep, Int_t &newReg,
-              Double_t &saf, Int_t &newLttc, Int_t &lttcFlag,
+          Double_t &saf, Int_t &newLttc, Int_t &lttcFlag,
           Double_t *sLt, Int_t *jrLt)
 
 {
    // Initialize FLUKa point and direction;
    gNstep++;
+
+//   if( (gNstep > 43912170 && gNstep < 43912196 ) ||
+//       (gNstep > 47424560 && gNstep < 47424581  ) ||
+//       (gNstep > 54388266 && gNstep < 54388319 )
+//       ) gMCGeom->SetDebugMode();
+//   else gMCGeom->SetDebugMode(kFALSE);
+   
    NORLAT.xn[0] = pSx;
    NORLAT.xn[1] = pSy;
    NORLAT.xn[2] = pSz;
@@ -1377,10 +1389,17 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
    if (oldLttc<=0) {
       gGeoManager->FindNode(pSx,pSy,pSz);
       olttc = gGeoManager->GetCurrentNodeId()+1;
-   }   
+      oldReg = gGeoManager->GetCurrentVolume()->GetNumber();
+   }
+
+   if (gMCGeom->IsDebugging()) cout << "g1wr gNstep=" << gNstep
+                                    << " oldReg="<< oldReg <<" olttc="<< olttc
+                                    << " track=" << TRACKR.ispusr[mkbmx2-1] << endl;
+
    Int_t ccreg,cclat;
    gMCGeom->GetCurrentRegion(ccreg,cclat);
-   Bool_t crossed = (ccreg==oldReg && cclat==oldLttc)?kFALSE:kTRUE;
+   Bool_t crossed = (ccreg==oldReg && cclat==olttc)?kFALSE:kTRUE;
+
    gMCGeom->SetCurrentRegion(oldReg, olttc);
    // Initialize default return values
    lttcFlag = 0;
@@ -1404,12 +1423,13 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
       sLt[lttcFlag] = TGeoShape::Tolerance(); // null step in next region
       jrLt[lttcFlag+1] = -1;
       sLt[lttcFlag+1] = 0.; // null step in next region;
+      if (gMCGeom->IsDebugging()) printf("=> crossed dummy!! newReg=%i newLttc=%i\n", newReg, newLttc);
       return;
-   }   
-      
+   }
+
    // Reset outside flag
    gGeoManager->SetOutside(kFALSE);
-   
+
    curLttc = gGeoManager->GetCurrentNodeId()+1;
    curReg = gGeoManager->GetCurrentVolume()->GetNumber();
    if (olttc != curLttc) {
@@ -1418,12 +1438,12 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
       gGeoManager->CdNode(olttc-1);
       curLttc = gGeoManager->GetCurrentNodeId()+1;
       curReg  = gGeoManager->GetCurrentVolume()->GetNumber();
-   }  
+   }
    // Now the current TGeo state reflects the FLUKA state 
-     
+
    gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
    gGeoManager->SetCurrentDirection(pV);
-   
+
    if (crossed) {
       gGeoManager->FindNextBoundaryAndStep(propStep);
       saf = 0.0;
@@ -1432,8 +1452,8 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
       saf = gGeoManager->GetSafeDistance();
       if (saf<0) saf=0.0;
       saf -= saf*3.0e-09;
-   }   
-      
+   }
+
    Double_t snext = gGeoManager->GetStep();
 
    if (snext<=0.0) snext = TGeoShape::Tolerance();
@@ -1445,7 +1465,7 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
    NORLAT.xn[2] += snext*pV[2];
    if (!gGeoManager->IsOnBoundary()) {
    // Next boundary further than proposed step, which is approved
-      if (saf>propStep) saf = propStep;       
+      if (saf>propStep) saf = propStep;
       retStep = propStep;
       sLt[lttcFlag] = propStep;
       return;
@@ -1454,51 +1474,61 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
    gGeoManager->SetCurrentPoint(pSx,pSy,pSz);
    newLttc = (gGeoManager->IsOutside())?(TFlukaMCGeometry::kLttcOutside):gGeoManager->GetCurrentNodeId()+1;
    newReg = (gGeoManager->IsOutside())?(gMCGeom->NofVolumes()+2):gGeoManager->GetCurrentVolume()->GetNumber();
-   if (gMCGeom->IsDebugging()) printf("   newReg=%i newLttc=%i\n", newReg, newLttc);
+   if (gMCGeom->IsDebugging()) printf("=> newReg=%i newLttc=%i\n", newReg, newLttc);
 
    // We really crossed the boundary, but is it the same region ?
    gMCGeom->SetNextRegion(newReg, newLttc);
 
    Int_t pid = TRACKR.jtrack;
-   if (newReg==oldReg && newLttc!=olttc && pid!=-1) {
+   if ( ((newReg==oldReg && newLttc!=olttc) || (oldReg!=newReg && olttc==newLttc) ) && pid!=-1) {
       // Virtual boundary between replicants
       newReg  = gFluka->GetDummyRegion();
       newLttc = TFlukaMCGeometry::kLttcVirtual;
-   } 
-   
+      if (gMCGeom->IsDebugging()) printf("=> virtual boundary!! newReg=%i newLttc=%i\n", newReg, newLttc);
+   }
+
    retStep = snext;
    sLt[lttcFlag] = snext;
    lttcFlag++;
    jrLt[lttcFlag] = newLttc;
    sLt[lttcFlag] = snext;
    jrLt[lttcFlag+1] = -1;
-   sLt[lttcFlag+1] = 0.;      
+
+   sLt[lttcFlag+1] = 0.;
    gGeoManager->SetOutside(kFALSE);
    gGeoManager->CdNode(olttc-1);
    if (gMCGeom->IsDebugging()) {
       printf("=> snext=%g safe=%g\n", snext, saf);
       for (Int_t i=0; i<lttcFlag+1; i++) printf("   jrLt[%i]=%i  sLt[%i]=%g\n", i,jrLt[i],i,sLt[i]);
-   }   
+   }
 }
 
 //_____________________________________________________________________________
 void g1rtwr()
 {
+
    if (gMCGeom->IsDebugging()) printf("========== Dummy G1RTWR\n");
 } 
 
 //_____________________________________________________________________________
-void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/)
+void conhwr(Int_t & intHist, Int_t & incrCount)
 {
-   if (gMCGeom->IsDebugging()) printf("========== Dummy CONHWR\n");
+   if (gMCGeom->IsDebugging()) printf("========== Dummy CONHWR intHist=%d incrCount=%d currentNodeId=%d\n",
+                                       intHist, incrCount, gGeoManager->GetCurrentNodeId()+1 );
+//   if( incrCount != -1 ) {
+//      if (intHist==0) gGeoManager->CdTop();
+//      else gGeoManager->CdNode(intHist-1);
+//   }
+//   intHist = gGeoManager->GetCurrentNodeId()+1;
 }
 
 //_____________________________________________________________________________
 void inihwr(Int_t &intHist)
 {
-   if (gMCGeom->IsDebugging()) printf("========== Inside INIHWR -> reinitializing history: %i\n", intHist);
+   if (gMCGeom->IsDebugging())
+      printf("========== Inside INIHWR -> reinitializing history: %i \n", intHist);
    if (gGeoManager->IsOutside()) gGeoManager->CdTop();
-   if (intHist<=0) {
+   if (intHist<0) {
 //      printf("=== wrong history number\n");
       return;
    }
@@ -1507,7 +1537,7 @@ void inihwr(Int_t &intHist)
    if (gMCGeom->IsDebugging()) {
       printf(" --- current path: %s\n", gGeoManager->GetPath());
       printf("<= INIHWR\n");
-   }   
+   }
 }
 
 //_____________________________________________________________________________
@@ -1551,7 +1581,7 @@ void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
 //_____________________________________________________________________________
 void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
             Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
-                     Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
+            Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
 {
    if (gMCGeom->IsDebugging()) {
       printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
@@ -1564,7 +1594,7 @@ void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
 //_____________________________________________________________________________
 void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
           Double_t *pV, const Int_t &oldReg, const Int_t &oldLttc,
-              Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
+          Int_t &flagErr, Int_t &newReg, Int_t &newLttc)
 {
    if (gMCGeom->IsDebugging()) {
       printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
@@ -1592,15 +1622,27 @@ void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
       if (newreg1==newReg && newlttc1==newLttc) {
          newReg = dummy;
          newLttc = TFlukaMCGeometry::kLttcVirtual;
-      }   
+         flagErr = newReg;
+         if (gMCGeom->IsDebugging()) printf("  virtual boundary (oldReg==dummy) !! newReg=%i newLttc=%i\n", newReg, newLttc);
+      }
       return;
-   }   
+   }
 
    if (oldReg==newReg && oldLttc!=newLttc) {
-      newReg  = gFluka->GetDummyRegion();
+      newReg  = dummy;
       newLttc = TFlukaMCGeometry::kLttcVirtual;
-   }   
-         
+      if (gMCGeom->IsDebugging()) printf("  virtual boundary!! newReg=%i newLttc=%i\n", newReg, newLttc);
+   }
+
+   if( oldReg!=newReg && oldLttc==newLttc ) {
+      // this should not happen!! ??? Ernesto
+       cout << "  lkwr      oldReg!=newReg ("<< oldReg <<"!="<< newReg
+            << ") && oldLttc==newLttc ("<< newLttc <<") !!!!!!!!!!!!!!!!!" << endl;
+      newReg  = dummy;
+      newLttc = TFlukaMCGeometry::kLttcVirtual;
+      flagErr = newReg;
+   }
+
    if (gMCGeom->IsDebugging()) {
       printf("  LKWR: newReg=%i newLttc=%i\n", newReg, newLttc);
    }   
@@ -1609,11 +1651,12 @@ void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
 //_____________________________________________________________________________
 void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
             Double_t &pVx, Double_t &pVy, Double_t &pVz,
-                Double_t *norml, const Int_t &oldReg, 
-                const Int_t &newReg, Int_t &flagErr)
+            Double_t *norml, const Int_t &oldReg,
+            const Int_t &newReg, Int_t &flagErr)
 {
    if (gMCGeom->IsDebugging()) {
       printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
+      printf("                         (%g, %g, %g)\n", NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
       printf("   oldReg=%i, newReg=%i\n", oldReg,newReg);
    }   
    gGeoManager->SetCurrentPoint(NORLAT.xn[0], NORLAT.xn[1], NORLAT.xn[2]);
@@ -1627,15 +1670,16 @@ void nrmlwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
       norml[1] = -pVy;   
       norml[2] = -pVz; 
    } else {
-      norml[0] = -dnorm[0];   
-      norml[1] = -dnorm[1];   
-      norml[2] = -dnorm[2]; 
+      norml[0] = -dnorm[0];
+      norml[1] = -dnorm[1];
+      norml[2] = -dnorm[2];
    }  
-       
+
    if (gMCGeom->IsDebugging()) {
       printf("   normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);  
       printf("<= NRMLWR\n");
-   }   
+   }
+
 }
 
 //_____________________________________________________________________________
@@ -1661,7 +1705,7 @@ Int_t isvhwr(const Int_t &check, const Int_t & intHist)
 
 // For TGeo, just return the current node ID. No copy need to be made.
 
-   if (gMCGeom->IsDebugging()) printf("=> Inside ISVHWR\n");
+   if (gMCGeom->IsDebugging()) printf("=> Inside ISVHWR check=%d intHist=%d\n", check, intHist);
    if (check<0) return intHist;
    Int_t histInt = gGeoManager->GetCurrentNodeId()+1;
    if (gMCGeom->IsDebugging()) printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath());
index 9588e1e..1ef55fa 100644 (file)
@@ -31,6 +31,7 @@ class TFlukaMCGeometry :public TNamed {
     virtual ~TFlukaMCGeometry();
   
     // get methods
+    Int_t         GetNstep();  // to be removed
     Int_t         GetMedium() const;
     Int_t        *GetRegionList(Int_t imed, Int_t &nreg);
     Int_t        *GetMaterialList(Int_t imat, Int_t &nreg);
@@ -46,6 +47,7 @@ class TFlukaMCGeometry :public TNamed {
     void          SetMreg(Int_t mreg, Int_t lttc);
     void          SetCurrentRegion(Int_t mreg, Int_t latt);
     void          GetCurrentRegion(Int_t &mreg, Int_t &latt) const {mreg=fCurrentRegion; latt=fCurrentLattice;}
+    Int_t         GetCurrentRegion() const {return fCurrentRegion;}
     Int_t         GetDummyRegion() const {return fDummyRegion;}
     Int_t         GetDummyLattice() const {return kLttcVirtual;}
     void          SetNextRegion(Int_t mreg, Int_t latt);
@@ -65,8 +67,9 @@ class TFlukaMCGeometry :public TNamed {
     void     Vname(const char *name, char *vname) const;
    
   private:
+    // Copy constructor and operator= declared but not implemented (-Weff++ flag)
     TFlukaMCGeometry(const TFlukaMCGeometry& rhs);
-    TFlukaMCGeometry& operator=(const TFlukaMCGeometry& /*rhs*/) {return (*this);}
+    TFlukaMCGeometry& operator=(const TFlukaMCGeometry& /*rhs*/); // {return (*this);}
 
     static TFlukaMCGeometry*  fgInstance; // singleton instance
     Bool_t       fDebug;                  // debug flag
@@ -80,6 +83,7 @@ class TFlukaMCGeometry :public TNamed {
     Int_t        fIndmat;                 // material index where pemf file creation starts
     TObjArray   *fMatList;                //! material list as known by FLUKA
     TObjArray   *fMatNames;               //! list of FLUKA material names
+    
   ClassDef(TFlukaMCGeometry,1)  //Virtual MonteCarlo Interface
 };
 
index e45b006..62a2290 100644 (file)
@@ -26,12 +26,19 @@ FILE*              TFlukaScoringOption::fgFile(0x0);
 TFlukaMCGeometry*  TFlukaScoringOption::fgGeom(0x0);
 
 TFlukaScoringOption::TFlukaScoringOption()
+   : fNopfl(0),
+     fOutFile(""),
+     fLun(0)
 {
+
     // Default constructor
 }
 
 TFlukaScoringOption::TFlukaScoringOption(const char* name, const char* sdum, Int_t nopfl, char* outfile, Float_t* what)
-    : TNamed(name, sdum)
+    : TNamed(name, sdum),
+      fNopfl(nopfl),
+      fOutFile(outfile),
+      fLun(0)
 {
     // Constructor
     fNopfl   = nopfl;
@@ -41,12 +48,15 @@ TFlukaScoringOption::TFlukaScoringOption(const char* name, const char* sdum, Int
 }
 
 TFlukaScoringOption::TFlukaScoringOption(const char* name, const char* sdum, Int_t nopfl, char* outfile, Float_t* what, 
-                                        const char* det1, const char* det2, const char* det3)
-    : TNamed(name, sdum)
+                                         const char* det1, const char* det2, const char* det3)
+    : TNamed(name, sdum),
+      fNopfl(nopfl + 2),
+      fOutFile(outfile),
+      fLun(0)
 {
     // Constructor
-    fNopfl   = nopfl + 2;
-    fOutFile = outfile;
+//    fNopfl   = nopfl + 2;
+//    fOutFile = outfile;
     Int_t npar = (nopfl == 0)? 6 : 12;
     for (Int_t i = 0; i < npar; i++)  fWhat[i] = what[i];
 
@@ -73,7 +83,7 @@ void TFlukaScoringOption::WriteOpenFlukaFile()
     // Write Fluka input card for output file opening
     // 
     fprintf(fgFile, "OPEN      %10.1f                                                    %s\n", 
-           GetLun(), "NEW");     
+            GetLun(), "NEW");
     fprintf(fgFile, "%s\n", GetFileName());
 }
 
@@ -97,63 +107,63 @@ void TFlukaScoringOption::WriteFlukaInputCards()
 // USRBIN
 //
     if(strncmp(GetName(), "USRBIN", 6) == 0){
-       if (Par() == 0) {
-           fprintf(fgFile, "USRBIN    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", 
-                   What(1), What(2), GetLun(), What(4), What(5), What(6));
-       } else if (Par() == 1) {
-           fprintf(fgFile, "USRBIN    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",  
-                   What(1), What(2), GetLun(), What(4), What(5),  What(6));
-           fprintf(fgFile, "USRBIN    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1f  %s\n",  
-                   What(7), What(8), What(9), What(10), What(11), What(12), cont_line);
-       } else if (Par() == 2) {
-           if(What(1) == 2.0 || What(1) == 12){
-               fprintf(fgFile, "USRBIN    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",  
-                       What(1), What(2), GetLun(), Float_t(GetRegionByName(GetRegName(1))),
-                       Float_t(GetRegionByName(GetRegName(2))), Float_t(GetRegionByName(GetRegName(3))));
-               fprintf(fgFile, "USRBIN    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1f  %s\n",  
-                       Float_t(GetRegionByName(GetRegName(1))), Float_t(GetRegionByName(GetRegName(2))),
-                       Float_t(GetRegionByName(GetRegName(3))), 1., 1., 1., cont_line);
-           } else {
-               printf("Check consistency of SetUserScoring values \n");
-           }
-       }
+        if (Par() == 0) {
+            fprintf(fgFile, "USRBIN    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",
+                    What(1), What(2), GetLun(), What(4), What(5), What(6));
+        } else if (Par() == 1) {
+            fprintf(fgFile, "USRBIN    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",
+                    What(1), What(2), GetLun(), What(4), What(5),  What(6));
+            fprintf(fgFile, "USRBIN    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1f  %s\n",
+                    What(7), What(8), What(9), What(10), What(11), What(12), cont_line);
+        } else if (Par() == 2) {
+            if(What(1) == 2.0 || What(1) == 12){
+                fprintf(fgFile, "USRBIN    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",
+                        What(1), What(2), GetLun(), Float_t(GetRegionByName(GetRegName(1))),
+                        Float_t(GetRegionByName(GetRegName(2))), Float_t(GetRegionByName(GetRegName(3))));
+                fprintf(fgFile, "USRBIN    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1f  %s\n",
+                        Float_t(GetRegionByName(GetRegName(1))), Float_t(GetRegionByName(GetRegName(2))),
+                        Float_t(GetRegionByName(GetRegName(3))), 1., 1., 1., cont_line);
+            } else {
+                printf("Check consistency of SetUserScoring values \n");
+            }
+        }
     }
    
 //
 // USRBDX
 //
     if(strncmp(GetName(), "USRBDX", 6) == 0){
-       if (Par() == 2) {
-           fprintf(fgFile, "USRBDX    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n", 
-                   What(1), What(2), GetLun(), Float_t(GetRegionByName(GetRegName(1))), 
-                   Float_t(GetRegionByName(GetRegName(2))), What(6));
-       } else if (Par() == 3) {
-           fprintf(fgFile, "USRBDX    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",  
-                   What(1), What(2), GetLun(), Float_t(GetRegionByName(GetRegName(1))),
-                   Float_t(GetRegionByName(GetRegName(2))), What(6));
-           fprintf(fgFile, "USRBDX    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1f  %s\n",  
-                   What(7), What(8), What(9), What(10), What(11), What(12), cont_line);              
-       }
+        if (Par() == 2) {
+            fprintf(fgFile, "USRBDX    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",
+                    What(1), What(2), GetLun(), Float_t(GetRegionByName(GetRegName(1))),
+                    Float_t(GetRegionByName(GetRegName(2))), What(6));
+        } else if (Par() == 3) {
+            fprintf(fgFile, "USRBDX    %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",
+                    What(1), What(2), GetLun(), Float_t(GetRegionByName(GetRegName(1))),
+                    Float_t(GetRegionByName(GetRegName(2))), What(6));
+            fprintf(fgFile, "USRBDX    %10.1f%10.4g%10.1f%10.1f%10.1f%10.1f  %s\n",
+                    What(7), What(8), What(9), What(10), What(11), What(12), cont_line);
+        }
     }
     
 //
 // USTRACK
 //
     if(strncmp(GetName(), "USRTRACK", 6) == 0){
-       fprintf(fgFile, "USRTRACK  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",  
-               What(1), What(2), GetLun(), Float_t(GetRegionByName(GetRegName(1))), What(4), What(5));
-       fprintf(fgFile, "USRTRACK  %10.1f%10.4g                                          %s\n",  
-               What(7), What(8), cont_line);
+        fprintf(fgFile, "USRTRACK  %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",
+                What(1), What(2), GetLun(), Float_t(GetRegionByName(GetRegName(1))), What(4), What(5));
+        fprintf(fgFile, "USRTRACK  %10.1f%10.4g                                          %s\n",
+                What(7), What(8), cont_line);
     }
     
 //
 // USRCOLL
 // 
     if(strncmp(GetName(), "USRCOLL", 6) == 0){
-       fprintf(fgFile, "USRCOLL   %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",  
-               What(1), What(2), GetLun(), Float_t(GetRegionByName(GetRegName(1))), What(4), What(5));
-       fprintf(fgFile, "USRCOLL   %10.1f%10.4g                                          %s\n",  
-               What(7), What(8), cont_line);         
+        fprintf(fgFile, "USRCOLL   %10.1f%10.1f%10.1f%10.1f%10.1f%10.1f\n",
+                What(1), What(2), GetLun(), Float_t(GetRegionByName(GetRegName(1))), What(4), What(5));
+        fprintf(fgFile, "USRCOLL   %10.1f%10.4g                                          %s\n",
+                What(7), What(8), cont_line);
     }
 }
 
index bc876a9..53fa537 100644 (file)
@@ -29,24 +29,25 @@ public:
     TFlukaScoringOption();
     TFlukaScoringOption(const char* name, const char* sdum, Int_t nopfl, char* outfile, Float_t* what);
     TFlukaScoringOption(const char* name, const char* sdum, Int_t nopfl, char* outfile, Float_t* what,
-                       const char* det1, const char* det2, const char* det3);
+                        const char* det1, const char* det2, const char* det3);
     // Getters
-    Float_t What(Int_t indi) const       {return fWhat[indi - 1];}
-    Int_t   Par() const                  {return fNopfl;}
-    char*   GetFileName() const          {return fOutFile;}
-    Float_t GetLun()  const              {return fLun;}
+    Float_t     What(Int_t indi) const       {return fWhat[indi - 1];}
+    Int_t       Par() const                  {return fNopfl;}
+    char*       GetFileName() const          {return fOutFile;}
+    Float_t     GetLun()  const              {return fLun;}
     const char* GetRegName(Int_t ndet);
-    void SetPar(Int_t val)                 {fNopfl   = val;}
-    void SetFileName(char* outfile)        {fOutFile = outfile;}
-    void SetLun(Float_t lun)               {fLun     = lun;}
-           
+    void        SetPar(Int_t val)            {fNopfl   = val;}
+    void        SetFileName(char* outfile)   {fOutFile = outfile;}
+    void        SetLun(Float_t lun)          {fLun     = lun;}
 //
-    void  WriteFlukaInputCards();
-    void  WriteOpenFlukaFile();
-    Int_t GetRegionByName(const char* detname);
-  
+    void        WriteFlukaInputCards();
+    void        WriteOpenFlukaFile();
+    Int_t       GetRegionByName(const char* detname);
+
     static void SetStaticInfo(FILE* file, TFlukaMCGeometry* geom)
-       {fgFile = file; fgGeom = geom;}
+                {fgFile = file; fgGeom = geom;}
+
  protected:
     Int_t         fNopfl;        // Number of paramters
     Float_t       fWhat[12];     // WHAT()
@@ -57,7 +58,12 @@ public:
     // Static
     static FILE*             fgFile;      // Output file
     static TFlukaMCGeometry* fgGeom;      // Pointer to geometry
-    
+ private:
+    // Copy constructor and operator= declared but not implemented (-Weff++ flag)
+    TFlukaScoringOption(const TFlukaScoringOption&);
+    TFlukaScoringOption& operator=(const TFlukaScoringOption&);
+
     ClassDef(TFlukaScoringOption, 1)  // Fluka Scoring Option
 };
        
index c9abc12..931b8dd 100644 (file)
 #else
 # define bxdraw BXDRAW
 #endif
+
+
+#include "TGeoManager.h"
+
+
 extern "C" {
 void bxdraw(Int_t& icode, Int_t& mreg, Int_t& newreg,
             Double_t& xsco, Double_t& ysco, Double_t& zsco)
 {
     TFluka* fluka = (TFluka*) gMC;
-    Int_t oldlttc = LTCLCM.mlatm1;
+    Int_t oldlttc = TRACKR.lt1trk; //LTCLCM.mlatm1;
     Int_t newlttc = LTCLCM.newlat;
     fluka->SetIcode((FlukaProcessCode_t)icode);
-    fluka->SetNewreg(newreg,newlttc);
+//    fluka->SetNewreg(newreg,newlttc);
     fluka->SetXsco(xsco);
     fluka->SetYsco(ysco);
     fluka->SetZsco(zsco);
     Int_t verbosityLevel = fluka->GetVerbosityLevel();
     Bool_t debug = (verbosityLevel>=3)?kTRUE:kFALSE;
+    // nothing to do if particle is crossing a dummy region
+    if (mreg   == fluka->GetDummyRegion() ||
+        newreg == fluka->GetDummyRegion() ||
+        oldlttc == TFlukaMCGeometry::kLttcVirtual ||
+        newlttc == TFlukaMCGeometry::kLttcVirtual
+        ) return;
+
 //
 // Double step for boundary crossing
 //
     fluka->SetTrackIsNew(kFALSE); // has to be called BEFORE Stepping()
-    if (mreg != fluka->GetDummyRegion()) {
+    if (mreg != fluka->GetDummyRegion() && newreg != fluka->GetDummyRegion()) {
        if (debug) printf("bxdraw (ex) \n");
        fluka->SetTrackIsExiting();
        fluka->SetCaller(kBXExiting);
        fluka->SetMreg(mreg,oldlttc);
+
+    // check region lattice consistency (debug Ernesto)
+    // *****************************************************
+    Int_t nodeId;
+    Int_t volId = fluka->CurrentVolID(nodeId);
+    Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
+    if( mreg != volId  && !gGeoManager->IsOutside()) {
+       cout << "  bxdraw:   track=" << TRACKR.ispusr[mkbmx2-1]<< " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+            << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
+            << "               fluka   mreg=" << mreg << " oldlttc=" << oldlttc << " newreg=" << newreg << " newlttc=" << newlttc << endl
+            << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
+            << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
+            << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
+            << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
+            << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
+        if( oldlttc == crtlttc ) cout << "   **************************** Exit *********************************" << endl;
+    }
+    // *****************************************************
+
+
+       
        TVirtualMCStack* cppstack = fluka->GetStack();
        cppstack->SetCurrentTrack( TRACKR.ispusr[mkbmx2-1] );
        (TVirtualMCApplication::Instance())->Stepping();
@@ -44,6 +77,25 @@ void bxdraw(Int_t& icode, Int_t& mreg, Int_t& newreg,
        fluka->SetTrackIsEntering();
        if (fluka->GetDummyBoundary() == 1) fluka->SetDummyBoundary(2);
        fluka->SetMreg(newreg,newlttc);
+
+    // check region lattice consistency (debug Ernesto)
+    // *****************************************************
+    Int_t nodeId;
+    Int_t volId = fluka->CurrentVolID(nodeId);
+    Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
+    if( newreg != volId  && !gGeoManager->IsOutside()) {
+       cout << "  bxdraw:   track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+            << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
+            << "               fluka   mreg=" << mreg << " oldlttc=" << oldlttc << " newreg=" << newreg << " newlttc=" << newlttc << endl
+            << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
+            << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
+            << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
+            << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
+            << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
+        if( newlttc == crtlttc ) cout << "   ******************************** Enter *****************************" << endl;
+    }
+    // *****************************************************
+
        TVirtualMCStack* cppstack = fluka->GetStack();
        cppstack->SetCurrentTrack( TRACKR.ispusr[mkbmx2-1] );
        (TVirtualMCApplication::Instance())->Stepping();
index a955075..593c122 100644 (file)
@@ -26,14 +26,14 @@ void endraw(Int_t& icode, Int_t& mreg, Double_t& rull, Double_t& xsco, Double_t&
   if (mreg == fluka->GetDummyRegion()) return;
   Int_t verbosityLevel = fluka->GetVerbosityLevel();
   Bool_t debug = (verbosityLevel >= 3)? kTRUE : kFALSE;
-  Int_t mlttc = LTCLCM.mlatm1;
+  Int_t mlttc = (icode==kKASKADinelarecoil) ? TRACKR.lt2trk : TRACKR.lt1trk; //LTCLCM.mlatm1;
   fluka->SetCaller(kENDRAW);
   fluka->SetRull(rull);
   fluka->SetXsco(xsco);
   fluka->SetYsco(ysco);
   fluka->SetZsco(zsco);
   fluka->SetMreg(mreg, mlttc);
-  
+
   Float_t edep = rull;
   
   if (TRACKR.jtrack == -1) {
@@ -50,29 +50,45 @@ void endraw(Int_t& icode, Int_t& mreg, Double_t& rull, Double_t& xsco, Double_t&
   }
 
   TVirtualMCStack* cppstack = fluka->GetStack();
-  Int_t saveTrackId = cppstack->GetCurrentTrackNumber();
 
   if (debug) {
-     cout << "ENDRAW For icode=" << icode << " stacktrack=" << saveTrackId
+     cout << "ENDRAW For icode=" << icode 
           << " track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
-          << " edep="<< edep <<endl;
+          << " edep="<< edep << " mreg=" << mreg << endl;
   }
 
+    // check region lattice consistency (debug Ernesto)
+    // *****************************************************
+    Int_t nodeId;
+    Int_t volId = fluka->CurrentVolID(nodeId);
+    Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
+    if( mreg != volId  && !gGeoManager->IsOutside() ) {
+       cout << "  endraw:   track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+            << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
+            << "               fluka   mreg=" << mreg << " mlttc=" << mlttc << endl
+            << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
+            << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
+            << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
+            << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
+            << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
+        if( mlttc == crtlttc ) cout << "   *************************************************************" << endl;
+    }
+    // *****************************************************
+
+
   if (icode != kEMFSCOstopping1 && icode != kEMFSCOstopping2) {
       fluka->SetIcode((FlukaProcessCode_t)icode);
       fluka->SetRull(edep);
       if (icode == kKASKADelarecoil && TRACKR.ispusr[mkbmx2-5]) {
-         //  Elastic recoil and in stuprf npprmr > 0,
-         //  the secondary being loaded is actually still the interacting particle
-         cppstack->SetCurrentTrack( TRACKR.ispusr[mkbmx2-4] );
-         //      cout << "endraw elastic recoil track=" << TRACKR.ispusr[mkbmx2-1] << " parent=" << TRACKR.ispusr[mkbmx2-4]
-         //           << endl;
+         //  Elastic recoil and in stuprf npprmr > 0,
+         //  the secondary being loaded is actually still the interacting particle
+         cppstack->SetCurrentTrack( TRACKR.ispusr[mkbmx2-4] );
+         //      cout << "endraw elastic recoil track=" << TRACKR.ispusr[mkbmx2-1] << " parent=" << TRACKR.ispusr[mkbmx2-4]
+         //           << endl;
       }
       else
-         cppstack->SetCurrentTrack(TRACKR.ispusr[mkbmx2-1] );
+          cppstack->SetCurrentTrack(TRACKR.ispusr[mkbmx2-1] );
       (TVirtualMCApplication::Instance())->Stepping();
-      
-//      cppstack->SetCurrentTrack( saveTrackId );
   } else {
   //
   // For icode 21,22 the particle has fallen below thresshold.
index 5e89412..74c9cd6 100644 (file)
@@ -19,6 +19,9 @@
 # define mgdraw MGDRAW
 #endif
 
+
+#include "TGeoManager.h" // <- delete
+
 extern "C" {
 void mgdraw(Int_t& icode, Int_t& mreg)
 {
@@ -42,11 +45,11 @@ void mgdraw(Int_t& icode, Int_t& mreg)
     Int_t verbosityLevel = fluka->GetVerbosityLevel();
 
     if (TRACKR.jtrack < -6) {
-       // (Unknow heavy Ion???)
-       // assing parent id ???
+       // from -7 to -12 = "heavy" fragment
+       // assing parent id
        // id < -6 was skipped in stuprf =>   if (kpart < -6) return;
        if (verbosityLevel >= 3) {
-          cout << "mgdraw: jtrack < -6 =" << TRACKR.jtrack
+          cout << "mgdraw: (heavy fragment) jtrack < -6 =" << TRACKR.jtrack
                << " assign parent pdg=" << fluka->PDGFromId(TRACKR.ispusr[mkbmx2 - 3]) << endl;
        }
        TRACKR.jtrack = TRACKR.ispusr[mkbmx2 - 3];
@@ -55,19 +58,47 @@ void mgdraw(Int_t& icode, Int_t& mreg)
     cppstack->SetCurrentTrack(trackId);
 //
 //    
-    Int_t mlttc = LTCLCM.mlatm1;
+    Int_t mlttc = TRACKR.lt1trk; // LTCLCM.mlatm1;
     fluka->SetMreg(mreg, mlttc);
-    fluka->SetNewreg(mreg, mlttc);
+//    fluka->SetNewreg(mreg, mlttc); // dont used!!
     fluka->SetIcode((FlukaProcessCode_t) icode);
     fluka->SetCaller(kMGDRAW);
 
+    Int_t nodeId;
+    Int_t volId = fluka->CurrentVolID(nodeId);
+    Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
+
+//
+//    if( (fluka->GetNstep() > 43912170 && fluka->GetNstep() < 43912196 ) ||
+//        (fluka->GetNstep() > 47424560 && fluka->GetNstep() < 47424581 ) ||
+//        (fluka->GetNstep() > 54388266 && fluka->GetNstep() < 54388319 )
+//      ) fluka->SetVerbosityLevel(3);
+//    else fluka->SetVerbosityLevel(0);
+
+    // check region lattice consistency (debug Ernesto)
+    // *****************************************************
+    if( mreg != volId  && !gGeoManager->IsOutside() ) {
+       cout << "  mgdraw:   track=" << trackId << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+            << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
+            << "               fluka   mreg=" << mreg << " mlttc=" << mlttc << endl
+            << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
+            << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
+            << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
+            << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
+            << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
+        if( mlttc == crtlttc ) cout << "   *************************************************************" << endl;
+    }
+    // *****************************************************
+
     if (!TRACKR.ispusr[mkbmx2 - 2]) {
         //
         // Single step
         if (verbosityLevel >= 3) {
            cout << endl << "mgdraw: energy deposition for:" << trackId
                  << " icode=" << icode
-                 << " pdg=" << fluka->PDGFromId(TRACKR.jtrack) << " flukaid="<< TRACKR.jtrack << endl;
+                 << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+                 << " flukaid="<< TRACKR.jtrack
+                 << " mreg=" << mreg << endl;
         }
         (TVirtualMCApplication::Instance())->Stepping();
         fluka->SetTrackIsNew(kFALSE);
index 1da1c03..6c6f014 100644 (file)
@@ -43,7 +43,7 @@ extern "C" {
   //
   void type_of_call geocrs(Double_t &, Double_t &, Double_t &);
   void type_of_call georeg(Double_t &, Double_t &, Double_t &, 
-                          Int_t &, Int_t &);
+                           Int_t &, Int_t &);
   void type_of_call geohsm(Int_t &, Int_t &, Int_t &, Int_t &);
   void type_of_call soevsv();
  /*
@@ -95,24 +95,26 @@ extern "C" {
     lastParticleWasPrimary = particleIsPrimary;
     
     if (itrack >= nprim) {
-       particleIsPrimary = kFALSE;
+        particleIsPrimary = kFALSE;
     } else {
-       particleIsPrimary = kTRUE;
+        particleIsPrimary = kTRUE;
     }
 
     if (lfirst) {
-       SOURCM.tkesum = zerzer;
-       lfirst = false;
-       SOURCM.lussrc = true;
+        SOURCM.tkesum = zerzer;
+        lfirst = false;
+        SOURCM.lussrc = true;
     } else {
 //
 // Post-track actions for primary track
 //
-       if (particleIsPrimary) {
-           TVirtualMCApplication::Instance()->PostTrack();
-           TVirtualMCApplication::Instance()->FinishPrimary();
-           if ((itrack%10)==0) printf("=== TRACKING PRIMARY %d ===\n", itrack);
-       }
+        if (particleIsPrimary) {
+            TVirtualMCApplication::Instance()->PostTrack();
+            TVirtualMCApplication::Instance()->FinishPrimary();
+            if ((itrack%10)==0)
+                cout << "=== TRACKING PRIMARY "<< itrack <<" ===" << endl;
+            //printf("=== TRACKING PRIMARY %d ===\n", itrack);
+        }
     }
 
     // Exit if itrack is negative (-1). Set lsouit to false to mark last track for this event
@@ -131,11 +133,11 @@ extern "C" {
     //
     // Handle user event abortion
     if (fluka->EventIsStopped()) {
-       printf("Event has been stopped by user !");
-       fluka->SetStopEvent(kFALSE);
-       nomore = 1;
-       SOURCM.lsouit = false;
-       return;
+        printf("Event has been stopped by user !");
+        fluka->SetStopEvent(kFALSE);
+        nomore = 1;
+        SOURCM.lsouit = false;
+        return;
     }
 
     //Get some info about the particle and print it
@@ -148,14 +150,14 @@ extern "C" {
        cout << "\t* Particle " << itrack << " retrieved..." << endl;
        cout << "\t\t+ Name = " << particle->GetName() << endl;
        cout << "\t\t+ PDG/Fluka code = " << pdg 
-           << " / " << fluka->IdFromPDG(pdg) << endl;
+            << " / " << fluka->IdFromPDG(pdg) << endl;
        cout << "\t\t+ P = (" 
-           << particle->Px() << " , "
-           << particle->Py() << " , "
-           << particle->Pz() << " ) --> "
-           << particle->P() << " GeV " 
-           << particle->Energy() << " GeV "  
-           << particle->GetMass() << " GeV " << endl;
+            << particle->Px() << " , "
+            << particle->Py() << " , "
+            << particle->Pz() << " ) --> "
+            << particle->P() << " GeV "
+            << particle->Energy() << " GeV "
+            << particle->GetMass() << " GeV " << endl;
     }   
     /* Npflka is the stack counter: of course any time source is called it
      * must be =0
@@ -167,163 +169,163 @@ extern "C" {
     if (particle->Pz() < 0.) cosz = -cosz;    
 
     if (pdg != 50000050 &&  pdg !=  50000051) {
-       FLKSTK.npflka++;
-       Int_t ifl =  fluka-> IdFromPDG(pdg);
-       FLKSTK.iloflk[FLKSTK.npflka] = ifl;
-       /* Wtflk is the weight of the particle*/
-       FLKSTK.wtflk[FLKSTK.npflka] = oneone;
-       SUMCOU.weipri += FLKSTK.wtflk[FLKSTK.npflka];
-       
-       FLKSTK.loflk[FLKSTK.npflka] = 1;
-       
-       /* User dependent flag:*/
-       FLKSTK.louse[FLKSTK.npflka] = 0;
-
-       /* User dependent spare variables:*/
-       Int_t ispr = 0;
-       for (ispr = 0; ispr < mkbmx1; ispr++)
-           FLKSTK.sparek[FLKSTK.npflka][ispr] = zerzer;
-
-       /* User dependent spare flags:*/
-       for (ispr = 0; ispr < mkbmx2; ispr++)
-           FLKSTK.ispark[FLKSTK.npflka][ispr] = 0;
-
-       /* Save the track number of the stack particle:*/
-       FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = itrack;
-       FLKSTK.nparma++;
-       FLKSTK.numpar[FLKSTK.npflka] = FLKSTK.nparma;
-       FLKSTK.nevent[FLKSTK.npflka] = 0;
-       FLKSTK.dfnear[FLKSTK.npflka] = +zerzer;
-       
-       /* Particle age (s)*/
-       FLKSTK.agestk[FLKSTK.npflka] = +zerzer;
-       FLKSTK.cmpath[FLKSTK.npflka] = +zerzer;
-       FLKSTK.aknshr[FLKSTK.npflka] = -twotwo;
-
-       /* Group number for "low" energy neutrons, set to 0 anyway*/
-       FLKSTK.igroup[FLKSTK.npflka] = 0;
-       
-       /* Kinetic energy */
-       Double_t p    = particle->P();
-       Double_t mass = PAPROP.am[ifl + 6];
-       FLKSTK.tkeflk[FLKSTK.npflka] =  TMath::Sqrt( p * p + mass * mass) - mass;
-       /* Particle momentum*/
-       FLKSTK.pmoflk [FLKSTK.npflka] = p;
-
-       FLKSTK.txflk [FLKSTK.npflka] = cosx;
-       FLKSTK.tyflk [FLKSTK.npflka] = cosy;
-       FLKSTK.tzflk [FLKSTK.npflka] = cosz;
-    
-       /* Polarization cosines:*/
-       if (polarisation.Mag()) {
-           Double_t cospolx = polarisation.Px() / polarisation.Mag();
-           Double_t cospoly = polarisation.Py() / polarisation.Mag();
-           Double_t cospolz = sqrt(oneone - cospolx * cospolx - cospoly * cospoly);
-           FLKSTK.txpol [FLKSTK.npflka] = cospolx;
-           FLKSTK.typol [FLKSTK.npflka] = cospoly;
-           FLKSTK.tzpol [FLKSTK.npflka] = cospolz;
-       }
-       else {
-           FLKSTK.txpol [FLKSTK.npflka] = -twotwo;
-           FLKSTK.typol [FLKSTK.npflka] = +zerzer;
-           FLKSTK.tzpol [FLKSTK.npflka] = +zerzer;
-       }
-       
-       /* Particle coordinates*/
-       // Vertext coordinates;
-       FLKSTK.xflk [FLKSTK.npflka] = particle->Vx();
-       FLKSTK.yflk [FLKSTK.npflka] = particle->Vy();
-       FLKSTK.zflk [FLKSTK.npflka] = particle->Vz();
+        FLKSTK.npflka++;
+        Int_t ifl =  fluka-> IdFromPDG(pdg);
+        FLKSTK.iloflk[FLKSTK.npflka] = ifl;
+        /* Wtflk is the weight of the particle*/
+        FLKSTK.wtflk[FLKSTK.npflka] = oneone;
+        SUMCOU.weipri += FLKSTK.wtflk[FLKSTK.npflka];
+        
+        FLKSTK.loflk[FLKSTK.npflka] = 1;
+        
+        /* User dependent flag:*/
+        FLKSTK.louse[FLKSTK.npflka] = 0;
+
+        /* User dependent spare variables:*/
+        Int_t ispr = 0;
+        for (ispr = 0; ispr < mkbmx1; ispr++)
+            FLKSTK.sparek[FLKSTK.npflka][ispr] = zerzer;
+
+        /* User dependent spare flags:*/
+        for (ispr = 0; ispr < mkbmx2; ispr++)
+            FLKSTK.ispark[FLKSTK.npflka][ispr] = 0;
+
+        /* Save the track number of the stack particle:*/
+        FLKSTK.ispark[FLKSTK.npflka][mkbmx2-1] = itrack;
+        FLKSTK.nparma++;
+        FLKSTK.numpar[FLKSTK.npflka] = FLKSTK.nparma;
+        FLKSTK.nevent[FLKSTK.npflka] = 0;
+        FLKSTK.dfnear[FLKSTK.npflka] = +zerzer;
+        
+        /* Particle age (s)*/
+        FLKSTK.agestk[FLKSTK.npflka] = +zerzer;
+        FLKSTK.cmpath[FLKSTK.npflka] = +zerzer;
+        FLKSTK.aknshr[FLKSTK.npflka] = -twotwo;
+
+        /* Group number for "low" energy neutrons, set to 0 anyway*/
+        FLKSTK.igroup[FLKSTK.npflka] = 0;
+        
+        /* Kinetic energy */
+        Double_t p    = particle->P();
+        Double_t mass = PAPROP.am[ifl + 6];
+        FLKSTK.tkeflk[FLKSTK.npflka] =  TMath::Sqrt( p * p + mass * mass) - mass;
+        /* Particle momentum*/
+        FLKSTK.pmoflk [FLKSTK.npflka] = p;
+
+        FLKSTK.txflk [FLKSTK.npflka] = cosx;
+        FLKSTK.tyflk [FLKSTK.npflka] = cosy;
+        FLKSTK.tzflk [FLKSTK.npflka] = cosz;
     
-       /*  Calculate the total kinetic energy of the primaries: don't change*/
-       Int_t st_ilo =  FLKSTK.iloflk[FLKSTK.npflka];
-       if ( st_ilo != 0 )
-           SOURCM.tkesum += 
-               ((FLKSTK.tkeflk[FLKSTK.npflka] + PAPROP.amdisc[st_ilo+6])
-                * FLKSTK.wtflk[FLKSTK.npflka]);
-       else
-           SOURCM.tkesum += (FLKSTK.tkeflk[FLKSTK.npflka] * FLKSTK.wtflk[FLKSTK.npflka]);
-       
-       /*  Here we ask for the region number of the hitting point.
-        *     NRGFLK (LFLKSTK) = ...
-        *  The following line makes the starting region search much more
-        *  robust if particles are starting very close to a boundary:
-        */
-       geocrs( FLKSTK.txflk[FLKSTK.npflka], 
-               FLKSTK.tyflk[FLKSTK.npflka], 
-               FLKSTK.tzflk[FLKSTK.npflka] );
+        /* Polarization cosines:*/
+        if (polarisation.Mag()) {
+            Double_t cospolx = polarisation.Px() / polarisation.Mag();
+            Double_t cospoly = polarisation.Py() / polarisation.Mag();
+            Double_t cospolz = sqrt(oneone - cospolx * cospolx - cospoly * cospoly);
+            FLKSTK.txpol [FLKSTK.npflka] = cospolx;
+            FLKSTK.typol [FLKSTK.npflka] = cospoly;
+            FLKSTK.tzpol [FLKSTK.npflka] = cospolz;
+        }
+        else {
+            FLKSTK.txpol [FLKSTK.npflka] = -twotwo;
+            FLKSTK.typol [FLKSTK.npflka] = +zerzer;
+            FLKSTK.tzpol [FLKSTK.npflka] = +zerzer;
+        }
+
+        /* Particle coordinates*/
+        // Vertext coordinates;
+        FLKSTK.xflk [FLKSTK.npflka] = particle->Vx();
+        FLKSTK.yflk [FLKSTK.npflka] = particle->Vy();
+        FLKSTK.zflk [FLKSTK.npflka] = particle->Vz();
+
+        /*  Calculate the total kinetic energy of the primaries: don't change*/
+        Int_t st_ilo =  FLKSTK.iloflk[FLKSTK.npflka];
+        if ( st_ilo != 0 )
+            SOURCM.tkesum +=
+                ((FLKSTK.tkeflk[FLKSTK.npflka] + PAPROP.amdisc[st_ilo+6])
+                 * FLKSTK.wtflk[FLKSTK.npflka]);
+        else
+            SOURCM.tkesum += (FLKSTK.tkeflk[FLKSTK.npflka] * FLKSTK.wtflk[FLKSTK.npflka]);
+
+        /*  Here we ask for the region number of the hitting point.
+         *     NRGFLK (LFLKSTK) = ...
+         *  The following line makes the starting region search much more
+         *  robust if particles are starting very close to a boundary:
+         */
+        geocrs( FLKSTK.txflk[FLKSTK.npflka],
+                FLKSTK.tyflk[FLKSTK.npflka],
+                FLKSTK.tzflk[FLKSTK.npflka] );
     
-       Int_t idisc;
-
-       georeg ( FLKSTK.xflk[FLKSTK.npflka], 
-                FLKSTK.yflk[FLKSTK.npflka], 
-                FLKSTK.zflk[FLKSTK.npflka],
-                FLKSTK.nrgflk[FLKSTK.npflka], 
-                idisc);//<-- dummy return variable not used
-       /*  Do not change these cards:*/
-       Int_t igeohsm1 = 1;
-       Int_t igeohsm2 = -11;
-       geohsm ( FLKSTK.nhspnt[FLKSTK.npflka], igeohsm1, igeohsm2, LTCLCM.mlattc );
-       FLKSTK.nlattc[FLKSTK.npflka] = LTCLCM.mlattc;
-       soevsv();
+        Int_t idisc;
+
+        georeg ( FLKSTK.xflk[FLKSTK.npflka],
+                 FLKSTK.yflk[FLKSTK.npflka],
+                 FLKSTK.zflk[FLKSTK.npflka],
+                 FLKSTK.nrgflk[FLKSTK.npflka],
+                 idisc);//<-- dummy return variable not used
+        /*  Do not change these cards:*/
+        Int_t igeohsm1 = 1;
+        Int_t igeohsm2 = -11;
+        geohsm ( FLKSTK.nhspnt[FLKSTK.npflka], igeohsm1, igeohsm2, LTCLCM.mlattc );
+        FLKSTK.nlattc[FLKSTK.npflka] = LTCLCM.mlattc;
+        soevsv();
     } else {
         //
-       // Next particle is optical photon
-       //
-       OPPHST.lstopp++;
-       OPPHST.donear [OPPHST.lstopp - 1] = 0.;
-
-       OPPHST.xoptph [OPPHST.lstopp - 1] = particle->Vx();
-       OPPHST.yoptph [OPPHST.lstopp - 1] = particle->Vy();
-       OPPHST.zoptph [OPPHST.lstopp - 1] = particle->Vz();
-
-       OPPHST.txopph [OPPHST.lstopp - 1] = cosx;
-       OPPHST.tyopph [OPPHST.lstopp - 1] = cosy;
-       OPPHST.tzopph [OPPHST.lstopp - 1] = cosz;
-
-
-       if (polarisation.Mag()) {
-           Double_t cospolx = polarisation.Px() / polarisation.Mag();
-           Double_t cospoly = polarisation.Py() / polarisation.Mag();
-           Double_t cospolz = sqrt(oneone - cospolx * cospolx - cospoly * cospoly);
-           OPPHST.txpopp [OPPHST.lstopp - 1] = cospolx;
-           OPPHST.typopp [OPPHST.lstopp - 1] = cospoly;
-           OPPHST.tzpopp [OPPHST.lstopp - 1] = cospolz;
-       }
-       else {
-           OPPHST.txpopp [OPPHST.lstopp - 1] = -twotwo;
-           OPPHST.typopp [OPPHST.lstopp - 1] = +zerzer;
-           OPPHST.tzpopp [OPPHST.lstopp - 1] = +zerzer;
-       }
-       
-       geocrs( OPPHST.txopph[OPPHST.lstopp - 1], 
-               OPPHST.tyopph[OPPHST.lstopp - 1], 
-               OPPHST.tzopph[OPPHST.lstopp - 1] );
-    
-       Int_t idisc;
-
-       georeg ( OPPHST.xoptph[OPPHST.lstopp - 1], 
-                OPPHST.yoptph[OPPHST.lstopp - 1], 
-                OPPHST.zoptph[OPPHST.lstopp - 1],
-                OPPHST.nregop[OPPHST.lstopp - 1], 
-                idisc);//<-- dummy return variable not used
-       
-       OPPHST.wtopph [OPPHST.lstopp - 1] = particle->GetWeight();
-       OPPHST.poptph [OPPHST.lstopp - 1] = particle->P();
-       OPPHST.agopph [OPPHST.lstopp - 1] = particle->T();      
-       OPPHST.cmpopp [OPPHST.lstopp - 1] = +zerzer;
-       OPPHST.loopph [OPPHST.lstopp - 1] = 0;
-       OPPHST.louopp [OPPHST.lstopp - 1] = itrack;
-    }
-    
+        // Next particle is optical photon
+        //
+        OPPHST.lstopp++;
+        OPPHST.donear [OPPHST.lstopp - 1] = 0.;
+
+        OPPHST.xoptph [OPPHST.lstopp - 1] = particle->Vx();
+        OPPHST.yoptph [OPPHST.lstopp - 1] = particle->Vy();
+        OPPHST.zoptph [OPPHST.lstopp - 1] = particle->Vz();
+
+        OPPHST.txopph [OPPHST.lstopp - 1] = cosx;
+        OPPHST.tyopph [OPPHST.lstopp - 1] = cosy;
+        OPPHST.tzopph [OPPHST.lstopp - 1] = cosz;
+
+
+        if (polarisation.Mag()) {
+            Double_t cospolx = polarisation.Px() / polarisation.Mag();
+            Double_t cospoly = polarisation.Py() / polarisation.Mag();
+            Double_t cospolz = sqrt(oneone - cospolx * cospolx - cospoly * cospoly);
+            OPPHST.txpopp [OPPHST.lstopp - 1] = cospolx;
+            OPPHST.typopp [OPPHST.lstopp - 1] = cospoly;
+            OPPHST.tzpopp [OPPHST.lstopp - 1] = cospolz;
+        }
+        else {
+            OPPHST.txpopp [OPPHST.lstopp - 1] = -twotwo;
+            OPPHST.typopp [OPPHST.lstopp - 1] = +zerzer;
+            OPPHST.tzpopp [OPPHST.lstopp - 1] = +zerzer;
+        }
+
+        geocrs( OPPHST.txopph[OPPHST.lstopp - 1],
+                OPPHST.tyopph[OPPHST.lstopp - 1],
+                OPPHST.tzopph[OPPHST.lstopp - 1] );
+
+        Int_t idisc;
+
+        georeg ( OPPHST.xoptph[OPPHST.lstopp - 1],
+                 OPPHST.yoptph[OPPHST.lstopp - 1],
+                 OPPHST.zoptph[OPPHST.lstopp - 1],
+                 OPPHST.nregop[OPPHST.lstopp - 1],
+                 idisc);//<-- dummy return variable not used
+
+        OPPHST.wtopph [OPPHST.lstopp - 1] = particle->GetWeight();
+        OPPHST.poptph [OPPHST.lstopp - 1] = particle->P();
+        OPPHST.agopph [OPPHST.lstopp - 1] = particle->T();
+        OPPHST.cmpopp [OPPHST.lstopp - 1] = +zerzer;
+        OPPHST.loopph [OPPHST.lstopp - 1] = 0;
+        OPPHST.louopp [OPPHST.lstopp - 1] = itrack;
+        OPPHST.nlatop [OPPHST.lstopp - 1] = LTCLCM.mlattc;
+     }
+
 //
 //  Pre-track actions at for primary tracks
 //
     if (particleIsPrimary) {
-       TVirtualMCApplication::Instance()->BeginPrimary();
-       TVirtualMCApplication::Instance()->PreTrack();
+        TVirtualMCApplication::Instance()->BeginPrimary();
+        TVirtualMCApplication::Instance()->PreTrack();
     }
-    
 //
     if (debug) cout << "<== source(" << nomore << ")" << endl;
   }
index 027497f..b1f5455 100644 (file)
@@ -52,13 +52,13 @@ void stupre()
 
   Int_t npnw, ispr;
   for (npnw = EMFSTK.npstrt-1; npnw <= EMFSTK.npemf-1; npnw++) {
-      
+
       for (ispr = 0; ispr <= mkbmx1-1; ispr++) 
-         EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
-      
+          EMFSTK.espark[npnw][ispr] = TRACKR.spausr[ispr];
+
       for (ispr = 0; ispr <= mkbmx2-1; ispr++) 
-         EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
-    
+          EMFSTK.iespak[npnw][ispr] = TRACKR.ispusr[ispr];
+
       EMFSTK.louemf[npnw] = TRACKR.llouse;
   }
 
@@ -77,35 +77,36 @@ void stupre()
   
   for (kp = EMFSTK.npstrt - 1; kp <= EMFSTK.npemf - 1; kp++) {
     
-// Ckeck transport cut first
+// Check transport cut first
     Int_t    ireg   = EMFSTK.iremf[kp];
     Double_t cut    = (TMath::Abs(EMFSTK.ichemf[kp]) == 1) ? EMFRGN.elethr[ireg-1] :  EMFRGN.phothr[ireg-1];
     Double_t e      = EMFSTK.etemf[kp];
 
     if ((e < cut) 
-       && ( 
-           (EMFSTK.ichemf[kp] ==  0) ||
-           (EMFSTK.ichemf[kp] == -1) ||
-           (EMFSTK.ichemf[kp] ==  1 &&  EMFRGN.phothr[ireg-1] > emassmev)
-           )
-       )
+        && (
+            (EMFSTK.ichemf[kp] ==  0) ||
+            (EMFSTK.ichemf[kp] == -1) ||
+            (EMFSTK.ichemf[kp] ==  1 &&  EMFRGN.phothr[ireg-1] > emassmev)
+            )
+        )
     {
-       EMFSTK.iespak[kp][mkbmx2-1] = TRACKR.ispusr[mkbmx2-1];
-       EMFSTK.iespak[kp][mkbmx2-2] =  0;
-       EMFSTK.iespak[kp][mkbmx2-3] = TRACKR.jtrack;
-       continue;
+        EMFSTK.iespak[kp][mkbmx2-1] = TRACKR.ispusr[mkbmx2-1];
+        EMFSTK.iespak[kp][mkbmx2-2] =  0;
+        EMFSTK.iespak[kp][mkbmx2-3] = TRACKR.jtrack;
+        EMFSTK.irlatt[kp] = TRACKR.lt1trk;
+        continue;
     }
 
 // Save the parent track number and reset it at each loop
     Int_t done = 0;
     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
     Int_t flukaid = 0;
-    
+
 // Identify particle type
     if      (EMFSTK.ichemf[kp] == -1)  flukaid = kFLUKAelectron;
     else if (EMFSTK.ichemf[kp] ==  0)  flukaid = kFLUKAphoton;
     else if (EMFSTK.ichemf[kp] ==  1)  flukaid = kFLUKApositron;
-    
+
 
     e *= emvgev;
     Int_t    pdg    = fluka->PDGFromId(flukaid);
@@ -130,73 +131,73 @@ void stupre()
 //* all secondaries are true
     if ((EVTFLG.lpairp == 1) || (EVTFLG.lphoel == 1) ||
         (EVTFLG.lannfl == 1) || (EVTFLG.lannrs == 1)) {
-       
-       if (EVTFLG.lpairp == 1) mech = kPPair;
-       else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
-       else mech = kPAnnihilation;
+
+        if (EVTFLG.lpairp == 1) mech = kPPair;
+        else if (EVTFLG.lphoel == 1) mech = kPPhotoelectric;
+        else mech = kPAnnihilation;
         cppstack->PushTrack(done, parent, pdg,
-                          px, py, pz, e, vx, vy, vz, tof,
-                          polx, poly, polz, mech, ntr, weight, is);
-       if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << "energy " << e-PAPROP.am[flukaid+6] << endl;
+                           px, py, pz, e, vx, vy, vz, tof,
+                           polx, poly, polz, mech, ntr, weight, is);
+        if (debug) cout << endl << " !!! stupre (PAIR, ..) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << " energy " << e-PAPROP.am[flukaid+6] << endl;
 
-       EMFSTK.iespak[kp][mkbmx2-1] = ntr;
-       EMFSTK.iespak[kp][mkbmx2-2] = 0;
+        EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+        EMFSTK.iespak[kp][mkbmx2-2] = 0;
     } // end of lpairp, lphoel, lannfl, lannrs
-    
+
 //* Compton: secondary is true only if charged (e+, e-)
     else if ((EVTFLG.lcmptn == 1)) {
 
-       if (EMFSTK.ichemf[kp] != 0) {
-           mech = kPCompton;
-           cppstack->PushTrack(done, parent, pdg,
-                              px, py, pz, e, vx, vy, vz, tof,
-                              polx, poly, polz, mech, ntr, weight, is);
-           if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
-           EMFSTK.iespak[kp][mkbmx2-1] = ntr;
-           EMFSTK.iespak[kp][mkbmx2-2] = 0;
-       }
+        if (EMFSTK.ichemf[kp] != 0) {
+            mech = kPCompton;
+            cppstack->PushTrack(done, parent, pdg,
+                               px, py, pz, e, vx, vy, vz, tof,
+                               polx, poly, polz, mech, ntr, weight, is);
+            if (debug) cout << endl << " !!! stupre (COMPTON) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
+            EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+            EMFSTK.iespak[kp][mkbmx2-2] = 0;
+        }
     } // end of lcmptn
-    
+
 //* Bremsstrahlung: true secondary only if charge = 0 (photon)
     else if ((EVTFLG.lbrmsp == 1)) {
-       if (EMFSTK.ichemf[kp] == 0) {
-           mech = kPBrem;
-           cppstack->PushTrack(done, parent, pdg,
-                              px, py, pz, e, vx, vy, vz, tof,
-                              polx, poly, polz, mech, ntr, weight, is);
-           if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
-           EMFSTK.iespak[kp][mkbmx2-1] = ntr;
-           EMFSTK.iespak[kp][mkbmx2-2] = 0;
-       }
+        if (EMFSTK.ichemf[kp] == 0) {
+            mech = kPBrem;
+            cppstack->PushTrack(done, parent, pdg,
+                               px, py, pz, e, vx, vy, vz, tof,
+                               polx, poly, polz, mech, ntr, weight, is);
+            if (debug) cout << endl << " !!! stupre (BREMS) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
+            EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+            EMFSTK.iespak[kp][mkbmx2-2] = 0;
+        }
     } // end of lbrmsp
-    
+
 //* Delta ray: If Bhabha, true secondary only if negative (electron)
     else if ((EVTFLG.ldltry == 1)) {
-       if (lbhabh == 1) {
-           if (EMFSTK.ichemf[kp] == -1) {
-               mech = kPDeltaRay;
-               cppstack->PushTrack(done, parent, pdg,
-                                  px, py, pz, e, vx, vy, vz, tof,
-                                  polx, poly, polz, mech, ntr, weight, is);
-               EMFSTK.iespak[kp][mkbmx2-1] = ntr;
-               EMFSTK.iespak[kp][mkbmx2-2] = 0;
-          if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
-           } // end of Bhabha
-       } // lbhabh == 1
-       
+        if (lbhabh == 1) {
+            if (EMFSTK.ichemf[kp] == -1) {
+                mech = kPDeltaRay;
+                cppstack->PushTrack(done, parent, pdg,
+                                   px, py, pz, e, vx, vy, vz, tof,
+                                   polx, poly, polz, mech, ntr, weight, is);
+                EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+                EMFSTK.iespak[kp][mkbmx2-2] = 0;
+           if (debug) cout << endl << " !!! stupre (BHABA) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
+            } // end of Bhabha
+        } // lbhabh == 1
+
 //* Delta ray: Otherwise Moller: true secondary is the electron with
 //*            lower energy, which has been put higher in the stack
-       else if (kp == EMFSTK.npemf-1) {
-           mech = kPDeltaRay;
-           cppstack->PushTrack(done, parent, pdg,
-                              px, py, pz, e, vx, vy, vz, tof,
-                              polx, poly, polz, mech, ntr, weight, is);
-           if (debug) cout << endl << " !!! stupre (Moller) : ntr=" << ntr << "pdg " << pdg << " parent=" << parent << endl;
-           EMFSTK.iespak[kp][mkbmx2-1] = ntr;
-           EMFSTK.iespak[kp][mkbmx2-2] = 0;
-       } // end of Delta ray
+        else if (kp == EMFSTK.npemf-1) {
+            mech = kPDeltaRay;
+            cppstack->PushTrack(done, parent, pdg,
+                               px, py, pz, e, vx, vy, vz, tof,
+                               polx, poly, polz, mech, ntr, weight, is);
+            if (debug) cout << endl << " !!! stupre (Moller) : ntr=" << ntr << " pdg " << pdg << " parent=" << parent << endl;
+            EMFSTK.iespak[kp][mkbmx2-1] = ntr;
+            EMFSTK.iespak[kp][mkbmx2-2] = 0;
+        } // end of Delta ray
     } // end of ldltry
-    
+
   } // end of loop
 } // end of stupre
 } // end of extern "C"
index e80f2d1..67cd585 100644 (file)
@@ -82,7 +82,7 @@ extern "C" {
 
     Int_t parent =  TRACKR.ispusr[mkbmx2-1];
     Int_t kpart  = GENSTK.kpart[numsec-1];
-    if (kpart < -6) return;
+    if (kpart < -6) return; // -7 to -12 = "heavy" fragment
 
     Int_t pdg = fluka->PDGFromId(kpart);
 
@@ -106,6 +106,7 @@ extern "C" {
     if (EVTFLG.ldecay == 1) {
         mech = kPDecay;
         if (debug) cout << endl << "Decay" << endl;
+        FLKSTK.nlattc[FLKSTK.npflka] = TRACKR.lt1trk;
     } else if (EVTFLG.ldltry == 1) {
         mech = kPDeltaRay;
         if( fluka->GetIcode() == kKASHEA ) {
@@ -144,7 +145,7 @@ extern "C" {
     if (debug)
        cout << endl << " !!! stuprf: ntr=" << ntr << " pdg " << pdg << " parent=" << parent
              << " parent_pdg="<< fluka->PDGFromId(TRACKR.jtrack) << " numsec "
-             << numsec << " npprmr " << npprmr << " icode=" << fluka->GetIcode()
+             << numsec << " npprmr " << npprmr << " icode=" << fluka->GetIcode() << endl
              << endl;
 
 //
index 3074d48..78ce06a 100644 (file)
@@ -12,6 +12,9 @@
 #else
 # define usdraw USDRAW
 #endif
+
+#include "TGeoManager.h" // <- delete
+
 extern "C" {
 void usdraw(Int_t& icode, Int_t& mreg, 
             Double_t& xsco, Double_t& ysco, Double_t& zsco)
@@ -26,32 +29,50 @@ void usdraw(Int_t& icode, Int_t& mreg,
 
   if (icode/100 == kEMFSCO) {
       for (Int_t npnw = EMFSTK.npstrt-1; npnw <= EMFSTK.npemf-1; npnw++) {
-         if (EMFSTK.iespak[npnw][mkbmx2-1] ==  TRACKR.ispusr[mkbmx2 - 1] ) {
-             EMFSTK.iespak[npnw][mkbmx2 - 2] = 1;
+          if (EMFSTK.iespak[npnw][mkbmx2-1] ==  TRACKR.ispusr[mkbmx2 - 1] ) {
+              EMFSTK.iespak[npnw][mkbmx2 - 2] = 1;
 // Save properties at point where particle disappears in case this is only an interruption
-             TLorentzVector p;
-             gMC->TrackMomentum(p);
-             EMFSTK.espark[npnw][0] = xsco;               // x
-             EMFSTK.espark[npnw][1] = ysco;               // y
-             EMFSTK.espark[npnw][2] = zsco;               // z
-             EMFSTK.espark[npnw][3] = gMC->TrackTime();   // t
-             EMFSTK.espark[npnw][4] = p[0];               // px
-             EMFSTK.espark[npnw][5] = p[1];               // py
-             EMFSTK.espark[npnw][6] = p[2];               // pz
-             EMFSTK.espark[npnw][7] = p[3];               // e
-             EMFSTK.espark[npnw][8] = gMC->TrackLength(); // Length 
-         } // Track found in stack
+              TLorentzVector p;
+              gMC->TrackMomentum(p);
+              EMFSTK.espark[npnw][0] = xsco;               // x
+              EMFSTK.espark[npnw][1] = ysco;               // y
+              EMFSTK.espark[npnw][2] = zsco;               // z
+              EMFSTK.espark[npnw][3] = gMC->TrackTime();   // t
+              EMFSTK.espark[npnw][4] = p[0];               // px
+              EMFSTK.espark[npnw][5] = p[1];               // py
+              EMFSTK.espark[npnw][6] = p[2];               // pz
+              EMFSTK.espark[npnw][7] = p[3];               // e
+              EMFSTK.espark[npnw][8] = gMC->TrackLength(); // Length
+          } // Track found in stack
       } // Loop over emf stack 
   } // Electromagnetic process
 
-  Int_t mlttc = LTCLCM.mlatm1;
+  Int_t mlttc = TRACKR.lt1trk; //LTCLCM.mlatm1;
   fluka->SetMreg(mreg, mlttc);
   fluka->SetXsco(xsco);
   fluka->SetYsco(ysco);
   fluka->SetZsco(zsco);
 
-  if (debug) printf("USDRAW: Number of track segments:%6d %6d icode=%d tof=%10.3e track=%d pdg=%d\n",
-  TRACKR.ntrack, TRACKR.mtrack, icode, TRACKR.atrack, TRACKR.ispusr[mkbmx2-1], fluka->PDGFromId(TRACKR.jtrack) );
+    // check region lattice consistency (debug Ernesto)
+    // *****************************************************
+    Int_t nodeId;
+    Int_t volId = fluka->CurrentVolID(nodeId);
+    Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
+    if( mreg != volId  && !gGeoManager->IsOutside() ) {
+       cout << "  usdraw:   track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
+            << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
+            << "               fluka   mreg=" << mreg << " mlttc=" << mlttc << endl
+            << "               TGeo   volId=" << volId << " crtlttc=" << crtlttc << endl
+            << "     common TRACKR   lt1trk=" << TRACKR.lt1trk << " lt2trk=" << TRACKR.lt2trk << endl
+            << "     common LTCLCM   newlat=" << LTCLCM.newlat << " mlatld=" <<  LTCLCM.mlatld << endl
+            << "                     mlatm1=" << LTCLCM.mlatm1 << " mltsen=" <<  LTCLCM.mltsen << endl
+            << "                     mltsm1=" << LTCLCM.mltsm1 << " mlattc=" << LTCLCM.mlattc << endl;
+        if( mlttc == crtlttc ) cout << "   *************************************************************" << endl;
+    }
+    // *****************************************************
+
+  if (debug) printf("USDRAW: Number of track segments:%6d %6d icode=%d tof=%10.3e track=%d pdg=%d mreg=%d\n",
+  TRACKR.ntrack, TRACKR.mtrack, icode, TRACKR.atrack, TRACKR.ispusr[mkbmx2-1], fluka->PDGFromId(TRACKR.jtrack), mreg );
 
   TVirtualMCStack* cppstack = fluka->GetStack();
   cppstack->SetCurrentTrack( TRACKR.ispusr[mkbmx2-1] );