Jpp test-rotations-new
the software that should make you happy
Loading...
Searching...
No Matches
Functions
JMarkovPathDisplayer.cc File Reference
#include <iostream>
#include <sstream>
#include <iomanip>
#include <vector>
#include <cmath>
#include <cstdlib>
#include "TRandom3.h"
#include "TFile.h"
#include "TPolyLine3D.h"
#include "TPolyMarker3D.h"
#include "TAxis3D.h"
#include "TView.h"
#include "TView3D.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TH3F.h"
#include "Jeep/JParser.hh"
#include "JMarkov/JPhotonPath.hh"
#include "JMarkov/JPhotonPathReader.hh"

Go to the source code of this file.

Functions

int main (int argc, char **argv)
 

Function Documentation

◆ main()

int main ( int argc,
char ** argv )

Definition at line 41 of file JMarkovPathDisplayer.cc.

41 {
42 cout << "JMarkovPathDisplayer" << endl
43 << "Written by Martijn Jongen" << endl
44 << endl ;
45 cout << "Type '" << argv[0] << " -h!' to display the command-line options." << endl ;
46 cout << endl ;
47
48 vector<string> ifnames ;
49 string ofname = "out.root" ;
50 int npaths_to_draw ;
51
52 try {
53 JParser<string> zap ; // this argument parser can handle strings
54 zap["f"] = make_field(ifnames,"input file name (binary file containing JPhotonPaths)") ;
55 zap["o"] = make_field(ofname,"output file name") ;
56 zap["n"] = make_field(npaths_to_draw,"max number of paths to draw") = 200 ;
57
58 if (zap.read(argc, argv) != 0) {
59 return 1 ;
60 }
61 }
62 catch(const exception &error) {
63 // ignore exceptions
64 }
65
66 // read the paths from the file
67 vector<JMARKOV::JPhotonPath> paths ;
68 cout << "Reading files." << endl ;
69 int nread_total = 0 ;
70 for( vector<string>::iterator it=ifnames.begin() ; it!=ifnames.end() ; ++it ) {
71 cout << "Reading '" << *it << "'." << endl ;
72 // read the paths from the file
73 JMARKOV::JPhotonPathReader reader ;
74 reader.open(it->c_str()) ;
75 if( !reader.is_open() ) {
76 cerr << "FATAL ERROR: unable to open input file '" << *it << "'." << endl ;
77 exit(1) ;
78 }
79 int nread = 0 ;
80 JMARKOV::JPhotonPath* p = NULL ;
81 while( reader.hasNext() && nread_total<=npaths_to_draw ) {
82 p = reader.next() ;
83 paths.push_back(*p) ;
84 ++nread ;
85 ++nread_total ;
86 }
87 cout << "--> read " << nread << " paths." << endl ;
88 reader.close() ;
89 }
90 cout << endl ;
91 cout << "Done reading paths." << endl ;
92 cout << endl ;
93
94 if( nread_total == 0 ) {
95 cerr << "FATAL ERROR: could not read any JPhotonPaths from the file(s)." << endl ;
96 exit(1) ;
97 }
98
99 cout << "Done reading files. Read " << nread_total << " paths from it." << endl ;
100 cout << endl ;
101
102 // variables to find out how large our view has to be
103 double xmin = 1.0/0.0 ;
104 double xmax = -1.0/0.0 ;
105 double ymin = 1.0/0.0 ;
106 double ymax = -1.0/0.0 ;
107 double zmin = 1.0/0.0 ;
108 double zmax = -1.0/0.0 ;
109
110 // location of the source and target
111 JGEOMETRY3D::JPosition3D xsource ;
112 JGEOMETRY3D::JPosition3D xtarget ;
113
114 TRandom3* ran = new TRandom3(0) ;
115
116 // limit number of paths to draw
117 npaths_to_draw = min(npaths_to_draw,nread_total) ;
118
119 vector<TPolyLine3D*> to_draw ;
120 for( int i=0 ; i<npaths_to_draw ; ++i ) {
121 int nvert = paths[i].size() ; // number of vertices
122
123 // create graphics object TPolyLine that will be used to draw the path
124 TPolyLine3D* pl = new TPolyLine3D(nvert) ;
125 pl->SetLineWidth(1) ;
126 pl->SetLineColor(kCyan) ;
127 to_draw.push_back(pl) ;
128
129 for( int j=0 ; j<nvert ; ++j ) {
130 double x = paths[i][j].getX() ;
131 double y = paths[i][j].getY() ;
132 double z = paths[i][j].getZ() ;
133 if( x>xmax ) xmax = x ;
134 if( y>ymax ) ymax = y ;
135 if( x<xmin ) xmin = x ;
136 if( y<ymin ) ymin = y ;
137 if( z<zmin ) zmin = z ;
138 if( z>zmax ) zmax = z ;
139 pl->SetPoint( j, x, y, z ) ;
140 }
141 xsource = paths[i][0] ;
142 xtarget = paths[i][nvert-1] ;
143 }
144
145 // adjust x,y, z min and max such that the length of the x, y and z
146 // ranges are equal
147 double range = 0 ;
148 if( xmax-xmin > range ) range = xmax-xmin ;
149 if( ymax-ymin > range ) range = ymax-ymin ;
150 if( zmax-zmin > range ) range = zmax-zmin ;
151 range *= 1.1 ;
152 double extra ;
153 // adjust x
154 extra = 0.5*(range-xmax-xmin) ;
155 xmax += extra ;
156 xmin -= extra ;
157 // adjust y
158 extra = 0.5*(range-ymax-ymin) ;
159 ymax += extra ;
160 ymin -= extra ;
161 // adjust z
162 extra = 0.5*(range-zmax-zmin) ;
163 zmax += extra ;
164 zmin -= extra ;
165
166 cout << "DIMENSIONS: " << endl
167 << "x = [ " << xmin << " , " << xmax << " ] [m]" << endl
168 << "y = [ " << ymin << " , " << xmax << " ] [m]" << endl
169 << "z = [ " << zmin << " , " << xmax << " ] [m]" << endl ;
170 cout << endl ;
171
172 // prepare TH3F to show the local path density
173 cout << "Creating and filling TH3F" << endl ;
174 int nbinsh3 = 200 ;
175 TH3F* h3 = new TH3F("hPathDensity_Z",
176 "Local path density;X [m];Y [m];Z [m]",
177 nbinsh3,xmin,xmax,
178 nbinsh3,ymin,ymax,
179 nbinsh3,zmin,zmax ) ;
180 h3->SetOption("colz") ;
181 for( Int_t zbin=1 ; zbin<=nbinsh3 ; ++zbin ) {
182 double z = h3->GetZaxis()->GetBinCenter(zbin) ;
183 for( int i=0 ; i<nread_total ; ++i ) {
184 vector<JGEOMETRY3D::JPosition3D> coords = paths[i].getPointsWithZ(z) ;
185 for( vector<JGEOMETRY3D::JPosition3D>::iterator it=coords.begin() ; it!=coords.end() ; ++it ) {
186 h3->Fill( it->getX(), it->getY(), it->getZ() ) ;
187 }
188 }
189 }
190 cout << "TH3F filled!" << endl ;
191 cout << endl ;
192
193 // prepare another TH3F to show the local path density
194 cout << "Creating and filling TH3F" << endl ;
195 TH3F* h3X = new TH3F("hPathDensity_X",
196 "Local path density;X [m];Y [m];Z [m]",
197 nbinsh3,xmin,xmax,
198 nbinsh3,ymin,ymax,
199 nbinsh3,zmin,zmax ) ;
200 h3X->SetOption("colz") ;
201 for( Int_t xbin=1 ; xbin<=nbinsh3 ; ++xbin ) {
202 double x = h3X->GetXaxis()->GetBinCenter(xbin) ;
203 for( int i=0 ; i<nread_total ; ++i ) {
204 vector<JGEOMETRY3D::JPosition3D> coords = paths[i].getPointsWithX(x) ;
205 for( vector<JGEOMETRY3D::JPosition3D>::iterator it=coords.begin() ; it!=coords.end() ; ++it ) {
206 h3X->Fill( it->getX(), it->getY(), it->getZ() ) ;
207 }
208 }
209 }
210 cout << "TH3F filled!" << endl ;
211 cout << endl ;
212
213 // prepare another TH3F to show the local path density
214 /*cout << "Creating and filling TH3F" << endl ;
215 nbinsh3 = 80 ;
216 TH3F* h3sphere = new TH3F("hPathDensity_sphere",
217 "Local path density;X [m];Y [m];Z [m]",
218 nbinsh3,xmin,xmax,
219 nbinsh3,ymin,ymax,
220 nbinsh3,zmin,zmax ) ;
221 h3sphere->SetOption("colz") ;
222 for( Int_t xbin=1 ; xbin<=nbinsh3 ; ++xbin ) {
223 double x = h3sphere->GetXaxis()->GetBinCenter(xbin) ;
224 cout << "x = " << x << endl ;
225 for( Int_t ybin=1 ; ybin<=nbinsh3 ; ++ybin ) {
226 double y = h3sphere->GetYaxis()->GetBinCenter(ybin) ;
227 for( Int_t zbin=1 ; zbin<=nbinsh3 ; ++zbin ) {
228 double z = h3sphere->GetZaxis()->GetBinCenter(zbin) ;
229 Int_t bin = h3sphere->GetBin(xbin,ybin,zbin) ;
230 for( int i=0 ; i<nread_total ; ++i ) {
231 JPosition3D pos(x,y,z) ;
232 if ( paths[i].hitsSphere( pos, range/nbinsh3 ) ) {
233 h3sphere->Fill( x,y,z ) ;
234 }
235 }
236 }
237 }
238 }
239 cout << "TH3F filled!" << endl ;
240 cout << endl ;*/
241
242 // make a canvas that shows the 3D paths
243 double csize = 2048 ;
244 TCanvas* cDisplay = new TCanvas("cDisplay","Display",10,10,csize,csize) ;
245 TPad* pad1 = new TPad("pad1","pad1",0,0,1,1,kBlue+4) ;
246 pad1->Draw() ;
247 pad1->cd() ;
248 TView3D* view = (TView3D*)TView::CreateView(1) ;
249 //double scale = 1.0/sqrt(2.0) ; // scale to use more of the visible area
250 view->SetRange(xmin,ymin,zmin,
251 xmax,ymax,zmax ) ;
252 //TAxis3D rulers ;
253 //rulers.Draw();
254
255 // draw source and target position
256 TPolyMarker3D* pm = new TPolyMarker3D(2,1) ;
257 pm->SetPoint(0, xsource.getX(), xsource.getY(), xsource.getZ() ) ;
258 pm->SetPoint(1, xtarget.getX(), xtarget.getY(), xtarget.getZ() ) ;
259 pm->SetMarkerSize(2) ;
260 pm->SetMarkerColor(kRed) ;
261 pm->SetMarkerStyle(20) ;
262 pm->Draw() ;
263
264 for( vector<TPolyLine3D*>::iterator it=to_draw.begin() ; it!=to_draw.end() ; ++it ) {
265 (*it)->Draw() ;
266 }
267
268 cDisplay->SaveAs("display.png") ;
269 cDisplay->SaveAs("display.pdf") ;
270 cout << "Output in 'display.png'." << endl ;
271 cout << endl ;
272
273 /*{
274 // delete the file if it is there (otherwise there will be complaints)
275 gSystem->Unlink("myGif.gif") ;
276 const int nrot = 40 ;
277 for( int i=0 ; i<nrot ; ++i ) {
278 view->RotateView( i*360.0/nrot,90 ) ;
279 if( i==nrot-1 ) {
280 cDisplay->Print("myGif.gif++10++") ;
281 } else {
282 cDisplay->Print("myGif.gif+10") ;
283 }
284 }
285
286 }
287\author mjongen
288 */
289 TFile* fout = new TFile(ofname.c_str(),"recreate") ;
290 cDisplay->Write() ;
291 h3->Write() ;
292 h3X->Write() ;
293 //h3sphere->Write() ;
294 fout->Close() ;
295 cout << "Output written to '" << ofname << "'." << endl ;
296 cout << endl ;
297 cout << "Done!" << endl ;
298}
make_field
#define make_field(A,...)
macro to convert parameter to JParserTemplateElement object
Definition JParser.hh:2142
JGEOMETRY3D::JPosition3D
Data structure for position in three dimensions.
Definition JPosition3D.hh:38
JGEOMETRY3D::JVector3D::getY
double getY() const
Get y position.
Definition JVector3D.hh:104
JGEOMETRY3D::JVector3D::getZ
double getZ() const
Get z position.
Definition JVector3D.hh:115
JGEOMETRY3D::JVector3D::getX
double getX() const
Get x position.
Definition JVector3D.hh:94
JMARKOV::JPhotonPath
A photon path.
Definition JPhotonPath.hh:38
JPARSER::JParser
Utility class to parse command line options.
Definition JParser.hh:1698
JPARSER::JParser::read
int read(const int argc, const char *const argv[])
Parse the program's command line options.
Definition JParser.hh:1992
GAUSS_LEGENDRE::xmax
const double xmax
Definition JQuadrature.cc:24
GAUSS_LEGENDRE::xmin
const double xmin
Definition JQuadrature.cc:23
JTOOLS::j
int j
Definition JPolint.hh:801
Generated by doxygen 1.12.0