Corrections for the TGeo interface. (A. Gheata)

diff --git a/TFluka/TFlukaMCGeometry.cxx b/TFluka/TFlukaMCGeometry.cxx
index d57facbe23c12a9f11ff9c25cc4e67a8a99aea4e..6b1e649a460727da5f83024e01fc00bfb391ed31 100644 (file)
--- a/TFluka/TFlukaMCGeometry.cxx
+++ b/TFluka/TFlukaMCGeometry.cxx
@@ -1374,13 +1374,9 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
 
 {
    // Initialize FLUKa point and direction;
 
 {
    // Initialize FLUKa point and direction;

+   static Int_t ierr = 0;

    gNstep++;
    gNstep++;

-
-//   if( (gNstep > 43912170 && gNstep < 43912196 ) ||
-//       (gNstep > 47424560 && gNstep < 47424581  ) ||
-//       (gNstep > 54388266 && gNstep < 54388319 )
-//       ) gMCGeom->SetDebugMode();
-//   else gMCGeom->SetDebugMode(kFALSE);
+//      gMCGeom->SetDebugMode(kTRUE);

    
    NORLAT.xn[0] = pSx;
    NORLAT.xn[1] = pSy;
    
    NORLAT.xn[0] = pSx;
    NORLAT.xn[1] = pSy;


@@ -1393,11 +1389,15 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
       oldReg = gGeoManager->GetCurrentVolume()->GetNumber();
    }
 
       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;
+   if (gMCGeom->IsDebugging()) {
+      cout << "g1wr gNstep=" << gNstep << " oldReg="<< oldReg <<" olttc="<< olttc
+           << " track=" << TRACKR.ispusr[mkbmx2-1] << endl;
+      cout << " point: (" << pSx << ", " << pSy << ", " << pSz << ")  dir: ("
+           << pV[0] << ", " << pV[1] << ", " << pV[2] << ")" << endl;
+   }        
+           
+
+   Int_t ccreg=0,cclat=0;

    gMCGeom->GetCurrentRegion(ccreg,cclat);
    Bool_t crossed = (ccreg==oldReg && cclat==olttc)?kFALSE:kTRUE;
 
    gMCGeom->GetCurrentRegion(ccreg,cclat);
    Bool_t crossed = (ccreg==oldReg && cclat==olttc)?kFALSE:kTRUE;
 


@@ -1444,21 +1444,68 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
 
    gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
    gGeoManager->SetCurrentDirection(pV);
 
    gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
    gGeoManager->SetCurrentDirection(pV);

-
+   Double_t pt[3], local[3], ldir[3];
+   Int_t pid = TRACKR.jtrack;
+   pt[0] = pSx;
+   pt[1] = pSy;
+   pt[2] = pSz;
+   gGeoManager->MasterToLocal(pt,local);
+   gGeoManager->MasterToLocalVect(pV,ldir);
+/*
+   Bool_t valid = gGeoManager->GetCurrentVolume()->Contains(local);
+   if (!valid) {
+      printf("location not valid in %s pid=%i\n", gGeoManager->GetPath(),pid);
+      printf("local:(%f, %f, %f)  ldir:(%f, %f, %f)\n", local[0],local[1],local[2],ldir[0],ldir[1],ldir[2]);
+//      gGeoManager->FindNode();
+//      printf("   -> actual location: %s\n", gGeoManager->GetPath());
+   }   
+*/
+   Double_t pstep = propStep;
+   Double_t snext = propStep;
+   const Double_t epsil = 0.9999999 * TGeoShape::Tolerance();
+   // This should never happen !!!
+   if (pstep<TGeoShape::Tolerance()) {
+      printf("Proposed step is 0 !!!\n");
+      pstep = 2.*TGeoShape::Tolerance();
+   }   

    if (crossed) {
    if (crossed) {

-      gGeoManager->FindNextBoundaryAndStep(propStep);
-      saf = 0.0;
+      gGeoManager->FindNextBoundaryAndStep(pstep);
+      snext = gGeoManager->GetStep();
+      saf = 0.;
+      if (snext <= TGeoShape::Tolerance()) {
+//         printf("FLUKA crossed boundary but snext=%g\n", snext);
+         ierr++;
+         snext = epsil;
+      } else {
+         snext += TGeoShape::Tolerance();
+         ierr = 0;
+      }     

    } else {
    } else {

-      gGeoManager->FindNextBoundaryAndStep(propStep, kTRUE);
+      gGeoManager->FindNextBoundaryAndStep(pstep, kTRUE);
+      snext = gGeoManager->GetStep();

       saf = gGeoManager->GetSafeDistance();
       saf = gGeoManager->GetSafeDistance();

-      if (saf<0) saf=0.0;
+      if (snext <= TGeoShape::Tolerance()) {
+//         printf("FLUKA put particle on bondary without crossing\n");
+         ierr++;
+         snext = epsil;
+         saf = 0.;
+      } else {
+         snext += TGeoShape::Tolerance();
+         ierr = 0;
+      }      
+      if (saf<0) saf=0.;

       saf -= saf*3.0e-09;
    }
       saf -= saf*3.0e-09;
    }

-
-   Double_t snext = gGeoManager->GetStep();
-
-   if (snext<=0.0) snext = TGeoShape::Tolerance();
-
+//   if (ierr>1) {
+//      printf("%d snext=%g\n", ierr, snext);
+//   }   
+   if (ierr == 10) {
+//      printf("Too many null steps - sending error code -33...\n");
+      newReg = -33;   // Error code
+      ierr = 0;
+      return;
+   }   
+   

    PAREM.dist = snext;
    NORLAT.distn = snext;
    NORLAT.xn[0] += snext*pV[0];
    PAREM.dist = snext;
    NORLAT.distn = snext;
    NORLAT.xn[0] += snext*pV[0];


@@ -1480,7 +1527,6 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
    // We really crossed the boundary, but is it the same region ?
    gMCGeom->SetNextRegion(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) || (oldReg!=newReg && olttc==newLttc) ) && pid!=-1) {
       // Virtual boundary between replicants
       newReg  = gFluka->GetDummyRegion();
    if ( ((newReg==oldReg && newLttc!=olttc) || (oldReg!=newReg && olttc==newLttc) ) && pid!=-1) {
       // Virtual boundary between replicants
       newReg  = gFluka->GetDummyRegion();


@@ -1603,7 +1649,8 @@ void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
    }   
    flagErr = 0;
       printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
    }   
    flagErr = 0;

-   TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
+   Double_t epsil = 1000.*TGeoShape::Tolerance();
+   TGeoNode *node = gGeoManager->FindNode(pSx+epsil*pV[0], pSy+epsil*pV[1], pSz+epsil*pV[2]);

    if (gGeoManager->IsOutside()) {
       newReg = gMCGeom->NofVolumes()+2;
       newLttc = TFlukaMCGeometry::kLttcOutside;
    if (gGeoManager->IsOutside()) {
       newReg = gMCGeom->NofVolumes()+2;
       newLttc = TFlukaMCGeometry::kLttcOutside;


@@ -1637,8 +1684,8 @@ void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
 
    if( oldReg!=newReg && oldLttc==newLttc ) {
       // this should not happen!! ??? Ernesto
 
    if( oldReg!=newReg && oldLttc==newLttc ) {
       // this should not happen!! ??? Ernesto

-       cout << "  lkwr      oldReg!=newReg ("<< oldReg <<"!="<< newReg
-            << ") && oldLttc==newLttc ("<< newLttc <<") !!!!!!!!!!!!!!!!!" << endl;
+//       cout << "  lkwr      oldReg!=newReg ("<< oldReg <<"!="<< newReg
+//            << ") && oldLttc==newLttc ("<< newLttc <<") !!!!!!!!!!!!!!!!!" << endl;

       newReg  = dummy;
       newLttc = TFlukaMCGeometry::kLttcVirtual;
       flagErr = newReg;
       newReg  = dummy;
       newLttc = TFlukaMCGeometry::kLttcVirtual;
       flagErr = newReg;

diff --git a/TFluka/bxdraw.cxx b/TFluka/bxdraw.cxx
index 931b8dda8ed3796d28b539801c63dc99d2bb7b91..f386743325dd4069fb9794db01bb57db8ed15821 100644 (file)
--- a/TFluka/bxdraw.cxx
+++ b/TFluka/bxdraw.cxx
@@ -52,7 +52,7 @@ void bxdraw(Int_t& icode, Int_t& mreg, Int_t& newreg,
     Int_t nodeId;
     Int_t volId = fluka->CurrentVolID(nodeId);
     Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
     Int_t nodeId;
     Int_t volId = fluka->CurrentVolID(nodeId);
     Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;

-    if( mreg != volId  && !gGeoManager->IsOutside()) {
+    if(debug && 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
        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


@@ -83,7 +83,7 @@ void bxdraw(Int_t& icode, Int_t& mreg, Int_t& newreg,
     Int_t nodeId;
     Int_t volId = fluka->CurrentVolID(nodeId);
     Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
     Int_t nodeId;
     Int_t volId = fluka->CurrentVolID(nodeId);
     Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;

-    if( newreg != volId  && !gGeoManager->IsOutside()) {
+    if(debug && 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
        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

diff --git a/TFluka/endraw.cxx b/TFluka/endraw.cxx
index 593c1229b9e2edf96374cdd6d37cd752b4f67994..81213d1a85d51dd7d70a3b83e5b8128f7ffad45e 100644 (file)
--- a/TFluka/endraw.cxx
+++ b/TFluka/endraw.cxx
@@ -62,7 +62,7 @@ void endraw(Int_t& icode, Int_t& mreg, Double_t& rull, Double_t& xsco, Double_t&
     Int_t nodeId;
     Int_t volId = fluka->CurrentVolID(nodeId);
     Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
     Int_t nodeId;
     Int_t volId = fluka->CurrentVolID(nodeId);
     Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;

-    if( mreg != volId  && !gGeoManager->IsOutside() ) {
+    if(debug && 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
        cout << "  endraw:   track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
             << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
             << "               fluka   mreg=" << mreg << " mlttc=" << mlttc << endl

diff --git a/TFluka/mgdraw.cxx b/TFluka/mgdraw.cxx
index d21ba91f1201fbedea580e6a9c8aa4080bece831..33accf8bf6a0155efe12321d9589a891eb092ba9 100644 (file)
--- a/TFluka/mgdraw.cxx
+++ b/TFluka/mgdraw.cxx
@@ -71,7 +71,7 @@ void mgdraw(Int_t& icode, Int_t& mreg)
 
     // check region lattice consistency (debug Ernesto)
     // *****************************************************
 
     // check region lattice consistency (debug Ernesto)
     // *****************************************************

-    if( mreg != volId  && !gGeoManager->IsOutside() ) {
+    if(verbosityLevel>=3 && 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
        cout << "  mgdraw:   track=" << trackId << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
             << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
             << "               fluka   mreg=" << mreg << " mlttc=" << mlttc << endl

diff --git a/TFluka/usdraw.cxx b/TFluka/usdraw.cxx
index 78ce06ae7d2ca52f3d576e5f7af76223a8e1742d..d8d68f06af66d5c37870b2ad87cacf625e39afa7 100644 (file)
--- a/TFluka/usdraw.cxx
+++ b/TFluka/usdraw.cxx
@@ -58,7 +58,7 @@ void usdraw(Int_t& icode, Int_t& mreg,
     Int_t nodeId;
     Int_t volId = fluka->CurrentVolID(nodeId);
     Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;
     Int_t nodeId;
     Int_t volId = fluka->CurrentVolID(nodeId);
     Int_t crtlttc = gGeoManager->GetCurrentNodeId()+1;

-    if( mreg != volId  && !gGeoManager->IsOutside() ) {
+    if(verbosityLevel>=3 && 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
        cout << "  usdraw:   track=" << TRACKR.ispusr[mkbmx2-1] << " pdg=" << fluka->PDGFromId(TRACKR.jtrack)
             << " icode=" << icode << " gNstep=" << fluka->GetNstep() << endl
             << "               fluka   mreg=" << mreg << " mlttc=" << mlttc << endl