#include <BesMdcDigitizer.hh>
Public Member Functions | |
| BesMdcDigitizer (G4String modName) | |
| BesMdcDigitizer (G4String modName) | |
| virtual void | Digitize () |
| virtual void | Digitize () |
| void | SetEff (G4int layer, G4double eff) |
| void | SetEff (G4int layer, G4double eff) |
| void | SetEffFlag (G4int flag) |
| void | SetEffFlag (G4int flag) |
| void | SetMdcDRes (G4double res) |
| void | SetMdcDRes (G4double res) |
| void | SetNoiseFlag (G4int flag) |
| void | SetNoiseFlag (G4int flag) |
| void | SetNoiseLevel (G4double level) |
| void | SetNoiseLevel (G4double level) |
| void | SetNoiseType (G4int type) |
| void | SetNoiseType (G4int type) |
| void | SetSmearFlag (G4int flag) |
| void | SetSmearFlag (G4int flag) |
| ~BesMdcDigitizer () | |
| ~BesMdcDigitizer () | |
Private Member Functions | |
| void | AddNoise (void) |
| void | AddNoise (void) |
| void | AddNoise2 (void) |
| void | AddNoise2 (void) |
| G4double | Smear (G4double, G4double, G4double, G4double, G4double, G4double) |
| G4double | Smear (G4double, G4double, G4double) |
| G4double | Smear (G4double) |
| G4double | Smear (G4double, G4double, G4double, G4double, G4double, G4double) |
| G4double | Smear (G4double, G4double, G4double) |
| G4double | Smear (G4double) |
Private Attributes | |
| G4int | digiPointer [43][288] |
| BesMdcDigisCollection * | digisCollection |
| BesMdcDigisCollection * | digisCollection |
| BesMdcDigitizerMessenger * | digitizerMessenger |
| BesMdcDigitizerMessenger * | digitizerMessenger |
| G4int | effFlag |
| TFile * | f |
| TFile * | f |
| TH1F * | h1 |
| TH1F * | h1 |
| TH1F * | h2 |
| TH1F * | h2 |
| TH1F * | h3 |
| TH1F * | h3 |
| vector< G4double > | layerEff |
| vector< G4double > | layerEff |
| NTuple::Item< long > | m_cellId |
| NTuple::Item< long > | m_cellId |
| NTuple::Item< double > | m_driftD |
| NTuple::Item< double > | m_driftD |
| NTuple::Item< double > | m_driftDNew |
| NTuple::Item< double > | m_driftDNew |
| NTuple::Item< double > | m_driftTNew |
| NTuple::Item< double > | m_driftTNew |
| NTuple::Item< double > | m_edep |
| NTuple::Item< double > | m_edep |
| NTuple::Item< double > | m_enterAngle |
| NTuple::Item< double > | m_enterAngle |
| G4Svc * | m_G4Svc |
| G4Svc * | m_G4Svc |
| NTuple::Item< double > | m_globalT |
| NTuple::Item< double > | m_globalT |
| NTuple::Item< long > | m_layerId |
| NTuple::Item< long > | m_layerId |
| NTuple::Item< long > | m_NHits |
| NTuple::Item< long > | m_NHits |
| NTuple::Item< double > | m_theta |
| NTuple::Item< double > | m_theta |
| NTuple::Tuple * | m_tupleMdc |
| NTuple::Tuple * | m_tupleMdc |
| G4double | maxNoiseT |
| BesMdcCalTransfer * | mdcCalPointer |
| BesMdcCalTransfer * | mdcCalPointer |
| G4double | mdcDRes |
| BesMdcGeoParameter * | mdcGeoPointer |
| BesMdcGeoParameter * | mdcGeoPointer |
| MdcTunningSvc * | mdcTunningSvc |
| MdcTunningSvc * | mdcTunningSvc |
| vector< G4double > | mixLevel |
| vector< G4double > | mixLevel |
| G4int | noiseFlag |
| G4double | noiseLevel |
| G4int | noiseType |
| G4int | smearFlag |
|
|
00029 :G4VDigitizerModule(modName){
00030 noiseFlag=0;
00031 noiseType=3;
00032 noiseLevel=0.1;//10%
00033 maxNoiseT=300.;//ns
00034 smearFlag=1;
00035 mdcDRes = 0.13; //mm
00036 effFlag = 0;
00037 for(G4int i=0; i<43;i++){
00038 layerEff.push_back(1.);
00039 }
00040 collectionName.push_back("BesMdcDigisCollection");
00041 digitizerMessenger = new BesMdcDigitizerMessenger(this);
00042 mdcGeoPointer=BesMdcGeoParameter::GetGeo();
00043 mdcCalPointer=new BesMdcCalTransfer;
00044
00045 // ISvcLocator* svcLocator = Gaudi::svcLocator();
00046 IG4Svc* tmpSvc;
00047 //G4Svc* m_G4Svc;
00048 StatusCode sc=Gaudi::svcLocator()->service("G4Svc", tmpSvc);
00049 if (!sc.isSuccess())
00050 G4cout <<" MdcDigitizer::Error,could not open G4Svc"<<G4endl;
00051
00052 m_G4Svc=dynamic_cast<G4Svc *>(tmpSvc);
00053
00054 IMdcTunningSvc* IMdcTunningSvc;
00055 sc= Gaudi::svcLocator()->service("MdcTunningSvc",IMdcTunningSvc);
00056 if (!sc.isSuccess()){
00057 G4cout <<" MdcDigitizer::Error,could not open Mdc Tunning Service"<<G4endl;
00058 }else{
00059 G4cout<<" MdcDigitizer:: Open Mdc Tunning Service"<<G4endl;
00060 }
00061 mdcTunningSvc=dynamic_cast<MdcTunningSvc *>(IMdcTunningSvc);
00062
00063 std::string noiseFile=m_G4Svc->GetMdcNoiseFile();
00064 f=new TFile(noiseFile.c_str());
00065 h1=(TH1F*)f->Get("h703");
00066 h2=(TH1F*)f->Get("h501");
00067 h3=(TH1F*)f->Get("h801");
00068 /*
00069 //get Mdc Ntuple from G4Svc
00070 if(m_G4Svc->MdcRootFlag())
00071 {
00072 m_tupleMdc = m_G4Svc->GetTupleMdc();
00073 sc = m_tupleMdc->addItem("NHits",m_NHits);
00074 sc = m_tupleMdc->addItem("LayerId",m_layerId);
00075 sc = m_tupleMdc->addItem("cellId",m_cellId);
00076 sc = m_tupleMdc->addItem("Edep",m_edep);
00077 sc = m_tupleMdc->addItem("driftD",m_driftD);
00078 // sc = m_tupleMdc->addItem("driftT",m_driftT);
00079 sc = m_tupleMdc->addItem("globalT",m_globalT);
00080 sc = m_tupleMdc->addItem("theta",m_theta);
00081 sc = m_tupleMdc->addItem("enterAngle",m_enterAngle);
00082 sc = m_tupleMdc->addItem("driftDNew",m_driftDNew);
00083 sc = m_tupleMdc->addItem("driftTNew",m_driftTNew);
00084 // sc = m_tupleMdc->addItem("adc",m_adc);
00085 // sc = m_tupleMdc->addItem("tdc",m_tdc);
00086 }
00087 */
00088 }
|
|
|
00090 {delete digitizerMessenger;}
|
|
|
|
|
|
|
|
|
|
|
|
00331 {
00332 G4double r0,r;
00333 vector<G4double> noise; //Noise level of each layer
00334
00335 G4int NLayer=mdcGeoPointer->SignalLayerNo();
00336 if(noiseType==0){
00337 for(G4int i=0;i<NLayer;i++){
00338 noise.push_back(noiseLevel);
00339 }
00340 }else if(noiseType==1){
00341 r0=mdcGeoPointer->SignalLayer(0).R();
00342 for(G4int i=0;i<NLayer;i++){
00343 r=mdcGeoPointer->SignalLayer(i).R();
00344 noise.push_back(noiseLevel * r0 / r);
00345 }
00346 }else if(noiseType==2){
00347 r0=mdcGeoPointer->SignalLayer(0).R();
00348 for(G4int i=0;i<NLayer;i++){
00349 r=mdcGeoPointer->SignalLayer(i).R();
00350 noise.push_back(noiseLevel * r0 * r0 / r / r);
00351 }
00352 }else if(noiseType==3){ // liugc add 22:11 4/14/06
00353 Int_t Nbins=(Int_t)h1->GetNbinsX();
00354
00355 double xmax=h1->GetXaxis()->GetXmax();
00356 double xmin=h1->GetXaxis()->GetXmin();
00357 double dx=(xmax-xmin)/Nbins;
00358 G4double y;
00359 for(G4int i=0;i<Nbins;i++){
00360 double x=double(i+1)*dx;
00361 y=(G4double)h1->GetBinContent(x);
00362 y=y*noiseLevel/0.05608559;// Normalize use noise level of 1st layer in "run23096noise.root"
00363 noise.push_back(y);
00364 }
00365 }
00366
00367 G4int wireNo;
00368 G4double random;
00369 G4double randomT;
00370 G4double randomQ;
00371 G4double T0=m_G4Svc->GetBeamStartTime();
00372 for(G4int i=0;i<NLayer;i++){
00373 wireNo=mdcGeoPointer->SignalLayer(i).WireNo()/2;
00374 for(G4int j=0;j<wireNo;j++){
00375 random=G4UniformRand();
00376 if(random < noise[i]){
00377 //randomT=G4UniformRand() * maxNoiseT;
00378 randomT=h2->GetRandom()+T0;
00379 // randomT=randomT;
00380 // randomQ=h3->GetRandom();
00381 // randomQ=randomQ*0.001*0.001; //Transfer from TDC channels to Mev, coef. is temporary one.
00382 if(isnan(randomT)){
00383 G4cout<<"MdcDigitizer: error, randomT is nan"<<G4endl;
00384 continue;
00385 }
00386
00387 randomQ=0.;
00388 if(digiPointer[i][j]!=-1){
00389 G4int pointer=digiPointer[i][j];
00390 G4double signalEdep=(*digisCollection)[pointer]->GetEdep();
00391 (*digisCollection)[pointer]->SetEdep(randomQ+signalEdep);
00392 G4double preDriftT=(*digisCollection)[pointer]->GetDriftT();
00393 if(preDriftT <= randomT)continue;
00394 (*digisCollection)[pointer]->SetDriftT(randomT);
00395 (*digisCollection)[pointer]->SetTrackID(-1);
00396 }else{
00397 BesMdcDigi* newDigi = new BesMdcDigi();
00398 newDigi->SetTrackID(-1);
00399 newDigi->SetLayerNo(i);
00400 newDigi->SetCellNo(j);
00401 newDigi->SetEdep(randomQ);
00402 newDigi->SetDriftT(randomT);
00403 digisCollection->insert(newDigi);
00404 }
00405 }
00406 }
00407 }
00408 }
|
|
|
|
|
|
00292 {
00293 G4int wireNo;
00294 G4double random;
00295 G4double randomT;
00296 G4double randomQ;
00297 G4int NLayer=mdcGeoPointer->SignalLayerNo();
00298 for(G4int i=0;i<NLayer;i++){
00299 wireNo=mdcGeoPointer->SignalLayer(i).WireNo()/2;
00300 for(G4int j=0;j<wireNo;j++){
00301 random=G4UniformRand();
00302 if(random < mixLevel[i]){
00303 randomT=G4UniformRand() * 2000;
00304 if(isnan(randomT)){
00305 G4cout<<"MdcDigitizer: error, randomT is nan"<<G4endl;
00306 continue;
00307 }
00308
00309 randomQ=200.;
00310 if(digiPointer[i][j]!=-1){
00311 G4int pointer=digiPointer[i][j];
00312 G4double preDriftT=(*digisCollection)[pointer]->GetDriftT();
00313 if(preDriftT <= randomT)continue;
00314 (*digisCollection)[pointer]->SetDriftT(randomT);
00315 (*digisCollection)[pointer]->SetTrackID(-1);
00316 }else{
00317 BesMdcDigi* newDigi = new BesMdcDigi();
00318 newDigi->SetTrackID(-1);
00319 newDigi->SetLayerNo(i);
00320 newDigi->SetCellNo(j);
00321 newDigi->SetEdep(randomQ);
00322 newDigi->SetDriftT(randomT);
00323 digisCollection->insert(newDigi);
00324 }
00325 }
00326 }
00327 }
00328 }
|
|
|
|
|
|
00102 {
00103
00104 //initialize
00105 for(G4int i=0; i<43;i++){
00106 for(G4int j=0;j<288;j++){
00107 digiPointer[i][j]=-1;
00108 }
00109 }
00110
00111 G4int NHits,layerId, cellId, posFlag;
00112 G4double edep,driftD,driftT, globalT, theta,cosTheta,enterAngle;
00113 G4double mean,sigma,mean1,mean2,sigma1, sigma2, f,sig,delSig, fRandom, driftDNew, driftTNew;
00114 G4double tempEff;
00115
00116 G4DigiManager* DigiMan = G4DigiManager::GetDMpointer();
00117
00118 //hits collection ID
00119 G4int THCID=-1;
00120 THCID = DigiMan->GetHitsCollectionID("BesMdcHitsCollection");
00121
00122 //hits collection
00123 BesMdcHitsCollection* THC = 0;
00124 THC = (BesMdcHitsCollection*) (DigiMan->GetHitsCollection(THCID));
00125
00126 if(THC){
00127 digisCollection=new BesMdcDigisCollection
00128 (moduleName, collectionName[0]);
00129 NHits=THC->entries();
00130 for(G4int i=0;i<NHits;i++){
00131 layerId = (*THC)[i]->GetLayerNo();
00132 cellId = (*THC)[i]->GetCellNo();
00133 edep = (*THC)[i]->GetEdep();
00134 driftD = (*THC)[i]->GetDriftD();
00135 globalT = (*THC)[i]->GetGlobalT();
00136 theta = (*THC)[i]->GetTheta();
00137 cosTheta = cos(theta);
00138 enterAngle = (*THC)[i]->GetEnterAngle();
00139 posFlag = (*THC)[i]->GetPosFlag();
00140
00141 //Transfer hit pointer to BesMdcCalTransfer
00142 mdcCalPointer->SetHitPointer((*THC)[i]);
00143
00144 //Filter with wire efficiency
00145 if(effFlag==0){
00146 //tempEff = mdcCalPointer->GetEff();
00147 tempEff=mdcTunningSvc->GetEff(layerId,cellId,driftD,cosTheta,posFlag);
00148 }else{
00149 tempEff = layerEff[layerId];
00150 }
00151 fRandom=G4UniformRand();
00152 if(fRandom>tempEff)continue;
00153
00154 //cout<<"layerid "<<layerId<<" cellid "<<cellId<<" theta "<<cosTheta<<" enterangle "<<enterAngle<<endl;
00155 //Drift distance smear
00156 if(smearFlag==0){ //No smear
00157 driftDNew = driftD;
00158 }else if(smearFlag==1){ //Smear from TuningSvc
00159 // mdcTunningSvc->GetRes(layerId,cellId,driftD,cosTheta,posFlag,enterAngle,mean,sigma);
00160 mdcTunningSvc->GetRes2(layerId,cellId,driftD,cosTheta,posFlag,enterAngle,f,mean1,sigma1,mean2,sigma2);
00161
00162 //driftDNew = Smear(driftD,f,mean1,sigma1,mean2,sigma2);
00163
00164 driftDNew = Smear(driftD-(f*mean1+(1-f)*mean2),f,mean1,sigma1,mean2,sigma2);//new method
00165
00166 }else if(smearFlag==2){ //Smear with fixed resolution
00167 driftDNew = Smear(driftD);
00168 }else{
00169 G4cerr<<"MdcDigitizer::worong smearFlag: "<<smearFlag<<G4endl;
00170 }
00171
00172 //Do X-T conversion
00173 driftTNew = mdcCalPointer->D2T(driftDNew);
00174
00175 //Do Q-T correct
00176 driftTNew += mdcCalPointer->GetTimeWalk();
00177
00178 //Add T0
00179 driftTNew += mdcCalPointer->GetT0();
00180
00181 //Add TOF
00182 driftTNew += globalT;
00183
00184 //Signal transfer time on wire
00185 // transferT=Transfer(layerId,cellId,hitPosition);
00186 //driftTNew+=transferT;
00187
00188 if(isnan(driftTNew)){
00189 G4cout<<"MdcDigitizer::error, driftT is nan"<<G4endl;
00190 continue;
00191 }
00192
00193 /*
00194 if(m_G4Svc->MdcRootFlag())
00195 {
00196 m_NHits= NHits;
00197 m_layerId= layerId;
00198 m_cellId= cellId;
00199 m_edep= edep;
00200 m_driftD= driftD;
00201 // m_driftT= driftT;
00202 m_globalT = globalT;
00203 m_enterAngle = enterAngle;
00204 m_driftDNew = driftDNew;
00205 m_driftTNew = driftTNew;
00206 m_theta = theta;
00207 m_tupleMdc ->write();
00208 }
00209 */
00210 BesMdcDigi* newDigi = new BesMdcDigi();
00211 newDigi->SetTrackID((*THC)[i]->GetTrackID());
00212 newDigi->SetLayerNo(layerId);
00213 newDigi->SetCellNo(cellId);
00214 newDigi->SetEdep(edep);
00215 newDigi->SetDriftT(driftTNew);
00216 G4int NbDigis = digisCollection->insert(newDigi);
00217 digiPointer[layerId][cellId]=NbDigis-1;
00218 }
00219
00220 if(noiseFlag==1)AddNoise();
00221 if(noiseFlag==2){
00222 ifstream readNoiseLevel("$MDCSIMROOT/share/noiselevel.txt");
00223 if(!readNoiseLevel.good()){
00224 std::cout<<" Error , noiselevel file not exist "<<std::endl;
00225 }else{
00226 std::cout<<" MdcDigitizer:: Open noiselevel file "<<std::endl;
00227 }
00228 G4int NLayer=mdcGeoPointer->SignalLayerNo();
00229 G4double level;
00230 for(G4int i=0;i<NLayer;i++){
00231 readNoiseLevel>>level;
00232 mixLevel.push_back(level);
00233 }
00234 AddNoise2();
00235 }
00236
00237 if (verboseLevel>0) {
00238 G4cout << "\n-------->digis Collection: in this event they are "
00239 << digisCollection->entries()
00240 << " digis in the MDC chambers: " << G4endl;
00241 digisCollection->PrintAllDigi();
00242 }
00243 StoreDigiCollection(digisCollection);
00244 }
00245
00246 }
|
|
||||||||||||
|
|
|
||||||||||||
|
00092 {
00093 if(layer==-1){
00094 for(G4int i=0; i<43;i++){
00095 layerEff[i]=eff;
00096 }
00097 }else{
00098 layerEff[layer]=eff;
00099 }
00100 }
|
|
|
00046 {effFlag=flag;}
|
|
|
00046 {effFlag=flag;}
|
|
|
00044 {mdcDRes=res;}
|
|
|
00044 {mdcDRes=res;}
|
|
|
00039 {noiseFlag=flag;}
|
|
|
00039 {noiseFlag=flag;}
|
|
|
00041 {noiseLevel=level;}
|
|
|
00041 {noiseLevel=level;}
|
|
|
00040 {noiseType=type;}
|
|
|
00040 {noiseType=type;}
|
|
|
00043 {smearFlag=flag;}
|
|
|
00043 {smearFlag=flag;}
|
|
||||||||||||||||||||||||||||
|
|
|
||||||||||||||||
|
|
|
|
|
|
||||||||||||||||||||||||||||
|
00270 {
00271 G4double r, driftDNew,mean,sigma;
00272 r = G4UniformRand();
00273 if(r<=f){
00274 mean=mean1;
00275 sigma=sigma1;
00276 }else{
00277 mean=mean2;
00278 sigma=sigma2;
00279 }
00280 int times=0;
00281 r = G4RandGauss::shoot();
00282 driftDNew = driftD + sigma * r+mean;
00283 while(driftDNew <= 0){
00284 r = G4RandGauss::shoot();
00285 driftDNew = driftD + sigma * r+mean;
00286 times++;
00287 if(times>10)driftDNew=0.01;
00288 }
00289 return driftDNew;
00290 }
|
|
||||||||||||||||
|
00259 {
00260 G4double r, driftDNew;
00261 r = G4RandGauss::shoot();
00262 driftDNew = driftD + sigma * r+mean;
00263 while(driftDNew <= 0){
00264 r = G4RandGauss::shoot();
00265 driftDNew = driftD + sigma * r+mean;
00266 }
00267 return driftDNew;
00268 }
|
|
|
00248 {
00249 G4double r, driftDNew;
00250 r = G4RandGauss::shoot();
00251 driftDNew = driftD + mdcDRes * r;
00252 while(driftDNew<=0){
00253 r = G4RandGauss::shoot();
00254 driftDNew = driftD + mdcDRes * r;
00255 }
00256 return driftDNew;
00257 }
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1