94{
96 vector<BesTruthTrack*>* trackList = sensitiveManager->
GetTrackList();
97 G4int nTrack = trackList->size();
99
100 vector<BesTruthVertex*>* vertexList = sensitiveManager->
GetVertexList();
101 G4int nVertex = vertexList->size();
103
104
105 for(int i=0;i<nTrack-1;i++)
106 for(int j=i+1;j<nTrack;j++)
107 if((*trackList)[i]->GetIndex()>(*trackList)[j]->GetIndex())
108 {
109 track=(*trackList)[i];
110 (*trackList)[i]=(*trackList)[j];
111 (*trackList)[j]=track;
112 }
113
115
116
117 for(int i=0;i<nTrack;i++)
118 {
119 track = (*trackList) [i];
120
121
122 bool isPrimary = false;
123 bool startPointFound = false;
125
126 for (int i=0;i<nVertex;i++)
127 {
128 vertex = (*vertexList) [i];
130 {
132 startPointFound = true;
133 startPoint = vertex;
134 break;
135 }
136 }
137
138 if (!startPointFound)
139 std::cout << "error in finding the start point of a track" <<std::endl;
140
141
142 bool endPointFound = false;
143
144 for (int i=0;i<nVertex;i++)
145 {
146 vertex = (*vertexList) [i];
148 {
150 {
151
153
154
155 HepLorentzVector initialMomentum(track->
GetP4().x()/1000., track->
GetP4().y()/1000., track->
GetP4().z()/1000., track->
GetP4().t()/1000.);
156
158
160
161
163
164
166
169
170
171
172
177 else
179
180
181
182
183
186
187
189 mcParticle->
finalize(finalPosition);
190
191
192 particleCol->push_back(mcParticle);
193
194 endPointFound = true;
195 }
196 }
197 }
198
199 if (!endPointFound)
200 {
201
203 {
204
206
207 HepLorentzVector initialMomentum( track->
GetP4().x()/1000., track->
GetP4().y()/1000., track->
GetP4().z()/1000., track->
GetP4().t()/1000.);
209
211
212
214
215
218
219
222
223
224
225
230 else
232
233
234
235
236
237
238 particleCol->push_back(mcParticle);
239 continue;
240 }
241 }
242 }
243
244
245 SmartRefVector<Event::McParticle> primaryParticleCol;
246 Event::McParticleCol::iterator iter_particle;
247 for ( iter_particle = particleCol->begin();
248 iter_particle != particleCol->end(); iter_particle++) {
249
250 if ((*iter_particle)->primaryParticle()) {
252 primaryParticleCol.push_back(mcPart);
253 }
254 }
255
256 if (primaryParticleCol.empty())
257 std::cout << "no primary particle found!" << std::endl;
258
259
260 SmartRefVector<Event::McParticle>::iterator iter_primary;
261 for ( iter_primary = primaryParticleCol.begin();
262 iter_primary != primaryParticleCol.end(); iter_primary++) {
263 if ( !((*iter_primary)->leafParticle()) )
265 }
266
267
268 StatusCode sc = m_evtSvc->registerObject("/Event/MC/McParticleCol", particleCol);
269 if(sc!=StatusCode::SUCCESS)
270 G4cout<< "Could not register McParticle collection" <<G4endl;
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295}
std::vector< BesTruthVertex * > * GetVertexList()
std::vector< BesTruthTrack * > * GetTrackList()
static BesSensitiveManager * GetSensitiveManager()
HepLorentzVector GetP4() const
BesTruthVertex * GetTerminalVertex() const
BesTruthVertex * GetVertex() const
vector< int > GetDaughterIndexes() const
BesTruthTrack * GetParentTrack() const
G4ThreeVector GetPosition() const
void setVertexIndex0(int index0)
methods for setting and getting vertex indexes
@ DECAYFLT
Decayed by generator.
@ PRIMARY
Decayed in flight.
@ ERROR
this particle is a leaf in the particle tree
void initialize(StdHepId id, unsigned int statusBits, const HepLorentzVector &initialMomentum, const HepLorentzVector &initialPosition, const std::string process="")
Set the initial attributes of the McParticle.
void addStatusFlag(unsigned int flag)
add a new flag to the status flags
void finalize(const HepLorentzVector &finalPosition)
Set the final attributes of the McParticle.
void setVertexIndex1(int index1)
void setTrackIndex(int trackIndex)
ObjectList< McParticle > McParticleCol