G4TheRayTracer.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 //
00030 //
00031 
00032 
00033 #include "G4TheRayTracer.hh"
00034 #include "G4SystemOfUnits.hh"
00035 #include "G4EventManager.hh"
00036 #include "G4RTMessenger.hh"
00037 #include "G4RayShooter.hh"
00038 #include "G4VFigureFileMaker.hh"
00039 #include "G4RTTrackingAction.hh"
00040 #include "G4RTSteppingAction.hh"
00041 #include "G4RayTrajectory.hh"
00042 #include "G4RayTrajectoryPoint.hh"
00043 #include "G4RTJpegMaker.hh"
00044 #include "G4RTSimpleScanner.hh"
00045 #include "G4GeometryManager.hh"
00046 #include "G4SDManager.hh"
00047 #include "G4StateManager.hh"
00048 #include "G4Event.hh"
00049 #include "G4TrajectoryContainer.hh"
00050 #include "G4Colour.hh"
00051 #include "G4VisAttributes.hh"
00052 #include "G4UImanager.hh"
00053 #include "G4TransportationManager.hh"
00054 #include "G4RegionStore.hh"
00055 #include "G4ProductionCutsTable.hh"
00056 
00057 G4TheRayTracer::G4TheRayTracer(G4VFigureFileMaker* figMaker,
00058                                G4VRTScanner* scanner)
00059 {
00060   theFigMaker = figMaker;
00061   if(!theFigMaker) theFigMaker = new G4RTJpegMaker;
00062   theScanner = scanner;
00063   if(!theScanner) theScanner = new G4RTSimpleScanner;
00064   theRayShooter = new G4RayShooter();
00065   theRayTracerEventAction = 0;
00066   theRayTracerStackingAction = 0;
00067   theRayTracerTrackingAction = 0;
00068   theRayTracerSteppingAction = 0;
00069   theMessenger = G4RTMessenger::GetInstance(this,theRayTracerSteppingAction);
00070   theEventManager = G4EventManager::GetEventManager();
00071 
00072   nColumn = 640;
00073   nRow = 640;
00074 
00075   eyePosition = G4ThreeVector(1.*m,1.*m,1.*m);
00076   targetPosition = G4ThreeVector(0.,0.,0.);
00077   lightDirection = G4ThreeVector(-0.1,-0.2,-0.3).unit();
00078   up = G4ThreeVector(0,1,0);
00079   viewSpan = 5.0*deg;
00080   headAngle = 0.;
00081   attenuationLength = 1.0*m;
00082 
00083   distortionOn = false;
00084   antialiasingOn = false;
00085 
00086   backgroundColour = G4Colour(1.,1.,1.);
00087 }
00088 
00089 G4TheRayTracer::~G4TheRayTracer()
00090 {
00091   delete theRayShooter;
00092   if(theRayTracerTrackingAction) delete theRayTracerTrackingAction;
00093   if(theRayTracerSteppingAction) delete theRayTracerSteppingAction;
00094   delete theMessenger;
00095   delete theScanner;
00096   delete theFigMaker;
00097 }
00098 
00099 void G4TheRayTracer::Trace(G4String fileName)
00100 {
00101   G4StateManager* theStateMan = G4StateManager::GetStateManager();
00102   G4ApplicationState currentState = theStateMan->GetCurrentState();
00103   if(currentState!=G4State_Idle)
00104   {
00105     G4cerr << "Illegal application state - Trace() ignored." << G4endl;
00106     return;
00107   }
00108 
00109   if(!theFigMaker)
00110   {
00111     G4cerr << "Figure file maker class is not specified - Trace() ignored." << G4endl;
00112     return;
00113   }
00114 
00115   G4UImanager* UI = G4UImanager::GetUIpointer();
00116   G4int storeTrajectory = UI->GetCurrentIntValue("/tracking/storeTrajectory");
00117   if(storeTrajectory==0) UI->ApplyCommand("/tracking/storeTrajectory 1");
00118 
00119 
00120   G4ThreeVector tmpVec = targetPosition - eyePosition;
00121   eyeDirection = tmpVec.unit();
00122   colorR = new unsigned char[nColumn*nRow];
00123   colorG = new unsigned char[nColumn*nRow];
00124   colorB = new unsigned char[nColumn*nRow];
00125 
00126   StoreUserActions();
00127   G4bool succeeded = CreateBitMap();
00128   if(succeeded)
00129   { CreateFigureFile(fileName); }
00130   else
00131   { G4cerr << "Could not create figure file" << G4endl;
00132     G4cerr << "You might set the eye position outside of the world volume" << G4endl; }
00133   RestoreUserActions();
00134 
00135   if(storeTrajectory==0) UI->ApplyCommand("/tracking/storeTrajectory 0");
00136 
00137   delete [] colorR;
00138   delete [] colorG;
00139   delete [] colorB;
00140 }
00141 
00142 void G4TheRayTracer::StoreUserActions()
00143 { 
00144   theUserEventAction = theEventManager->GetUserEventAction();
00145   theUserStackingAction = theEventManager->GetUserStackingAction();
00146   theUserTrackingAction = theEventManager->GetUserTrackingAction();
00147   theUserSteppingAction = theEventManager->GetUserSteppingAction();
00148 
00149   if(!theRayTracerTrackingAction) theRayTracerTrackingAction = new G4RTTrackingAction();
00150   if(!theRayTracerSteppingAction) theRayTracerSteppingAction = new G4RTSteppingAction();
00151 
00152   theEventManager->SetUserAction(theRayTracerEventAction);
00153   theEventManager->SetUserAction(theRayTracerStackingAction);
00154   theEventManager->SetUserAction(theRayTracerTrackingAction);
00155   theEventManager->SetUserAction(theRayTracerSteppingAction);
00156 
00157   G4SDManager* theSDMan = G4SDManager::GetSDMpointerIfExist();
00158   if(theSDMan)
00159   { theSDMan->Activate("/",false); }
00160 
00161   G4GeometryManager* theGeomMan = G4GeometryManager::GetInstance();
00162   theGeomMan->OpenGeometry();
00163   theGeomMan->CloseGeometry(true);
00164 }
00165 
00166 void G4TheRayTracer::RestoreUserActions()
00167 {
00168   theEventManager->SetUserAction(theUserEventAction);
00169   theEventManager->SetUserAction(theUserStackingAction);
00170   theEventManager->SetUserAction(theUserTrackingAction);
00171   theEventManager->SetUserAction(theUserSteppingAction);
00172 
00173   G4SDManager* theSDMan = G4SDManager::GetSDMpointerIfExist();
00174   if(theSDMan)
00175   { theSDMan->Activate("/",true); }
00176 }
00177 
00178 #include "G4ProcessManager.hh"
00179 #include "G4ProcessVector.hh"
00180 #include "G4Geantino.hh"
00181 
00182 G4bool G4TheRayTracer::CreateBitMap()
00183 {
00184   G4int iEvent = 0;
00185   G4double stepAngle = viewSpan/100.;
00186   G4double viewSpanX = stepAngle*nColumn;
00187   G4double viewSpanY = stepAngle*nRow;
00188   G4bool succeeded;
00189 
00190 // Confirm process(es) of Geantino is initialized
00191   G4VPhysicalVolume* pWorld =
00192         G4TransportationManager::GetTransportationManager()->
00193         GetNavigatorForTracking()->GetWorldVolume();
00194   G4RegionStore::GetInstance()->UpdateMaterialList(pWorld);
00195   G4ProductionCutsTable::GetProductionCutsTable()->UpdateCoupleTable(pWorld);
00196   G4ProcessVector* pVector
00197     = G4Geantino::GeantinoDefinition()->GetProcessManager()->GetProcessList();
00198   for (G4int j=0; j < pVector->size(); ++j) {
00199       (*pVector)[j]->BuildPhysicsTable(*(G4Geantino::GeantinoDefinition()));
00200   }
00201 
00202 // Close geometry and set the application state
00203   G4GeometryManager* geomManager = G4GeometryManager::GetInstance();
00204   geomManager->OpenGeometry();
00205   geomManager->CloseGeometry(1,0);
00206   
00207   G4ThreeVector center(0,0,0);
00208   G4Navigator* navigator =
00209       G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking();
00210   navigator->LocateGlobalPointAndSetup(center,0,false);
00211 
00212   G4StateManager* theStateMan = G4StateManager::GetStateManager();
00213   theStateMan->SetNewState(G4State_GeomClosed); 
00214 
00215 // Event loop
00216   theScanner->Initialize(nRow,nColumn);
00217   G4int iRow, iColumn;
00218   while (theScanner->Coords(iRow,iColumn)) {
00219       G4int iCoord = iRow * nColumn + iColumn;
00220       G4double dRow = 0, dColumn = 0;  // Antialiasing increments.
00221       G4Event* anEvent = new G4Event(iEvent++);
00222       G4double angleX = -(viewSpanX/2. - (iColumn+dColumn)*stepAngle);
00223       G4double angleY = viewSpanY/2. - (iRow+dRow)*stepAngle;
00224       G4ThreeVector rayDirection;
00225       if(distortionOn)
00226       {
00227         rayDirection = G4ThreeVector(-std::tan(angleX)/std::cos(angleY),std::tan(angleY)/std::cos(angleX),1.0);
00228       }
00229       else
00230       {
00231         rayDirection = G4ThreeVector(-std::tan(angleX),std::tan(angleY),1.0);
00232       }
00233       G4double cp = std::cos(eyeDirection.phi());
00234       G4double sp = std::sqrt(1.-cp*cp);
00235       G4double ct = std::cos(eyeDirection.theta());
00236       G4double st = std::sqrt(1.-ct*ct);
00237       G4double gamma = std::atan2(ct*cp*up.x()+ct*sp*up.y()-st*up.z(), -sp*up.x()+cp*up.y());
00238       rayDirection.rotateZ(-gamma);
00239       rayDirection.rotateZ(headAngle);
00240       rayDirection.rotateUz(eyeDirection);
00241       G4ThreeVector rayPosition(eyePosition);
00242       G4bool interceptable = true;
00243       // Check if rayPosition is in the world.
00244       EInside whereisit =
00245         pWorld->GetLogicalVolume()->GetSolid()->Inside(rayPosition);
00246       if (whereisit != kInside) {
00247         // It's outside the world, so move it inside.
00248         G4double outsideDistance =
00249           pWorld->GetLogicalVolume()->GetSolid()->
00250           DistanceToIn(rayPosition,rayDirection);  
00251         if (outsideDistance != kInfinity) {
00252           // Borrowing from geometry, where 1e-8 < epsilon < 1e-3, in
00253           // absolute/internal length units, is used for ensuring good
00254           // behaviour, choose to add 0.001 to ensure rayPosition is
00255           // definitely inside the world volume (JA 16/9/2005)...
00256           rayPosition = rayPosition+(outsideDistance+0.001)*rayDirection;
00257         }
00258         else {
00259           interceptable = false;
00260         }
00261       }
00262       if (interceptable) {
00263         theRayShooter->Shoot(anEvent,rayPosition,rayDirection.unit());
00264         theEventManager->ProcessOneEvent(anEvent);
00265         succeeded = GenerateColour(anEvent);
00266         colorR[iCoord] = (unsigned char)(int(255*rayColour.GetRed()));
00267         colorG[iCoord] = (unsigned char)(int(255*rayColour.GetGreen()));
00268         colorB[iCoord] = (unsigned char)(int(255*rayColour.GetBlue()));
00269       } else {  // Ray does not intercept world at all.
00270         // Store background colour...
00271         colorR[iCoord] = (unsigned char)(int(255*backgroundColour.GetRed()));
00272         colorG[iCoord] = (unsigned char)(int(255*backgroundColour.GetGreen()));
00273         colorB[iCoord] = (unsigned char)(int(255*backgroundColour.GetBlue()));
00274         succeeded = true;
00275       }
00276 
00277       theScanner->Draw(colorR[iCoord],colorG[iCoord],colorB[iCoord]);
00278 
00279       delete anEvent;
00280       if(!succeeded) return false;
00281   }
00282 
00283   theStateMan->SetNewState(G4State_Idle); 
00284   return true;
00285 }
00286 
00287 void G4TheRayTracer::CreateFigureFile(G4String fileName)
00288 {
00289   //G4cout << nColumn << " " << nRow << G4endl;
00290   theFigMaker->CreateFigureFile(fileName,nColumn,nRow,colorR,colorG,colorB);
00291 }
00292 
00293 G4bool G4TheRayTracer::GenerateColour(G4Event* anEvent)
00294 {
00295   G4TrajectoryContainer * trajectoryContainer = anEvent->GetTrajectoryContainer();
00296   
00297   G4RayTrajectory* trajectory = (G4RayTrajectory*)( (*trajectoryContainer)[0] );
00298   if(!trajectory) return false;
00299 
00300   G4int nPoint = trajectory->GetPointEntries();
00301   if(nPoint==0) return false;
00302 
00303   G4Colour initialColour(backgroundColour);
00304   if( trajectory->GetPointC(nPoint-1)->GetPostStepAtt() )
00305   { initialColour = GetSurfaceColour(trajectory->GetPointC(nPoint-1)); }
00306   rayColour = Attenuate(trajectory->GetPointC(nPoint-1),initialColour);
00307 
00308   for(int i=nPoint-2;i>=0;i--)
00309   {
00310     G4Colour surfaceColour = GetSurfaceColour(trajectory->GetPointC(i));
00311     G4double weight = 1.0 - surfaceColour.GetAlpha();
00312     G4Colour mixedColour = GetMixedColour(rayColour,surfaceColour,weight);
00313     rayColour = Attenuate(trajectory->GetPointC(i),mixedColour);
00314   }
00315     
00316   return true;
00317 }
00318 
00319 G4Colour G4TheRayTracer::GetMixedColour(G4Colour surfCol,G4Colour transCol,G4double weight)
00320 {
00321   G4double red   = weight*surfCol.GetRed() + (1.-weight)*transCol.GetRed();
00322   G4double green = weight*surfCol.GetGreen() + (1.-weight)*transCol.GetGreen();
00323   G4double blue  = weight*surfCol.GetBlue() + (1.-weight)*transCol.GetBlue();
00324   G4double alpha = weight*surfCol.GetAlpha() + (1.-weight)*transCol.GetAlpha();
00325   return G4Colour(red,green,blue,alpha);
00326 }
00327 
00328 G4Colour G4TheRayTracer::GetSurfaceColour(G4RayTrajectoryPoint* point)
00329 {
00330   const G4VisAttributes* preAtt = point->GetPreStepAtt();
00331   const G4VisAttributes* postAtt = point->GetPostStepAtt();
00332 
00333   G4bool preVis = ValidColour(preAtt);
00334   G4bool postVis = ValidColour(postAtt);
00335 
00336   G4Colour transparent(1.,1.,1.,0.);
00337 
00338   if(!preVis&&!postVis) return transparent;
00339 
00340   G4ThreeVector normal = point->GetSurfaceNormal();
00341 
00342   G4Colour preCol(1.,1.,1.);
00343   G4Colour postCol(1.,1.,1.);
00344 
00345   if(preVis)
00346   {
00347     G4double brill = (1.0-(-lightDirection).dot(normal))/2.0;
00348     G4double red   = preAtt->GetColour().GetRed();
00349     G4double green = preAtt->GetColour().GetGreen();
00350     G4double blue  = preAtt->GetColour().GetBlue();
00351     preCol = G4Colour
00352       (red*brill,green*brill,blue*brill,preAtt->GetColour().GetAlpha());
00353   }
00354   else
00355   { preCol = transparent; }
00356 
00357   if(postVis)
00358   {
00359     G4double brill = (1.0-(-lightDirection).dot(-normal))/2.0;
00360     G4double red   = postAtt->GetColour().GetRed();
00361     G4double green = postAtt->GetColour().GetGreen();
00362     G4double blue  = postAtt->GetColour().GetBlue();
00363     postCol = G4Colour
00364       (red*brill,green*brill,blue*brill,postAtt->GetColour().GetAlpha());
00365   }
00366   else
00367   { postCol = transparent; }
00368     
00369   if(!preVis) return postCol;
00370   if(!postVis) return preCol;
00371 
00372   G4double weight = 0.5;
00373   return GetMixedColour(preCol,postCol,weight);
00374 }
00375 
00376 G4Colour G4TheRayTracer::Attenuate(G4RayTrajectoryPoint* point, G4Colour sourceCol)
00377 {
00378   const G4VisAttributes* preAtt = point->GetPreStepAtt();
00379 
00380   G4bool visible = ValidColour(preAtt);
00381   if(!visible) return sourceCol;
00382 
00383   G4Colour objCol = preAtt->GetColour();
00384   G4double stepRed = objCol.GetRed();
00385   G4double stepGreen = objCol.GetGreen();
00386   G4double stepBlue = objCol.GetBlue();
00387   G4double stepAlpha = objCol.GetAlpha();
00388   G4double stepLength = point->GetStepLength();
00389 
00390   G4double attenuationFuctor;
00391   if(stepAlpha > 0.9999999){ stepAlpha = 0.9999999; } // patch to the next line
00392     attenuationFuctor = -stepAlpha/(1.0-stepAlpha)*stepLength/attenuationLength;
00393  
00394   G4double KtRed = std::exp((1.0-stepRed)*attenuationFuctor);
00395   G4double KtGreen = std::exp((1.0-stepGreen)*attenuationFuctor);
00396   G4double KtBlue = std::exp((1.0-stepBlue)*attenuationFuctor);
00397   if(KtRed>1.0){KtRed=1.0;}
00398   if(KtGreen>1.0){KtGreen=1.0;}
00399   if(KtBlue>1.0){KtBlue=1.0;}
00400   return G4Colour(sourceCol.GetRed()*KtRed,
00401     sourceCol.GetGreen()*KtGreen,sourceCol.GetBlue()*KtBlue);
00402 }
00403 
00404 G4bool G4TheRayTracer::ValidColour(const G4VisAttributes* visAtt)
00405 {
00406   G4bool val = true;
00407   if(!visAtt)
00408   { val = false; }
00409   else if(!(visAtt->IsVisible()))
00410   { val = false; }
00411   else if(visAtt->IsForceDrawingStyle()
00412     &&(visAtt->GetForcedDrawingStyle()==G4VisAttributes::wireframe))
00413   { val = false; }
00414   return val;
00415 }
00416 

Generated on Mon May 27 17:50:00 2013 for Geant4 by  doxygen 1.4.7