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

updated kinematic scrubbing, added histcomp

Browse Source
This commit is contained in:
James Szalkie 2026-04-09 13:03:14 -04:00
parent 1339218088
commit ab06f42956
4 changed files with 136 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
Armory/ClassAnasen.h
View File

@ -108,17 +108,21 @@ inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1,
geom = new TGeoManager("Detector", "ANASEN");
//--- define some materials
TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0);
TGeoMaterial *matVacuum = new TGeoMaterial("Vacuum", 0,0,0); //name, A, Z, density
TGeoMaterial *matHe = new TGeoMaterial("He", 4.0026, 2, 0.000861);
TGeoMaterial *matAl = new TGeoMaterial("Al", 26.98,13,2.7);
TGeoMaterial *matSi = new TGeoMaterial("Si", 28.085,14,2.33);
//--- define some media
TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum);
TGeoMedium *Vacuum = new TGeoMedium("Vacuum",1, matVacuum); //name, number of materials, material
TGeoMedium *He = new TGeoMedium("He",2, matHe);
TGeoMedium *Al = new TGeoMedium("Root Material",2, matAl);
TGeoMedium *Si = new TGeoMedium("Si",3, matSi);
//--- make the top container volume
Double_t worldx = 200.; //mm
Double_t worldy = 200.; //mm
Double_t worldz = 200.; //mm
worldBox = geom->MakeBox("ROOT", Vacuum, worldx, worldy, worldz);
worldBox = geom->MakeBox("ROOT", He, worldx, worldy, worldz); // name, medium, x half-length, y half-length, z half-length
geom->SetTopVolume(worldBox);
//--- making axis
@ -183,7 +187,7 @@ inline void ANASEN::Construct3DModel(int anodeID1, int anodeID2, int cathodeID1,
new TGeoRotation("rot1", wirePhi , wireTheta, 0.)));
}
TGeoVolume * sx3Det = geom->MakeBox("box", Al, 0.1, sx3->GetWidth()/2, sx3->GetLength()/2);
TGeoVolume * sx3Det = geom->MakeBox("box", Si, 0.1, sx3->GetWidth()/2, sx3->GetLength()/2);
sx3Det->SetLineColor(kGreen+3);
for( int i = 0; i < sx3->GetNumDet(); i++){
@ -243,7 +247,7 @@ inline void ANASEN::DrawTrack(TVector3 pos, TVector3 direction, bool drawEstima
Track->SetLineColor(kRed);
worldBox->AddNode(Track, 1, new TGeoCombiTrans( pos.X(), pos.Y(), pos.Z(), new TGeoRotation("rotA", phi + 90, theta, 0.)));
TGeoVolume * startPos = geom->MakeSphere("startPos", 0, 0, 3);
TGeoVolume * startPos = geom->MakeSphere("startPos", 0, 0, 3);
startPos->SetLineColor(kBlack);
worldBox->AddNode(startPos, 3, new TGeoCombiTrans( pos.X(), pos.Y(), pos.Z(), new TGeoRotation("rotA", 0, 0, 0.)));

2
Armory/ClassTransfer.h
View File

@ -40,7 +40,7 @@ public:
std::string beamStoppingPowerFile; ///stopping_power_for_beam
std::string recoilLightStoppingPowerFile; ///stopping_power_for_light_recoil
std::string recoilHeavyStoppingPowerFile; ///stopping_power_for_heavy_recoil
bool isDecay; ///isDacay
bool isDecay; ///isDecay
int heavyDecayA; ///decayNucleus_A
int heavyDecayZ; ///decayNucleus_Z
bool isRedo; ///isReDo

58
Armory/anasenMS.cpp
View File

@ -34,15 +34,15 @@ int main(int argc, char **argv){
// load energy loss tables (assume units: E in MeV, dE/dx in MeV/(mg/cm²), density in mg/cm³)
TGraph* elossLight = LoadELoss("../ELoss/Eloss_alpha"); // for light particle (alpha)
TGraph* elossHeavy = LoadELoss("../ELoss/Eloss_p"); // for heavy particle (proton)
double density = 2.7; // example for aluminum target, adjust as needed
double density = 0.0000876; // example for aluminum target, adjust as needed
// Reaction setup: projectile + target configuration, energy, and product IDs
TransferReaction transfer;
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
transfer.Seta( 4, 2); // identify reaction product a in internal indexing e.g., 4He (alpha)
transfer.Setb( 1, 1); // identify reaction product b e.g., 1H (proton)
// TODO add alpha source or alternative reaction channel selection
@ -100,7 +100,7 @@ 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
// beam and CM variables saved in tree
double KEA;
@ -175,6 +175,8 @@ int main(int argc, char **argv){
double z0;
tree->Branch("z0", &z0, "reconstucted_Z/D");
TTree* tree2 = tree->CloneTree(0);
tree2->SetName("tree2");
//========timer
TBenchmark clock;
@ -268,10 +270,34 @@ int main(int argc, char **argv){
visTrackHitPos.push_back(hitPos);
visTrackWires.push_back({anodeID[0], cathodeID[0]});
}
/*
// 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());
}
// update kinetic energy after loss
Tb = Pb.E() - Pb.M();*/
// fill tree with original data before energy loss
// 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();
z0 = pw->GetZ0();
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();
@ -286,19 +312,6 @@ int main(int argc, char **argv){
}
// update kinetic energy after loss
Tb = Pb.E() - Pb.M();
// 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();
z0 = pw->GetZ0();
// fill tree2 with energy loss adjusted data
tree2->Fill();
@ -321,6 +334,7 @@ int main(int argc, char **argv){
// fill tree with original data (no energy loss for these events)
tree->Fill();
tree2->Fill();
}
//#################################################################### Timer
@ -343,8 +357,10 @@ int main(int argc, char **argv){
}
// write results to ROOT file and close
tree->Write();
tree2->Write();
//tree->Write();
//tree2->Write();
tree->Write("", TObject::kOverwrite);
tree2->Write("", TObject::kOverwrite);
int count = tree->GetEntries();
int count2 = tree2->GetEntries();
saveFile->Close();

89
Armory/histcomp.C Normal file
View File

@ -0,0 +1,89 @@
void histcomp() {
// Open file
TFile *f = new TFile("SimAnasen1.root");
// Get trees (MAKE SURE names are correct)
TTree *tree1 = (TTree*)f->Get("tree");
TTree *tree2 = (TTree*)f->Get("tree2");
if (!tree1 || !tree2) {
printf("Error: could not find trees. Check names!\n");
return;
}
// Create output directory (overwrite-safe)
gSystem->Exec("mkdir -p plots");
// Get list of branches
TObjArray *branches = tree1->GetListOfBranches();
int nBranches = branches->GetEntries();
//int nBranches = 1;
// Loop over branches
for (int i = 0; i < nBranches; i++) {
TBranch *br = (TBranch*)branches->At(i);
TString name = br->GetName();
printf("Processing branch: %s\n", name.Data());
// Create histograms (auto-range using Draw first)
TString h1name = "h1_" + name;
TString h2name = "h2_" + name;
// Temporary draw to get range
tree1->Draw(name, "", "goff");
double min = tree1->GetMinimum(name);
double max = tree1->GetMaximum(name);
//if (min == max) continue; // skip constant branches
// Expand range slightly
double margin = 0.1 * (max - min);
min -= margin;
max += margin;
TH1D *h1 = new TH1D(h1name, name, 100, min, max);
TH1D *h2 = new TH1D(h2name, name, 100, min, max);
// Fill histograms
tree1->Draw(name + ">>" + h1name, "", "goff");
tree2->Draw(name + ">>" + h2name, "", "goff");
// Style
h1->SetLineColor(kRed);
h1->SetLineWidth(2);
h2->SetLineColor(kBlue);
h2->SetLineWidth(2);
// Normalize (optional but useful)
//if (h1->GetEntries() > 0) h1->Scale(1.0 / h1->GetEntries());
//if (h2->GetEntries() > 0) h2->Scale(1.0 / h2->GetEntries());
// Canvas
TCanvas *c = new TCanvas("c", name, 800, 600);
h1->SetTitle(name + ";"+name+";Counts");
h1->Draw("HIST");
h2->Draw("HIST SAME");
// Legend
TLegend *leg = new TLegend(0.65,0.75,0.88,0.88);
leg->AddEntry(h1, "tree1", "l");
leg->AddEntry(h2, "tree2", "l");
leg->Draw();
// Save plot (overwrite each run)
TString filename = "plots/" + name + ".png";
c->SaveAs(filename);
// Clean up
delete c;
delete h1;
delete h2;
}
printf("Done! Plots saved in ./plots/\n");
}