00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011 #include "BesEmcDigitizer.hh"
00012 #include "BesEmcDigi.hh"
00013 #include "BesEmcHit.hh"
00014 #include "BesEmcWaveform.hh"
00015 #include "G4DigiManager.hh"
00016 #include "BesEmcParameter.hh"
00017 #include "Randomize.hh"
00018
00019 #include "GaudiKernel/Bootstrap.h"
00020 #include "GaudiKernel/ISvcLocator.h"
00021 #include "G4Svc/G4Svc.h"
00022 #include "EmcCalibConstSvc/EmcCalibConstSvc.h"
00023
00024 BesEmcDigitizer::BesEmcDigitizer(G4String modName)
00025 :G4VDigitizerModule(modName),m_emcCalibConstSvc(0)
00026 {
00027 collectionName.push_back("BesEmcDigitsCollection");
00028 m_besEmcDigitsCollection = 0;
00029
00030
00031 ISvcLocator* svcLocator = Gaudi::svcLocator();
00032 IG4Svc* iG4Svc;
00033 StatusCode sc=svcLocator->service("G4Svc", iG4Svc);
00034 m_G4Svc=dynamic_cast<G4Svc *>(iG4Svc);
00035
00036
00037 if(m_G4Svc->EmcRootFlag())
00038 {
00039 m_tupleEmc1 = m_G4Svc->GetTupleEmc1();
00040 sc = m_tupleEmc1->addItem("partId",m_partId);
00041 sc = m_tupleEmc1->addItem("nTheta",m_nTheta);
00042 sc = m_tupleEmc1->addItem("nPhi",m_nPhi);
00043 sc = m_tupleEmc1->addItem("edep",m_eDep);
00044 sc = m_tupleEmc1->addItem("nHits",m_nHits);
00045 sc = m_tupleEmc1->addItem("adc",m_adc);
00046 sc = m_tupleEmc1->addItem("tdc",m_tdc);
00047
00048 m_tupleEmc2 = m_G4Svc->GetTupleEmc2();
00049 sc = m_tupleEmc2->addItem("etot",m_eTot);
00050 sc = m_tupleEmc2->addItem("nDigi",m_nDigi);
00051 }
00052
00053
00054 sc = svcLocator->service("EmcCalibConstSvc", m_emcCalibConstSvc);
00055 if(sc != StatusCode::SUCCESS) {
00056 G4cout << "BesEmcDigitizer Error: Can't get EmcCalibConstSvc." << G4endl;
00057 }
00058
00059 }
00060
00061 BesEmcDigitizer::~BesEmcDigitizer()
00062 {
00063 }
00064
00065 void BesEmcDigitizer::Initialize()
00066 {
00067 }
00068
00069 void BesEmcDigitizer::Digitize()
00070 {
00071 Initialize();
00072
00073
00074 m_besEmcDigitsCollection = new BesEmcDigitsCollection
00075 ("BesEmcDigitizer","BesEmcDigitsCollection");
00076 G4DigiManager* DigiMan = G4DigiManager::GetDMpointer();
00077
00078
00079
00080 G4int EHCID;
00081 EHCID = DigiMan->GetHitsCollectionID("BesEmcHitsCollection");
00082
00083
00084 BesEmcHitsCollection* EHC = 0;
00085 EHC = (BesEmcHitsCollection*) (DigiMan->GetHitsCollection(EHCID));
00086
00087 if (EHC)
00088 {
00089
00090 m_crystalGroup = new vector<CrystalSingle*>;
00091 GroupHits(EHC);
00092 G4int size=m_crystalGroup->size();
00093 CrystalSingle* cryst;
00094 G4int partId, nTheta, nPhi, nHits;
00095 G4double eTot=0, eDigi;
00096 BesEmcHit* hit;
00097
00098 G4double coherentNoise = RandGauss::shoot()*m_G4Svc->EmcCoherentNoise();
00099
00100 for(G4int i=0;i<size;i++)
00101 {
00102 cryst = (*m_crystalGroup)[i];
00103 partId = cryst->GetPartId();
00104 nTheta = cryst->GetNTheta();
00105 nPhi = cryst->GetNPhi();
00106 nHits= cryst->GetHitIndexes()->size();
00107 eDigi = cryst->GetEdep();
00108 eTot += eDigi;
00109
00110 BesEmcDigi *digi = new BesEmcDigi;
00111 BesEmcWaveform *wave = digi->GetWaveform();
00112 G4long bin = 0;
00113
00114 const int indexSize = 200;
00115 G4double e[indexSize];
00116 for(G4int i=0;i<indexSize;i++)
00117 e[i]=0;
00118 G4int index=0;
00119 G4double energy=0;
00120
00121 for(G4int j=0;j<nHits;j++)
00122 {
00123 hit= (*EHC)[( *(cryst->GetHitIndexes()) )[j]];
00124 energy = hit->GetEdepCrystal();
00125 index = hit->GetTrackIndex();
00126 if(index<indexSize&&index>=0)
00127 e[index]+=energy;
00128 else
00129 G4cout<<"Track index overload!"<<G4endl;
00130 }
00131
00132 G4double maxi=e[0];
00133 for(G4int i=1;i<indexSize;i++)
00134 {
00135 if(e[i]>maxi)
00136 {
00137 maxi = e[i];
00138 index = i;
00139 }
00140 }
00141
00142 if(eDigi>0)
00143 {
00144 digi->SetPartId(partId);
00145 digi->SetThetaNb(nTheta);
00146 digi->SetPhiNb(nPhi);
00147 digi->SetEnergy(eDigi);
00148 digi->SetTime(m_G4Svc->EmcTime());
00149 digi->SetTrackIndex(index);
00150
00151 wave->updateWaveform(digi);
00152 if(m_G4Svc->EmcNoiseLevel()>0)
00153 wave->addElecNoise(m_G4Svc->EmcIncoherentNoise(),coherentNoise);
00154
00155
00156 m_energy = wave->max(bin);
00157
00158 m_energy -= 0.46*MeV;
00159 wave->digitize();
00160 wave->max(bin);
00161
00162 if(m_G4Svc->EmcLightOutput())
00163 {
00164 G4int index = m_emcCalibConstSvc->getIndex(partId,nTheta,nPhi);
00165 G4double adc2e = m_emcCalibConstSvc->getDigiCalibConst(index);
00166
00167
00168 if (m_G4Svc->EmcElecSaturation()==1){
00169 G4double emaxData = m_emcCalibConstSvc->getCrystalEmaxData(index);
00170
00171 if (emaxData>0) {
00172
00173 adc2e=emaxData/2.5;
00174 }
00175 }
00176
00177 if(adc2e<=1e-5)
00178 {
00179 m_energy = 0;
00180 }
00181 else
00182 {
00183
00184 m_energy /= adc2e;
00185
00186
00187
00188 }
00189 }
00190
00191
00192 if(m_G4Svc->EmcRootFlag())
00193 {
00194 m_partId = partId;
00195 m_nTheta = nTheta;
00196 m_nPhi = nPhi;
00197 m_eDep = eDigi;
00198 m_nHits = nHits;
00199 m_adc = m_energy;
00200 m_tdc = bin;
00201 m_tupleEmc1->write();
00202 }
00203
00204 digi->SetEnergy(m_energy);
00205 digi->SetTime(bin);
00206 digi->SetWaveform(wave);
00207 m_besEmcDigitsCollection->insert(digi);
00208 }
00209 }
00210
00211
00212 if(m_G4Svc->EmcNoiseLevel()==2)
00213 AddNoise5x5(coherentNoise);
00214 else if(m_G4Svc->EmcNoiseLevel()==3)
00215 ;
00216
00217
00218
00219 if(m_G4Svc->EmcRootFlag())
00220 {
00221 m_eTot = eTot;
00222 m_nDigi = size;
00223 m_tupleEmc2->write();
00224 }
00225
00226 StoreDigiCollection(m_besEmcDigitsCollection);
00227
00228 for(size_t i=0;i<m_crystalGroup->size();i++)
00229 {
00230 delete (*m_crystalGroup)[i];
00231 }
00232 m_crystalGroup->clear();
00233 delete m_crystalGroup;
00234 }
00235 }
00236
00237 void BesEmcDigitizer::GroupHits(BesEmcHitsCollection* m_EHC)
00238 {
00239 G4int partId, nTheta, nPhi, size, flag;
00240 G4double edep;
00241 BesEmcHit* hit;
00242 G4int nHits = m_EHC->entries();
00243
00244
00245 for(G4int i=0;i<nHits;i++)
00246 {
00247 hit = (*m_EHC)[i];
00248 partId = hit->GetPartId();
00249 nTheta = hit->GetNumThetaCrystal();
00250 nPhi = hit->GetNumPhiCrystal();
00251 edep = hit->GetEdepCrystal();
00252 size = m_crystalGroup->size();
00253 flag=0;
00254
00255 if(size>0)
00256 {
00257 CrystalSingle* oldCryst;
00258 for(G4int j=0; j<size;j++)
00259 {
00260 oldCryst = (*m_crystalGroup)[j];
00261 if((oldCryst->GetNTheta()==nTheta)&&(oldCryst->GetNPhi()==nPhi)&&(oldCryst->GetPartId()==partId))
00262 {
00263 oldCryst->GetHitIndexes()->push_back(i);
00264 oldCryst->AddEdep(edep);
00265 flag=1;
00266 break;
00267 }
00268 }
00269 }
00270
00271 if(flag==0)
00272 {
00273 CrystalSingle* newCryst = new CrystalSingle;
00274 newCryst->SetPartId(partId);
00275 newCryst->SetNTheta(nTheta);
00276 newCryst->SetNPhi(nPhi);
00277 newCryst->SetEdep(edep);
00278 newCryst->GetHitIndexes()->push_back(i);
00279 m_crystalGroup->push_back(newCryst);
00280 }
00281
00282 }
00283 }
00284
00285 void BesEmcDigitizer::AddNoise5x5(G4double coherentNoise)
00286 {
00287 BesEmcParameter& emcPara = BesEmcParameter::GetInstance();
00288 vector<BesEmcDigi*>* vecDC = m_besEmcDigitsCollection->GetVector();
00289 G4int nDigi = m_besEmcDigitsCollection->entries();
00290 G4int partMax,thetaMax,phiMax;
00291 partMax=thetaMax=phiMax=-99;
00292 G4double eMax = 0;
00293
00294 for(G4int i=0;i<nDigi;i++) {
00295 BesEmcDigi *digi = (*vecDC)[i];
00296 double eDigi = digi->GetEnergy();
00297 if(eDigi>eMax) {
00298 eMax = eDigi;
00299 partMax = digi->GetPartId();
00300 thetaMax = digi->GetThetaNb();
00301 phiMax = digi->GetPhiNb();
00302 }
00303 }
00304
00305 if(partMax==1) {
00306 G4int thetan,thetap,phin,phip;
00307 thetan = thetaMax-2;
00308 thetap = thetaMax+2;
00309 phin = phiMax-2;
00310 phip = phiMax+2;
00311
00312 if(thetaMax==0) {
00313 thetan = thetaMax;
00314 } else if(thetaMax==1) {
00315 thetan = thetaMax-1;
00316 } else if(thetaMax==emcPara.GetBSCNbTheta()*2-1) {
00317 thetap = thetaMax;
00318 } else if(thetaMax==emcPara.GetBSCNbTheta()*2-2) {
00319 thetap = thetaMax+1;
00320 }
00321
00322 if(phiMax==0) {
00323 phin = emcPara.GetBSCNbPhi()-2;
00324 } else if(phiMax==1) {
00325 phin = emcPara.GetBSCNbPhi()-2;
00326 } else if(phiMax==emcPara.GetBSCNbPhi()-1) {
00327 phip = 1;
00328 } else if(phiMax==emcPara.GetBSCNbPhi()-2) {
00329 phip = 0;
00330 }
00331
00332 for(G4int theta=thetan;theta<=thetap;theta++) {
00333 for(G4int phi=phin;phi<=phip;phi++) {
00334 G4bool flag = true;
00335
00336 if(nDigi>0) {
00337 for(G4int i=0;i<nDigi;i++) {
00338 BesEmcDigi *digi = (*vecDC)[i];
00339 if( partMax == digi->GetPartId()
00340 && theta == digi->GetThetaNb()
00341 && phi == digi->GetPhiNb() ) {
00342 flag=false;
00343 break;
00344 }
00345 }
00346 }
00347
00348 if(flag) {
00349 BesEmcDigi *digi = new BesEmcDigi;
00350 BesEmcWaveform *wave = digi->GetWaveform();
00351 digi->SetPartId(partMax);
00352 digi->SetThetaNb(theta);
00353 digi->SetPhiNb(phi);
00354 digi->SetEnergy(0);
00355 digi->SetTime(m_G4Svc->EmcTime());
00356 digi->SetTrackIndex(-9);
00357
00358 wave->updateWaveform(digi);
00359 wave->addElecNoise(m_G4Svc->EmcIncoherentNoise(),coherentNoise);
00360 wave->digitize();
00361
00362 G4long bin = 0;
00363 m_energy = wave->max(bin);
00364
00365 if(m_G4Svc->EmcLightOutput()) {
00366 m_energy *= emcPara.GetLightOutput(partMax,theta,phi);
00367 }
00368
00369 digi->SetEnergy(m_energy);
00370 digi->SetTime(bin);
00371 digi->SetWaveform(wave);
00372
00373
00374 G4double ecorr;
00375 if(m_energy<625.) {
00376 ecorr = -0.1285*log(m_energy/6805.);
00377 } else {
00378 ecorr = -2.418+9.513e-4*m_energy;
00379 }
00380
00381 if(m_energy-ecorr>m_G4Svc->EmcNoiseThreshold()) {
00382 m_besEmcDigitsCollection->insert(digi);
00383 } else {
00384 delete digi;
00385 }
00386 }
00387 }
00388 }
00389
00390 }
00391 }
00392
00393 void BesEmcDigitizer::AddNoiseAll(G4double coherentNoise)
00394 {
00395 BesEmcParameter& emcPara = BesEmcParameter::GetInstance();
00396 vector<BesEmcDigi*>* vecDC = m_besEmcDigitsCollection->GetVector();
00397 G4int nDigi = m_besEmcDigitsCollection->entries();
00398
00399
00400 for(G4int part=0;part<3;part++) {
00401
00402 G4int thetaNb;
00403 if(part == 1) {
00404 thetaNb = emcPara.GetBSCNbTheta()*2;
00405 } else {
00406 thetaNb = 6;
00407 }
00408
00409 for(G4int theta=0;theta<thetaNb;theta++) {
00410
00411 G4int phiNb;
00412 if(part == 1) {
00413 phiNb = emcPara.GetBSCNbPhi();
00414 } else {
00415 phiNb = emcPara.GetCryInOneLayer(theta);
00416 }
00417
00418 for(G4int phi=0;phi<phiNb;phi++) {
00419
00420 G4bool flag = true;
00421
00422 if(nDigi>0) {
00423
00424
00425 for(G4int i=0;i<nDigi;i++) {
00426 BesEmcDigi *digi = (*vecDC)[i];
00427 if( part == digi->GetPartId()
00428 && theta == digi->GetThetaNb()
00429 && phi == digi->GetPhiNb() ) {
00430
00431 flag=false;
00432 break;
00433 }
00434 }
00435 }
00436
00437 if(flag) {
00438 BesEmcDigi *digi = new BesEmcDigi;
00439 digi->SetTrackIndex(-9);
00440 digi->SetPartId(part);
00441 digi->SetThetaNb(theta);
00442 digi->SetPhiNb(phi);
00443
00444 bool fastSimulation = true;
00445 if(fastSimulation) {
00446
00447 m_energy = RandGauss::shoot()*m_G4Svc->EmcNoiseSigma();
00448 m_energy += m_G4Svc->EmcNoiseMean();
00449 digi->SetTime((G4int)(G4UniformRand()*60));
00450
00451 } else {
00452
00453 BesEmcWaveform *wave = digi->GetWaveform();
00454 digi->SetTime(m_G4Svc->EmcTime());
00455
00456 wave->updateWaveform(digi);
00457 wave->addElecNoise(m_G4Svc->EmcIncoherentNoise(),coherentNoise);
00458 wave->digitize();
00459
00460 G4long bin = 0;
00461 m_energy = wave->max(bin);
00462 digi->SetTime(bin);
00463 digi->SetWaveform(wave);
00464 }
00465
00466 if(m_G4Svc->EmcLightOutput()) {
00467 m_energy *= emcPara.GetLightOutput(part,theta,phi);
00468 }
00469 digi->SetEnergy(m_energy);
00470
00471
00472 G4double ecorr;
00473 if(m_energy<625.) {
00474 ecorr = -0.1285*log(m_energy/6805.);
00475 } else {
00476 ecorr = -2.418+9.513e-4*m_energy;
00477 }
00478
00479 if(m_energy-ecorr>m_G4Svc->EmcNoiseThreshold()) {
00480 m_besEmcDigitsCollection->insert(digi);
00481 } else {
00482 delete digi;
00483 }
00484 }
00485
00486 }
00487 }
00488 }
00489 }