Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/visualization/RayTracer/src/G4TheRayTracer.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 //
 28 //
 29 //
 30 
 31 
 32 #include "G4TheRayTracer.hh"
 33 #include "G4SystemOfUnits.hh"
 34 #include "G4EventManager.hh"
 35 #include "G4RTMessenger.hh"
 36 #include "G4RayShooter.hh"
 37 #include "G4VFigureFileMaker.hh"
 38 #include "G4RTTrackingAction.hh"
 39 #include "G4RTSteppingAction.hh"
 40 #include "G4RayTrajectory.hh"
 41 #include "G4RayTrajectoryPoint.hh"
 42 #include "G4RTJpegMaker.hh"
 43 #include "G4RTSimpleScanner.hh"
 44 #include "G4GeometryManager.hh"
 45 #include "G4SDManager.hh"
 46 #include "G4StateManager.hh"
 47 #include "G4Event.hh"
 48 #include "G4TrajectoryContainer.hh"
 49 #include "G4Colour.hh"
 50 #include "G4VisAttributes.hh"
 51 #include "G4UImanager.hh"
 52 #include "G4TransportationManager.hh"
 53 #include "G4RegionStore.hh"
 54 #include "G4ProductionCutsTable.hh"
 55 #include "G4VVisManager.hh"
 56 
 57 #define G4warn G4cout
 58 
 59 G4VFigureFileMaker * G4TheRayTracer::theFigMaker = 0;
 60 G4VRTScanner * G4TheRayTracer::theScanner = 0;
 61 
 62 G4TheRayTracer::G4TheRayTracer(G4VFigureFileMaker* figMaker,
 63              G4VRTScanner* scanner)
 64 {
 65   theFigMaker = figMaker;
 66   if(!theFigMaker) theFigMaker = new G4RTJpegMaker;
 67   theScanner = scanner;
 68   if(!theScanner) theScanner = new G4RTSimpleScanner;
 69   theRayShooter = new G4RayShooter();
 70   theUserEventAction = 0;
 71   theUserStackingAction = 0;
 72   theUserTrackingAction = 0;
 73   theUserSteppingAction = 0;
 74   theRayTracerEventAction = 0;
 75   theRayTracerStackingAction = 0;
 76   theRayTracerTrackingAction = 0;
 77   theRayTracerSteppingAction = 0;
 78   colorR = 0;
 79   colorG = 0;
 80   colorB = 0;
 81 
 82   theMessenger = G4RTMessenger::GetInstance(this);
 83   theEventManager = G4EventManager::GetEventManager();
 84 
 85   nColumn = 640;
 86   nRow = 640;
 87 
 88   eyePosition = G4ThreeVector(1.*m,1.*m,1.*m);
 89   targetPosition = G4ThreeVector(0.,0.,0.);
 90   lightDirection = G4ThreeVector(-0.1,-0.2,-0.3).unit();
 91   up = G4ThreeVector(0,1,0);
 92   viewSpan = 5.0*deg;
 93   headAngle = 0.;
 94   attenuationLength = 1.0*m;
 95 
 96   distortionOn = false;
 97   antialiasingOn = false;
 98 
 99   backgroundColour = G4Colour(1.,1.,1.);
100 }
101 
102 G4TheRayTracer::~G4TheRayTracer()
103 {
104   delete theRayShooter;
105   if(theRayTracerTrackingAction) delete theRayTracerTrackingAction;
106   if(theRayTracerSteppingAction) delete theRayTracerSteppingAction;
107   delete theMessenger;
108   delete theScanner;
109   delete theFigMaker;
110 }
111 
112 void G4TheRayTracer::Trace(const G4String& fileName)
113 {
114   G4StateManager* theStateMan = G4StateManager::GetStateManager();
115   G4ApplicationState currentState = theStateMan->GetCurrentState();
116   if(currentState!=G4State_Idle)
117   {
118     G4warn << "Illegal application state - Trace() ignored." << G4endl;
119     return;
120   }
121 
122   if(!theFigMaker)
123   {
124     G4warn << "Figure file maker class is not specified - Trace() ignored." << G4endl;
125     return;
126   }
127 
128   G4UImanager* UI = G4UImanager::GetUIpointer();
129   G4int storeTrajectory = UI->GetCurrentIntValue("/tracking/storeTrajectory");
130   if(storeTrajectory==0) UI->ApplyCommand("/tracking/storeTrajectory 1");
131 
132 
133   G4ThreeVector tmpVec = targetPosition - eyePosition;
134   eyeDirection = tmpVec.unit();
135   colorR = new unsigned char[nColumn*nRow];
136   colorG = new unsigned char[nColumn*nRow];
137   colorB = new unsigned char[nColumn*nRow];
138 
139   StoreUserActions();
140   G4bool succeeded = CreateBitMap();
141   if(succeeded)
142   { CreateFigureFile(fileName); }
143   else
144   { G4warn << "Could not create figure file" << G4endl;
145     G4warn << "You might set the eye position outside of the world volume" << G4endl; }
146   RestoreUserActions();
147 
148   if(storeTrajectory==0) UI->ApplyCommand("/tracking/storeTrajectory 0");
149 
150   delete [] colorR;
151   delete [] colorG;
152   delete [] colorB;
153 }
154 
155 void G4TheRayTracer::StoreUserActions()
156 { 
157   theUserEventAction = theEventManager->GetUserEventAction();
158   theUserStackingAction = theEventManager->GetUserStackingAction();
159   theUserTrackingAction = theEventManager->GetUserTrackingAction();
160   theUserSteppingAction = theEventManager->GetUserSteppingAction();
161 
162   if(!theRayTracerTrackingAction) theRayTracerTrackingAction = new G4RTTrackingAction();
163   if(!theRayTracerSteppingAction) theRayTracerSteppingAction = new G4RTSteppingAction();
164 
165   theEventManager->SetUserAction(theRayTracerEventAction);
166   theEventManager->SetUserAction(theRayTracerStackingAction);
167   theEventManager->SetUserAction(theRayTracerTrackingAction);
168   theEventManager->SetUserAction(theRayTracerSteppingAction);
169 
170   G4SDManager* theSDMan = G4SDManager::GetSDMpointerIfExist();
171   if(theSDMan)
172   { theSDMan->Activate("/",false); }
173 
174   G4GeometryManager* theGeomMan = G4GeometryManager::GetInstance();
175   theGeomMan->OpenGeometry();
176   theGeomMan->CloseGeometry(true);
177 }
178 
179 void G4TheRayTracer::RestoreUserActions()
180 {
181   theEventManager->SetUserAction(theUserEventAction);
182   theEventManager->SetUserAction(theUserStackingAction);
183   theEventManager->SetUserAction(theUserTrackingAction);
184   theEventManager->SetUserAction(theUserSteppingAction);
185 
186   G4SDManager* theSDMan = G4SDManager::GetSDMpointerIfExist();
187   if(theSDMan)
188   { theSDMan->Activate("/",true); }
189 }
190 
191 #include "G4ProcessManager.hh"
192 #include "G4ProcessVector.hh"
193 #include "G4Geantino.hh"
194 
195 G4bool G4TheRayTracer::CreateBitMap()
196 {
197   G4int iEvent = 0;
198   G4double stepAngle = viewSpan/100.;
199   G4double viewSpanX = stepAngle*nColumn;
200   G4double viewSpanY = stepAngle*nRow;
201   G4bool succeeded;
202 
203   G4VVisManager* visMan = G4VVisManager::GetConcreteInstance();
204   visMan->IgnoreStateChanges(true);
205 
206 // Confirm process(es) of Geantino is initialized
207   G4VPhysicalVolume* pWorld =
208   G4TransportationManager::GetTransportationManager()->
209   GetNavigatorForTracking()->GetWorldVolume();
210   G4RegionStore::GetInstance()->UpdateMaterialList(pWorld);
211   G4ProductionCutsTable::GetProductionCutsTable()->UpdateCoupleTable(pWorld);
212   G4ProcessVector* pVector
213     = G4Geantino::GeantinoDefinition()->GetProcessManager()->GetProcessList();
214   for (G4int j=0; j < (G4int)pVector->size(); ++j) {
215       (*pVector)[j]->BuildPhysicsTable(*(G4Geantino::GeantinoDefinition()));
216   }
217 
218 // Close geometry and set the application state
219   G4GeometryManager* geomManager = G4GeometryManager::GetInstance();
220   geomManager->OpenGeometry();
221   geomManager->CloseGeometry(1,0);
222   
223   G4ThreeVector center(0,0,0);
224   G4Navigator* navigator =
225       G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
226   navigator->LocateGlobalPointAndSetup(center,0,false);
227 
228   G4StateManager* theStateMan = G4StateManager::GetStateManager();
229   theStateMan->SetNewState(G4State_GeomClosed); 
230 
231 // Event loop
232   theScanner->Initialize(nRow,nColumn);
233   G4int iRow, iColumn;
234   while (theScanner->Coords(iRow,iColumn)) {
235       G4int iCoord = iRow * nColumn + iColumn;
236       G4double dRow = 0, dColumn = 0;  // Antialiasing increments.
237       G4Event* anEvent = new G4Event(iEvent++);
238       G4double angleX = -(viewSpanX/2. - (iColumn+dColumn)*stepAngle);
239       G4double angleY = viewSpanY/2. - (iRow+dRow)*stepAngle;
240       G4ThreeVector rayDirection;
241       if(distortionOn)
242       {
243   rayDirection = G4ThreeVector(-std::tan(angleX)/std::cos(angleY),std::tan(angleY)/std::cos(angleX),1.0);
244       }
245       else
246       {
247   rayDirection = G4ThreeVector(-std::tan(angleX),std::tan(angleY),1.0);
248       }
249       G4double cp = std::cos(eyeDirection.phi());
250       G4double sp = std::sqrt(1.-cp*cp);
251       G4double ct = std::cos(eyeDirection.theta());
252       G4double st = std::sqrt(1.-ct*ct);
253       G4double gamma = std::atan2(ct*cp*up.x()+ct*sp*up.y()-st*up.z(), -sp*up.x()+cp*up.y());
254       rayDirection.rotateZ(-gamma);
255       rayDirection.rotateZ(headAngle);
256       rayDirection.rotateUz(eyeDirection);
257       G4ThreeVector rayPosition(eyePosition);
258       G4bool interceptable = true;
259       // Check if rayPosition is in the world.
260       EInside whereisit =
261   pWorld->GetLogicalVolume()->GetSolid()->Inside(rayPosition);
262       if (whereisit != kInside) {
263   // It's outside the world, so move it inside.
264   G4double outsideDistance =
265     pWorld->GetLogicalVolume()->GetSolid()->
266     DistanceToIn(rayPosition,rayDirection);  
267   if (outsideDistance != kInfinity) {
268     // Borrowing from geometry, where 1e-8 < epsilon < 1e-3, in
269     // absolute/internal length units, is used for ensuring good
270     // behaviour, choose to add 0.001 to ensure rayPosition is
271     // definitely inside the world volume (JA 16/9/2005)...
272     rayPosition = rayPosition+(outsideDistance+0.001)*rayDirection;
273   }
274   else {
275     interceptable = false;
276   }
277       }
278       if (interceptable) {
279   theRayShooter->Shoot(anEvent,rayPosition,rayDirection.unit());
280   theEventManager->ProcessOneEvent(anEvent);
281   succeeded = GenerateColour(anEvent);
282   colorR[iCoord] = (unsigned char)(G4int(255*rayColour.GetRed()));
283   colorG[iCoord] = (unsigned char)(G4int(255*rayColour.GetGreen()));
284   colorB[iCoord] = (unsigned char)(G4int(255*rayColour.GetBlue()));
285       } else {  // Ray does not intercept world at all.
286   // Store background colour...
287   colorR[iCoord] = (unsigned char)(G4int(255*backgroundColour.GetRed()));
288   colorG[iCoord] = (unsigned char)(G4int(255*backgroundColour.GetGreen()));
289   colorB[iCoord] = (unsigned char)(G4int(255*backgroundColour.GetBlue()));
290   succeeded = true;
291       }
292 
293       theScanner->Draw(colorR[iCoord],colorG[iCoord],colorB[iCoord]);
294 
295       delete anEvent;
296       if(!succeeded) return false;
297   }
298 
299   theStateMan->SetNewState(G4State_Idle); 
300   visMan->IgnoreStateChanges(false);
301   return true;
302 }
303 
304 void G4TheRayTracer::CreateFigureFile(const G4String& fileName)
305 {
306   //G4cout << nColumn << " " << nRow << G4endl;
307   theFigMaker->CreateFigureFile(fileName,nColumn,nRow,colorR,colorG,colorB);
308 }
309 
310 G4bool G4TheRayTracer::GenerateColour(G4Event* anEvent)
311 {
312   G4TrajectoryContainer * trajectoryContainer = anEvent->GetTrajectoryContainer();
313   
314   G4RayTrajectory* trajectory = (G4RayTrajectory*)( (*trajectoryContainer)[0] );
315   if(!trajectory) return false;
316 
317   G4int nPoint = trajectory->GetPointEntries();
318   if(nPoint==0) return false;
319 
320   G4Colour initialColour(backgroundColour);
321   if( trajectory->GetPointC(nPoint-1)->GetPostStepAtt() )
322   { initialColour = GetSurfaceColour(trajectory->GetPointC(nPoint-1)); }
323   rayColour = Attenuate(trajectory->GetPointC(nPoint-1),initialColour);
324 
325   for(G4int i=nPoint-2;i>=0;--i)
326   {
327     G4Colour surfaceColour = GetSurfaceColour(trajectory->GetPointC(i));
328     G4double weight = 1.0 - surfaceColour.GetAlpha();
329     G4Colour mixedColour = GetMixedColour(rayColour,surfaceColour,weight);
330     rayColour = Attenuate(trajectory->GetPointC(i),mixedColour);
331   }
332     
333   return true;
334 }
335 
336 G4Colour G4TheRayTracer::GetMixedColour
337 (const G4Colour& surfCol,const G4Colour& transCol,G4double weight)
338 {
339   G4double red   = weight*surfCol.GetRed() + (1.-weight)*transCol.GetRed();
340   G4double green = weight*surfCol.GetGreen() + (1.-weight)*transCol.GetGreen();
341   G4double blue  = weight*surfCol.GetBlue() + (1.-weight)*transCol.GetBlue();
342   G4double alpha = weight*surfCol.GetAlpha() + (1.-weight)*transCol.GetAlpha();
343   return G4Colour(red,green,blue,alpha);
344 }
345 
346 G4Colour G4TheRayTracer::GetSurfaceColour(G4RayTrajectoryPoint* point)
347 {
348   const G4VisAttributes* preAtt = point->GetPreStepAtt();
349   const G4VisAttributes* postAtt = point->GetPostStepAtt();
350 
351   G4bool preVis = ValidColour(preAtt);
352   G4bool postVis = ValidColour(postAtt);
353 
354   G4Colour transparent(1.,1.,1.,0.);
355 
356   if(!preVis&&!postVis) return transparent;
357 
358   G4ThreeVector normal = point->GetSurfaceNormal();
359 
360   G4Colour preCol(1.,1.,1.);
361   G4Colour postCol(1.,1.,1.);
362 
363   if(preVis)
364   {
365     const G4Colour& preAttColour = preAtt->GetColour();
366     G4double brill = (1.0-(-lightDirection).dot(normal))/2.0;
367     G4double red   = preAttColour.GetRed();
368     G4double green = preAttColour.GetGreen();
369     G4double blue  = preAttColour.GetBlue();
370     preCol = G4Colour
371       (red*brill,green*brill,blue*brill,preAttColour.GetAlpha());
372   }
373   else
374   { preCol = transparent; }
375 
376   if(postVis)
377   {
378     const G4Colour& postAttColour = postAtt->GetColour();
379     G4double brill = (1.0-(-lightDirection).dot(-normal))/2.0;
380     G4double red   = postAttColour.GetRed();
381     G4double green = postAttColour.GetGreen();
382     G4double blue  = postAttColour.GetBlue();
383     postCol = G4Colour
384       (red*brill,green*brill,blue*brill,postAttColour.GetAlpha());
385   }
386   else
387   { postCol = transparent; }
388     
389   if(!preVis) return postCol;
390   if(!postVis) return preCol;
391 
392   G4double weight = 0.5;
393   return GetMixedColour(preCol,postCol,weight);
394 }
395 
396 G4Colour G4TheRayTracer::Attenuate
397 (G4RayTrajectoryPoint* point,const G4Colour& sourceCol)
398 {
399   const G4VisAttributes* preAtt = point->GetPreStepAtt();
400 
401   G4bool visible = ValidColour(preAtt);
402   if(!visible) return sourceCol;
403 
404   G4Colour objCol = preAtt->GetColour();
405   G4double stepRed = objCol.GetRed();
406   G4double stepGreen = objCol.GetGreen();
407   G4double stepBlue = objCol.GetBlue();
408   G4double stepAlpha = objCol.GetAlpha();
409   G4double stepLength = point->GetStepLength();
410 
411   G4double attenuationFuctor;
412   if(stepAlpha > 0.9999999){ stepAlpha = 0.9999999; } // patch to the next line
413     attenuationFuctor = -stepAlpha/(1.0-stepAlpha)*stepLength/attenuationLength;
414  
415   G4double KtRed = std::exp((1.0-stepRed)*attenuationFuctor);
416   G4double KtGreen = std::exp((1.0-stepGreen)*attenuationFuctor);
417   G4double KtBlue = std::exp((1.0-stepBlue)*attenuationFuctor);
418   if(KtRed>1.0){KtRed=1.0;}
419   if(KtGreen>1.0){KtGreen=1.0;}
420   if(KtBlue>1.0){KtBlue=1.0;}
421   return G4Colour(sourceCol.GetRed()*KtRed,
422     sourceCol.GetGreen()*KtGreen,sourceCol.GetBlue()*KtBlue);
423 }
424 
425 G4bool G4TheRayTracer::ValidColour(const G4VisAttributes* visAtt)
426 {
427   G4bool val = true;
428   if(!visAtt)
429   { val = false; }
430   else if(!(visAtt->IsVisible()))
431   { val = false; }
432   else if(visAtt->IsForceDrawingStyle()
433     &&(visAtt->GetForcedDrawingStyle()==G4VisAttributes::wireframe))
434   { val = false; }
435   return val;
436 }
437 
438