CGEM BOSS 6.6.5.h
BESIII Offline Software System
Loading...
Searching...
No Matches
CgemSim-01-00-40/src/BesCgemSD.cc
Go to the documentation of this file.
1//-------------------------------------------------------------------------------------//
2// BOOST --- BESIII Object_Oriented Simulation Tool //
3//-------------------------------------------------------------------------------------//
4/*
5 * =====================================================================================
6 *
7 * Filename: BesCgemSD.cc
8 * Description:
9 * Conventions:
10 * gv_ : global variable
11 * lv_ : local variable used in function
12 * lvd_ : local variable double
13 * lv3_ : local variable ThreeVector
14 * m_* : normal member of class
15 * sm_* : static member of class
16 * m_M_* : class data member, Material of each layer
17 * m_N_* : class data member, Number of layers (CgemLayer,GemFoil)
18 * m_R_* : class data member, Radius of each (material) layer, and so on
19 * m_L_* : class data member, Length of each layer or hole pitch
20 * m_T_* : class data member, Thickness of each (material) layer
21 * m_W* : calss data member, Width of strip
22 * m_A_* : class data member, Angle of anode VStrip
23 * Version: CgemSim-01-00-00
24 * Created: 01/02/2014 08:59:21 AM
25 * Revision: CgemSim-00-00-01(in CMT version CgemBoss-0.0.1 written by xiuql)
26 * Compiler: gcc
27 * Author: [email protected]
28 * Organization: DG1,EPC,IHEP
29 * History:
30 * <Num> <Author> <Time> <Version> <remark>
31 * 0 juxd 20140102 01-00-00 created,
32 *
33 * =====================================================================================
34 */
35
36//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
37/* Header file: BOSS */
38#include "BesCgemSD.hh"
39#include "ReadBoostRoot.hh"
40
41//#include "GaudiKernel/SmartIF.h"
42//#include "GaudiKernel/Bootstrap.h"
43//#include "GaudiKernel/IProperty.h"
44//#include "GaudiKernel/ServiceLocatorHelper.h"
45//#include "GaudiKernel/ISvcLocator.h"
46
47/* Header file: Geant4 */
48#include "G4RunManager.hh"
49#include "G4SDManager.hh"
50#include "G4HCofThisEvent.hh"
51#include "G4Step.hh"
52#include "G4ThreeVector.hh"
53#include "G4UnitsTable.hh"
54#include "G4ios.hh"
55#include "Randomize.hh"
56#include "globals.hh"
57
58/* Header file: ROOT */
59
60/* Header file: C++ */
61#include <iostream>
62#include <iomanip>
63#include <fstream>
64#include <vector>
65#include <string>
66#include <cfloat>
67
68// Gaudi
69#include "GaudiKernel/IMessageSvc.h"
70#include "GaudiKernel/MsgStream.h"
71#include "GaudiKernel/ISvcLocator.h"
72#include "GaudiKernel/Bootstrap.h"
73#include "GaudiKernel/SmartIF.h"
74#include "GaudiKernel/Property.h"
75#include "GaudiKernel/IJobOptionsSvc.h"
76
77using namespace std;
78
79//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
80 BesCgemSD::BesCgemSD(G4String name)
81:BesSensitiveDetector(name)
82{
83 m_F_printTruth = 0;
84 m_F_outputTruth = 0;
85 myCgemTruth=NULL;
86
87 /* Define hit collection name handled by SD */
88 collectionName.insert("BesCgemHitsCollection");
89 collectionName.insert("BesCgemTruthCollection");
90
91 // get info from JobOptionsSvc.
92 ISvcLocator* iface = Gaudi::svcLocator();
93 SmartIF<IJobOptionsSvc> jos(iface->service("JobOptionsSvc"));
94 const std::vector<const Property*>* ppp = jos->getProperties("BesCgemSD");
95
96 string outputFile;
97 //outputFile += "/dat/outputMcTruth.txt";
98
99 if (ppp) {
100 //std::cout << "BesCgemSD size: " << ppp->size() << std::endl;
101 for (std::vector<const Property*>::const_iterator it = ppp->begin();
102 it != ppp->end(); ++it) {
103 if ( (*it)->name() == "outputFile") {
104 outputFile = (*it)->toString();
105 }
106 }
107 }
108 IMessageSvc* msgSvc;
109 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
110 MsgStream log(msgSvc, "BesCgemSD::BesCgemSD():BesSensitiveDetector()");
111
112 /* Save McTruth information to compare with Root RawData Information */
113 if (m_F_outputTruth == 1)
114 {
115 log<< MSG::INFO << "CgemSim::BesCgemSD::Initialize(), McTruth Hit will be output to a txt file!"
116 << endreq
117 << "outputFile Path: " << outputFile << endreq;
118
119 m_outputTruth.open(outputFile.c_str(), ios::out);
120
121 if (!m_outputTruth.good())
122 {
123 log<< MSG::ERROR << "CgemSim::BesCgemSD::Initialize(), Fail to open McTruth output file: "
124 << outputFile << endreq;
125 return ;
126 }
127 } /* End of 'if (m_outputTruth == 1)' */
128}
129
130//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131BesCgemSD::~BesCgemSD()
132{
133 m_outputTruth.close();
134}
135
136//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
137/* Invoked at the beginning of each event */
138void BesCgemSD::Initialize(G4HCofThisEvent* HCE)
139{
140 /* Instantiate hit collection */
141 m_collection_hit = new BesCgemHitsCollection(SensitiveDetectorName,collectionName[0]);
142
143 //Add by sunxh (lvi_ID_hit)
144 lvi_ID_hit = 0;
145 //end add
146 //IMessageSvc* msgSvc;
147 //Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
148 //MsgStream log(msgSvc, "BesCgemSD::Initialize()");
149 //log<< MSG::INFO << "|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||" << endreq;
150 //log<< MSG::INFO << "CgemSim::BesCgemSD::Initialize(), Begin !!!" << endreq;
151}
152
153//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
154/* MC Truth */
155void BesCgemSD::BeginOfTruthEvent(const G4Event* f_event)
156{
157 m_collection_truth = new BesCgemHitsCollection(SensitiveDetectorName,collectionName[1]);
158 lvi_ID_hit=0;
159
160}
161
162//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
163/* MC Truth */
164void BesCgemSD::EndOfTruthEvent(const G4Event* f_event)
165{
166 static G4int HLID=-1;
167 if (HLID < 0)
168 {
169 HLID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[1]);
170 }
171
172 G4HCofThisEvent *HCE = f_event->GetHCofThisEvent();
173
174 IMessageSvc* msgSvc;
175 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
176 MsgStream log(msgSvc, "BesCgemSD::EndOfTruthEvent()");
177 //log<< MSG::INFO << "CgemSim::BesCgemSD::EndOfTruthEvent(), Entries of truth collection is "
178 // << m_collection_truth->entries() << endreq;
179
180 HCE->AddHitsCollection(HLID, m_collection_truth);
181
182 /* Print and output McTruth information */
183 if (m_F_printTruth==1 or m_F_outputTruth==1)
184 {
185 /* Print McTruth information */
186 if (m_F_printTruth == 1)
187 {
188 log<< MSG::INFO << "-------------------------------------------------------------------------------"
189 << endreq;
190 log<< MSG::INFO << "CgemSim::BesCgemSD::EndOfTruthEvent(), McTruth information: " << endreq;
191 log<< MSG::INFO << left << "Track " << left << "Layer " << left << "PDG "
192 << left << setw(12) << "R_Pre "
193 << left << setw(12) << "R_Post "
194 << left << setw(40) << "pos_Pre "
195 << left << setw(40) << "pos_Post "
196 << left << setw(12) << "P_Pre "
197 << left << setw(12) << "P_Post "
198 << left << setw(16) << "StepL "
199 << left << setw(16) <<"Time "
200 << left << "Edep "
201 << endreq;
202 }
203
204
205 for (G4int ii=0; ii<m_collection_truth->entries(); ii++)
206 {
207 BesCgemHit *lv_truth = (*m_collection_truth)[ii];
208
209 /* Print McTruth information */
210 if (m_F_printTruth == 1)
211 {
212 G4cout << left << setw( 6) << lv_truth-> GetTrackID ()
213 << left << setw( 6) << lv_truth-> GetLayerID ()
214 << left << setw( 6) << lv_truth-> GetPDGCode ()
215 << left << setw(12) << lv_truth-> GetPositionOfPrePoint ().rho()
216 << left << setw(12) << lv_truth-> GetPositionOfPostPoint ().rho()
217 << left << lv_truth-> GetPositionOfPrePoint()
218 << left << lv_truth-> GetPositionOfPostPoint()
219 << left << setw(12) << lv_truth-> GetMomentumOfPrePoint().mag()
220 << left << setw(12) << lv_truth-> GetMomentumOfPostPoint().mag()
221 << left << setw(16) << lv_truth-> GetStepLength ()
222 << left << setw(16) << lv_truth-> GetGlobalTime ()
223 << left << setw(16) << lv_truth-> GetTotalEnergyDeposit ()
224 << G4endl;
225 }
226
227 /* Save McTruth information to compare with Root RawData Information */
228 if (m_F_outputTruth == 1)
229 {
230 m_outputTruth << left << setw( 5) << lv_truth->GetTrackID()
231 << left << setw( 3) << lv_truth->GetLayerID()
232 << left << setw(12) << lv_truth->GetTotalEnergyDeposit()
233 << left << setw(12) << lv_truth->GetPositionOfPrePoint().x()
234 << left << setw(12) << lv_truth->GetPositionOfPrePoint().y()
235 << left << setw(12) << lv_truth->GetPositionOfPrePoint().z()
236 << left << setw(12) << lv_truth->GetPositionOfPostPoint().x()
237 << left << setw(12) << lv_truth->GetPositionOfPostPoint().y()
238 << left << setw(12) << lv_truth->GetPositionOfPostPoint().z()
239 << left << setw(12) << lv_truth->GetMomentumOfPrePoint().x()
240 << left << setw(12) << lv_truth->GetMomentumOfPrePoint().y()
241 << left << setw(12) << lv_truth->GetMomentumOfPrePoint().z()
242 << left << setw(12) << lv_truth->GetMomentumOfPostPoint().x()
243 << left << setw(12) << lv_truth->GetMomentumOfPostPoint().y()
244 << left << setw(12) << lv_truth->GetMomentumOfPostPoint().z()
245 << G4endl;
246 }
247 } /* End of 'for (G4int ii=0; ii<m_collection_truth->entries(); ii++)' */
248
249 //m_outputTruth.close();
250
251 } /* End of 'if (m_F_printTruth==1 or m_F_outputTruth==1)' */
252}
253//************************************************************************
254//************************************************************************
255//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
256/* Invoked at the beginning of every step in the SD volume */
257G4bool BesCgemSD::ProcessHits(G4Step* f_step, G4TouchableHistory* f_touchable)
258{
259 bool printFlag=false;
260
261 IMessageSvc* msgSvc;
262 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
263 MsgStream log(msgSvc, "BesCgemSD::ProcessHits()");
264
265 /* Event -> Track -> Step */
266 G4Track *lv_track = f_step->GetTrack();
267
268 /* Check tracks */
269 /* Skip neutral tracks */
270 G4double charge = lv_track->GetDefinition()->GetPDGCharge();
271 if (charge == 0.)
272 {
273 if(printFlag)
274 log<< MSG::WARNING << "CgemSim::BesCgemSD::ProcessHits(), Track is neutral, skip the track" << endreq;
275 return false;
276 }
277 G4int lvi_pdg_code = lv_track->GetDefinition()->GetPDGEncoding();
278
279 /* Global Time (time since the current event began ) */
280 G4double lvd_global_time = lv_track->GetGlobalTime();
281 if (isnan(lvd_global_time))
282 {
283 if(printFlag)
284 log<< MSG::WARNING << "CgemSim::BesCgemSD::ProcessHits(), Gloabal time is nan" << endreq;
285 return false;
286 }
287 /* MDC T window is 2 microsecond */
288 if (lvd_global_time > 2000)
289 {
290 if(printFlag)
291 log<< MSG::WARNING << "CgemSim::BesCgemSD::ProcessHits(), Time window is >2000" << endreq;
292 return false;
293 }
294
295 /* Skip no energy deposit tracks */
296 /* if (lvd_E_deposit == 0.)
297 {
298 log<< MSG::WARNING << "CgemSim::BesCgemSD::ProcessHits(), Energy deposition is 0,lose track="
299 << primTrackIdx << endreq;
300 return false;
301 }*/
302
303 /* Skip no length tracks */
304 G4double lvd_L_step = f_step->GetStepLength();
305 if (lvd_L_step == 0.)
306 {
307 if(printFlag)
308 log<< MSG::WARNING << "CgemSim::BesCgemSD::ProcessHits(), Step length is 0, skip track" << endreq;
309 return false;
310 }
311
312 /* A G4Step object consists of two points */
313 G4StepPoint *lv_step_pre = f_step->GetPreStepPoint();
314 G4StepPoint *lv_step_post = f_step->GetPostStepPoint();
315 G4TrackStatus currentTrackStatus = lv_track->GetTrackStatus();
316
317 /* Skip delta electron, added by sunxh */
318 G4int lvi_RdtElectron = 0;
319 G4String process_name="Generator";
320 if (NULL != lv_track->GetCreatorProcess()){
321 process_name = lv_track->GetCreatorProcess()->GetProcessName();
322 if (process_name == "hIoni" || process_name == "ionIoni" || process_name == "conv" || process_name == "compt" || process_name == "muIoni" || process_name == "eIoni" || process_name == "phot" || process_name == "muMinusCaptureAtRest" ) lvi_RdtElectron = 1;
323 //if(lvi_RdtElectron != 1&&lv_track->GetDefinition()->GetPDGEncoding()!=22){
324 // G4cout << "process_name = " << process_name << G4endl;
325 // G4cout << "PDG = " << lv_track->GetDefinition()->GetPDGEncoding()<<"\tTrack ID = "<<lv_track->GetTrackID() <<"\tParent ID = "<<lv_track->GetParentID() << G4endl;
326 //}
327 //cout<<"particle pdg="<<lvi_pdg_code<<" with CreatorProcess "<<process_name<<endl;
328 }
329 if(lvi_RdtElectron) return false;
330
331 //else {
332 // G4cout << "process_name = null" << G4endl;
333 // G4cout << "PDG = " << lv_track->GetDefinition()->GetPDGEncoding()<<"\tTrack ID = "<<lv_track->GetTrackID() <<"\tParent ID = "<<lv_track->GetParentID() << G4endl;
334 //}
335 // G4cout << "PDG = " << lv_track->GetDefinition()->GetPDGEncoding() << G4endl;
336 // else G4cout << "No process_name PDG = " << lv_track->GetDefinition()->GetPDGEncoding() << G4endl;
337 // if (abs(lv_track->GetDefinition()->GetPDGEncoding()) == 11 ) lvi_RdtElectron = 1;
338 // else{
339 // if (abs(lv_track->GetDefinition()->GetPDGEncoding()) == 11) {
340 // G4cout << "No process_name PDG = " << lv_track->GetDefinition()->GetPDGEncoding() << G4endl;
341 // }
342 // }
343 //End Add
344
345 /* Local variables used to get Hit information */
346 /* Get track ID */
347 G4int primTrackIdx = 0;// trk index for the primary track, comments by llwang
348 G4int G4TrackId = 0;// G4 trk id for the current track
349 G4int primG4TrackId = 0;// G4 trk id for the primary track
350
351 /* In the following we use 'SetTrackID(primTrackIdx)' */
352 //primTrackIdx = lv_track->GetTrackID();
353 //GetCurrentTrackIndex(G4TrackId, primG4TrackId);
354
355 //Add by sunxh
356 G4TrackId = lv_track->GetTrackID();
357 GetCurrentTrackIndex(primTrackIdx, primG4TrackId);
358 // G4cout<<"PDG = "<<lv_track->GetDefinition()->GetPDGEncoding()<<"\tprocess_name = "<<process_name<<"\tlvi_ID_truth = " << G4TrackId<<"\tprimTrackIdx = "<<primTrackIdx<<"\tlvi_ID_G4Track = "<<primG4TrackId<<G4endl;
359 //end
360
361 G4int lvi_ID_layer = 0;
362 //G4int lvi_pdg_code = lv_track->GetDefinition()->GetPDGEncoding();
363 G4int lvi_ID_parent = lv_track->GetParentID();// G4 trk id for the parent
364 G4double lvd_E_deposit = f_step->GetTotalEnergyDeposit();
365 /* Positions in the global coordinate system */
366 G4ThreeVector lv3_XYZ_pre = lv_step_pre ->GetPosition();
367 G4ThreeVector lv3_XYZ_post = lv_step_post->GetPosition();
368
369 G4ThreeVector lv3_P_pre = lv_step_pre ->GetMomentum();
370 G4ThreeVector lv3_P_post = lv_step_post->GetMomentum();
371
372 if(printFlag) {
373 log<< MSG::INFO << "..............................................................................." << endreq;
374 log<< MSG::INFO << "CgemSim::BesCgemSD::ProcessHits(), Process a hit at track="
375 << primTrackIdx << endreq;
376 }
377
378
379
380 /* Check and Save Hits and MCTruth Hits */
381 /* Get pre step point status and its volume */
382 /* Volume depth: 0-Gap_D, 1-CgemLayer, 2-Cgem, 3-BES */
383 G4StepStatus lv_status_pre = lv_step_pre ->GetStepStatus();
384 G4StepStatus lv_status_post = lv_step_post->GetStepStatus();
385
386 /* Retrieve a touchable by creating a handle for it */
387 const G4VTouchable *lv_TOUCHABLE = lv_step_pre ->GetTouchable();
388 /* Get the current volume(where the step has just gone through) and its logical name */
389 G4String lvs_volume_name = lv_TOUCHABLE ->GetVolume(0)->GetLogicalVolume()->GetName();
390
391 /* Save Hits in SD Geometry only! */
392 //if ((lv_status_pre == fGeomBoundary) && (lvs_volume_name == "Gap_D_Logic"))
393
394
395 //if (lvs_volume_name == "Gap_D_logic")
396 if (lvs_volume_name == "Gap_D_logic0" or lvs_volume_name == "Gap_D_logic1" or lvs_volume_name == "Gap_D_logic2")
397 {
398 G4double trkLen = lv_track->GetTrackLength();// track length at post point?
399
400 lvi_ID_layer = lv_TOUCHABLE->GetVolume(1)->GetCopyNo();
401
402 if(printFlag) {
403
404 log<< MSG::INFO << "CgemSim::BesCgemSD::ProcessHits(), Track enters the SD volume: "
405 << lvs_volume_name << endreq;
406
407 /* CgemLayer copy number */
408
409 // KeyType lv_key(primTrackIdx, lvi_ID_layer);
410
411 /* Create Hits */
412 log<< MSG::INFO << "CgemSim::BesCgemSD::ProcessHits(), Create a hit at layer="
413 << lvi_ID_layer << " track=" << primTrackIdx << endreq;
414 }
415
416 BesCgemHit *l_new_hit = new BesCgemHit();
417
418 //Add by sunxh
419 l_new_hit->SetRdtEl ( lvi_RdtElectron );
420 l_new_hit->SetHitID ( lvi_ID_hit );
421 //end
422 l_new_hit->SetTrackID ( primTrackIdx );// primary track idx in McParticle Collection
423 G4int isSecondary = (G4TrackId==primG4TrackId)? 0:1;
424 //l_new_hit->SetCurrentTrackID ( G4TrackId );
425 l_new_hit->SetIsSecondary ( isSecondary );
426 l_new_hit->SetLayerID ( lvi_ID_layer );
427 l_new_hit->SetPDGCode ( lvi_pdg_code );
428 l_new_hit->SetParentID ( lvi_ID_parent );// G4 trk id for the parent
429 l_new_hit->SetGlobalTime ( lvd_global_time );
430 l_new_hit->SetTotalEnergyDeposit ( lvd_E_deposit );
431 l_new_hit->SetStepLength ( lvd_L_step );
432
433 // --> modified by wulh for alignment on Apr 5, 2019 begin:
434 double x_pre_glo = lv3_XYZ_pre.x();
435 double y_pre_glo = lv3_XYZ_pre.y();
436 double z_pre_glo = lv3_XYZ_pre.z();
437 double x_post_glo = lv3_XYZ_post.x();
438 double y_post_glo = lv3_XYZ_post.y();
439 double z_post_glo = lv3_XYZ_post.z();
440 double r_pre = sqrt(x_pre_glo*x_pre_glo+y_pre_glo*y_pre_glo);
441 double r_post = sqrt(x_post_glo*x_post_glo+y_post_glo*y_post_glo);
442
443 CgemGeoAlign* cgemAlign = m_cgemGeomSvc->getAlignPtr();
444 HepPoint3D pos_pre(x_pre_glo, y_pre_glo, z_pre_glo);
445 HepPoint3D pos_post(x_post_glo, y_post_glo, z_post_glo);
446 StraightLine line_glo(pos_pre, pos_post);
447 StraightLine line_loc = cgemAlign->StraightLineConversion(lvi_ID_layer, line_glo);
448 HepPoint3D crossPoint_pre;
449 HepPoint3D crossPoint_post;
450
451 double rDrift_inn = m_cgemGeomSvc->getCgemLayer(lvi_ID_layer)->getInnerROfGapD();
452 double rDrift_out = m_cgemGeomSvc->getCgemLayer(lvi_ID_layer)->getOuterROfGapD();
453
454 double distLineGlo = fabs(line_glo.dr());
455 if(fabs(distLineGlo) > rDrift_inn){
456 crossPoint_pre = line_loc.xAtR(r_post, -1);
457 crossPoint_post = line_loc.xAtR(r_post, 1);
458 } else if(r_pre < r_post){
459 crossPoint_pre = line_loc.xAtR(r_pre, 1);
460 crossPoint_post = line_loc.xAtR(r_post, 1);
461 } else{
462 crossPoint_pre = line_loc.xAtR(r_pre, -1);
463 crossPoint_post = line_loc.xAtR(r_post, -1);
464 }
465
466 G4ThreeVector lv3_XYZ_pre_align(crossPoint_pre.x(), crossPoint_pre.y(), crossPoint_pre.z());
467 G4ThreeVector lv3_XYZ_post_align(crossPoint_post.x(), crossPoint_post.y(), crossPoint_post.z());
468
469 l_new_hit->SetPositionOfPrePointAlign ( lv3_XYZ_pre_align );
470 l_new_hit->SetPositionOfPostPointAlign ( lv3_XYZ_post_align );
471 // <-- modified by wulh for alignment on Apr 5, 2019 --- end.
472
473
474 l_new_hit->SetPositionOfPrePoint ( lv3_XYZ_pre );
475 l_new_hit->SetPositionOfPostPoint ( lv3_XYZ_post );
476 l_new_hit->SetMomentumOfPrePoint ( lv3_P_pre );
477 l_new_hit->SetMomentumOfPostPoint ( lv3_P_post );
478 l_new_hit->SetCreatorProcess ( process_name );
479 l_new_hit->SetFlightLengthPostPoint(trkLen);
480
481 m_collection_hit->insert(l_new_hit);
482
483 //Add by sunxh
484 hit_ID_vector.push_back(lvi_ID_hit);
485 lvi_ID_hit++;
486 //end
487
488
489
490 /* Create Truth Hit */
491 /* MCTruth Hit only contains the point input and output the volume (not step by step) */
492 /* Meaning of 'MCTruth hit' parameters is not the same as 'hit' */
493 /* All steps in SD volume will be merged to 1 step */
494 /* Check over that the particle has just entered in the current volume */
495 /* It is at the first step in the volume; the preStepPoint is at the boundary */
496
497 KeyType lv_key(primTrackIdx, lvi_ID_layer);
498
499 if (m_collection_truth)
500 {
501 int newVer=1;
502
503 if(newVer==0) {
504 m_map_truth_it = m_map_truth.find(lv_key);
505 if (m_map_truth_it == m_map_truth.end())
506 {
507 log<< MSG::INFO << "CgemSim::BesCgemSD::ProcessHits(), Create a Truth at layer="
508 << lvi_ID_layer << " track=" << primTrackIdx << endreq;
509
510 BesCgemHit *l_new_truth = new BesCgemHit();
511
512 l_new_truth->SetTrackID ( primTrackIdx ); /* SAME */
513 l_new_truth->SetLayerID ( lvi_ID_layer ); /* SAME */
514 l_new_truth->SetPDGCode ( lvi_pdg_code ); /* SAME */
515 l_new_truth->SetParentID ( lvi_ID_parent ); /* SAME */
516 //l_new_truth->SetCurrentTrackID ( G4TrackId );
517 l_new_truth->SetGlobalTime ( lvd_global_time ); /* DIFF */
518 l_new_truth->SetTotalEnergyDeposit ( lvd_E_deposit ); /* SAME */
519 l_new_truth->SetStepLength ( lvd_L_step ); /* SAME */
520 l_new_truth->SetPositionOfPrePoint ( lv3_XYZ_pre ); /* SAME */
521 l_new_truth->SetPositionOfPostPoint ( lv3_XYZ_post ); /* SAME */
522 l_new_truth->SetPositionOfPrePointAlign(lv3_XYZ_pre_align); /* SAME */
523 l_new_truth->SetPositionOfPostPointAlign(lv3_XYZ_post_align); /* SAME */
524 l_new_truth->SetMomentumOfPrePoint ( lv3_P_pre ); /* SAME */
525 l_new_truth->SetMomentumOfPostPoint ( lv3_P_post ); /* SAME */
526
527 //Add by sunxh
528 int tmp[2000];
529 for(int ii =0; ii < hit_ID_vector.size();ii++){
530 tmp[ii] = hit_ID_vector[ii];
531 }
532 l_new_truth->AddIdentifier(tmp,hit_ID_vector.size());
533 hit_ID_vector.clear();
534 //end
535
536
537 m_collection_truth->insert(l_new_truth);
538
539 m_map_truth[lv_key] = m_collection_truth->entries() - 1;
540
541 if (lv_step_pre->GetStepStatus() == fGeomBoundary)
542 {
543 log<< MSG::INFO << "CgemSim::BesCgemSD::ProcessHits(), Truth starts from boundry!!!"
544 << endreq;
545 }
546 else
547 {
548 log<< MSG::INFO << "CgemSim::BesCgemSD::ProcessHits(), Truth doesn't start from boundry!!!"
549 << endreq;
550 }
551 }//if (m_map_truth_it == m_map_truth.end())
552 else
553 {
554 if(printFlag)
555 log<< MSG::INFO << "CgemSim::BesCgemSD::ProcessHits(), Add a step to Truth at layer="
556 << lvi_ID_layer << " track=" << primTrackIdx << endreq;
557
558 G4int lvi_ID_truth = (*m_map_truth_it).second;
559 G4int primTrackIdx_pre = (*m_collection_truth)[lvi_ID_truth]->GetTrackID();
560 G4double lvd_E_pre = (*m_collection_truth)[lvi_ID_truth]->GetTotalEnergyDeposit();
561 G4double lvd_L_step_pre = (*m_collection_truth)[lvi_ID_truth]->GetStepLength();
562
563 (*m_collection_truth)[lvi_ID_truth]->SetTotalEnergyDeposit(lvd_E_deposit + lvd_E_pre);
564 (*m_collection_truth)[lvi_ID_truth]->SetStepLength(lvd_L_step + lvd_L_step_pre);
565 (*m_collection_truth)[lvi_ID_truth]->SetPositionOfPostPoint(lv3_XYZ_post);
566 (*m_collection_truth)[lvi_ID_truth]->SetPositionOfPostPointAlign(lv3_XYZ_post_align);
567 (*m_collection_truth)[lvi_ID_truth]->SetMomentumOfPostPoint(lv3_P_post);
568 }
569
570 /* Check over that the particle is leaving the current volume */
571 /* It is at the last step in the volume; the postStepPoint is at the boundary */
572 if(printFlag)
573 if (lv_step_post->GetStepStatus() == fGeomBoundary)
574 {
575 log<< MSG::INFO << "CgemSim::BesCgemSD::ProcessHits(), Finish a Truth at layer="
576 << lvi_ID_layer << " track=" << primTrackIdx << endreq;
577 }
578 }// old version
579 else
580 {
581 if(myCgemTruth==NULL)
582 {
583 myCgemTruth= new BesCgemHit();
584 if(printFlag)
585 cout<<"creat a new myCgemTruth, prePos is Boundary: "<<(lv_status_pre==fGeomBoundary)<<" trkID="<<primTrackIdx<<", LayerID="<<lvi_ID_layer<<endl;
586 // myCgemTruth->SetTrackID ( G4TrackId); /* SAME */
587 myCgemTruth->SetTrackID ( primTrackIdx); /* SAME */
588 myCgemTruth->SetLayerID ( lvi_ID_layer ); /* SAME */
589 myCgemTruth->SetPDGCode ( lvi_pdg_code ); /* SAME */
590 myCgemTruth->SetParentID ( lvi_ID_parent ); /* SAME */
591 myCgemTruth->SetGlobalTime ( lvd_global_time ); /* DIFF */
592 myCgemTruth->SetTotalEnergyDeposit ( lvd_E_deposit ); /* SAME */
593 myCgemTruth->SetStepLength ( lvd_L_step ); /* SAME */
594 myCgemTruth->SetPositionOfPrePoint ( lv3_XYZ_pre ); /* SAME */
595 myCgemTruth->SetPositionOfPostPoint ( lv3_XYZ_post ); /* SAME */
596 myCgemTruth->SetPositionOfPrePointAlign(lv3_XYZ_pre_align); /* SAME */
597 myCgemTruth->SetPositionOfPostPointAlign(lv3_XYZ_post_align); /* SAME */
598 myCgemTruth->SetMomentumOfPrePoint ( lv3_P_pre ); /* SAME */
599 myCgemTruth->SetMomentumOfPostPoint ( lv3_P_post ); /* SAME */
600 myCgemTruth->SetCreatorProcess ( process_name );
601 myCgemTruth->SetRdtEl ( lvi_RdtElectron );
602 myCgemTruth->SetIsSecondary (isSecondary );
603 myCgemTruth->SetFlightLengthPostPoint(trkLen);
604 //double dphi=lv3_XYZ_pre.phi()-lv3_XYZ_post.phi();
605 //double dX=dphi*87.5;
606 //if(lvi_ID_layer==0) cout<<"dX="<<dX<<" , nXStrips="<<dX/0.65<<endl;
607 }
608 else
609 {
610 G4double lvd_E_pre = myCgemTruth->GetTotalEnergyDeposit();
611 G4double lvd_L_step_pre = myCgemTruth->GetStepLength();
612 G4ThreeVector pos_pre_last = myCgemTruth->GetPositionOfPrePoint();
613 G4ThreeVector pos_post_last = myCgemTruth->GetPositionOfPostPoint();
614 //cout<<"add CGEM truth, last hit: "<<pos_pre_last<<" to "<<pos_post_last<<endl;
615 //cout<<" new hit: "<<lv3_XYZ_pre <<" to "<<lv3_XYZ_post<<endl;
616 myCgemTruth->SetTotalEnergyDeposit(lvd_E_deposit + lvd_E_pre);
617 myCgemTruth->SetStepLength(lvd_L_step + lvd_L_step_pre);
618 myCgemTruth->SetPositionOfPostPoint(lv3_XYZ_post);
619 myCgemTruth->SetPositionOfPostPointAlign(lv3_XYZ_post_align);
620 myCgemTruth->SetMomentumOfPostPoint(lv3_P_post);
621 myCgemTruth->SetFlightLengthPostPoint(trkLen);
622 if(printFlag)
623 cout<<"update myCgemTruth: prePos is Boundary: "<<(lv_status_pre==fGeomBoundary)<<", postPos is Boundary:"<<(lv_status_post==fGeomBoundary)<<", trkID="<<primTrackIdx<<", LayerID="<<lvi_ID_layer<<endl;
624 }
625 if( myCgemTruth && (currentTrackStatus==fStopAndKill || lv_status_post==fGeomBoundary) )
626 {
627 //Add by sunxh
628 int tmp[2000];
629 for(int ii =0; ii < hit_ID_vector.size();ii++){
630 tmp[ii] = hit_ID_vector[ii];
631 }
632 myCgemTruth->AddIdentifier(tmp,hit_ID_vector.size());
633 //TArrayI hit_tmp=myCgemTruth->GetIdentifier();
634 //end
635
636 m_collection_truth->insert(myCgemTruth);
637 myCgemTruth=NULL;
638
639 //Add by sunxh
640 hit_ID_vector.clear();
641 //end
642
643 if(printFlag)
644 cout<<"save a myCgemTruth: prePos is Boundary: "<<(lv_status_pre==fGeomBoundary)<<", postPos is Boundary:"<<(lv_status_post==fGeomBoundary)<<", trkID="<<primTrackIdx<<", LayerID="<<lvi_ID_layer<<endl;
645 }
646 }//if(newVer==0)
647 } /* End of 'if (m_collection_truth)' */
648 } /* End of 'if (lvs_volume_name == "Gap_D_Logic")' */
649
650 return true;
651}
652
653//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
654/* Invoked at the end of processing an event (even the event is aborted) */
655void BesCgemSD::EndOfEvent(G4HCofThisEvent* HCE)
656{
657 IMessageSvc* msgSvc;
658 Gaudi::svcLocator() -> service("MessageSvc", msgSvc);
659 MsgStream log(msgSvc, "BesCgemSD::EndOfEvent()");
660
661 bool printFlag=false;
662 m_map_truth.clear();
663
664 /* Get the unique ID for this collection */
665 /* GetCollectionID() is a heavy operation, shouldn't be called for every event. */
666 /* It's available after SD is constructed and registed to G4SDManager. */
667 /* It can't be invoked in the constructor. */
668 static G4int HCID = -1;
669 if (HCID < 0)
670 {
671 HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]);
672 }
673
674 /* Attach hit collection to G4HCofThisEvent object */
675 if(printFlag)
676 {
677 log<< MSG::INFO << "-------------------------------------------------------------------------------" << endreq;
678 log<< MSG::INFO << "CgemSim::BesCgemSD::EndOfEvent(), Entries of hits collection is "
679 << m_collection_hit->entries() << endreq;
680 }
681
682 HCE->AddHitsCollection(HCID, m_collection_hit);
683
684 if (verboseLevel > 0)
685 {
686 m_collection_hit->PrintAllHits();
687 }
688}
689
690//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
BesCgemHitsCollection
G4THitsCollection< BesCgemHit > BesCgemHitsCollection
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:177
msgSvc
IMessageSvc * msgSvc()
Definition ServiceAccessor.h:13
BesCgemHit::SetTotalEnergyDeposit
void SetTotalEnergyDeposit(G4double f_E_deposit)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:82
BesCgemHit::SetGlobalTime
void SetGlobalTime(G4double f_global_time)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:81
BesCgemHit::GetPositionOfPrePoint
G4ThreeVector GetPositionOfPrePoint() const
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:118
BesCgemHit::AddIdentifier
void AddIdentifier(G4int f_ID_Identifier[2000], G4int N_dim)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:95
BesCgemHit::GetMomentumOfPostPoint
G4ThreeVector GetMomentumOfPostPoint() const
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:123
BesCgemHit::SetMomentumOfPostPoint
void SetMomentumOfPostPoint(G4ThreeVector f_P_post)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:89
BesCgemHit::SetPDGCode
void SetPDGCode(G4int f_pdg_code)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:79
BesCgemHit::SetLayerID
void SetLayerID(G4int f_ID_layer)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:78
BesCgemHit::SetRdtEl
void SetRdtEl(G4int f_RdtElectron)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:73
BesCgemHit::SetParentID
void SetParentID(G4int f_ID_parent)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:80
BesCgemHit::SetFlightLengthPostPoint
void SetFlightLengthPostPoint(G4double len)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:101
BesCgemHit::SetTrackID
void SetTrackID(G4int f_ID_track)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:76
BesCgemHit::GetPositionOfPostPoint
G4ThreeVector GetPositionOfPostPoint() const
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:119
BesCgemHit::SetPositionOfPrePoint
void SetPositionOfPrePoint(G4ThreeVector f_XYZ_pre)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:84
BesCgemHit::SetMomentumOfPrePoint
void SetMomentumOfPrePoint(G4ThreeVector f_P_pre)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:88
BesCgemHit::GetMomentumOfPrePoint
G4ThreeVector GetMomentumOfPrePoint() const
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:122
BesCgemHit::SetCreatorProcess
void SetCreatorProcess(G4String proName)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:97
BesCgemHit::SetPositionOfPostPointAlign
void SetPositionOfPostPointAlign(G4ThreeVector f_XYZ_post)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:87
BesCgemHit::SetPositionOfPrePointAlign
void SetPositionOfPrePointAlign(G4ThreeVector f_XYZ_pre)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:86
BesCgemHit::SetPositionOfPostPoint
void SetPositionOfPostPoint(G4ThreeVector f_XYZ_post)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:85
BesCgemHit::SetStepLength
void SetStepLength(G4double f_L_step)
Definition CgemSim-01-00-40/CgemSim/BesCgemHit.hh:83
BesCgemSD::EndOfTruthEvent
void EndOfTruthEvent(const G4Event *)
Definition CgemSim-01-00-40/src/BesCgemSD.cc:164
BesCgemSD::EndOfEvent
void EndOfEvent(G4HCofThisEvent *)
Definition CgemSim-01-00-40/src/BesCgemSD.cc:655
BesCgemSD::ProcessHits
G4bool ProcessHits(G4Step *, G4TouchableHistory *)
Definition CgemSim-01-00-40/src/BesCgemSD.cc:257
BesCgemSD::BeginOfTruthEvent
void BeginOfTruthEvent(const G4Event *)
Definition CgemSim-01-00-40/src/BesCgemSD.cc:155
BesCgemSD::Initialize
void Initialize(G4HCofThisEvent *)
Definition CgemSim-01-00-40/src/BesCgemSD.cc:138
BesSensitiveDetector::GetCurrentTrackIndex
void GetCurrentTrackIndex(G4int &trackIndex, G4int &g4TrackId) const
Definition BesSensitiveDetector.cc:30
CgemGeoAlign::StraightLineConversion
StraightLine StraightLineConversion(int layer_vir, StraightLine lineOriginal)
Definition CgemGeoAlign.cxx:59
CgemGeoLayer::getOuterROfGapD
double getOuterROfGapD() const
Definition CgemGeoLayer.h:104
CgemGeoLayer::getInnerROfGapD
double getInnerROfGapD() const
Definition CgemGeoLayer.h:103
CgemGeomSvc::getCgemLayer
CgemGeoLayer * getCgemLayer(int i) const
Definition CgemGeomSvc.h:48
CgemGeomSvc::getAlignPtr
CgemGeoAlign * getAlignPtr() const
Definition CgemGeomSvc.h:69
StraightLine::xAtR
HepPoint3D xAtR(double R, int direction=1) const
Definition StraightLine.cxx:100
StraightLine::dr
double dr(void) const
returns an element of parameters.
Definition StraightLine.h:45