rtang/ANASEN_analysis
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

redo

Browse Source
This commit is contained in:
James Szalkie 2026-03-30 10:49:41 -04:00
parent e3ffe56351
commit a0f0a7083d
6 changed files with 305 additions and 62 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
.gitignore vendored
View File

@ -10,6 +10,7 @@ EventBuilder*
*.png
Mapper
AnasenMS
Armory/anasenMS
data/
data_proton/

9
Armory/Makefile
View File

@ -8,7 +8,11 @@ CC = g++
#COPTS = -fPIC -DLINUX -O2 -std=c++17 -lpthread
COPTS = -fPIC -DLINUX -g -O0 -Wall -std=c++17 -lpthread
ROOTLIBS = `root-config --cflags --glibs`
ROOTCFLAGS := $(shell root-config --cflags)
ROOTLIBS := $(shell root-config --libs)
CXX := g++
CXXFLAGS := -O2 $(ROOTCFLAGS)
ALL = Mapper EventBuilder#AnasenMS
@ -31,3 +35,6 @@ 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)

63
Armory/README_anasenMS.md Normal file
View File

@ -0,0 +1,63 @@
# ANASEN Monte Carlo (anasenMS.cpp)
## Overview
`anasenMS.cpp` is a standalone Monte Carlo simulation for an ANASEN-style detector setup. It generates transfer reaction kinematics, propagates products to a wire chamber (PW) and a silicon array (SX3), reconstructs tracks, and writes output to a ROOT tree.
## Requirements
- ROOT (e.g. `root-config` for compile flags)
- C++ compiler (gcc/g++)
- Project includes: `ClassTransfer.h`, `ClassAnasen.h`, plus their dependent implementation files.
## Build
In `Armory` directory:
```bash
g++ -O2 -o anasenMS anasenMS.cpp ClassTransfer.cpp ClassAnasen.cpp ... `root-config --cflags --libs`
```
(Adjust source file list based on actual project layout.)
## Run
```bash
./anasenMS [numEvents]
```
- `numEvents`: optional integer, default `1000000`
- Outputs: `SimAnasen1.root` containing `tree`
## What the code does
- Initializes reaction: `TransferReaction transfer`
- `SetA(24,12,0)` target
- `SetIncidentEnergyAngle(10,0,0)` beam energy and direction
- `Seta`, `Setb` reaction fragment indices
- Sets excitation lists: `ExAList`, `ExList`
- Vertex and resolution settings:
- `vertexX/Y/Z` ranges
- `sigmaSX3_W`, `sigmaSX3_L`, `sigmaPW_A`, `sigmaPW_C`
- Prepares ROOT output tree and branches for truth/reconstructed
- Loop over events:
- Sample excitation and CM direction
- `transfer.Event(thetaCM, phiCM)` outputs `TLorentzVector` products
- Compute lab angles/energies
- Random vertex inside target volume
- Run detector response:
- `pw->FindWireID(...)`
- `sx3->FindSX3Pos(...)`
- Read out wire hits and SX3 channel + depth
- Apply position smearing for SX3
- Reconstruct track via `pw->CalTrack` and `pw->CalTrack2`
- Fill ROOT tree
- At end: write tree, close file, clean up
## Notes
- The code now has expanded inline comments explaining each stage.
- Important methods are from:
- `ClassTransfer` (`SetA`, `SetIncidentEnergyAngle`, `Seta`, `Setb`, `SetExA`, `SetExB`, `CalReactionConstant`, `Event`)
- `ClassAnasen` / `SX3` / `PW` (`FindWireID`, `FindSX3Pos`, `GetHitInfo`, `CalTrack`, `CalTrack2`, `GetTrackTheta`, `GetTrackPhi`, `GetZ0`, `GetHitPosWithSigma`, `GetID`, etc.)
- Optional: change excitation lists, vertex spread, and sigma values to mimic different beam/target conditions.

107
Armory/anasenMS.cpp
View File

@ -1,18 +1,17 @@
#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 "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 "TH2.h" // 2D histograms
#include "TStyle.h" // ROOT plotting style controls
#include "TCanvas.h" // ROOT canvas drawing
#include "TBenchmark.h" // timing measurement
#include "ClassTransfer.h"
#include "ClassAnasen.h"
#include "ClassTransfer.h" // Reaction kinematics and MC event generation
#include "ClassAnasen.h" // ANASEN detector model classes (SX3, PW, etc.)
//======== Gerneate light particle based on reaction
// find out the CalTrack and the real track
// find out the Q-value uncertaintly
//======== Generate light particle based on reaction
// calculate real and reconstructed tracks and Q-value uncertainty
int main(int argc, char **argv){
@ -20,33 +19,36 @@ int main(int argc, char **argv){
printf("=== ANASEN Monte Carlo ===\n");
printf("=========================================\n");
// number of events can be overridden from command line
int numEvent = 1000000;
if( argc >= 2 ) numEvent = atoi(argv[1]);
//Reaction
// Reaction setup: projectile + target configuration, energy, and product IDs
TransferReaction transfer;
transfer.SetA(24,12, 0);
transfer.SetIncidentEnergyAngle(10, 0, 0);
transfer.Seta( 4, 2);
transfer.Setb( 1, 1);
transfer.SetA(24,12, 0); // e.g., 24Mg (Z=12) with 0 excitation
transfer.SetIncidentEnergyAngle(10, 0, 0); // 10 MeV beam, 0 polar and azimuthal angle
transfer.Seta( 4, 2); // identify reaction product a in internal indexing
transfer.Setb( 1, 1); // identify reaction product b
//TODO add alpha source
// TODO add alpha source or alternative reaction channel selection
// Excited state lists (target and projectile/excited products)
std::vector<float> ExAList = {0};
std::vector<float> ExList = {0, 1, 2};
// define vertex position uniform distribution ranges (mm)
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.
//###################################################
// detector resolution / uncertainty parameters
double sigmaSX3_W = -1; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension)
double sigmaSX3_L = 3; // mm, vertical spread for SX3
double sigmaPW_A = 0; // normalized anode uncertainty term (0-1)
double sigmaPW_C = 0; // normalized cathode uncertainty term (0-1)
// status printout
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]);
@ -57,20 +59,25 @@ int main(int argc, char **argv){
printf(" Anode : %.1f mm\n", sigmaPW_A);
printf(" Cathode : %.1f mm\n", sigmaPW_C);
printf(" num_eve : %d \n",numEvent);
// calculates energy/momentum/kinematics constants for transfer reaction
transfer.CalReactionConstant();
int nExA = ExAList.size();
int nEx = ExList.size();
ANASEN * anasen = new ANASEN();
SX3 * sx3 = anasen->GetSX3();
PW * pw = anasen->GetPW();
// create detector representation in memory
ANASEN * anasen = new ANASEN(); // top-level detector object
SX3 * sx3 = anasen->GetSX3(); // silicon array part
PW * pw = anasen->GetPW(); // proportional wire chamber part
// output file + tree
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");
// beam and CM variables saved in tree
double KEA;
tree->Branch("beamKEA", &KEA, "beamKEA/D");
@ -78,6 +85,7 @@ int main(int argc, char **argv){
tree->Branch("thetaCM", &thetaCM, "thetaCM/D");
tree->Branch("phiCM", &phiCM, "phiCM/D");
// outgoing particles in lab frame (light/heavy)
double thetab, phib, Tb;
double thetaB, phiB, TB;
tree->Branch("thetab", &thetab, "thetab/D");
@ -87,6 +95,7 @@ int main(int argc, char **argv){
tree->Branch("phiB", &phiB, "phiB/D");
tree->Branch("TB", &TB, "TB/D");
// excitation state identifiers
int ExAID;
double ExA;
tree->Branch("ExAID", &ExAID, "ExAID/I");
@ -97,24 +106,29 @@ int main(int argc, char **argv){
tree->Branch("ExID", &ExID, "ExID/I");
tree->Branch("Ex", &Ex, "Ex/D");
// true vertex position in target volume
double vertexX, vertexY, vertexZ;
tree->Branch("vX", &vertexX, "VertexX/D");
tree->Branch("vY", &vertexY, "VertexY/D");
tree->Branch("vZ", &vertexZ, "VertexZ/D");
// reconstructed SX3 hit position
double sx3X, sx3Y, sx3Z;
tree->Branch("sx3X", &sx3X, "sx3X/D");
tree->Branch("sx3Y", &sx3Y, "sx3Y/D");
tree->Branch("sx3Z", &sx3Z, "sx3Z/D");
// PW nearest and next nearest wires
int anodeID[2], cathodeID[2];
tree->Branch("aID", anodeID, "anodeID/I");
tree->Branch("cID", cathodeID, "cathodeID/I");
// distances to nearest wires
double anodeDist[2], cathodeDist[2];
tree->Branch("aDist", anodeDist, "anodeDist/D");
tree->Branch("cDist", cathodeDist, "cathodeDist/D");
// SX3 channel assignment and Z fraction (depth) information
int sx3ID, sx3Up, sx3Dn, sx3Bk;
double sx3ZFrac;
tree->Branch("sx3ID", &sx3ID, "sx3ID/I");
@ -123,6 +137,7 @@ int main(int argc, char **argv){
tree->Branch("sx3Bk", &sx3Bk, "sx3Bk/I");
tree->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D");
// reconstructed angles from PW track fit, method 1 and 2
double reTheta, rePhi;
tree->Branch("reTheta", &reTheta, "reconstucted_theta/D");
tree->Branch("rePhi", &rePhi, "reconstucted_phi/D");
@ -131,6 +146,7 @@ int main(int argc, char **argv){
tree->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D");
tree->Branch("rePhi1", &rePhi1, "reconstucted_phi1/D");
// reconstructed vertex Z from PW fit
double z0;
tree->Branch("z0", &z0, "reconstucted_Z/D");
@ -142,9 +158,10 @@ int main(int argc, char **argv){
clock.Start("timer");
shown = false;
//================================= Calculate event
//================================= Calculate event loop
for( int i = 0; i < numEvent ; i++){
// randomly sample target/projectile excitations
ExAID = gRandom->Integer(nExA);
ExA = ExAList[ExAID];
transfer.SetExA(ExA);
@ -153,15 +170,17 @@ int main(int argc, char **argv){
Ex = ExList[ExID];
transfer.SetExB(Ex);
// recalc kinematic constants for chosen states
transfer.CalReactionConstant();
// isotropic CM direction
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];
//==== Calculate reaction kinematics in lab frame
TLorentzVector * output = transfer.Event(thetaCM, phiCM); // returns array of outputs
TLorentzVector Pb = output[2]; // light particle or product A
TLorentzVector PB = output[3]; // heavy particle or product B
thetab = Pb.Theta() * TMath::RadToDeg();
thetaB = PB.Theta() * TMath::RadToDeg();
@ -172,22 +191,25 @@ int main(int argc, char **argv){
phib = Pb.Phi() * TMath::RadToDeg();
phiB = PB.Phi() * TMath::RadToDeg();
// vertex position in target volume
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);
// set direction vector from lab angle
TVector3 dir(1, 0, 0);
dir.SetTheta(thetab * TMath::DegToRad());
dir.SetPhi(phib * TMath::DegToRad());
// run detector response models for PW and SX3
pw->FindWireID(vertex, dir, false);
sx3->FindSX3Pos(vertex, dir, false);
PWHitInfo hitInfo = pw->GetHitInfo();
// store nearest/next closest wire data
anodeID[0] = hitInfo.nearestWire.first;
cathodeID[0] = hitInfo.nearestWire.second;
anodeID[1] = hitInfo.nextNearestWire.first;
@ -198,6 +220,7 @@ int main(int argc, char **argv){
anodeDist[1] = hitInfo.nextNearestDist.first;
cathodeDist[1] = hitInfo.nextNearestDist.second;
// SX3 hit channel info and depth fraction
sx3ID = sx3->GetID();
if( sx3ID >= 0 ){
sx3Up = sx3->GetChUp();
@ -205,18 +228,20 @@ int main(int argc, char **argv){
sx3Bk = sx3->GetChBk();
sx3ZFrac = sx3->GetZFrac();
//Introduce uncertaity
// TVector3 hitPos = sx3->GetHitPos();
// apply intrinsic detector resolution to true SX3 hit position
// for no smearing comment out and use GetHitPos();
TVector3 hitPos = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
sx3X = hitPos.X();
sx3Y = hitPos.Y();
sx3Z = hitPos.Z();
// reconstruct track from PW readings + SX3 hit
pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);
reTheta = pw->GetTrackTheta() * TMath::RadToDeg();
rePhi = pw->GetTrackPhi() * TMath::RadToDeg();
// alternative track algorithm with uncertainty parameters
pw->CalTrack2(hitPos, hitInfo, sigmaPW_A, sigmaPW_C, false);
reTheta1 = pw->GetTrackTheta() * TMath::RadToDeg();
rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg();
@ -224,6 +249,7 @@ int main(int argc, char **argv){
z0 = pw->GetZ0();
}else{
// no valid SX3 hit: mark clearly invalid
sx3Up = -1;
sx3Dn = -1;
sx3Bk = -1;
@ -233,23 +259,17 @@ int main(int argc, char **argv){
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();
}
tree->Fill();
//#################################################################### Timer
// measure elapsed real time and print progress roughly every 10 sec
clock.Stop("timer");
Double_t time = clock.GetRealTime("timer");
clock.Start("timer");
@ -267,6 +287,7 @@ int main(int argc, char **argv){
}
// write results to ROOT file and close
tree->Write();
int count = tree->GetEntries();
saveFile->Close();

131
DetectorConstruction.cc Normal file
View File

@ -0,0 +1,131 @@
#include "DetectorConstruction.hh"
#include "G4RunManager.hh"
#include "G4NistManager.hh"
#include "G4Box.hh"
#include "G4Tubs.hh"
#include "G4LogicalVolume.hh"
#include "G4PVPlacement.hh"
#include "G4SystemOfUnits.hh"
#include "G4RotationMatrix.hh"
#include "G4ThreeVector.hh"
DetectorConstruction::DetectorConstruction()
: G4VUserDetectorConstruction(),
fScoringVolume(nullptr)
{}
DetectorConstruction::~DetectorConstruction()
{}
G4VPhysicalVolume* DetectorConstruction::Construct()
{
// Get nist material manager
G4NistManager* nist = G4NistManager::Instance();
// Materials
G4Material* vacuum = nist->FindOrBuildMaterial("G4_Galactic");
G4Material* al = nist->FindOrBuildMaterial("G4_Al");
// World
G4double worldx = 200.*mm;
G4double worldy = 200.*mm;
G4double worldz = 200.*mm;
G4Box* solidWorld = new G4Box("World", worldx, worldy, worldz);
G4LogicalVolume* logicWorld = new G4LogicalVolume(solidWorld, vacuum, "World");
G4VPhysicalVolume* physWorld = new G4PVPlacement(0, G4ThreeVector(), logicWorld, "World", 0, false, 0, true);
// Axes (optional, for visualization)
G4Tubs* solidAxis = new G4Tubs("Axis", 0, 0.1*mm, 5.*mm, 0, 360*deg);
G4LogicalVolume* logicAxisX = new G4LogicalVolume(solidAxis, al, "AxisX");
G4LogicalVolume* logicAxisY = new G4LogicalVolume(solidAxis, al, "AxisY");
G4LogicalVolume* logicAxisZ = new G4LogicalVolume(solidAxis, al, "AxisZ");
G4RotationMatrix* rotX = new G4RotationMatrix();
rotX->rotateY(90*deg);
rotX->rotateZ(90*deg);
new G4PVPlacement(rotX, G4ThreeVector(5*mm, 0, 0), logicAxisX, "AxisX", logicWorld, false, 0);
G4RotationMatrix* rotY = new G4RotationMatrix();
rotY->rotateX(90*deg);
new G4PVPlacement(rotY, G4ThreeVector(0, 5*mm, 0), logicAxisY, "AxisY", logicWorld, false, 0);
new G4PVPlacement(0, G4ThreeVector(0, 0, 5*mm), logicAxisZ, "AxisZ", logicWorld, false, 0);
// ANASEN geometry
const int nWire = 24;
const int wireShift = 3;
const G4double zLen = 350.*mm;
const G4double radiusA = 38.*mm;
const G4double radiusC = 43.*mm;
G4double dAngle = wireShift * 2 * CLHEP::pi / nWire;
G4double radiusAnew = radiusA * cos(dAngle / 2.);
G4double wireALength = sqrt(zLen*zLen + pow(2 * radiusA * sin(dAngle/2.), 2));
G4double wireATheta = atan2(2 * radiusA * sin(dAngle / 2.), zLen);
G4Tubs* solidPC_A = new G4Tubs("PC_A", 0, 0.01*mm, wireALength/2., 0, 360*deg);
G4LogicalVolume* logicPC_A = new G4LogicalVolume(solidPC_A, al, "PC_A");
for(int i = 0; i < nWire; i++){
G4double phi = 2 * CLHEP::pi / nWire * i + dAngle / 2.;
G4ThreeVector pos(radiusAnew * cos(phi), radiusAnew * sin(phi), 0);
G4RotationMatrix* rot = new G4RotationMatrix();
rot->rotateZ(360./nWire * (i + wireShift/2.) * deg);
rot->rotateY(wireATheta);
new G4PVPlacement(rot, pos, logicPC_A, "PC_A", logicWorld, false, i);
}
G4double radiusCnew = radiusC * cos(dAngle / 2.);
G4double wireCLength = sqrt(zLen*zLen + pow(2 * radiusC * sin(dAngle/2.), 2));
G4double wireCTheta = atan2(2 * radiusC * sin(dAngle / 2.), zLen);
G4Tubs* solidPC_C = new G4Tubs("PC_C", 0, 0.01*mm, wireCLength/2., 0, 360*deg);
G4LogicalVolume* logicPC_C = new G4LogicalVolume(solidPC_C, al, "PC_C");
for(int i = 0; i < nWire; i++){
G4double phi = 2 * CLHEP::pi / nWire * i - dAngle/2.;
G4ThreeVector pos(radiusCnew * cos(phi), radiusCnew * sin(phi), 0);
G4RotationMatrix* rot = new G4RotationMatrix();
rot->rotateZ(360./nWire * (i - wireShift/2.) * deg);
rot->rotateY(-wireCTheta);
new G4PVPlacement(rot, pos, logicPC_C, "PC_C", logicWorld, false, i);
}
const int nSX3 = 12;
const G4double sx3Radius = 88.*mm;
const G4double sx3Width = 40.*mm;
const G4double sx3Length = 75.*mm;
const G4double sx3Gap = 5.*mm;
G4Box* solidSX3 = new G4Box("SX3", 0.1*mm, sx3Width/2., sx3Length/2.);
G4LogicalVolume* logicSX3 = new G4LogicalVolume(solidSX3, al, "SX3");
fScoringVolume = logicSX3;
for(int i = 0; i < nSX3; i++){
G4double phi = 2 * CLHEP::pi / nSX3 * (i + 0.5);
G4ThreeVector pos1(sx3Radius * cos(phi), sx3Radius * sin(phi), sx3Length/2. + sx3Gap);
G4RotationMatrix* rot1 = new G4RotationMatrix();
rot1->rotateZ(360./nSX3 * (i + 0.5) * deg);
new G4PVPlacement(rot1, pos1, logicSX3, "SX3_front", logicWorld, false, 2*i);
G4ThreeVector pos2(sx3Radius * cos(phi), sx3Radius * sin(phi), -sx3Length/2. - sx3Gap);
G4RotationMatrix* rot2 = new G4RotationMatrix();
rot2->rotateZ(360./nSX3 * (i + 0.5) * deg);
new G4PVPlacement(rot2, pos2, logicSX3, "SX3_back", logicWorld, false, 2*i+1);
}
const G4double qqqR1 = 50.*mm;
const G4double qqqR2 = 100.*mm;
G4Tubs* solidQQQ = new G4Tubs("QQQ", qqqR1, qqqR2, 0.5*mm, 5*deg, 85*deg);
G4LogicalVolume* logicQQQ = new G4LogicalVolume(solidQQQ, al, "QQQ");
for(int i = 0; i < 4; i++){
G4ThreeVector pos(0, 0, 100.*mm);
G4RotationMatrix* rot = new G4RotationMatrix();
rot->rotateZ(360./4 * i * deg);
new G4PVPlacement(rot, pos, logicQQQ, "QQQ", logicWorld, false, i);
}
return physWorld;
}

20
DetectorConstruction.hh Normal file
View File

@ -0,0 +1,20 @@
#ifndef DetectorConstruction_h
#define DetectorConstruction_h 1
#include "G4VUserDetectorConstruction.hh"
class DetectorConstruction : public G4VUserDetectorConstruction
{
public:
DetectorConstruction();
virtual ~DetectorConstruction();
virtual G4VPhysicalVolume* Construct();
G4LogicalVolume* GetScoringVolume() const { return fScoringVolume; }
private:
G4LogicalVolume* fScoringVolume;
};
#endif