Improvements in Wrappers. (A. Gheata)
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 10 Feb 2004 09:24:35 +0000 (09:24 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 10 Feb 2004 09:24:35 +0000 (09:24 +0000)
TFluka/TFlukaMCGeometry.cxx

index 1fb30e8..6c8fcb9 100644 (file)
@@ -1180,7 +1180,7 @@ Int_t idnrwr(const Int_t & /*nreg*/, const Int_t & /*mlat*/)
 // card in fluka input), returns 1 if user wants Fluka always to 
 // use DNEAR (in this case, be sure that GEANT4 DNEAR is unique, 
 // coming from all directions!!!)
-   printf("=> Inside IDNRWR\n");
+   printf("========== Dummy IDNRWR\n");
    return 0;
 }
 
@@ -1188,7 +1188,7 @@ 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)
 {
 //   from FLUGG:
@@ -1196,113 +1196,129 @@ void g1wr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
 // particle and all variables that fluka G1 computes.
 
    // Initialize current point/direction
-   printf("=> Inside G1WR\n");
+   printf("========== Inside G1WR\n");
+   printf("   point/dir:(%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pV[0],pV[1],pV[2]);
    gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
    gGeoManager->SetCurrentDirection(pV);
-   
-   // Initialize old path (FLUKA lattice history)
-   if (oldLttc != jrLt[lttcFlag])
-      printf("Woops: old history not matching jrLt[%i]. Checking other histories.\n",lttcFlag);
-   
-   gGeoManager->CdNode(oldLttc);
-   TGeoVolume *oldvol = (gGeoManager->IsOutside())?0:gGeoManager->GetCurrentVolume();
-   Int_t oldregion = (oldvol)?(mcgeom->NofVolumes()+1):oldvol->GetNumber(); // should it be 0?
-   if (oldregion != oldReg) {
-      while (lttcFlag>=0) {
-         gGeoManager->CdNode(jrLt[lttcFlag]-1);
-         oldvol = (gGeoManager->IsOutside())?0:gGeoManager->GetCurrentVolume();
-         oldregion = (oldvol)?(mcgeom->NofVolumes()+1):oldvol->GetNumber();
-         if (oldregion == oldReg) break;
-         // bad history -> clean up jrLt[lttcFlag], sLt[lttcFlag]
-         sLt[lttcFlag] = 0.;
-         jrLt[lttcFlag] = -1;
-         lttcFlag--;
-      }         
-         
-      if (oldregion != oldReg) {   
-         printf("Error: g1wr: history not found\n");
-         printf("   relocating current point (%f, %f, %f)\n", pSx, pSy, pSz);
-         gGeoManager->FindNode();
-         oldvol = (gGeoManager->IsOutside())?0:gGeoManager->GetCurrentVolume();
-         oldregion = (oldvol)?(mcgeom->NofVolumes()+1):oldvol->GetNumber();
-         lttcFlag = 0;
-         jrLt[lttcFlag]=isvhwr(0,0);
-      }   
-   }
+   printf("   oldReg=%i  oldLttc=%i  pstep=%f\n",oldReg, oldLttc, propStep);
+   Int_t curLttc = gGeoManager->GetCurrentNodeId()+1;
+   Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
+   printf("   curReg=%i  curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
+   Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE;
+   if (!regsame) printf("   REGIONS DOES NOT MATCH\n");
+   if (oldLttc != curLttc) {
+      printf("   HISTORIES DOES NOT MATCH\n");
+      gGeoManager->CdNode(oldLttc-1);
+      curLttc = gGeoManager->GetCurrentNodeId()+1;
+      curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
+      printf("   re-initialized point: curReg=%i  curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
+   }         
+   lttcFlag = 0;
    sLt[lttcFlag] = 0.;   
+   jrLt[lttcFlag] = curLttc;
    // now 'oldregion' contains the real region, matching or not the old history
    
    // Compute geometry step/safety within physical step limit
-   newReg = oldregion;
+//   newReg = oldregion;
+   Double_t *point = gGeoManager->GetCurrentPoint();
+   Double_t *dir = gGeoManager->GetCurrentDirection();
    Double_t steptot = 0.;
    Double_t snext = 0.;
    Int_t istep = 0;
    Bool_t done = kFALSE;
+   Int_t i;
    while (!done) {
-      gGeoManager->FindNextBoundary(propStep-steptot);
+      gGeoManager->FindNextBoundary(-propStep);
       snext = gGeoManager->GetStep();
-      if (steptot == 0) saf = gGeoManager->GetSafeDistance();
+      printf("   FindNextBoundary(%g) snext=%g\n", propStep, snext);
+      if (steptot == 0) {
+         saf = gGeoManager->GetSafeDistance();
+         printf("   Safety: %g\n", saf);
+      }   
+      sLt[lttcFlag] = propStep;
+      jrLt[lttcFlag] = gGeoManager->GetCurrentNodeId()+1;     
+      lttcFlag++; //1
+      newReg = curreg;
+      newLttc = oldLttc;
       if (snext<propStep) {
          // There is a boundary on the way.
          // Make a step=snext+1E-6 to force boundary crossing
+         lttcFlag--; // 0
          steptot += snext;
          sLt[lttcFlag] = snext;
          retStep = snext;
-         gGeoManager->Step();
-         // Hopefully we end-up in a new region, else we do few small steps
-         if (!gGeoManager->IsEntering()) {
-//            sameregion = kTRUE;
-            gGeoManager->SetStep(0.);
-            istep = 0;
-         }   
-         while (!gGeoManager->IsEntering() && steptot<propStep) {
-            gGeoManager->Step();
+//         lttcFlag++;
+         // make the step to get into the next region
+         for (i=0;i<3;i++) point[i]+=(snext+1E-6)*dir[i];
+         gGeoManager->FindNode();
+         istep = 0;
+         printf("   boundary: step made %g\n", snext);
+         while (gGeoManager->IsSameLocation() && steptot<propStep) {
+            if (istep>1E3) {
+               printf("Geometry error: could not cross boundary after extra 10 microns\n");
+               return;
+            }   
+            for (i=0;i<3;i++) point[i]+=1E-6*dir[i];
+            gGeoManager->FindNode();
             sLt[lttcFlag] += 1E-6;
             retStep = sLt[lttcFlag];
             steptot += 1E-6;
             istep++;
-            if (istep>1000) {
-            // we already made 10 extra microns and nothing
-               printf("Woops: g1wr: extra 10 microns and no boundary...\n");
-               gGeoManager->SetStep(propStep-steptot-1E-6);
-               gGeoManager->Step();
-               if (gGeoManager->IsEntering()) {
-                  retStep = sLt[lttcFlag];
-                  lttcFlag++;
-                  sLt[lttcFlag] = propStep-steptot;
-                  newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
-               } else {
-                  sLt[lttcFlag] += propStep-steptot;
-               }         
-               return;
-            }   
-         }   
-         if (steptot>propStep) return;
+         }            
+         lttcFlag++; //1                
+         if (steptot>propStep) {printf("Error\n");return;}
          // we managed to cross the boundary -> in which region
          newReg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
-         lttcFlag++;
-         if (gGeoManager->IsOutside()) return;
-         
-      }      
+         sLt[lttcFlag] = propStep-steptot;
+         newLttc = (gGeoManager->IsOutside())?999999999:gGeoManager->GetCurrentNodeId()+1;
+         jrLt[lttcFlag] = newLttc;
+         if (!gGeoManager->IsOutside()) {
+            lttcFlag++; //2
+            printf("   ENTERED region %i, newLttc=%i in: %s\n", newReg,newLttc,gGeoManager->GetPath());
+         } else printf("   EXIT GEOMETRY: BLKHOLE reg=%i\n", newReg);
+      } 
+      // no boundary within proposed step
+      lttcFlag--;
+      done = kTRUE;
+   }   
+   printf("=> newReg=%i newLttc=%i lttcFlag=%i\n", newReg, newLttc, lttcFlag);
+   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]);
+   if (newLttc!=oldLttc) {
+      if (gGeoManager->IsOutside()) {
+         gGeoManager->SetOutside(kFALSE);
+         gGeoManager->CdTop();
+      }   
+      gGeoManager->CdNode(oldLttc-1);
    }   
+   printf("<= G1WR (in: %s)\n", gGeoManager->GetPath());
 }
 
 //_____________________________________________________________________________
 void g1rtwr()
 {
-   printf("=> Inside G1RTWR\n");
+   printf("========== Dummy G1RTWR\n");
 } 
 
 //_____________________________________________________________________________
 void conhwr(Int_t & /*intHist*/, Int_t * /*incrCount*/)
 {
-   printf("=> Inside CONHWR\n");
+   printf("========== Dummy CONHWR\n");
 }
 
 //_____________________________________________________________________________
-void inihwr(Int_t & /*intHist*/)
+void inihwr(Int_t &intHist)
 {
-   printf("=> Inside INIHWR\n");
+   printf("========== Inside INIHWR -> reinitializing history: %i\n", intHist);
+   if (gGeoManager->IsOutside()) gGeoManager->CdTop();
+   if (intHist<=0) {
+//      printf("=== wrong history number\n");
+      return;
+   }
+   if (intHist==0) gGeoManager->CdTop();
+   else gGeoManager->CdNode(intHist-1);
+   printf(" --- current path: %s\n", gGeoManager->GetPath());
+   printf("<= INIHWR\n");
 }
 
 //_____________________________________________________________________________
@@ -1312,8 +1328,9 @@ void  jomiwr(const Int_t & /*nge*/, const Int_t & /*lin*/, const Int_t & /*lou*/
 // Geometry initialization wrapper called by FLUKAM. Provides to FLUKA the
 // number of regions (volumes in TGeo)
    // build application geometry
-   printf("=> Inside JOMIWR\n");
-   flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast()+1;
+   printf("========== Inside JOMIWR\n");
+   flukaReg = gGeoManager->GetListOfUVolumes()->GetEntriesFast();
+   printf("<= JOMIWR: last region=%i\n", flukaReg);
 }   
 
 //_____________________________________________________________________________
@@ -1321,8 +1338,7 @@ void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
             Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
             Int_t &newReg, Int_t &flagErr, Int_t &newLttc)             
 {
-   printf("=> Inside LKDBWR\n");
-   printf("Locate :(%f, %f, %f)\n", pSx, pSy, pSz);
+   printf("========== Inside LKDBWR (%f, %f, %f)\n",pSx, pSy, pSz);
 //   printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
    printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
    TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
@@ -1334,6 +1350,7 @@ void lkdbwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
    newLttc = gGeoManager->GetCurrentNodeId()+1; 
    flagErr = newReg;
    printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
+   printf("<= LKDBWR\n");
 }
 
 //_____________________________________________________________________________
@@ -1341,8 +1358,7 @@ void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
             Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
             Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
 {
-   printf("=> Inside LKFXWR\n");
-   printf("Locate :(%f, %f, %f)\n", pSx, pSy, pSz);
+   printf("========== Inside LKFXWR (%f, %f, %f)\n",pSx, pSy, pSz);
 //   printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
    printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
    TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
@@ -1354,6 +1370,7 @@ void lkfxwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
    newLttc = gGeoManager->GetCurrentNodeId()+1; 
    flagErr = newReg;
    printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
+   printf("<= LKFXWR\n");
 }
 
 //_____________________________________________________________________________
@@ -1361,8 +1378,7 @@ 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)
 {
-   printf("=> Inside LKMGWR\n");
-   printf("Locate :(%f, %f, %f)\n", pSx, pSy, pSz);
+   printf("========== Inside LKMGWR (%f, %f, %f)\n",pSx, pSy, pSz);
 //   printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
    printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
    TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
@@ -1374,6 +1390,7 @@ void lkmgwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
    newLttc = gGeoManager->GetCurrentNodeId()+1; 
    flagErr = newReg;
    printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
+   printf("<= LKMGWR\n");
 }
 
 //_____________________________________________________________________________
@@ -1381,8 +1398,7 @@ void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
           Double_t * /*pV*/, const Int_t &oldReg, const Int_t &oldLttc,
               Int_t &newReg, Int_t &flagErr, Int_t &newLttc)
 {
-   printf("=> Inside LKWR\n");
-   printf("Locate :(%f, %f, %f)\n", pSx, pSy, pSz);
+   printf("========== Inside LKWR (%f, %f, %f)\n",pSx, pSy, pSz);
 //   printf("   in: pV=(%f, %f, %f)\n", pV[0], pV[1], pV[2]);
    printf("   in: oldReg=%i oldLttc=%i\n", oldReg, oldLttc);
    TGeoNode *node = gGeoManager->FindNode(pSx, pSy, pSz);
@@ -1393,23 +1409,55 @@ void lkwr(Double_t &pSx, Double_t &pSy, Double_t &pSz,
    newReg = node->GetVolume()->GetNumber();
    newLttc = gGeoManager->GetCurrentNodeId()+1; 
    flagErr = newReg;
-   printf("  out: newReg=%i newLttc=%i\n", newReg, newLttc);
+   printf("  out: newReg=%i newLttc=%i in %s\n", newReg, newLttc, gGeoManager->GetPath());
+   printf("<= LKWR\n");
 }
 
 //_____________________________________________________________________________
-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*/)
+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)
 {
-   printf("=> Inside NRMLWR\n");
+   printf("========== Inside NRMLWR (%g, %g, %g, %g, %g, %g)\n", pSx,pSy,pSz,pVx,pVy,pVz);
+   printf("   oldReg=%i, newReg=%i\n", oldReg,newReg);
+   Int_t curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
+   Int_t curLttc = gGeoManager->GetCurrentNodeId()+1;
+   printf("   curReg=%i, curLttc=%i in: %s\n", curreg, curLttc, gGeoManager->GetPath());
+   Bool_t regsame = (curreg==oldReg)?kTRUE:kFALSE;
+   gGeoManager->SetCurrentPoint(pSx, pSy, pSz);
+   gGeoManager->SetCurrentDirection(pVx,pVy,pVz);
+   if (!regsame) {
+      printf("   REGIONS DOEN NOT MATCH\n");
+      gGeoManager->FindNode();
+      curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
+      curLttc = gGeoManager->GetCurrentNodeId()+1;
+      printf("   re-initialized point: curReg=%i  curLttc=%i curPath=%s\n", curreg, curLttc, gGeoManager->GetPath());
+   }
+   Double_t *dnorm = gGeoManager->FindNormalFast();
+   flagErr = 0;
+   if (!dnorm) {
+      printf("   ERROR: Cannot compute fast normal\n");
+      flagErr = 1;
+      norml[0] = -pVx;   
+      norml[1] = -pVy;   
+      norml[2] = -pVz; 
+   }
+   norml[0] = -dnorm[0];   
+   norml[1] = -dnorm[1];   
+   norml[2] = -dnorm[2]; 
+   printf("   normal to boundary: (%g, %g, %g)\n", norml[0], norml[1], norml[2]);  
+   curreg = (gGeoManager->IsOutside())?(mcgeom->NofVolumes()+1):gGeoManager->GetCurrentVolume()->GetNumber();
+   curLttc = gGeoManager->GetCurrentNodeId()+1;
+   printf("   final location: curReg=%i, curLttc=%i in %s\n", curreg,curLttc,gGeoManager->GetPath());
+   printf("<= NRMLWR\n");
 }
 
 //_____________________________________________________________________________
 void rgrpwr(const Int_t & /*flukaReg*/, const Int_t & /*ptrLttc*/, Int_t & /*g4Reg*/,
             Int_t * /*indMother*/, Int_t * /*repMother*/, Int_t & /*depthFluka*/)
 {
-   printf("=> Inside RGRPWR\n");
+   printf("=> Dummy RGRPWR\n");
 }
 
 //_____________________________________________________________________________
@@ -1430,7 +1478,8 @@ Int_t isvhwr(const Int_t &check, const Int_t & intHist)
 
    printf("=> Inside ISVHWR\n");
    if (check<0) return intHist;
-   Int_t histInt = gGeoManager->GetCurrentNodeId();
+   Int_t histInt = gGeoManager->GetCurrentNodeId()+1;
+   printf("<= ISVHWR: history is: %i in: %s\n", histInt, gGeoManager->GetPath());
    return histInt;
 }