Skip to content

Commit

Permalink
Update
Browse files Browse the repository at this point in the history
  • Loading branch information
@bschmookler
bschmookler committed Jan 26, 2024
1 parent 3453e22 commit fb5b278
Show file tree
Hide file tree
Showing 5 changed files with 52 additions and 869 deletions.
92 changes: 52 additions & 40 deletions ...emp/analysis/HEXPLIT_EICRecon/DEMP_reco.C → benchmarks/demp/analysis/DEMP_reco.C
View file
Original file line number Diff line number Diff line change
Expand Up @@ -122,10 +122,10 @@ void DEMP_reco(){
TTreeReaderArray<float> rec_mass(tr, "ReconstructedChargedParticles.mass");

//HEXPLIT Clusters
TTreeReaderArray<float> zdc_cluster_energy(tr, "ZDC_HEXPLITClusters.energy");
TTreeReaderArray<float> zdc_cluster_x(tr, "ZDC_HEXPLITClusters.position.x");
TTreeReaderArray<float> zdc_cluster_y(tr, "ZDC_HEXPLITClusters.position.y");
TTreeReaderArray<float> zdc_cluster_z(tr, "ZDC_HEXPLITClusters.position.z");
TTreeReaderArray<float> zdc_cluster_energy(tr, "HcalFarForwardZDCClusters.energy");
TTreeReaderArray<float> zdc_cluster_x(tr, "HcalFarForwardZDCClusters.position.x");
TTreeReaderArray<float> zdc_cluster_y(tr, "HcalFarForwardZDCClusters.position.y");
TTreeReaderArray<float> zdc_cluster_z(tr, "HcalFarForwardZDCClusters.position.z");

TTreeReaderArray<pair<string,vector<string>>> weight_map(tr,"_stringMap");

Expand Down Expand Up @@ -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

Expand Down
Loading

0 comments on commit fb5b278

Please sign in to comment.