88{
90 vector<BesTruthTrack*>* trackList = sensitiveManager->
GetTrackList();
91 G4int nTrack = trackList->size();
93
94 vector<BesTruthVertex*>* vertexList = sensitiveManager->
GetVertexList();
95 G4int nVertex = vertexList->size();
97
98
99 for(int i=0;i<nTrack-1;i++)
100 for(int j=i+1;j<nTrack;j++)
101 if((*trackList)[i]->GetIndex()>(*trackList)[j]->GetIndex())
102 {
103 track=(*trackList)[i];
104 (*trackList)[i]=(*trackList)[j];
105 (*trackList)[j]=track;
106 }
107
109
110
111 for(int i=0;i<nTrack;i++)
112 {
113 track = (*trackList) [i];
114
115
116 bool isPrimary = false;
117 bool startPointFound = false;
119
120 for (int i=0;i<nVertex;i++)
121 {
122 vertex = (*vertexList) [i];
124 {
126 startPointFound = true;
127 startPoint = vertex;
128 break;
129 }
130 }
131
132 if (!startPointFound)
133 std::cout << "error in finding the start point of a track" <<std::endl;
134
135
136 bool endPointFound = false;
137
138 for (int i=0;i<nVertex;i++)
139 {
140 vertex = (*vertexList) [i];
142 {
144 {
145
147
148
149 HepLorentzVector initialMomentum(track->
GetP4().x()/1000., track->
GetP4().y()/1000., track->
GetP4().z()/1000., track->
GetP4().t()/1000.);
150
152
154
155
157
158
160
163
164
165
166
171 else
173
174
175
176
177
180
181
183 mcParticle->
finalize(finalPosition);
184
185
186 particleCol->push_back(mcParticle);
187
188 endPointFound = true;
189 }
190 }
191 }
192
193 if (!endPointFound)
194 {
195
197 {
198
200
201 HepLorentzVector initialMomentum( track->
GetP4().x()/1000., track->
GetP4().y()/1000., track->
GetP4().z()/1000., track->
GetP4().t()/1000.);
203
205
206
208
209
212
213
216
217
218
219
224 else
226
227
228
229
230
231
232 particleCol->push_back(mcParticle);
233 continue;
234 }
235 }
236 }
237
238
239 SmartRefVector<Event::McParticle> primaryParticleCol;
240 Event::McParticleCol::iterator iter_particle;
241 for ( iter_particle = particleCol->begin();
242 iter_particle != particleCol->end(); iter_particle++) {
243
244 if ((*iter_particle)->primaryParticle()) {
246 primaryParticleCol.push_back(mcPart);
247 }
248 }
249
250 if (primaryParticleCol.empty())
251 std::cout << "no primary particle found!" << std::endl;
252
253
254 SmartRefVector<Event::McParticle>::iterator iter_primary;
255 for ( iter_primary = primaryParticleCol.begin();
256 iter_primary != primaryParticleCol.end(); iter_primary++) {
257 if ( !((*iter_primary)->leafParticle()) )
259 }
260
261
262 StatusCode sc = m_evtSvc->registerObject("/Event/MC/McParticleCol", particleCol);
263 if(sc!=StatusCode::SUCCESS)
264 G4cout<< "Could not register McParticle collection" <<G4endl;
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289}
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
@ PRIMARY
Decayed in flight.
@ ERROR
this particle is a leaf in the particle tree
@ DECAYFLT
Decayed by generator.
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