energy dep comments

This commit is contained in:
James Szalkie 2026-03-30 15:21:44 -04:00
parent e78b378e05
commit 6e969434da
2 changed files with 35 additions and 23 deletions

View File

@ -26,18 +26,30 @@ void ANASEN_model(int anodeID1 = -1, int anodeID2 = -1, int cathodeID1 = -1, int
TGeoVolume *worldBox = geom->MakeBox("ROOT", Vacuum, worldx, worldy, worldz);
geom->SetTopVolume(worldBox);
//--- making axis
TGeoVolume *axisX = geom->MakeTube("axisX", Al, 0, 0.1, 5.);
axisX->SetLineColor(1);
worldBox->AddNode(axisX, 1, new TGeoCombiTrans(5, 0, 0., new TGeoRotation("rotA", 90., 90., 0.)));
//--- making axis (X=red, Y=green, Z=blue)
Double_t axisLen = 100.;
TGeoVolume *axisX = geom->MakeTube("axisX", Al, 0, 0.1, axisLen/2);
axisX->SetLineColor(kRed);
worldBox->AddNode(axisX, 1, new TGeoCombiTrans(axisLen/2, 0, 0., new TGeoRotation("rotA", 90., 90., 0.)));
TGeoVolume *axisY = geom->MakeTube("axisY", Al, 0, 0.1, 5.);
axisY->SetLineColor(1);
worldBox->AddNode(axisY, 1, new TGeoCombiTrans(0, 5, 0., new TGeoRotation("rotB", 0., 90., 0.)));
TGeoVolume *axisY = geom->MakeTube("axisY", Al, 0, 0.1, axisLen/2);
axisY->SetLineColor(kGreen);
worldBox->AddNode(axisY, 1, new TGeoCombiTrans(0, axisLen/2, 0., new TGeoRotation("rotB", 0., 90., 0.)));
TGeoVolume *axisZ = geom->MakeTube("axisZ", Al, 0, 0.1, 5.);
axisZ->SetLineColor(1);
worldBox->AddNode(axisZ, 1, new TGeoTranslation(0, 0, 5));
TGeoVolume *axisZ = geom->MakeTube("axisZ", Al, 0, 0.1, axisLen/2);
axisZ->SetLineColor(kBlue);
worldBox->AddNode(axisZ, 1, new TGeoTranslation(0, 0, axisLen/2));
//--- axis labels (draw as TPolyMarker3D + text because TGeo does not label directly)
TPolyMarker3D *marker = new TPolyMarker3D();
marker->SetMarkerSize(1.2);
marker->SetMarkerColor(kRed);
marker->SetPoint(0, axisLen, 0, 0); // X
marker->SetMarkerColor(kGreen);
marker->SetPoint(1, 0, axisLen, 0); // Y
marker->SetMarkerColor(kBlue);
marker->SetPoint(2, 0, 0, axisLen); // Z
marker->Draw();
//--- making ANASEN
const int nWire = 24;

View File

@ -237,7 +237,7 @@ int main(int argc, char **argv){
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);
@ -247,19 +247,19 @@ int main(int argc, char **argv){
sx3Z = hitPos.Z();
// apply energy loss from vertex to SX3 hit position (for light particle)
double dl = (hitPos - vertex).Mag() / 10.0; // path length in cm (positions in mm)
double EkinLight = Pb.E() - Pb.M();
double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx
double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling)
if (dE_light < EkinLight) {
Pb.SetE(Pb.E() - dE_light);
// update momentum to conserve direction
double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M());
TVector3 dir_new = Pb.Vect().Unit() * p_new;
Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E());
}
//double dl = (hitPos - vertex).Mag() / 10.0; // path length in cm (positions in mm)
//double EkinLight = Pb.E() - Pb.M();
//double dedxLight = elossLight->Eval(EkinLight); // interpolate dE/dx
//double dE_light = dedxLight * dl * density / 1000.0; // adjust for units (example scaling)
//if (dE_light < EkinLight) {
// Pb.SetE(Pb.E() - dE_light);
// // update momentum to conserve direction
// double p_new = TMath::Sqrt(Pb.E()*Pb.E() - Pb.M()*Pb.M());
// TVector3 dir_new = Pb.Vect().Unit() * p_new;
// Pb.SetPxPyPzE(dir_new.X(), dir_new.Y(), dir_new.Z(), Pb.E());
//}
// update kinetic energy after loss
Tb = Pb.E() - Pb.M();
//Tb = Pb.E() - Pb.M();
// reconstruct track from PW readings + SX3 hit
pw->CalTrack(hitPos, anodeID[0], cathodeID[0], false);