BOSS 7.0.9
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | List of all members
BesMdcSD Class Reference

#include <BesMdcSD.hh>

+ Inheritance diagram for BesMdcSD:

Public Member Functions

 BesMdcSD (G4String)
 
 ~BesMdcSD ()
 
void Initialize (G4HCofThisEvent *)
 
G4bool ProcessHits (G4Step *, G4TouchableHistory *)
 
void EndOfEvent (G4HCofThisEvent *)
 
void BeginOfTruthEvent (const G4Event *)
 
void EndOfTruthEvent (const G4Event *)
 
G4double Distance (G4int, G4int, G4ThreeVector, G4ThreeVector, G4ThreeVector &, G4double &)
 
void dedxFuncInti (void)
 
- Public Member Functions inherited from BesSensitiveDetector
 BesSensitiveDetector (const G4String name)
 
virtual ~BesSensitiveDetector ()
 
virtual void BeginOfTruthEvent (const G4Event *)
 
virtual void EndOfTruthEvent (const G4Event *)
 
virtual void BeginOfTrack (const G4Track *)
 
virtual void EndOfTrack (const G4Track *)
 

Additional Inherited Members

- Protected Member Functions inherited from BesSensitiveDetector
void GetCurrentTrackIndex (G4int &trackIndex, G4int &g4TrackId) const
 

Detailed Description

Definition at line 36 of file BesMdcSD.hh.

Constructor & Destructor Documentation

◆ BesMdcSD()

BesMdcSD::BesMdcSD ( G4String  name)

Definition at line 28 of file BesMdcSD.cc.

29:BesSensitiveDetector(name)
30{
31 collectionName.insert("BesMdcHitsCollection");
32 collectionName.insert("BesMdcTruthCollection");
33
34 mdcGeoPointer=BesMdcGeoParameter::GetGeo();
35 mdcCalPointer=new BesMdcCalTransfer;
36
37 IMdcGeomSvc* ISvc;
38 StatusCode sc=Gaudi::svcLocator()->service("MdcGeomSvc", ISvc);
39 if (!sc.isSuccess())
40 std::cout<<"BesMdcSD::Could not open MdcGeomSvc"<<std::endl;
41 mdcGeomSvc=dynamic_cast<MdcGeomSvc *>(ISvc);
42
43 IG4Svc* tmpSvc;
44 sc=Gaudi::svcLocator()->service("G4Svc", tmpSvc);
45 if (!sc.isSuccess())
46 G4cout <<" MdcSD::Error,could not open G4Svc"<<G4endl;
47 m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc);
48
49 if(m_G4Svc->GetMdcDedxFlag()==1){
50 G4cout <<" MdcSD: Use sampled dedx instead of Geant4 value"<<G4endl;
51 dedxFuncInti();
52 }
53
54 //get DedxSimSvc
55 IDedxSimSvc* tmp_dedxSimSvc;
56 sc = Gaudi::svcLocator()->service("DedxSimSvc", tmp_dedxSimSvc, true);
57 if (!sc.isSuccess())
58 {
59 std::cout << "Could not get DedxSimSvc" << std::endl;
60 }
61 m_pDedxSimSvc = dynamic_cast<DedxSimSvc*>(tmp_dedxSimSvc);
62 //get CalibCurSvc
63 IDedxCurSvc* tmp_dedxCurSvc;
64 sc = Gaudi::svcLocator()->service("DedxCurSvc", tmp_dedxCurSvc, true);
65 if (!sc.isSuccess())
66 {
67 std::cout << "Could not get DedxCurSvc"
68 << std::endl;
69 }
70 m_pDedxCurSvc = dynamic_cast<DedxCurSvc*>(tmp_dedxCurSvc);
71
72 if(m_G4Svc->MdcRootFlag())
73 {
74 m_tupleMdc = m_G4Svc->GetTupleMdc();
75 sc = m_tupleMdc->addItem("betaGamma",m_betaGamma);
76 sc = m_tupleMdc->addItem("fitval",m_fitval);
77 sc = m_tupleMdc->addItem("dedx",m_dedx);
78 sc = m_tupleMdc->addItem("de",m_de);
79 //sc = m_tupleMdc->addItem("length",m_length);
80 //sc = m_tupleMdc->addItem("random",m_random);
81 sc = m_tupleMdc->addItem("charge", m_charge);
82 sc = m_tupleMdc->addItem("costheta", m_costheta);
83 }
84}
BesMdcGeoParameter::GetGeo
static BesMdcGeoParameter * GetGeo(void)
Definition: BesMdcGeoParameter.cc:34
BesMdcSD::dedxFuncInti
void dedxFuncInti(void)
Definition: BesMdcSD.cc:412
G4Svc
Definition: G4Svc.h:32
G4Svc::MdcRootFlag
bool MdcRootFlag()
Definition: G4Svc.h:117
G4Svc::GetTupleMdc
NTuple::Tuple * GetTupleMdc()
Definition: G4Svc.h:96
G4Svc::GetMdcDedxFlag
int GetMdcDedxFlag()
Definition: G4Svc.h:94
IG4Svc
Definition: IG4Svc.h:30

◆ ~BesMdcSD()

BesMdcSD::~BesMdcSD ( )

Definition at line 86 of file BesMdcSD.cc.

86 {
87}

Member Function Documentation

◆ BeginOfTruthEvent()

void BesMdcSD::BeginOfTruthEvent ( const G4Event *  evt)
virtual

Reimplemented from BesSensitiveDetector.

Definition at line 107 of file BesMdcSD.cc.

108{
109 truthCollection = new BesMdcHitsCollection
110 (SensitiveDetectorName,collectionName[1]);
111 // G4cout<<" begin event "<<evt->GetEventID()<<G4endl;
112}
BesMdcHitsCollection
G4THitsCollection< BesMdcHit > BesMdcHitsCollection
Definition: BesMdcHit.hh:78

◆ dedxFuncInti()

void BesMdcSD::dedxFuncInti ( void  )

Definition at line 412 of file BesMdcSD.cc.

413{
414 dEdE_mylanfunc = new TF1("dEdE_mylanfunc",
415 "[3]*TMath::Exp([2]*((x[0]-[0])/[1]+TMath::Exp(-1*((x[0]-[0])/[1]))))",
416 0,
417 7500);
418 //dEdE_mylanfunc->SetParameters(2009.35,559.776,-1.0932,6327.38);
419 dEdE_mylanfunc->SetParNames("MPV","Sigma","constant1","constant2");
420}

Referenced by BesMdcSD().

◆ Distance()

G4double BesMdcSD::Distance ( G4int  layerId,
G4int  cellId,
G4ThreeVector  pointIn,
G4ThreeVector  pointOut,
G4ThreeVector &  hitPosition,
G4double &  transferT 
)

Definition at line 351 of file BesMdcSD.cc.

352{
353 //For two lines r=r1+t1.v1 & r=r2+t2.v2
354 //the closest approach is d=|(r2-r1).(v1 x v2)|/|v1 x v2|
355 //the point where closest approach are
356 //t1=(v1 x v2).[(r2-r1) x v2]/[(v1 x v2).(v1 x v2)]
357 //t2=(v1 x v2).[(r2-r1) x v1]/[(v1 x v2).(v1 x v2)]
358 //if v1 x v2=0 means two lines are parallel
359 //d=|(r2-r1) x v1|/|v1|
360
361 G4double t1,distance,dInOut,dHitIn,dHitOut;
362 //Get wirepoint @ endplate
363 G4ThreeVector east=mdcGeomSvc->Wire(layerId,cellId)->Backward();
364 G4ThreeVector west=mdcGeomSvc->Wire(layerId,cellId)->Forward();
365 G4ThreeVector wireLine=east-west;
366 G4ThreeVector hitLine=pointOut-pointIn;
367
368 G4ThreeVector hitXwire=hitLine.cross(wireLine);
369 G4ThreeVector wire2hit=east-pointOut;
370 //Hitposition is the position on hit line where closest approach
371 //of two lines, but it may out the area from pointIn to pointOut
372 if(hitXwire.mag()==0){
373 distance=wireLine.cross(wire2hit).mag()/wireLine.mag();
374 hitPosition=pointIn;
375 }else{
376 t1=hitXwire.dot(wire2hit.cross(wireLine))/hitXwire.mag2();
377 hitPosition=pointOut+t1*hitLine;
378
379 dInOut=(pointOut-pointIn).mag();
380 dHitIn=(hitPosition-pointIn).mag();
381 dHitOut=(hitPosition-pointOut).mag();
382 if(dHitIn<=dInOut && dHitOut<=dInOut){ //Between point in & out
383 distance=fabs(wire2hit.dot(hitXwire)/hitXwire.mag());
384 }else if(dHitOut>dHitIn){ // out pointIn
385 distance=wireLine.cross(pointIn-east).mag()/wireLine.mag();
386 hitPosition=pointIn;
387 }else{ // out pointOut
388 distance=wireLine.cross(pointOut-east).mag()/wireLine.mag();
389 hitPosition=pointOut;
390 }
391 }
392
393 //Calculate signal transferT on wire
394 G4double halfLayerLength=mdcGeomSvc->Layer(layerId)->Length()/2.;
395 G4double halfWireLength=wireLine.mag()/2.;
396 G4double transferZ=0;
397 if(layerId%2==0){
398 transferZ=halfLayerLength+hitPosition.z(); //West readout
399 }else{
400 transferZ=halfLayerLength-hitPosition.z(); //East readout
401 }
402 if(layerId<8){
403 transferT=transferZ*halfWireLength/halfLayerLength/220;
404 }else{
405 transferT=transferZ*halfWireLength/halfLayerLength/240;
406 }
407
408 return distance;
409
410}
MdcGeoLayer::Length
double Length(void) const
Definition: MdcGeoLayer.h:161
MdcGeoWire::Forward
HepPoint3D Forward(void) const
Definition: MdcGeoWire.h:129
MdcGeoWire::Backward
HepPoint3D Backward(void) const
Definition: MdcGeoWire.h:128
MdcGeomSvc::Wire
const MdcGeoWire *const Wire(unsigned id)
Definition: MdcGeomSvc.cxx:768
MdcGeomSvc::Layer
const MdcGeoLayer *const Layer(unsigned id)
Definition: MdcGeomSvc.cxx:784

Referenced by ProcessHits().

◆ EndOfEvent()

void BesMdcSD::EndOfEvent ( G4HCofThisEvent *  )

Definition at line 338 of file BesMdcSD.cc.

339{
340 if (verboseLevel>0) {
341 hitsCollection->PrintAllHits();
342 /*
343 G4int NbHits = hitsCollection->entries();
344 G4cout << "\n-------->Hits Collection: in this event they are " << NbHits
345 << " hits in the MDC chambers: " << G4endl;
346 for (G4int i=0;i<NbHits;i++) (*hitsCollection)[i]->Print();
347 */
348 }
349}

◆ EndOfTruthEvent()

void BesMdcSD::EndOfTruthEvent ( const G4Event *  evt)
virtual

Reimplemented from BesSensitiveDetector.

Definition at line 114 of file BesMdcSD.cc.

115{ static G4int HLID=-1;
116 if(HLID<0)
117 {
118 HLID = G4SDManager::GetSDMpointer()->
119 GetCollectionID(collectionName[1]);
120 }
121 G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
122 HCE->AddHitsCollection(HLID,truthCollection);
123}

◆ Initialize()

void BesMdcSD::Initialize ( G4HCofThisEvent *  HCE)

Definition at line 89 of file BesMdcSD.cc.

90{
91 hitsCollection = new BesMdcHitsCollection
92 (SensitiveDetectorName,collectionName[0]);
93 static G4int HCID = -1;
94 if(HCID<0)
95 { HCID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); }
96 HCE->AddHitsCollection( HCID, hitsCollection );
97 G4int i,j;
98 for(i=0; i<43;i++){
99 for(j=0;j<288;j++){
100 hitPointer[i][j]=-1;
101 truthPointer[i][j]=-1;
102 }
103 }
104}

◆ ProcessHits()

G4bool BesMdcSD::ProcessHits ( G4Step *  aStep,
G4TouchableHistory *   
)

Definition at line 125 of file BesMdcSD.cc.

126{
127 G4Track* gTrack = aStep->GetTrack() ;
128
129 G4double globalT=gTrack->GetGlobalTime();//Time since the event in which the track belongs is created
130 if(std::isnan(globalT)){
131 G4cout<<"MdcSD:error, globalT is nan "<<G4endl;
132 return false;
133 }
134 if(globalT > 2000)return false; //MDC T window is 2 microsecond
135
136 //skip neutral tracks
137 G4int charge = gTrack->GetDefinition()->GetPDGCharge();
138 if (charge == 0) return false;
139
140 //skip no energy deposit tracks
141 G4double stepLength=aStep->GetStepLength();
142 if(stepLength==0){
143 // G4cout<<"step length equal 0!!"<<G4endl;
144 return false;
145 }
146
147 G4double edep = aStep->GetTotalEnergyDeposit() / stepLength;
148 if(edep==0.) return false;
149
150 // get position of the track at the beginning and at the end of step
151 G4StepPoint* prePoint = aStep->GetPreStepPoint() ;
152 G4StepPoint* postPoint = aStep->GetPostStepPoint() ;
153
154 //get position coordinate
155 G4ThreeVector pointIn = prePoint->GetPosition();
156 G4ThreeVector pointOut = postPoint->GetPosition();
157
158 // get physical volumes
159 const G4VTouchable *touchable = prePoint->GetTouchable();
160 G4VPhysicalVolume *volume = touchable->GetVolume(0);
161
162 G4double driftD = 0;
163 G4double driftT = 0;
164 G4double edepTemp = 0;
165 G4double lengthTemp = 0;
166 G4int cellId=0;
167 G4int layerId = touchable->GetVolume(1)->GetCopyNo();
168 if(volume->IsReplicated()){
169 cellId = touchable->GetReplicaNumber();
170 }else{
171 cellId=touchable->GetVolume(0)->GetCopyNo();
172 }
173 if(layerId==18&&(cellId==27||cellId==42))return false; // no sense wire
174
175 if(ReadBoostRoot::GetMdc() == 2) { //gdml
176 //layerId 0-35 -> CopyNo 0-35 in gdml
177 //layerId 36-42 -> CopyNo (36,37),(38,39),...(48,49)
178 if(layerId >= 36) layerId = 36 + (layerId-36)/2;
179 }
180
181 G4double halfLayerLength=mdcGeomSvc->Layer(layerId)->Length()/2.;
182 if(((fabs(pointIn.z())-halfLayerLength)>-7.)
183 &&((fabs(pointOut.z())-halfLayerLength)>-7.))return false;//Out sensitive area
184
185 G4int trackID= gTrack->GetTrackID(); //G4 track ID of current track.
186 G4int truthID, g4TrackID;
187 GetCurrentTrackIndex(truthID, g4TrackID); //ID of current primary track.
188
189 G4double theta=gTrack->GetMomentumDirection().theta();
190
191 G4ThreeVector hitPosition(0,0,0);
192 G4double transferT=0;
193 driftD = Distance(layerId,cellId,pointIn,pointOut,hitPosition,transferT);
194
195 G4double posPhi, wirePhi;
196 posPhi=hitPosition.phi();//from -pi to pi
197 if(posPhi<0)posPhi += 2*pi;
198 wirePhi=mdcGeoPointer->SignalWire(layerId,cellId).Phi(hitPosition.z());//from 0 to 2pi
199
200 G4int posFlag;
201 if(posPhi<=wirePhi){
202 posFlag = 0;
203 }else{
204 posFlag = 1;
205 }
206 // if x axis is between pos and wire, phi will has a jump of one of them.
207 if(posPhi < 1. && wirePhi > 5.)posFlag = 1;
208 if(posPhi > 5. && wirePhi < 1.)posFlag = 0;
209
210 G4ThreeVector hitLine=pointOut-pointIn;
211 G4double enterAngle=hitLine.phi()-wirePhi;
212 while(enterAngle<-pi/2.)enterAngle+=pi;
213 while(enterAngle>pi/2.)enterAngle-=pi;
214
215 if(m_G4Svc->GetMdcDedxFlag()==1){
216 G4double betaGamma=aStep->GetPreStepPoint()->GetBeta() * aStep->GetPreStepPoint()->GetGamma();
217 if(betaGamma<0.01)return false;//too low momentum
218 //if (betaGamma < 10.0) betaGamma = 10.0;
219
220 G4double eCount=dedxSample(betaGamma, charge, theta);
221 edep=eCount;
222 }
223
224 BesMdcHit* newHit = new BesMdcHit();
225 newHit->SetTrackID(truthID);
226 //newHit->SetTrkID(trackID);
227 newHit->SetLayerNo(layerId);
228 newHit->SetCellNo(cellId);
229 newHit->SetEdep(edep);
230 newHit->SetPos(hitPosition);
231 newHit->SetDriftD(driftD);
232 newHit->SetTheta(theta);
233 newHit->SetPosFlag(posFlag);
234 newHit->SetEnterAngle(enterAngle);
235
236 //Transfer hit pointer to BesMdcCalTransfer
237 mdcCalPointer->SetHitPointer(newHit);
238
239 driftT=mdcCalPointer->D2T(driftD);
240 globalT+=transferT;
241 driftT+=globalT;
242
243 newHit->SetDriftT (driftT);
244 newHit->SetGlobalT(globalT);
245
246 if (hitPointer[layerId][cellId] == -1) {
247 hitsCollection->insert(newHit);
248 G4int NbHits = hitsCollection->entries();
249 hitPointer[layerId][cellId]=NbHits-1;
250 }
251 else
252 {
253 G4int pointer=hitPointer[layerId][cellId];
254 if (g4TrackID == trackID) {
255 G4double preDriftT=(*hitsCollection)[pointer]->GetDriftT();
256 }
257
258 G4double preDriftT = (*hitsCollection)[pointer]->GetDriftT();
259 if (driftT < preDriftT) {
260 (*hitsCollection)[pointer]->SetTrackID(truthID);
261 //(*hitsCollection)[pointer]->SetTrkID(trackID);
262 (*hitsCollection)[pointer]->SetDriftD(driftD);
263 (*hitsCollection)[pointer]->SetDriftT(driftT);
264 (*hitsCollection)[pointer]->SetPos(hitPosition);
265 (*hitsCollection)[pointer]->SetGlobalT(globalT);
266 (*hitsCollection)[pointer]->SetTheta(theta);
267 (*hitsCollection)[pointer]->SetPosFlag(posFlag);
268 (*hitsCollection)[pointer]->SetEnterAngle(enterAngle);
269 }
270
271 delete newHit;
272 }
273
274 //for mc truth
275 if(truthCollection){
276 if(g4TrackID==trackID){ //This track is the primary track & will appear in MC truth
277 G4int pointer=truthPointer[layerId][cellId];
278 if(pointer==-1){
279 BesMdcHit* truthHit = new BesMdcHit();
280 truthHit->SetTrackID (truthID);
281 truthHit->SetLayerNo(layerId);
282 truthHit->SetCellNo(cellId);
283 truthHit->SetEdep (edep);
284 truthHit->SetPos (hitPosition);
285 truthHit->SetDriftD (driftD);
286 truthHit->SetDriftT (driftT);
287 truthHit->SetGlobalT(globalT);
288 truthHit->SetTheta(theta);
289 truthHit->SetPosFlag(posFlag);
290 truthHit->SetEnterAngle(enterAngle);
291
292 truthCollection->insert(truthHit);
293 G4int NbHits = truthCollection->entries();
294 truthPointer[layerId][cellId]=NbHits-1;
295 }
296 else {
297 if(truthID == (*truthCollection)[pointer]->GetTrackID()){
298 G4double preDriftT=(*truthCollection)[pointer]->GetDriftT();
299 if(driftT<preDriftT){
300 (*truthCollection)[pointer]->SetDriftD(driftD);
301 (*truthCollection)[pointer]->SetDriftT(driftT);
302 (*truthCollection)[pointer]->SetPos(hitPosition);
303 (*truthCollection)[pointer]->SetGlobalT(globalT);
304 (*truthCollection)[pointer]->SetTheta(theta);
305 (*truthCollection)[pointer]->SetPosFlag(posFlag);
306 (*truthCollection)[pointer]->SetEnterAngle(enterAngle);
307 }
308 edepTemp=(*truthCollection)[pointer]->GetEdep();
309 (*truthCollection)[pointer]->SetEdep(edepTemp+edep);
310 } else {
311 BesMdcHit* truthHit = new BesMdcHit();
312 truthHit->SetTrackID (truthID);
313 truthHit->SetLayerNo(layerId);
314 truthHit->SetCellNo(cellId);
315 truthHit->SetEdep(edep);
316 truthHit->SetPos(hitPosition);
317 truthHit->SetDriftD (driftD);
318 truthHit->SetDriftT (driftT);
319 truthHit->SetGlobalT(globalT);
320 truthHit->SetTheta(theta);
321 truthHit->SetPosFlag(posFlag);
322 truthHit->SetEnterAngle(enterAngle);
323
324 truthCollection->insert(truthHit);
325 G4int NbHits = truthCollection->entries();
326 truthPointer[layerId][cellId]=NbHits-1;
327 }
328 }
329 }
330 }
331
332 //newHit->Print();
333// newHit->Draw();
334
335 return true;
336}
BesMdcCalTransfer::D2T
double D2T(double driftDNew)
Definition: BesMdcCalTransfer.cc:109
BesMdcCalTransfer::SetHitPointer
void SetHitPointer(BesMdcHit *hit)
Definition: BesMdcCalTransfer.cc:25
BesMdcGeoParameter::SignalWire
BesMdcWire SignalWire(int, int)
Definition: BesMdcGeoParameter.cc:118
BesMdcHit::SetEdep
void SetEdep(G4double de)
Definition: BesMdcHit.hh:41
BesMdcHit::SetDriftT
void SetDriftT(G4double time)
Definition: BesMdcHit.hh:44
BesMdcHit::SetEnterAngle
void SetEnterAngle(G4double angle)
Definition: BesMdcHit.hh:47
BesMdcHit::SetCellNo
void SetCellNo(G4int cell)
Definition: BesMdcHit.hh:40
BesMdcHit::SetPos
void SetPos(G4ThreeVector xyz)
Definition: BesMdcHit.hh:42
BesMdcHit::SetTrackID
void SetTrackID(G4int track)
Definition: BesMdcHit.hh:38
BesMdcHit::SetLayerNo
void SetLayerNo(G4int layer)
Definition: BesMdcHit.hh:39
BesMdcHit::SetTheta
void SetTheta(G4double angle)
Definition: BesMdcHit.hh:46
BesMdcHit::SetDriftD
void SetDriftD(G4double distance)
Definition: BesMdcHit.hh:43
BesMdcHit::SetGlobalT
void SetGlobalT(G4double time)
Definition: BesMdcHit.hh:45
BesMdcHit::SetPosFlag
void SetPosFlag(G4int flag)
Definition: BesMdcHit.hh:48
BesMdcSD::Distance
G4double Distance(G4int, G4int, G4ThreeVector, G4ThreeVector, G4ThreeVector &, G4double &)
Definition: BesMdcSD.cc:351
BesMdcWire::Phi
double Phi(void) const
Definition: BesMdcGeoParameter.hh:33
BesSensitiveDetector::GetCurrentTrackIndex
void GetCurrentTrackIndex(G4int &trackIndex, G4int &g4TrackId) const
Definition: BesSensitiveDetector.cc:30
ReadBoostRoot::GetMdc
static G4int GetMdc()
Definition: ReadBoostRoot.hh:13
charge
float charge
Definition: prepare_data.cxx:28
pi
const float pi
Definition: vector3.h:133
The documentation for this class was generated from the following files: