From fb5b2787e2861cc666e0d54c3fee4f738571db91 Mon Sep 17 00:00:00 2001 From: Barak Schmookler Date: Fri, 26 Jan 2024 17:41:46 -0500 Subject: [PATCH] Update --- .../{HEXPLIT_EICRecon => }/DEMP_reco.C | 92 ++--- .../analysis/HEXPLIT_standalone/DEMP_reco.C | 325 ------------------ .../HEXPLIT_standalone/ZDC_neutron_recon.h | 65 ---- .../analysis/HEXPLIT_standalone/demp_tests.C | 249 -------------- .../analysis/HEXPLIT_standalone/hexplit.h | 190 ---------- 5 files changed, 52 insertions(+), 869 deletions(-) rename benchmarks/demp/analysis/{HEXPLIT_EICRecon => }/DEMP_reco.C (83%) delete mode 100644 benchmarks/demp/analysis/HEXPLIT_standalone/DEMP_reco.C delete mode 100644 benchmarks/demp/analysis/HEXPLIT_standalone/ZDC_neutron_recon.h delete mode 100644 benchmarks/demp/analysis/HEXPLIT_standalone/demp_tests.C delete mode 100644 benchmarks/demp/analysis/HEXPLIT_standalone/hexplit.h diff --git a/benchmarks/demp/analysis/HEXPLIT_EICRecon/DEMP_reco.C b/benchmarks/demp/analysis/DEMP_reco.C similarity index 83% rename from benchmarks/demp/analysis/HEXPLIT_EICRecon/DEMP_reco.C rename to benchmarks/demp/analysis/DEMP_reco.C index fd6414e9..5149ed1f 100644 --- a/benchmarks/demp/analysis/HEXPLIT_EICRecon/DEMP_reco.C +++ b/benchmarks/demp/analysis/DEMP_reco.C @@ -122,10 +122,10 @@ void DEMP_reco(){ TTreeReaderArray rec_mass(tr, "ReconstructedChargedParticles.mass"); //HEXPLIT Clusters - TTreeReaderArray zdc_cluster_energy(tr, "ZDC_HEXPLITClusters.energy"); - TTreeReaderArray zdc_cluster_x(tr, "ZDC_HEXPLITClusters.position.x"); - TTreeReaderArray zdc_cluster_y(tr, "ZDC_HEXPLITClusters.position.y"); - TTreeReaderArray zdc_cluster_z(tr, "ZDC_HEXPLITClusters.position.z"); + TTreeReaderArray zdc_cluster_energy(tr, "HcalFarForwardZDCClusters.energy"); + TTreeReaderArray zdc_cluster_x(tr, "HcalFarForwardZDCClusters.position.x"); + TTreeReaderArray zdc_cluster_y(tr, "HcalFarForwardZDCClusters.position.y"); + TTreeReaderArray zdc_cluster_z(tr, "HcalFarForwardZDCClusters.position.z"); TTreeReaderArray>> weight_map(tr,"_stringMap"); @@ -204,62 +204,74 @@ void DEMP_reco(){ }//End loop over reconstructed particles //Neutron Reconstruction using HEXSPLIT algorithm - auto ZDC_theta = std::asin( std::hypot(zdc_cluster_x[0],zdc_cluster_y[0])/zdc_cluster_z[0] ); - auto ZDC_phi = std::atan2(zdc_cluster_y[0],zdc_cluster_x[0]); - - auto neut_rec_E = zdc_cluster_energy[0]; + + double ZDC_theta(0),ZDC_phi(0),neut_rec_E(0); + double neut_rec_p(0),neut_rec_px(0),neut_rec_py(0),neut_rec_pz(0); double neut_mass = 0.9396; - auto neut_rec_p = std::sqrt( neut_rec_E*neut_rec_E - neut_mass*neut_mass ); - auto neut_rec_px = neut_rec_p * sin(ZDC_theta) * cos(ZDC_phi); - auto neut_rec_py = neut_rec_p * sin(ZDC_theta) * sin(ZDC_phi); - auto neut_rec_pz = neut_rec_p * cos(ZDC_theta); + TLorentzVector neut_rec; + + if(zdc_cluster_x.GetSize()>0){ + //Just take first cluster for now + ZDC_theta = std::atan( std::hypot(zdc_cluster_x[0],zdc_cluster_y[0])/zdc_cluster_z[0] ); + ZDC_phi = std::atan2(zdc_cluster_y[0],zdc_cluster_x[0]); + neut_rec_E = zdc_cluster_energy[0]; + neut_rec_p = std::sqrt( neut_rec_E*neut_rec_E - neut_mass*neut_mass ); + neut_rec_px = neut_rec_p * sin(ZDC_theta) * cos(ZDC_phi); + neut_rec_py = neut_rec_p * sin(ZDC_theta) * sin(ZDC_phi); + neut_rec_pz = neut_rec_p * cos(ZDC_theta); + - TLorentzVector neut_rec(neut_rec_px,neut_rec_py,neut_rec_pz,neut_rec_E); + neut_rec.SetPxPyPzE(neut_rec_px,neut_rec_py,neut_rec_pz,neut_rec_E); - h2_neut->Fill(neut_rec.Theta()*TMath::RadToDeg(),neut_rec.E()); + h2_neut->Fill(neut_rec.Theta()*TMath::RadToDeg(),neut_rec.E()); - //W.r.t proton beam direction - TLorentzVector rotated = neut_rec; - rotated.RotateY(0.025); - h4_neut->Fill(rotated.Theta()*TMath::RadToDeg(),rotated.Phi()*TMath::RadToDeg()); + //W.r.t proton beam direction + TLorentzVector rotated = neut_rec; + rotated.RotateY(0.025); + h4_neut->Fill(rotated.Theta()*TMath::RadToDeg(),rotated.Phi()*TMath::RadToDeg()); + } - //Calculate reconstructed t + //Calculate reconstructed t and fill histograms + + //Truth + ht_true->Fill(-1.*t_true); + htw_true->Fill(-1.*t_true,weight); //Method 1 auto t_rec1 = (e_beam - e_s_rec - pi_rec).Mag2(); + ht_rec1->Fill(-1.*t_rec1); + htw_rec1->Fill(-1.*t_rec1,weight); //Method 2 - auto p_miss = (e_beam + p_beam - e_s_rec - pi_rec); + if(zdc_cluster_x.GetSize()>0){ + auto p_miss = (e_beam + p_beam - e_s_rec - pi_rec); - TLorentzVector neut_opt; - auto neut_opt_px = p_miss.P()*sin(neut_rec.Theta())*cos(neut_rec.Phi()); - auto neut_opt_py = p_miss.P()*sin(neut_rec.Theta())*sin(neut_rec.Phi()); - auto neut_opt_pz = p_miss.P()*cos(neut_rec.Theta()); + TLorentzVector neut_opt; + auto neut_opt_px = p_miss.P()*sin(neut_rec.Theta())*cos(neut_rec.Phi()); + auto neut_opt_py = p_miss.P()*sin(neut_rec.Theta())*sin(neut_rec.Phi()); + auto neut_opt_pz = p_miss.P()*cos(neut_rec.Theta()); + + neut_opt.SetXYZM(neut_opt_px,neut_opt_py,neut_opt_pz,neut_mass); - neut_opt.SetXYZM(neut_opt_px,neut_opt_py,neut_opt_pz,neut_mass); + auto t_rec2 = (p_beam - neut_opt).Mag2(); - auto t_rec2 = (p_beam - neut_opt).Mag2(); + ht_rec2->Fill(-1.*t_rec2); + htw_rec2->Fill(-1.*t_rec2,weight); + } //Method 3 auto sum_epi = e_s_rec + pi_rec; sum_epi.RotateY(0.025); auto t_rec3 = -1.*(sum_epi).Perp2(); //Make sure t is negative - - //Method 4 - auto t_rec4 = (p_beam - neut_rec).Mag2(); - - //Fill additional histograms - ht_true->Fill(-1.*t_true); - ht_rec1->Fill(-1.*t_rec1); - ht_rec2->Fill(-1.*t_rec2); ht_rec3->Fill(-1.*t_rec3); - ht_rec4->Fill(-1.*t_rec4); + htw_rec3->Fill(-1.*t_rec3,weight); - htw_true->Fill(-1.*t_true,weight); - htw_rec1->Fill(-1.*t_rec1,weight); - htw_rec2->Fill(-1.*t_rec2,weight); - htw_rec3->Fill(-1.*t_rec3,weight); - htw_rec4->Fill(-1.*t_rec4,weight); + //Method 4 + if(zdc_cluster_x.GetSize()>0){ + auto t_rec4 = (p_beam - neut_rec).Mag2(); + ht_rec4->Fill(-1.*t_rec4); + htw_rec4->Fill(-1.*t_rec4,weight); + } } //End loop over events diff --git a/benchmarks/demp/analysis/HEXPLIT_standalone/DEMP_reco.C b/benchmarks/demp/analysis/HEXPLIT_standalone/DEMP_reco.C deleted file mode 100644 index f366095b..00000000 --- a/benchmarks/demp/analysis/HEXPLIT_standalone/DEMP_reco.C +++ /dev/null @@ -1,325 +0,0 @@ -#include "ZDC_neutron_recon.h" - -//------------------ -void DEMP_reco(){ - - //Define Style - gStyle->SetOptStat(0); - gStyle->SetPadBorderMode(0); - gStyle->SetFrameBorderMode(0); - gStyle->SetFrameLineWidth(2); - gStyle->SetLabelSize(0.035,"X"); - gStyle->SetLabelSize(0.035,"Y"); - //gStyle->SetLabelOffset(0.01,"X"); - //gStyle->SetLabelOffset(0.01,"Y"); - gStyle->SetTitleXSize(0.04); - gStyle->SetTitleXOffset(0.9); - gStyle->SetTitleYSize(0.04); - gStyle->SetTitleYOffset(0.9); - - //Define histograms - TH2 *h1_elec = new TH2D("h1_elec","Scattered electron true momentum vs. polar angle",100,110,170,100,4.6,6.8); - h1_elec->GetXaxis()->SetTitle("#theta [deg]"); h1_elec->GetXaxis()->CenterTitle(); - h1_elec->GetYaxis()->SetTitle("p [GeV/c]"); h1_elec->GetYaxis()->CenterTitle(); - - TH2 *h1_pion = new TH2D("h1_pion","#pi^{+} true momentum vs. polar angle",100,0,50,100,0,50); - h1_pion->GetXaxis()->SetTitle("#theta [deg]"); h1_pion->GetXaxis()->CenterTitle(); - h1_pion->GetYaxis()->SetTitle("p [GeV/c]"); h1_pion->GetYaxis()->CenterTitle(); - - TH2 *h1_neut = new TH2D("h1_neut","Neutron true energy vs. polar angle",100,0.4,2.4,100,40,100); - h1_neut->GetXaxis()->SetTitle("#theta [deg]"); h1_neut->GetXaxis()->CenterTitle(); - h1_neut->GetYaxis()->SetTitle("E [GeV]"); h1_neut->GetYaxis()->CenterTitle(); - - TH2 *h2_elec = new TH2D("h2_elec","Scattered electron reconstructed momentum vs. polar angle",100,110,170,100,4.6,6.8); - h2_elec->GetXaxis()->SetTitle("#theta [deg]"); h2_elec->GetXaxis()->CenterTitle(); - h2_elec->GetYaxis()->SetTitle("p [GeV/c]"); h2_elec->GetYaxis()->CenterTitle(); - - TH2 *h2_pion = new TH2D("h2_pion","#pi^{+} reconstructed momentum vs. polar angle",100,0,50,100,0,50); - h2_pion->GetXaxis()->SetTitle("#theta [deg]"); h2_pion->GetXaxis()->CenterTitle(); - h2_pion->GetYaxis()->SetTitle("p [GeV/c]"); h2_pion->GetYaxis()->CenterTitle(); - - TH2 *h2_neut = new TH2D("h2_neut","Neutron reconstructed energy vs. polar angle",100,0.4,2.4,100,40,100); - h2_neut->GetXaxis()->SetTitle("#theta [deg]"); h2_neut->GetXaxis()->CenterTitle(); - h2_neut->GetYaxis()->SetTitle("E [GeV]"); h2_neut->GetYaxis()->CenterTitle(); - - TH2 *h3_neut = new TH2D("h3_neut","Neutron truth azimuth vs. polar angle around p axis",100,0,0.6,100,-200,200); - h3_neut->GetXaxis()->SetTitle("#theta^{*} [deg]"); h3_neut->GetXaxis()->CenterTitle(); - h3_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h3_neut->GetYaxis()->CenterTitle(); - - TH2 *h4_neut = new TH2D("h4_neut","Neutron reconstructed azimuth vs. polar angle around p axis",100,0,0.6,100,-200,200); - h4_neut->GetXaxis()->SetTitle("#theta^{*} [deg]"); h4_neut->GetXaxis()->CenterTitle(); - h4_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h4_neut->GetYaxis()->CenterTitle(); - - TH1 *ht_true = new TH1D("ht_true","",100,-0.1,2); //True t - ht_true->GetXaxis()->SetTitle("-t [GeV^{2}]"); ht_true->GetXaxis()->CenterTitle(); - ht_true->GetYaxis()->SetTitle("Yield"); ht_true->GetYaxis()->CenterTitle(); - ht_true->SetLineWidth(2);ht_true->SetLineColor(kBlue); - - TH1 *ht_rec1 = new TH1D("ht_rec1","",100,-0.1,2); //Reconstructed t -- 4 vector reconstruction using only electron and pion - ht_rec1->GetXaxis()->SetTitle("-t [GeV^{2}]"); ht_rec1->GetXaxis()->CenterTitle(); - ht_rec1->GetYaxis()->SetTitle("Yield"); ht_rec1->GetYaxis()->CenterTitle(); - ht_rec1->SetLineWidth(2);ht_rec1->SetLineColor(kRed); - - TH1 *ht_rec2 = new TH1D("ht_rec2","",100,-0.1,2); //Reconstructed t -- ECCE method using electron and pion + neutron angle - ht_rec2->GetXaxis()->SetTitle("-t [GeV^{2}]"); ht_rec2->GetXaxis()->CenterTitle(); - ht_rec2->GetYaxis()->SetTitle("Yield"); ht_rec2->GetYaxis()->CenterTitle(); - ht_rec2->SetLineWidth(2);ht_rec2->SetLineColor(kGreen); - - TH1 *ht_rec3 = new TH1D("ht_rec3","",100,-0.1,2); //Reconstructed t -- Pt-based reconstruction using only electron and pion - ht_rec3->GetXaxis()->SetTitle("-t [GeV^{2}]"); ht_rec3->GetXaxis()->CenterTitle(); - ht_rec3->GetYaxis()->SetTitle("Yield"); ht_rec3->GetYaxis()->CenterTitle(); - ht_rec3->SetLineWidth(2);ht_rec3->SetLineColor(kOrange); - - TH1 *ht_rec4 = new TH1D("ht_rec4","",100,-0.1,2); //Reconstructed t -- 4 vector reconstruction using only neutron - ht_rec4->GetXaxis()->SetTitle("-t [GeV^{2}]"); ht_rec4->GetXaxis()->CenterTitle(); - ht_rec4->GetYaxis()->SetTitle("Yield"); ht_rec4->GetYaxis()->CenterTitle(); - ht_rec4->SetLineWidth(2);ht_rec4->SetLineColor(kMagenta); - - TH1 *htw_true = new TH1D("htw_true","",100,-0.1,2); //True t weighted - htw_true->GetXaxis()->SetTitle("-t [GeV^{2}]"); htw_true->GetXaxis()->CenterTitle(); - htw_true->GetYaxis()->SetTitle("Yield [arb.]"); htw_true->GetYaxis()->CenterTitle(); - htw_true->SetLineWidth(2);htw_true->SetLineColor(kBlue); - - TH1 *htw_rec1 = new TH1D("htw_rec1","",100,-0.1,2); //Reconstructed t weighted - htw_rec1->GetXaxis()->SetTitle("-t [GeV^{2}]"); htw_rec1->GetXaxis()->CenterTitle(); - htw_rec1->GetYaxis()->SetTitle("Yield [arb.]"); htw_rec1->GetYaxis()->CenterTitle(); - htw_rec1->SetLineWidth(2);htw_rec1->SetLineColor(kRed); - - TH1 *htw_rec2 = new TH1D("htw_rec2","",100,-0.1,2); //Reconstructed t weighted - htw_rec2->GetXaxis()->SetTitle("-t [GeV^{2}]"); htw_rec2->GetXaxis()->CenterTitle(); - htw_rec2->GetYaxis()->SetTitle("Yield [arb.]"); htw_rec2->GetYaxis()->CenterTitle(); - htw_rec2->SetLineWidth(2);htw_rec2->SetLineColor(kGreen); - - TH1 *htw_rec3 = new TH1D("htw_rec3","",100,-0.1,2); //Reconstructed t - htw_rec3->GetXaxis()->SetTitle("-t [GeV^{2}]"); htw_rec3->GetXaxis()->CenterTitle(); - htw_rec3->GetYaxis()->SetTitle("Yield"); htw_rec3->GetYaxis()->CenterTitle(); - htw_rec3->SetLineWidth(2);htw_rec3->SetLineColor(kOrange); - - TH1 *htw_rec4 = new TH1D("htw_rec4","",100,-0.1,2); //Reconstructed t - htw_rec4->GetXaxis()->SetTitle("-t [GeV^{2}]"); htw_rec4->GetXaxis()->CenterTitle(); - htw_rec4->GetYaxis()->SetTitle("Yield"); htw_rec4->GetYaxis()->CenterTitle(); - htw_rec4->SetLineWidth(2);htw_rec4->SetLineColor(kMagenta); - - TFile *f = new TFile("eicrecon_out.root"); - TTree *tree = (TTree*) f->Get("events"); - - //Create Array Reader - TTreeReader tr(tree); - - TTreeReaderArray gen_status(tr, "MCParticles.generatorStatus"); - TTreeReaderArray gen_pid(tr, "MCParticles.PDG"); - TTreeReaderArray gen_px(tr, "MCParticles.momentum.x"); - TTreeReaderArray gen_py(tr, "MCParticles.momentum.y"); - TTreeReaderArray gen_pz(tr, "MCParticles.momentum.z"); - TTreeReaderArray gen_mass(tr, "MCParticles.mass"); - TTreeReaderArray gen_charge(tr, "MCParticles.charge"); - TTreeReaderArray gen_vx(tr, "MCParticles.vertex.x"); - TTreeReaderArray gen_vy(tr, "MCParticles.vertex.y"); - TTreeReaderArray gen_vz(tr, "MCParticles.vertex.z"); - - TTreeReaderArray rec_pid(tr, "ReconstructedChargedParticles.PDG"); - TTreeReaderArray rec_px(tr, "ReconstructedChargedParticles.momentum.x"); - TTreeReaderArray rec_py(tr, "ReconstructedChargedParticles.momentum.y"); - TTreeReaderArray rec_pz(tr, "ReconstructedChargedParticles.momentum.z"); - TTreeReaderArray rec_mass(tr, "ReconstructedChargedParticles.mass"); - - //local positions and energies in the ZDC - TTreeReaderArray zdc_hit_energies(tr, "ZDCRecHits.energy"); - TTreeReaderArray zdc_hit_times(tr, "ZDCRecHits.time"); - TTreeReaderArray zdc_hit_local_x(tr, "ZDCRecHits.local.x"); - TTreeReaderArray zdc_hit_local_y(tr, "ZDCRecHits.local.y"); - TTreeReaderArray zdc_hit_local_z(tr, "ZDCRecHits.local.z"); - - TTreeReaderArray>> weight_map(tr,"_stringMap"); - - //Other variables - TLorentzVector gen_vec; - TVector3 gen_vertex; - TLorentzVector rec_vec; - int counter(0); - - TLorentzVector p_beam_real; //Proton beam (incorporating event-by-event fluctuations) for true t calculation - TLorentzVector neut_true; //True neutron for true t calculation - TLorentzVector e_beam(0.,0.,-5.,5.); //Electron beam for t reconstruction (ignoring event-by-event fluctuations) - TLorentzVector p_beam(100.*sin(-25./1000),0,100*cos(-25./1000),100); //Proton beam for t reconstruction (ignoring event-by-event fluctuations) - TLorentzVector e_s_rec; //Reconstructed scattered electron from tracking detector - TLorentzVector pi_rec; //Reconstructed pi+ from tracking detector - - //Loop over events - while (tr.Next()) { - - if(counter%100==0) cout<<"Analyzing event "< weight_value = weight_map[0].second; - - double weight = std::stod( *(weight_value.begin()) ); - //cout<<"Weight = "<Fill(gen_vec.Theta()*TMath::RadToDeg(),gen_vec.P()); - if(gen_pid[igen]==211) h1_pion->Fill(gen_vec.Theta()*TMath::RadToDeg(),gen_vec.P()); - - if(gen_pid[igen]==2112){ - h1_neut->Fill(gen_vec.Theta()*TMath::RadToDeg(),gen_vec.E()); - neut_true = gen_vec; - - //W.r.t proton beam direction - TLorentzVector rotated = gen_vec; - rotated.RotateY(0.025); - h3_neut->Fill(rotated.Theta()*TMath::RadToDeg(),rotated.Phi()*TMath::RadToDeg()); - - } - - } - } //End loop over generated particles - - //Calculate true t - auto t_true = (p_beam_real - neut_true).Mag2(); - - //Loop over reconstructed tracks for electron and pion - for(int irec=0;irecFill(rec_vec.Theta()*TMath::RadToDeg(),rec_vec.P()); - e_s_rec = rec_vec; - } - - if(rec_pid[irec]==211){ - h2_pion->Fill(rec_vec.Theta()*TMath::RadToDeg(),rec_vec.P()); - pi_rec = rec_vec; - } - - }//End loop over reconstructed particles - - //Neutron Reconstruction using HEXSPLIT algorithm - TLorentzVector neut_rec = ZDC_neutron_recon(zdc_hit_energies, zdc_hit_times, - zdc_hit_local_x, zdc_hit_local_y,zdc_hit_local_z); - - h2_neut->Fill(neut_rec.Theta()*TMath::RadToDeg(),neut_rec.E()); - - //W.r.t proton beam direction - TLorentzVector rotated = neut_rec; - rotated.RotateY(0.025); - h4_neut->Fill(rotated.Theta()*TMath::RadToDeg(),rotated.Phi()*TMath::RadToDeg()); - - //Calculate reconstructed t - - //Method 1 - auto t_rec1 = (e_beam - e_s_rec - pi_rec).Mag2(); - - //Method 2 - auto p_miss = (e_beam + p_beam - e_s_rec - pi_rec); - - TLorentzVector neut_opt; - auto neut_opt_px = p_miss.P()*sin(neut_rec.Theta())*cos(neut_rec.Phi()); - auto neut_opt_py = p_miss.P()*sin(neut_rec.Theta())*sin(neut_rec.Phi()); - auto neut_opt_pz = p_miss.P()*cos(neut_rec.Theta()); - double neut_mass = 0.9396; - - neut_opt.SetXYZM(neut_opt_px,neut_opt_py,neut_opt_pz,neut_mass); - - auto t_rec2 = (p_beam - neut_opt).Mag2(); - - //Method 3 - auto sum_epi = e_s_rec + pi_rec; - sum_epi.RotateY(0.025); - auto t_rec3 = -1.*(sum_epi).Perp2(); //Make sure t is negative - - //Method 4 - auto t_rec4 = (p_beam - neut_rec).Mag2(); - - //Fill additional histograms - ht_true->Fill(-1.*t_true); - ht_rec1->Fill(-1.*t_rec1); - ht_rec2->Fill(-1.*t_rec2); - ht_rec3->Fill(-1.*t_rec3); - ht_rec4->Fill(-1.*t_rec4); - - htw_true->Fill(-1.*t_true,weight); - htw_rec1->Fill(-1.*t_rec1,weight); - htw_rec2->Fill(-1.*t_rec2,weight); - htw_rec3->Fill(-1.*t_rec3,weight); - htw_rec4->Fill(-1.*t_rec4,weight); - - } //End loop over events - - //Make plots - TCanvas *c1 = new TCanvas("c1"); - h1_elec->Draw("colz"); - - TCanvas *c2 = new TCanvas("c2"); - h1_pion->Draw("colz"); - - TCanvas *c3 = new TCanvas("c3"); - h1_neut->Draw("colz"); - - TCanvas *c4 = new TCanvas("c4"); - h2_elec->Draw("colz"); - - TCanvas *c5 = new TCanvas("c5"); - h2_pion->Draw("colz"); - - TCanvas *c6 = new TCanvas("c6"); - h2_neut->Draw("colz"); - - TCanvas *c7 = new TCanvas("c7"); - h3_neut->Draw("colz"); - - TCanvas *c8 = new TCanvas("c8"); - h4_neut->Draw("colz"); - - TCanvas *c9 = new TCanvas("c9"); - ht_true->Draw("hist"); - ht_rec1->Draw("hist same"); - ht_rec2->Draw("hist same"); - ht_rec3->Draw("hist same"); - ht_rec4->Draw("hist same"); - - TLegend *leg1 = new TLegend(0.5,0.6,0.85,0.85); - leg1->SetBorderSize(0);leg1->SetTextSize(0.03); - leg1->SetFillStyle(0); - leg1->AddEntry(ht_true,"Truth (after beam-effects)","l"); - leg1->AddEntry(ht_rec1,"Electron+Pion 4-vector","l"); - leg1->AddEntry(ht_rec4,"Neutron only","l"); - leg1->AddEntry(ht_rec3,"Electron+Pion P_{T}-based","l"); - leg1->AddEntry(ht_rec2,"ECCE paper method","l"); - leg1->Draw(); - - TCanvas *c10 = new TCanvas("c10"); - htw_true->Draw("hist"); - htw_rec1->Draw("hist same"); - htw_rec2->Draw("hist same"); - htw_rec3->Draw("hist same"); - htw_rec4->Draw("hist same"); - leg1->Draw(); - - //Print plots to file - c1->Print("DEMP_reco.pdf["); - c1->Print("DEMP_reco.pdf"); - c2->Print("DEMP_reco.pdf"); - c3->Print("DEMP_reco.pdf"); - c4->Print("DEMP_reco.pdf"); - c5->Print("DEMP_reco.pdf"); - c6->Print("DEMP_reco.pdf"); - c7->Print("DEMP_reco.pdf"); - c8->Print("DEMP_reco.pdf"); - c9->Print("DEMP_reco.pdf"); - c10->Print("DEMP_reco.pdf"); - c10->Print("DEMP_reco.pdf]"); - -} diff --git a/benchmarks/demp/analysis/HEXPLIT_standalone/ZDC_neutron_recon.h b/benchmarks/demp/analysis/HEXPLIT_standalone/ZDC_neutron_recon.h deleted file mode 100644 index 889040d7..00000000 --- a/benchmarks/demp/analysis/HEXPLIT_standalone/ZDC_neutron_recon.h +++ /dev/null @@ -1,65 +0,0 @@ -#include "hexplit.h" -TLorentzVector ZDC_neutron_recon(TTreeReaderArray&E,TTreeReaderArray&t,TTreeReaderArray&x,TTreeReaderArray&y, - TTreeReaderArray&z){ - //cout << "hit count" << E.GetSize() << endl; - vector subcellE, subcellx, subcelly, subcellz; - hexplit_opts opts1; - opts1.MIP=0.000470; - opts1.side_length=31.3; - opts1.layer_spacing=25.1; - opts1.Emin=opts1.MIP*0.1; - opts1.tmax=320; - - subcell_reweight(E, t, x, y, z, subcellE, subcellx, subcelly, subcellz, opts1); - - //now determine the position from the subcells - position_recon_opts opts2; - opts2.sampling_fraction=0.0203; - - - /*vector flatE;flatE.reserve(2048*10); - vector flatx;flatx.reserve(2048*10); - vector flaty;flaty.reserve(2048*10); - vector flatz;flatz.reserve(2048*10); - - flatten(subcellE, subcellx, subcelly, subcellz, flatE, flatx, flaty, flatz); - subcellE=flatE; - subcellx=flatx; - subcelly=flaty; - subcellz=flatz;*/ - - - // parameterization of values to determine which - // value of w0 to use - opts2.w0_a=5.0; - opts2.w0_b=0.65; - opts2.w0_c=0.31; - opts2.E0=50; - - /*opts2.w0_a=3.0; - opts2.w0_b=0; - opts2.w0_c=0;*/ - - //obtains the local position - double E_recon, x_recon, y_recon, z_recon; - - recon_position(subcellE, subcellx, subcelly, subcellz, E_recon, x_recon, - y_recon, z_recon, opts2); - - //translate and then rotate - double z_offset=36601; //z' position of the center of the ZDC SiPM-on-tile detector - double cross=-0.025; - z_recon=z_recon+z_offset; - double z_recon_new=z_recon*cos(cross)-x_recon*sin(cross); - double x_recon_new=x_recon*cos(cross)+z_recon*sin(cross); - - x_recon=x_recon_new; - z_recon=z_recon_new; - double r_recon=sqrt(x_recon*x_recon+y_recon*y_recon+z_recon*z_recon); - double m_n=0.939565; - double p_recon=sqrt(E_recon*E_recon-m_n*m_n); - //cout << "x, y, z recon " <SetTextColor(kBlack); - title->SetNDC(); - title->SetTextAlign(22); - //title->SetTextSize(0.05); // Adjust the size as needed - title->Draw(); - - gPad->Update(); -} -//------------------ -void demp_tests(const char* fname = "rec_demp.root"){ - - //Define Style - gStyle->SetOptStat(0); - //titles in ROOT don't do what you want them to do. - //use TLatex instead - gStyle->SetOptTitle(0); - gStyle->SetPadTopMargin(0.1); - gStyle->SetOptFit(1); - //gStyle->SetTitleSize(1,"T"); - //gStyle->SetTitleOffset(1, "T"); - //gStyle->SetTitleSize(14, "T"); - /*gStyle->SetPadBorderMode(0); - gStyle->SetFrameBorderMode(0); - gStyle->SetFrameLineWidth(2); - gStyle->SetLabelSize(0.035,"X"); - gStyle->SetLabelSize(0.035,"Y"); - //gStyle->SetLabelOffset(0.01,"X"); - //gStyle->SetLabelOffset(0.01,"Y"); - - gStyle->SetTitleXSize(0.04); - gStyle->SetTitleXOffset(0.9); - gStyle->SetTitleYSize(0.04); - gStyle->SetTitleYOffset(0.9);//*/ - - //Define histograms - TH2 *h1_elec = new TH2D("h1_elec","Scattered electron true momentum vs polar angle",100,110,170,100,4.6,6.8); - h1_elec->GetXaxis()->SetTitle("#theta [deg]"); h1_elec->GetXaxis()->CenterTitle(); - h1_elec->GetYaxis()->SetTitle("p [GeV/c]"); h1_elec->GetYaxis()->CenterTitle(); - - TH2 *h1_pion = new TH2D("h1_pion","#pi^{+} true momentum vs. polar angle",100,0,50,100,0,50); - h1_pion->GetXaxis()->SetTitle("#theta [deg]"); h1_pion->GetXaxis()->CenterTitle(); - h1_pion->GetYaxis()->SetTitle("p [GeV/c]"); h1_pion->GetYaxis()->CenterTitle(); - - TH2 *h1_neut = new TH2D("h1_neut","Neutron true momentum vs. polar angle",100,0.4,2.4,100,40,100); - h1_neut->GetXaxis()->SetTitle("#theta [deg]"); h1_neut->GetXaxis()->CenterTitle(); - h1_neut->GetYaxis()->SetTitle("p [GeV/c]"); h1_neut->GetYaxis()->CenterTitle(); - - TH2 *h2_elec = new TH2D("h2_elec","Scattered electron reconstructed momentum vs. polar angle",100,110,170,100,4.6,6.8); - h2_elec->GetXaxis()->SetTitle("#theta [deg]"); h2_elec->GetXaxis()->CenterTitle(); - h2_elec->GetYaxis()->SetTitle("p [GeV/c]"); h2_elec->GetYaxis()->CenterTitle(); - - TH2 *h2_pion = new TH2D("h2_pion","#pi^{+} reconstructed momentum vs. polar angle",100,0,50,100,0,50); - h2_pion->GetXaxis()->SetTitle("#theta [deg]"); h2_pion->GetXaxis()->CenterTitle(); - h2_pion->GetYaxis()->SetTitle("p [GeV/c]"); h2_pion->GetYaxis()->CenterTitle(); - - TH2 *h2_neut = new TH2D("h2_neut","Neutron reconstructed momentum vs. polar angle",100,0.4,2.4,100,40,100); - h2_neut->GetXaxis()->SetTitle("#theta [deg]"); h2_neut->GetXaxis()->CenterTitle(); - h2_neut->GetYaxis()->SetTitle("p [GeV/c]"); h2_neut->GetYaxis()->CenterTitle(); - - TH1 *h3_neut = new TH1D("h3_neut","Neutron theta resolution",100,-1,1); - h3_neut->GetXaxis()->SetTitle("#Delta#theta [mrad]"); h3_neut->GetXaxis()->CenterTitle(); - h3_neut->GetYaxis()->SetTitle("events"); h3_neut->GetYaxis()->CenterTitle(); - - TH1 *h4_neut = new TH1D("h4_neut","Neutron energy resolution",100,-30,10); - h4_neut->GetXaxis()->SetTitle("#Delta E/E [%]"); h4_neut->GetXaxis()->CenterTitle(); - h4_neut->GetYaxis()->SetTitle("events"); h4_neut->GetYaxis()->CenterTitle(); - - TH2 *h5_neut = new TH2D("h5_neut","Neutron truth azimuth vs. polar angle",100,0.4,2.4,100,160,200); - h5_neut->GetXaxis()->SetTitle("#theta [deg]"); h5_neut->GetXaxis()->CenterTitle(); - h5_neut->GetYaxis()->SetTitle("#phi [deg]"); h5_neut->GetYaxis()->CenterTitle(); - - TH2 *h6_neut = new TH2D("h6_neut","Neutron reconstructed azimuth vs. polar angle",100,0.4,2.4,100,160,200); - h6_neut->GetXaxis()->SetTitle("#theta [deg]"); h6_neut->GetXaxis()->CenterTitle(); - h6_neut->GetYaxis()->SetTitle("#phi [deg]"); h6_neut->GetYaxis()->CenterTitle(); - - TH2 *h7_neut = new TH2D("h7_neut","Neutron truth azimuth vs. polar angle around p axis",100,0,0.5,100,-180,180); - h7_neut->GetXaxis()->SetTitle("#theta^{*} [deg]"); h7_neut->GetXaxis()->CenterTitle(); - h7_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h7_neut->GetYaxis()->CenterTitle(); - - TH2 *h8_neut = new TH2D("h8_neut","Neutron reconstructed azimuth vs. polar angle around p axis",100,0,0.5,100,-180,180); - h8_neut->GetXaxis()->SetTitle("#theta^{*} [deg]"); h8_neut->GetXaxis()->CenterTitle(); - h8_neut->GetYaxis()->SetTitle("#phi^{*} [deg]"); h8_neut->GetYaxis()->CenterTitle(); - - - TFile *f = new TFile(fname); - TTree *tree = (TTree*) f->Get("events"); - - //Create Array Reader - TTreeReader tr(tree); - - TTreeReaderArray gen_status(tr, "MCParticles.generatorStatus"); - TTreeReaderArray gen_pid(tr, "MCParticles.PDG"); - TTreeReaderArray gen_px(tr, "MCParticles.momentum.x"); - TTreeReaderArray gen_py(tr, "MCParticles.momentum.y"); - TTreeReaderArray gen_pz(tr, "MCParticles.momentum.z"); - TTreeReaderArray gen_mass(tr, "MCParticles.mass"); - TTreeReaderArray gen_charge(tr, "MCParticles.charge"); - TTreeReaderArray gen_vx(tr, "MCParticles.vertex.x"); - TTreeReaderArray gen_vy(tr, "MCParticles.vertex.y"); - TTreeReaderArray gen_vz(tr, "MCParticles.vertex.z"); - - TTreeReaderArray rec_pid(tr, "ReconstructedChargedParticles.PDG"); - TTreeReaderArray rec_px(tr, "ReconstructedChargedParticles.momentum.x"); - TTreeReaderArray rec_py(tr, "ReconstructedChargedParticles.momentum.y"); - TTreeReaderArray rec_pz(tr, "ReconstructedChargedParticles.momentum.z"); - TTreeReaderArray rec_mass(tr, "ReconstructedChargedParticles.mass"); - - //local positions and energies in the ZDC - TTreeReaderArray zdc_hit_energies(tr, "ZDCRecHits.energy"); - TTreeReaderArray zdc_hit_times(tr, "ZDCRecHits.time"); - TTreeReaderArray zdc_hit_local_x(tr, "ZDCRecHits.local.x"); - TTreeReaderArray zdc_hit_local_y(tr, "ZDCRecHits.local.y"); - TTreeReaderArray zdc_hit_local_z(tr, "ZDCRecHits.local.z"); - - //Other variables - TLorentzVector gen_vec; - TVector3 gen_vertex; - TLorentzVector rec_vec; - int counter(0); - - //Loop over events - while (tr.Next()) { - - if(counter%100==0) cout<<"Analyzing event "<Fill(gen_vec.Theta()*TMath::RadToDeg(),gen_vec.P()); - if(gen_pid[igen]==211) h1_pion->Fill(gen_vec.Theta()*TMath::RadToDeg(),gen_vec.P()); - if(gen_pid[igen]==2112) { - h1_neut->Fill(gen_vec.Theta()*TMath::RadToDeg(),gen_vec.P()); - h5_neut->Fill(gen_vec.Theta()*TMath::RadToDeg(),(gen_vec.Phi()+(gen_vec.Phi()<0)*2*M_PI)*TMath::RadToDeg()); - theta_neut_truth=gen_vec.Theta(); - E_neut_truth=gen_vec.E(); - - TLorentzVector rotated = gen_vec; - rotated.RotateY(0.025); - h7_neut->Fill(rotated.Theta()*TMath::RadToDeg(),rotated.Phi()*TMath::RadToDeg()); - } - - } - } //End loop over generated particles - - //Loop over reconstructed tracks for electron and pion - for(int irec=0;irecFill(rec_vec.Theta()*TMath::RadToDeg(),rec_vec.P()); - if(rec_pid[irec]==211) h2_pion->Fill(rec_vec.Theta()*TMath::RadToDeg(),rec_vec.P()); - - }//End loop over reconstructed particles - - TLorentzVector neutron = ZDC_neutron_recon(zdc_hit_energies, zdc_hit_times, zdc_hit_local_x, zdc_hit_local_y, - zdc_hit_local_z); - //cout << "neutron theta " << neutron.Theta() << endl; - h2_neut->Fill(neutron.Theta()*TMath::RadToDeg(),neutron.P()); - h3_neut->Fill((neutron.Theta()-theta_neut_truth)*1000); - h4_neut->Fill((neutron.E()/E_neut_truth-1)*100); - h6_neut->Fill(neutron.Theta()*TMath::RadToDeg(),(neutron.Phi()+(neutron.Phi()<0)*2*M_PI)*TMath::RadToDeg()); - - TLorentzVector rotated = neutron; - rotated.RotateY(0.025); - h8_neut->Fill(rotated.Theta()*TMath::RadToDeg(),rotated.Phi()*TMath::RadToDeg()); - - } //End loop over events - - //Make plots - int w=800, h=600; - TCanvas *c1 = new TCanvas("c1", "", w,h); - h1_elec->Draw("colz"); - draw_title(h1_elec->GetTitle()); - - TCanvas *c2 = new TCanvas("c2", "", w,h); - h1_pion->Draw("colz"); - draw_title(h1_pion->GetTitle()); - - TCanvas *c3 = new TCanvas("c3", "", w,h); - h1_neut->Draw("colz"); - draw_title(h1_neut->GetTitle()); - cout << "generated neutron count " << h1_neut->Integral(); - - TCanvas *c4 = new TCanvas("c4", "", w,h); - h2_elec->Draw("colz"); - draw_title(h2_elec->GetTitle()); - - TCanvas *c5 = new TCanvas("c5", "", w,h); - h2_pion->Draw("colz"); - draw_title(h2_pion->GetTitle()); - - TCanvas *c6 = new TCanvas("c6", "", w,h); - h2_neut->Draw("colz"); - draw_title(h2_neut->GetTitle()); - cout << "reconstructed neutron count " << h2_neut->Integral(); - - - TCanvas *c7 = new TCanvas("c7", "", w,h); - TF1* fnc = new TF1("f1", "gaus(0)", -1,1); - h3_neut->Draw("colz"); - h3_neut->Fit(fnc); - draw_title(h3_neut->GetTitle()); - - TCanvas *c8 = new TCanvas("c8", "", w,h); - h4_neut->Draw("colz"); - h4_neut->Fit(fnc); - draw_title(h4_neut->GetTitle()); - - TCanvas *c9 = new TCanvas("c9", "", w,h); - h5_neut->Draw("colz"); - draw_title(h5_neut->GetTitle()); - - TCanvas *c10 = new TCanvas("c10", "", w,h); - h6_neut->Draw("colz"); - draw_title(h6_neut->GetTitle()); - - TCanvas *c11 = new TCanvas("c11", "", w,h); - h7_neut->Draw("colz"); - draw_title(h7_neut->GetTitle()); - - TCanvas *c12 = new TCanvas("c12", "", w,h); - h8_neut->Draw("colz"); - draw_title(h8_neut->GetTitle()); - - //Print plots to file - c1->Print("results/demp/DEMP_reco.pdf["); - c1->Print("results/demp/DEMP_reco.pdf"); - c2->Print("results/demp/DEMP_reco.pdf"); - c3->Print("results/demp/DEMP_reco.pdf"); - c4->Print("results/demp/DEMP_reco.pdf"); - c5->Print("results/demp/DEMP_reco.pdf"); - c6->Print("results/demp/DEMP_reco.pdf"); - c7->Print("results/demp/DEMP_reco.pdf"); - c8->Print("results/demp/DEMP_reco.pdf"); - c9->Print("results/demp/DEMP_reco.pdf"); - c10->Print("results/demp/DEMP_reco.pdf"); - c11->Print("results/demp/DEMP_reco.pdf"); - c12->Print("results/demp/DEMP_reco.pdf"); - c12->Print("results/demp/DEMP_reco.pdf]"); -} diff --git a/benchmarks/demp/analysis/HEXPLIT_standalone/hexplit.h b/benchmarks/demp/analysis/HEXPLIT_standalone/hexplit.h deleted file mode 100644 index 3b370567..00000000 --- a/benchmarks/demp/analysis/HEXPLIT_standalone/hexplit.h +++ /dev/null @@ -1,190 +0,0 @@ -#include -//#include -#define SUBCELLS 12 -struct hexplit_opts { - // longitudinal space between layers - double layer_spacing; - // side length of a cell - double side_length; - // MIP value for the detector - double MIP; - // minimum energy cut - double Emin=0; - // maximum time cut - double tmax=1000; -}; -struct position_recon_opts{ - //use the sampling fraction to get the recon energy, - // which is used in the parameterization for the w0. - double sampling_fraction; - - // parameterization of values to determine which - // value of w0 to use - double w0_a; - double w0_b; - double w0_c; - double E0=50; - -}; -const double pi=M_PI; -/* -phi=np.linspace(0, np.pi*5/3, 6) - cph=np.cos(phi) - sph=np.sin(phi) - - offsetsx=np.concatenate((1.5*cph, -sqrt3/2*sph))*sl - offsetsy=np.concatenate((1.5*sph, sqrt3/2*cph))*sl - - #determine the positions of the vertices - subcell_offsetsx=[] - subcell_offsetsy=[] - subcell_offsetsx+=[[sl/2*cph[k], sl/2*(cph[k]+cph[(k+1)%6]), - sl*cph[k], sl/2*(cph[k]+cph[(k+5)%6])] for k in range(6)] - subcell_offsetsx+=[[0, sl/4*cph[k]-sqrt3*sl*sph[k]/4, - -sqrt3*sl*sph[k]/2, -sl/4*cph[k]-sqrt3*sl*sph[k]/4] for k in range(6)] - - subcell_offsetsy+=[[sl/2*sph[k], sl/2*(sph[k]+sph[(k+1)%6]), - sl*sph[k], sl/2*(sph[k]+sph[(k+5)%6])] for k in range(6)] - subcell_offsetsy+=[[0, sqrt3*sl*cph[k]/4+sl/4*sph[k], - sqrt3*sl*cph[k]/2, sqrt3*sl*cph[k]/4-sl/4*sph[k]] for k in range(6)]*/ - -//positions where the overlapping cells are relative to a given cell (in units of hexagon side length) -double neighbor_offsets_x[]={1.5*cos(0), 1.5*cos(pi/3), 1.5*cos(2*pi/3),1.5*cos(3*pi/3), 1.5*cos(4*pi/3), 1.5*cos(5*pi/3), - -sqrt(3)/2.*sin(0),-sqrt(3)/2.*sin(pi/3),-sqrt(3)/2.*sin(2*pi/3),-sqrt(3)/2.*sin(3*pi/3),-sqrt(3)/2.*sin(4*pi/3),-sqrt(3)/2.*sin(5*pi/3)}; -double neighbor_offsets_y[]={1.5*sin(0), 1.5*sin(pi/3), 1.5*sin(2*pi/3),1.5*sin(3*pi/3), 1.5*sin(4*pi/3), 1.5*sin(5*pi/3), - sqrt(3)/2.*cos(0), sqrt(3)/2.*cos(pi/3), sqrt(3)/2.*cos(2*pi/3), sqrt(3)/2.*cos(3*pi/3), sqrt(3)/2.*cos(4*pi/3), sqrt(3)/2.*cos(5*pi/3)}; - -//indices of the neighboring cells which overlap to produce a given subcell -int neighbor_indices[12][3]={{0, 11,10}, {1, 6, 11},{2, 7, 6}, {3,8,7}, {4,9,8}, {5,10,9}, - {6, 11, 7}, {7, 6, 8}, {8, 7, 9}, {9,8,10},{10,9,11},{11,10,6}}; - -//positions of the centers of subcells -double subcell_offsets_x[]={0.75*cos(0), 0.75*cos(pi/3), 0.75*cos(2*pi/3), 0.75*cos(3*pi/3), 0.75*cos(4*pi/3), 0.75*cos(5*pi/3), - -sqrt(3)/4*cos(0),-sqrt(3)/4*cos(pi/3),-sqrt(3)/4*cos(2*pi/3),-sqrt(3)/4*cos(3*pi/3),-sqrt(3)/4*cos(4*pi/3),-sqrt(3)/4*cos(5*pi/3)}; -double subcell_offsets_y[]={0.75*sin(0), 0.75*sin(pi/3), 0.75*sin(2*pi/3), 0.75*sin(3*pi/3), 0.75*sin(4*pi/3), 0.75*sin(5*pi/3), - sqrt(3)/4*sin(0), sqrt(3)/4*sin(pi/3), sqrt(3)/4*sin(2*pi/3), sqrt(3)/4*sin(3*pi/3), sqrt(3)/4*sin(4*pi/3), sqrt(3)/4*sin(5*pi/3)}; - -void subcell_reweight(TTreeReaderArray&E,TTreeReaderArray&t,TTreeReaderArray&x,TTreeReaderArray&y,TTreeReaderArray&z, //number of hits, energies, time, and local x,y,z positions of the hits - vector & subcellE, vector& subcellx, vector & subcelly, vector & subcellz, //returned position of the cluster - hexplit_opts opts) { - int nhits=E.GetSize(); - double sl=opts.side_length; - double layer_spacing=opts.layer_spacing; - double Emin=opts.Emin; - double tmax=opts.tmax; - double MIP=opts.MIP; - - double Esum=0; - for(int i=0; itmax) - continue; - - //keep track of the energy in each neighboring cell - double Eneighbors[SUBCELLS]; - for (int j=0; j2.5*layer_spacing || z[i]==z[j]) - continue; - //cout << "z cut" << endl; - if (E[j]tmax) - continue; - //difference in transverse position (in units of side lengths) - double dx=(x[j]-x[i])/sl; - double dy=(y[j]-y[i])/sl; - if (abs(dx)>2 || abs(dy)>sqrt(3)) - continue; - //cout << "dx, dy cut passed: " << dx << " " << dy << endl; - - //loop over locations of the neighboring cells - //and check if the jth hit matches this location - double tol=0.01; //tolerance for rounding errors - for(int k=0;k&E, vector& x, vector& y, vector& z, - vector&Enew,vector& xnew, vector& ynew, vector& znew){ - Enew.clear(); - xnew.clear(); - ynew.clear(); - znew.clear(); - for(int i = 0; i& E, vector& x, vector& y, vector& z, double & E_recon, double & x_recon, - double & y_recon, double & z_recon, position_recon_opts& opts){ - - //first determine the total energy of the particle and of the shower - E_recon=0; - double E_shower=0; - for (int i = 0; i