| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 #include "Par04SensitiveDetector.hh" 26 #include "Par04SensitiveDetector.hh"
27 << 27 #include <CLHEP/Vector/Rotation.h> // for HepRotation
28 #include "Par04EventInformation.hh" // for Pa <<
29 #include "Par04Hit.hh" // for Par04Hit, Par04 <<
30 <<
31 #include "G4Event.hh" // for G4Event <<
32 #include "G4EventManager.hh" // for G4EventMa <<
33 #include "G4HCofThisEvent.hh" // for G4HCofTh <<
34 #include "G4SDManager.hh" // for G4SDManager <<
35 #include "G4Step.hh" // for G4Step <<
36 #include "G4Track.hh" // for G4Track <<
37 <<
38 #include <CLHEP/Vector/Rotation.h> // for Hep <<
39 #include <CLHEP/Vector/ThreeVector.h> // for 28 #include <CLHEP/Vector/ThreeVector.h> // for Hep3Vector
40 #include <G4CollectionNameVector.hh> // for G << 29 #include <cmath> // for floor
41 #include <G4FastHit.hh> // for G4FastHit << 30 #include <G4CollectionNameVector.hh> // for G4CollectionNameVector
42 #include <G4FastTrack.hh> // for G4FastTrack << 31 #include <G4FastHit.hh> // for G4FastHit
43 #include <G4RotationMatrix.hh> // for G4Rotat << 32 #include <G4FastTrack.hh> // for G4FastTrack
44 #include <G4StepPoint.hh> // for G4StepPoint << 33 #include <G4RotationMatrix.hh> // for G4RotationMatrix
45 #include <G4THitsCollection.hh> // for G4THit << 34 #include <G4StepPoint.hh> // for G4StepPoint
46 #include <G4ThreeVector.hh> // for G4ThreeVec << 35 #include <G4THitsCollection.hh> // for G4THitsCollection
47 #include <G4VSensitiveDetector.hh> // for G4V << 36 #include <G4ThreeVector.hh> // for G4ThreeVector
>> 37 #include <G4VSensitiveDetector.hh> // for G4VSensitiveDetector
48 #include <G4VUserEventInformation.hh> // for 38 #include <G4VUserEventInformation.hh> // for G4VUserEventInformation
49 #include <cmath> // for floor << 39 #include <cstddef> // for size_t
50 #include <cstddef> // for size_t << 40 #include <vector> // for vector
51 #include <vector> // for vector << 41 #include "G4Event.hh" // for G4Event
>> 42 #include "G4EventManager.hh" // for G4EventManager
>> 43 #include "G4HCofThisEvent.hh" // for G4HCofThisEvent
>> 44 #include "G4SDManager.hh" // for G4SDManager
>> 45 #include "G4Step.hh" // for G4Step
>> 46 #include "G4Track.hh" // for G4Track
>> 47 #include "Par04EventInformation.hh" // for Par04EventInformation
>> 48 #include "Par04Hit.hh" // for Par04Hit, Par04HitsCollection
52 49
53 //....oooOO0OOooo........oooOO0OOooo........oo 50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
54 51
55 Par04SensitiveDetector::Par04SensitiveDetector << 52 Par04SensitiveDetector::Par04SensitiveDetector(G4String aName)
>> 53 : G4VSensitiveDetector(aName)
56 { 54 {
57 collectionName.insert("hits"); 55 collectionName.insert("hits");
58 } 56 }
59 //....oooOO0OOooo........oooOO0OOooo........oo 57 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
60 58
61 Par04SensitiveDetector::Par04SensitiveDetector 59 Par04SensitiveDetector::Par04SensitiveDetector(G4String aName, G4ThreeVector aNb,
62 60 G4ThreeVector aSize)
63 : G4VSensitiveDetector(aName), fMeshNbOfCell << 61 : G4VSensitiveDetector(aName)
>> 62 , fMeshNbOfCells(aNb)
>> 63 , fMeshSizeOfCells(aSize)
64 { 64 {
65 collectionName.insert("hits"); 65 collectionName.insert("hits");
66 } 66 }
67 67
68 //....oooOO0OOooo........oooOO0OOooo........oo 68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
69 69
70 Par04SensitiveDetector::~Par04SensitiveDetecto << 70 Par04SensitiveDetector::~Par04SensitiveDetector() {}
71 71
72 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73 73
74 void Par04SensitiveDetector::Initialize(G4HCof 74 void Par04SensitiveDetector::Initialize(G4HCofThisEvent* aHCE)
75 { 75 {
76 fHitsCollection = new Par04HitsCollection(Se 76 fHitsCollection = new Par04HitsCollection(SensitiveDetectorName, collectionName[0]);
77 if (fHitCollectionID < 0) { << 77 if(fHitCollectionID < 0)
>> 78 {
78 fHitCollectionID = G4SDManager::GetSDMpoin 79 fHitCollectionID = G4SDManager::GetSDMpointer()->GetCollectionID(fHitsCollection);
79 } 80 }
80 aHCE->AddHitsCollection(fHitCollectionID, fH 81 aHCE->AddHitsCollection(fHitCollectionID, fHitsCollection);
81 82
>> 83 // fill calorimeter hits with zero energy deposition
>> 84 for(G4int iphi = 0; iphi < fMeshNbOfCells.y(); iphi++)
>> 85 for(G4int irho = 0; irho < fMeshNbOfCells.x(); irho++)
>> 86 for(G4int iz = 0; iz < fMeshNbOfCells.z(); iz++)
>> 87 {
>> 88 Par04Hit* hit = new Par04Hit();
>> 89 fHitsCollection->insert(hit);
>> 90 }
82 // reset entrance position 91 // reset entrance position
83 fEntrancePosition.set(-1, -1, -1); 92 fEntrancePosition.set(-1, -1, -1);
84 fEntranceDirection.set(-1, -1, -1); 93 fEntranceDirection.set(-1, -1, -1);
85 } 94 }
86 95
87 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
88 97
89 G4bool Par04SensitiveDetector::ProcessHits(G4S 98 G4bool Par04SensitiveDetector::ProcessHits(G4Step* aStep, G4TouchableHistory*)
90 { 99 {
91 G4double edep = aStep->GetTotalEnergyDeposit 100 G4double edep = aStep->GetTotalEnergyDeposit();
92 if (edep == 0.) return true; << 101 if(edep == 0.)
>> 102 return true;
93 103
94 auto hit = RetrieveAndSetupHit(aStep->GetPos 104 auto hit = RetrieveAndSetupHit(aStep->GetPostStepPoint()->GetPosition());
95 if (hit == nullptr) return true; << 105 if(hit == nullptr)
>> 106 return true;
96 107
97 // Add energy deposit from G4Step 108 // Add energy deposit from G4Step
98 hit->AddEdep(edep); 109 hit->AddEdep(edep);
99 // Increment the counter <<
100 hit->AddNdep(1); <<
101 110
102 // Fill time information from G4Step 111 // Fill time information from G4Step
103 // If it's already filled, choose hit with e 112 // If it's already filled, choose hit with earliest global time
104 if (hit->GetTime() == -1 || hit->GetTime() > << 113 if(hit->GetTime() == -1 || hit->GetTime() > aStep->GetTrack()->GetGlobalTime())
105 hit->SetTime(aStep->GetTrack()->GetGlobalT 114 hit->SetTime(aStep->GetTrack()->GetGlobalTime());
106 115
107 // Set hit type to full simulation (only if 116 // Set hit type to full simulation (only if hit is not already marked as fast
108 // sim) 117 // sim)
109 if (hit->GetType() != 1) hit->SetType(0); << 118 if(hit->GetType() != 1)
>> 119 hit->SetType(0);
110 120
111 return true; 121 return true;
112 } 122 }
113 123
114 //....oooOO0OOooo........oooOO0OOooo........oo 124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
115 125
116 G4bool Par04SensitiveDetector::ProcessHits(con 126 G4bool Par04SensitiveDetector::ProcessHits(const G4FastHit* aHit, const G4FastTrack* aTrack,
117 G4T 127 G4TouchableHistory*)
118 { 128 {
119 G4double edep = aHit->GetEnergy(); 129 G4double edep = aHit->GetEnergy();
120 if (edep == 0.) return true; << 130 if(edep == 0.)
>> 131 return true;
121 132
122 auto hit = RetrieveAndSetupHit(aHit->GetPosi 133 auto hit = RetrieveAndSetupHit(aHit->GetPosition());
123 if (hit == nullptr) return true; << 134 if(hit == nullptr)
>> 135 return true;
124 136
125 // Add energy deposit from G4FastHit 137 // Add energy deposit from G4FastHit
126 hit->AddEdep(edep); 138 hit->AddEdep(edep);
127 // Increment the counter <<
128 hit->AddNdep(1); <<
129 139
130 // Fill time information from G4FastTrack 140 // Fill time information from G4FastTrack
131 // If it's already filled, choose hit with e 141 // If it's already filled, choose hit with earliest global time
132 if (hit->GetTime() == -1 || hit->GetTime() > << 142 if(hit->GetTime() == -1 || hit->GetTime() > aTrack->GetPrimaryTrack()->GetGlobalTime())
>> 143 {
133 hit->SetTime(aTrack->GetPrimaryTrack()->Ge 144 hit->SetTime(aTrack->GetPrimaryTrack()->GetGlobalTime());
134 } 145 }
135 146
136 // Set hit type to fast simulation (even if 147 // Set hit type to fast simulation (even if hit was already marked as full
137 // sim, overwrite it) 148 // sim, overwrite it)
138 hit->SetType(1); 149 hit->SetType(1);
139 150
140 return true; 151 return true;
141 } 152 }
142 153
143 //....oooOO0OOooo........oooOO0OOooo........oo 154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
144 155
145 Par04Hit* Par04SensitiveDetector::RetrieveAndS 156 Par04Hit* Par04SensitiveDetector::RetrieveAndSetupHit(G4ThreeVector aGlobalPosition)
146 { 157 {
147 if (fEntrancePosition.x() == -1) { << 158 if(fEntrancePosition.x() == -1)
148 auto info = dynamic_cast<Par04EventInforma << 159 {
>> 160 Par04EventInformation* info = dynamic_cast<Par04EventInformation*>(
149 G4EventManager::GetEventManager()->GetCo 161 G4EventManager::GetEventManager()->GetConstCurrentEvent()->GetUserInformation());
150 if (info == nullptr) return nullptr; << 162 if(info == nullptr)
151 fEntrancePosition = info->GetPosition(); << 163 return nullptr;
>> 164 fEntrancePosition = info->GetPosition();
152 fEntranceDirection = info->GetDirection(); 165 fEntranceDirection = info->GetDirection();
153 } 166 }
154 167
155 auto delta = aGlobalPosition - fEntrancePosi 168 auto delta = aGlobalPosition - fEntrancePosition;
156 169
157 // Calculate rotation matrix along the parti 170 // Calculate rotation matrix along the particle momentum direction
158 // It will rotate the shower axes to match t 171 // It will rotate the shower axes to match the incoming particle direction
159 G4RotationMatrix rotMatrix = G4RotationMatri 172 G4RotationMatrix rotMatrix = G4RotationMatrix();
160 double particleTheta = fEntranceDirection.th << 173 double particleTheta = fEntranceDirection.theta();
161 double particlePhi = fEntranceDirection.phi( << 174 double particlePhi = fEntranceDirection.phi();
162 rotMatrix.rotateZ(-particlePhi); 175 rotMatrix.rotateZ(-particlePhi);
163 rotMatrix.rotateY(-particleTheta); 176 rotMatrix.rotateY(-particleTheta);
164 G4RotationMatrix rotMatrixInv = CLHEP::inver 177 G4RotationMatrix rotMatrixInv = CLHEP::inverseOf(rotMatrix);
165 178
166 delta = rotMatrix * delta; 179 delta = rotMatrix * delta;
167 180
168 G4int rhoNo = std::floor(delta.perp() / fMes 181 G4int rhoNo = std::floor(delta.perp() / fMeshSizeOfCells.x());
169 G4int phiNo = std::floor((CLHEP::pi + delta. 182 G4int phiNo = std::floor((CLHEP::pi + delta.phi()) / fMeshSizeOfCells.y());
170 G4int zNo = std::floor(delta.z() / fMeshSize << 183 G4int zNo = std::floor(delta.z() / fMeshSizeOfCells.z());
171 184
172 std::size_t hitID = 185 std::size_t hitID =
173 fMeshNbOfCells.x() * fMeshNbOfCells.z() * 186 fMeshNbOfCells.x() * fMeshNbOfCells.z() * phiNo + fMeshNbOfCells.z() * rhoNo + zNo;
174 187
175 if (zNo >= fMeshNbOfCells.z() || rhoNo >= fM << 188 if(hitID >= fHitsCollection->entries() || zNo > fMeshNbOfCells.z() ||
>> 189 rhoNo > fMeshNbOfCells.x() || zNo < 0)
>> 190 {
176 return nullptr; 191 return nullptr;
177 } 192 }
178 193
179 auto hit = fHitsMap[hitID].get(); << 194 Par04Hit* hit = (*fHitsCollection)[hitID];
180 if (hit == nullptr) { << 195
181 fHitsMap[hitID] = std::unique_ptr<Par04Hit << 196 if(hit->GetRhoId() < 0)
182 hit = fHitsMap[hitID].get(); << 197 {
183 hit->SetPhiId(phiNo); <<
184 hit->SetRhoId(rhoNo); 198 hit->SetRhoId(rhoNo);
>> 199 hit->SetPhiId(phiNo);
185 hit->SetZid(zNo); 200 hit->SetZid(zNo);
186 hit->SetRot(rotMatrixInv); 201 hit->SetRot(rotMatrixInv);
187 hit->SetPos(fEntrancePosition << 202 hit->SetPos(fEntrancePosition +
188 + rotMatrixInv * G4ThreeVector << 203 rotMatrixInv * G4ThreeVector(0, 0, (zNo + 0.5) * fMeshSizeOfCells.z()));
189 } 204 }
190 return hit; 205 return hit;
191 } <<
192 <<
193 //....oooOO0OOooo........oooOO0OOooo........oo <<
194 <<
195 void Par04SensitiveDetector::EndOfEvent(G4HCof <<
196 { <<
197 for (const auto& hits : fHitsMap) { <<
198 fHitsCollection->insert(new Par04Hit(*hits <<
199 } <<
200 fHitsMap.clear(); <<
201 } 206 }
202 207