[807] | 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 | |
---|