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

Last change on this file since 1309 was 807, checked in by garnier, 16 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.