Flugg/WrapG1.cxx

   1
   2 // Flugg tag
   3
   4 ///////////////////////////////////////////////////////////////////
   5 //
   6 // WrapG1.hh - Sara Vanini
   7 //
   8 // Wrapper for geometry tracking: returns approved step of
   9 // particle and alla variables that fluka G1 computes.
  10 //
  11 // modified 11/III/99 : included jrLtGeant array for storing
  12 // lattice histories; fixed jump-on-boundaries and back-scattering
  13 // modified 20/IV/99 : history counters of jrLt array are
  14 // incremented when created and decremented when check is required:
  15 // this for not deleting histories twice.
  16 // modified 18/X/99 : LocateGlobalPointWithinVolume used when position
  17 // is changed from last located point.
  18 // modified 1/III/00 : update utilities histories when crossing
  19 // identical volume boundaries.
  20 //
  21 // modified 22/III/00 : fixed LttcFlag and jrLt return values.
  22 // modified 12/VI/00 : end-step on Boundary bug fixed.
  23 // modified 5/VII/00 : boundary not seen by G4 geometry bug fixed.
  24 // modified 24.10.01: by I. Hrivnacova
  25 //   functions declarations separated from implementation
  26 //   (moved to Wrappers.hh);
  27 //
  28 /////////////////////////////////////////////////////////////////////
  29
  30 #include "Wrappers.hh"
  31 #include "WrapUtils.hh"
  32 #include "FGeometryInit.hh"
  33 #include "NavHistWithCount.hh"
  34 #include "G4VPhysicalVolume.hh"
  35 #include "G4NavigationHistory.hh"
  36 #include "FluggNavigator.hh"
  37 #include "G4PhysicalVolumeStore.hh"
  38 #include "globals.hh"
  39
  40
  41 void g1wr(G4double& pSx, G4double& pSy, G4double& pSz, G4double* pV,
  42           G4int& oldReg, const G4int& oldLttc, G4double& propStep,
  43           G4int& nascFlag, G4double& retStep, G4int& newReg,
  44           G4double& saf, G4int& newLttc, G4int& LttcFlag,
  45           G4double* sLt, G4int* jrLt)
  46 {
  47   //flag
  48 #ifdef G4GEOMETRY_DEBUG
  49   G4cout<<"============= G1WR =============="<<G4endl;
  50 #endif
  51
  52   //static G4int count=0;
  53   //count+=1;
  54   //G4cout<<"contatore G1="<<count<<G4endl;
  55
  56   ///////////////////////// INPUT ///////////////////////////
  57   //Geoinit, Navigator, TouchableHistory pointers
  58   static FGeometryInit * ptrGeoInit = FGeometryInit::GetInstance();
  59   FluggNavigator * ptrNavig = ptrGeoInit->getNavigatorForTracking();
  60   G4TouchableHistory * ptrTouchHist = ptrGeoInit->GetTouchableHistory();
  61   G4PhysicalVolumeStore * pVolStore = G4PhysicalVolumeStore::GetInstance();
  62
  63   //setting variables (and dimension: Fluka uses cm.!)
  64   G4ThreeVector partLoc(pSx,pSy,pSz);
  65   partLoc *= 10.0; // in millimeters!
  66   static G4ThreeVector partLocOld = partLoc;
  67   static G4ThreeVector oldLocalPoint =
  68     ptrNavig->ComputeLocalPoint(partLocOld);
  69
  70   G4ThreeVector pVec(pV[0],pV[1],pV[2]);
  71   const G4double physStep=G4double(propStep*10.);
  72
  73   G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
  74
  75
  76 #ifdef G4GEOMETRY_DEBUG
  77   G4cout.precision(10);
  78   G4cout << "Position (cm):" << pSx << "," << pSy << "," << pSz << G4endl;
  79   G4cout << "Direction: "    << pVec << G4endl;
  80   G4cout << "Proposed step :"<< propStep << G4endl;
  81 #endif
  82
  83
  84
  85   ///////////////////// FLUKA/G4 REGIONS COMPARISON //////////////////////
  86   //get oldReg pointer and G4volume of previous step end-point pointer
  87   G4VPhysicalVolume * ptrOldReg = (*pVolStore)[oldReg-1];
  88   G4VPhysicalVolume * lastVolume = ptrTouchHist->GetVolume();
  89   G4int touVolCpNb = 0;
  90   if (lastVolume)
  91     touVolCpNb = lastVolume->GetCopyNo();
  92
  93 #ifdef  G4GEOMETRY_DEBUG
  94   G4int fluVolCpNb = ptrOldReg->GetCopyNo();
  95   G4cout << "Fluka volume before step: " << ptrOldReg->GetName()
  96          << "," << fluVolCpNb<< G4endl;
  97   if(lastVolume)
  98     G4cout << "G4 Touch Hist volume: "
  99            << lastVolume->GetName() << "," << touVolCpNb << G4endl;
 100   G4cout << "------------------------------------------------" << G4endl;
 101 #endif
 102
 103
 104   //if volume is changed, this is a new particle tracking, or fluka tries
 105   //to reach a boundary more softly with the same particle. In this case
 106   //fluka restart tracking from old history, in general. For tracking in
 107   //lattice-volume fluka goes back to one of the previous lattice volumes.
 108   //Check if ptrOldReg is equal to old history, or to one of the N lattice
 109   //volumes stored in jrLt. Then reinitialise step-histories! Otherwise
 110   //must relocate.
 111   //NB. jrLtGeant stores lattice volume histories until LttcFlag==-1,
 112   //then all histories are checked and deleted, and array is reinitialised
 113
 114
 115   G4int haveHistNb = -1;
 116   G4int newRegErr=0;
 117   G4int indHist = LttcFlagGeant;
 118   G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
 119
 120
 121   G4NavigationHistory* ptrLttcHist=0;
 122   if(oldLttc)
 123     ptrLttcHist = reinterpret_cast<NavHistWithCount*>(oldLttc)->GetNavHistPtr();
 124
 125
 126   while(indHist>=0 && haveHistNb==-1) {
 127 #ifdef  G4GEOMETRY_DEBUG
 128     G4cout<<"Searching in jrLtc...."<<G4endl;
 129 #endif
 130     if(oldLttc==jrLtGeant[indHist])
 131       haveHistNb=indHist;
 132     indHist-=1;
 133   }
 134
 135   if(haveHistNb!=-1) {
 136     //fluka history found in jrLtGeant
 137     if(haveHistNb<LttcFlagGeant) {
 138 #ifdef  G4GEOMETRY_DEBUG
 139       G4cout<<"* Fluka reaches boundary more softly..."<<G4endl;
 140       G4cout<<"* Re-initializing FluggNavigator history"<<G4endl;
 141       G4cout<<"* and updating step-histories"<<G4endl;
 142 #endif
 143
 144       ptrNavig->UpdateNavigatorHistory(ptrLttcHist);
 145       ptrGeoInit->UpdateHistories(ptrLttcHist,0);
 146     }
 147 #ifdef  G4GEOMETRY_DEBUG
 148     if(haveHistNb==LttcFlagGeant)
 149       G4cout << "Continuing step...." << G4endl;
 150 #endif
 151     jrLt[0]=oldLttc;
 152   }
 153   else {
 154     //not found fluka history in jrLttGeant!
 155     G4cout << "* ERROR! Geometry not correctly initialised in fluka history!"
 156            << G4endl;
 157
 158     //relocation!
 159     ptrNavig->LocateGlobalPointAndUpdateTouchable(partLoc,0,ptrTouchHist,true);
 160
 161     G4cout << "* ATTENTION: point relocation in: "
 162            << ptrTouchHist->GetVolume()->GetName() << G4endl;
 163
 164     ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),1);
 165
 166     if(ptrTouchHist->GetVolume() != ptrOldReg) {
 167       G4cout << "* ERROR! Point not in fluka volume!" << G4endl;
 168       newRegErr=-3;
 169     }
 170
 171     //save new history in jrLt[0] and increment its counter
 172 #ifdef G4GEOMETRY_DEBUG
 173     G4cout << "* ISVHWR call to store new NavHistWithCount in jrLt[0]"
 174            << G4endl;
 175 #endif
 176     jrLt[0]=isvhwr(0,0);
 177
 178 #ifdef G4GEOMETRY_DEBUG
 179     G4cout << "* CONHWR call to increment counter" << G4endl;
 180 #endif
 181     G4int incrCount2=1;
 182     conhwr(jrLt[0],&incrCount2);
 183   }
 184
 185
 186   //jrLtGeant - history check: decrement counter and delete
 187   //histories, if LttcFlag=-1, then reinitialise array with -1.
 188   if(LttcFlag==-1 && LttcFlagGeant>=0) {
 189 #ifdef G4GEOMETRY_DEBUG
 190     G4cout << "* CONHWR call to check and delete histories in jrLtGeant[]:"
 191            << G4endl;
 192 #endif
 193
 194     for(G4int ind=0; ind<=LttcFlagGeant; ind++) {
 195       G4int incrCount1=-1;
 196       if(jrLtGeant[ind]!=jrLt[0]) conhwr(jrLtGeant[ind],&incrCount1);
 197       jrLtGeant[ind]=-1;
 198     }
 199     LttcFlagGeant=-1;
 200   }
 201
 202
 203   //update jrLt and sLt arrays
 204   G4int end = 99;
 205   if(LttcFlag>=0)
 206     end=LttcFlag;
 207
 208   // Added by A.Solodkov
 209   if (end>=100) {
 210     G4cout << "Problems in WrapG1 routine" << G4endl;
 211     G4cout << "Index LttcFlag=" << end << " is outside array bounds" << G4endl;
 212     G4cout << "Better to stop immediately !" << G4endl;
 213     exit(1);
 214   }
 215   //jrLt re-initialization with -1 (jrLt[0] is already set)
 216   for(G4int vv=1;vv<=end;vv++)
 217     jrLt[vv]=-1;
 218   //sLt re-initialization
 219   for(G4int vs=0;vs<=end;vs++)
 220     sLt[vs]=0;
 221
 222   LttcFlag=0;
 223
 224
 225   /////////////////////  COMPUTE STEP  ////////////////////////
 226   //update Navigator private flags and voxel stack if point is
 227   //changed from last located point (otherwise troubles come
 228   //when fluka changes location or particle because G4 computes
 229   //from last located point).
 230   G4ThreeVector newLocalPoint = ptrNavig->ComputeLocalPoint(partLoc);
 231   G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
 232   if(moveLenSq>=kCarTolerance*kCarTolerance)
 233     ptrNavig->LocateGlobalPointWithinVolume(partLoc);
 234
 235
 236   //compute step and new location
 237   newReg = oldReg;
 238   G4double appStep = 0;
 239   G4double safety = 0;
 240   G4bool onBoundaryRound = false;
 241   G4bool crossBound = false;
 242   G4double physStepTmp = physStep;
 243   G4bool flagError = false;
 244
 245   while( (oldReg==newReg && appStep<physStepTmp) || onBoundaryRound ) {
 246 #ifdef G4GEOMETRY_DEBUG
 247     G4cout << "* Computing step..." << G4endl;
 248 #endif
 249
 250     //update variables
 251     oldReg = newReg;
 252
 253     if(onBoundaryRound) {
 254       physStepTmp=10.e-10;
 255       //compute step and location: returns newReg
 256       newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
 257                              safety,onBoundaryRound,flagError,oldReg);
 258       physStepTmp=0.;
 259       crossBound=true;
 260     }
 261     else {
 262       physStepTmp -= appStep;
 263       //compute step and location: returns newReg
 264       newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
 265                              safety,onBoundaryRound,flagError,oldReg);
 266     }
 267
 268     G4bool EqualHist;
 269     EqualHist = EqualHistories(ptrTouchHist->GetHistory(),
 270                                ptrGeoInit->GetTempNavHist()->GetHistory());
 271
 272     if(!EqualHist && flagError) {
 273       pVec=-pVec;
 274       newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
 275                              safety,onBoundaryRound,flagError,oldReg);
 276       pVec=-pVec;
 277       physStepTmp+=1;
 278       newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
 279                              safety,onBoundaryRound,flagError,oldReg);
 280
 281       EqualHist = EqualHistories(ptrTouchHist->GetHistory(),
 282                                  ptrGeoInit->GetTempNavHist()->GetHistory());
 283     }
 284
 285     //update sLt
 286     G4double pas = G4double(appStep);
 287     sLt[LttcFlag] += pas/cm;
 288     safety = (oldReg!=newReg)?0.:safety;
 289
 290     if(!EqualHist) {
 291       //end-step point is on boundary between volumes;
 292       //==> update step-histories, save new NavHistWithCount
 293       //and save its pointer in jrLt; save step in sLt.
 294
 295       //set onBoundaryRound=false to avoid re-compute step!
 296       onBoundaryRound=false;
 297
 298 #ifdef G4GEOMETRY_DEBUG
 299       G4cout << "* History is changed!" << G4endl;
 300       G4cout << "* updating step-histories, jrLt, LttcFlag" << G4endl;
 301 #endif
 302
 303       ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
 304
 305       LttcFlag += 1;
 306
 307 #ifdef G4GEOMETRY_DEBUG
 308       G4cout << "* ISVHWR call to store new NavHistWithCount in jrLt"
 309              << G4endl;
 310 #endif
 311
 312       jrLt[LttcFlag] = isvhwr(0,0);
 313
 314 #ifdef G4GEOMETRY_DEBUG
 315       G4cout << "* CONHWR call to increment counter" << G4endl;
 316 #endif
 317       G4int incrCount3=1;
 318       conhwr(jrLt[LttcFlag],&incrCount3);
 319
 320       sLt[LttcFlag] = sLt[LttcFlag-1];
 321     }
 322   }
 323
 324   //////////////////////   OUTPUT   //////////////////////////
 325   //If back-scattering occured, and fluka is in the wrong region, return -3.
 326   //(N. B. Boundary between replicans are not seen when physStep=distance
 327   //particle-boundary: in this case step=kInfinity and history is unchanged.
 328   //Following step is =0 and then history changes.)
 329   if(nascFlag<0 && !appStep && physStep && newReg!=oldReg && !crossBound) {
 330     //don't need to compare histories because boundary between
 331     //identical volumes in different replicans are not seen
 332 #ifdef G4GEOMETRY_DEBUG
 333     G4cout << "* Back-scattering!" << G4endl;
 334 #endif
 335
 336     newReg=-3;
 337   }
 338   else if(newRegErr<0)
 339     newReg=newRegErr;
 340
 341
 342   //compute output variables (in cm.!)
 343   //final step
 344   retStep = sLt[LttcFlag];
 345
 346   //safety (Fluka sottracts a bit to safety to be sure
 347   //not to jump on a boundary)
 348   G4double s = G4double(safety);
 349   s -= s*3.0e-09;
 350   saf=s/cm;
 351
 352   //update wrapper utility variables
 353   //copy jrLt in jrLtGeant
 354   G4int start=0;
 355   if(haveHistNb!=-1 && LttcFlagGeant!=-1)
 356     start=1;
 357   for(G4int lt=start;lt<=LttcFlag;lt++)
 358     jrLtGeant[LttcFlagGeant+1+lt-start]=jrLt[lt];
 359   LttcFlagGeant+=(1+LttcFlag-start);
 360   newLttc = jrLt[LttcFlag];
 361   ptrGeoInit->SetLttcFlagGeant(LttcFlagGeant);
 362
 363   partLocOld=partLoc;
 364   oldLocalPoint = ptrNavig->ComputeLocalPoint(partLocOld);
 365
 366   //compute new position
 367   G4ThreeVector oldPos = G4ThreeVector(pSx,pSy,pSz);
 368   G4ThreeVector newPos = oldPos + retStep*pVec;
 369
 370 #ifdef G4GEOMETRY_DEBUG
 371   G4cout << "New position: " << newPos << G4endl;
 372   G4cout << "Output region: " << newReg << G4endl;
 373   G4cout << "G4 safety (cm): " << (safety*0.1) << G4endl;
 374   G4cout << "Fluka safety (cm): " << saf << G4endl;
 375   G4cout << "Approved step: " << retStep << G4endl;
 376   G4cout << "LttcFlag = " << LttcFlag << G4endl;
 377   for(G4int i=0;i<=LttcFlag+1;i++) {
 378     G4cout << "jrLt[" << i <<"]=" << jrLt[i] << G4endl;
 379     G4cout << "sLt[" << i <<"]=" << sLt[i] << G4endl;
 380   }
 381
 382   G4cout << "LttcFlagGeant =" << LttcFlagGeant << G4endl;
 383   for(G4int ib=0;ib<=LttcFlagGeant+1;ib++) {
 384     G4cout << "jrLtGeant[" << ib <<"]=" << jrLtGeant[ib] << G4endl;
 385   }
 386   G4cout << "newLttc=" << newLttc << G4endl;
 387 #endif
 388 }
 389
 390
 391
 392