00001
00011
00012 #include "BesTofDigitizerBrV2.hh"
00013 #include "BesTofHit.hh"
00014 #include "G4DigiManager.hh"
00015 #include "G4RunManager.hh"
00016 #include "Randomize.hh"
00017 #include "G4Poisson.hh"
00018 #include "BesTofDigi.hh"
00019 #include "BesTofGeoParameter.hh"
00020 #include "GaudiKernel/ISvcLocator.h"
00021 #include "GaudiKernel/Bootstrap.h"
00022 #include "GaudiKernel/IDataProviderSvc.h"
00023 #include "G4Svc/IG4Svc.h"
00024 #include "G4Svc/G4Svc.h"
00025 #include "TMath.h"
00026
00027 BesTofDigitizerBrV2::BesTofDigitizerBrV2()
00028 {
00029 ReadData();
00030 m_timeBinSize=0.005;
00031
00032
00033 ISvcLocator* svcLocator = Gaudi::svcLocator();
00034 IG4Svc* tmpSvc;
00035 StatusCode sc = svcLocator->service("G4Svc", tmpSvc);
00036 if(!sc.isSuccess())
00037 {
00038 std::cout << " Could not initialize Realization Service in BesTofDigitizerBrV2" << std::endl;
00039 }
00040 else
00041 {
00042 m_G4Svc = dynamic_cast<G4Svc *>(tmpSvc);
00043 }
00044
00045
00046 IRealizationSvc *tmpReal;
00047 StatusCode scReal = svcLocator->service("RealizationSvc",tmpReal);
00048 if (!scReal.isSuccess())
00049 {
00050 std::cout << " Could not initialize Realization Service in BesTofDigitizerBrV2" << std::endl;
00051 }
00052 else
00053 {
00054 m_RealizationSvc = dynamic_cast<RealizationSvc*>(tmpReal);
00055 }
00056
00057
00058 }
00059
00060 void BesTofDigitizerBrV2::ReadData()
00061 {
00062 BesTofGeoParameter* tofPara = BesTofGeoParameter::GetInstance();
00063
00064 m_scinLength = tofPara->GetBr1L();
00065 m_tau1 = tofPara->GetTau1();
00066 m_tau2 = tofPara->GetTau2();
00067 m_tau3 = tofPara->GetTau3();
00068 m_tauRatio = tofPara->GetTauRatio();
00069 m_refIndex = tofPara->GetRefIndex();
00070 m_phNConst = tofPara->GetPhNConst();
00071 m_Cpe2pmt = tofPara->GetCpe2pmt();
00072 m_rAngle = tofPara->GetRAngle();
00073 m_QE = tofPara->GetQE();
00074 m_CE = tofPara->GetCE();
00075 m_peCorFac = tofPara->GetPeCorFac();
00076
00077 m_ttsMean = tofPara->GetTTSmean();
00078 m_ttsSigma = tofPara->GetTTSsigma();
00079 m_Ce = tofPara->GetCe();
00080 m_LLthresh = tofPara->GetLLthresh();
00081 m_HLthresh = tofPara->GetHLthresh();
00082 m_preGain = tofPara->GetPreGain();
00083 m_noiseSigma = tofPara->GetNoiseSigma();
00084
00085
00086 }
00087
00088 BesTofDigitizerBrV2::~BesTofDigitizerBrV2()
00089 {
00090 ;
00091 }
00092
00093 void BesTofDigitizerBrV2::Digitize(ScintSingle* scint, BesTofDigitsCollection* DC)
00094 {
00095 m_beamTime = m_G4Svc->GetBeamTime() * ns;
00096 m_besTofDigitsCollection = DC;
00097
00098 G4DigiManager* digiManager = G4DigiManager::GetDMpointer();
00099 G4int THCID = digiManager->GetHitsCollectionID("BesTofHitsCollection");
00100 m_THC = (BesTofHitsCollection*) (digiManager->GetHitsCollection(THCID));
00101
00102 if (m_G4Svc->TofRootFlag())
00103 {
00104 m_eTotal = 0;
00105 m_nDigi = 0;
00106 m_partIdMPV = -9;
00107 m_scinNbMPV = -9;
00108 m_edepMPV = 0;
00109 m_nDigiOut = 0;
00110 }
00111
00112 if (m_THC)
00113 {
00114
00115 G4int partId, scinNb, nHits;
00116 G4double edep;
00117 BesTofHit* hit;
00118 partId=scint->GetPartId();
00119 scinNb=scint->GetScinNb();
00120 edep = scint->GetEdep();
00121 nHits=scint->GetHitIndexes()->size();
00122
00123
00124
00125
00126
00127
00128
00129 TofPmtInit();
00130
00131
00132 if (m_G4Svc->TofRootFlag())
00133 {
00134 if (edep>m_edepMPV)
00135 {
00136 m_partIdMPV = partId;
00137 m_scinNbMPV = scinNb;
00138 m_edepMPV = edep;
00139 }
00140 m_eTotal += edep;
00141 m_nDigi ++;
00142
00143 m_partId = partId;
00144 m_scinNb = scinNb;
00145 m_edep = edep;
00146 m_nHits = nHits;
00147 }
00148
00149 if (edep>0.01)
00150 {
00151 for (G4int j=0;j<nHits;j++)
00152 {
00153 hit= (*m_THC)[( *(scint->GetHitIndexes()) )[j]];
00154 TofPmtAccum(hit, scinNb);
00155 }
00156 if (m_G4Svc->TofRootFlag())
00157 {
00158 m_time1st0=m_t1st[0];
00159 m_time1st1=m_t1st[1];
00160 m_timelast0=m_tLast[0];
00161 m_timelast1=m_tLast[1];
00162 m_totalPhot0=m_totalPhot[0];
00163 m_totalPhot1=m_totalPhot[1];
00164 }
00165
00166 TofPmtRspns(scinNb);
00167
00168
00169
00170
00171 G4double temp0 = m_ADC[0]+m_TDC[0];
00172 G4double temp1 = m_ADC[1]+m_TDC[1];
00173
00174
00175 if ( partId==1 && (temp0 > -1998. || temp1 > -1998.))
00176 {
00177
00178 BesTofDigi* digi = new BesTofDigi;
00179 digi->SetTrackIndex(m_trackIndex);
00180 digi->SetPartId(partId);
00181 digi->SetScinNb(scinNb);
00182 digi->SetForwADC( m_ADC[0]) ;
00183 digi->SetBackADC( m_ADC[1]) ;
00184 if (m_TDC[0]>0)
00185 m_TDC[0] = m_TDC[0]+m_beamTime;
00186 if (m_TDC[1]>0)
00187 m_TDC[1] = m_TDC[1]+m_beamTime;
00188 digi->SetForwTDC( m_TDC[0]) ;
00189 digi->SetBackTDC( m_TDC[1]) ;
00190
00191
00192
00193
00194 m_besTofDigitsCollection->insert(digi);
00195
00196 if (m_G4Svc->TofRootFlag())
00197 m_nDigiOut++;
00198
00199 }
00200 if (m_G4Svc->TofRootFlag())
00201 m_tupleTof1->write();
00202
00203 }
00204
00205 if (m_G4Svc->TofRootFlag())
00206 m_tupleTof2->write();
00207 }
00208 }
00209
00210 void BesTofDigitizerBrV2::TofPmtInit()
00211 {
00212 Initialize();
00213
00214 m_t1st[0]=100;
00215 m_t1st[1]=100;
00216 m_tLast[0]=0.;
00217 m_tLast[1]=0;
00218 m_totalPhot[0]=0;
00219 m_totalPhot[1]=0;
00220 for (G4int i=0;i<2;i++)
00221 for (G4int j=0;j<m_profBinN;j++)
00222 m_nPhot[j][i]=0;
00223
00224 if (m_G4Svc->TofRootFlag())
00225 {
00226 m_partId = -9;
00227 m_scinNb = -9;
00228 m_edep = 0;
00229 m_nHits = 0;
00230 m_time1st0 = 100;
00231 m_time1st1 = 100;
00232 m_timelast0 = 0;
00233 m_timelast1 = 0;
00234 m_totalPhot0 = 0;
00235 m_totalPhot1 = 0;
00236 m_NphAllSteps = 0;
00237 m_max0 = 0;
00238 m_max1 = 0;
00239 m_tdc0 = -999;
00240 m_adc0 = -999;
00241 m_tdc1 = -999;
00242 m_adc1 = -999;
00243 }
00244 }
00245
00246 void BesTofDigitizerBrV2::TofPmtAccum(BesTofHit* hit, G4int scinNb)
00247 {
00248 BesTofGeoParameter* tofPara = BesTofGeoParameter::GetInstance();
00249
00250 G4double cvelScint=c_light/m_refIndex/1.07;
00251
00252 G4ThreeVector pos=hit->GetPos();
00253 G4int trackIndex = hit->GetTrackIndex();
00254 G4int partId =hit->GetPartId();
00255 G4double edep=hit->GetEdep();
00256 G4double stepL=hit->GetStepL();
00257 G4double deltaT=hit->GetDeltaT();
00258 G4double timeFlight=hit->GetTime()-m_beamTime;
00259
00260 if (timeFlight < m_globalTime)
00261 {
00262 m_globalTime = timeFlight;
00263 m_trackIndex = trackIndex;
00264 }
00265
00266 G4ThreeVector pDirection=hit->GetPDirection();
00267 G4double nx=pDirection.x();
00268 G4double ny=pDirection.y();
00269 G4double nz=pDirection.z();
00270
00271
00272
00273
00274
00275
00276
00277
00278
00279
00280
00281 G4double thetaProb=1-cos( asin(1.43/m_refIndex));
00282
00283
00284
00285 G4double nMean, nPhoton;
00286
00287 nMean = m_phNConst*BirksLaw(hit);
00288 G4int runId = m_RealizationSvc->getRunId();
00289 if( ( runId>=-11396 && runId<=-8093 ) || ( runId>-80000 && runId<=-23463 ) ) {
00290 nMean = 9000.0*BirksLaw(hit);
00291 }
00292
00293
00294 if(nMean>10)
00295 {
00296 G4double resolutionScale=1.;
00297 G4double sigma=resolutionScale*sqrt(nMean);
00298 nPhoton=G4int(G4RandGauss::shoot(nMean,sigma)+0.5);
00299 }
00300 else
00301 nPhoton=G4int(G4Poisson(nMean));
00302
00303
00304 G4int NphStep;
00305 if (partId==1)
00306 NphStep=G4int(nPhoton*thetaProb*m_rAngle*m_QE*m_CE*m_peCorFac);
00307
00308
00309
00310
00311
00312
00313
00314
00315 if (m_G4Svc->TofRootFlag())
00316 m_NphAllSteps += NphStep;
00317
00318
00319 G4double ddS, ddT;
00320 if (NphStep>0)
00321 {
00322 for (G4int i=0;i<NphStep;i++)
00323 {
00324
00325 ddS=stepL*G4UniformRand();
00326 ddT=deltaT*G4UniformRand();
00327 G4ThreeVector emtPos;
00328 emtPos.setX(pos.x() + nx*ddS);
00329 emtPos.setY(pos.y() + ny*ddS);
00330 emtPos.setZ(pos.z() + nz*ddS);
00331
00332
00333
00334 G4double pathL=0;
00335 G4int forb;
00336 DirectPh(partId, emtPos, pathL, forb);
00337
00338
00339
00340
00341 G4double ran = G4UniformRand();
00342
00343
00344 G4double attenL = m_tofSimSvc->BarAttenLength(scinNb);
00345 attenL = 10.*attenL/0.75;
00346
00347 if (pathL>0 && exp(-pathL/attenL) > ran)
00348 {
00349
00350 G4double scinSwim=pathL/cvelScint;
00351
00352
00353 G4double scinTime=Scintillation(partId);
00354
00355
00356 G4double transitTime=TransitTime();
00357
00358 G4double endTime= timeFlight + ddT + scinSwim + scinTime + transitTime;
00359
00360
00361 if (m_G4Svc->TofRootFlag())
00362 {
00363
00364
00365
00366 m_scinSwim = scinSwim;
00367
00368
00369 m_endTime = endTime;
00370 m_tupleTof3->write();
00371 }
00372
00373 AccuSignal(endTime, forb);
00374
00375
00376 if (m_t1st[forb]>endTime) m_t1st[forb] = endTime;
00377 if (m_tLast[forb]<endTime) m_tLast[forb]= endTime;
00378 }
00379 }
00380 }
00381 }
00382
00383 G4double BesTofDigitizerBrV2::BirksLaw(BesTofHit* hit)
00384 {
00385 const G4double kappa = 0.015*cm/MeV;
00386 const G4String brMaterial = "BC408";
00387 G4double dE = hit->GetEdep();
00388
00389 G4double dX = hit->GetStepL();
00390
00391 G4double charge = hit->GetCharge();
00392 G4double cor_dE = dE;
00393
00394 if(charge!=0.&& dX!=0.)
00395 {
00396 cor_dE = dE/(1+kappa*dE/dX);
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 }
00410 return cor_dE;
00411
00412 }
00413
00414 G4double BesTofDigitizerBrV2::Reflectivity(G4double n1,G4double n2,G4double n3,G4double theta)
00415 {
00416 G4double I1,I2,I3,rp1,rs1,rp2,rs2,Rp1,Rs1,Rp2,Rs2,Rp,Rs;
00417 G4double R=1.0;
00418
00419
00420 I1=theta;
00421 if(I1<asin(n2/n1))
00422 {
00423 I2=asin(sin(I1)*(n1/n2));
00424 rp1=(n1/cos(I1)-n2/cos(I2))/(n1/cos(I1)+n2/cos(I2));
00425 rs1=(n1*cos(I1)-n2*cos(I2))/(n1*cos(I1)+n2*cos(I2));
00426 Rp1=rp1*rp1;
00427 Rs1=rs1*rs1;
00428
00429 I3=asin(sin(I1)*(n1/n3));
00430 rp2=(n2/cos(I2)-n3/cos(I3))/(n2/cos(I2)+n3/cos(I3));
00431 rs2=(n2*cos(I2)-n3*cos(I3))/(n2*cos(I2)+n3*cos(I3));
00432 Rp2=rp2*rp2;
00433 Rs2=rs2*rs2;
00434 Rp=(Rp1+Rp2-2*Rp1*Rp2)/(1-Rp1*Rp2);
00435 Rs=(Rs1+Rs2-2*Rs1*Rs2)/(1-Rs1*Rs2);
00436 R=(Rp+Rs)/2.;
00437 }
00438 return R;
00439 }
00440
00441
00442 G4double BesTofDigitizerBrV2::Reflectivity(G4double n1,G4double n2,G4double theta)
00443 {
00444 G4double I1,I2,rp,rs,Rp,Rs;
00445 G4double R=1.0;
00446
00447
00448 I1=theta;
00449 if (I1<asin(n2/n1))
00450 {
00451 I2=asin(sin(I1)*(n1/n2));
00452 rp=(n1/cos(I1)-n2/cos(I2))/(n1/cos(I1)+n2/cos(I2));
00453 rs=(n1*cos(I1)-n2*cos(I2))/(n1*cos(I1)+n2*cos(I2));
00454 Rp=rp*rp;
00455 Rs=rs*rs;
00456 R=(Rp+Rs)/2.;
00457 }
00458 return R;
00459 }
00460
00461 void BesTofDigitizerBrV2::DirectPh(G4int partId, G4ThreeVector emtPos, G4double& pathL, G4int& forb)
00462 {
00463
00464 const static G4double silicon_oil_index = 1.43;
00465 const static G4double glass_index = 1.532;
00466
00467
00468 G4double cos_span=1-cos( asin(silicon_oil_index/m_refIndex) );
00469
00470 G4double dcos, ran;
00471 ran=G4UniformRand();
00472
00473 dcos=cos_span*(ran*2.0-1.0);
00474
00475 G4double r2=sqrt(emtPos.x()*emtPos.x()+emtPos.y()*emtPos.y());
00476 G4int nRef=0;
00477 G4double costheta,sintheta;
00478 G4double theta,thetaR;
00479 costheta=dcos>0?(1-dcos):(1+dcos);
00480 theta=acos(costheta);
00481 sintheta=sin(theta);
00482 thetaR=asin(1/m_refIndex);
00483 G4double R1;
00484 R1=Reflectivity(m_refIndex,1.0,theta);
00485 G4double ratio1Mean=(CLHEP::pi)*25.5*25.5/(57.1+60.7)/25.0;
00486 G4double ratio2Mean=(19.5/25.5)*(19.5/25.5);
00487 G4double ratio1=G4RandGauss::shoot(ratio1Mean,0.015);
00488 G4double ratio2=G4RandGauss::shoot(ratio2Mean,0.015);
00489 G4double ratio3Mean=0.945*ratio1Mean*ratio2Mean;
00490 G4double ratio3=G4RandGauss::shoot(ratio3Mean,0.015);
00491
00492 G4double R2=1-Reflectivity(m_refIndex,silicon_oil_index, glass_index, theta);
00493 if (dcos > 0 && dcos != 1)
00494 {
00495 if (theta < thetaR)
00496 {
00497 if (G4UniformRand()<ratio1)
00498 {
00499 if (G4UniformRand()<R2)
00500 {
00501 if (G4UniformRand()<ratio2)
00502 {
00503 forb=0;
00504 pathL=(m_scinLength/2-emtPos.z())/costheta;
00505 }
00506 }
00507 else
00508 {
00509 if (G4UniformRand()<ratio3)
00510 {
00511 forb=1;
00512 pathL=(1.5*m_scinLength-emtPos.z())/costheta;
00513 }
00514 }
00515 }
00516 else
00517 {
00518 if (G4UniformRand()<R1)
00519 {
00520 G4double tempran=G4UniformRand();
00521 if (tempran<ratio3)
00522 {
00523 forb=1;
00524 pathL=(1.5*m_scinLength-emtPos.z())/costheta;
00525 }
00526 else if (tempran>ratio1&&G4UniformRand()<R1&&G4UniformRand()<ratio3)
00527 {
00528 forb=0;
00529 pathL=(2.5*m_scinLength-emtPos.z())/costheta;
00530 }
00531 }
00532 }
00533 }
00534 else
00535 {
00536 if (G4UniformRand()<ratio1)
00537 {
00538 if (G4UniformRand()<R2)
00539 {
00540 if (G4UniformRand()<ratio2)
00541 {
00542 forb=0;
00543 pathL=(m_scinLength/2-emtPos.z())/costheta;
00544 }
00545 }
00546 else
00547 {
00548 if (G4UniformRand()<ratio3)
00549 {
00550 forb=1;
00551 pathL=(1.5*m_scinLength-emtPos.z())/costheta;
00552 }
00553 }
00554 }
00555 else
00556 {
00557 G4double tempran=G4UniformRand();
00558 if (tempran<ratio3)
00559 {
00560 forb=1;
00561 pathL=(1.5*m_scinLength-emtPos.z())/costheta;
00562 }
00563 else if (tempran>ratio1 && G4UniformRand()<ratio3)
00564 {
00565 forb=0;
00566 pathL=(2.5*m_scinLength-emtPos.z())/costheta;
00567 }
00568 }
00569 }
00570 }
00571 else if (dcos < 0 && dcos != -1)
00572 {
00573 if (theta < thetaR)
00574 {
00575 if (G4UniformRand()<ratio1)
00576 {
00577 if (G4UniformRand()<R2)
00578 {
00579 if (G4UniformRand()<ratio2)
00580 {
00581 forb=1;
00582 pathL=(m_scinLength/2+emtPos.z())/costheta;
00583 }
00584 }
00585 else
00586 {
00587 if (G4UniformRand()<ratio3)
00588 {
00589 forb=0;
00590 pathL=(1.5*m_scinLength+emtPos.z())/costheta;
00591 }
00592 }
00593 }
00594 else
00595 {
00596 if (G4UniformRand()<R1)
00597 {
00598 G4double tempran=G4UniformRand();
00599 if (tempran<ratio3)
00600 {
00601 forb=0;
00602 pathL=(1.5*m_scinLength+emtPos.z())/costheta;
00603 }
00604 else if (tempran>ratio1&&G4UniformRand()<R1&&G4UniformRand()<ratio3)
00605 {
00606 forb=1;
00607 pathL=(2.5*m_scinLength+emtPos.z())/costheta;
00608 }
00609 }
00610 }
00611 }
00612 else
00613 {
00614 if (G4UniformRand()<ratio1)
00615 {
00616 if (G4UniformRand()<R2)
00617 {
00618 if (G4UniformRand()<ratio2)
00619 {
00620 forb=1;
00621 pathL=(m_scinLength/2+emtPos.z())/costheta;
00622 }
00623 }
00624 else
00625 {
00626 if (G4UniformRand()<ratio3)
00627 {
00628 forb=0;
00629 pathL=(1.5*m_scinLength+emtPos.z())/costheta;
00630 }
00631 }
00632 }
00633 else
00634 {
00635 G4double tempran=G4UniformRand();
00636 if (tempran<ratio3)
00637 {
00638 forb=0;
00639 pathL=(1.5*m_scinLength+emtPos.z())/costheta;
00640 }
00641 else if (tempran>ratio1 && G4UniformRand()<ratio3)
00642 {
00643 forb=1;
00644 pathL=(2.5*m_scinLength+emtPos.z())/costheta;
00645 }
00646 }
00647 }
00648 }
00649
00650 G4double convFactor=180./3.1415926;
00651 if (theta>asin(1)-thetaR)
00652 {
00653 G4double alpha = G4UniformRand()*asin(1.);
00654 G4int nRef1 = pathL*sintheta*cos(alpha)/50.0+0.5;
00655 G4int nRef2 = pathL*sintheta*sin(alpha)/58.9+0.5;
00656 G4double beta1=acos(sintheta*cos(alpha));
00657 G4double beta2=acos(sintheta*sin(alpha));
00658 beta2 -= 3.75/convFactor;
00659 G4double R21,R22;
00660 R21=Reflectivity(m_refIndex,1.0,beta1);
00661 R22=Reflectivity(m_refIndex,1.0,beta2);
00662 for (G4int i=0;i<nRef1;i++)
00663 {
00664 if (G4UniformRand()<(1-R21) && G4UniformRand()<0.15)
00665 pathL=-9;
00666 }
00667 for (G4int i=0;i<nRef2;i++)
00668 {
00669 if (G4UniformRand()<(1-R22) && G4UniformRand()<0.15)
00670 pathL=-9;
00671 }
00672 }
00673 }
00674
00675
00676
00677
00678
00679
00680
00681
00682
00683
00684
00685
00686
00687
00688
00689
00690
00691
00692
00693
00694
00695
00696
00697 G4double BesTofDigitizerBrV2::Scintillation(G4int partId)
00698 {
00699 G4double tmp_tauRatio,tmp_tau1,tmp_tau2,tmp_tau3;
00700 tmp_tauRatio = m_tauRatio;
00701 tmp_tau1 = m_tau1;
00702 tmp_tau2 = m_tau2;
00703 tmp_tau3 = m_tau3;
00704 G4double UniformR = tmp_tauRatio/(1+tmp_tauRatio);
00705 G4double EmissionTime;
00706 if (G4UniformRand()>UniformR) {
00707 while (1) {
00708 EmissionTime = -tmp_tau2*log( G4UniformRand() );
00709 if (G4UniformRand()-exp(EmissionTime/tmp_tau2-EmissionTime/tmp_tau1)>1.E-8)
00710 break;
00711 }
00712 }
00713 else EmissionTime = -tmp_tau3*log( G4UniformRand() );
00714 return EmissionTime;
00715 }
00716
00717
00718
00719
00720
00721
00722
00723
00724
00725
00726
00727
00728
00729
00730
00731
00732
00733
00734
00735
00736
00737
00738
00739
00740
00741
00742
00743
00744 G4double BesTofDigitizerBrV2::TransitTime()
00745 {
00746
00747 return G4RandGauss::shoot(m_ttsMean,m_ttsSigma);
00748 }
00749
00750 void BesTofDigitizerBrV2::AccuSignal(G4double endTime, G4int forb)
00751 {
00752 G4int ihst;
00753 ihst=G4int(endTime/m_timeBinSize);
00754 if (ihst>0 &&ihst<m_profBinN)
00755 {
00756 m_nPhot[ihst][forb]=m_nPhot[ihst][forb]+1;
00757 m_totalPhot[forb]=m_totalPhot[forb]+1;
00758 }
00759 }
00760
00761 void BesTofDigitizerBrV2::TofPmtRspns(G4int scinNb)
00762 {
00763
00764 BesTofGeoParameter* tofPara = BesTofGeoParameter::GetInstance();
00765
00766
00767 G4double snpe[m_snpeBinN][2];
00768
00769
00770
00771
00772 G4double norma_const;
00773 G4double echarge=1.6e-7;
00774
00775
00776 G4double profPmt[m_profBinN][2];
00777
00778 G4double t;
00779 G4double tau;
00780 G4int n1, n2, ii;
00781 G4int phtn;
00782
00783 for (G4int i=0;i<m_snpeBinN;i++)
00784 {
00785 tau = tofPara->GetBrERiseTime(scinNb);
00786 norma_const = sqrt(CLHEP::pi)*tau*tau*tau/4.0;
00787 t = (i+1)*m_timeBinSize;
00788 snpe[i][0] = m_Ce*t*t*exp(- (t/tau) * (t/tau) )/norma_const;
00789 }
00790
00791 for (G4int i=0;i<m_snpeBinN;i++)
00792 {
00793 tau = tofPara->GetBrWRiseTime(scinNb);
00794 norma_const = sqrt(CLHEP::pi)*tau*tau*tau/4.0;
00795 t = (i+1)*m_timeBinSize;
00796 snpe[i][1] = m_Ce*t*t*exp(- (t/tau) * (t/tau) )/norma_const;
00797 }
00798
00799 G4int fb=2;
00800 G4double Npoisson;
00801 G4double pmtGain0,pmtGain,relativeGain,smearPMTgain;
00802 G4double smearADC[2] = {0.,0.};
00803 G4int runId = m_RealizationSvc->getRunId();
00804 pmtGain0 = m_tofSimSvc->BarPMTGain();
00805
00806
00807
00808
00809
00810
00811
00812
00813
00814
00815
00816
00817
00818
00819 G4double timeSmear = G4RandGauss::shoot(0,0.020);
00820 if (runId>=-22913 && runId<=-20448) {
00821 timeSmear = G4RandGauss::shoot(0,0.040);
00822 }
00823 else if( ( runId>=-11396 && runId<=-8093 ) || ( runId>-80000 && runId<=-23463 ) ) {
00824 timeSmear = G4RandGauss::shoot(0,0.025);
00825 }
00826
00827 for (G4int j=0; j<fb; j++)
00828 {
00829 if (m_totalPhot[j] > 0)
00830 {
00831 n1=G4int(m_t1st[j]/m_timeBinSize);
00832 n2=G4int(m_tLast[j]/m_timeBinSize);
00833
00834
00835 for (G4int i=0;i<m_profBinN;i++)
00836 profPmt[i][j]=0.0;
00837
00838
00839 n2 = n2<m_profBinN ? n2:m_profBinN;
00840 for (G4int i=n1;i<n2;i++)
00841 {
00842 phtn=m_nPhot[i][j];
00843 if (phtn>0)
00844 {
00845 while(1) {
00846 Npoisson = G4Poisson(10.0);
00847 if(Npoisson>0.) break;
00848 }
00849 while(1) {
00850
00851
00852 relativeGain = j ? m_tofSimSvc->BarGain2(scinNb) : m_tofSimSvc->BarGain1(scinNb);
00853 pmtGain = pmtGain0*relativeGain;
00854 smearPMTgain = G4RandGauss::shoot(pmtGain,pmtGain/sqrt(Npoisson));
00855
00856 if(smearPMTgain>0) break;
00857 }
00858 smearADC[j] += phtn*smearPMTgain;
00859
00860 for (G4int ihst=0; ihst<m_snpeBinN; ihst++)
00861 {
00862 ii=i+ihst;
00863 if (ii<m_profBinN)
00864 profPmt[ii][j] += smearPMTgain*phtn*snpe[ihst][j];
00865 else
00866 break;
00867 }
00868 }
00869 }
00870
00871
00872 for (int i=0;i<m_profBinN;i++)
00873 {
00874 if (profPmt[i][j]>0)
00875 profPmt[i][j] = m_preGain*profPmt[i][j]+G4RandGauss::shoot(0,m_noiseSigma);
00876 }
00877
00878
00879 G4double max=0;
00880 for (int i=n1;i<m_profBinN;i++)
00881 {
00882 if (profPmt[i][j]>max)
00883 max=profPmt[i][j];
00884 }
00885 if (m_G4Svc->TofRootFlag())
00886 {
00887 if (j==0) m_max0=max;
00888 else m_max1=max;
00889 }
00890
00891 G4double tmp_HLthresh, tmp_LLthresh, adcFactor;
00892 G4double ratio;
00893
00894 if(runId>=-80000 && runId<=-9484)
00895 {
00896
00897
00898 tmp_HLthresh = m_HLthresh;
00899 tmp_LLthresh = m_LLthresh;
00900 adcFactor = 5.89;
00901 }
00902 else
00903 {
00904
00905
00906 tmp_HLthresh = 600;
00907 tmp_LLthresh = 150;
00908 adcFactor = 4.8;
00909 }
00910
00911
00912
00913
00914
00915
00916
00917
00918 tmp_HLthresh = m_tofSimSvc->BarHighThres();
00919 tmp_LLthresh = m_tofSimSvc->BarLowThres();
00920 ratio = m_tofSimSvc->BarConstant();
00921
00922
00923 if (max>=tmp_HLthresh)
00924 {
00925 for (int i=0;i<m_profBinN;i++)
00926 {
00927 if ( profPmt[i][j] >= tmp_LLthresh )
00928 {
00929 m_TDC[j] = i*m_timeBinSize + timeSmear +G4RandGauss::shoot(0,0.025);
00930
00931 if( ( runId>=-11396 && runId<=-8093 ) || ( runId>-80000 && runId<=-23463 ) ) {
00932 if( scinNb<88 ) {
00933 m_TDC[j] = i*m_timeBinSize + timeSmear +G4RandGauss::shoot(0,0.055);
00934 }
00935 else {
00936 m_TDC[j] = i*m_timeBinSize + timeSmear +G4RandGauss::shoot(0,0.062);
00937 }
00938 }
00939
00940 if( m_G4Svc->TofSaturationFlag())
00941 {
00942
00943 double x = m_preGain*smearADC[j]*echarge*ratio;
00944 if (j==0)
00945 {
00946 m_ADC[j] = m_tofQElecSvc->BQChannel1(scinNb,x);
00947 }
00948 else
00949 {
00950 m_ADC[j] = m_tofQElecSvc->BQChannel2(scinNb,x);
00951 }
00952 }
00953 else
00954 m_ADC[j] = m_preGain*smearADC[j]*echarge*adcFactor;
00955
00956
00957 if (m_G4Svc->TofRootFlag())
00958 {
00959 if (j==0) {
00960 m_tdc0 = m_TDC[0];
00961 m_adc0 = m_ADC[0];
00962 }
00963 else {
00964 m_tdc1 = m_TDC[1];
00965 m_adc1 = m_ADC[1];
00966 }
00967 }
00968 break;
00969 }
00970 }
00971 }
00972 }
00973 }
00974 }