/* This file is part of MAUS: http://micewww.pp.rl.ac.uk:8080/projects/maus
*
* MAUS is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* MAUS is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MAUS. If not, see .
*
*/
#include "src/utilities/event-viewer/EVHepRepExporter.hh"
#include
#include
#include
#include
#include "src/legacy/Config/MiceModule.hh"
namespace EventViewer {
EVHepRepExporter::EVHepRepExporter(EVEvent evEvent) {
event = evEvent;
}
void EVHepRepExporter::Export(int display) {
HepRepXMLWriter* writer = new HepRepXMLWriter();
std::stringstream ssOutFileName;
// - retrieve destination directory
char *destDir = getenv("EV_DEST_DIR");
if (destDir != NULL)
ssOutFileName << destDir << "/";
else
ssOutFileName << "./";
// - open output heprep file
ssOutFileName << "MICEEvent_" << event.runNumber << "_" << event.spillNumber
<< "_" << event.eventNumber << ".heprep";
std::string dataFileName = "data.heprep";
writer->open(dataFileName.c_str());
// - define tof points
writer->addType("TOFPoints", 0);
writer->addInstance();
writer->addAttValue("DrawAs", "Point");
writer->addAttValue("MarkName", "Box");
writer->addAttValue("MarkColor", 0., 1., 1.);
for (int point = 0; point < 3; ++point) {
double x = event.tofPoints[point][0];
double y = event.tofPoints[point][1];
double z = event.tofPoints[point][2];
if (z != 5673390) {
writer->addPrimitive();
writer->addPoint(x, y, z);
}
}
// - define track points
writer->addType("TrackPoints", 0);
writer->addInstance();
writer->addAttValue("DrawAs", "Point");
// - plot TKU track points in HepRApp
writer->addType("TKUTrackPoints", 1);
writer->addInstance();
writer->addAttValue("MarkName", "Cross");
writer->addAttValue("MarkColor", 1., 1., 0.);
for (int point = 0; point < 5; ++point) {
writer->addType("TKUTrackPoint", 2); // add layer of depth to apply attibutes
writer->addInstance();
double x = event.scifiUSTrackerPoints[point][0];
double y = event.scifiUSTrackerPoints[point][1];
double z = event.scifiUSTrackerPoints[point][2];
// attributes - momentum components
writer->addAttDef("Px", "Track momentum x component", "Physics", "GeV");
writer->addAttDef("Py", "Track momentum y component", "Physics", "GeV");
writer->addAttDef("Pz", "Track momentum z component", "Physics", "GeV");
writer->addAttValue("Px", event.scifiUSTrackerPointsMomenta[point][0]);
writer->addAttValue("Py", event.scifiUSTrackerPointsMomenta[point][1]);
writer->addAttValue("Pz", event.scifiUSTrackerPointsMomenta[point][2]);
if (z != 5673390) {
writer->addPrimitive();
writer->addPoint(x, y, z);
}
}
// - plot TKD track points in HepRApp
writer->addType("TKDTrackPoints", 1);
writer->addInstance();
writer->addAttValue("MarkName", "Cross");
writer->addAttValue("MarkColor", 1., 1., 0.);
for (int point = 0; point < 5; ++point) {
writer->addType("TKDTrackPoint", 2); // add layer of depth to apply attibutes
writer->addInstance();
double x = event.scifiDSTrackerPoints[point][0];
double y = event.scifiDSTrackerPoints[point][1];
double z = event.scifiDSTrackerPoints[point][2];
// attributes - momentum components
writer->addAttDef("Px", "Track momentum x component", "Physics", "GeV");
writer->addAttDef("Py", "Track momentum y component", "Physics", "GeV");
writer->addAttDef("Pz", "Track momentum z component", "Physics", "GeV");
writer->addAttValue("Px", event.scifiUSTrackerPointsMomenta[point][0]);
writer->addAttValue("Py", event.scifiUSTrackerPointsMomenta[point][1]);
writer->addAttValue("Pz", event.scifiUSTrackerPointsMomenta[point][2]);
if (z != 5673390) {
writer->addPrimitive();
writer->addPoint(x, y, z);
}
}
// - define space points
writer->addType("SpacePoints", 0);
writer->addInstance();
writer->addAttValue("DrawAs", "Point");
// - plot TKU space points in HepRApp
writer->addType("TKUSpacePoints", 1);
writer->addInstance();
writer->addAttValue("MarkName", "Cross");
writer->addAttValue("MarkColor", 1., 0.3, 0.);
for (std::vector::iterator it = event.scifiUSTrackerSpacePoints.begin(),
end = event.scifiUSTrackerSpacePoints.end(); it != end; ++it) {
writer->addType("TKUSpacePoint", 2); // add layer of depth to apply attibutes
writer->addInstance();
double x = it->x();
double y = it->y();
double z = it->z();
writer->addAttDef("alpha P0", "Cluster info for plane 0 - alpha", "Physics", "channel");
writer->addAttDef("z P0", "Cluster info for plane 0 - z", "Physics", "mm");
writer->addAttDef("alpha P1", "Cluster info for plane 1", "Physics", "channel");
writer->addAttDef("z P1", "Cluster info for plane 1 - z", "Physics", "mm");
writer->addAttDef("alpha P2", "Cluster info for plane 2", "Physics", "channel");
writer->addAttDef("z P2", "Cluster info for plane 2 - z", "Physics", "mm");
unsigned int station = it - event.scifiUSTrackerSpacePoints.begin();
if (station < event.scifiUSTrackerClusters[0][0].size()) {
writer->addAttValue("alpha P0", event.scifiUSTrackerClusters[0][0][station]);
writer->addAttValue("z P0", event.scifiUSTrackerClusters[0][1][station]);
}
if (station < event.scifiUSTrackerClusters[1][0].size()) {
writer->addAttValue("alpha P1", event.scifiUSTrackerClusters[1][0][station]);
writer->addAttValue("z P1", event.scifiUSTrackerClusters[1][1][station]);
}
if (station < event.scifiUSTrackerClusters[2][0].size()) {
writer->addAttValue("alpha P2", event.scifiUSTrackerClusters[2][0][station]);
writer->addAttValue("z P2", event.scifiUSTrackerClusters[2][1][station]);
}
if (z != 5673390) { // should be 5673390
writer->addPrimitive();
writer->addPoint(x, y, z);
}
}
// - plot TKD space points in HepRApp
writer->addType("TKDSpacePoints", 1);
writer->addInstance();
writer->addAttValue("MarkName", "Cross");
writer->addAttValue("MarkColor", 1., 0.3, 0.);
for (std::vector::iterator it = event.scifiDSTrackerSpacePoints.begin(),
end = event.scifiDSTrackerSpacePoints.end(); it != end; ++it) {
writer->addType("TKDSpacePoint", 2); // add layer of depth to apply attibutes
writer->addInstance();
double x = it->x();
double y = it->y();
double z = it->z();
writer->addAttDef("alpha P0", "Cluster info for plane 0 - alpha", "Physics", "channel");
writer->addAttDef("z P0", "Cluster info for plane 0 - z", "Physics", "mm");
writer->addAttDef("alpha P1", "Cluster info for plane 1", "Physics", "channel");
writer->addAttDef("z P1", "Cluster info for plane 1 - z", "Physics", "mm");
writer->addAttDef("alpha P2", "Cluster info for plane 2", "Physics", "channel");
writer->addAttDef("z P2", "Cluster info for plane 2 - z", "Physics", "mm");
unsigned int station = it - event.scifiDSTrackerSpacePoints.begin();
if (station < event.scifiDSTrackerClusters[0][0].size()) {
writer->addAttValue("alpha P0", event.scifiDSTrackerClusters[0][0][station]);
writer->addAttValue("z P0", event.scifiDSTrackerClusters[0][1][station]);
}
if (station < event.scifiDSTrackerClusters[1][0].size()) {
writer->addAttValue("alpha P1", event.scifiDSTrackerClusters[1][0][station]);
writer->addAttValue("z P1", event.scifiDSTrackerClusters[1][1][station]);
}
if (station < event.scifiDSTrackerClusters[2][0].size()) {
writer->addAttValue("alpha P2", event.scifiDSTrackerClusters[2][0][station]);
writer->addAttValue("z P2", event.scifiDSTrackerClusters[2][1][station]);
}
if (z != 5673390) { // should be 5673390
writer->addPrimitive();
writer->addPoint(x, y, z);
}
}
// - define straight track
writer->addType("StraightTrack", 0);
writer->addInstance();
writer->addAttValue("DrawAs", "Line");
writer->addAttValue("LineColor", 1., 0., 0.);
// - plot TKU straight track(s)
writer->addType("TKUStraightTrack", 1);
writer->addInstance();
writer->addPrimitive();
for (int point = 0; point < 2; ++point) {
double x = event.scifiUSTrackerStraightTrackPoints[point][0];
double y = event.scifiUSTrackerStraightTrackPoints[point][1];
double z = event.scifiUSTrackerStraightTrackPoints[point][2];
if (z != 5673390) { // ToDo: define dummy value at one place and use through out
writer->addPoint(x, y, z);
}
}
// - load parent geometry mice module (needed for helical track transformations)
std::stringstream modulePath;
if (getenv("EV_GEOM_HOME"))
modulePath << getenv("EV_GEOM_HOME") << "/ParentGeometryFile.dat";
else
modulePath << getenv("MAUS_ROOT_DIR") << "/files/geometry/download/ParentGeometryFile.dat";
// - Adam's fix (redirecting stdout/stderr)
std::ostream* stdout = new std::ofstream();
stdout->rdbuf(std::cout.rdbuf());
std::ostream* stderr = new std::ofstream();
stderr->rdbuf(std::cerr.rdbuf());
try {
mm = new MiceModule(modulePath.str());
std::cout.rdbuf(stdout->rdbuf());
std::cerr.rdbuf(stderr->rdbuf());
} catch (MAUS::Exceptions::Exception e) {
std::cerr << "Caught exception while building MiceModules" << std::endl;
e.what();
e.Print();
modulePath.str("");
// this file should be included in the release so no try/catch neccessary
modulePath << getenv("MAUS_ROOT_DIR") << "/src/utilities/event-viewer/ParentGeometryFile.dat";
mm = new MiceModule(modulePath.str());
std::cout.rdbuf(stdout->rdbuf());
std::cerr.rdbuf(stderr->rdbuf());
}
// --- fetch US tracker global geometry
mmTrackerUS = mm->findModulesByName("ParentGeometryFile.dat/Tracker0.dat");
// - get TrackerRef0 position
double trackerRef0posX = mmTrackerUS[0]->daughter(0)->globalPosition().x();
double trackerRef0posY = mmTrackerUS[0]->daughter(0)->globalPosition().y();
double trackerRef0posZ = mmTrackerUS[0]->daughter(0)->globalPosition().z();
// - get rotation matrix
double xx = mmTrackerUS[0]->globalRotation().xx();
double xy = mmTrackerUS[0]->globalRotation().xy();
double xz = mmTrackerUS[0]->globalRotation().xz();
double yx = mmTrackerUS[0]->globalRotation().yx();
double yy = mmTrackerUS[0]->globalRotation().yy();
double yz = mmTrackerUS[0]->globalRotation().yz();
double zx = mmTrackerUS[0]->globalRotation().zx();
double zy = mmTrackerUS[0]->globalRotation().zy();
double zz = mmTrackerUS[0]->globalRotation().zz();
// - fill matrix elements (for VectorTransform)
std::vector matrixElements;
matrixElements.push_back(xx);
matrixElements.push_back(xy);
matrixElements.push_back(xz);
matrixElements.push_back(yx);
matrixElements.push_back(yy);
matrixElements.push_back(yz);
matrixElements.push_back(zx);
matrixElements.push_back(zy);
matrixElements.push_back(zz);
// - dimensions (y coordinate is length)
Hep3Vector dimensions = mmTrackerUS[0]->dimensions();
double trackerLength = dimensions.y();
// - plot TKU helical track(s)
writer->addType("TKUHelicalTrackTest", 1);
writer->addInstance();
writer->addPrimitive();
double circle_x0 = event.scifiUSTrackerHelicalTrackParameters[0];
double circle_y0 = event.scifiUSTrackerHelicalTrackParameters[1];
double radius = event.scifiUSTrackerHelicalTrackParameters[2];
double dsdz = event.scifiUSTrackerHelicalTrackParameters[3];
double line_dz_c = event.scifiUSTrackerHelicalTrackParameters[4];
double posZ = event.scifiUSTrackerHelicalTrackParameters[7];
double zLocal = posZ;
int zStep = trackerLength/300;
while (zLocal < trackerLength) {
zLocal += zStep; // z value in tracker local reference system
// x value in tracker local reference system
double upstream_helix_x_value
= circle_x0 + (radius*cos((1/radius)*(dsdz*zLocal+line_dz_c)));
// y value in tracker local reference system
double upstream_helix_y_value
= circle_y0 + (radius*sin((1/radius)*(dsdz*zLocal+line_dz_c)));
// transformation to global coordinate system
std::vector vSpacePoint
= {upstream_helix_x_value, upstream_helix_y_value, zLocal};
VectorTransform trans(vSpacePoint);
trans.Rotate(matrixElements, "reg");
trans.Translate(trackerRef0posX, trackerRef0posY, trackerRef0posZ, "reg");
double x = trans.GetX();
double y = trans.GetY();
double z = trans.GetZ();
writer->addPoint(x, y, z);
}
// - plot TKD straight track(s)
writer->addType("TKDStraightTrack", 1);
writer->addInstance();
writer->addPrimitive();
for (int point = 0; point < 2; ++point) {
double x = event.scifiDSTrackerStraightTrackPoints[point][0];
double y = event.scifiDSTrackerStraightTrackPoints[point][1];
double z = event.scifiDSTrackerStraightTrackPoints[point][2];
if (z != 5673390) {
writer->addPoint(x, y, z);
}
}
// --- fetch DS tracker global geometry
mmTrackerDS = mm->findModulesByName("ParentGeometryFile.dat/Tracker1.dat");
// - get TrackerRef1 position
double trackerRef1posX = mmTrackerDS[0]->daughter(0)->globalPosition().x();
double trackerRef1posY = mmTrackerDS[0]->daughter(0)->globalPosition().y();
double trackerRef1posZ = mmTrackerDS[0]->daughter(0)->globalPosition().z();
// - get rotation matrix
xx = mmTrackerDS[0]->globalRotation().xx();
xy = mmTrackerDS[0]->globalRotation().xy();
xz = mmTrackerDS[0]->globalRotation().xz();
yx = mmTrackerDS[0]->globalRotation().yx();
yy = mmTrackerDS[0]->globalRotation().yy();
yz = mmTrackerDS[0]->globalRotation().yz();
zx = mmTrackerDS[0]->globalRotation().zx();
zy = mmTrackerDS[0]->globalRotation().zy();
zz = mmTrackerDS[0]->globalRotation().zz();
// - fill matrix elements (for VectorTransform)
matrixElements.clear();
matrixElements.push_back(xx);
matrixElements.push_back(xy);
matrixElements.push_back(xz);
matrixElements.push_back(yx);
matrixElements.push_back(yy);
matrixElements.push_back(yz);
matrixElements.push_back(zx);
matrixElements.push_back(zy);
matrixElements.push_back(zz);
// - plot TKD helical track(s)
writer->addType("TKDHelicalTrackTest", 1);
writer->addInstance();
writer->addPrimitive();
circle_x0 = event.scifiDSTrackerHelicalTrackParameters[0];
circle_y0 = event.scifiDSTrackerHelicalTrackParameters[1];
radius = event.scifiDSTrackerHelicalTrackParameters[2];
dsdz = event.scifiDSTrackerHelicalTrackParameters[3];
line_dz_c = event.scifiDSTrackerHelicalTrackParameters[4];
posZ = event.scifiDSTrackerHelicalTrackParameters[7];
zLocal = posZ;
while (zLocal < trackerLength) {
zLocal += zStep; // z value in tracker local reference system
// x value in tracker local reference system
double downstream_helix_x_value
= circle_x0 + (radius*cos((1/radius)*(dsdz*zLocal+line_dz_c)));
// y value in tracker local reference system
double downstream_helix_y_value
= circle_y0 + (radius*sin((1/radius)*(dsdz*zLocal+line_dz_c)));
// transformation to global coordinate system
std::vector vSpacePoint
= {downstream_helix_x_value, downstream_helix_y_value, zLocal};
VectorTransform trans(vSpacePoint);
trans.Rotate(matrixElements, "reg");
trans.Translate(trackerRef1posX, trackerRef1posY, trackerRef1posZ, "reg");
double x = trans.GetX();
double y = trans.GetY();
double z = trans.GetZ();
writer->addPoint(x, y, z);
}
writer->close();
char *detGeomFileName = getenv("EV_DETGEOM_FILE");
// disable writing to HepRep options if geometry file is not selected
if (detGeomFileName == NULL) {
std::cerr << "WARNING: Detector geometry file not defined! Modify and source configure.sh or "
<< "choose geometry file using selection button (if using GUI). "
<< "Not writting output heprep file!" << std::endl;
} else {
std::string sOutFileName = ssOutFileName.str();
std::string sDetGeomFileName(detGeomFileName);
Concatenate(sDetGeomFileName, dataFileName, sOutFileName);
}
// - remove MAUS data heprep file
std::stringstream cmdRemove;
cmdRemove << "rm " << dataFileName;
system(cmdRemove.str().c_str());
// - open the file for viewing
if (display) {
std::cout << display << std::endl;
std::stringstream cmdDisplay;
int displayCheck = 1;
char *javaDir = getenv("EV_JAVA_DIR");
if (javaDir != NULL) {
cmdDisplay << javaDir << "/bin/java -jar ";
} else {
std::cerr << "WARNING: No java directory selected! HepRApp display unavailable! "
<< "Select directory in which your JRE can be found!" << std::endl;
displayCheck = 0;
}
char *heprappDir = getenv("EV_HEPRAPP_DIR");
if (heprappDir != NULL) {
cmdDisplay << heprappDir << "/HepRApp.jar -file ";
} else {
std::cerr << "WARNING: No HepRapp directory selected! HepRApp display unavailable! "
<< "Select directory in which your HepRApp.jar can be found!" << std::endl;
displayCheck = 0;
}
cmdDisplay << ssOutFileName.str();
std::cout << cmdDisplay.str() << std::endl;
if (displayCheck)
system(cmdDisplay.str().c_str());
}
}
int EVHepRepExporter::Concatenate(std::string& geometryFileName, std::string& dataFileName,
std::string& outFileName) {
std::ifstream geomFile;
geomFile.open(geometryFileName.c_str(), std::ifstream::in);
if (!geomFile.is_open()) {
std::cerr << "ERROR: Could not open geometry file needed for output to heprep!" << std::endl;
return 0;
}
std::ifstream dataFile;
dataFile.open(dataFileName.c_str(), std::fstream::in);
if (!dataFile.is_open()) {
std::cerr << "ERROR: Could not open data file needed for output to heprep!" << std::endl;
return 0;
}
std::ofstream outFile;
outFile.open(outFileName.c_str(), std::ifstream::out);
if (!outFile.is_open()) {
std::cerr << "ERROR: Could not open temporary file needed for output to heprep!" << std::endl;
return 0;
}
while (!geomFile.eof()) {
std::string line, target;
std::getline(geomFile, line);
target = line.substr(0, 16);
if (target != "")
outFile << line << std::endl;
}
int counter = 0;
while (!dataFile.eof()) {
++counter;
std::string line;
std::getline(dataFile, line);
if (counter > 3)
outFile << line << std::endl;
}
return 1;
}
} // ~namespace EventViewer