1 | // %%%%%%%%%% |
---|
2 | // G4 headers |
---|
3 | // %%%%%%%%%% |
---|
4 | #include "G4SDManager.hh" |
---|
5 | #include "G4UnitsTable.hh" |
---|
6 | #include "G4ParticleTypes.hh" |
---|
7 | |
---|
8 | |
---|
9 | // %%%%%%%%%%%%% |
---|
10 | // gemc headers |
---|
11 | // %%%%%%%%%%%%% |
---|
12 | #include "identifier.h" |
---|
13 | #include "MSensitiveDetector.h" |
---|
14 | |
---|
15 | |
---|
16 | MSensitiveDetector::MSensitiveDetector(G4String name, gemc_opts opt):G4VSensitiveDetector(name), HCID(-1), gemcOpt(opt) |
---|
17 | { |
---|
18 | HCname = name; |
---|
19 | collectionName.insert(HCname); |
---|
20 | SDID = get_SDId(HCname, gemcOpt); |
---|
21 | SDID.IDshifts = set_IDshifts(SDID.IDmaxs); |
---|
22 | minEnergy = SDID.minEnergy; |
---|
23 | } |
---|
24 | |
---|
25 | MSensitiveDetector::~MSensitiveDetector(){;} |
---|
26 | |
---|
27 | |
---|
28 | void MSensitiveDetector::Initialize(G4HCofThisEvent* HCE) |
---|
29 | { |
---|
30 | Id_Set.clear(); |
---|
31 | Tr_Map.clear(); |
---|
32 | hitCollection = new MHitCollection(HCname, collectionName[0]); |
---|
33 | if(HCID < 0) HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); |
---|
34 | HCE->AddHitsCollection( HCID, hitCollection ); |
---|
35 | } |
---|
36 | |
---|
37 | |
---|
38 | |
---|
39 | G4bool MSensitiveDetector::ProcessHits(G4Step* aStep, G4TouchableHistory*) |
---|
40 | { |
---|
41 | // First check on energy deposited |
---|
42 | double depe = aStep->GetTotalEnergyDeposit(); |
---|
43 | /* if(depe < minEnergy && aStep->GetTrack()->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) return false;*/ |
---|
44 | if(depe == 0 && aStep->GetTrack()->GetDefinition() != G4OpticalPhoton::OpticalPhotonDefinition()) return false; |
---|
45 | |
---|
46 | string hd_msg1 = gemcOpt.args["LOG_MSG"].args + " New Hit: <<< "; |
---|
47 | string hd_msg2 = gemcOpt.args["LOG_MSG"].args + " > "; |
---|
48 | string catch_v = gemcOpt.args["CATCH"].args; |
---|
49 | double HIT_VERBOSITY = gemcOpt.args["HIT_VERBOSITY"].arg ; |
---|
50 | |
---|
51 | G4StepPoint *prestep = aStep->GetPreStepPoint(); |
---|
52 | G4StepPoint *poststep = aStep->GetPostStepPoint(); |
---|
53 | G4VTouchable* TH = (G4VTouchable*) aStep->GetPreStepPoint()->GetTouchable(); |
---|
54 | |
---|
55 | string name = TH->GetVolume(0)->GetName(); |
---|
56 | G4ThreeVector xyz = poststep->GetPosition(); |
---|
57 | G4ThreeVector Lxyz = prestep->GetTouchableHandle()->GetHistory() |
---|
58 | ->GetTopTransform().TransformPoint(xyz); |
---|
59 | double ctime = poststep->GetGlobalTime(); |
---|
60 | G4ThreeVector pxyz = prestep->GetMomentum(); |
---|
61 | double ene = prestep->GetTotalEnergy(); |
---|
62 | G4ThreeVector ver = aStep->GetTrack()->GetVertexPosition(); |
---|
63 | string pid = aStep->GetTrack()->GetDefinition()->GetParticleName(); |
---|
64 | vector<identifier> VID = SetId(GetDetectorIdentifier(name), TH, ctime, SDID.TimeWindow); |
---|
65 | int tid = aStep->GetTrack()->GetTrackID(); |
---|
66 | int ptid = aStep->GetTrack()->GetParentID(); |
---|
67 | |
---|
68 | if(Tr_Map.find(tid) == Tr_Map.end()) |
---|
69 | { |
---|
70 | MTrackInfo trinf; |
---|
71 | trinf.pid = pid; |
---|
72 | trinf.v = ver; |
---|
73 | Tr_Map[tid] = trinf; |
---|
74 | } |
---|
75 | |
---|
76 | MPHBaseClass *ProcessHitRoutine = GetMPHClass((*MProcessHit_Map), GetDetectorHitType(name)); |
---|
77 | if(!ProcessHitRoutine) |
---|
78 | { |
---|
79 | cout << endl << hd_msg2 << " Hit Process Routine <" << GetDetectorHitType(name) << "> not found. Exiting. " << endl; |
---|
80 | exit(0); |
---|
81 | } |
---|
82 | |
---|
83 | vector<identifier> PID = ProcessHitRoutine->ProcessID(VID, aStep, (*Hall_Map)[name]); |
---|
84 | |
---|
85 | if(HIT_VERBOSITY > 11 || name.find(catch_v) != string::npos) |
---|
86 | cout << endl << hd_msg2 << " Before hit Process Identification:" << endl << VID |
---|
87 | << hd_msg2 << " After: hit Process Identification:" << endl << PID << endl; |
---|
88 | |
---|
89 | |
---|
90 | set<vector<identifier> > :: iterator itid; |
---|
91 | int hit_found = 0; |
---|
92 | |
---|
93 | if(HIT_VERBOSITY > 5 || name.find(catch_v) != string::npos) cout << endl << endl << " BEGIN SEARCH for same hit in Identifier Set..." << endl; |
---|
94 | |
---|
95 | for(itid = Id_Set.begin(); itid!= Id_Set.end() && !hit_found; itid++) |
---|
96 | { |
---|
97 | if(*itid == PID) hit_found=1; |
---|
98 | if(HIT_VERBOSITY > 5 || name.find(catch_v) != string::npos) |
---|
99 | cout << " >> Current Step: " << PID |
---|
100 | << " >> Set Hit Index: " << *itid |
---|
101 | << (hit_found ? " >> FOUND at this Set Entry. " : " >> Not found yet. ") << endl; |
---|
102 | } |
---|
103 | if(HIT_VERBOSITY > 5 || name.find(catch_v) != string::npos) cout << " SEARCH ENDED." << (hit_found ? " 1 " : " No ") << "hit found in the Set." << endl << endl; |
---|
104 | |
---|
105 | // New Hit |
---|
106 | if(!hit_found) |
---|
107 | { |
---|
108 | MHit *thisHit = new MHit(); |
---|
109 | thisHit->SetPos(xyz); |
---|
110 | thisHit->SetLPos(Lxyz); |
---|
111 | thisHit->SetEdep(depe); |
---|
112 | thisHit->SetTime(ctime); |
---|
113 | thisHit->SetMom(pxyz); |
---|
114 | thisHit->SetE(ene); |
---|
115 | thisHit->SetVert(ver); |
---|
116 | thisHit->SetPID(pid); |
---|
117 | thisHit->SetId(PID); |
---|
118 | thisHit->SetTrackId(tid); |
---|
119 | thisHit->SetDetector((*Hall_Map)[name]); |
---|
120 | thisHit->SetThreshold(SDID.minEnergy); |
---|
121 | if(Tr_Map.count(ptid) > 0) |
---|
122 | { |
---|
123 | thisHit->SetmPID(Tr_Map[ptid].pid); |
---|
124 | thisHit->SetmVert(Tr_Map[ptid].v); |
---|
125 | } |
---|
126 | else |
---|
127 | { |
---|
128 | thisHit->SetmPID("unknown"); |
---|
129 | thisHit->SetmVert(G4ThreeVector(0, 0, 0)); |
---|
130 | } |
---|
131 | |
---|
132 | hitCollection->insert(thisHit); |
---|
133 | Id_Set.insert(PID); |
---|
134 | |
---|
135 | if(HIT_VERBOSITY > 3 || name.find(catch_v) != string::npos) |
---|
136 | cout << endl << hd_msg1 << endl |
---|
137 | << " > This element was not hit yet in this event. Identity:" << endl << thisHit->GetId() |
---|
138 | << " > Creating new hit by a E=" << ene/MeV << ", p=" << pxyz.mag()/MeV << " MeV " << pid |
---|
139 | << ", track ID = " << tid << ", inside Hit Collection <" << HCname << ">." << endl |
---|
140 | << " > Energy Deposited this step: " << depe/MeV << " MeV" << endl |
---|
141 | << " > Time of this step: " << G4BestUnit(ctime, "Time") << endl |
---|
142 | << " > Position of this step: " << xyz/cm << " cm" << endl |
---|
143 | << " > Local Position in volume: " << Lxyz/cm << " cm" << endl; |
---|
144 | } |
---|
145 | else |
---|
146 | { |
---|
147 | // Adding hit info only if the poststeppint remains in the volume? |
---|
148 | //if( aStep->GetPreStepPoint()->GetTouchable()->GetVolume(0) == aStep->GetPostStepPoint()->GetTouchable()->GetVolume(0)) |
---|
149 | { |
---|
150 | MHit *thisHit = find_existing_hit(PID); |
---|
151 | if(!thisHit) |
---|
152 | { |
---|
153 | cout << " Hit not found in collection but found in PID. This should never happen. Exiting." << endl; |
---|
154 | exit(0); |
---|
155 | } |
---|
156 | else |
---|
157 | { |
---|
158 | thisHit->SetPos(xyz); |
---|
159 | thisHit->SetLPos(Lxyz); |
---|
160 | thisHit->SetEdep(depe); |
---|
161 | thisHit->SetTime(ctime); |
---|
162 | thisHit->SetTrackId(tid); |
---|
163 | |
---|
164 | if(HIT_VERBOSITY > 3 || name.find(catch_v) != string::npos) |
---|
165 | cout << hd_msg2 << " Step Number " << thisHit->GetPos().size() |
---|
166 | << " inside Identity: " << endl << thisHit->GetId() |
---|
167 | << " > Adding hit inside Hit Collection <" << HCname << ">." |
---|
168 | << " by a E=" << ene/MeV << ", p=" << pxyz.mag()/MeV << " MeV " |
---|
169 | << pid << ", track ID = " << tid << endl |
---|
170 | << " > Energy Deposited this step: " << depe/MeV << " MeV" << endl |
---|
171 | << " > Time of this step: " << G4BestUnit(ctime, "Time") |
---|
172 | << " is within this element Time Window of " << SDID.TimeWindow/ns << " ns. " << endl |
---|
173 | << " > Position of this step: " << xyz/cm << " cm" << endl |
---|
174 | << " > Local Position in volume: " << Lxyz/cm << " cm" << endl; |
---|
175 | } |
---|
176 | } |
---|
177 | } |
---|
178 | return true; |
---|
179 | } |
---|
180 | |
---|
181 | |
---|
182 | void MSensitiveDetector::EndOfEvent(G4HCofThisEvent *HCE) |
---|
183 | { |
---|
184 | string hd_msg = gemcOpt.args["LOG_MSG"].args + " End of Hit Collections >> "; |
---|
185 | double HIT_VERBOSITY = gemcOpt.args["HIT_VERBOSITY"].arg ; |
---|
186 | string catch_v = gemcOpt.args["CATCH"].args; |
---|
187 | |
---|
188 | int nhitC = hitCollection->GetSize(); |
---|
189 | if(!nhitC) return; |
---|
190 | MHit *aHit; |
---|
191 | double Etot; |
---|
192 | if(HIT_VERBOSITY > 2 && nhitC) |
---|
193 | { |
---|
194 | cout << endl; |
---|
195 | cout << hd_msg << " Hit Collections <" << HCname << ">: " << nhitC << " hits." << endl; |
---|
196 | |
---|
197 | for (int i=0; i<nhitC; i++) |
---|
198 | { |
---|
199 | aHit = (*hitCollection)[i]; |
---|
200 | string vname = aHit->GetId()[aHit->GetId().size()-1].name; |
---|
201 | if(HIT_VERBOSITY > 4 || vname.find(catch_v) != string::npos) |
---|
202 | { |
---|
203 | cout << hd_msg << " Hit " << i + 1 << " -- total number of steps this hit: " << aHit->GetPos().size() << endl; |
---|
204 | cout << aHit->GetId(); |
---|
205 | Etot = 0; |
---|
206 | for(int e=0; e<aHit->GetPos().size(); e++) Etot = Etot + aHit->GetEdep()[e]; |
---|
207 | cout << " Total energy deposited: " << Etot/MeV << " MeV" << endl; |
---|
208 | } |
---|
209 | } |
---|
210 | } |
---|
211 | |
---|
212 | |
---|
213 | } |
---|
214 | |
---|
215 | int MSensitiveDetector::get_bit_compressed_id(vector<identifier> id) |
---|
216 | { |
---|
217 | string hd_msg = gemcOpt.args["LOG_MSG"].args + " Bit Compressing >> "; |
---|
218 | double HIT_VERBOSITY = gemcOpt.args["HIT_VERBOSITY"].arg ; |
---|
219 | |
---|
220 | int bitc = 0; |
---|
221 | int depth = id.size(); |
---|
222 | // check for consistency |
---|
223 | if(depth != SDID.IDnames.size()) |
---|
224 | { |
---|
225 | cout << hd_msg << " Something is wrong. Identifier size is different than SDId sizes: " << endl; |
---|
226 | cout << " id size: " << id.size() << " SDID size: " << endl; |
---|
227 | exit(0); |
---|
228 | } |
---|
229 | else |
---|
230 | { |
---|
231 | for(int i=0; i<depth; i++) |
---|
232 | { |
---|
233 | // cout << SDID.IDnames[i] << " " << id[i].name << " " << id[i].id << endl; |
---|
234 | if(i==0) bitc = id[i].id; |
---|
235 | else bitc = bitc + (id[i].id << SDID.IDshifts[i]) ; |
---|
236 | } |
---|
237 | } |
---|
238 | |
---|
239 | return bitc; |
---|
240 | } |
---|
241 | |
---|
242 | |
---|
243 | MHit* MSensitiveDetector::find_existing_hit(vector<identifier> PID) |
---|
244 | { |
---|
245 | for(int i=0; i<hitCollection->GetSize(); i++) |
---|
246 | if((*hitCollection)[i]->GetId() == PID) return (*hitCollection)[i]; |
---|
247 | |
---|
248 | return NULL; |
---|
249 | } |
---|
250 | |
---|
251 | |
---|
252 | |
---|
253 | |
---|
254 | |
---|
255 | |
---|
256 | |
---|
257 | |
---|
258 | |
---|