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
  68 // change because of geant4 6.0
  69 //static G4ThreeVector oldLocalPoint =
  70 //  ptrNavig->ComputeLocalPoint(partLocOld);
  71  static G4ThreeVector oldLocalPoint =
  72     ptrNavig->GetGlobalToLocalTransform().TransformPoint(partLocOld);
  73
  74   G4ThreeVector pVec(pV[0],pV[1],pV[2]);
  75   const G4double physStep=G4double(propStep*10.);
  76
  77   G4int LttcFlagGeant = ptrGeoInit->GetLttcFlagGeant();
  78
  79
  80 #ifdef G4GEOMETRY_DEBUG
  81   G4cout.precision(10);
  82   G4cout << "Position (cm):" << pSx << "," << pSy << "," << pSz << G4endl;
  83   G4cout << "Direction: "    << pVec << G4endl;
  84   G4cout << "Proposed step :"<< propStep << G4endl;
  85 #endif
  86
  87
  88
  89   ///////////////////// FLUKA/G4 REGIONS COMPARISON //////////////////////
  90   //get oldReg pointer and G4volume of previous step end-point pointer
  91   G4VPhysicalVolume * ptrOldReg = (*pVolStore)[oldReg-1];
  92   G4VPhysicalVolume * lastVolume = ptrTouchHist->GetVolume();
  93   G4int touVolCpNb = 0;
  94   if (lastVolume)
  95     touVolCpNb = lastVolume->GetCopyNo();
  96
  97 #ifdef  G4GEOMETRY_DEBUG
  98   G4int fluVolCpNb = ptrOldReg->GetCopyNo();
  99   G4cout << "Fluka volume before step: " << ptrOldReg->GetName()
 100          << "," << fluVolCpNb<< G4endl;
 101   if(lastVolume)
 102     G4cout << "G4 Touch Hist volume: "
 103            << lastVolume->GetName() << "," << touVolCpNb << G4endl;
 104   G4cout << "------------------------------------------------" << G4endl;
 105 #endif
 106
 107
 108   //if volume is changed, this is a new particle tracking, or fluka tries
 109   //to reach a boundary more softly with the same particle. In this case
 110   //fluka restart tracking from old history, in general. For tracking in
 111   //lattice-volume fluka goes back to one of the previous lattice volumes.
 112   //Check if ptrOldReg is equal to old history, or to one of the N lattice
 113   //volumes stored in jrLt. Then reinitialise step-histories! Otherwise
 114   //must relocate.
 115   //NB. jrLtGeant stores lattice volume histories until LttcFlag==-1,
 116   //then all histories are checked and deleted, and array is reinitialised
 117
 118
 119   G4int haveHistNb = -1;
 120   G4int newRegErr=0;
 121   G4int indHist = LttcFlagGeant;
 122   G4int * jrLtGeant = ptrGeoInit->GetJrLtGeantArray();
 123
 124
 125   G4NavigationHistory* ptrLttcHist=0;
 126   if(oldLttc)
 127     ptrLttcHist = reinterpret_cast<NavHistWithCount*>(oldLttc)->GetNavHistPtr();
 128
 129
 130   while(indHist>=0 && haveHistNb==-1) {
 131 #ifdef  G4GEOMETRY_DEBUG
 132     G4cout<<"Searching in jrLtc...."<<G4endl;
 133 #endif
 134     if(oldLttc==jrLtGeant[indHist])
 135       haveHistNb=indHist;
 136     indHist-=1;
 137   }
 138
 139   if(haveHistNb!=-1) {
 140     //fluka history found in jrLtGeant
 141     if(haveHistNb<LttcFlagGeant) {
 142 #ifdef  G4GEOMETRY_DEBUG
 143       G4cout<<"* Fluka reaches boundary more softly..."<<G4endl;
 144       G4cout<<"* Re-initializing FluggNavigator history"<<G4endl;
 145       G4cout<<"* and updating step-histories"<<G4endl;
 146 #endif
 147
 148       ptrNavig->UpdateNavigatorHistory(ptrLttcHist);
 149       ptrGeoInit->UpdateHistories(ptrLttcHist,0);
 150     }
 151 #ifdef  G4GEOMETRY_DEBUG
 152     if(haveHistNb==LttcFlagGeant)
 153       G4cout << "Continuing step...." << G4endl;
 154 #endif
 155     jrLt[0]=oldLttc;
 156   }
 157   else {
 158     //not found fluka history in jrLttGeant!
 159     G4cout << "* ERROR! Geometry not correctly initialised in fluka history!"
 160            << G4endl;
 161
 162     //relocation!
 163     ptrNavig->LocateGlobalPointAndUpdateTouchable(partLoc,0,ptrTouchHist,true);
 164
 165     G4cout << "* ATTENTION: point relocation in: "
 166            << ptrTouchHist->GetVolume()->GetName() << G4endl;
 167
 168     ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),1);
 169
 170     if(ptrTouchHist->GetVolume() != ptrOldReg) {
 171       G4cout << "* ERROR! Point not in fluka volume!" << G4endl;
 172       newRegErr=-3;
 173     }
 174
 175     //save new history in jrLt[0] and increment its counter
 176 #ifdef G4GEOMETRY_DEBUG
 177     G4cout << "* ISVHWR call to store new NavHistWithCount in jrLt[0]"
 178            << G4endl;
 179 #endif
 180     jrLt[0]=isvhwr(0,0);
 181
 182 #ifdef G4GEOMETRY_DEBUG
 183     G4cout << "* CONHWR call to increment counter" << G4endl;
 184 #endif
 185     G4int incrCount2=1;
 186     conhwr(jrLt[0],&incrCount2);
 187   }
 188
 189
 190   //jrLtGeant - history check: decrement counter and delete
 191   //histories, if LttcFlag=-1, then reinitialise array with -1.
 192   if(LttcFlag==-1 && LttcFlagGeant>=0) {
 193 #ifdef G4GEOMETRY_DEBUG
 194     G4cout << "* CONHWR call to check and delete histories in jrLtGeant[]:"
 195            << G4endl;
 196 #endif
 197
 198     for(G4int ind=0; ind<=LttcFlagGeant; ind++) {
 199       G4int incrCount1=-1;
 200       if(jrLtGeant[ind]!=jrLt[0]) conhwr(jrLtGeant[ind],&incrCount1);
 201       jrLtGeant[ind]=-1;
 202     }
 203     LttcFlagGeant=-1;
 204   }
 205
 206
 207   //update jrLt and sLt arrays
 208   G4int end = 99;
 209   if(LttcFlag>=0)
 210     end=LttcFlag;
 211
 212   // Added by A.Solodkov
 213   if (end>=100) {
 214     G4cout << "Problems in WrapG1 routine" << G4endl;
 215     G4cout << "Index LttcFlag=" << end << " is outside array bounds" << G4endl;
 216     G4cout << "Better to stop immediately !" << G4endl;
 217     exit(1);
 218   }
 219   //jrLt re-initialization with -1 (jrLt[0] is already set)
 220   for(G4int vv=1;vv<=end;vv++)
 221     jrLt[vv]=-1;
 222   //sLt re-initialization
 223   for(G4int vs=0;vs<=end;vs++)
 224     sLt[vs]=0;
 225
 226   LttcFlag=0;
 227
 228
 229   /////////////////////  COMPUTE STEP  ////////////////////////
 230   //update Navigator private flags and voxel stack if point is
 231   //changed from last located point (otherwise troubles come
 232   //when fluka changes location or particle because G4 computes
 233   //from last located point).
 234
 235 // change because of geant4 6.0
 236 //G4ThreeVector newLocalPoint = ptrNavig->ComputeLocalPoint(partLoc);
 237   G4ThreeVector newLocalPoint =
 238     ptrNavig->GetGlobalToLocalTransform().TransformPoint(partLoc);
 239
 240
 241
 242   G4double moveLenSq = (newLocalPoint-oldLocalPoint).mag2();
 243   if(moveLenSq>=kCarTolerance*kCarTolerance)
 244     ptrNavig->LocateGlobalPointWithinVolume(partLoc);
 245
 246
 247   //compute step and new location
 248   newReg = oldReg;
 249   G4double appStep = 0;
 250   G4double safety = 0;
 251   G4bool onBoundaryRound = false;
 252   G4bool crossBound = false;
 253   G4double physStepTmp = physStep;
 254   G4bool flagError = false;
 255
 256   while( (oldReg==newReg && appStep<physStepTmp) || onBoundaryRound ) {
 257 #ifdef G4GEOMETRY_DEBUG
 258     G4cout << "* Computing step..." << G4endl;
 259 #endif
 260
 261     //update variables
 262     oldReg = newReg;
 263
 264     if(onBoundaryRound) {
 265       physStepTmp=10.e-10;
 266       //compute step and location: returns newReg
 267       newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
 268                              safety,onBoundaryRound,flagError,oldReg);
 269       physStepTmp=0.;
 270       crossBound=true;
 271     }
 272     else {
 273       physStepTmp -= appStep;
 274       //compute step and location: returns newReg
 275       newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
 276                              safety,onBoundaryRound,flagError,oldReg);
 277     }
 278
 279     G4bool EqualHist;
 280     EqualHist = EqualHistories(ptrTouchHist->GetHistory(),
 281                                ptrGeoInit->GetTempNavHist()->GetHistory());
 282
 283     if(!EqualHist && flagError) {
 284       pVec=-pVec;
 285       newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
 286                              safety,onBoundaryRound,flagError,oldReg);
 287       pVec=-pVec;
 288       physStepTmp+=1;
 289       newReg=StepAndLocation(partLoc,pVec,physStepTmp,appStep,
 290                              safety,onBoundaryRound,flagError,oldReg);
 291
 292       EqualHist = EqualHistories(ptrTouchHist->GetHistory(),
 293                                  ptrGeoInit->GetTempNavHist()->GetHistory());
 294     }
 295
 296     //update sLt
 297     G4double pas = G4double(appStep);
 298     sLt[LttcFlag] += pas/cm;
 299     safety = (oldReg!=newReg)?0.:safety;
 300
 301     if(!EqualHist) {
 302       //end-step point is on boundary between volumes;
 303       //==> update step-histories, save new NavHistWithCount
 304       //and save its pointer in jrLt; save step in sLt.
 305
 306       //set onBoundaryRound=false to avoid re-compute step!
 307       onBoundaryRound=false;
 308
 309 #ifdef G4GEOMETRY_DEBUG
 310       G4cout << "* History is changed!" << G4endl;
 311       G4cout << "* updating step-histories, jrLt, LttcFlag" << G4endl;
 312 #endif
 313
 314       ptrGeoInit->UpdateHistories(ptrTouchHist->GetHistory(),2);
 315
 316       LttcFlag += 1;
 317
 318 #ifdef G4GEOMETRY_DEBUG
 319       G4cout << "* ISVHWR call to store new NavHistWithCount in jrLt"
 320              << G4endl;
 321 #endif
 322
 323       jrLt[LttcFlag] = isvhwr(0,0);
 324
 325 #ifdef G4GEOMETRY_DEBUG
 326       G4cout << "* CONHWR call to increment counter" << G4endl;
 327 #endif
 328       G4int incrCount3=1;
 329       conhwr(jrLt[LttcFlag],&incrCount3);
 330
 331       sLt[LttcFlag] = sLt[LttcFlag-1];
 332     }
 333   }
 334
 335   //////////////////////   OUTPUT   //////////////////////////
 336   //If back-scattering occured, and fluka is in the wrong region, return -3.
 337   //(N. B. Boundary between replicans are not seen when physStep=distance
 338   //particle-boundary: in this case step=kInfinity and history is unchanged.
 339   //Following step is =0 and then history changes.)
 340   if(nascFlag<0 && !appStep && physStep && newReg!=oldReg && !crossBound) {
 341     //don't need to compare histories because boundary between
 342     //identical volumes in different replicans are not seen
 343 #ifdef G4GEOMETRY_DEBUG
 344     G4cout << "* Back-scattering!" << G4endl;
 345 #endif
 346
 347     newReg=-3;
 348   }
 349   else if(newRegErr<0)
 350     newReg=newRegErr;
 351 //======================= Endre print ===================
 352 printf("\nWrapG1 mreg=%d",newReg);
 353 //======================= Endre print ===================
 354
 355   //compute output variables (in cm.!)
 356   //final step
 357   retStep = sLt[LttcFlag];
 358
 359   //safety (Fluka sottracts a bit to safety to be sure
 360   //not to jump on a boundary)
 361   G4double s = G4double(safety);
 362   s -= s*3.0e-09;
 363   saf=s/cm;
 364
 365   //update wrapper utility variables
 366   //copy jrLt in jrLtGeant
 367   G4int start=0;
 368   if(haveHistNb!=-1 && LttcFlagGeant!=-1)
 369     start=1;
 370   for(G4int lt=start;lt<=LttcFlag;lt++)
 371     jrLtGeant[LttcFlagGeant+1+lt-start]=jrLt[lt];
 372   LttcFlagGeant+=(1+LttcFlag-start);
 373   newLttc = jrLt[LttcFlag];
 374   ptrGeoInit->SetLttcFlagGeant(LttcFlagGeant);
 375
 376   partLocOld=partLoc;
 377
 378 // change because of geant4 6.0
 379 //oldLocalPoint = ptrNavig->ComputeLocalPoint(partLocOld);
 380   oldLocalPoint =
 381     ptrNavig->GetGlobalToLocalTransform().TransformPoint(partLocOld);
 382
 383   //compute new position
 384   G4ThreeVector oldPos = G4ThreeVector(pSx,pSy,pSz);
 385   G4ThreeVector newPos = oldPos + retStep*pVec;
 386
 387 #ifdef G4GEOMETRY_DEBUG
 388   G4cout << "New position: " << newPos << G4endl;
 389   G4cout << "Output region: " << newReg << G4endl;
 390   G4cout << "G4 safety (cm): " << (safety*0.1) << G4endl;
 391   G4cout << "Fluka safety (cm): " << saf << G4endl;
 392   G4cout << "Approved step: " << retStep << G4endl;
 393   G4cout << "LttcFlag = " << LttcFlag << G4endl;
 394   for(G4int i=0;i<=LttcFlag+1;i++) {
 395     G4cout << "jrLt[" << i <<"]=" << jrLt[i] << G4endl;
 396     G4cout << "sLt[" << i <<"]=" << sLt[i] << G4endl;
 397   }
 398
 399   G4cout << "LttcFlagGeant =" << LttcFlagGeant << G4endl;
 400   for(G4int ib=0;ib<=LttcFlagGeant+1;ib++) {
 401     G4cout << "jrLtGeant[" << ib <<"]=" << jrLtGeant[ib] << G4endl;
 402   }
 403   G4cout << "newLttc=" << newLttc << G4endl;
 404 #endif
 405 }
 406
 407
 408
 409