Skip to content

add FILTER_BEAM to cut away events only recording beam particles #205

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 3 commits into
base: main
Choose a base branch
from
Open

add FILTER_BEAM to cut away events only recording beam particles #205

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
26ad7a3
add FILTER_BEAM to cut away events only recording beam particles
Aug 2, 2021
9c20243
Merge branch 'gemc:master' into master
zhaozhiwen Sep 22, 2021
10246e3
change procID from int to string
Sep 22, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion output/gbank.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -119,7 +119,7 @@ map<string, gBank> read_banks(goptions gemcOpt, map<string, string> allSystems)
abank.load_variable("mvz", 22, "Rd", "z component of point of origin of the mother of the particle entering the sensitive volume");
abank.load_variable("avg_t", 23, "Rd", "Average time");
abank.load_variable("nsteps", 24, "Ri", "Number of geant4 steps");
abank.load_variable("procID", 25, "Ri", "Process that created the particle. It's an integer described at gemc.jlab.org");
abank.load_variable("procID", 25, "Rs", "Process that created the particle. It's an integer described at gemc.jlab.org");
abank.load_variable("hitn", 99, "Ri", "Hit Number");
abank.orderNames();
banks["raws"] = abank;
Expand Down
4 changes: 2 additions & 2 deletions sensitivity/Hit.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@ MHit::MHit(double energy, double tim, vector<identifier> vid, int pid)
otrackID.push_back(-1);
mvert.push_back(G4ThreeVector(0,0,0));
materialName.push_back("noise");
processID.push_back(999);
processID.push_back("999");
mgnf.push_back(0);

identity = vid;
Expand Down Expand Up @@ -111,7 +111,7 @@ MHit::MHit(double energy, double tim, int nphe, vector<identifier> vid)
otrackID.push_back(-1);
mvert.push_back(G4ThreeVector(0,0,0));
materialName.push_back("backgroundHit");
processID.push_back(999);
processID.push_back("999");
mgnf.push_back(0);

identity = vid;
Expand Down
10 changes: 5 additions & 5 deletions sensitivity/Hit.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,7 @@ class MHit : public G4VHit
vector<int> otrackID; ///< Original G4Track ID in each step
vector<G4ThreeVector> mvert; ///< Primary Vertex of the track's mother
vector<string> materialName; ///< Material name
vector<int> processID; ///< Process that originated this step
vector<string> processID; ///< Process that originated this step
vector<double> mgnf; ///< magnetic field

vector<detector> Detectors; ///< Detectors Hit. It might be a vector if multiple detectors have the same identifier
Expand Down Expand Up @@ -156,10 +156,10 @@ class MHit : public G4VHit
inline string GetMatName() { return materialName[0]; }
inline vector<string> GetMatNames() { return materialName; }

inline void SetProcID(int procID) { processID.push_back(procID); }
inline void SetProcID(vector<int> procIDs) { processID = procIDs; }
inline int GetProcID() { return processID[0]; }
inline vector<int> GetProcIDs() { return processID; }
inline void SetProcID(string procID) { processID.push_back(procID); }
inline void SetProcID(vector<string> procIDs) { processID = procIDs; }
inline string GetProcID() { return processID[0]; }
inline vector<string> GetProcIDs() { return processID; }

inline void SetSDID(sensitiveID s) { SID = s; }
inline sensitiveID GetSDID() { return SID; }
Expand Down
2 changes: 1 addition & 1 deletion sensitivity/HitProcess.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@ map<string, double> HitProcess::integrateRaw(MHit* aHit, int hitn, bool WRITEBAN
raws["hitn"] = hitn;
raws["totEdep"] = aHit->GetEdep().front();
raws["avg_t"] = aHit->GetTime().front();
raws["procID"] = -1;
raws["procID"] = "-1";
raws["nsteps"] = 1;

if(filterDummyBanks == false) {
Expand Down
113 changes: 58 additions & 55 deletions sensitivity/sensitiveDetector.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -371,62 +371,65 @@ MHit* sensitiveDetector::find_existing_hit(vector<identifier> PID) ///< return
// to check process name go to $G4ROOT/$GEANT4_VERSION/source/geant$GEANT4_VERSION/source/processes/
// mgrep "const G4String&" | grep process
// notice: this map should be written out in the output stream
int sensitiveDetector::processID(string procName)
string sensitiveDetector::processID(string procName)
{
if(procName == "eIoni") return 1;
if(procName == "compt") return 2;
if(procName == "eBrem") return 3;
if(procName == "phot") return 4;
if(procName == "conv") return 5;
if(procName == "annihil") return 6;
if(procName == "photonNuclear") return 7;
if(procName == "electronNuclear") return 8;
if(procName == "positronNuclear") return 9;
if(procName == "CoulombScat") return 10;
if(procName == "Cerenkov") return 11;
if(procName == "Scintillation") return 12;
if(procName == "SynRad") return 13;
if(procName == "SynchrotronRadiation") return 14;
if(procName == "hadElastic") return 20;
if(procName == "hBrems") return 21;
if(procName == "hIoni") return 22;
if(procName == "hPairProd") return 23;
if(procName == "protonInelastic") return 30;
if(procName == "neutronInelastic") return 31;
if(procName == "nCapture") return 32;
if(procName == "anti_neutronInelastic") return 33;
if(procName == "anti_protonInelastic") return 34;
if(procName == "hBertiniCaptureAtRest") return 35;
if(procName == "pi-Inelastic") return 40;
if(procName == "pi+Inelastic") return 41;
if(procName == "Decay") return 50;
if(procName == "DecayWithSpin") return 51;
if(procName == "muIoni") return 60;
if(procName == "muPairProd") return 61;
if(procName == "muBrems") return 62;
if(procName == "muonNuclear") return 63;
if(procName == "muMinusCaptureAtRest") return 64;
if(procName == "kaon-Inelastic") return 70;
if(procName == "kaon+Inelastic") return 71;
if(procName == "kaon0Inelastic") return 72;
if(procName == "kaon0LInelastic") return 73;
if(procName == "kaon0SInelastic") return 74;
if(procName == "alphaInelastic") return 80;
if(procName == "lambdaInelastic") return 90;
if(procName == "sigma-Inelastic") return 100;
if(procName == "dInelastic") return 110;
if(procName == "ionIoni") return 120;
if(procName == "ionInelastic") return 121;
if(procName == "tInelastic") return 130;
if(procName == "xi0Inelastic") return 140;
if(procName == "omega-Inelastic") return 150;

if(procName == "na") return 999;

if(verbosity > 0)
cout << " Warning: process name " << procName << " not catalogued." << endl;

return 999;
// cout << " Warning: process name " << procName << " alllll." << endl;
// if(procName == "eIoni") return 1;
// if(procName == "compt") return 2;
// if(procName == "eBrem") return 3;
// if(procName == "phot") return 4;
// if(procName == "conv") return 5;
// if(procName == "annihil") return 6;
// if(procName == "photonNuclear") return 7;
// if(procName == "electronNuclear") return 8;
// if(procName == "positronNuclear") return 9;
// if(procName == "CoulombScat") return 10;
// if(procName == "Cerenkov") return 11;
// if(procName == "Scintillation") return 12;
// if(procName == "SynRad") return 13;
// if(procName == "SynchrotronRadiation") return 14;
// if(procName == "hadElastic") return 20;
// if(procName == "hBrems") return 21;
// if(procName == "hIoni") return 22;
// if(procName == "hPairProd") return 23;
// if(procName == "protonInelastic") return 30;
// if(procName == "neutronInelastic") return 31;
// if(procName == "nCapture") return 32;
// if(procName == "anti_neutronInelastic") return 33;
// if(procName == "anti_protonInelastic") return 34;
// if(procName == "hBertiniCaptureAtRest") return 35;
// if(procName == "pi-Inelastic") return 40;
// if(procName == "pi+Inelastic") return 41;
// if(procName == "Decay") return 50;
// if(procName == "DecayWithSpin") return 51;
// if(procName == "muIoni") return 60;
// if(procName == "muPairProd") return 61;
// if(procName == "muBrems") return 62;
// if(procName == "muonNuclear") return 63;
// if(procName == "muMinusCaptureAtRest") return 64;
// if(procName == "kaon-Inelastic") return 70;
// if(procName == "kaon+Inelastic") return 71;
// if(procName == "kaon0Inelastic") return 72;
// if(procName == "kaon0LInelastic") return 73;
// if(procName == "kaon0SInelastic") return 74;
// if(procName == "alphaInelastic") return 80;
// if(procName == "lambdaInelastic") return 90;
// if(procName == "sigma-Inelastic") return 100;
// if(procName == "dInelastic") return 110;
// if(procName == "ionIoni") return 120;
// if(procName == "ionInelastic") return 121;
// if(procName == "tInelastic") return 130;
// if(procName == "xi0Inelastic") return 140;
// if(procName == "omega-Inelastic") return 150;
//
// // if(procName == "na") return 999;
// cout << " Warning: process name " << procName << " not catalogued." << endl;
// if(verbosity > 0)
// cout << " Warning: process name " << procName << " not catalogued." << endl;
//
// return 999;

return procName;
}


2 changes: 1 addition & 1 deletion sensitivity/sensitiveDetector.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -77,7 +77,7 @@ class sensitiveDetector : public G4VSensitiveDetector
MHitCollection* GetMHitCollection() {if(hitCollection) return hitCollection; else return NULL;} ///< returns hit collection
MHit* find_existing_hit(vector<identifier>); ///< returns hit collection hit inside identifer

int processID(string procName); // return an ID from a process name.
string processID(string procName); // return an ID from a process name.
};


Expand Down
64 changes: 63 additions & 1 deletion src/MEventAction.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,8 +92,10 @@ MEventAction::MEventAction(goptions opts, map<string, double> gpars)
gPars = gpars;
MAXP = (int) gemcOpt.optMap["NGENP"].arg;
FILTER_HITS = (int) gemcOpt.optMap["FILTER_HITS"].arg;
FILTER_BEAM = (int) gemcOpt.optMap["FILTER_BEAM"].arg;
FILTER_HADRONS = (int) gemcOpt.optMap["FILTER_HADRONS"].arg;
FILTER_HIGHMOM = (int) gemcOpt.optMap["FILTER_HIGHMOM"].arg;
FILTER_HIGHMOM = (int) gemcOpt.optMap["FILTER_HIGHMOM"].arg;
FILTER_HIGHTHETA = (int) gemcOpt.optMap["FILTER_HIGHTHETA"].arg;
SKIPREJECTEDHITS = (int) gemcOpt.optMap["SKIPREJECTEDHITS"].arg;
rw = runWeights(opts);

Expand Down Expand Up @@ -260,6 +262,38 @@ void MEventAction::EndOfEventAction(const G4Event* evt)
if(anyHit==0) return;
}

if(FILTER_BEAM) {
int foundBeam = 0;
// cout << SeDe_Map.size() << " SeDe_Map.size() " << endl; // always 1
for(map<string, sensitiveDetector*>::iterator it = SeDe_Map.begin(); it!= SeDe_Map.end(); it++) {
MHC = it->second->GetMHitCollection();
if (MHC) nhits = MHC->GetSize();
else nhits = 0;
if (nhits==1){
for (int h=0; h<nhits; h++)
{
vector<int> tids = (*MHC)[h]->GetTIds();
vector<G4ThreeVector> mmts = (*MHC)[h]->GetMoms();
int thiscounter=0;
for (vector<int>::const_iterator tit = tids.begin(); tit != tids.end(); tit++) {

// cout << nhits << " " << tids.size() << " " << mmts.size() << " " << *tit << " " << mmts[thiscounter].z() << endl;
if (*tit ==1 && mmts[thiscounter].z()==FILTER_BEAM) {
foundBeam = 1;
break;
}

thiscounter++;
}

}
}

}

// stop here if there are no hits and FILTER_HITS is set
if(foundBeam) return;
}

// if FILTER_HADRONS is set, checking if there are any (matching) hadrons
if (FILTER_HADRONS == 1 || abs(FILTER_HADRONS) > 99) {
Expand Down Expand Up @@ -315,6 +349,34 @@ void MEventAction::EndOfEventAction(const G4Event* evt)
if (not foundHighmom) return;
}

// if FILTER_HIGHTHETA is set, checking if there is any high mom hits
if (FILTER_HIGHTHETA != 0)
{
int foundHightheta = 0;
for (map<string, sensitiveDetector*>::iterator it = SeDe_Map.begin(); it!= SeDe_Map.end(); it++)
{
MHC = it->second->GetMHitCollection();
if (MHC) nhits = MHC->GetSize();
else nhits = 0;
for (int h=0; h<nhits; h++)
{
vector<G4ThreeVector> mmts = (*MHC)[h]->GetMoms();
for(unsigned int t=0; t<mmts.size(); t++){
cout << "theta " << mmts[t].theta() << " mom " << mmts[t].mag() << endl;
if (mmts[t].theta()/3.1416*180. > FILTER_HIGHTHETA)
{
foundHightheta = 1;
break;
}
}

}
}

// stop here if there are no any high mom hits
if (not foundHightheta) return;
}

if(evtN%Modulo == 0 )
cout << hd_msg << " Starting Event Action Routine " << evtN << " Run Number: " << rw.runNo << endl;

Expand Down
4 changes: 3 additions & 1 deletion src/MEventAction.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -153,8 +153,10 @@ class MEventAction : public G4UserEventAction
int SAVE_ALL_ANCESTORS; ///< Outputs info on all ancestors of tracks with hits
int MAXP; ///< Max number of generated particles to save on output stream
int FILTER_HITS; ///< If set to 1, do not write any output unless there is a hit somewhere
int FILTER_HADRONS; ///< If set to 1, do not write any output unless there is a hadron somewhere
int FILTER_BEAM; ///< If set to non-0, do not write any output when there is only a hit from beam
int FILTER_HADRONS; ///< If set to non-0, do not write any output unless there is a hadron somewhere
int FILTER_HIGHMOM; ///< If set to non-0, do not write any output unless there is high mom hit
int FILTER_HIGHTHETA; ///< If set to non-0, do not write any output unless there is high theta angle hit
int SKIPREJECTEDHITS; ///< Skips hits that are rejected by digitization. Default: yes
string WRITE_ALLRAW; ///< List of detectors for which geant4 all raw info need to be saved
string WRITE_INTRAW; ///< List of detectors for which geant4 raw integrated info need to be saved
Expand Down
15 changes: 15 additions & 0 deletions src/gemc_options.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -776,6 +776,12 @@ void goptions::setGoptions()
optMap["FILTER_HITS"].type = 0;
optMap["FILTER_HITS"].ctgr = "output";

optMap["FILTER_BEAM"].arg = 0;
optMap["FILTER_BEAM"].help = "If set to non-0, but as beam Mom in MeV, do not write output if there are only 1 hit from the beam";
optMap["FILTER_BEAM"].name = "If set to non-0, but as beam Mom in MeV, do not write output if there are only 1 hit from the beam";
optMap["FILTER_BEAM"].type = 0;
optMap["FILTER_BEAM"].ctgr = "output";

optMap["FILTER_HADRONS"].arg = 0;
optMap["FILTER_HADRONS"].help = "If set to 1, do not write events if there are no hadrons. Otherwise if \n";
optMap["FILTER_HADRONS"].help += "nonzero write only events having a hadron with matching ID. For example\n";
Expand All @@ -793,6 +799,15 @@ void goptions::setGoptions()
optMap["FILTER_HIGHMOM"].type = 0;
optMap["FILTER_HIGHMOM"].ctgr = "output";

// No output if no high theta hit
optMap["FILTER_HIGHTHETA"].arg = 0;
optMap["FILTER_HIGHTHETA"].help = "If set to non-0, do not write events if there are no high theta hit. Otherwise if \n";
optMap["FILTER_HIGHTHETA"].help += "nonzero write only events having a hit with theta > FILTER_HIGHTHETA. For example\n";
optMap["FILTER_HIGHTHETA"].help += " -FILTER_HIGHTHETA=1 for theta > 1deg";
optMap["FILTER_HIGHTHETA"].name = "If set to non-0 (or >1), do not write output if there are no high theta hit";
optMap["FILTER_HIGHTHETA"].type = 0;
optMap["FILTER_HIGHTHETA"].ctgr = "output";

// sampling time of electronics (typically FADC), and number of sampling / event
// the VT output is sampled every TSAMPLING nanoseconds to produce a ADC
// the default number of samples is 500 ADC points, at 4ns intervals (total electronic event time = 2 microseconds)
Expand Down