2 #include "GenFit/FitStatus.h"
3 #include "GenFit/GFRaveVertexFactory.h"
6 #include "StEvent/StEvent.h"
7 #include "StEvent/StFttCollection.h"
8 #include "StEvent/StFttCluster.h"
12 float DEGS_TO_RAD = 3.14159f/180.0f;
19 const static bool verbose =
false;
22 void vert( ofstream &of,
float x,
float y,
float z ){
23 of <<
"v " << x <<
" " << y <<
" " << z << endl;
26 void vert( ofstream &of, TVector3 v ){
27 of <<
"v " << v.X() <<
" " << v.Y() <<
" " << v.Z() << endl;
32 void sphere(TVector3 pt,
float radius,
int nLatitude,
int nLongitude, ofstream &fout){
35 int nPitch = nLongitude + 1;
37 float pitchInc = (180. / (float)nPitch) * DEGS_TO_RAD;
38 float rotInc = (360. / (float)nLatitude) * DEGS_TO_RAD;
41 vert(fout, pt.X(), pt.Y()+radius, pt.Z());
42 vert(fout, pt.X(), pt.Y()-radius, pt.Z());
44 int fVert = numVertices;
45 for(p=1; p<nPitch; p++) {
46 out = fabs(radius * sin((
float)p * pitchInc));
47 y = radius * cos(p * pitchInc);
48 for(s=0; s<nLatitude; s++) {
49 x = out * cos(s * rotInc);
50 z = out * sin(s * rotInc);
51 fout << TString::Format(
"v %g %g %g\n", x+pt.X(), y+pt.Y(), z+pt.Z()).Data();
58 for(p=1; p<nPitch-1; p++) {
59 for(s=0; s<nLatitude; s++) {
61 j = (s==nLatitude-1) ? i-nLatitude : i;
62 fout << TString::Format(
"f %d %d %d %d\n",
63 (i+1-nLatitude)+fVert, (j+2-nLatitude)+fVert, (j+2)+fVert, (i+1)+fVert);
69 int offLastVerts = fVert + (nLatitude * (nLongitude-1));
70 for(s=0; s<nLatitude; s++)
72 j = (s==nLatitude-1) ? -1 : s;
73 fout << TString::Format(
"f %d %d %d\n", fVert-1, (j+2)+fVert, (s+1)+fVert ).Data();
74 fout << TString::Format(
"f %d %d %d\n", fVert, (s+1)+offLastVerts, (j+2)+offLastVerts).Data();
79 void face( ofstream &fout,
int anchor,
int v0,
int v1,
int v2 ){
80 fout << TString::Format(
"f %d %d %d\n", anchor+v0, anchor+v1, anchor+v2).Data();
82 void face( ofstream &fout,
int anchor,
int v0,
int v1,
int v2,
int v3 ){
83 fout << TString::Format(
"f %d %d %d %d\n", anchor+v0, anchor+v1, anchor+v2, anchor+v3).Data();
86 void prism( ofstream &fout, TVector3 po, TVector3 ds ){
87 int fVert = numVertices;
89 vert( fout, po.X(), po.Y(), po.Z() );
90 vert( fout, po.X() + ds.X(), po.Y(), po.Z() );
91 vert( fout, po.X() + ds.X(), po.Y() + ds.Y(), po.Z() );
92 vert( fout, po.X(), po.Y() + ds.Y(), po.Z() );
95 vert( fout, po.X(), po.Y(), po.Z() + ds.Z() );
96 vert( fout, po.X() + ds.X(), po.Y(), po.Z() + ds.Z() );
97 vert( fout, po.X() + ds.X(), po.Y() + ds.Y(), po.Z() + ds.Z() );
98 vert( fout, po.X(), po.Y() + ds.Y(), po.Z() + ds.Z() );
100 face( fout, fVert, 1, 2, 3, 4 );
101 face( fout, fVert, 5, 6, 7, 8 );
102 face( fout, fVert, 1, 2, 6, 5 );
103 face( fout, fVert, 2, 3, 7, 6 );
104 face( fout, fVert, 3, 4, 8, 7 );
105 face( fout, fVert, 1, 4, 8, 5 );
109 void tri( ofstream &fout, TVector3 v0,
float dz,
float w,
float h,
float phi ) {
111 int fVert = numVertices;
113 float yp = v0.Y() + (sin( phi + dphi ) + cos( phi + dphi ) ) * h;
114 float yn = v0.Y() + (sin( phi - dphi ) + cos( phi - dphi ) ) * h;
115 float xp = v0.X() + (cos( phi + dphi ) - sin( phi + dphi)) * h;
116 float xn = v0.X() + (cos( phi - dphi ) - sin( phi - dphi)) * h;
118 vert( fout, v0.X(), v0.Y(), v0.Z() );
119 vert( fout, xp, yp, v0.Z() );
120 vert( fout, xn, yn, v0.Z() );
123 vert( fout, v0.X(), v0.Y(), v0.Z() - dz );
124 vert( fout, xp, yp, v0.Z() - dz );
125 vert( fout, xn, yn, v0.Z() - dz );
127 face( fout, fVert, 1, 2, 3 );
128 face( fout, fVert, 4, 5, 6 );
129 face( fout, fVert, 1, 2, 5, 4 );
130 face( fout, fVert, 2, 3, 6, 5 );
131 face( fout, fVert, 3, 1, 4, 6 );
135 static TVector3 projectAsStraightLine( genfit::Track * t,
float z0,
float z1,
float zf,
float * cov, TVector3 &mom ) {
136 TVector3 tv3A = trackPosition( t, z0, cov, mom );
137 TVector3 tv3B = trackPosition( t, z1, cov, mom );
140 LOG_INFO <<
"Straight Line Projection using" << endm;
141 LOG_INFO <<
"A.x = " << tv3A.X() << endm;
142 LOG_INFO <<
"B.x = " << tv3B.X() << endm;
145 double dxdz = ( tv3B.X() - tv3A.X() ) / ( tv3B.Z() - tv3A.Z() );
146 double dydz = ( tv3B.Y() - tv3A.Y() ) / ( tv3B.Z() - tv3A.Z() );
148 double dx = dxdz * ( zf - z1 );
149 double dy = dydz * ( zf - z1 );
150 TVector3 r( tv3B.X() + dx, tv3B.Y() + dy, zf );
155 static TVector3 trackPosition( genfit::Track * t,
float z,
float * cov, TVector3 &mom ){
159 auto plane = genfit::SharedPlanePtr(
161 new genfit::DetPlane(TVector3(0, 0, z), TVector3(1, 0, 0), TVector3(0, 1, 0) )
164 genfit::MeasuredStateOnPlane tst = t->getFittedState(iPoint);
165 auto TCM = t->getCardinalRep()->get6DCov(tst);
167 t->getCardinalRep()->extrapolateToPlane(tst, plane,
false,
true);
169 TCM = t->getCardinalRep()->get6DCov(tst);
172 float x = tst.getPos().X();
173 float y = tst.getPos().Y();
174 float _z = tst.getPos().Z();
176 mom.SetXYZ( tst.getMom().X(), tst.getMom().Y(), tst.getMom().Z() );
179 cov[0] = TCM(0,0); cov[1] = TCM(1,0); cov[2] = TCM(2,0);
180 cov[3] = TCM(0,1); cov[4] = TCM(1,1); cov[5] = TCM(2,1);
181 cov[6] = TCM(0,2); cov[7] = TCM(1,2); cov[8] = TCM(2,2);
185 return TVector3( x, y, _z );
186 }
catch ( genfit::Exception e ){
187 LOG_INFO <<
"Track projection Failed from trackPoint " << iPoint <<
" E: " << e.what() << endm;
188 return TVector3( -990, -990, -990 );
192 return TVector3( -990, -990, -990 );
196 static TVector3 trackPosition( genfit::Track * t,
float z ){
199 return trackPosition( t, z, cov, mom);
203 void output_ftt_strips(
207 fout <<
"\n" << endl;
208 fout <<
"o fttStrips" << endl;
209 fout <<
"usemtl stgc_hits\n" << endl;
210 float pz[] = {280.90499, 303.70498, 326.60501, 349.40499};
214 for (
size_t i = 0; i <
event->fttCollection()->numberOfClusters(); i++ ){
215 StFttCluster* c =
event->fttCollection()->clusters()[i];
216 if ( c->nStrips() < 2 )
continue;
217 float dw = 0.05, dlh = 60.0, dlv = 60.0;
218 float mx = 0.0, my = 0.0;
219 float sx = 1.0, sy = 1.0;
222 if ( c->quadrant() == kFttQuadrantA ){
225 }
else if ( c->quadrant() == kFttQuadrantB ){
226 mx = 10.16*SCALE; my = 0.0*SCALE;
230 }
else if ( c->quadrant() == kFttQuadrantC ){
231 mx = -10.16*SCALE ; my = -00.0*SCALE;
232 sx = -1.0; sy = -1.0;
233 dlh = -dlh; dlv = -dlv;
235 }
else if ( c->quadrant() == kFttQuadrantD ){
240 cp.SetZ( -pz[ c->plane() ] * SCALE );
241 if ( c->orientation() == kFttHorizontal ){
242 cp.SetY( my + sy * c->x()/10.0 * SCALE );
244 prism( fout, cp, TVector3( dlh * SCALE, dw, dw ) );
245 }
else if ( c->orientation() == kFttVertical ){
246 cp.SetX( mx + sx * c->x()/10.0 * SCALE );
248 prism( fout, cp, TVector3( dw, dlv * SCALE, dw ) );
254 void output( std::string filename,
256 std::vector< Seed_t> seeds,
257 std::vector< genfit::Track *> tracks,
258 const std::vector< genfit::GFRaveVertex *> &vertices,
259 std::vector<TVector3> &fttHits,
260 std::vector<TVector3> &fstHits,
261 std::vector<TVector3> &fcsPreHits,
262 std::vector<TVector3> &fcsClusters,
263 std::vector<float> &fcsClusterEnergy
266 LOG_INFO <<
"Writing: " << filename << endm;
269 ofstream ofile( (filename +
".obj" ).c_str() );
271 ofile <<
"\nmtllib materials.mtl\n\n" << endl;
274 output_ftt_strips( ofile, event );
278 if ( vertices.size() > 0 ){
279 ofile <<
"o FwdVertices" << endl;
280 for (
auto v : vertices ) {
281 startPos.SetXYZ( v->getPos().X(), v->getPos().Y(), v->getPos().Z() );
282 sphere( TVector3( v->getPos().X() * SCALE, v->getPos().Y() * SCALE, -v->getPos().Z() * SCALE ), 0.5, 10, 10, ofile );
288 LOG_INFO <<
"Viz has " << fttHits.size() <<
" FTT Hits" << endm;
290 if ( fttHits.size() > 0 ){
291 ofile <<
"\n" << endl;
292 ofile <<
"o fttHits" << endl;
293 ofile <<
"usemtl stgc_hits\n" << endl;
294 for (
auto p : fttHits ){
295 sphere( TVector3( p.X() * SCALE, p.Y() * SCALE, -p.Z() * SCALE ), 0.15, 12, 12, ofile );
301 LOG_INFO <<
"Viz has " << fstHits.size() <<
" FST Hits" << endm;
303 if ( fstHits.size() > 0 ) {
304 ofile <<
"\n" << endl;
305 ofile <<
"o fstHits" << endl;
306 ofile <<
"usemtl fst_hits\n" << endl;
307 for (
auto p : fstHits ){
309 float fstphi = TMath::ATan2( p.Y(), p.X() );
310 printf(
"FST PHI: %f \n", fstphi );
312 sphere( TVector3( p.X() * SCALE, p.Y() * SCALE, -p.Z() * SCALE ), 0.3, 10, 10, ofile );
318 LOG_INFO <<
"Viz has " << fcsPreHits.size() <<
" EPD Hits" << endm;
320 if ( fcsPreHits.size() > 0 ){
321 ofile <<
"\n" << endl;
322 ofile <<
"o epd" << endl;
323 ofile <<
"usemtl fcs_hits\n" << endl;
324 for (
auto p : fcsPreHits ){
326 sphere( TVector3( p.X() * SCALE, p.Y() * SCALE, -p.Z() * SCALE ), 0.25, 10, 10, ofile );
331 LOG_INFO <<
"Viz has " << fcsClusters.size() <<
" FCS Hits" << endm;
333 if ( fcsClusters.size() > 0 ){
334 ofile <<
"\n" << endl;
335 ofile <<
"o fcs" << endl;
336 ofile <<
"usemtl fcs_hits\n" << endl;
338 for (
auto p : fcsClusters ){
340 TVector3 ds = TVector3( 1, 1, fcsClusterEnergy[i] );
341 prism( ofile, TVector3( p.X() * SCALE, p.Y() * SCALE, -p.Z() * SCALE ), ds );
348 LOG_INFO <<
"Viz has " << seeds.size() <<
" seeds" << endm;
350 if ( seeds.size() > 0 ){
351 ofile <<
"\n\no FwdSeeds\n" << endl;
352 ofile <<
"usemtl seeds\n" << endl;
354 for (
auto s : seeds ) {
355 size_t vStart = numVertices;
358 vert( ofile, h->getX() * SCALE, h->getY() * SCALE, -h->getZ() * SCALE );
363 for (
size_t i = vStart; i < numVertices; i++){
376 LOG_INFO <<
"Viz has " << tracks.size() <<
" tracks Hits" << endm;
378 if ( tracks.size() > 0 ){
379 ofile <<
"\n\no FwdTracks\n" << endl;
380 ofile <<
"usemtl tracks\n" << endl;
382 for (
auto t : tracks ) {
383 size_t vStart = numVertices;
387 for (
float z = startPos.Z(); z < 875; z += zStep ){
388 TVector3 point = trackPosition( t, z );
389 if ( point.X() < -900 && point.Y() < -900 )
break;
390 if ( point.X() < -90 && point.Y() < -90 ) { z+= 50;
continue;}
392 vert( ofile, point.X() * SCALE, point.Y() * SCALE, -point.Z() * SCALE );
397 for (
size_t i = vStart; i < numVertices; i++){