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