source: trunk/examples/novice/gemc/src/MSensitiveDetector.cc@ 1199

Last change on this file since 1199 was 807, checked in by garnier, 17 years ago

update
File size: 9.1 KB
Line 
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
16MSensitiveDetector::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
25MSensitiveDetector::~MSensitiveDetector(){;}
26
27
28void 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
39G4bool 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
182void 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
215int 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
243MHit* 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
Note: See TracBrowser for help on using the repository browser.