diff --git a/.gitignore b/.gitignore index 193cfae..bc47dfe 100644 --- a/.gitignore +++ b/.gitignore @@ -8,6 +8,7 @@ EventBuilder* *.err *.seq *.png +*.pdf Mapper AnasenMS Armory/anasenMS diff --git a/.vscode/c_cpp_properties.json b/.vscode/c_cpp_properties.json index 744c2d6..1f22739 100644 --- a/.vscode/c_cpp_properties.json +++ b/.vscode/c_cpp_properties.json @@ -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": [ diff --git a/Armory/Makefile b/Armory/Makefile index 14092ca..31a4ecf 100644 --- a/Armory/Makefile +++ b/Armory/Makefile @@ -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) diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index fa3d31b..c3e7129 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -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 +#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 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(); - } - tree->Fill(); + // 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; } diff --git a/Armory/anasenMS_root b/Armory/anasenMS_root deleted file mode 100755 index cdb335e..0000000 Binary files a/Armory/anasenMS_root and /dev/null differ diff --git a/Armory/anasenMS_root.cpp b/Armory/anasenMS_root.cpp deleted file mode 100644 index 5a8ad39..0000000 --- a/Armory/anasenMS_root.cpp +++ /dev/null @@ -1,309 +0,0 @@ -// add includes at top -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#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 ExAList = {0}; - std::vector 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 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;iSetNextPoint(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 diff --git a/Armory/vis_helpers.h b/Armory/vis_helpers.h index b6c01c8..af16edc 100644 --- a/Armory/vis_helpers.h +++ b/Armory/vis_helpers.h @@ -1,6 +1,8 @@ -// vis_helpers.h (or paste into anasenMS.cpp) +#ifndef VIS_HELPERS_H +#define VIS_HELPERS_H + +#include #include -#include #include #include #include @@ -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& x, - const std::vector& y, - const std::vector& z, - TTree* outTree = nullptr) +inline void UpdateVisPointSet(const std::vector& x, + const std::vector& y, + const std::vector& z) { - if(outTree){ - outTree->SetBranchAddress("x",(void*)&x); - outTree->SetBranchAddress("y",(void*)&y); - outTree->SetBranchAddress("z",(void*)&z); - outTree->Fill(); - outTree->GetCurrentFile()->Flush(); - } - if(!gVisPts) return; std::lock_guard lk(gVisMutex); gVisPts->Reset(); - for(size_t i=0;iSetNextPoint(x[i], y[i], z[i]); - gEve->Redraw3D(); - gSystem->ProcessEvents(); + size_t n = std::min({x.size(), y.size(), z.size()}); + for(size_t i=0; iSetNextPoint(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& x, + const std::vector& y, + const std::vector& z) +{ + if(!outTree) return; + static std::vector 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(); +} + +inline void PushEventAndRecord(const std::vector& x, + const std::vector& y, + const std::vector& z, + TTree* outTree = nullptr) +{ + if(outTree) RecordTreeXYZ(outTree, x, y, z); + UpdateVisPointSet(x,y,z); +} + +#endif // VIS_HELPERS_H