Commit 4892e0db authored by Wouter Deconinck's avatar Wouter Deconinck
Browse files

mcID has value and source

parent ffade60c
......@@ -67,8 +67,8 @@ public:
bool ebeam_found = false;
bool pbeam_found = false;
eic::Index scatID;
int32_t mcscatID = -1;
for (const auto& p : parts) {
if (p.genStatus() == 4 && p.pdgID() == 11) { // Incoming electron
......@@ -89,15 +89,15 @@ public:
// which seems to be correct based on a cursory glance at the Pythia8 output. In the future,
// it may be better to trace back each final-state electron and see which one originates from
// the beam.
if (p.genStatus() == 1 && p.pdgID() == 11 && !scatID) {
if (p.genStatus() == 1 && p.pdgID() == 11 && mcscatID == -1) {
ef.setPx(p.ps().x);
ef.setPy(p.ps().y);
ef.setPz(p.ps().z);
ef.setE(p.energy());
scatID = p.ID();
mcscatID = p.ID();
}
if (ebeam_found && pbeam_found && scatID) {
if (ebeam_found && pbeam_found && mcscatID != -1) {
// all done!
break;
}
......@@ -116,7 +116,7 @@ public:
}
return StatusCode::SUCCESS;
}
if (scatID == 0) {
if (mcscatID == -1) {
if (msgLevel(MSG::DEBUG)) {
debug() << "No scattered electron found" << endmsg;
}
......@@ -131,7 +131,7 @@ public:
kin.nu(q * pi / m_proton);
kin.x(kin.Q2() / (2. * q * pi));
kin.W(sqrt((pi + q).m2()));
kin.scatID(scatID);
kin.scatID(mcscatID);
// Debugging output
if (msgLevel(MSG::DEBUG)) {
......
......@@ -77,7 +77,7 @@ public:
bool found_proton = false;
//Also get the true scattered electron, which will not be included in the sum
//over final-state particles for the JB reconstruction
eic::Index mcscatID;
int32_t mcscatID = -1;
for (const auto& p : mcparts) {
if (p.genStatus() == 4 && p.pdgID() == 11) {
// Incoming electron
......@@ -131,11 +131,11 @@ public:
found_proton = true;
}
// Index of true Scattered electron. Currently taken as first status==1 electron in HEPMC record.
if (p.genStatus() == 1 && p.pdgID() == 11 && !mcscatID) {
if (p.genStatus() == 1 && p.pdgID() == 11 && mcscatID == -1) {
mcscatID = p.ID();
}
if (found_electron && found_proton && mcscatID) {
if (found_electron && found_proton && mcscatID != -1) {
break;
}
}
......@@ -152,6 +152,12 @@ public:
}
return StatusCode::SUCCESS;
}
if (mcscatID == -1) {
if (msgLevel(MSG::DEBUG)) {
debug() << "No truth scattered electron found" << endmsg;
}
return StatusCode::SUCCESS;
}
// Loop over reconstructed particles to get all outgoing particles
// -----------------------------------------------------------------
......@@ -183,7 +189,7 @@ public:
double E_boosted = p.energy() + (-m_crossingAngle)/2.*p.p().x;
//Get the scattered electron index and angle
if(p.mcID() == mcscatID){
if(p.mcID().value == mcscatID){
scatID = p.ID();
theta_e = acos( pz_boosted /sqrt(px_boosted*px_boosted + py_boosted*py_boosted + pz_boosted*pz_boosted) );
}
......
......@@ -77,7 +77,7 @@ public:
bool found_proton = false;
//Also get the true scattered electron, which will not be included in the sum
//over final-state particles for the JB reconstruction
eic::Index mcscatID;
int32_t mcscatID = -1;
for (const auto& p : mcparts) {
if (p.genStatus() == 4 && p.pdgID() == 11) {
// Incoming electron
......@@ -131,11 +131,11 @@ public:
found_proton = true;
}
// Index of true Scattered electron. Currently taken as first status==1 electron in HEPMC record.
if (p.genStatus() == 1 && p.pdgID() == 11 && !mcscatID) {
if (p.genStatus() == 1 && p.pdgID() == 11 && mcscatID == -1) {
mcscatID = p.ID();
}
if (found_electron && found_proton && mcscatID) {
if (found_electron && found_proton && mcscatID != -1) {
break;
}
}
......@@ -152,13 +152,19 @@ public:
}
return StatusCode::SUCCESS;
}
if (mcscatID == -1) {
if (msgLevel(MSG::DEBUG)) {
debug() << "No truth scattered electron found" << endmsg;
}
return StatusCode::SUCCESS;
}
// Loop over reconstructed particles to get outgoing scattered electron
// Use the true scattered electron from the MC information
typedef std::pair<eic::VectorXYZT, eic::Index> t_electron;
std::vector<t_electron> electrons;
for (const auto& p : parts) {
if (p.mcID() == mcscatID) {
if (p.mcID().value == mcscatID) {
// Outgoing electron
electrons.push_back(t_electron(eic::VectorXYZT(p.p().x, p.p().y, p.p().z, p.energy()), eic::Index(p.ID())));
}
......
......@@ -77,7 +77,7 @@ public:
bool found_proton = false;
//Also get the true scattered electron, which will not be included in the sum
//over final-state particles for the JB reconstruction
eic::Index mcscatID;
int32_t mcscatID = -1;
for (const auto& p : mcparts) {
if (p.genStatus() == 4 && p.pdgID() == 11) {
// Incoming electron
......@@ -131,11 +131,11 @@ public:
found_proton = true;
}
// Index of true Scattered electron. Currently taken as first status==1 electron in HEPMC record.
if (p.genStatus() == 1 && p.pdgID() == 11 && !mcscatID) {
if (p.genStatus() == 1 && p.pdgID() == 11 && mcscatID == -1) {
mcscatID = p.ID();
}
if (found_electron && found_proton && mcscatID) {
if (found_electron && found_proton && mcscatID != -1) {
break;
}
}
......@@ -151,6 +151,12 @@ public:
}
return StatusCode::SUCCESS;
}
if (mcscatID == -1) {
if (msgLevel(MSG::DEBUG)) {
debug() << "No truth scattered electron found" << endmsg;
}
return StatusCode::SUCCESS;
}
// Loop over reconstructed particles to get all outgoing particles other than the scattered electron
// -----------------------------------------------------------------
......@@ -176,7 +182,7 @@ public:
for (const auto& p : parts) {
//Get Reconstructed Index of scattered electron
if (p.mcID() == mcscatID){
if (p.mcID().value == mcscatID){
scatID = p.ID();
}
//Sum over all particles other than scattered electron
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment