diff --git a/output/gbank.cc b/output/gbank.cc index 006eb579..96164e6a 100644 --- a/output/gbank.cc +++ b/output/gbank.cc @@ -119,7 +119,7 @@ map read_banks(goptions gemcOpt, map 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; diff --git a/sensitivity/Hit.cc b/sensitivity/Hit.cc index 655921a9..f5a3c51a 100644 --- a/sensitivity/Hit.cc +++ b/sensitivity/Hit.cc @@ -82,7 +82,7 @@ MHit::MHit(double energy, double tim, vector 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; @@ -111,7 +111,7 @@ MHit::MHit(double energy, double tim, int nphe, vector 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; diff --git a/sensitivity/Hit.h b/sensitivity/Hit.h index 27a98ea8..b86eb019 100644 --- a/sensitivity/Hit.h +++ b/sensitivity/Hit.h @@ -61,7 +61,7 @@ class MHit : public G4VHit vector otrackID; ///< Original G4Track ID in each step vector mvert; ///< Primary Vertex of the track's mother vector materialName; ///< Material name - vector processID; ///< Process that originated this step + vector processID; ///< Process that originated this step vector mgnf; ///< magnetic field vector Detectors; ///< Detectors Hit. It might be a vector if multiple detectors have the same identifier @@ -156,10 +156,10 @@ class MHit : public G4VHit inline string GetMatName() { return materialName[0]; } inline vector GetMatNames() { return materialName; } - inline void SetProcID(int procID) { processID.push_back(procID); } - inline void SetProcID(vector procIDs) { processID = procIDs; } - inline int GetProcID() { return processID[0]; } - inline vector GetProcIDs() { return processID; } + inline void SetProcID(string procID) { processID.push_back(procID); } + inline void SetProcID(vector procIDs) { processID = procIDs; } + inline string GetProcID() { return processID[0]; } + inline vector GetProcIDs() { return processID; } inline void SetSDID(sensitiveID s) { SID = s; } inline sensitiveID GetSDID() { return SID; } diff --git a/sensitivity/HitProcess.cc b/sensitivity/HitProcess.cc index 2f0e65be..f47ac66e 100644 --- a/sensitivity/HitProcess.cc +++ b/sensitivity/HitProcess.cc @@ -42,7 +42,7 @@ map 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) { diff --git a/sensitivity/sensitiveDetector.cc b/sensitivity/sensitiveDetector.cc index 7bfd9afd..c8a5fb9a 100644 --- a/sensitivity/sensitiveDetector.cc +++ b/sensitivity/sensitiveDetector.cc @@ -371,62 +371,65 @@ MHit* sensitiveDetector::find_existing_hit(vector 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; } diff --git a/sensitivity/sensitiveDetector.h b/sensitivity/sensitiveDetector.h index 2125d9f7..a904c4e9 100644 --- a/sensitivity/sensitiveDetector.h +++ b/sensitivity/sensitiveDetector.h @@ -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); ///< 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. }; diff --git a/src/MEventAction.cc b/src/MEventAction.cc index 22bdb442..fb9fd7c7 100644 --- a/src/MEventAction.cc +++ b/src/MEventAction.cc @@ -92,8 +92,10 @@ MEventAction::MEventAction(goptions opts, map 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); @@ -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::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 tids = (*MHC)[h]->GetTIds(); + vector mmts = (*MHC)[h]->GetMoms(); + int thiscounter=0; + for (vector::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) { @@ -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::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 mmts = (*MHC)[h]->GetMoms(); + for(unsigned int t=0; t 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; diff --git a/src/MEventAction.h b/src/MEventAction.h index 5033a0ed..8f21861f 100644 --- a/src/MEventAction.h +++ b/src/MEventAction.h @@ -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 diff --git a/src/gemc_options.cc b/src/gemc_options.cc index 87938c1b..9e811df7 100644 --- a/src/gemc_options.cc +++ b/src/gemc_options.cc @@ -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"; @@ -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)