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

catches all events

Browse Source
This commit is contained in:
James Szalkie 2026-08-25 12:01:40 -04:00
parent 8240310039
commit 0a9dc19904
1 changed files with 48 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

79
Armory/anasenMS.cpp
View File

@ -264,6 +264,7 @@ int main(int argc, char **argv){
TTree * tree1 = new TTree("tree1", "tree1");
// beam and CM variables saved in tree
double eventID;
double KEA;
double beamPath_cm;
int MBeamOut;
@ -275,6 +276,7 @@ int main(int argc, char **argv){
int ZLightOut;
int ZHeavyOut;
double beamEnergy;
tree1->Branch("eventID", &eventID, "eventID/D");
tree1->Branch("beamKEA", &KEA, "beamKEA/D");
tree1->Branch("beamPath_cm", &beamPath_cm, "beamPath_cm/D");
tree1->Branch("beamEnergy", &beamEnergy, "beamEnergy/D");
@ -414,6 +416,42 @@ int main(int argc, char **argv){
for( int i = 0; i < numEvent ; i++){
if(quit) break; // exit gracefully if signal Ctrl+C received
// randomly sample target/projectile excitations
eventID = i;
thetaCM = TMath::QuietNaN();
phiCM = TMath::QuietNaN();
thetab = TMath::QuietNaN();
phib = TMath::QuietNaN();
Tb = TMath::QuietNaN();
thetaB = TMath::QuietNaN();
phiB = TMath::QuietNaN();
TB = TMath::QuietNaN();
T = {TMath::QuietNaN(), TMath::QuietNaN()};
qqqTb = TMath::QuietNaN();
qqqTB = TMath::QuietNaN();
sx3Tb = TMath::QuietNaN();
sx3TB = TMath::QuietNaN();
Esx3 = TMath::QuietNaN();
Eqqq = TMath::QuietNaN();
Eanode = TMath::QuietNaN();
Ecathode = TMath::QuietNaN();
EPC = TMath::QuietNaN();
sx3X = TMath::QuietNaN();
sx3Y = TMath::QuietNaN();
sx3Z = TMath::QuietNaN();
qqqX = TMath::QuietNaN();
qqqY = TMath::QuietNaN();
qqqZ = TMath::QuietNaN();
detX = TMath::QuietNaN();
detY = TMath::QuietNaN();
detZ = TMath::QuietNaN();
sx3ID = -1;
sx3Up = -1;
sx3Dn = -1;
sx3Bk = -1;
qqqID = -1;
ExAID = gRandom->Integer(nExA);
ExA = ExAList[ExAID];
transfer.SetExA(ExA);
@ -491,24 +529,25 @@ int main(int argc, char **argv){
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;
//if(IsDeadAnode(anodeID[0])) continue;
//if(IsDeadCathode(cathodeID[0])) continue;
// SX3 hit channel info and depth fraction
sx3ID = sx3->GetID();
qqqID = qqq->GetID();
if(IsDeadSX3(sx3ID)) continue;
//if(IsDeadSX3(sx3ID)) continue;
anodeDist[0] = hitInfo.nearestDist.first; // distance to nearest anode wire
cathodeDist[0] = hitInfo.nearestDist.second; // distance to nearest cathode wire
//start HERE
if( sx3ID >= 0 ){
sx3Up = sx3->GetChUp();
sx3Dn = sx3->GetChDn();
sx3Bk = sx3->GetChBk();
if(IsDeadSX3ChannelCombo(sx3ID, sx3Up, sx3Dn, sx3Bk)) continue;
//if(IsDeadSX3ChannelCombo(sx3ID, sx3Up, sx3Dn, sx3Bk)) continue;
sx3ZFrac = sx3->GetZFrac();
// apply intrinsic detector resolution to true SX3 hit position
@ -586,41 +625,21 @@ int main(int argc, char **argv){
b, "He", Tb,
distance_C);
if (Esx3 <= 0 || Eanode <= 0 || Ecathode <= 0) {
Esx3 = NAN;
beamEnergy = NAN;
Ex = NAN;
continue;
//continue;
}
sx3Tb = Tb; // for simplicity, using the same kinetic energy for SX3 hit events, can be modified to simulate energy loss if desired
sx3TB = TB;
qqqTb = TMath::QuietNaN(); // mark kinetic energy as invalid for SX3 hit case
qqqTB = TMath::QuietNaN();
Eqqq = TMath::QuietNaN(); // mark QQQ energy as invalid for SX3 hit case
Edet = Esx3;
EPC = Eanode - Ecathode;
AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)");
AutoHist2D::Fill("EPC x sin(theta) vs Esx3", Esx3, EPC * sin(thetab * TMath::DegToRad()), "Esx3 (MeV)", "EPC x sin(theta) (MeV)");
tree1->Fill();
//tree1->Fill();
}else if (qqqID >= 0){
// handle QQQ hit case
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();
TVector3 hitPos = qqq->GetHitPos();
@ -666,21 +685,18 @@ int main(int argc, char **argv){
if (Eqqq <= 0 || Eanode <= 0 || Ecathode <= 0) {
Eqqq = NAN;
Ex = NAN;
continue;
//continue;
}
qqqTb = Tb; // for simplicity, using the same kinetic energy for QQQ hit events, can be modified to simulate energy loss if desired
qqqTB = TB;
sx3Tb = TMath::QuietNaN(); // mark kinetic energy as invalid for QQQ hit case
sx3TB = TMath::QuietNaN();
Esx3 = TMath::QuietNaN(); // mark SX3 energy as invalid for QQQ hit case
Edet = Eqqq;
EPC = Eanode - Ecathode;
AutoHist2D::Fill("beamEnergy_vs_vZ", vertexZ / 10, beamEnergy, "vZ (cm)", "beamEnergy (MeV)");
AutoHist2D::Fill("EPC x sin(theta) vs Eqqq", Eqqq, EPC * sin(thetab * TMath::DegToRad()), "Eqqq (MeV)", "EPC x sin(theta) (MeV)");
beamEnergy = TMath::QuietNaN(); // mark beam energy as invalid for QQQ hit case
tree1->Fill();
//tree1->Fill();
}else{
// no valid SX3 hit: mark clearly invalid
@ -708,6 +724,7 @@ int main(int argc, char **argv){
//comment out tree fill for no hit case
//tree1->Fill();
}
tree1->Fill();
//#################################################################### Timer
// measure elapsed real time and print progress roughly every 10 sec