two trees, cleaned up
This commit is contained in:
James Szalkie
parent
6e969434da
commit
4bbb1399cc
1
.gitignore
vendored
1
.gitignore
vendored
|
|
@ -8,6 +8,7 @@ EventBuilder*
|
|||
*.err
|
||||
*.seq
|
||||
*.png
|
||||
*.pdf
|
||||
Mapper
|
||||
AnasenMS
|
||||
Armory/anasenMS
|
||||
|
|
|
|||
13
.vscode/c_cpp_properties.json
vendored
13
.vscode/c_cpp_properties.json
vendored
|
|
@ -1,5 +1,18 @@
|
|||
{
|
||||
"configurations": [
|
||||
{
|
||||
"name": "Linux",
|
||||
"includePath": [
|
||||
"${workspaceFolder}/**",
|
||||
"/opt/root-6.36.06/include",
|
||||
"/home/jamesszalkie/anasen/Armory"
|
||||
],
|
||||
"defines": [],
|
||||
"compilerPath": "/usr/bin/g++",
|
||||
"cStandard": "c11",
|
||||
"cppStandard": "c++17",
|
||||
"intelliSenseMode": "gcc-x64"
|
||||
},
|
||||
{
|
||||
"name": "Hades",
|
||||
"includePath": [
|
||||
|
|
|
|||
|
|
@ -25,16 +25,16 @@ clean :
|
|||
|
||||
Mapper : Mapper.cpp ../mapping.h ClassDet.h
|
||||
@echo "--------- making Mapper"
|
||||
$(CC) $(COPTS) -o Mapper Mapper.cpp $(ROOTLIBS)
|
||||
$(CC) $(COPTS) $(ROOTCFLAGS) -o Mapper Mapper.cpp $(ROOTLIBS)
|
||||
|
||||
AnasenMS : constant.h Isotope.h ClassTransfer.h ClassSX3.h ClassPW.h ClassAnasen.h anasenMS.cpp
|
||||
@echo "--------- making ANASEN Monte Carlo"
|
||||
$(CC) $(COPTS) -o AnasenMS anasenMS.cpp $(ROOTLIBS)
|
||||
$(CC) $(COPTS) $(ROOTCFLAGS) -o AnasenMS anasenMS.cpp $(ROOTLIBS) -lEve -lGui -lGeom
|
||||
|
||||
EventBuilder : EventBuilder.cpp ClassData.h fsuReader.h Hit.h
|
||||
@echo "--------- making EventBuilder"
|
||||
$(CC) $(COPTS) -o EventBuilder EventBuilder.cpp $(ROOTLIBS)
|
||||
|
||||
anasenMS: anasenMS.cpp
|
||||
$(CXX) $(CXXFLAGS) anasenMS.cpp -o anasenMS $(ROOTLIBS)
|
||||
#anasenMS: anasenMS.cpp
|
||||
# $(CXX) $(CXXFLAGS) anasenMS.cpp -o anasenMS $(ROOTLIBS)
|
||||
|
||||
|
|
|
|||
|
|
@ -1,13 +1,19 @@
|
|||
#include "TRandom.h" // ROOT random number generators, gRandom
|
||||
#include "TFile.h" // ROOT file I/O
|
||||
#include "TTree.h" // ROOT tree storage
|
||||
#include "TH1.h" // 1D histograms (not directly used here but common in analyzers)
|
||||
#include "TH1.h" // 1D histograms
|
||||
#include "TH2.h" // 2D histograms
|
||||
#include "TStyle.h" // ROOT plotting style controls
|
||||
#include "TCanvas.h" // ROOT canvas drawing
|
||||
#include "TBenchmark.h" // timing measurement#include "TGraph.h" // for energy loss interpolation
|
||||
#include "TBenchmark.h" // timing measurement
|
||||
#include "TGraph.h" // for energy loss interpolation
|
||||
#include <cstring>
|
||||
#include "TApplication.h" // ROOT app loop for TEve
|
||||
#include "TEveManager.h"
|
||||
#include "TEvePointSet.h"
|
||||
#include "ClassTransfer.h" // Reaction kinematics and MC event generation
|
||||
#include "ClassAnasen.h" // ANASEN detector model classes (SX3, PW, etc.)
|
||||
#include "vis_helpers.h" // Visualization utilities for TEve
|
||||
|
||||
//======== Generate light particle based on reaction
|
||||
// calculate real and reconstructed tracks and Q-value uncertainty
|
||||
|
|
@ -76,6 +82,20 @@ int main(int argc, char **argv){
|
|||
int nExA = ExAList.size();
|
||||
int nEx = ExList.size();
|
||||
|
||||
// optional visualization control: pass "vis" as 3rd arg
|
||||
bool enableVis = (argc >= 3 && strcmp(argv[2], "vis") == 0);
|
||||
TApplication *app = nullptr;
|
||||
if(enableVis){
|
||||
app = new TApplication("anasenVis", &argc, argv);
|
||||
TEveManager::Create();
|
||||
TEvePointSet* pts = new TEvePointSet("hits");
|
||||
pts->SetMarkerStyle(20);
|
||||
pts->SetMarkerSize(1.4);
|
||||
pts->SetMarkerColor(kRed);
|
||||
gEve->AddElement(pts);
|
||||
SetVisPointSet(pts);
|
||||
}
|
||||
|
||||
// create detector representation in memory
|
||||
ANASEN * anasen = new ANASEN(); // top-level detector object
|
||||
SX3 * sx3 = anasen->GetSX3(); // silicon array part
|
||||
|
|
@ -86,6 +106,15 @@ int main(int argc, char **argv){
|
|||
printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
|
||||
TFile * saveFile = new TFile(saveFileName, "recreate");
|
||||
TTree * tree = new TTree("tree", "tree");
|
||||
TTree* tree2 = tree->CloneTree(0); // for 2D histograms or alternative data structure if needed
|
||||
TTree* visTree = nullptr;
|
||||
std::vector<double> visX, visY, visZ;
|
||||
if(enableVis){
|
||||
visTree = new TTree("visTree", "vis points");
|
||||
visTree->Branch("x", &visX);
|
||||
visTree->Branch("y", &visY);
|
||||
visTree->Branch("z", &visZ);
|
||||
}
|
||||
|
||||
// beam and CM variables saved in tree
|
||||
double KEA;
|
||||
|
|
@ -246,20 +275,32 @@ int main(int argc, char **argv){
|
|||
sx3Y = hitPos.Y();
|
||||
sx3Z = hitPos.Z();
|
||||
|
||||
// visualization point list
|
||||
if(enableVis) {
|
||||
visX.clear(); visY.clear(); visZ.clear();
|
||||
visX.push_back(sx3X);
|
||||
visY.push_back(sx3Y);
|
||||
visZ.push_back(sx3Z);
|
||||
PushEventAndRecord(visX, visY, visZ, visTree);
|
||||
}
|
||||
|
||||
// fill tree with original data before energy loss
|
||||
tree->Fill();
|
||||
|
||||
// apply energy loss from vertex to SX3 hit position (for light particle)
|
||||
//double dl = (hitPos - vertex).Mag() / 10.0; // path length in cm (positions in mm)
|
||||
//double EkinLight = Pb.E() - Pb.M();
|
||||
//double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx
|
||||
//double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling)
|
||||
//if (dE_light < EkinLight) {
|
||||
// Pb.SetE(Pb.E() - dE_light);
|
||||
// // update momentum to conserve direction
|
||||
// double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M());
|
||||
// TVector3 dir_new = Pb.Vect().Unit() * p_new;
|
||||
// Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E());
|
||||
//}
|
||||
double dl = (hitPos - vertex).Mag() / 10.0; // path length in cm (positions in mm)
|
||||
double EkinLight = Pb.E() - Pb.M();
|
||||
double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx
|
||||
double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling)
|
||||
if (dE_light < EkinLight) {
|
||||
Pb.SetE(Pb.E() - dE_light);
|
||||
// update momentum to conserve direction
|
||||
double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M());
|
||||
TVector3 dir_new = Pb.Vect().Unit() * p_new;
|
||||
Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E());
|
||||
}
|
||||
// update kinetic energy after loss
|
||||
//Tb = Pb.E() - Pb.M();
|
||||
Tb = Pb.E() - Pb.M();
|
||||
|
||||
// reconstruct track from PW readings + SX3 hit
|
||||
pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
|
||||
|
|
@ -273,6 +314,9 @@ int main(int argc, char **argv){
|
|||
|
||||
z0 = pw->GetZ0();
|
||||
|
||||
// fill tree2 with energy loss adjusted data
|
||||
tree2->Fill();
|
||||
|
||||
}else{
|
||||
// no valid SX3 hit: mark clearly invalid
|
||||
sx3Up = -1;
|
||||
|
|
@ -289,9 +333,10 @@ int main(int argc, char **argv){
|
|||
reTheta1 = TMath::QuietNaN();
|
||||
rePhi1 = TMath::QuietNaN();
|
||||
z0 = TMath::QuietNaN();
|
||||
}
|
||||
|
||||
// fill tree with original data (no energy loss for these events)
|
||||
tree->Fill();
|
||||
}
|
||||
|
||||
//#################################################################### Timer
|
||||
// measure elapsed real time and print progress roughly every 10 sec
|
||||
|
|
@ -314,15 +359,23 @@ int main(int argc, char **argv){
|
|||
|
||||
// write results to ROOT file and close
|
||||
tree->Write();
|
||||
tree2->Write();
|
||||
if(visTree) visTree->Write();
|
||||
int count = tree->GetEntries();
|
||||
int count2 = tree2->GetEntries();
|
||||
saveFile->Close();
|
||||
|
||||
printf("=============== done. saved as %s. count(hit==1) : %d\n", saveFileName.Data(), count);
|
||||
printf("=============== done. saved as %s. tree entries: %d, tree2 entries: %d\n", saveFileName.Data(), count, count2);
|
||||
|
||||
delete anasen;
|
||||
delete elossLight;
|
||||
delete elossHeavy;
|
||||
|
||||
if(enableVis && app){
|
||||
printf("Entering TEve GUI event loop (close window to finish)\n");
|
||||
app->Run();
|
||||
}
|
||||
|
||||
return 0;
|
||||
|
||||
}
|
||||
|
|
|
|||
Binary file not shown.
|
|
@ -1,309 +0,0 @@
|
|||
// add includes at top
|
||||
#include <TApplication.h>
|
||||
#include <TSystem.h>
|
||||
#include <TFile.h>
|
||||
#include <TTree.h>
|
||||
#include <TEveManager.h>
|
||||
#include <TEvePointSet.h>
|
||||
#include <TEveGeoNode.h>
|
||||
#include <TGeoManager.h>
|
||||
#include <vector>
|
||||
#include <mutex>
|
||||
|
||||
#include "TRandom.h"
|
||||
#include "TFile.h"
|
||||
#include "TTree.h"
|
||||
#include "TH1.h"
|
||||
#include "TH2.h"
|
||||
#include "TStyle.h"
|
||||
#include "TCanvas.h"
|
||||
#include "TBenchmark.h"
|
||||
|
||||
#include "ClassTransfer.h"
|
||||
#include "ClassAnasen.h"
|
||||
|
||||
// expose to ROOT
|
||||
int Run(int nEvents=1000, const char* outFile=nullptr){
|
||||
// Ensure TEve exists (create after geometry has been built if possible)
|
||||
if(!gEve) TEveManager::Create();
|
||||
|
||||
// if a geometry has already been loaded by ANASEN_model.C, make sure it
|
||||
// shows up in the TEve scene. TEveManager::Create() normally pulls in
|
||||
// gGeoManager, but we do it explicitly to be safe. We must wrap the
|
||||
// top node/volume in a TEveGeoTopNode (not pass a raw TGeoVolume).
|
||||
if(gGeoManager){
|
||||
// create a TEve wrapper around the top node
|
||||
TEveGeoTopNode *top = new TEveGeoTopNode(gGeoManager, gGeoManager->GetTopNode());
|
||||
gEve->AddElement(top);
|
||||
}
|
||||
|
||||
// Reaction
|
||||
TransferReaction transfer;
|
||||
|
||||
transfer.SetA(24,12, 0);
|
||||
transfer.SetIncidentEnergyAngle(10, 0, 0);
|
||||
transfer.Seta( 4, 2);
|
||||
transfer.Setb( 1, 1);
|
||||
|
||||
//TODO add alpha source
|
||||
|
||||
std::vector<float> ExAList = {0};
|
||||
std::vector<float> ExList = {0, 1, 2};
|
||||
|
||||
double vertexXRange[2] = { -5, 5}; // mm
|
||||
double vertexYRange[2] = { -5, 5};
|
||||
double vertexZRange[2] = { -100, 100};
|
||||
|
||||
double sigmaSX3_W = -1; // mm, < 0 use mid-point
|
||||
double sigmaSX3_L = 3; // mm, < 0 use mid-point
|
||||
double sigmaPW_A = 0; // from 0 to 1.
|
||||
double sigmaPW_C = 0; // from 0 to 1.
|
||||
|
||||
//###################################################
|
||||
|
||||
printf("------------ Vertex :\n");
|
||||
printf("X : %7.2f - %7.2f mm\n", vertexXRange[0], vertexXRange[1]);
|
||||
printf("Y : %7.2f - %7.2f mm\n", vertexYRange[0], vertexYRange[1]);
|
||||
printf("Z : %7.2f - %7.2f mm\n", vertexZRange[0], vertexZRange[1]);
|
||||
printf("------------ Uncertainty :\n");
|
||||
printf(" SX3 horizontal : %.1f\n", sigmaSX3_W);
|
||||
printf(" SX3 vertical : %.1f\n", sigmaSX3_L);
|
||||
printf(" Anode : %.1f mm\n", sigmaPW_A);
|
||||
printf(" Cathode : %.1f mm\n", sigmaPW_C);
|
||||
printf(" num_eve : %d \n",nEvents);
|
||||
transfer.CalReactionConstant();
|
||||
|
||||
int nExA = ExAList.size();
|
||||
int nEx = ExList.size();
|
||||
|
||||
ANASEN * anasen = new ANASEN();
|
||||
SX3 * sx3 = anasen->GetSX3();
|
||||
PW * pw = anasen->GetPW();
|
||||
|
||||
TString saveFileName = "SimAnasen1.root";
|
||||
printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
|
||||
TFile * saveFile = new TFile(saveFileName, "recreate");
|
||||
TTree * tree = new TTree("tree", "tree");
|
||||
|
||||
double KEA;
|
||||
tree->Branch("beamKEA", &KEA, "beamKEA/D");
|
||||
|
||||
double thetaCM, phiCM;
|
||||
tree->Branch("thetaCM", &thetaCM, "thetaCM/D");
|
||||
tree->Branch("phiCM", &phiCM, "phiCM/D");
|
||||
|
||||
double thetab, phib, Tb;
|
||||
double thetaB, phiB, TB;
|
||||
tree->Branch("thetab", &thetab, "thetab/D");
|
||||
tree->Branch("phib", &phib, "phib/D");
|
||||
tree->Branch("Tb", &Tb, "Tb/D");
|
||||
tree->Branch("thetaB", &thetaB, "thetaB/D");
|
||||
tree->Branch("phiB", &phiB, "phiB/D");
|
||||
tree->Branch("TB", &TB, "TB/D");
|
||||
|
||||
int ExAID;
|
||||
double ExA;
|
||||
tree->Branch("ExAID", &ExAID, "ExAID/I");
|
||||
tree->Branch("ExA", &ExA, "ExA/D");
|
||||
|
||||
int ExID;
|
||||
double Ex;
|
||||
tree->Branch("ExID", &ExID, "ExID/I");
|
||||
tree->Branch("Ex", &Ex, "Ex/D");
|
||||
|
||||
double vertexX, vertexY, vertexZ;
|
||||
tree->Branch("vX", &vertexX, "VertexX/D");
|
||||
tree->Branch("vY", &vertexY, "VertexY/D");
|
||||
tree->Branch("vZ", &vertexZ, "VertexZ/D");
|
||||
|
||||
double sx3X, sx3Y, sx3Z;
|
||||
tree->Branch("sx3X", &sx3X, "sx3X/D");
|
||||
tree->Branch("sx3Y", &sx3Y, "sx3Y/D");
|
||||
tree->Branch("sx3Z", &sx3Z, "sx3Z/D");
|
||||
|
||||
int anodeID[2], cathodeID[2];
|
||||
tree->Branch("aID", anodeID, "anodeID/I");
|
||||
tree->Branch("cID", cathodeID, "cathodeID/I");
|
||||
|
||||
double anodeDist[2], cathodeDist[2];
|
||||
tree->Branch("aDist", anodeDist, "anodeDist/D");
|
||||
tree->Branch("cDist", cathodeDist, "cathodeDist/D");
|
||||
|
||||
int sx3ID, sx3Up, sx3Dn, sx3Bk;
|
||||
double sx3ZFrac;
|
||||
tree->Branch("sx3ID", &sx3ID, "sx3ID/I");
|
||||
tree->Branch("sx3Up", &sx3Up, "sx3Up/I");
|
||||
tree->Branch("sx3Dn", &sx3Dn, "sx3Dn/I");
|
||||
tree->Branch("sx3Bk", &sx3Bk, "sx3Bk/I");
|
||||
tree->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D");
|
||||
|
||||
double reTheta, rePhi;
|
||||
tree->Branch("reTheta", &reTheta, "reconstucted_theta/D");
|
||||
tree->Branch("rePhi", &rePhi, "reconstucted_phi/D");
|
||||
|
||||
double reTheta1, rePhi1;
|
||||
tree->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D");
|
||||
tree->Branch("rePhi1", &rePhi1, "reconstucted_phi1/D");
|
||||
|
||||
double z0;
|
||||
tree->Branch("z0", &z0, "reconstucted_Z/D");
|
||||
|
||||
|
||||
//========timer
|
||||
TBenchmark clock;
|
||||
bool shown ;
|
||||
clock.Reset();
|
||||
clock.Start("timer");
|
||||
shown = false;
|
||||
|
||||
// Create a point set to show hits
|
||||
TEvePointSet *pts = new TEvePointSet("hits");
|
||||
pts->SetMarkerStyle(20);
|
||||
pts->SetMarkerColor(kRed);
|
||||
gEve->AddElement(pts);
|
||||
|
||||
// Optionally open output file/tree
|
||||
TFile *fout = nullptr;
|
||||
TTree *tout = nullptr;
|
||||
std::vector<double> vx, vy, vz;
|
||||
if(outFile){
|
||||
fout = TFile::Open(outFile,"RECREATE");
|
||||
tout = new TTree("evt","events");
|
||||
tout->Branch("x",&vx);
|
||||
tout->Branch("y",&vy);
|
||||
tout->Branch("z",&vz);
|
||||
}
|
||||
|
||||
// Simulation loop (replace with your sim code that fills vx,vy,vz per event)
|
||||
|
||||
for( int i = 0; i < nEvents ; i++){
|
||||
|
||||
ExAID = gRandom->Integer(nExA);
|
||||
ExA = ExAList[ExAID];
|
||||
transfer.SetExA(ExA);
|
||||
|
||||
ExID = gRandom->Integer(nEx);
|
||||
Ex = ExList[ExID];
|
||||
transfer.SetExB(Ex);
|
||||
|
||||
transfer.CalReactionConstant();
|
||||
|
||||
thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ;
|
||||
phiCM = (gRandom->Rndm() - 0.5) * TMath::TwoPi();
|
||||
|
||||
//==== Calculate reaction
|
||||
TLorentzVector * output = transfer.Event(thetaCM, phiCM);
|
||||
TLorentzVector Pb = output[2];
|
||||
TLorentzVector PB = output[3];
|
||||
|
||||
thetab = Pb.Theta() * TMath::RadToDeg();
|
||||
thetaB = PB.Theta() * TMath::RadToDeg();
|
||||
|
||||
Tb = Pb.E() - Pb.M();
|
||||
TB = PB.E() - PB.M();
|
||||
|
||||
phib = Pb.Phi() * TMath::RadToDeg();
|
||||
phiB = PB.Phi() * TMath::RadToDeg();
|
||||
|
||||
vertexX = (vertexXRange[1]- vertexXRange[0])*gRandom->Rndm() + vertexXRange[0];
|
||||
vertexY = (vertexYRange[1]- vertexYRange[0])*gRandom->Rndm() + vertexYRange[0];
|
||||
vertexZ = (vertexZRange[1]- vertexZRange[0])*gRandom->Rndm() + vertexZRange[0];
|
||||
|
||||
TVector3 vertex(vertexX, vertexY, vertexZ);
|
||||
|
||||
TVector3 dir(1, 0, 0);
|
||||
dir.SetTheta(thetab * TMath::DegToRad());
|
||||
dir.SetPhi(phib * TMath::DegToRad());
|
||||
|
||||
|
||||
pw->FindWireID(vertex, dir, false);
|
||||
sx3->FindSX3Pos(vertex, dir, false);
|
||||
|
||||
PWHitInfo hitInfo = pw->GetHitInfo();
|
||||
|
||||
anodeID[0] = hitInfo.nearestWire.first;
|
||||
cathodeID[0] = hitInfo.nearestWire.second;
|
||||
anodeID[1] = hitInfo.nextNearestWire.first;
|
||||
cathodeID[1] = hitInfo.nextNearestWire.second;
|
||||
|
||||
anodeDist[0] = hitInfo.nearestDist.first;
|
||||
cathodeDist[0] = hitInfo.nearestDist.second;
|
||||
anodeDist[1] = hitInfo.nextNearestDist.first;
|
||||
cathodeDist[1] = hitInfo.nextNearestDist.second;
|
||||
|
||||
sx3ID = sx3->GetID();
|
||||
if( sx3ID >= 0 ){
|
||||
sx3Up = sx3->GetChUp();
|
||||
sx3Dn = sx3->GetChDn();
|
||||
sx3Bk = sx3->GetChBk();
|
||||
sx3ZFrac = sx3->GetZFrac();
|
||||
|
||||
//Introduce uncertaity
|
||||
// TVector3 hitPos = sx3->GetHitPos();
|
||||
TVector3 hitPos = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
|
||||
|
||||
sx3X = hitPos.X();
|
||||
sx3Y = hitPos.Y();
|
||||
sx3Z = hitPos.Z();
|
||||
|
||||
pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
|
||||
reTheta = pw->GetTrackTheta() * TMath::RadToDeg();
|
||||
rePhi = pw->GetTrackPhi() * TMath::RadToDeg();
|
||||
|
||||
pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false);
|
||||
reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg();
|
||||
rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg();
|
||||
|
||||
z0 = pw->GetZ0();
|
||||
|
||||
}else{
|
||||
sx3Up = -1;
|
||||
sx3Dn = -1;
|
||||
sx3Bk = -1;
|
||||
sx3ZFrac = TMath::QuietNaN();
|
||||
|
||||
sx3X = TMath::QuietNaN();
|
||||
sx3Y = TMath::QuietNaN();
|
||||
sx3Z = TMath::QuietNaN();
|
||||
|
||||
// for( int i = 0; i < 12; i++){
|
||||
// sx3Index[i] = -1;
|
||||
// }
|
||||
|
||||
reTheta = TMath::QuietNaN();
|
||||
rePhi = TMath::QuietNaN();
|
||||
|
||||
reTheta1 = TMath::QuietNaN();
|
||||
rePhi1 = TMath::QuietNaN();
|
||||
|
||||
z0 = TMath::QuietNaN();
|
||||
|
||||
}
|
||||
// -----------------------------------------------------------------------
|
||||
|
||||
|
||||
}
|
||||
|
||||
// update TEve
|
||||
pts->Reset();
|
||||
for(size_t i=0;i<vx.size(); ++i) pts->SetNextPoint(vx[i], vy[i], vz[i]);
|
||||
gEve->Redraw3D();
|
||||
gSystem->ProcessEvents();
|
||||
|
||||
// write to tree
|
||||
if(tout){ tout->Fill(); fout->Flush(); }
|
||||
|
||||
if(fout) fout->Close();
|
||||
return 0;
|
||||
}
|
||||
|
||||
// optional main to keep standalone build working
|
||||
#ifndef __CLING__
|
||||
int main(int argc, char** argv){
|
||||
TApplication app("app",&argc,argv);
|
||||
// if you want to import geometry here when running standalone:
|
||||
// TGeoManager::Import("yourGeom.root");
|
||||
Run(500, "sim_out.root");
|
||||
return 0;
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1,6 +1,8 @@
|
|||
// vis_helpers.h (or paste into anasenMS.cpp)
|
||||
#ifndef VIS_HELPERS_H
|
||||
#define VIS_HELPERS_H
|
||||
|
||||
#include <TSystem.h>
|
||||
#include <TEvePointSet.h>
|
||||
#include <TFile.h>
|
||||
#include <TTree.h>
|
||||
#include <vector>
|
||||
#include <mutex>
|
||||
|
|
@ -8,27 +10,52 @@
|
|||
static TEvePointSet* gVisPts = nullptr;
|
||||
static std::mutex gVisMutex;
|
||||
|
||||
// Call from your ROOT session after creating TEve objects:
|
||||
void SetVisPointSet(TEvePointSet* pts){ gVisPts = pts; }
|
||||
// Recommended: call once after opening TEve and adding a point set to gEve
|
||||
inline void SetVisPointSet(TEvePointSet* pts) { gVisPts = pts; }
|
||||
|
||||
// Call this from your sim loop to update visualization and optionally write data:
|
||||
void PushEventAndRecord(const std::vector<double>& x,
|
||||
inline void UpdateVisPointSet(const std::vector<double>& x,
|
||||
const std::vector<double>& y,
|
||||
const std::vector<double>& z,
|
||||
TTree* outTree = nullptr)
|
||||
const std::vector<double>& z)
|
||||
{
|
||||
if(outTree){
|
||||
outTree->SetBranchAddress("x",(void*)&x);
|
||||
outTree->SetBranchAddress("y",(void*)&y);
|
||||
outTree->SetBranchAddress("z",(void*)&z);
|
||||
if(!gVisPts) return;
|
||||
std::lock_guard<std::mutex> lk(gVisMutex);
|
||||
gVisPts->Reset();
|
||||
size_t n = std::min({x.size(), y.size(), z.size()});
|
||||
for(size_t i=0; i<n; ++i) gVisPts->SetNextPoint(x[i], y[i], z[i]);
|
||||
if(gEve) {
|
||||
gEve->Redraw3D();
|
||||
gSystem->ProcessEvents();
|
||||
}
|
||||
}
|
||||
|
||||
// Fill a tree with pointlists (one entry per event); must have branches defined once by caller
|
||||
inline void RecordTreeXYZ(TTree* outTree,
|
||||
const std::vector<double>& x,
|
||||
const std::vector<double>& y,
|
||||
const std::vector<double>& z)
|
||||
{
|
||||
if(!outTree) return;
|
||||
static std::vector<double> tx, ty, tz;
|
||||
tx = x;
|
||||
ty = y;
|
||||
tz = z;
|
||||
|
||||
if(outTree->GetBranch("x") == nullptr) outTree->Branch("x", &tx);
|
||||
if(outTree->GetBranch("y") == nullptr) outTree->Branch("y", &ty);
|
||||
if(outTree->GetBranch("z") == nullptr) outTree->Branch("z", &tz);
|
||||
|
||||
// Do NOT call SetBranchAddress() for the branch we are filling.
|
||||
outTree->Fill();
|
||||
outTree->GetCurrentFile()->Flush();
|
||||
}
|
||||
|
||||
if(!gVisPts) return;
|
||||
std::lock_guard<std::mutex> lk(gVisMutex);
|
||||
gVisPts->Reset();
|
||||
for(size_t i=0;i<x.size(); ++i) gVisPts->SetNextPoint(x[i], y[i], z[i]);
|
||||
gEve->Redraw3D();
|
||||
gSystem->ProcessEvents();
|
||||
inline void PushEventAndRecord(const std::vector<double>& x,
|
||||
const std::vector<double>& y,
|
||||
const std::vector<double>& z,
|
||||
TTree* outTree = nullptr)
|
||||
{
|
||||
if(outTree) RecordTreeXYZ(outTree, x, y, z);
|
||||
UpdateVisPointSet(x,y,z);
|
||||
}
|
||||
|
||||
#endif // VIS_HELPERS_H
|
||||
|
|
|
|||
Loading…
Reference in New Issue
Block a user