sequential decay
This commit is contained in:
James Szalkie
parent
d2953531ef
commit
4feca6c104
|
|
@ -175,7 +175,7 @@ public:
|
||||||
void SetA(int A, int Z, double Ex);
|
void SetA(int A, int Z, double Ex);
|
||||||
void Seta(int A, int Z);
|
void Seta(int A, int Z);
|
||||||
void Setb(int A, int Z);
|
void Setb(int A, int Z);
|
||||||
void SetB(int A, int Z);
|
void SetB(int A, int Z, double Ex);
|
||||||
void SetIncidentEnergyAngle(double KEA, double theta, double phi);
|
void SetIncidentEnergyAngle(double KEA, double theta, double phi);
|
||||||
void SetExA(double Ex); // excitation energy of A in MeV
|
void SetExA(double Ex); // excitation energy of A in MeV
|
||||||
void SetExB(double Ex); // excitation energy of B in MeV
|
void SetExB(double Ex); // excitation energy of B in MeV
|
||||||
|
|
@ -245,7 +245,7 @@ TransferReaction::TransferReaction(){
|
||||||
SetA(24, 12, 0);
|
SetA(24, 12, 0);
|
||||||
Seta(4,2);
|
Seta(4,2);
|
||||||
Setb(1,1);
|
Setb(1,1);
|
||||||
SetB(27,13);
|
SetB(27,13, 0);
|
||||||
TA = 2.5; // MeV/u
|
TA = 2.5; // MeV/u
|
||||||
T = TA * reaction.beamA;
|
T = TA * reaction.beamA;
|
||||||
|
|
||||||
|
|
@ -300,12 +300,14 @@ void TransferReaction::Setb(int A, int Z){
|
||||||
isReady = false;
|
isReady = false;
|
||||||
isBSet = false;
|
isBSet = false;
|
||||||
}
|
}
|
||||||
void TransferReaction::SetB(int A, int Z){
|
void TransferReaction::SetB(int A, int Z, double Ex = 0){
|
||||||
Isotope temp (A, Z);
|
Isotope temp (A, Z);
|
||||||
mB = temp.Mass;
|
double mB0 = temp.Mass; // ground state mass
|
||||||
|
mB = mB0;
|
||||||
reaction.recoilHeavyA = A;
|
reaction.recoilHeavyA = A;
|
||||||
reaction.recoilHeavyZ = Z;
|
reaction.recoilHeavyZ = Z;
|
||||||
nameB = temp.Name;
|
nameB = temp.Name;
|
||||||
|
ExB = Ex;
|
||||||
isReady = false;
|
isReady = false;
|
||||||
isBSet = true;
|
isBSet = true;
|
||||||
}
|
}
|
||||||
|
|
@ -389,7 +391,8 @@ void TransferReaction::CalReactionConstant(){
|
||||||
beta = k / (mA + ExA + ma + T);
|
beta = k / (mA + ExA + ma + T);
|
||||||
gamma = 1 / TMath::Sqrt(1- beta * beta);
|
gamma = 1 / TMath::Sqrt(1- beta * beta);
|
||||||
Etot = TMath::Sqrt(TMath::Power(mA + ExA + ma + T,2) - k * k);
|
Etot = TMath::Sqrt(TMath::Power(mA + ExA + ma + T,2) - k * k);
|
||||||
p = TMath::Sqrt( (Etot*Etot - TMath::Power(mb + mB + ExB,2)) * (Etot*Etot - TMath::Power(mb - mB - ExB,2)) ) / 2 / Etot;
|
double mBtot = mB + ExB;
|
||||||
|
p = TMath::Sqrt( (Etot*Etot - TMath::Power(mb + mBtot,2)) * (Etot*Etot - TMath::Power(mb - mBtot,2)) ) / 2 / Etot;
|
||||||
|
|
||||||
PA.SetXYZM(0, 0, k, mA + ExA);
|
PA.SetXYZM(0, 0, k, mA + ExA);
|
||||||
PA.RotateY(thetaIN);
|
PA.RotateY(thetaIN);
|
||||||
|
|
|
||||||
632
Armory/anasenMS Decay Version.cpp
Normal file
632
Armory/anasenMS Decay Version.cpp
Normal file
|
|
@ -0,0 +1,632 @@
|
||||||
|
#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
|
||||||
|
#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 <cstring>
|
||||||
|
#include "TApplication.h" // ROOT app loop
|
||||||
|
#include "ClassTransfer.h" // Reaction kinematics and MC event generation
|
||||||
|
#include "ClassAnasen.h" // ANASEN detector model classes (SX3, PW, etc.)
|
||||||
|
#include <stdio.h>
|
||||||
|
#include <stdlib.h>
|
||||||
|
#include <set>
|
||||||
|
#include "TLegend.h"
|
||||||
|
#include "TH1D.h"
|
||||||
|
#include "TObjArray.h"
|
||||||
|
#include "TBranch.h"
|
||||||
|
#include <iostream>
|
||||||
|
#include <fstream>
|
||||||
|
|
||||||
|
//======== Generate light particle based on reaction
|
||||||
|
// calculate real and reconstructed tracks and Q-value uncertainty
|
||||||
|
|
||||||
|
// Function to load energy loss table from file
|
||||||
|
TGraph* LoadELoss(const char* filename) {
|
||||||
|
TGraph* g = new TGraph(filename, "%lg %lg");
|
||||||
|
return g;
|
||||||
|
}
|
||||||
|
|
||||||
|
bool IsDeadAnode(int id){
|
||||||
|
static std::set<int> dead = {}; // add dead anode IDs here, 0-23
|
||||||
|
return dead.count(id);
|
||||||
|
}
|
||||||
|
|
||||||
|
bool IsDeadCathode(int id){
|
||||||
|
static std::set<int> dead = {}; // add dead cathode IDs here, 0-23
|
||||||
|
return dead.count(id);
|
||||||
|
}
|
||||||
|
|
||||||
|
bool IsDeadSX3(int id){
|
||||||
|
static std::set<int> dead = {}; // add dead SX3 IDs here, 0-23 1,7,9,3
|
||||||
|
return dead.count(id);
|
||||||
|
}
|
||||||
|
|
||||||
|
// Simulate sequential two-body decay of an unstable parent in its rest frame.
|
||||||
|
// The parent is boosted from the lab frame, the daughter (A1,Z1) is returned in lab frame,
|
||||||
|
// and the emitted ejectile (A2,Z2) is written to ejectileOut.
|
||||||
|
TLorentzVector SimulateSequentialDecay(const TLorentzVector &parent,
|
||||||
|
int daughterA, int daughterZ,
|
||||||
|
int ejectA, int ejectZ,
|
||||||
|
TLorentzVector &ejectileOut){
|
||||||
|
Isotope daughter(daughterA, daughterZ);
|
||||||
|
Isotope ejectile(ejectA, ejectZ);
|
||||||
|
|
||||||
|
double M = parent.M();
|
||||||
|
double mD = daughter.Mass;
|
||||||
|
double mE = ejectile.Mass;
|
||||||
|
|
||||||
|
double sqM = M * M;
|
||||||
|
double sum = mD + mE;
|
||||||
|
double diff = mD - mE;
|
||||||
|
double p2 = (sqM - sum*sum) * (sqM - diff*diff) / (4.0 * sqM);
|
||||||
|
if( p2 < 0 ) p2 = 0;
|
||||||
|
double p = TMath::Sqrt(p2);
|
||||||
|
|
||||||
|
double cosTheta = 2.0 * gRandom->Rndm() - 1.0;
|
||||||
|
double theta = TMath::ACos(cosTheta);
|
||||||
|
double phi = gRandom->Rndm() * TMath::TwoPi();
|
||||||
|
|
||||||
|
TVector3 v;
|
||||||
|
v.SetMagThetaPhi(p, theta, phi);
|
||||||
|
|
||||||
|
TLorentzVector daughterLab;
|
||||||
|
daughterLab.SetVectM(v, mD);
|
||||||
|
|
||||||
|
TLorentzVector ejectileLab;
|
||||||
|
ejectileLab.SetVectM(-v, mE);
|
||||||
|
|
||||||
|
TVector3 boost = parent.BoostVector();
|
||||||
|
daughterLab.Boost(boost);
|
||||||
|
ejectileLab.Boost(boost);
|
||||||
|
|
||||||
|
ejectileOut = ejectileLab;
|
||||||
|
return daughterLab;
|
||||||
|
}
|
||||||
|
|
||||||
|
int main(int argc, char **argv){
|
||||||
|
|
||||||
|
printf("=========================================\n");
|
||||||
|
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 setup for 18Ne + 4He -> p + 21Na*.
|
||||||
|
// The heavy product 21Na* is then decayed to 20Ne + p in the simulation.
|
||||||
|
TransferReaction transfer;
|
||||||
|
|
||||||
|
transfer.SetA(18, 10, 0); // 18Ne
|
||||||
|
transfer.SetIncidentEnergyAngle(4.4, 0, 0); // KEA in MeV/u, theta and phi in degree
|
||||||
|
transfer.Seta(4, 2); // 4He target
|
||||||
|
transfer.Setb(1, 1); // outgoing proton from the primary transfer
|
||||||
|
transfer.SetB(21, 11); // 21Na* heavy product
|
||||||
|
|
||||||
|
bool enableSequentialDecay = true;
|
||||||
|
const int decayDaughterA = 20;
|
||||||
|
const int decayDaughterZ = 10;
|
||||||
|
const int decayEjectA = 1;
|
||||||
|
const int decayEjectZ = 1;
|
||||||
|
|
||||||
|
// Excited state lists (projectile and heavy-product excitation states)
|
||||||
|
std::vector<float> ExAList = {0}; // 18Ne projectile excitations in MeV
|
||||||
|
std::vector<float> ExList = {2.5}; // 21Na* excitation in MeV (above proton separation threshold)
|
||||||
|
|
||||||
|
// define vertex position uniform distribution ranges (mm)
|
||||||
|
double vertexXRange[2] = { -5, 5}; // mm
|
||||||
|
double vertexYRange[2] = { -5, 5};
|
||||||
|
double vertexZRange[2] = { -100, 100};
|
||||||
|
|
||||||
|
// 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]);
|
||||||
|
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",numEvent);
|
||||||
|
|
||||||
|
// calculates energy/momentum/kinematics constants for transfer reaction
|
||||||
|
transfer.CalReactionConstant();
|
||||||
|
printf("Primary reaction: %s at %.2f MeV/u\n", transfer.GetReactionName().Data(), 4.4);
|
||||||
|
printf("Sequential decay enabled: %s\n", enableSequentialDecay ? "yes" : "no");
|
||||||
|
|
||||||
|
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);
|
||||||
|
}
|
||||||
|
|
||||||
|
// storage for tracks during simulation (for visualization)
|
||||||
|
std::vector<TVector3> visTrackVertex, visTrackDir, visTrackHitPos;
|
||||||
|
std::vector<std::pair<int,int>> visTrackWires; // {anodeID, cathodeID}
|
||||||
|
|
||||||
|
// 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 + trees
|
||||||
|
TString saveFileName = "SimAnasen1.root";
|
||||||
|
printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
|
||||||
|
TFile * saveFile = new TFile(saveFileName, "recreate");
|
||||||
|
TTree * tree1 = new TTree("tree1", "tree1");
|
||||||
|
TTree * tree2 = new TTree("tree2", "tree2");
|
||||||
|
|
||||||
|
// beam and CM variables saved in tree
|
||||||
|
double KEA;
|
||||||
|
double KEA2;
|
||||||
|
tree1->Branch("beamKEA", &KEA, "beamKEA/D");
|
||||||
|
tree2->Branch("beamKEA", &KEA2, "beamKEA/D");
|
||||||
|
|
||||||
|
double thetaCM, phiCM;
|
||||||
|
double thetaCM2, phiCM2;
|
||||||
|
tree1->Branch("thetaCM", &thetaCM, "thetaCM/D");
|
||||||
|
tree1->Branch("phiCM", &phiCM, "phiCM/D");
|
||||||
|
tree2->Branch("thetaCM", &thetaCM2, "thetaCM/D");
|
||||||
|
tree2->Branch("phiCM", &phiCM2, "phiCM/D");
|
||||||
|
|
||||||
|
// outgoing particles in lab frame (light/heavy)
|
||||||
|
double thetab, phib;
|
||||||
|
double Tb;
|
||||||
|
double thetaB, phiB, TB;
|
||||||
|
std::array<double, 2> T;
|
||||||
|
tree1->Branch("thetab", &thetab, "thetab/D"); // polar angle of light particle in lab frame
|
||||||
|
tree1->Branch("phib", &phib, "phib/D"); // azimuthal angle of light particle in lab frame
|
||||||
|
tree1->Branch("Tb", &Tb, "Tb/D"); // kinetic energy of light particle at vertex (before energy loss)
|
||||||
|
tree1->Branch("thetaB", &thetaB, "thetaB/D");
|
||||||
|
tree1->Branch("phiB", &phiB, "phiB/D");
|
||||||
|
tree1->Branch("TB", &TB, "TB/D"); // kinetic energy of heavy particle at vertex
|
||||||
|
tree1->Branch("T", &T, "T/D"); // placeholder for true Q-value, currently set to 0 for simplicity
|
||||||
|
|
||||||
|
double thetab2, phib2;
|
||||||
|
double Tb2;
|
||||||
|
double thetaB2, phiB2, TB2;
|
||||||
|
std::array<double, 2> T2;
|
||||||
|
tree2->Branch("thetab", &thetab2, "thetab/D");
|
||||||
|
tree2->Branch("phib", &phib2, "phib/D");
|
||||||
|
tree2->Branch("Tb", &Tb2, "Tb/D");
|
||||||
|
tree2->Branch("thetaB", &thetaB2, "thetaB/D");
|
||||||
|
tree2->Branch("phiB", &phiB2, "phiB/D");
|
||||||
|
tree2->Branch("TB", &TB2, "TB/D");
|
||||||
|
tree2->Branch("T", &T2, "T/D");
|
||||||
|
|
||||||
|
// excitation state identifiers
|
||||||
|
int ExAID;
|
||||||
|
double ExA;
|
||||||
|
tree1->Branch("ExAID", &ExAID, "ExAID/I"); // projectile excitation state ID
|
||||||
|
tree1->Branch("ExA", &ExA, "ExA/D"); // projectile excitation energy in MeV
|
||||||
|
|
||||||
|
int ExAID2;
|
||||||
|
double ExA2;
|
||||||
|
tree2->Branch("ExAID", &ExAID2, "ExAID/I");
|
||||||
|
tree2->Branch("ExA", &ExA2, "ExA/D");
|
||||||
|
|
||||||
|
int ExID;
|
||||||
|
double Ex;
|
||||||
|
tree1->Branch("ExID", &ExID, "ExID/I"); // target excitation state ID
|
||||||
|
tree1->Branch("Ex", &Ex, "Ex/D"); // target excitation energy in MeV
|
||||||
|
|
||||||
|
int ExID2;
|
||||||
|
double Ex2;
|
||||||
|
tree2->Branch("ExID", &ExID2, "ExID/I");
|
||||||
|
tree2->Branch("Ex", &Ex2, "Ex/D");
|
||||||
|
|
||||||
|
// true vertex position in target volume
|
||||||
|
double vertexX, vertexY, vertexZ;
|
||||||
|
tree1->Branch("vX", &vertexX, "VertexX/D"); // true vertex X position in mm
|
||||||
|
tree1->Branch("vY", &vertexY, "VertexY/D"); // true vertex Y position in mm
|
||||||
|
tree1->Branch("vZ", &vertexZ, "VertexZ/D"); // true vertex Z position in mm
|
||||||
|
|
||||||
|
double vertexX2, vertexY2, vertexZ2;
|
||||||
|
tree2->Branch("vX", &vertexX2, "VertexX/D");
|
||||||
|
tree2->Branch("vY", &vertexY2, "VertexY/D");
|
||||||
|
tree2->Branch("vZ", &vertexZ2, "VertexZ/D");
|
||||||
|
|
||||||
|
// reconstructed SX3 hit position
|
||||||
|
double sx3X, sx3Y, sx3Z;
|
||||||
|
tree1->Branch("sx3X", &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing)
|
||||||
|
tree1->Branch("sx3Y", &sx3Y, "sx3Y/D"); // reconstructed Y position from SX3 (with optional smearing)
|
||||||
|
tree1->Branch("sx3Z", &sx3Z, "sx3Z/D"); // reconstructed Z position from SX3 (with optional smearing)
|
||||||
|
|
||||||
|
double sx3X2, sx3Y2, sx3Z2;
|
||||||
|
tree2->Branch("sx3X", &sx3X2, "sx3X/D");
|
||||||
|
tree2->Branch("sx3Y", &sx3Y2, "sx3Y/D");
|
||||||
|
tree2->Branch("sx3Z", &sx3Z2, "sx3Z/D");
|
||||||
|
|
||||||
|
// PW nearest and next nearest wires
|
||||||
|
int anodeID[2], cathodeID[2];
|
||||||
|
int anodeID2[2], cathodeID2[2];
|
||||||
|
tree1->Branch("aID", anodeID, "anodeID/I"); // anodeID[0] is nearest anode wire, anodeID[1] is next nearest anode wire
|
||||||
|
tree1->Branch("cID", cathodeID, "cathodeID/I"); // cathodeID[0] is nearest cathode wire, cathodeID[1] is next nearest cathode wire
|
||||||
|
tree2->Branch("aID", anodeID2, "anodeID/I");
|
||||||
|
tree2->Branch("cID", cathodeID2, "cathodeID/I");
|
||||||
|
|
||||||
|
// distances to nearest wires
|
||||||
|
double anodeDist[2], cathodeDist[2];
|
||||||
|
double anodeDist2[2], cathodeDist2[2];
|
||||||
|
tree1->Branch("aDist", anodeDist, "anodeDist/D");
|
||||||
|
tree1->Branch("cDist", cathodeDist, "cathodeDist/D");
|
||||||
|
tree2->Branch("aDist", anodeDist2, "anodeDist/D");
|
||||||
|
tree2->Branch("cDist", cathodeDist2, "cathodeDist/D");
|
||||||
|
|
||||||
|
// SX3 channel assignment and Z fraction (depth) information
|
||||||
|
int sx3ID, sx3Up, sx3Dn, sx3Bk;
|
||||||
|
double sx3ZFrac;
|
||||||
|
int sx3ID2, sx3Up2, sx3Dn2, sx3Bk2;
|
||||||
|
double sx3ZFrac2;
|
||||||
|
tree1->Branch("sx3ID", &sx3ID, "sx3ID/I");
|
||||||
|
tree1->Branch("sx3Up", &sx3Up, "sx3Up/I");
|
||||||
|
tree1->Branch("sx3Dn", &sx3Dn, "sx3Dn/I");
|
||||||
|
tree1->Branch("sx3Bk", &sx3Bk, "sx3Bk/I");
|
||||||
|
tree1->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D");
|
||||||
|
tree2->Branch("sx3ID", &sx3ID2, "sx3ID/I");
|
||||||
|
tree2->Branch("sx3Up", &sx3Up2, "sx3Up/I");
|
||||||
|
tree2->Branch("sx3Dn", &sx3Dn2, "sx3Dn/I");
|
||||||
|
tree2->Branch("sx3Bk", &sx3Bk2, "sx3Bk/I");
|
||||||
|
tree2->Branch("sx3ZFrac", &sx3ZFrac2, "sx3ZFrac/D");
|
||||||
|
|
||||||
|
// reconstructed angles from PW track fit, method 1 and 2
|
||||||
|
double reTheta, rePhi;
|
||||||
|
double reTheta2, rePhi2;
|
||||||
|
tree1->Branch("reTheta", &reTheta, "reconstucted_theta/D");
|
||||||
|
tree1->Branch("rePhi", &rePhi, "reconstucted_phi/D");
|
||||||
|
tree2->Branch("reTheta", &reTheta2, "reconstucted_theta/D");
|
||||||
|
tree2->Branch("rePhi", &rePhi2, "reconstucted_phi/D");
|
||||||
|
|
||||||
|
double reTheta1, rePhi1;
|
||||||
|
double reTheta12, rePhi12;
|
||||||
|
tree1->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D");
|
||||||
|
tree1->Branch("rePhi1", &rePhi1, "reconstucted_phi1/D");
|
||||||
|
tree2->Branch("reTheta1", &reTheta12, "reconstucted_theta1/D");
|
||||||
|
tree2->Branch("rePhi1", &rePhi12, "reconstucted_phi1/D");
|
||||||
|
|
||||||
|
// reconstructed vertex Z from PW fit
|
||||||
|
double z0;
|
||||||
|
double z02;
|
||||||
|
tree1->Branch("z0", &z0, "reconstucted_Z/D");
|
||||||
|
tree2->Branch("z0", &z02, "reconstucted_Z/D");
|
||||||
|
|
||||||
|
//========timer
|
||||||
|
TBenchmark clock;
|
||||||
|
bool shown ;
|
||||||
|
clock.Reset();
|
||||||
|
clock.Start("timer");
|
||||||
|
shown = false;
|
||||||
|
|
||||||
|
//================================= 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);
|
||||||
|
|
||||||
|
ExID = gRandom->Integer(nEx);
|
||||||
|
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 kinematics in lab frame for the primary transfer
|
||||||
|
TLorentzVector * output = transfer.Event(thetaCM, phiCM); // returns array of outputs
|
||||||
|
TLorentzVector Pb = output[2]; // primary proton from transfer
|
||||||
|
TLorentzVector PB = output[3]; // excited 21Na* heavy product
|
||||||
|
|
||||||
|
thetab = Pb.Theta() * TMath::RadToDeg();
|
||||||
|
Tb = (Pb.E() - Pb.M()); // kinetic energy of the light proton from the primary transfer
|
||||||
|
thetaB = PB.Theta() * TMath::RadToDeg();
|
||||||
|
TB = (PB.E() - PB.M());
|
||||||
|
phib = Pb.Phi() * TMath::RadToDeg();
|
||||||
|
phiB = PB.Phi() * TMath::RadToDeg();
|
||||||
|
T[0] = Tb;
|
||||||
|
T[1] = TB;
|
||||||
|
|
||||||
|
// prepare secondary proton from 21Na* sequential decay
|
||||||
|
TLorentzVector decayProton;
|
||||||
|
TLorentzVector heavy20;
|
||||||
|
if(enableSequentialDecay){
|
||||||
|
heavy20 = SimulateSequentialDecay(PB, decayDaughterA, decayDaughterZ,
|
||||||
|
decayEjectA, decayEjectZ, decayProton);
|
||||||
|
thetab2 = decayProton.Theta() * TMath::RadToDeg();
|
||||||
|
phib2 = decayProton.Phi() * TMath::RadToDeg();
|
||||||
|
Tb2 = decayProton.E() - decayProton.M();
|
||||||
|
thetaB2 = heavy20.Theta() * TMath::RadToDeg();
|
||||||
|
phiB2 = heavy20.Phi() * TMath::RadToDeg();
|
||||||
|
TB2 = heavy20.E() - heavy20.M();
|
||||||
|
T2[0] = Tb2;
|
||||||
|
T2[1] = TB2;
|
||||||
|
} else {
|
||||||
|
thetab2 = TMath::QuietNaN();
|
||||||
|
phib2 = TMath::QuietNaN();
|
||||||
|
Tb2 = TMath::QuietNaN();
|
||||||
|
thetaB2 = TMath::QuietNaN();
|
||||||
|
phiB2 = TMath::QuietNaN();
|
||||||
|
TB2 = TMath::QuietNaN();
|
||||||
|
T2[0] = TMath::QuietNaN();
|
||||||
|
T2[1] = TMath::QuietNaN();
|
||||||
|
}
|
||||||
|
|
||||||
|
delete [] output;
|
||||||
|
|
||||||
|
// 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();
|
||||||
|
|
||||||
|
anodeID[0] = hitInfo.nearestWire.first; // nearest anode wire ID
|
||||||
|
cathodeID[0] = hitInfo.nearestWire.second; // nearest cathode wire ID
|
||||||
|
anodeID[1] = hitInfo.nextNearestWire.first; // next nearest anode wire ID
|
||||||
|
cathodeID[1] = hitInfo.nextNearestWire.second; // next nearest cathode wire ID
|
||||||
|
|
||||||
|
anodeDist[1] = hitInfo.nextNearestDist.first; // distance to next nearest anode wire
|
||||||
|
cathodeDist[1] = hitInfo.nextNearestDist.second; // distance to next nearest cathode wire
|
||||||
|
|
||||||
|
if(IsDeadAnode(anodeID[0])) continue;
|
||||||
|
if(IsDeadCathode(cathodeID[0])) continue;
|
||||||
|
|
||||||
|
// SX3 hit channel info and depth fraction
|
||||||
|
sx3ID = sx3->GetID();
|
||||||
|
|
||||||
|
if(IsDeadSX3(sx3ID)) continue;
|
||||||
|
|
||||||
|
anodeDist[0] = hitInfo.nearestDist.first; // distance to nearest anode wire
|
||||||
|
cathodeDist[0] = hitInfo.nearestDist.second; // distance to nearest cathode wire
|
||||||
|
|
||||||
|
if( sx3ID >= 0 ){
|
||||||
|
sx3Up = sx3->GetChUp();
|
||||||
|
sx3Dn = sx3->GetChDn();
|
||||||
|
sx3Bk = sx3->GetChBk();
|
||||||
|
sx3ZFrac = sx3->GetZFrac();
|
||||||
|
|
||||||
|
// 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();
|
||||||
|
|
||||||
|
// store track data for visualization if enabled
|
||||||
|
if(enableVis){
|
||||||
|
visTrackVertex.push_back(vertex);
|
||||||
|
visTrackDir.push_back(dir);
|
||||||
|
visTrackHitPos.push_back(hitPos);
|
||||||
|
visTrackWires.push_back({anodeID[0], cathodeID[0]});
|
||||||
|
}
|
||||||
|
// 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();
|
||||||
|
tree1->Fill();
|
||||||
|
|
||||||
|
// fill tree2 using the secondary proton (proton 2) track
|
||||||
|
TVector3 dir2(1, 0, 0);
|
||||||
|
dir2.SetTheta(thetab2 * TMath::DegToRad());
|
||||||
|
dir2.SetPhi(phib2 * TMath::DegToRad());
|
||||||
|
|
||||||
|
pw->FindWireID(vertex, dir2, false);
|
||||||
|
sx3->FindSX3Pos(vertex, dir2, false);
|
||||||
|
PWHitInfo hitInfo2 = pw->GetHitInfo();
|
||||||
|
|
||||||
|
anodeID2[0] = hitInfo2.nearestWire.first;
|
||||||
|
cathodeID2[0] = hitInfo2.nearestWire.second;
|
||||||
|
anodeID2[1] = hitInfo2.nextNearestWire.first;
|
||||||
|
cathodeID2[1] = hitInfo2.nextNearestWire.second;
|
||||||
|
|
||||||
|
anodeDist2[1] = hitInfo2.nextNearestDist.first;
|
||||||
|
cathodeDist2[1] = hitInfo2.nextNearestDist.second;
|
||||||
|
|
||||||
|
if(IsDeadAnode(anodeID2[0]) || IsDeadCathode(cathodeID2[0])){
|
||||||
|
sx3ID2 = -1;
|
||||||
|
} else {
|
||||||
|
sx3ID2 = sx3->GetID();
|
||||||
|
}
|
||||||
|
|
||||||
|
if(sx3ID2 < 0 || IsDeadSX3(sx3ID2)){
|
||||||
|
sx3ID2 = -1;
|
||||||
|
sx3Up2 = -1;
|
||||||
|
sx3Dn2 = -1;
|
||||||
|
sx3Bk2 = -1;
|
||||||
|
sx3ZFrac2 = TMath::QuietNaN();
|
||||||
|
sx3X2 = TMath::QuietNaN();
|
||||||
|
sx3Y2 = TMath::QuietNaN();
|
||||||
|
sx3Z2 = TMath::QuietNaN();
|
||||||
|
anodeDist2[0] = TMath::QuietNaN();
|
||||||
|
cathodeDist2[0] = TMath::QuietNaN();
|
||||||
|
reTheta2 = TMath::QuietNaN();
|
||||||
|
rePhi2 = TMath::QuietNaN();
|
||||||
|
reTheta12 = TMath::QuietNaN();
|
||||||
|
rePhi12 = TMath::QuietNaN();
|
||||||
|
z02 = TMath::QuietNaN();
|
||||||
|
} else {
|
||||||
|
anodeDist2[0] = hitInfo2.nearestDist.first;
|
||||||
|
cathodeDist2[0] = hitInfo2.nearestDist.second;
|
||||||
|
sx3Up2 = sx3->GetChUp();
|
||||||
|
sx3Dn2 = sx3->GetChDn();
|
||||||
|
sx3Bk2 = sx3->GetChBk();
|
||||||
|
sx3ZFrac2 = sx3->GetZFrac();
|
||||||
|
TVector3 hitPos2 = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
|
||||||
|
sx3X2 = hitPos2.X();
|
||||||
|
sx3Y2 = hitPos2.Y();
|
||||||
|
sx3Z2 = hitPos2.Z();
|
||||||
|
pw->CalTrack(hitPos2, anodeID2[0], cathodeID2[0], false);
|
||||||
|
reTheta2 = pw->GetTrackTheta() * TMath::RadToDeg();
|
||||||
|
rePhi2 = pw->GetTrackPhi() * TMath::RadToDeg();
|
||||||
|
pw->CalTrack2(hitPos2, hitInfo2, sigmaPW_A, sigmaPW_C, false);
|
||||||
|
reTheta12 = pw->GetTrackTheta() * TMath::RadToDeg();
|
||||||
|
rePhi12 = pw->GetTrackPhi() * TMath::RadToDeg();
|
||||||
|
z02 = pw->GetZ0();
|
||||||
|
}
|
||||||
|
|
||||||
|
// copy common event info to tree2
|
||||||
|
KEA2 = KEA;
|
||||||
|
thetaCM2 = thetaCM;
|
||||||
|
phiCM2 = phiCM;
|
||||||
|
ExAID2 = ExAID;
|
||||||
|
ExA2 = ExA;
|
||||||
|
ExID2 = ExID;
|
||||||
|
Ex2 = Ex;
|
||||||
|
vertexX2 = vertexX;
|
||||||
|
vertexY2 = vertexY;
|
||||||
|
vertexZ2 = vertexZ;
|
||||||
|
|
||||||
|
tree2->Fill();
|
||||||
|
|
||||||
|
}else{
|
||||||
|
// no valid SX3 hit: mark clearly invalid
|
||||||
|
sx3Up = -1;
|
||||||
|
sx3Dn = -1;
|
||||||
|
sx3Bk = -1;
|
||||||
|
sx3ZFrac = TMath::QuietNaN();
|
||||||
|
|
||||||
|
sx3X = TMath::QuietNaN();
|
||||||
|
sx3Y = TMath::QuietNaN();
|
||||||
|
sx3Z = TMath::QuietNaN();
|
||||||
|
|
||||||
|
reTheta = TMath::QuietNaN();
|
||||||
|
rePhi = TMath::QuietNaN();
|
||||||
|
reTheta1 = TMath::QuietNaN();
|
||||||
|
rePhi1 = TMath::QuietNaN();
|
||||||
|
z0 = TMath::QuietNaN();
|
||||||
|
//Tb = -12354567; // mark kinetic energy as invalid for no hit case
|
||||||
|
// fill tree with original data (no energy loss for these events)
|
||||||
|
//comment out tree fill for no hit case
|
||||||
|
//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");
|
||||||
|
|
||||||
|
if ( !shown ) {
|
||||||
|
if (fmod(time, 10) < 1 ){
|
||||||
|
printf( "%10d[%2d%%]| %8.2f sec | expect: %5.1f min \n", i, TMath::Nint((i+1)*100./numEvent), time , numEvent*time/(i+1)/60);
|
||||||
|
shown = 1;
|
||||||
|
}
|
||||||
|
} else {
|
||||||
|
if (fmod(time, 10) > 9 ){
|
||||||
|
shown = 0;
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
// write results to ROOT file and close
|
||||||
|
tree1->Write("", TObject::kOverwrite);
|
||||||
|
tree2->Write("", TObject::kOverwrite);
|
||||||
|
int count1 = tree1->GetEntries();
|
||||||
|
int count2 = tree2->GetEntries();
|
||||||
|
saveFile->Close();
|
||||||
|
|
||||||
|
printf("=============== done. saved as %s. tree1 entries: %d, tree2 entries: %d\n", saveFileName.Data(), count1, count2);
|
||||||
|
|
||||||
|
if(enableVis){ // to enable visualization, run with 3rd argument "vis", e.g. "./anasenMC 1000 vis"
|
||||||
|
printf("Displaying geometry with %zu tracks from simulation\n", visTrackVertex.size());
|
||||||
|
|
||||||
|
// Build full geometry with all wires
|
||||||
|
anasen->DrawAnasen(0, 23, 0, 23, -1, true);
|
||||||
|
|
||||||
|
// Add all stored tracks to the geometry
|
||||||
|
TGeoManager *geom = anasen->GetGeoManager();
|
||||||
|
TGeoVolume *worldBox = anasen->GetWorldBox();
|
||||||
|
|
||||||
|
if(geom && worldBox && visTrackVertex.size() > 0){
|
||||||
|
int trackNodeID = 500; // start node IDs for tracks
|
||||||
|
|
||||||
|
for(size_t iTrack = 0; iTrack < visTrackVertex.size(); ++iTrack){
|
||||||
|
TVector3 vertex = visTrackVertex[iTrack];
|
||||||
|
TVector3 dir = visTrackDir[iTrack];
|
||||||
|
TVector3 hitPos = visTrackHitPos[iTrack];
|
||||||
|
|
||||||
|
double theta = dir.Theta() * TMath::RadToDeg();
|
||||||
|
double phi = dir.Phi() * TMath::RadToDeg();
|
||||||
|
|
||||||
|
// Add a line marker at the vertex
|
||||||
|
TGeoVolume *startMarker = geom->MakeSphere("startMarker", 0, 0, 2.0);
|
||||||
|
startMarker->SetLineColor(kBlack);
|
||||||
|
worldBox->AddNode(startMarker, trackNodeID,
|
||||||
|
new TGeoCombiTrans(vertex.X(), vertex.Y(), vertex.Z(),
|
||||||
|
new TGeoRotation("rot", 0, 0, 0)));
|
||||||
|
trackNodeID++;
|
||||||
|
|
||||||
|
// Add track line from vertex toward hit position
|
||||||
|
TGeoVolume *trackLine = geom->MakeTube("trackLine", 0, 0, 0.08, 150.0);
|
||||||
|
trackLine->SetLineColor(kBlue);
|
||||||
|
worldBox->AddNode(trackLine, trackNodeID,
|
||||||
|
new TGeoCombiTrans(vertex.X(), vertex.Y(), vertex.Z(),
|
||||||
|
new TGeoRotation("rotTrack", phi + 90, theta, 0)));
|
||||||
|
trackNodeID++;
|
||||||
|
|
||||||
|
// Add hit position marker
|
||||||
|
TGeoVolume *hitMarker = geom->MakeSphere("hitMarker", 0, 0, 2.0);
|
||||||
|
hitMarker->SetLineColor(kRed);
|
||||||
|
worldBox->AddNode(hitMarker, trackNodeID,
|
||||||
|
new TGeoCombiTrans(hitPos.X(), hitPos.Y(), hitPos.Z(),
|
||||||
|
new TGeoRotation("rotHit", 0, 0, 0)));
|
||||||
|
trackNodeID++;
|
||||||
|
}
|
||||||
|
|
||||||
|
// Redraw geometry with all tracks
|
||||||
|
geom->CloseGeometry();
|
||||||
|
geom->SetVisLevel(4);
|
||||||
|
worldBox->Draw("ogle");
|
||||||
|
}
|
||||||
|
|
||||||
|
if(app){
|
||||||
|
printf("Entering ROOT event loop\n");
|
||||||
|
app->Run();
|
||||||
|
}
|
||||||
|
}
|
||||||
|
|
||||||
|
delete anasen;
|
||||||
|
|
||||||
|
|
||||||
|
return 0;
|
||||||
|
|
||||||
|
}
|
||||||
|
|
@ -45,6 +45,48 @@ bool IsDeadSX3(int id){
|
||||||
return dead.count(id);
|
return dead.count(id);
|
||||||
}
|
}
|
||||||
|
|
||||||
|
// Simulate sequential two-body decay of an unstable parent in its rest frame.
|
||||||
|
// The parent is boosted from the lab frame, the daughter (A1,Z1) is returned in lab frame,
|
||||||
|
// and the emitted ejectile (A2,Z2) is written to ejectileOut.
|
||||||
|
TLorentzVector SimulateSequentialDecay(const TLorentzVector &parent,
|
||||||
|
int daughterA, int daughterZ,
|
||||||
|
int ejectA, int ejectZ,
|
||||||
|
TLorentzVector &ejectileOut){
|
||||||
|
Isotope daughter(daughterA, daughterZ);
|
||||||
|
Isotope ejectile(ejectA, ejectZ);
|
||||||
|
|
||||||
|
double M = parent.M();
|
||||||
|
double mD = daughter.Mass;
|
||||||
|
double mE = ejectile.Mass;
|
||||||
|
|
||||||
|
double sqM = M * M;
|
||||||
|
double sum = mD + mE;
|
||||||
|
double diff = mD - mE;
|
||||||
|
double p2 = (sqM - sum*sum) * (sqM - diff*diff) / (4.0 * sqM);
|
||||||
|
if( p2 < 0 ) p2 = 0;
|
||||||
|
double p = TMath::Sqrt(p2);
|
||||||
|
|
||||||
|
double cosTheta = 2.0 * gRandom->Rndm() - 1.0;
|
||||||
|
double theta = TMath::ACos(cosTheta);
|
||||||
|
double phi = gRandom->Rndm() * TMath::TwoPi();
|
||||||
|
|
||||||
|
TVector3 v;
|
||||||
|
v.SetMagThetaPhi(p, theta, phi);
|
||||||
|
|
||||||
|
TLorentzVector daughterLab;
|
||||||
|
daughterLab.SetVectM(v, mD);
|
||||||
|
|
||||||
|
TLorentzVector ejectileLab;
|
||||||
|
ejectileLab.SetVectM(-v, mE);
|
||||||
|
|
||||||
|
TVector3 boost = parent.BoostVector();
|
||||||
|
daughterLab.Boost(boost);
|
||||||
|
ejectileLab.Boost(boost);
|
||||||
|
|
||||||
|
ejectileOut = ejectileLab;
|
||||||
|
return daughterLab;
|
||||||
|
}
|
||||||
|
|
||||||
int main(int argc, char **argv){
|
int main(int argc, char **argv){
|
||||||
|
|
||||||
printf("=========================================\n");
|
printf("=========================================\n");
|
||||||
|
|
@ -54,20 +96,25 @@ int main(int argc, char **argv){
|
||||||
// number of events can be overridden from command line
|
// number of events can be overridden from command line
|
||||||
int numEvent = 1000000;
|
int numEvent = 1000000;
|
||||||
if( argc >= 2 ) numEvent = atoi(argv[1]);
|
if( argc >= 2 ) numEvent = atoi(argv[1]);
|
||||||
// Reaction setup: projectile + target configuration, energy, and product IDs
|
// Reaction setup: 18Ne + 4He -> p + 21Na*
|
||||||
|
// then sequential decay of 21Na* -> p + 20Ne
|
||||||
TransferReaction transfer;
|
TransferReaction transfer;
|
||||||
|
|
||||||
transfer.SetA(14, 7, 0); // e.g., 24Mg (Z=12) with 0 excitation
|
transfer.SetA(27, 13, 0); // 18Ne projectile
|
||||||
transfer.SetIncidentEnergyAngle((42.82/14.0), 0, 0); // arguments are KEA in MeV/u, theta and phi in degree
|
transfer.SetIncidentEnergyAngle((72/27.0), 0, 0); // KEA in MeV/u, theta and phi in degree
|
||||||
transfer.Seta( 4, 2); // identify reaction product a in internal indexing e.g., 4He (alpha)
|
transfer.Seta(4, 2); // 4He target
|
||||||
transfer.Setb(1, 1); // identify reaction product b e.g., 1H (proton)
|
transfer.Setb(4, 2); // outgoing proton from the primary transfer
|
||||||
transfer.SetB(14, 7); // identify reaction product B e.g., 23Na (Z=11)
|
transfer.SetB(27, 13); // 21Na* heavy product
|
||||||
|
|
||||||
// TODO add alpha source or alternative reaction channel selection
|
bool enableSequentialDecay = false; // turning to false to disable sequential decay for now, can be set to true to enable
|
||||||
|
const int decayDaughterA = 20;
|
||||||
|
const int decayDaughterZ = 10;
|
||||||
|
const int decayEjectA = 1;
|
||||||
|
const int decayEjectZ = 1;
|
||||||
|
|
||||||
// Excited state lists (target and projectile/excited products)
|
// Excited state lists (projectile and heavy-product excitation states)
|
||||||
std::vector<float> ExAList = {0}; // projectile excitation states in MeV
|
std::vector<float> ExAList = {0}; // 18Ne projectile excitations in MeV
|
||||||
std::vector<float> ExList = {0}; // target excitation states in MeV
|
std::vector<float> ExList = {2.5}; // 21Na* excitation in MeV
|
||||||
|
|
||||||
// define vertex position uniform distribution ranges (mm)
|
// define vertex position uniform distribution ranges (mm)
|
||||||
double vertexXRange[2] = { -5, 5}; // mm
|
double vertexXRange[2] = { -5, 5}; // mm
|
||||||
|
|
@ -114,87 +161,145 @@ int main(int argc, char **argv){
|
||||||
SX3 * sx3 = anasen->GetSX3(); // silicon array part
|
SX3 * sx3 = anasen->GetSX3(); // silicon array part
|
||||||
PW * pw = anasen->GetPW(); // proportional wire chamber part
|
PW * pw = anasen->GetPW(); // proportional wire chamber part
|
||||||
|
|
||||||
// output file + tree
|
// output file + trees
|
||||||
TString saveFileName = "SimAnasen1.root";
|
TString saveFileName = "SimAnasen1.root";
|
||||||
printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
|
printf("\e[32m#################################### building Tree in %s\e[0m\n", saveFileName.Data());
|
||||||
TFile * saveFile = new TFile(saveFileName, "recreate");
|
TFile * saveFile = new TFile(saveFileName, "recreate");
|
||||||
TTree * tree = new TTree("tree", "tree");
|
TTree * tree1 = new TTree("tree1", "tree1");
|
||||||
|
TTree * tree2 = new TTree("tree2", "tree2");
|
||||||
|
|
||||||
|
|
||||||
// beam and CM variables saved in tree
|
// beam and CM variables saved in tree
|
||||||
double KEA;
|
double KEA;
|
||||||
tree->Branch("beamKEA", &KEA, "beamKEA/D");
|
double KEA2;
|
||||||
|
tree1->Branch("beamKEA", &KEA, "beamKEA/D");
|
||||||
|
tree2->Branch("beamKEA", &KEA2, "beamKEA/D");
|
||||||
|
|
||||||
double thetaCM, phiCM;
|
double thetaCM, phiCM;
|
||||||
tree->Branch("thetaCM", &thetaCM, "thetaCM/D");
|
double thetaCM2, phiCM2;
|
||||||
tree->Branch("phiCM", &phiCM, "phiCM/D");
|
tree1->Branch("thetaCM", &thetaCM, "thetaCM/D");
|
||||||
|
tree1->Branch("phiCM", &phiCM, "phiCM/D");
|
||||||
|
tree2->Branch("thetaCM", &thetaCM2, "thetaCM/D");
|
||||||
|
tree2->Branch("phiCM", &phiCM2, "phiCM/D");
|
||||||
|
|
||||||
// outgoing particles in lab frame (light/heavy)
|
// outgoing particles in lab frame (light/heavy)
|
||||||
double thetab, phib, Tb;
|
double thetab, phib, Tb;
|
||||||
double thetaB, phiB, TB;
|
double thetaB, phiB, TB;
|
||||||
std::array<double, 2> T;
|
std::array<double, 2> T;
|
||||||
tree->Branch("thetab", &thetab, "thetab/D"); // polar angle of light particle in lab frame
|
tree1->Branch("thetab", &thetab, "thetab/D"); // polar angle of light particle in lab frame
|
||||||
tree->Branch("phib", &phib, "phib/D"); // azimuthal angle of light particle in lab frame
|
tree1->Branch("phib", &phib, "phib/D"); // azimuthal angle of light particle in lab frame
|
||||||
tree->Branch("Tb", &Tb, "Tb/D"); // kinetic energy of light particle at vertex (before energy loss)
|
tree1->Branch("Tb", &Tb, "Tb/D"); // kinetic energy of light particle at vertex (before energy loss)
|
||||||
tree->Branch("thetaB", &thetaB, "thetaB/D");
|
tree1->Branch("thetaB", &thetaB, "thetaB/D");
|
||||||
tree->Branch("phiB", &phiB, "phiB/D");
|
tree1->Branch("phiB", &phiB, "phiB/D");
|
||||||
tree->Branch("TB", &TB, "TB/D"); // kinetic energy of heavy particle at vertex (before energy loss)
|
tree1->Branch("TB", &TB, "TB/D"); // kinetic energy of heavy particle at vertex (before energy loss)
|
||||||
tree->Branch("T", &T, "T/D"); // placeholder for true Q-value, currently set to 0 for simplicity
|
tree1->Branch("T", &T, "T/D"); // placeholder for true Q-value, currently set to 0 for simplicity
|
||||||
|
|
||||||
|
double thetab2, phib2, Tb2;
|
||||||
|
double thetaB2, phiB2, TB2;
|
||||||
|
std::array<double, 2> T2;
|
||||||
|
tree2->Branch("thetab", &thetab2, "thetab/D");
|
||||||
|
tree2->Branch("phib", &phib2, "phib/D");
|
||||||
|
tree2->Branch("Tb", &Tb2, "Tb/D");
|
||||||
|
tree2->Branch("thetaB", &thetaB2, "thetaB/D");
|
||||||
|
tree2->Branch("phiB", &phiB2, "phiB/D");
|
||||||
|
tree2->Branch("TB", &TB2, "TB/D");
|
||||||
|
tree2->Branch("T", &T2, "T/D");
|
||||||
|
|
||||||
// excitation state identifiers
|
// excitation state identifiers
|
||||||
int ExAID;
|
int ExAID;
|
||||||
double ExA;
|
double ExA;
|
||||||
tree->Branch("ExAID", &ExAID, "ExAID/I"); // projectile excitation state ID
|
tree1->Branch("ExAID", &ExAID, "ExAID/I"); // projectile excitation state ID
|
||||||
tree->Branch("ExA", &ExA, "ExA/D"); // projectile excitation energy in MeV
|
tree1->Branch("ExA", &ExA, "ExA/D"); // projectile excitation energy in MeV
|
||||||
|
|
||||||
|
int ExAID2;
|
||||||
|
double ExA2;
|
||||||
|
tree2->Branch("ExAID", &ExAID2, "ExAID/I");
|
||||||
|
tree2->Branch("ExA", &ExA2, "ExA/D");
|
||||||
|
|
||||||
int ExID;
|
int ExID;
|
||||||
double Ex;
|
double Ex;
|
||||||
tree->Branch("ExID", &ExID, "ExID/I"); // target excitation state ID
|
tree1->Branch("ExID", &ExID, "ExID/I"); // target excitation state ID
|
||||||
tree->Branch("Ex", &Ex, "Ex/D"); // target excitation energy in MeV
|
tree1->Branch("Ex", &Ex, "Ex/D"); // target excitation energy in MeV
|
||||||
|
|
||||||
|
int ExID2;
|
||||||
|
double Ex2;
|
||||||
|
tree2->Branch("ExID", &ExID2, "ExID/I");
|
||||||
|
tree2->Branch("Ex", &Ex2, "Ex/D");
|
||||||
|
|
||||||
// true vertex position in target volume
|
// true vertex position in target volume
|
||||||
double vertexX, vertexY, vertexZ;
|
double vertexX, vertexY, vertexZ;
|
||||||
tree->Branch("vX", &vertexX, "VertexX/D"); // true vertex X position in mm
|
tree1->Branch("vX", &vertexX, "VertexX/D"); // true vertex X position in mm
|
||||||
tree->Branch("vY", &vertexY, "VertexY/D"); // true vertex Y position in mm
|
tree1->Branch("vY", &vertexY, "VertexY/D"); // true vertex Y position in mm
|
||||||
tree->Branch("vZ", &vertexZ, "VertexZ/D"); // true vertex Z position in mm
|
tree1->Branch("vZ", &vertexZ, "VertexZ/D"); // true vertex Z position in mm
|
||||||
|
|
||||||
|
double vertexX2, vertexY2, vertexZ2;
|
||||||
|
tree2->Branch("vX", &vertexX2, "VertexX/D");
|
||||||
|
tree2->Branch("vY", &vertexY2, "VertexY/D");
|
||||||
|
tree2->Branch("vZ", &vertexZ2, "VertexZ/D");
|
||||||
|
|
||||||
// reconstructed SX3 hit position
|
// reconstructed SX3 hit position
|
||||||
double sx3X, sx3Y, sx3Z;
|
double sx3X, sx3Y, sx3Z;
|
||||||
tree->Branch("sx3X", &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing)
|
tree1->Branch("sx3X", &sx3X, "sx3X/D"); // reconstructed X position from SX3 (with optional smearing)
|
||||||
tree->Branch("sx3Y", &sx3Y, "sx3Y/D"); // reconstructed Y position from SX3 (with optional smearing)
|
tree1->Branch("sx3Y", &sx3Y, "sx3Y/D"); // reconstructed Y position from SX3 (with optional smearing)
|
||||||
tree->Branch("sx3Z", &sx3Z, "sx3Z/D"); // reconstructed Z position from SX3 (with optional smearing)
|
tree1->Branch("sx3Z", &sx3Z, "sx3Z/D"); // reconstructed Z position from SX3 (with optional smearing)
|
||||||
|
|
||||||
|
double sx3X2, sx3Y2, sx3Z2;
|
||||||
|
tree2->Branch("sx3X", &sx3X2, "sx3X/D");
|
||||||
|
tree2->Branch("sx3Y", &sx3Y2, "sx3Y/D");
|
||||||
|
tree2->Branch("sx3Z", &sx3Z2, "sx3Z/D");
|
||||||
|
|
||||||
// PW nearest and next nearest wires
|
// PW nearest and next nearest wires
|
||||||
int anodeID[2], cathodeID[2];
|
int anodeID[2], cathodeID[2];
|
||||||
tree->Branch("aID", anodeID, "anodeID/I"); // anodeID[0] is nearest anode wire, anodeID[1] is next nearest anode wire
|
int anodeID2[2], cathodeID2[2];
|
||||||
tree->Branch("cID", cathodeID, "cathodeID/I"); // cathodeID[0] is nearest cathode wire, cathodeID[1] is next nearest cathode wire
|
tree1->Branch("aID", anodeID, "anodeID/I"); // anodeID[0] is nearest anode wire, anodeID[1] is next nearest anode wire
|
||||||
|
tree1->Branch("cID", cathodeID, "cathodeID/I"); // cathodeID[0] is nearest cathode wire, cathodeID[1] is next nearest cathode wire
|
||||||
|
tree2->Branch("aID", anodeID2, "anodeID/I");
|
||||||
|
tree2->Branch("cID", cathodeID2, "cathodeID/I");
|
||||||
|
|
||||||
// distances to nearest wires
|
// distances to nearest wires
|
||||||
double anodeDist[2], cathodeDist[2];
|
double anodeDist[2], cathodeDist[2];
|
||||||
tree->Branch("aDist", anodeDist, "anodeDist/D");
|
double anodeDist2[2], cathodeDist2[2];
|
||||||
tree->Branch("cDist", cathodeDist, "cathodeDist/D");
|
tree1->Branch("aDist", anodeDist, "anodeDist/D");
|
||||||
|
tree1->Branch("cDist", cathodeDist, "cathodeDist/D");
|
||||||
|
tree2->Branch("aDist", anodeDist2, "anodeDist/D");
|
||||||
|
tree2->Branch("cDist", cathodeDist2, "cathodeDist/D");
|
||||||
|
|
||||||
// SX3 channel assignment and Z fraction (depth) information
|
// SX3 channel assignment and Z fraction (depth) information
|
||||||
int sx3ID, sx3Up, sx3Dn, sx3Bk;
|
int sx3ID, sx3Up, sx3Dn, sx3Bk;
|
||||||
double sx3ZFrac;
|
double sx3ZFrac;
|
||||||
tree->Branch("sx3ID", &sx3ID, "sx3ID/I");
|
int sx3ID2, sx3Up2, sx3Dn2, sx3Bk2;
|
||||||
tree->Branch("sx3Up", &sx3Up, "sx3Up/I");
|
double sx3ZFrac2;
|
||||||
tree->Branch("sx3Dn", &sx3Dn, "sx3Dn/I");
|
tree1->Branch("sx3ID", &sx3ID, "sx3ID/I");
|
||||||
tree->Branch("sx3Bk", &sx3Bk, "sx3Bk/I");
|
tree1->Branch("sx3Up", &sx3Up, "sx3Up/I");
|
||||||
tree->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D");
|
tree1->Branch("sx3Dn", &sx3Dn, "sx3Dn/I");
|
||||||
|
tree1->Branch("sx3Bk", &sx3Bk, "sx3Bk/I");
|
||||||
|
tree1->Branch("sx3ZFrac", &sx3ZFrac, "sx3ZFrac/D");
|
||||||
|
tree2->Branch("sx3ID", &sx3ID2, "sx3ID/I");
|
||||||
|
tree2->Branch("sx3Up", &sx3Up2, "sx3Up/I");
|
||||||
|
tree2->Branch("sx3Dn", &sx3Dn2, "sx3Dn/I");
|
||||||
|
tree2->Branch("sx3Bk", &sx3Bk2, "sx3Bk/I");
|
||||||
|
tree2->Branch("sx3ZFrac", &sx3ZFrac2, "sx3ZFrac/D");
|
||||||
|
|
||||||
// reconstructed angles from PW track fit, method 1 and 2
|
// reconstructed angles from PW track fit, method 1 and 2
|
||||||
double reTheta, rePhi;
|
double reTheta, rePhi;
|
||||||
tree->Branch("reTheta", &reTheta, "reconstucted_theta/D");
|
double reTheta2, rePhi2;
|
||||||
tree->Branch("rePhi", &rePhi, "reconstucted_phi/D");
|
tree1->Branch("reTheta", &reTheta, "reconstucted_theta/D");
|
||||||
|
tree1->Branch("rePhi", &rePhi, "reconstucted_phi/D");
|
||||||
|
tree2->Branch("reTheta", &reTheta2, "reconstucted_theta/D");
|
||||||
|
tree2->Branch("rePhi", &rePhi2, "reconstucted_phi/D");
|
||||||
|
|
||||||
double reTheta1, rePhi1;
|
double reTheta1, rePhi1;
|
||||||
tree->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D");
|
double reTheta12, rePhi12;
|
||||||
tree->Branch("rePhi1", &rePhi1, "reconstucted_phi1/D");
|
tree1->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D");
|
||||||
|
tree1->Branch("rePhi1", &rePhi1, "reconstucted_phi1/D");
|
||||||
|
tree2->Branch("reTheta1", &reTheta12, "reconstucted_theta1/D");
|
||||||
|
tree2->Branch("rePhi1", &rePhi12, "reconstucted_phi1/D");
|
||||||
|
|
||||||
// reconstructed vertex Z from PW fit
|
// reconstructed vertex Z from PW fit
|
||||||
double z0;
|
double z0;
|
||||||
tree->Branch("z0", &z0, "reconstucted_Z/D");
|
double z02;
|
||||||
|
tree1->Branch("z0", &z0, "reconstucted_Z/D");
|
||||||
|
tree2->Branch("z0", &z02, "reconstucted_Z/D");
|
||||||
|
|
||||||
//========timer
|
//========timer
|
||||||
TBenchmark clock;
|
TBenchmark clock;
|
||||||
|
|
@ -222,25 +327,46 @@ int main(int argc, char **argv){
|
||||||
thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ;
|
thetaCM = TMath::ACos(2 * gRandom->Rndm() - 1) ;
|
||||||
phiCM = (gRandom->Rndm() - 0.5) * TMath::TwoPi();
|
phiCM = (gRandom->Rndm() - 0.5) * TMath::TwoPi();
|
||||||
|
|
||||||
//==== Calculate reaction kinematics in lab frame
|
//==== Calculate reaction kinematics in lab frame for the primary transfer
|
||||||
TLorentzVector * output = transfer.Event(thetaCM, phiCM); // returns array of outputs
|
TLorentzVector * output = transfer.Event(thetaCM, phiCM); // returns array of outputs
|
||||||
TLorentzVector Pb = output[2]; // light particle or product A
|
TLorentzVector Pb = output[2]; // primary proton from transfer
|
||||||
TLorentzVector PB = output[3]; // heavy particle or product B
|
TLorentzVector PB = output[3]; // excited 21Na* heavy product
|
||||||
|
|
||||||
thetab = Pb.Theta() * TMath::RadToDeg();
|
thetab = Pb.Theta() * TMath::RadToDeg();
|
||||||
|
Tb = (Pb.E() - Pb.M()); // kinetic energy of the light proton from the primary transfer
|
||||||
thetaB = PB.Theta() * TMath::RadToDeg();
|
thetaB = PB.Theta() * TMath::RadToDeg();
|
||||||
|
|
||||||
Tb = (Pb.E() - Pb.M()); // kinetic energy of light particle at vertex (before energy loss) units of MeV
|
|
||||||
TB = (PB.E() - PB.M());
|
TB = (PB.E() - PB.M());
|
||||||
T[0] = Tb;
|
|
||||||
T[1] = TB;
|
|
||||||
//if (Tb < 1.5) {
|
|
||||||
// //skip event if light particle energy after loss is below detection threshold of 1.5 MeV
|
|
||||||
// continue;
|
|
||||||
//}
|
|
||||||
|
|
||||||
phib = Pb.Phi() * TMath::RadToDeg();
|
phib = Pb.Phi() * TMath::RadToDeg();
|
||||||
phiB = PB.Phi() * TMath::RadToDeg();
|
phiB = PB.Phi() * TMath::RadToDeg();
|
||||||
|
T[0] = Tb;
|
||||||
|
T[1] = TB;
|
||||||
|
|
||||||
|
// prepare secondary proton from 21Na* sequential decay
|
||||||
|
TLorentzVector decayProton;
|
||||||
|
TLorentzVector heavy20;
|
||||||
|
if(enableSequentialDecay){
|
||||||
|
heavy20 = SimulateSequentialDecay(PB, decayDaughterA, decayDaughterZ,
|
||||||
|
decayEjectA, decayEjectZ, decayProton);
|
||||||
|
thetab2 = decayProton.Theta() * TMath::RadToDeg();
|
||||||
|
phib2 = decayProton.Phi() * TMath::RadToDeg();
|
||||||
|
Tb2 = decayProton.E() - decayProton.M();
|
||||||
|
thetaB2 = heavy20.Theta() * TMath::RadToDeg();
|
||||||
|
phiB2 = heavy20.Phi() * TMath::RadToDeg();
|
||||||
|
TB2 = heavy20.E() - heavy20.M();
|
||||||
|
T2[0] = Tb2;
|
||||||
|
T2[1] = TB2;
|
||||||
|
} else {
|
||||||
|
thetab2 = TMath::QuietNaN();
|
||||||
|
phib2 = TMath::QuietNaN();
|
||||||
|
Tb2 = TMath::QuietNaN();
|
||||||
|
thetaB2 = TMath::QuietNaN();
|
||||||
|
phiB2 = TMath::QuietNaN();
|
||||||
|
TB2 = TMath::QuietNaN();
|
||||||
|
T2[0] = TMath::QuietNaN();
|
||||||
|
T2[1] = TMath::QuietNaN();
|
||||||
|
}
|
||||||
|
|
||||||
|
delete [] output;
|
||||||
|
|
||||||
// vertex position in target volume
|
// vertex position in target volume
|
||||||
vertexX = (vertexXRange[1]- vertexXRange[0])*gRandom->Rndm() + vertexXRange[0];
|
vertexX = (vertexXRange[1]- vertexXRange[0])*gRandom->Rndm() + vertexXRange[0];
|
||||||
|
|
@ -311,7 +437,79 @@ int main(int argc, char **argv){
|
||||||
rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg();
|
rePhi1 = pw->GetTrackPhi() * TMath::RadToDeg();
|
||||||
|
|
||||||
z0 = pw->GetZ0();
|
z0 = pw->GetZ0();
|
||||||
tree->Fill();
|
tree1->Fill();
|
||||||
|
|
||||||
|
// fill tree2 using the secondary proton from 21Na* decay
|
||||||
|
TVector3 dir2(1, 0, 0);
|
||||||
|
dir2.SetTheta(thetab2 * TMath::DegToRad());
|
||||||
|
dir2.SetPhi(phib2 * TMath::DegToRad());
|
||||||
|
|
||||||
|
pw->FindWireID(vertex, dir2, false);
|
||||||
|
sx3->FindSX3Pos(vertex, dir2, false);
|
||||||
|
PWHitInfo hitInfo2 = pw->GetHitInfo();
|
||||||
|
|
||||||
|
anodeID2[0] = hitInfo2.nearestWire.first;
|
||||||
|
cathodeID2[0] = hitInfo2.nearestWire.second;
|
||||||
|
anodeID2[1] = hitInfo2.nextNearestWire.first;
|
||||||
|
cathodeID2[1] = hitInfo2.nextNearestWire.second;
|
||||||
|
|
||||||
|
anodeDist2[1] = hitInfo2.nextNearestDist.first;
|
||||||
|
cathodeDist2[1] = hitInfo2.nextNearestDist.second;
|
||||||
|
|
||||||
|
if(IsDeadAnode(anodeID2[0]) || IsDeadCathode(cathodeID2[0])){
|
||||||
|
sx3ID2 = -1;
|
||||||
|
} else {
|
||||||
|
sx3ID2 = sx3->GetID();
|
||||||
|
}
|
||||||
|
|
||||||
|
if(sx3ID2 < 0 || IsDeadSX3(sx3ID2)){
|
||||||
|
sx3ID2 = -1;
|
||||||
|
sx3Up2 = -1;
|
||||||
|
sx3Dn2 = -1;
|
||||||
|
sx3Bk2 = -1;
|
||||||
|
sx3ZFrac2 = TMath::QuietNaN();
|
||||||
|
sx3X2 = TMath::QuietNaN();
|
||||||
|
sx3Y2 = TMath::QuietNaN();
|
||||||
|
sx3Z2 = TMath::QuietNaN();
|
||||||
|
anodeDist2[0] = TMath::QuietNaN();
|
||||||
|
cathodeDist2[0] = TMath::QuietNaN();
|
||||||
|
reTheta2 = TMath::QuietNaN();
|
||||||
|
rePhi2 = TMath::QuietNaN();
|
||||||
|
reTheta12 = TMath::QuietNaN();
|
||||||
|
rePhi12 = TMath::QuietNaN();
|
||||||
|
z02 = TMath::QuietNaN();
|
||||||
|
} else {
|
||||||
|
anodeDist2[0] = hitInfo2.nearestDist.first;
|
||||||
|
cathodeDist2[0] = hitInfo2.nearestDist.second;
|
||||||
|
sx3Up2 = sx3->GetChUp();
|
||||||
|
sx3Dn2 = sx3->GetChDn();
|
||||||
|
sx3Bk2 = sx3->GetChBk();
|
||||||
|
sx3ZFrac2 = sx3->GetZFrac();
|
||||||
|
TVector3 hitPos2 = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L);
|
||||||
|
sx3X2 = hitPos2.X();
|
||||||
|
sx3Y2 = hitPos2.Y();
|
||||||
|
sx3Z2 = hitPos2.Z();
|
||||||
|
pw->CalTrack(hitPos2, anodeID2[0], cathodeID2[0], false);
|
||||||
|
reTheta2 = pw->GetTrackTheta() * TMath::RadToDeg();
|
||||||
|
rePhi2 = pw->GetTrackPhi() * TMath::RadToDeg();
|
||||||
|
pw->CalTrack2(hitPos2, hitInfo2, sigmaPW_A, sigmaPW_C, false);
|
||||||
|
reTheta12 = pw->GetTrackTheta() * TMath::RadToDeg();
|
||||||
|
rePhi12 = pw->GetTrackPhi() * TMath::RadToDeg();
|
||||||
|
z02 = pw->GetZ0();
|
||||||
|
}
|
||||||
|
|
||||||
|
KEA2 = KEA;
|
||||||
|
thetaCM2 = thetaCM;
|
||||||
|
phiCM2 = phiCM;
|
||||||
|
ExAID2 = ExAID;
|
||||||
|
ExA2 = ExA;
|
||||||
|
ExID2 = ExID;
|
||||||
|
Ex2 = Ex;
|
||||||
|
vertexX2 = vertexX;
|
||||||
|
vertexY2 = vertexY;
|
||||||
|
vertexZ2 = vertexZ;
|
||||||
|
|
||||||
|
tree2->Fill();
|
||||||
|
|
||||||
}else{
|
}else{
|
||||||
// no valid SX3 hit: mark clearly invalid
|
// no valid SX3 hit: mark clearly invalid
|
||||||
|
|
@ -332,7 +530,7 @@ int main(int argc, char **argv){
|
||||||
//Tb = -12354567; // mark kinetic energy as invalid for no hit case
|
//Tb = -12354567; // mark kinetic energy as invalid for no hit case
|
||||||
// fill tree with original data (no energy loss for these events)
|
// fill tree with original data (no energy loss for these events)
|
||||||
//comment out tree fill for no hit case
|
//comment out tree fill for no hit case
|
||||||
//tree->Fill();
|
//tree1->Fill();
|
||||||
}
|
}
|
||||||
|
|
||||||
//#################################################################### Timer
|
//#################################################################### Timer
|
||||||
|
|
@ -355,12 +553,13 @@ int main(int argc, char **argv){
|
||||||
}
|
}
|
||||||
|
|
||||||
// write results to ROOT file and close
|
// write results to ROOT file and close
|
||||||
//tree->Write();
|
tree1->Write("", TObject::kOverwrite);
|
||||||
tree->Write("", TObject::kOverwrite);
|
tree2->Write("", TObject::kOverwrite);
|
||||||
int count = tree->GetEntries();
|
int count1 = tree1->GetEntries();
|
||||||
|
int count2 = tree2->GetEntries();
|
||||||
saveFile->Close();
|
saveFile->Close();
|
||||||
|
|
||||||
printf("=============== done. saved as %s. tree entries: %d\n", saveFileName.Data(), count);
|
printf("=============== done. saved as %s. tree1 entries: %d, tree2 entries: %d\n", saveFileName.Data(), count1, count2);
|
||||||
|
|
||||||
if(enableVis){ // to enable visualization, run with 3rd argument "vis", e.g. "./anasenMC 1000 vis"
|
if(enableVis){ // to enable visualization, run with 3rd argument "vis", e.g. "./anasenMC 1000 vis"
|
||||||
printf("Displaying geometry with %zu tracks from simulation\n", visTrackVertex.size());
|
printf("Displaying geometry with %zu tracks from simulation\n", visTrackVertex.size());
|
||||||
|
|
|
||||||
|
|
@ -1,795 +0,0 @@
|
||||||
#include "ClassData.h"
|
|
||||||
#include "Hit.h"
|
|
||||||
#include <algorithm>
|
|
||||||
#include <filesystem>
|
|
||||||
|
|
||||||
// #include "AggSeparator.h"
|
|
||||||
|
|
||||||
class FSUReader{
|
|
||||||
|
|
||||||
public:
|
|
||||||
FSUReader();
|
|
||||||
FSUReader(std::string fileName, uInt dataSize = 100, int verbose = 1);
|
|
||||||
FSUReader(std::vector<std::string> fileList, uInt dataSize = 100, int verbose = 1);
|
|
||||||
~FSUReader();
|
|
||||||
|
|
||||||
void OpenFile(std::string fileName, uInt dataSize, int verbose = 1);
|
|
||||||
bool IsOpen() const{return inFile == nullptr ? false : true;}
|
|
||||||
bool IsEndOfFile() const {
|
|
||||||
// printf("%s : %d | %ld |%ld\n", __func__, feof(inFile), ftell(inFile), inFileSize);
|
|
||||||
if(fileList.empty() ) {
|
|
||||||
if( (uLong )ftell(inFile) >= inFileSize){
|
|
||||||
return true;
|
|
||||||
}else{
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
}else{
|
|
||||||
if( fileID + 1 == (int) fileList.size() && ((uLong)ftell(inFile) >= inFileSize) ) {
|
|
||||||
return true;
|
|
||||||
}else{
|
|
||||||
return false;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void ScanNumBlock(int verbose = 1, uShort saveData = 0); // saveData = 0 (no save), 1 (no trace), 2 (with trace);
|
|
||||||
int ReadNextBlock(bool traceON = false, int verbose = 0, uShort saveData = 0); // saveData = 0 (no save), 1 (no trace), 2 (with trace);
|
|
||||||
int ReadBlock(unsigned int ID, int verbose = 0);
|
|
||||||
|
|
||||||
unsigned int GetFilePos() const {return filePos;}
|
|
||||||
unsigned long GetTotNumBlock() const{ return totNumBlock;}
|
|
||||||
std::vector<unsigned int> GetBlockTimestamp() const {return blockTimeStamp;}
|
|
||||||
|
|
||||||
Data * GetData() const{return data;}
|
|
||||||
|
|
||||||
std::string GetFileName() const{return fileName;}
|
|
||||||
int GetDPPType() const{return DPPType;}
|
|
||||||
int GetSN() const{return sn;}
|
|
||||||
int GetTick2ns() const{return tick2ns;}
|
|
||||||
int GetNumCh() const{return numCh;}
|
|
||||||
int GetFileOrder() const{return order;}
|
|
||||||
int GetChMask() const{return chMask;}
|
|
||||||
unsigned long GetFileByteSize() const {return inFileSize;}
|
|
||||||
|
|
||||||
void ClearHitList() { hit.clear();}
|
|
||||||
ulong GetHitListLength() const {return hit.size();}
|
|
||||||
std::vector<Hit> GetHitVector() const {return hit;}
|
|
||||||
void SortHit(int verbose = false);
|
|
||||||
Hit GetHit(int id) const {
|
|
||||||
if( id < 0 ) id = hit.size() + id;
|
|
||||||
return hit[id];
|
|
||||||
}
|
|
||||||
|
|
||||||
void ClearHitCount() {hitCount = 0;}
|
|
||||||
ulong GetHitCount() const{return hitCount;}
|
|
||||||
|
|
||||||
std::vector<Hit> ReadBatch(unsigned int batchSize = 1000000, bool verbose = false); // output the sorted Hit
|
|
||||||
|
|
||||||
// std::string SaveHit(std::vector<Hit> hitList, bool isAppend = false);
|
|
||||||
// std::string SaveHit2NewFile(std::string saveFolder = "./", std::string indexStr = "");
|
|
||||||
// void SortAndSaveTS(unsigned int batchSize = 1000000, bool verbose = false);
|
|
||||||
// off_t GetTSFileSize() const {return tsFileSize;}
|
|
||||||
|
|
||||||
//TODO
|
|
||||||
//void SplitFile(unsigned long hitSizePreFile);
|
|
||||||
|
|
||||||
void PrintHit(ulong numHit = -1, ulong startIndex = 0) {
|
|
||||||
for( ulong i = startIndex; i < std::min(numHit, hitCount); i++){
|
|
||||||
printf("%10zu ", i); hit[i].Print();
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
static void PrintHitListInfo(std::vector<Hit> * hitList, std::string name){
|
|
||||||
size_t n = hitList->size();
|
|
||||||
size_t s = sizeof(Hit);
|
|
||||||
printf("============== %s, size : %zu | %.2f MByte\n", name.c_str(), n, n*s/1024./1024.);
|
|
||||||
if( n > 0 ){
|
|
||||||
printf("t0 : %15llu ns\n", hitList->front().timestamp);
|
|
||||||
printf("t1 : %15llu ns\n", hitList->back().timestamp);
|
|
||||||
printf("dt : %15.3f ms\n", (hitList->back().timestamp - hitList->front().timestamp)/1e6);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
void PrintHitListInfo(){
|
|
||||||
size_t n = hit.size();
|
|
||||||
size_t s = sizeof(Hit);
|
|
||||||
printf("============== reader.hit, size : %zu | %.2f MByte\n", n, n*s/1024./1024.);
|
|
||||||
if( n > 0 ){
|
|
||||||
printf("t0 : %15llu ns\n", hit.at(0).timestamp);
|
|
||||||
printf("t1 : %15llu ns\n", hit.back().timestamp);
|
|
||||||
printf("dt : %15.3f ms\n", (hit.back().timestamp - hit.front().timestamp)/1e6);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
//void SaveAsCAENCoMPASSFormat();
|
|
||||||
|
|
||||||
private:
|
|
||||||
|
|
||||||
FILE * inFile;
|
|
||||||
Data * data;
|
|
||||||
|
|
||||||
std::string fileName;
|
|
||||||
std::vector<std::string> fileList;
|
|
||||||
short fileID;
|
|
||||||
unsigned long inFileSize;
|
|
||||||
unsigned int filePos;
|
|
||||||
unsigned long totNumBlock;
|
|
||||||
unsigned int blockID;
|
|
||||||
|
|
||||||
bool isDualBlock;
|
|
||||||
|
|
||||||
uShort sn;
|
|
||||||
uShort DPPType;
|
|
||||||
uShort tick2ns;
|
|
||||||
uShort order;
|
|
||||||
uShort chMask;
|
|
||||||
uShort numCh;
|
|
||||||
|
|
||||||
std::vector<unsigned int> blockPos;
|
|
||||||
std::vector<unsigned int > blockTimeStamp;
|
|
||||||
|
|
||||||
unsigned long hitCount;
|
|
||||||
|
|
||||||
std::vector<Hit> hit;
|
|
||||||
|
|
||||||
unsigned int word[1]; /// 4 byte
|
|
||||||
size_t dummy;
|
|
||||||
char * buffer;
|
|
||||||
|
|
||||||
off_t tsFileSize;
|
|
||||||
|
|
||||||
};
|
|
||||||
|
|
||||||
inline FSUReader::~FSUReader(){
|
|
||||||
delete data;
|
|
||||||
|
|
||||||
if( inFile ) fclose(inFile);
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
inline FSUReader::FSUReader(){
|
|
||||||
inFile = nullptr;
|
|
||||||
data = nullptr;
|
|
||||||
|
|
||||||
blockPos.clear();
|
|
||||||
blockTimeStamp.clear();
|
|
||||||
hit.clear();
|
|
||||||
|
|
||||||
fileList.clear();
|
|
||||||
fileID = -1;
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
inline FSUReader::FSUReader(std::string fileName, uInt dataSize, int verbose){
|
|
||||||
inFile = nullptr;
|
|
||||||
data = nullptr;
|
|
||||||
|
|
||||||
blockPos.clear();
|
|
||||||
blockTimeStamp.clear();
|
|
||||||
hit.clear();
|
|
||||||
|
|
||||||
fileList.clear();
|
|
||||||
fileID = -1;
|
|
||||||
OpenFile(fileName, dataSize, verbose);
|
|
||||||
}
|
|
||||||
|
|
||||||
inline FSUReader::FSUReader(std::vector<std::string> fileList, uInt dataSize, int verbose){
|
|
||||||
inFile = nullptr;
|
|
||||||
data = nullptr;
|
|
||||||
|
|
||||||
blockPos.clear();
|
|
||||||
blockTimeStamp.clear();
|
|
||||||
hit.clear();
|
|
||||||
//The files are the same DPPType and sn
|
|
||||||
this->fileList = fileList;
|
|
||||||
fileID = 0;
|
|
||||||
OpenFile(fileList[fileID], dataSize, verbose);
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
inline void FSUReader::OpenFile(std::string fileName, uInt dataSize, int verbose){
|
|
||||||
|
|
||||||
/// File format must be YYY...Y_runXXX_AAA_BBB_TT_CCC.fsu
|
|
||||||
/// YYY...Y = prefix
|
|
||||||
/// XXX = runID, 3 digits
|
|
||||||
/// AAA = board Serial Number, 3 digits
|
|
||||||
/// BBB = DPPtype, 3 digits
|
|
||||||
/// TT = tick2ns, any digits
|
|
||||||
/// CCC = over size index, 3 digits
|
|
||||||
|
|
||||||
if( inFile != nullptr ) fclose(inFile);
|
|
||||||
|
|
||||||
inFile = fopen(fileName.c_str(), "r");
|
|
||||||
|
|
||||||
if( inFile == NULL ){
|
|
||||||
printf("FSUReader::Cannot open file : %s \n", fileName.c_str());
|
|
||||||
this->fileName = "";
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
this->fileName = fileName;
|
|
||||||
|
|
||||||
fseek(inFile, 0L, SEEK_END);
|
|
||||||
inFileSize = ftell(inFile);
|
|
||||||
if(verbose) printf("%s | file size : %ld Byte = %.2f MB\n", fileName.c_str() , inFileSize, inFileSize/1024./1024.);
|
|
||||||
fseek(inFile, 0L, SEEK_SET);
|
|
||||||
filePos = 0;
|
|
||||||
|
|
||||||
if( fileID > 0 ) return;
|
|
||||||
|
|
||||||
totNumBlock = 0;
|
|
||||||
blockID = 0;
|
|
||||||
blockPos.clear();
|
|
||||||
blockTimeStamp.clear();
|
|
||||||
|
|
||||||
hitCount = 0;
|
|
||||||
hit.clear();
|
|
||||||
|
|
||||||
//check is the file is *.fsu or *.fsu.X
|
|
||||||
size_t found = fileName.find_last_of('.');
|
|
||||||
std::string ext = fileName.substr(found + 1);
|
|
||||||
|
|
||||||
if( ext.find("fsu") != std::string::npos ) {
|
|
||||||
if(verbose > 1) printf("It is an raw data *.fsu format\n");
|
|
||||||
isDualBlock = false;
|
|
||||||
chMask = -1;
|
|
||||||
}else{
|
|
||||||
chMask = atoi(ext.c_str());
|
|
||||||
isDualBlock = true;
|
|
||||||
if(verbose > 1) printf("It is a splitted dual block data *.fsu.X format, dual channel mask : %d \n", chMask);
|
|
||||||
}
|
|
||||||
|
|
||||||
std::string fileNameNoExt;
|
|
||||||
found = fileName.find_last_of(".fsu");
|
|
||||||
size_t found2 = fileName.find_last_of('/');
|
|
||||||
if( found2 == std::string::npos ){
|
|
||||||
fileNameNoExt = fileName.substr(0, found-4);
|
|
||||||
}else{
|
|
||||||
fileNameNoExt = fileName.substr(found2+1, found-4);
|
|
||||||
}
|
|
||||||
|
|
||||||
// Split the string by underscores
|
|
||||||
std::istringstream iss(fileNameNoExt);
|
|
||||||
std::vector<std::string> tokens;
|
|
||||||
std::string token;
|
|
||||||
|
|
||||||
while (std::getline(iss, token, '_')) { tokens.push_back(token); }
|
|
||||||
sn = atoi(tokens[2].c_str());
|
|
||||||
tick2ns = atoi(tokens[4].c_str());
|
|
||||||
order = atoi(tokens[5].c_str());
|
|
||||||
|
|
||||||
DPPType = 0;
|
|
||||||
if( fileName.find("PHA") != std::string::npos ) DPPType = DPPTypeCode::DPP_PHA_CODE;
|
|
||||||
if( fileName.find("PSD") != std::string::npos ) DPPType = DPPTypeCode::DPP_PSD_CODE;
|
|
||||||
if( fileName.find("QDC") != std::string::npos ) DPPType = DPPTypeCode::DPP_QDC_CODE;
|
|
||||||
if( DPPType == 0 ){
|
|
||||||
fclose(inFile);
|
|
||||||
inFile = nullptr;
|
|
||||||
printf("Cannot find DPPType in the filename. Abort.");
|
|
||||||
return ;
|
|
||||||
}
|
|
||||||
|
|
||||||
numCh = (DPPType == DPPTypeCode::DPP_QDC_CODE ? 64 : 16);
|
|
||||||
|
|
||||||
data = new Data(numCh, dataSize);
|
|
||||||
data->tick2ns = tick2ns;
|
|
||||||
data->boardSN = sn;
|
|
||||||
data->DPPType = DPPType;
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
inline int FSUReader::ReadNextBlock(bool traceON, int verbose, uShort saveData){
|
|
||||||
if( inFile == NULL ) return -1;
|
|
||||||
if( feof(inFile) || filePos >= inFileSize) {
|
|
||||||
if( fileID >= 0 && fileID + 1 < (short) fileList.size() ){
|
|
||||||
printf("-------------- next file\n");
|
|
||||||
fileID ++;
|
|
||||||
OpenFile(fileList[fileID], data->GetDataSize(), 1 );
|
|
||||||
}else{
|
|
||||||
return -1;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
dummy = fread(word, 4, 1, inFile);
|
|
||||||
fseek(inFile, -4, SEEK_CUR);
|
|
||||||
|
|
||||||
if( dummy != 1) {
|
|
||||||
printf("fread error, should read 4 bytes, but read %ld x 4 byte, file pos: %ld / %ld byte\n",
|
|
||||||
dummy, ftell(inFile), inFileSize);
|
|
||||||
return -10;
|
|
||||||
}
|
|
||||||
short header = ((word[0] >> 28 ) & 0xF);
|
|
||||||
|
|
||||||
Hit temp;
|
|
||||||
|
|
||||||
if( header == 0xA ) { ///normal header
|
|
||||||
|
|
||||||
unsigned int aggSize = (word[0] & 0x0FFFFFFF) * 4; ///byte
|
|
||||||
if( aggSize > inFileSize - ftell(inFile)) aggSize = inFileSize - ftell(inFile);
|
|
||||||
buffer = new char[aggSize];
|
|
||||||
dummy = fread(buffer, aggSize, 1, inFile);
|
|
||||||
filePos = ftell(inFile);
|
|
||||||
if( dummy != 1) {
|
|
||||||
printf("fread error, should read %d bytes, but read %ld x %d byte, file pos: %ld / %ld byte \n",
|
|
||||||
aggSize, dummy, aggSize, ftell(inFile), inFileSize);
|
|
||||||
return -30;
|
|
||||||
}
|
|
||||||
|
|
||||||
data->DecodeBuffer(buffer, aggSize, !traceON, verbose); // data will own the buffer
|
|
||||||
//printf(" word Index = %u | filePos : %u | ", data->GetWordIndex(), filePos);
|
|
||||||
|
|
||||||
}else if( (header & 0xF ) == 0x8 ) { /// dual channel header
|
|
||||||
|
|
||||||
unsigned int dualSize = (word[0] & 0x7FFFFFFF) * 4; ///byte
|
|
||||||
buffer = new char[dualSize];
|
|
||||||
dummy = fread(buffer, dualSize, 1, inFile);
|
|
||||||
filePos = ftell(inFile);
|
|
||||||
|
|
||||||
data->buffer = buffer;
|
|
||||||
data->DecodeDualBlock(buffer, dualSize, DPPType, chMask, !traceON, verbose);
|
|
||||||
|
|
||||||
}else{
|
|
||||||
printf("incorrect header.\n trminate.");
|
|
||||||
return -20;
|
|
||||||
}
|
|
||||||
|
|
||||||
unsigned int eventCout = 0;
|
|
||||||
|
|
||||||
for( int ch = 0; ch < data->GetNChannel(); ch++){
|
|
||||||
if( data->NumEventsDecoded[ch] == 0 ) continue;
|
|
||||||
|
|
||||||
hitCount += data->NumEventsDecoded[ch];
|
|
||||||
eventCout += data->NumEventsDecoded[ch];
|
|
||||||
|
|
||||||
if( saveData ){
|
|
||||||
int start = data->GetDataIndex(ch) - data->NumEventsDecoded[ch] + 1;
|
|
||||||
if( start < 0 ) start = start + data->GetDataSize();
|
|
||||||
|
|
||||||
for( int i = start; i < start + data->NumEventsDecoded[ch]; i++ ){
|
|
||||||
int k = i % data->GetDataSize();
|
|
||||||
|
|
||||||
temp.sn = sn;
|
|
||||||
temp.ch = ch;
|
|
||||||
temp.energy = data->GetEnergy(ch, k);
|
|
||||||
temp.energy2 = data->GetEnergy2(ch, k);
|
|
||||||
temp.timestamp = data->GetTimestamp(ch, k);
|
|
||||||
temp.fineTime = data->GetFineTime(ch, k);
|
|
||||||
temp.pileUp = data->GetPileUp(ch, k);
|
|
||||||
if( saveData > 1 ) {
|
|
||||||
temp.traceLength = data->Waveform1[ch][k].size();
|
|
||||||
temp.trace = data->Waveform1[ch][k];
|
|
||||||
}else{
|
|
||||||
temp.traceLength = 0;
|
|
||||||
if( temp.trace.size() > 0 ) temp.trace.clear();
|
|
||||||
}
|
|
||||||
|
|
||||||
hit.push_back(temp);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
data->ClearTriggerRate();
|
|
||||||
data->ClearNumEventsDecoded();
|
|
||||||
data->ClearBuffer(); // this will clear the buffer.
|
|
||||||
return 0;
|
|
||||||
}
|
|
||||||
|
|
||||||
inline int FSUReader::ReadBlock(unsigned int ID, int verbose){
|
|
||||||
if( totNumBlock == 0 )return -1;
|
|
||||||
if( ID >= totNumBlock )return -1;
|
|
||||||
|
|
||||||
data->ClearData();
|
|
||||||
|
|
||||||
fseek( inFile, 0L, SEEK_SET);
|
|
||||||
|
|
||||||
if( verbose ) printf("Block index: %u, File Pos: %u byte\n", ID, blockPos[ID]);
|
|
||||||
|
|
||||||
fseek(inFile, blockPos[ID], SEEK_CUR);
|
|
||||||
filePos = blockPos[ID];
|
|
||||||
blockID = ID;
|
|
||||||
return ReadNextBlock(false, verbose, false);
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
inline void FSUReader::SortHit(int verbose){
|
|
||||||
if( verbose) printf("\nQuick Sort hit array according to time...");
|
|
||||||
std::sort(hit.begin(), hit.end(), [](const Hit& a, const Hit& b) {
|
|
||||||
return a.timestamp < b.timestamp;
|
|
||||||
});
|
|
||||||
if( verbose) printf(".......done.\n");
|
|
||||||
}
|
|
||||||
|
|
||||||
inline void FSUReader::ScanNumBlock(int verbose, uShort saveData){
|
|
||||||
if( inFile == nullptr ) return;
|
|
||||||
if( feof(inFile) ) return;
|
|
||||||
|
|
||||||
blockID = 0;
|
|
||||||
blockPos.push_back(0);
|
|
||||||
|
|
||||||
data->ClearData();
|
|
||||||
rewind(inFile);
|
|
||||||
filePos = 0;
|
|
||||||
|
|
||||||
bool isTraceOn = saveData < 2 ? false : true;
|
|
||||||
|
|
||||||
while( ReadNextBlock(isTraceOn, verbose - 1, saveData) == 0 ){
|
|
||||||
blockPos.push_back(filePos);
|
|
||||||
blockTimeStamp.push_back(data->aggTime);
|
|
||||||
blockID ++;
|
|
||||||
if(verbose && blockID % 10000 == 0) printf("%u, %.2f%% %u/%lu\n\033[A\r", blockID, filePos*100./inFileSize, filePos, inFileSize);
|
|
||||||
}
|
|
||||||
|
|
||||||
totNumBlock = blockID;
|
|
||||||
if(verbose) {
|
|
||||||
printf("\nScan complete: number of data Block : %lu\n", totNumBlock);
|
|
||||||
printf( " number of hit : %lu", hitCount);
|
|
||||||
if( hitCount > 1e6 ) printf(" = %.3f million", hitCount/1e6);
|
|
||||||
printf("\n");
|
|
||||||
if( saveData )printf( " size of the hit array : %lu\n", hit.size());
|
|
||||||
|
|
||||||
if( saveData ){
|
|
||||||
size_t sizeT = sizeof(hit[0]) * hit.size();
|
|
||||||
printf("size of hit array : %lu byte = %.2f kByte, = %.2f MByte\n", sizeT, sizeT/1024., sizeT/1024./1024.);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
if( fileList.size() > 0 ){
|
|
||||||
fileID = 0;
|
|
||||||
OpenFile(fileList[fileID], data->GetDataSize(), 0);
|
|
||||||
}
|
|
||||||
|
|
||||||
rewind(inFile);
|
|
||||||
blockID = 0;
|
|
||||||
filePos = 0;
|
|
||||||
|
|
||||||
//check is the hitCount == hit.size();
|
|
||||||
if( saveData ){
|
|
||||||
if( hitCount != hit.size()){
|
|
||||||
printf("!!!!!! the Data::dataSize is not big enough. !!!!!!!!!!!!!!!\n");
|
|
||||||
}else{
|
|
||||||
SortHit(verbose+1);
|
|
||||||
}
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
inline std::vector<Hit> FSUReader::ReadBatch(unsigned int batchSize, bool verbose){
|
|
||||||
|
|
||||||
// printf("%s sn:%d. filePos : %lu\n", __func__, sn, ftell(inFile));
|
|
||||||
|
|
||||||
std::vector<Hit> hitList_A;
|
|
||||||
if( IsEndOfFile() ) {
|
|
||||||
hitList_A = hit;
|
|
||||||
hit.clear();
|
|
||||||
return hitList_A;
|
|
||||||
}
|
|
||||||
|
|
||||||
if( hit.size() == 0 ){
|
|
||||||
int res = 0;
|
|
||||||
do{
|
|
||||||
res = ReadNextBlock(true, 0, 3);
|
|
||||||
}while ( hit.size() < batchSize && res == 0);
|
|
||||||
SortHit();
|
|
||||||
uLong t0_B = hit.at(0).timestamp;
|
|
||||||
uLong t1_B = hit.back().timestamp;
|
|
||||||
if( verbose ) {
|
|
||||||
printf(" hit in memeory : %7zu | %u | %lu \n", hit.size(), filePos, inFileSize);
|
|
||||||
printf("t0 : %15lu ns\n", t0_B);
|
|
||||||
printf("t1 : %15lu ns\n", t1_B);
|
|
||||||
printf("dt : %15.3f ms\n", (t1_B - t0_B)/1e6);
|
|
||||||
}
|
|
||||||
|
|
||||||
hitList_A = hit;
|
|
||||||
hit.clear();
|
|
||||||
|
|
||||||
}else{
|
|
||||||
|
|
||||||
hitList_A = hit;
|
|
||||||
hit.clear();
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
if( IsEndOfFile() ) return hitList_A; // when file finished for 1st batch read
|
|
||||||
|
|
||||||
int res = 0;
|
|
||||||
do{
|
|
||||||
res = ReadNextBlock(true, 0, 3);
|
|
||||||
}while ( hit.size() < batchSize && res == 0);
|
|
||||||
SortHit();
|
|
||||||
uLong t0_B = hit.at(0).timestamp;
|
|
||||||
uLong t1_B = hit.back().timestamp;
|
|
||||||
|
|
||||||
if( verbose ) {
|
|
||||||
printf(" hit in memeory : %7zu | %u | %lu \n", hit.size(), filePos, inFileSize);
|
|
||||||
printf("t0 : %15lu\n", t0_B);
|
|
||||||
printf("t1 : %15lu\n", t1_B);
|
|
||||||
printf("dt : %15.3f ms\n", (t1_B - t0_B)/1e6);
|
|
||||||
}
|
|
||||||
|
|
||||||
uLong t0_A = hitList_A.at(0).timestamp;
|
|
||||||
uLong t1_A = hitList_A.back().timestamp;
|
|
||||||
ulong ID_A = 0;
|
|
||||||
ulong ID_B = 0;
|
|
||||||
|
|
||||||
if( t0_A >= t0_B) {
|
|
||||||
printf("\033[0;31m!!!!!!!!!!!!!!!!! %s | Need to increase the batch size. \033[0m\n", __func__);
|
|
||||||
return std::vector<Hit> ();
|
|
||||||
}
|
|
||||||
|
|
||||||
if( t1_A > t0_B) { // need to sort between two hitList
|
|
||||||
|
|
||||||
if( verbose ) {
|
|
||||||
printf("############# need to sort \n");
|
|
||||||
printf("=========== sume of A + B : %zu \n", hitList_A.size() + hit.size());
|
|
||||||
}
|
|
||||||
|
|
||||||
std::vector<Hit> hitTemp;
|
|
||||||
|
|
||||||
// find the hit that is >= t0_B, save them to hitTemp
|
|
||||||
for( size_t j = 0; j < hitList_A.size() ; j++){
|
|
||||||
if( hitList_A[j].timestamp < t0_B ) continue;;
|
|
||||||
if( ID_A == 0 ) ID_A = j;
|
|
||||||
hitTemp.push_back(hitList_A[j]);
|
|
||||||
}
|
|
||||||
|
|
||||||
// remove hitList_A element that is >= t0_B
|
|
||||||
hitList_A.erase(hitList_A.begin() + ID_A, hitList_A.end() );
|
|
||||||
|
|
||||||
// find the hit that is <= t1_A, save them to hitTemp
|
|
||||||
for( size_t j = 0; j < hit.size(); j++){
|
|
||||||
if( hit[j].timestamp > t1_A ) {
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
hitTemp.push_back(hit[j]);
|
|
||||||
ID_B = j + 1;
|
|
||||||
}
|
|
||||||
|
|
||||||
// remove hit elements that is <= t1_A
|
|
||||||
hit.erase(hit.begin(), hit.begin() + ID_B );
|
|
||||||
|
|
||||||
// sort hitTemp
|
|
||||||
std::sort(hitTemp.begin(), hitTemp.end(), [](const Hit& a, const Hit& b) {
|
|
||||||
return a.timestamp < b.timestamp;
|
|
||||||
});
|
|
||||||
|
|
||||||
|
|
||||||
if( verbose ) {
|
|
||||||
printf("----------------- ID_A : %lu, Drop\n", ID_A);
|
|
||||||
printf("----------------- ID_B : %lu, Drop\n", ID_B);
|
|
||||||
PrintHitListInfo(&hitList_A, "hitList_A");
|
|
||||||
PrintHitListInfo(&hitTemp, "hitTemp");
|
|
||||||
PrintHitListInfo();
|
|
||||||
printf("=========== sume of A + B + Temp : %zu \n", hitList_A.size() + hit.size() + hitTemp.size());
|
|
||||||
printf("----------------- refill hitList_A \n");
|
|
||||||
}
|
|
||||||
|
|
||||||
for( size_t j = 0; j < hitTemp.size(); j++){
|
|
||||||
hitList_A.push_back(hitTemp[j]);
|
|
||||||
}
|
|
||||||
hitTemp.clear();
|
|
||||||
|
|
||||||
if( verbose ) {
|
|
||||||
PrintHitListInfo(&hitList_A, "hitList_A");
|
|
||||||
PrintHitListInfo();
|
|
||||||
printf("=========== sume of A + B : %zu \n", hitList_A.size() + hit.size());
|
|
||||||
}
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
return hitList_A;
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
/*
|
|
||||||
inline void FSUReader::SortAndSaveTS(unsigned int batchSize, bool verbose){
|
|
||||||
|
|
||||||
int count = 0;
|
|
||||||
std::vector<Hit> hitList_A ;
|
|
||||||
|
|
||||||
do{
|
|
||||||
|
|
||||||
if( verbose ) printf("***************************************************\n");
|
|
||||||
|
|
||||||
int res = 0;
|
|
||||||
do{
|
|
||||||
res = ReadNextBlock(true, 0, 3);
|
|
||||||
}while ( hit.size() < batchSize && res == 0);
|
|
||||||
|
|
||||||
SortHit();
|
|
||||||
uLong t0_B = hit.at(0).timestamp;
|
|
||||||
uLong t1_B = hit.back().timestamp;
|
|
||||||
|
|
||||||
if( verbose ) {
|
|
||||||
printf(" hit in memeory : %7zu | %u | %lu \n", hit.size(), filePos, inFileSize);
|
|
||||||
printf("t0 : %15lu\n", t0_B);
|
|
||||||
printf("t1 : %15lu\n", t1_B);
|
|
||||||
}
|
|
||||||
|
|
||||||
if( count == 0 ) {
|
|
||||||
hitList_A = hit; // copy hit
|
|
||||||
}else{
|
|
||||||
|
|
||||||
uLong t0_A = hitList_A.at(0).timestamp;
|
|
||||||
uLong t1_A = hitList_A.back().timestamp;
|
|
||||||
ulong ID_A = 0;
|
|
||||||
ulong ID_B = 0;
|
|
||||||
|
|
||||||
if( t0_A > t0_B) {
|
|
||||||
printf("Need to increase the batch size. \n");
|
|
||||||
return;
|
|
||||||
}
|
|
||||||
|
|
||||||
if( t1_A > t0_B) { // need to sort between two hitList
|
|
||||||
|
|
||||||
if( verbose ) {
|
|
||||||
printf("############# need to sort \n");
|
|
||||||
printf("=========== sume of A + B : %zu \n", hitList_A.size() + hit.size());
|
|
||||||
}
|
|
||||||
|
|
||||||
std::vector<Hit> hitTemp;
|
|
||||||
|
|
||||||
for( size_t j = 0; j < hitList_A.size() ; j++){
|
|
||||||
if( hitList_A[j].timestamp < t0_B ) continue;
|
|
||||||
if( ID_A == 0 ) ID_A = j;
|
|
||||||
hitTemp.push_back(hitList_A[j]);
|
|
||||||
}
|
|
||||||
|
|
||||||
hitList_A.erase(hitList_A.begin() + ID_A, hitList_A.end() );
|
|
||||||
if( verbose ) {
|
|
||||||
printf("----------------- ID_A : %lu, Drop\n", ID_A);
|
|
||||||
PrintHitListInfo(hitList_A, "hitList_A");
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
for( size_t j = 0; j < hit.size(); j++){
|
|
||||||
if( hit[j].timestamp > t1_A ) {
|
|
||||||
ID_B = j;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
hitTemp.push_back(hit[j]);
|
|
||||||
}
|
|
||||||
|
|
||||||
std::sort(hitTemp.begin(), hitTemp.end(), [](const Hit& a, const Hit& b) {
|
|
||||||
return a.timestamp < b.timestamp;
|
|
||||||
});
|
|
||||||
|
|
||||||
hit.erase(hit.begin(), hit.begin() + ID_B );
|
|
||||||
|
|
||||||
if( verbose ) {
|
|
||||||
PrintHitListInfo(hitTemp, "hitTemp");
|
|
||||||
printf("----------------- ID_B : %lu, Drop\n", ID_B);
|
|
||||||
PrintHitListInfo(hit, "hit");
|
|
||||||
printf("=========== sume of A + B + Temp : %zu \n", hitList_A.size() + hit.size() + hitTemp.size());
|
|
||||||
printf("----------------- refill hitList_A \n");
|
|
||||||
}
|
|
||||||
ulong ID_Temp = 0;
|
|
||||||
for( size_t j = 0; j < hitTemp.size(); j++){
|
|
||||||
hitList_A.push_back(hitTemp[j]);
|
|
||||||
if( hitList_A.size() >= batchSize ) {
|
|
||||||
ID_Temp = j+1;
|
|
||||||
break;
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
hitTemp.erase(hitTemp.begin(), hitTemp.begin() + ID_Temp );
|
|
||||||
for( size_t j = 0 ; j < hit.size(); j ++){
|
|
||||||
hitTemp.push_back(hit[j]);
|
|
||||||
}
|
|
||||||
SaveHit(hitList_A, count <= 1 ? false : true);
|
|
||||||
|
|
||||||
if( verbose ) {
|
|
||||||
PrintHitListInfo(hitList_A, "hitList_A");
|
|
||||||
PrintHitListInfo(hitTemp, "hitTemp");
|
|
||||||
printf("----------------- replace hitList_A by hitTemp \n");
|
|
||||||
}
|
|
||||||
|
|
||||||
hitList_A.clear();
|
|
||||||
hitList_A = hitTemp;
|
|
||||||
hit.clear();
|
|
||||||
|
|
||||||
if( verbose ) {
|
|
||||||
PrintHitListInfo(hitList_A, "hitList_A");
|
|
||||||
printf("===========================================\n");
|
|
||||||
}
|
|
||||||
|
|
||||||
}else{ // save hitList_A, replace hitList_A
|
|
||||||
|
|
||||||
SaveHit(hitList_A, count <= 1? false : true);
|
|
||||||
hitList_A.clear();
|
|
||||||
hitList_A = hit;
|
|
||||||
if( verbose ) PrintHitListInfo(hitList_A, "hitList_A");
|
|
||||||
|
|
||||||
}
|
|
||||||
}
|
|
||||||
|
|
||||||
ClearHitList();
|
|
||||||
count ++;
|
|
||||||
}while(filePos < inFileSize);
|
|
||||||
|
|
||||||
SaveHit(hitList_A, count <= 1 ? false : true);
|
|
||||||
|
|
||||||
printf("================= finished.\n");
|
|
||||||
}
|
|
||||||
*/
|
|
||||||
|
|
||||||
/*
|
|
||||||
inline std::string FSUReader::SaveHit(std::vector<Hit> hitList, bool isAppend){
|
|
||||||
|
|
||||||
std::string outFileName;
|
|
||||||
if( fileList.empty() ) {
|
|
||||||
outFileName = fileName + ".ts" ;
|
|
||||||
}else{
|
|
||||||
outFileName = fileList[0] + ".ts" ;
|
|
||||||
}
|
|
||||||
uint64_t hitSize = hitList.size();
|
|
||||||
|
|
||||||
FILE * outFile ;
|
|
||||||
if( isAppend ) {
|
|
||||||
outFile = fopen(outFileName.c_str(), "rb+"); //read/write bineary
|
|
||||||
|
|
||||||
rewind(outFile);
|
|
||||||
fseek( outFile, 4, SEEK_CUR);
|
|
||||||
uint64_t org_hitSize;
|
|
||||||
fread(&org_hitSize, 8, 1, outFile);
|
|
||||||
|
|
||||||
rewind(outFile);
|
|
||||||
fseek( outFile, 4, SEEK_CUR);
|
|
||||||
|
|
||||||
org_hitSize += hitSize;
|
|
||||||
|
|
||||||
fwrite(&org_hitSize, 8, 1, outFile);
|
|
||||||
fseek(outFile, 0, SEEK_END);
|
|
||||||
|
|
||||||
}else{
|
|
||||||
outFile = fopen(outFileName.c_str(), "wb"); //overwrite binary
|
|
||||||
uint32_t header = 0xAA000000;
|
|
||||||
header += sn;
|
|
||||||
fwrite( &header, 4, 1, outFile );
|
|
||||||
fwrite( &hitSize, 8, 1, outFile);
|
|
||||||
}
|
|
||||||
|
|
||||||
|
|
||||||
for( ulong i = 0; i < hitSize; i++){
|
|
||||||
|
|
||||||
if( i% 10000 == 0 ) printf("Saving %lu/%lu Hit (%.2f%%)\n\033[A\r", i, hitSize, i*100./hitSize);
|
|
||||||
|
|
||||||
uint16_t flag = hitList[i].ch + (hitList[i].pileUp << 8) ;
|
|
||||||
|
|
||||||
if( DPPType == DPPTypeCode::DPP_PSD_CODE ) flag += ( 1 << 15);
|
|
||||||
if( hitList[i].traceLength > 0 ) flag += (1 << 14);
|
|
||||||
|
|
||||||
// fwrite( &(hit[i].ch), 1, 1, outFile);
|
|
||||||
fwrite( &flag, 2, 1, outFile);
|
|
||||||
fwrite( &(hitList[i].energy), 2, 1, outFile);
|
|
||||||
if( DPPType == DPPTypeCode::DPP_PSD_CODE ) fwrite( &(hitList[i].energy2), 2, 1, outFile);
|
|
||||||
fwrite( &(hitList[i].timestamp), 6, 1, outFile);
|
|
||||||
fwrite( &(hitList[i].fineTime), 2, 1, outFile);
|
|
||||||
if( hitList[i].traceLength > 0 ) fwrite( &(hitList[i].traceLength), 2, 1, outFile);
|
|
||||||
|
|
||||||
for( uShort j = 0; j < hitList[i].traceLength; j++){
|
|
||||||
fwrite( &(hitList[i].trace[j]), 2, 1, outFile);
|
|
||||||
}
|
|
||||||
|
|
||||||
}
|
|
||||||
|
|
||||||
off_t tsFileSize = ftello(outFile); // unsigned int = Max ~4GB
|
|
||||||
fclose(outFile);
|
|
||||||
|
|
||||||
printf("Saved to %s, size: ", outFileName.c_str());
|
|
||||||
if( tsFileSize < 1024 ) {
|
|
||||||
printf(" %ld Byte", tsFileSize);
|
|
||||||
}else if( tsFileSize < 1024*1024 ) {
|
|
||||||
printf(" %.2f kB", tsFileSize/1024.);
|
|
||||||
}else if( tsFileSize < 1024*1024*1024){
|
|
||||||
printf(" %.2f MB", tsFileSize/1024./1024.);
|
|
||||||
}else{
|
|
||||||
printf(" %.2f GB", tsFileSize/1024./1024./1024.);
|
|
||||||
}
|
|
||||||
printf("\n");
|
|
||||||
|
|
||||||
return outFileName;
|
|
||||||
|
|
||||||
}
|
|
||||||
*/
|
|
||||||
|
|
||||||
|
|
||||||
3594
Armory/mass20.txt
3594
Armory/mass20.txt
File diff suppressed because it is too large
Load Diff
|
|
@ -1,7 +0,0 @@
|
||||||
.L ANASEN_model.C
|
|
||||||
.L anasenMS_root.cpp+
|
|
||||||
|
|
||||||
void test_run(){
|
|
||||||
ANASEN_model();
|
|
||||||
Run(10);
|
|
||||||
}
|
|
||||||
Loading…
Reference in New Issue
Block a user