diff --git a/Armory/EnergyLoss.h b/Armory/EnergyLoss.h index ab9ea7f..1b441bb 100644 --- a/Armory/EnergyLoss.h +++ b/Armory/EnergyLoss.h @@ -99,4 +99,17 @@ inline double CalculateEnergyLoss(double vx, double vy, double vz, return table.EnergyAfterDistance(initialEnergy_MeV, distance_cm); } +inline double CalculateOriginalEnergy(double fx, double fy, double fz, + double vx, double vy, double vz, + const std::string& particle, + const std::string& medium, + double initialEnergy_MeV, + double& distance_cm) { + const EnergyLossTable& table = EnergyLossTable::Get(particle, medium); + distance_cm = - CalcPathLength_cm(vx, vy, vz, fx, fy, fz); + return table.EnergyAfterDistance(initialEnergy_MeV, distance_cm); +} + + + #endif diff --git a/Armory/Kinematics.h b/Armory/Kinematics.h index 768dcdc..5d47d28 100644 --- a/Armory/Kinematics.h +++ b/Armory/Kinematics.h @@ -5,6 +5,7 @@ #include #include #include +#include "EnergyLoss.h" const double u_MeV = 931.49410372; //u in MeV class Kinematics { diff --git a/Armory/anasenMS.cpp b/Armory/anasenMS.cpp index e8cb372..9d92c58 100644 --- a/Armory/anasenMS.cpp +++ b/Armory/anasenMS.cpp @@ -63,7 +63,7 @@ bool IsDeadCathode(int id){ bool IsDeadSX3(int id){ static std::set dead = { - //0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23 + 0, 2, 4, 5, 6, 8, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23 }; // add dead SX3 IDs here, 0-23 1,7,9,3 return dead.count(id); } @@ -84,7 +84,7 @@ bool IsDeadSX3FrontDnChannel(int sx3ID, int chDn){ bool IsDeadSX3BackChannel(int sx3ID, int chBk){ static std::set> dead = { - //{1, 10} + {1, 10} // {sx3ID, back-channel} }; return dead.count({sx3ID, chBk}); @@ -161,7 +161,7 @@ int main(int argc, char **argv){ // Excited state lists (projectile and heavy-product excitation states) std::vector ExAList = {0}; // Beam excited energy - std::vector ExList = {0.0, 2.2, 3.4, 6.0}; // Heavy product excited energy + std::vector ExList = {0, 2.2, 3.4, 6.0}; // Heavy product excited energy const int kMBeam = reactionConfig.beamA; // mass number of beam const int kMTarget = reactionConfig.targetA; // mass number of target @@ -229,8 +229,8 @@ int main(int argc, char **argv){ elossBeamInverse->Sort(); // TGraph::Eval requires ascending x // detector resolution / uncertainty parameters - double sigmaSX3_W = 0; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension) - double sigmaSX3_L = 0; // mm, vertical spread for SX3 + double sigmaSX3_W = -1; // mm, if < 0 use mid-point (no spread in SX3 horizontal dimension) + double sigmaSX3_L = -1; // 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) @@ -407,8 +407,7 @@ int main(int argc, char **argv){ tree1->Branch("qqqUp", &qqqUp, "qqqUp/I"); tree1->Branch("qqqBk", &qqqBk, "qqqBk/I"); - double EBeam_Kin_gs=NAN, EBeam_Kin_2_2=NAN, EBeam_Kin_3_4=NAN, Ex_recon=NAN - ; + double EBeam_Kin_gs=NAN, EBeam_Kin_2_2=NAN, EBeam_Kin_3_4=NAN, Ex_recon=NAN; tree1->Branch("EBeam_Kin", &EBeam_Kin_gs, "EBeam_Kin/D"); tree1->Branch("EBeam_Kin_2.2", &EBeam_Kin_2_2, "EBeam_Kin_2.2/D"); tree1->Branch("EBeam_Kin_3.4", &EBeam_Kin_3_4, "EBeam_Kin_3.4/D"); @@ -423,6 +422,10 @@ int main(int argc, char **argv){ tree1->Branch("reTheta1", &reTheta1, "reconstucted_theta1/D"); tree1->Branch("rePhi1", &rePhi1, "reconstucted_phi1/D"); + double hitTheta, originalEnergy; + tree1->Branch("hitTheta", &hitTheta, "hitTheta/D"); + tree1->Branch("originalEnergy", &originalEnergy, "originalEnergy/D"); + // reconstructed vertex Z from PW fit double z0; tree1->Branch("z0", &z0, "reconstucted_Z/D"); @@ -490,7 +493,7 @@ int main(int argc, char **argv){ KEA = beamEnergy / beamA; beamPath_cm = elossBeamInverse->Eval(beamEnergy); // beamE maps to x=0 after path shift vertexZ = beamEntranceZ + beamPath_cm * 10.0; // cm -> mm - + // transverse sampling range from the beam's distance straggle at this energy const double sigmaX_mm = std::max(0.0, sigmaXBeam->Eval(beamEnergy)) * 100.0; // cm -> mm //vertexX = 2.0 * sigmaX_mm * gRandom->Rndm() - sigmaX_mm; @@ -501,21 +504,19 @@ int main(int argc, char **argv){ vertexX = gRandom->Gaus(0, vertexRangeX); // mean and standard deviation vertexY = gRandom->Gaus(0, vertexRangeY); // mean and standard deviation - //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); // compute beam energy at the event vertex from the gas path length - //beamPath_cm = TVector3(vertexZ - beamEntranceZ, vertexX, vertexY).Mag() * 0.1; - //beamDistance = vertexZ - beamEntranceZ; - /* + beamPath_cm = TVector3(vertexZ - beamEntranceZ, vertexX, vertexY).Mag() * 0.1; + beamDistance = vertexZ - beamEntranceZ; + if( beamPath_cm < 0 ) beamPath_cm = 0; - beamEnergy = elossBeam->Eval(beamPath_cm); // MeV - double beamEnergyLoss = elossBeam->Eval(0.0) - beamEnergy; - KEA = beamEnergy / beamA;*/ + //beamEnergy = elossBeam->Eval(beamPath_cm); // MeV + //double beamEnergyLoss = elossBeam->Eval(0.0) - beamEnergy; + //KEA = beamEnergy / beamA; //KEA = gRandom->Uniform(0, beamE); transfer.SetIncidentEnergyAngle(KEA, 0, 0); @@ -586,7 +587,7 @@ int main(int argc, char **argv){ // 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); + TVector3 hitPos = sx3->GetHitPos(); sx3X = hitPos.X(); sx3Y = hitPos.Y(); @@ -647,6 +648,18 @@ int main(int argc, char **argv){ distance_sx3); dl = distance_sx3; + //std::cout << std::lround(sx3X / 10) * 10 << " " << std::lround(sx3Y / 10) * 10 << " " << std::lround(sx3Z / 10) * 10 << std::endl; + TVector3 hitSigma = sx3->GetHitPosWithSigma(sigmaSX3_W, sigmaSX3_L); + + double hitX = hitSigma.X(); + double hitY = hitSigma.Y(); + double hitZ = hitSigma.Z(); + + originalEnergy = CalculateOriginalEnergy(hitX, hitY, hitZ, + 0, 0, std::lround(vertexZ / 10) * 10, + b, "He", Esx3, + distance_sx3); + double distance_A; Eanode = CalculateEnergyLoss(vertexX, vertexY, vertexZ, aX, aY, aZ, @@ -673,14 +686,13 @@ int main(int argc, char **argv){ AutoHist2D::Fill("EPC x sin(theta) vs Esx3", Esx3, EPC * sin(thetab * TMath::DegToRad()), "Esx3 (MeV)", "EPC x sin(theta) (MeV)"); //tree1->Fill(); - //Kinematics aakin_27Al(26.981538408,4.00260325413,4.0026035413,26.981538408,beam_energy_at_vertex/26.981538408); //m3 is alpha Kinematics apkin_27Al(26.981538408,4.00260325413,1.00782503224,29.973770136,beamEnergy/26.981538408); //m3 is proton - //Kinematics apkin_27Al(1.00782503224,4.00260325413,4.00260325413,1.00782503224,beamEnergy/1.00782503224); //m3 is proton - Ex_recon = apkin_27Al.getExc(Esx3, thetab); - EBeam_Kin_gs = apkin_27Al.getEbeam_givenQ(Esx3, 0.0, thetab); - EBeam_Kin_2_2 = apkin_27Al.getEbeam_givenQ(Esx3, 2.2, thetab); - EBeam_Kin_3_4 = apkin_27Al.getEbeam_givenQ(Esx3, 3.4, thetab); - //std::cout << EBeam_Kin << std::endl; + Ex_recon = apkin_27Al.getExc(originalEnergy, std::lround(thetab / 10) * 10); + //EBeam_Kin_gs = apkin_27Al.getEbeam_givenQ(Esx3, 0.0, thetab); + //EBeam_Kin_2_2 = apkin_27Al.getEbeam_givenQ(Esx3, 2.2, thetab); + //EBeam_Kin_3_4 = apkin_27Al.getEbeam_givenQ(Esx3, 3.4, thetab); + //std::cout << std::lround(thetab / 10) * 10 << std::endl; + hitTheta = hitPos.Theta() * 180.0 / TMath::Pi(); }else if (qqqID >= 0){ diff --git a/Armory/efficiency study 8-27.zip b/Armory/efficiency study 8-27.zip new file mode 100644 index 0000000..5951916 Binary files /dev/null and b/Armory/efficiency study 8-27.zip differ diff --git a/Armory/testMacro.C b/Armory/testMacro.C index febf3d7..736db1d 100644 --- a/Armory/testMacro.C +++ b/Armory/testMacro.C @@ -2,10 +2,35 @@ TFile *file0 = TFile::Open("SimAnasen1.root"); file0->ls(); TTree *tree1 = (TTree*)(file0->Get("tree1")); + new TBrowser(); tree1->Draw("Tb:thetab","","col"); tree1->Draw("Tb:thetab","vZ>-140 && vZ<-130", "box same"); - TH2F *h2 = new TH2F("EPCvEsx3", "EPC x sin(thetab) vs Esx3;Esx3;EPC * sin(thetab * TMath::DegToRad())", 200, 0, 0, 200, 0, 0); + TH2F *h2 = new TH2F("EPCvEsx3", "EPC x sin(thetab) vs Esx3;Esx3;EPC * sin(thetab * TMath::DegToRad())", 200, 0, 0, 200, 0, 0); //arguments are (name, title, nbinsX, xlow, xup, nbinsY, ylow, yup) tree1->Draw("EPC*sin(thetab * TMath::DegToRad()):Esx3 >> EPCvEsx3", "Esx3 > 0", "colz"); - new TBrowser(); + + new TCanvas(); + TH2F *h3 = new TH2F( + "hExReconVsVZ", + "Ex_recon vs vZ;vZ (mm);Excitation energy (MeV)", + 100, 0, 0, + 100, -1, 0 + ); + tree1->Draw("Ex_recon:vZ>>hExReconVsVZ", "sx3ID >= 0", "colz"); + + TH2F *h4 = new TH2F( + "hExVsVZ", + "", + 100, 0, 0, + 100, -1, 0 + ); + tree1->Draw("Ex:vZ>>hExVsVZ", "sx3ID >= 0", "box same"); + + new TCanvas(); + tree1->Draw("originalEnergy:Tb", "sx3ID >=0", ""); + TF1 *f_diagonal = new TF1("f_diagonal", "x", 0, 200); // assuming the range for Tb is 0 to 200 + f_diagonal->SetLineColor(kRed); + f_diagonal->SetLineStyle(2); + f_diagonal->Draw("same"); + }