CGEM BOSS 6.6.5.h
BESIII Offline Software System
Loading...
Searching...
No Matches
BesMdcDigitizer.cc
Go to the documentation of this file.
1//---------------------------------------------------------------------------//
2// BOOST --- BESIII Object_Oriented Simulation Tool //
3//---------------------------------------------------------------------------//
4//Description:
5//Author: Yuan Ye([email protected])
6//Created: Oct. 26, 2004
7//Modified:
8//Comment:
9//---------------------------------------------------------------------------//
10
11#include "BesMdcDigitizer.hh"
12#include "BesMdcDigitizerMessenger.hh"
13#include "BesMdcHit.hh"
14
15#include "G4DigiManager.hh"
16#include "Randomize.hh"
17#include "G4ios.hh"
18#include <string>
19
20#include "TFile.h"
21#include "TH1F.h"
22
23#include "G4Svc/IG4Svc.h"
24#include "G4Svc/G4Svc.h"
25#include "GaudiKernel/ISvcLocator.h"
26#include "GaudiKernel/Bootstrap.h"
27#include "GaudiKernel/IDataProviderSvc.h"
28
29BesMdcDigitizer::BesMdcDigitizer(G4String modName):G4VDigitizerModule(modName){
30 noiseFlag=0;
31 noiseType=3;
32 noiseLevel=0.1;//10%
33 maxNoiseT=300.;//ns
34 smearFlag=1;
35 mdcDRes = 0.13; //mm
36 effFlag = 0;
37 for(G4int i=0; i<43;i++){
38 layerEff.push_back(1.);
39 }
40 collectionName.push_back("BesMdcDigisCollection");
41 digitizerMessenger = new BesMdcDigitizerMessenger(this);
42 mdcGeoPointer=BesMdcGeoParameter::GetGeo();
43 mdcCalPointer=new BesMdcCalTransfer;
44
45 // ISvcLocator* svcLocator = Gaudi::svcLocator();
46 IG4Svc* tmpSvc;
47 //G4Svc* m_G4Svc;
48 StatusCode sc=Gaudi::svcLocator()->service("G4Svc", tmpSvc);
49 if (!sc.isSuccess())
50 G4cout <<" MdcDigitizer::Error,could not open G4Svc"<<G4endl;
51
52 m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc);
53
54 IMdcTunningSvc* IMdcTunningSvc;
55 sc= Gaudi::svcLocator()->service("MdcTunningSvc",IMdcTunningSvc);
56 if (!sc.isSuccess()){
57 G4cout <<" MdcDigitizer::Error,could not open Mdc Tunning Service"<<G4endl;
58 }else{
59 G4cout<<" MdcDigitizer:: Open Mdc Tunning Service"<<G4endl;
60 }
61 mdcTunningSvc=dynamic_cast<MdcTunningSvc *>(IMdcTunningSvc);
62
63 std::string noiseFile=m_G4Svc->GetMdcNoiseFile();
64 f=new TFile(noiseFile.c_str());
65 h1=(TH1F*)f->Get("h703");
66 h2=(TH1F*)f->Get("h501");
67 h3=(TH1F*)f->Get("h801");
68 /*
69 //get Mdc Ntuple from G4Svc
70 if(m_G4Svc->MdcRootFlag())
71 {
72 m_tupleMdc = m_G4Svc->GetTupleMdc();
73 sc = m_tupleMdc->addItem("NHits",m_NHits);
74 sc = m_tupleMdc->addItem("LayerId",m_layerId);
75 sc = m_tupleMdc->addItem("cellId",m_cellId);
76 sc = m_tupleMdc->addItem("Edep",m_edep);
77 sc = m_tupleMdc->addItem("driftD",m_driftD);
78 // sc = m_tupleMdc->addItem("driftT",m_driftT);
79 sc = m_tupleMdc->addItem("globalT",m_globalT);
80 sc = m_tupleMdc->addItem("theta",m_theta);
81 sc = m_tupleMdc->addItem("enterAngle",m_enterAngle);
82 sc = m_tupleMdc->addItem("driftDNew",m_driftDNew);
83 sc = m_tupleMdc->addItem("driftTNew",m_driftTNew);
84 // sc = m_tupleMdc->addItem("adc",m_adc);
85 // sc = m_tupleMdc->addItem("tdc",m_tdc);
86 }
87 */
88}
89
90BesMdcDigitizer::~BesMdcDigitizer(){delete digitizerMessenger;}
91
92void BesMdcDigitizer::SetEff(G4int layer, G4double eff){
93 if(layer==-1){
94 for(G4int i=0; i<43;i++){
95 layerEff[i]=eff;
96 }
97 }else{
98 layerEff[layer]=eff;
99 }
100}
101
102void BesMdcDigitizer::Digitize(){
103
104 //initialize
105 for(G4int i=0; i<43;i++){
106 for(G4int j=0;j<288;j++){
107 digiPointer[i][j]=-1;
108 }
109 }
110
111 G4int NHits,layerId, cellId, posFlag;
112 G4double edep,driftD,driftT, globalT, theta,cosTheta,enterAngle;
113 G4double mean,sigma,mean1,mean2,sigma1, sigma2, f,sig,delSig, fRandom, driftDNew, driftTNew;
114 G4double tempEff;
115 G4double resLargest,resSmallest,resRatio;//added by liukai
116
117 G4DigiManager* DigiMan = G4DigiManager::GetDMpointer();
118
119 //hits collection ID
120 G4int THCID=-1;
121 THCID = DigiMan->GetHitsCollectionID("BesMdcHitsCollection");
122
123 //hits collection
124 BesMdcHitsCollection* THC = 0;
125 THC = (BesMdcHitsCollection*) (DigiMan->GetHitsCollection(THCID));
126
127 if(THC){
128 digisCollection=new BesMdcDigisCollection
129 (moduleName, collectionName[0]);
130 NHits=THC->entries();
131 for(G4int i=0;i<NHits;i++){
132 layerId = (*THC)[i]->GetLayerNo();
133 cellId = (*THC)[i]->GetCellNo();
134 edep = (*THC)[i]->GetEdep();
135 driftD = (*THC)[i]->GetDriftD();
136 globalT = (*THC)[i]->GetGlobalT();
137 theta = (*THC)[i]->GetTheta();
138 cosTheta = cos(theta);
139 enterAngle = (*THC)[i]->GetEnterAngle();
140 posFlag = (*THC)[i]->GetPosFlag();
141
142 //Transfer hit pointer to BesMdcCalTransfer
143 mdcCalPointer->SetHitPointer((*THC)[i]);
144
145 //Filter with wire efficiency
146 if(effFlag==0){
147 //tempEff = mdcCalPointer->GetEff();
148 tempEff=mdcTunningSvc->GetEff(layerId,cellId,driftD,cosTheta,posFlag);
149 }else{
150 tempEff = layerEff[layerId];
151 }
152 fRandom=G4UniformRand();
153 if(fRandom>tempEff)continue;
154
155 //cout<<"layerid "<<layerId<<" cellid "<<cellId<<" theta "<<cosTheta<<" enterangle "<<enterAngle<<endl;
156 //Drift distance smear
157 if(smearFlag==0){ //No smear
158 driftDNew = driftD;
159 }else if(smearFlag==1){ //Smear from TuningSvc
160 // mdcTunningSvc->GetRes(layerId,cellId,driftD,cosTheta,posFlag,enterAngle,mean,sigma);
161 //mdcTunningSvc->GetRes2(layerId,cellId,driftD,cosTheta,posFlag,enterAngle,f,mean1,sigma1,mean2,sigma2);
162 mdcTunningSvc->GetRes3(layerId,cellId,driftD,cosTheta,posFlag,enterAngle,f,mean1,sigma1,mean2,sigma2,resLargest,resSmallest,resRatio);
163
164 //driftDNew = Smear(driftD,f,mean1,sigma1,mean2,sigma2);
165 //driftDNew = Smear(driftD-(f*mean1+(1-f)*mean2),f,mean1,sigma1,mean2,sigma2);//new method
166
167 driftDNew = Smear(driftD-(f*mean1+(1-f)*mean2),f,mean1,sigma1,mean2,sigma2,resLargest,resSmallest,resRatio);//----added by liukai 2012-6-4
168
169
170 }else if(smearFlag==2){ //Smear with fixed resolution
171 driftDNew = Smear(driftD);
172 }else{
173 G4cerr<<"MdcDigitizer::worong smearFlag: "<<smearFlag<<G4endl;
174 }
175
176 //Do X-T conversion
177 driftTNew = mdcCalPointer->D2T(driftDNew);
178
179 //Do Q-T correct
180 driftTNew += mdcCalPointer->GetTimeWalk();
181
182 //Add T0
183 driftTNew += mdcCalPointer->GetT0();
184
185 //Add TOF
186 driftTNew += globalT;
187
188 //Signal transfer time on wire
189 // transferT=Transfer(layerId,cellId,hitPosition);
190 //driftTNew+=transferT;
191
192 if(isnan(driftTNew)){
193 G4cout<<"MdcDigitizer::error, driftT is nan"<<G4endl;
194 continue;
195 }
196
197 /*
198 if(m_G4Svc->MdcRootFlag())
199 {
200 m_NHits= NHits;
201 m_layerId= layerId;
202 m_cellId= cellId;
203 m_edep= edep;
204 m_driftD= driftD;
205 // m_driftT= driftT;
206 m_globalT = globalT;
207 m_enterAngle = enterAngle;
208 m_driftDNew = driftDNew;
209 m_driftTNew = driftTNew;
210 m_theta = theta;
211 m_tupleMdc ->write();
212 }
213 */
214 BesMdcDigi* newDigi = new BesMdcDigi();
215 newDigi->SetTrackID((*THC)[i]->GetTrackID());
216 newDigi->SetLayerNo(layerId);
217 newDigi->SetCellNo(cellId);
218 newDigi->SetEdep(edep);
219 newDigi->SetDriftT(driftTNew);
220 G4int NbDigis = digisCollection->insert(newDigi);
221 digiPointer[layerId][cellId]=NbDigis-1;
222 (*THC)[i]->SetDigiId(NbDigis-1);
223 }
224
225 if(noiseFlag==1)AddNoise();
226 if(noiseFlag==2){
227 ifstream readNoiseLevel("$MDCSIMROOT/share/noiselevel.txt");
228 if(!readNoiseLevel.good()){
229 std::cout<<" Error , noiselevel file not exist "<<std::endl;
230 }else{
231 std::cout<<" MdcDigitizer:: Open noiselevel file "<<std::endl;
232 }
233 G4int NLayer=mdcGeoPointer->SignalLayerNo();
234 G4double level;
235 for(G4int i=0;i<NLayer;i++){
236 readNoiseLevel>>level;
237 mixLevel.push_back(level);
238 }
239 AddNoise2();
240 }
241
242 if (verboseLevel>0) {
243 G4cout << "\n-------->digis Collection: in this event they are "
244 << digisCollection->entries()
245 << " digis in the MDC chambers: " << G4endl;
246 digisCollection->PrintAllDigi();
247 }
248 StoreDigiCollection(digisCollection);
249 }
250
251}
252
253G4double BesMdcDigitizer::Smear(G4double driftD){
254 G4double r, driftDNew;
255 r = G4RandGauss::shoot();
256 driftDNew = driftD + mdcDRes * r;
257 while(driftDNew<=0){
258 r = G4RandGauss::shoot();
259 driftDNew = driftD + mdcDRes * r;
260 }
261 return driftDNew;
262}
263
264G4double BesMdcDigitizer::Smear(G4double driftD,G4double sigma,G4double mean){
265 G4double r, driftDNew;
266 r = G4RandGauss::shoot();
267 driftDNew = driftD + sigma * r+mean;
268 while(driftDNew <= 0){
269 r = G4RandGauss::shoot();
270 driftDNew = driftD + sigma * r+mean;
271 }
272 return driftDNew;
273}
274
275G4double BesMdcDigitizer::Smear(G4double driftD,G4double f,G4double mean1,G4double sigma1,G4double mean2,G4double sigma2){
276 G4double r, driftDNew,mean,sigma;
277 r = G4UniformRand();
278 if(r<=f){
279 mean=mean1;
280 sigma=sigma1;
281 }else{
282 mean=mean2;
283 sigma=sigma2;
284 }
285 int times=0;
286 r = G4RandGauss::shoot();
287 driftDNew = driftD + sigma * r+mean;
288 while(driftDNew <= 0){
289 r = G4RandGauss::shoot();
290 driftDNew = driftD + sigma * r+mean;
291 times++;
292 if(times>10)driftDNew=0.01;
293 }
294 return driftDNew;
295}
296
297G4double BesMdcDigitizer::Smear(G4double driftD,G4double f,G4double mean1,G4double sigma1,G4double mean2,G4double sigma2,G4double resLargest,G4double resSmallest,G4double resRatio){
298 G4double r, driftDNew,mean,sigma;
299 G4double ratio,tempd;
300 ratio=G4UniformRand();
301 int times;
302 if(ratio<=resRatio)
303 {
304 //for hitOnTrk distribution
305 r = G4UniformRand();
306 if(r<=f){
307 mean=mean1;
308 sigma=sigma1;
309 }else{
310 mean=mean2;
311 sigma=sigma2;
312 }
313 times=0;
314 r = G4RandGauss::shoot();
315 driftDNew = driftD + sigma * r+mean;
316 }
317 else//for hitNotOnTrk
318 {
319 tempd=G4UniformRand()*2.0+resLargest;
320 times=0;
321 driftDNew = driftD + tempd;
322 }
323 while(driftDNew <= 0){
324 r = G4RandGauss::shoot();
325 driftDNew = driftD + sigma * r+mean;
326 times++;
327 if(times>10)driftDNew=0.01;
328 }
329 return driftDNew;
330}
331
332void BesMdcDigitizer::AddNoise2(void){
333 G4int wireNo;
334 G4double random;
335 G4double randomT;
336 G4double randomQ;
337 G4int NLayer=mdcGeoPointer->SignalLayerNo();
338 for(G4int i=0;i<NLayer;i++){
339 wireNo=mdcGeoPointer->SignalLayer(i).WireNo()/2;
340 for(G4int j=0;j<wireNo;j++){
341 random=G4UniformRand();
342 if(random < mixLevel[i]){
343 randomT=G4UniformRand() * 2000;
344 if(isnan(randomT)){
345 G4cout<<"MdcDigitizer: error, randomT is nan"<<G4endl;
346 continue;
347 }
348
349 randomQ=200.;
350 if(digiPointer[i][j]!=-1){
351 G4int pointer=digiPointer[i][j];
352 G4double preDriftT=(*digisCollection)[pointer]->GetDriftT();
353 if(preDriftT <= randomT)continue;
354 (*digisCollection)[pointer]->SetDriftT(randomT);
355 (*digisCollection)[pointer]->SetTrackID(-1);
356 }else{
357 BesMdcDigi* newDigi = new BesMdcDigi();
358 newDigi->SetTrackID(-1);
359 newDigi->SetLayerNo(i);
360 newDigi->SetCellNo(j);
361 newDigi->SetEdep(randomQ);
362 newDigi->SetDriftT(randomT);
363 digisCollection->insert(newDigi);
364 }
365 }
366 }
367 }
368}
369
370
371void BesMdcDigitizer::AddNoise(void){
372 G4double r0,r;
373 vector<G4double> noise; //Noise level of each layer
374
375 G4int NLayer=mdcGeoPointer->SignalLayerNo();
376 if(noiseType==0){
377 for(G4int i=0;i<NLayer;i++){
378 noise.push_back(noiseLevel);
379 }
380 }else if(noiseType==1){
381 r0=mdcGeoPointer->SignalLayer(0).R();
382 for(G4int i=0;i<NLayer;i++){
383 r=mdcGeoPointer->SignalLayer(i).R();
384 noise.push_back(noiseLevel * r0 / r);
385 }
386 }else if(noiseType==2){
387 r0=mdcGeoPointer->SignalLayer(0).R();
388 for(G4int i=0;i<NLayer;i++){
389 r=mdcGeoPointer->SignalLayer(i).R();
390 noise.push_back(noiseLevel * r0 * r0 / r / r);
391 }
392 }else if(noiseType==3){ // liugc add 22:11 4/14/06
393 Int_t Nbins=(Int_t)h1->GetNbinsX();
394
395 double xmax=h1->GetXaxis()->GetXmax();
396 double xmin=h1->GetXaxis()->GetXmin();
397 double dx=(xmax-xmin)/Nbins;
398 G4double y;
399 for(G4int i=0;i<Nbins;i++){
400 double x=double(i+1)*dx;
401 y=(G4double)h1->GetBinContent(x);
402 y=y*noiseLevel/0.05608559;// Normalize use noise level of 1st layer in "run23096noise.root"
403 noise.push_back(y);
404 }
405 }
406
407 G4int wireNo;
408 G4double random;
409 G4double randomT;
410 G4double randomQ;
411 G4double T0=m_G4Svc->GetBeamStartTime();
412 for(G4int i=0;i<NLayer;i++){
413 wireNo=mdcGeoPointer->SignalLayer(i).WireNo()/2;
414 for(G4int j=0;j<wireNo;j++){
415 random=G4UniformRand();
416 if(random < noise[i]){
417 //randomT=G4UniformRand() * maxNoiseT;
418 randomT=h2->GetRandom()+T0;
419 // randomT=randomT;
420 // randomQ=h3->GetRandom();
421 // randomQ=randomQ*0.001*0.001; //Transfer from TDC channels to Mev, coef. is temporary one.
422 if(isnan(randomT)){
423 G4cout<<"MdcDigitizer: error, randomT is nan"<<G4endl;
424 continue;
425 }
426
427 randomQ=0.;
428 if(digiPointer[i][j]!=-1){
429 G4int pointer=digiPointer[i][j];
430 G4double signalEdep=(*digisCollection)[pointer]->GetEdep();
431 (*digisCollection)[pointer]->SetEdep(randomQ+signalEdep);
432 G4double preDriftT=(*digisCollection)[pointer]->GetDriftT();
433 if(preDriftT <= randomT)continue;
434 (*digisCollection)[pointer]->SetDriftT(randomT);
435 (*digisCollection)[pointer]->SetTrackID(-1);
436 }else{
437 BesMdcDigi* newDigi = new BesMdcDigi();
438 newDigi->SetTrackID(-1);
439 newDigi->SetLayerNo(i);
440 newDigi->SetCellNo(j);
441 newDigi->SetEdep(randomQ);
442 newDigi->SetDriftT(randomT);
443 digisCollection->insert(newDigi);
444 }
445 }
446 }
447 }
448}
cos
double cos(const BesAngle a)
Definition BesAngle.h:213
BesMdcDigisCollection
G4TDigiCollection< BesMdcDigi > BesMdcDigisCollection
Definition BesMdcDigi.hh:59
BesMdcHitsCollection
G4THitsCollection< BesMdcHit > BesMdcHitsCollection
Definition BesMdcHit.hh:97
sigma2
double sigma2(0)
x
Double_t x[10]
Definition DataBase/tau_mode.c:57
wireNo
const int wireNo
Definition DedxCalibWireGain.cxx:18
BesMdcCalTransfer::D2T
double D2T(double driftDNew)
Definition BesMdcCalTransfer.cc:109
BesMdcCalTransfer::SetHitPointer
void SetHitPointer(BesMdcHit *hit)
Definition BesMdcCalTransfer.cc:25
BesMdcDigi::SetEdep
void SetEdep(G4double de)
Definition BesMdcDigi.hh:40
BesMdcDigi::SetCellNo
void SetCellNo(G4int cell)
Definition BesMdcDigi.hh:39
BesMdcDigi::SetDriftT
void SetDriftT(G4double time)
Definition BesMdcDigi.hh:41
BesMdcDigi::SetTrackID
void SetTrackID(G4int track)
Definition BesMdcDigi.hh:37
BesMdcDigi::SetLayerNo
void SetLayerNo(G4int layer)
Definition BesMdcDigi.hh:38
BesMdcDigitizer::Digitize
virtual void Digitize()
Definition BesMdcDigitizer.cc:102
BesMdcDigitizer::BesMdcDigitizer
BesMdcDigitizer(G4String modName)
Definition BesMdcDigitizer.cc:29
BesMdcDigitizer::SetEff
void SetEff(G4int layer, G4double eff)
Definition BesMdcDigitizer.cc:92
BesMdcGeoParameter::GetGeo
static BesMdcGeoParameter * GetGeo(void)
Definition BesMdcGeoParameter.cc:34
BesMdcGeoParameter::SignalLayer
const BesMdcLayer & SignalLayer(int) const
Definition BesMdcGeoParameter.cc:140
BesMdcLayer::WireNo
int WireNo(void) const
Definition BesMdcGeoParameter.hh:68
BesMdcWire::R
double R(void) const
Definition BesMdcGeoParameter.hh:35
G4Svc
Definition G4Svc.h:32
G4Svc::GetMdcNoiseFile
std::string GetMdcNoiseFile()
Definition G4Svc.h:87
G4Svc::GetBeamStartTime
double GetBeamStartTime()
Definition G4Svc.h:80
MdcTunningSvc::GetEff
double GetEff(int layerId, int cellId, double driftD, double cosTheta, int posFlag)
Definition MdcTunningSvc.cc:495
MdcTunningSvc::GetRes3
double GetRes3(int layerId, int cellId, double driftD, double cosTheta, int posFlag, double entranceAngle, double &f, double &mean1, double &sigma1, double &mean2, double &sigma2, double &ResLargest, double &ResSmallest, double &ResRatio)
Definition MdcTunningSvc.cc:640