#include <BesTofDigitizerBrV2.hh>
Inheritance diagram for BesTofDigitizerBrV2:
Public Member Functions | |
| void | AccuSignal (G4double, G4int) |
| void | AccuSignal (G4double, G4int) |
| BesTofDigitizerBrV2 () | |
| BesTofDigitizerBrV2 () | |
| G4double | BirksLaw (BesTofHit *hit) |
| G4double | BirksLaw (BesTofHit *hit) |
| virtual void | Digitize (ScintSingle *, BesTofDigitsCollection *) |
| virtual void | Digitize (ScintSingle *, BesTofDigitsCollection *) |
| void | DirectPh (G4int, G4ThreeVector, G4double &, G4int &) |
| void | DirectPh (G4int, G4ThreeVector, G4double &, G4int &) |
| void | Initialize () |
| void | Initialize () |
| void | ReadData () |
| void | ReadData () |
| G4double | Reflectivity (G4double n1, G4double n2, G4double n3, G4double theta) |
| G4double | Reflectivity (G4double n1, G4double n2, G4double theta) |
| G4double | Reflectivity (G4double n1, G4double n2, G4double n3, G4double theta) |
| G4double | Reflectivity (G4double n1, G4double n2, G4double theta) |
| G4double | Scintillation (G4int) |
| G4double | Scintillation (G4int) |
| void | TofPmtAccum (BesTofHit *, G4int) |
| void | TofPmtAccum (BesTofHit *, G4int) |
| void | TofPmtInit () |
| void | TofPmtInit () |
| void | TofPmtRspns (G4int) |
| void | TofPmtRspns (G4int) |
| G4double | TransitTime () |
| G4double | TransitTime () |
| ~BesTofDigitizerBrV2 () | |
| ~BesTofDigitizerBrV2 () | |
Protected Attributes | |
| G4double | m_ADC [2] |
| BesTofDigitsCollection * | m_besTofDigitsCollection |
| BesTofDigitsCollection * | m_besTofDigitsCollection |
| G4double | m_globalTime |
| G4double | m_TDC [2] |
| BesTofHitsCollection * | m_THC |
| BesTofHitsCollection * | m_THC |
| ITofCaliSvc * | m_tofCaliSvc |
| ITofCaliSvc * | m_tofCaliSvc |
| ITofQElecSvc * | m_tofQElecSvc |
| ITofQElecSvc * | m_tofQElecSvc |
| ITofSimSvc * | m_tofSimSvc |
| ITofSimSvc * | m_tofSimSvc |
| G4int | m_trackIndex |
Static Protected Attributes | |
| NTuple::Item< double > | m_adc0 |
| NTuple::Item< double > | m_adc0 |
| NTuple::Item< double > | m_adc1 |
| NTuple::Item< double > | m_adc1 |
| bool | m_booked = false |
| NTuple::Item< double > | m_ddT |
| NTuple::Item< double > | m_ddT |
| NTuple::Item< double > | m_edep |
| NTuple::Item< double > | m_edep |
| NTuple::Item< double > | m_edepHit |
| NTuple::Item< double > | m_edepHit |
| NTuple::Item< double > | m_edepMPV |
| NTuple::Item< double > | m_edepMPV |
| NTuple::Item< double > | m_endTime |
| NTuple::Item< double > | m_endTime |
| NTuple::Item< double > | m_eTotal |
| NTuple::Item< double > | m_eTotal |
| NTuple::Item< double > | m_forb |
| NTuple::Item< double > | m_forb |
| NTuple::Item< double > | m_max0 |
| NTuple::Item< double > | m_max0 |
| NTuple::Item< double > | m_max1 |
| NTuple::Item< double > | m_max1 |
| NTuple::Item< double > | m_nDigi |
| NTuple::Item< double > | m_nDigi |
| NTuple::Item< double > | m_nDigiOut |
| NTuple::Item< double > | m_nDigiOut |
| NTuple::Item< double > | m_nHits |
| NTuple::Item< double > | m_nHits |
| NTuple::Item< double > | m_NphAllSteps |
| NTuple::Item< double > | m_NphAllSteps |
| NTuple::Item< double > | m_partId |
| NTuple::Item< double > | m_partId |
| NTuple::Item< double > | m_partIdMPV |
| NTuple::Item< double > | m_partIdMPV |
| NTuple::Item< double > | m_scinNb |
| NTuple::Item< double > | m_scinNb |
| NTuple::Item< double > | m_scinNbMPV |
| NTuple::Item< double > | m_scinNbMPV |
| NTuple::Item< double > | m_scinSwim |
| NTuple::Item< double > | m_scinSwim |
| NTuple::Item< double > | m_scinTime |
| NTuple::Item< double > | m_scinTime |
| NTuple::Item< double > | m_tdc0 |
| NTuple::Item< double > | m_tdc0 |
| NTuple::Item< double > | m_tdc1 |
| NTuple::Item< double > | m_tdc1 |
| NTuple::Item< double > | m_time1st0 |
| NTuple::Item< double > | m_time1st0 |
| NTuple::Item< double > | m_time1st1 |
| NTuple::Item< double > | m_time1st1 |
| NTuple::Item< double > | m_timeFlight |
| NTuple::Item< double > | m_timeFlight |
| NTuple::Item< double > | m_timelast0 |
| NTuple::Item< double > | m_timelast0 |
| NTuple::Item< double > | m_timelast1 |
| NTuple::Item< double > | m_timelast1 |
| NTuple::Item< double > | m_totalPhot0 |
| NTuple::Item< double > | m_totalPhot0 |
| NTuple::Item< double > | m_totalPhot1 |
| NTuple::Item< double > | m_totalPhot1 |
| NTuple::Item< double > | m_transitTime |
| NTuple::Item< double > | m_transitTime |
| NTuple::Tuple * | m_tupleTof1 |
| NTuple::Tuple * | m_tupleTof1 = 0 |
| NTuple::Tuple * | m_tupleTof2 |
| NTuple::Tuple * | m_tupleTof2 = 0 |
| NTuple::Tuple * | m_tupleTof3 |
| NTuple::Tuple * | m_tupleTof3 = 0 |
Private Attributes | |
| G4double | m_beamTime |
| G4double | m_Ce |
| G4double | m_CE |
| G4double | m_Cpe2pmt |
| G4Svc * | m_G4Svc |
| G4Svc * | m_G4Svc |
| G4double | m_HLthresh |
| G4double | m_LLthresh |
| G4double | m_noiseSigma |
| G4int | m_nPhot [m_profBinN][2] |
| G4double | m_peCorFac |
| G4double | m_phNConst |
| G4double | m_preGain |
| G4double | m_QE |
| G4double | m_rAngle |
| RealizationSvc * | m_RealizationSvc |
| RealizationSvc * | m_RealizationSvc |
| G4double | m_refIndex |
| G4double | m_scinLength |
| G4double | m_t1st [2] |
| G4double | m_tau1 |
| G4double | m_tau2 |
| G4double | m_tau3 |
| G4double | m_tauRatio |
| G4double | m_timeBinSize |
| G4double | m_tLast [2] |
| G4int | m_totalPhot [2] |
| G4double | m_ttsMean |
| G4double | m_ttsSigma |
|
|
00028 {
00029 ReadData();
00030 m_timeBinSize=0.005;
00031
00032 //retrieve G4Svc
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 //retrieve RealizationSvc
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 }
|
|
|
00089 {
00090 ;
00091 }
|
|
|
|
|
|
|
|
||||||||||||
|
|
|
||||||||||||
|
00745 {
00746 G4int ihst;
00747 ihst=G4int(endTime/m_timeBinSize);
00748 if (ihst>0 &&ihst<m_profBinN)
00749 {
00750 m_nPhot[ihst][forb]=m_nPhot[ihst][forb]+1;
00751 m_totalPhot[forb]=m_totalPhot[forb]+1;
00752 }
00753 }
00754
|
|
|
|
|
|
00378 {
00379 const G4double kappa = 0.015*cm/MeV;
00380 const G4String brMaterial = "BC408";
00381 G4double dE = hit->GetEdep();
00382 //G4cout << "The edep is "<< dE << G4endl;
00383 G4double dX = hit->GetStepL();
00384 //G4Material* materiral = hit->GetMaterial();
00385 G4double charge = hit->GetCharge();
00386 G4double cor_dE = dE;
00387 //if((materiral->GetName()==brMaterial) && charge!=0.&& dX!=0.)
00388 if(charge!=0.&& dX!=0.)
00389 {
00390 cor_dE = dE/(1+kappa*dE/dX);
00391 //if(dE>20)
00392 //{
00393 // G4cout << "\n dE > 20. Details are below:" << G4endl;
00394 // G4cout << "dE/dx:" << dE/dX << G4endl;
00395 // G4cout << "dE:" << dE << "; dX:" << dX << G4endl;
00396 // G4cout << "It is BC408. cor_dE is " << cor_dE << G4endl;
00397 // G4double ratio = cor_dE/dE;
00398 // G4cout << "The ratio cor_dE/edep is "<< ratio << G4endl;
00399 //}
00400 //G4cout << "It is BC408. cor_dE is " << cor_dE << G4endl;
00401 //G4double ratio = cor_dE/dE;
00402 //G4cout << "The ratio cor_dE/edep is "<< ratio << G4endl;
00403 }
00404 return cor_dE;
00405
00406 }
00407
|
|
||||||||||||
|
Reimplemented from BesTofDigitizerV. |
|
||||||||||||
|
Reimplemented from BesTofDigitizerV. 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 //for each digi, compute TDC and ADC
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 //cout << "*** scinNb:" << scinNb << " *** " << m_tofCaliSvc->BAtten(scinNb) << "***" << endl;
00124 //cout << "*****scinNb:"<< scinNb << " ***** A2:" << m_tofCaliSvc->BGainBackward(scinNb)
00125 // << " ***** A1:" << m_tofCaliSvc->BGainForward(scinNb) << endl;
00126
00127 TofPmtInit();
00128
00129 //fill tof Ntuple
00130 if (m_G4Svc->TofRootFlag())
00131 {
00132 if (edep>m_edepMPV)
00133 {
00134 m_partIdMPV = partId;
00135 m_scinNbMPV = scinNb;
00136 m_edepMPV = edep;
00137 }
00138 m_eTotal += edep;
00139 m_nDigi ++;
00140
00141 m_partId = partId;
00142 m_scinNb = scinNb;
00143 m_edep = edep;
00144 m_nHits = nHits;
00145 }
00146
00147 if (edep>0.01)
00148 {
00149 for (G4int j=0;j<nHits;j++)
00150 {
00151 hit= (*m_THC)[( *(scint->GetHitIndexes()) )[j]];
00152 TofPmtAccum(hit, scinNb);
00153 }
00154 if (m_G4Svc->TofRootFlag())
00155 {
00156 m_time1st0=m_t1st[0];
00157 m_time1st1=m_t1st[1];
00158 m_timelast0=m_tLast[0];
00159 m_timelast1=m_tLast[1];
00160 m_totalPhot0=m_totalPhot[0];
00161 m_totalPhot1=m_totalPhot[1];
00162 }
00163 //get final tdc and adc
00164 TofPmtRspns(scinNb);
00165 //G4cout<<"pre-cut " << partId << "\nadc0:"<<m_ADC[0]<<"; adc1:"
00166 // <<m_ADC[1]<<"; tdc0:"<<m_TDC[0]<<"; tdc1:"
00167 // <<m_TDC[1]<<G4endl;
00168
00169 G4double temp0 = m_ADC[0]+m_TDC[0];
00170 G4double temp1 = m_ADC[1]+m_TDC[1];
00171 //cout << "partid: " << partId << " temp0: " << temp0 << " temp1: " << temp1 << endl;
00172 //if ( partId==1 && m_ADC[0]>255 && m_ADC[1]>255 && m_TDC[0]>0. && m_TDC[1]>0.)
00173 if ( partId==1 && (temp0 > -1998. || temp1 > -1998.))
00174 {
00175 //const double MAX_ADC = 8191 * 0.3; // channel set up to 8192 will lead to overflow.
00176 BesTofDigi* digi = new BesTofDigi;
00177 digi->SetTrackIndex(m_trackIndex);
00178 digi->SetPartId(partId);
00179 digi->SetScinNb(scinNb);
00180 digi->SetForwADC( m_ADC[0]) ;
00181 digi->SetBackADC( m_ADC[1]) ;
00182 if (m_TDC[0]>0)
00183 m_TDC[0] = m_TDC[0]+m_beamTime;
00184 if (m_TDC[1]>0)
00185 m_TDC[1] = m_TDC[1]+m_beamTime;
00186 digi->SetForwTDC( m_TDC[0]) ;
00187 digi->SetBackTDC( m_TDC[1]) ;
00188 //G4cout<<"+++++++++++++++++++++++++++++barrel\nadc0:"<<m_ADC[0]<<"; adc1:"
00189 // <<m_ADC[1]<<"; tdc0:"<<m_TDC[0]<<"; tdc1:"
00190 // <<m_TDC[1]<<G4endl;
00191
00192 m_besTofDigitsCollection->insert(digi);
00193
00194 if (m_G4Svc->TofRootFlag())
00195 m_nDigiOut++;
00196
00197 }
00198 if (m_G4Svc->TofRootFlag())
00199 m_tupleTof1->write();
00200
00201 }
00202
00203 if (m_G4Svc->TofRootFlag())
00204 m_tupleTof2->write();
00205 }
00206 }
|
|
||||||||||||||||||||
|
|
|
||||||||||||||||||||
|
00456 {
00457 //std::cout << "BesTofDigitizerBrV2::DirectPh()" << std::endl;
00458 const static G4double silicon_oil_index = 1.43;
00459 const static G4double glass_index = 1.532;
00460 //generation photon have random direction
00461 //optical parameters of scintillation simulation
00462 G4double cos_span=1-cos( asin(silicon_oil_index/m_refIndex) );
00463 //G4double cos_spanEc = 1;
00464 G4double dcos, ran;
00465 ran=G4UniformRand();
00466 //assuming uniform phi distribution, simulate cos distr only
00467 dcos=cos_span*(ran*2.0-1.0);
00468 //G4double dcosEc = cos_spanEc*(ran*2.0-1.0);
00469 G4double r2=sqrt(emtPos.x()*emtPos.x()+emtPos.y()*emtPos.y());
00470 G4int nRef=0;
00471 G4double costheta,sintheta;
00472 G4double theta,thetaR; // thetaR is scin to air full ref angle. about 39.265 degree.
00473 costheta=dcos>0?(1-dcos):(1+dcos);
00474 theta=acos(costheta);
00475 sintheta=sin(theta);
00476 thetaR=asin(1/m_refIndex);
00477 G4double R1;
00478 R1=Reflectivity(m_refIndex,1.0,theta);
00479 G4double ratio1Mean=(CLHEP::pi)*25.5*25.5/(57.1+60.7)/25.0; //0.693657
00480 G4double ratio2Mean=(19.5/25.5)*(19.5/25.5); //0.584775
00481 G4double ratio1=G4RandGauss::shoot(ratio1Mean,0.015);
00482 G4double ratio2=G4RandGauss::shoot(ratio2Mean,0.015);
00483 G4double ratio3Mean=0.945*ratio1Mean*ratio2Mean;
00484 G4double ratio3=G4RandGauss::shoot(ratio3Mean,0.015);
00485
00486 G4double R2=1-Reflectivity(m_refIndex,silicon_oil_index, glass_index, theta);
00487 if (dcos > 0 && dcos != 1)
00488 {
00489 if (theta < thetaR) // theta < 39.265 degree
00490 {
00491 if (G4UniformRand()<ratio1) //coup1
00492 {
00493 if (G4UniformRand()<R2)
00494 {
00495 if (G4UniformRand()<ratio2) //PMT 39mm
00496 {
00497 forb=0;
00498 pathL=(m_scinLength/2-emtPos.z())/costheta;
00499 }
00500 }
00501 else
00502 {
00503 if (G4UniformRand()<ratio3)
00504 {
00505 forb=1;
00506 pathL=(1.5*m_scinLength-emtPos.z())/costheta;
00507 }
00508 }
00509 }
00510 else //Air
00511 {
00512 if (G4UniformRand()<R1)
00513 {
00514 G4double tempran=G4UniformRand();
00515 if (tempran<ratio3)
00516 {
00517 forb=1;
00518 pathL=(1.5*m_scinLength-emtPos.z())/costheta;
00519 }
00520 else if (tempran>ratio1&&G4UniformRand()<R1&&G4UniformRand()<ratio3)
00521 {
00522 forb=0;
00523 pathL=(2.5*m_scinLength-emtPos.z())/costheta;
00524 }
00525 }
00526 }
00527 }
00528 else // 39.265 <= theta < 64.832
00529 {
00530 if (G4UniformRand()<ratio1) //coup1
00531 {
00532 if (G4UniformRand()<R2)
00533 {
00534 if (G4UniformRand()<ratio2) //PMT 39mm
00535 {
00536 forb=0;
00537 pathL=(m_scinLength/2-emtPos.z())/costheta;
00538 }
00539 }
00540 else
00541 {
00542 if (G4UniformRand()<ratio3)
00543 {
00544 forb=1;
00545 pathL=(1.5*m_scinLength-emtPos.z())/costheta;
00546 }
00547 }
00548 }
00549 else //Air
00550 {
00551 G4double tempran=G4UniformRand();
00552 if (tempran<ratio3)
00553 {
00554 forb=1;
00555 pathL=(1.5*m_scinLength-emtPos.z())/costheta;
00556 }
00557 else if (tempran>ratio1 && G4UniformRand()<ratio3)
00558 {
00559 forb=0;
00560 pathL=(2.5*m_scinLength-emtPos.z())/costheta;
00561 }
00562 }
00563 }
00564 }
00565 else if (dcos < 0 && dcos != -1)
00566 {
00567 if (theta < thetaR) // theta < 39.265 degree
00568 {
00569 if (G4UniformRand()<ratio1) //coup1
00570 {
00571 if (G4UniformRand()<R2)
00572 {
00573 if (G4UniformRand()<ratio2) //PMT 39mm
00574 {
00575 forb=1;
00576 pathL=(m_scinLength/2+emtPos.z())/costheta;
00577 }
00578 }
00579 else
00580 {
00581 if (G4UniformRand()<ratio3)
00582 {
00583 forb=0;
00584 pathL=(1.5*m_scinLength+emtPos.z())/costheta;
00585 }
00586 }
00587 }
00588 else //Air
00589 {
00590 if (G4UniformRand()<R1)
00591 {
00592 G4double tempran=G4UniformRand();
00593 if (tempran<ratio3)
00594 {
00595 forb=0;
00596 pathL=(1.5*m_scinLength+emtPos.z())/costheta;
00597 }
00598 else if (tempran>ratio1&&G4UniformRand()<R1&&G4UniformRand()<ratio3)
00599 {
00600 forb=1;
00601 pathL=(2.5*m_scinLength+emtPos.z())/costheta;
00602 }
00603 }
00604 }
00605 }
00606 else // 39.265 <= theta < 64.832
00607 {
00608 if (G4UniformRand()<ratio1) //coup1
00609 {
00610 if (G4UniformRand()<R2)
00611 {
00612 if (G4UniformRand()<ratio2) //PMT 39mm
00613 {
00614 forb=1;
00615 pathL=(m_scinLength/2+emtPos.z())/costheta;
00616 }
00617 }
00618 else
00619 {
00620 if (G4UniformRand()<ratio3)
00621 {
00622 forb=0;
00623 pathL=(1.5*m_scinLength+emtPos.z())/costheta;
00624 }
00625 }
00626 }
00627 else //Air
00628 {
00629 G4double tempran=G4UniformRand();
00630 if (tempran<ratio3)
00631 {
00632 forb=0;
00633 pathL=(1.5*m_scinLength+emtPos.z())/costheta;
00634 }
00635 else if (tempran>ratio1 && G4UniformRand()<ratio3)
00636 {
00637 forb=1;
00638 pathL=(2.5*m_scinLength+emtPos.z())/costheta;
00639 }
00640 }
00641 }
00642 }
00643
00644 G4double convFactor=180./3.1415926;
00645 if (theta>asin(1)-thetaR)
00646 {
00647 G4double alpha = G4UniformRand()*asin(1.);
00648 G4int nRef1 = pathL*sintheta*cos(alpha)/50.0+0.5;
00649 G4int nRef2 = pathL*sintheta*sin(alpha)/58.9+0.5;
00650 G4double beta1=acos(sintheta*cos(alpha));
00651 G4double beta2=acos(sintheta*sin(alpha));
00652 beta2 -= 3.75/convFactor;
00653 G4double R21,R22;
00654 R21=Reflectivity(m_refIndex,1.0,beta1);
00655 R22=Reflectivity(m_refIndex,1.0,beta2);
00656 for (G4int i=0;i<nRef1;i++)
00657 {
00658 if (G4UniformRand()<(1-R21) && G4UniformRand()<0.15)
00659 pathL=-9;
00660 }
00661 for (G4int i=0;i<nRef2;i++)
00662 {
00663 if (G4UniformRand()<(1-R22) && G4UniformRand()<0.15)
00664 pathL=-9;
00665 }
00666 }
00667 }
00668
|
|
|
|
|
|
00149 {
00150 for (G4int i=0;i<2;i++)
00151 {
00152 m_ADC[i] = -999;
00153 m_TDC[i] = -999;
00154 }
00155 m_trackIndex = -999;
00156 m_globalTime = 9999;
00157 }
|
|
|
|
|
|
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 }
|
|
||||||||||||||||||||
|
|
|
||||||||||||||||
|
|
|
||||||||||||||||||||
|
00409 {
00410 G4double I1,I2,I3,rp1,rs1,rp2,rs2,Rp1,Rs1,Rp2,Rs2,Rp,Rs;
00411 G4double R=1.0;
00412 //n1=m_refIndex;
00413 //n2=1.0;
00414 I1=theta;
00415 if(I1<asin(n2/n1))
00416 {
00417 I2=asin(sin(I1)*(n1/n2));
00418 rp1=(n1/cos(I1)-n2/cos(I2))/(n1/cos(I1)+n2/cos(I2));
00419 rs1=(n1*cos(I1)-n2*cos(I2))/(n1*cos(I1)+n2*cos(I2));
00420 Rp1=rp1*rp1;
00421 Rs1=rs1*rs1;
00422
00423 I3=asin(sin(I1)*(n1/n3));
00424 rp2=(n2/cos(I2)-n3/cos(I3))/(n2/cos(I2)+n3/cos(I3));
00425 rs2=(n2*cos(I2)-n3*cos(I3))/(n2*cos(I2)+n3*cos(I3));
00426 Rp2=rp2*rp2;
00427 Rs2=rs2*rs2;
00428 Rp=(Rp1+Rp2-2*Rp1*Rp2)/(1-Rp1*Rp2);
00429 Rs=(Rs1+Rs2-2*Rs1*Rs2)/(1-Rs1*Rs2);
00430 R=(Rp+Rs)/2.;
00431 }
00432 return R;
00433 }
00434
|
|
||||||||||||||||
|
00437 {
00438 G4double I1,I2,rp,rs,Rp,Rs;
00439 G4double R=1.0;
00440 //n1=m_refIndex;
00441 //n2=1.0;
00442 I1=theta;
00443 if (I1<asin(n2/n1))
00444 {
00445 I2=asin(sin(I1)*(n1/n2));
00446 rp=(n1/cos(I1)-n2/cos(I2))/(n1/cos(I1)+n2/cos(I2));
00447 rs=(n1*cos(I1)-n2*cos(I2))/(n1*cos(I1)+n2*cos(I2));
00448 Rp=rp*rp;
00449 Rs=rs*rs;
00450 R=(Rp+Rs)/2.;
00451 }
00452 return R;
00453 }
00454
|
|
|
|
|
|
00692 {
00693 G4double tmp_tauRatio,tmp_tau1,tmp_tau2,tmp_tau3;
00694 tmp_tauRatio = m_tauRatio;
00695 tmp_tau1 = m_tau1;
00696 tmp_tau2 = m_tau2;
00697 tmp_tau3 = m_tau3;
00698 G4double UniformR = tmp_tauRatio/(1+tmp_tauRatio);
00699 G4double EmissionTime;
00700 if (G4UniformRand()>UniformR) {
00701 while (1) {
00702 EmissionTime = -tmp_tau2*log( G4UniformRand() );
00703 if (G4UniformRand()-exp(EmissionTime/tmp_tau2-EmissionTime/tmp_tau1)>1.E-8)
00704 break;
00705 }
00706 }
00707 else EmissionTime = -tmp_tau3*log( G4UniformRand() );
00708 return EmissionTime;
00709 }
00710
|
|
||||||||||||
|
|
|
||||||||||||
|
00245 {
00246 BesTofGeoParameter* tofPara = BesTofGeoParameter::GetInstance();
00247 //std::cout << "BesTofDigitizerBrV2::TofPmtAccum()" << std::endl;
00248 G4double cvelScint=c_light/m_refIndex/1.07;
00249 //Get information of this step
00250 G4ThreeVector pos=hit->GetPos();
00251 G4int trackIndex = hit->GetTrackIndex();
00252 G4int partId =hit->GetPartId();
00253 G4double edep=hit->GetEdep();
00254 G4double stepL=hit->GetStepL();
00255 G4double deltaT=hit->GetDeltaT();
00256 G4double timeFlight=hit->GetTime()-m_beamTime;
00257 //std::cout << "timeFlight: " << timeFlight << std::endl;
00258 if (timeFlight < m_globalTime)
00259 {
00260 m_globalTime = timeFlight;
00261 m_trackIndex = trackIndex;
00262 }
00263
00264 G4ThreeVector pDirection=hit->GetPDirection();
00265 G4double nx=pDirection.x();
00266 G4double ny=pDirection.y();
00267 G4double nz=pDirection.z();
00268
00269 //phNConst=(Nph/MeV)*(QE)*(CE)*(1-cos(crit));
00270 // =10000 * 0.2 * 0.6 * (1-cos(39.265))=270.931
00271 //asin(1/1.58) = 39.265248
00272 //only the light has theta=0---39.265 can go out to PMT, the probability is computed with solid angle
00273 //solid angle = phi*(1-cos(theta)), phi=2pi
00274
00275 //Cpe2pmt=cathode area/scint area
00276 //peCorFac=correction factor for eff. of available Npe
00277
00278 //G4double thetaProb = 1-sqrt(m_refIndex*m_refIndex-1)/m_refIndex;
00279 G4double thetaProb=1-cos( asin(1.43/m_refIndex));
00280 //G4double thetaProbEc = 1-1/m_refIndex;
00281
00282 //number of photons generated in this step
00283 G4double nMean, nPhoton;
00284 //std::cout << "0 BirksLaw(): " << std::endl;
00285 nMean = m_phNConst*BirksLaw(hit);
00286 //std::cout << "1 BirksLaw(): " << std::endl;
00287
00288 if(nMean>10)
00289 {
00290 G4double resolutionScale=1.;
00291 G4double sigma=resolutionScale*sqrt(nMean);
00292 nPhoton=G4int(G4RandGauss::shoot(nMean,sigma)+0.5);
00293 }
00294 else
00295 nPhoton=G4int(G4Poisson(nMean));
00296 //G4cout<<"nPhoton:"<<nPhoton<<G4endl;
00297
00298 G4int NphStep;
00299 if (partId==1)
00300 NphStep=G4int(nPhoton*thetaProb*m_rAngle*m_QE*m_CE*m_peCorFac);
00301 //else
00302 //NphStep=G4int(edep*m_phNConstEc*m_Cpe2pmtEc*0.66*m_QEEc*m_CE*m_peCorFac);
00303 //introduce poission distribution of Npe
00304 //G4double Navr = NphStep;
00305 //G4double NpePoiss = G4double(G4Poisson(Navr));
00306 //G4double adcW_corr =5.0;
00307 //NphStep=G4int( (NpePoiss - Navr)*adcW_corr + Navr );
00308
00309 if (m_G4Svc->TofRootFlag())
00310 m_NphAllSteps += NphStep;
00311 //std::cout << "m_G4Svc->TofRootFlag(): " << m_G4Svc->TofRootFlag() << std::endl;
00312
00313 G4double ddS, ddT;
00314 if (NphStep>0)
00315 {
00316 for (G4int i=0;i<NphStep;i++)
00317 {
00318 //uniform scintilation in each step
00319 ddS=stepL*G4UniformRand();
00320 ddT=deltaT*G4UniformRand();
00321 G4ThreeVector emtPos;
00322 emtPos.setX(pos.x() + nx*ddS);
00323 emtPos.setY(pos.y() + ny*ddS);
00324 emtPos.setZ(pos.z() + nz*ddS);
00325
00326 //check scinillation light whether to hit the pmt or not
00327 //forb=0/1 for forward(z>0, east) and backward(z<0, west)
00328 G4double pathL=0;
00329 G4int forb;
00330 DirectPh(partId, emtPos, pathL, forb);
00331
00332
00333 //check if photon can reach PMT or not, after attenuation
00334
00335 G4double ran = G4UniformRand();
00336 //G4double attenL = tofPara->GetAtten(scinNb);
00337 //G4double attenL = 10.*(m_tofCaliSvc->BAtten(scinNb))/0.75; // cm into mm
00338 G4double attenL = m_tofSimSvc->BarAttenLength(scinNb);
00339 attenL = 10.*attenL/0.75; // cm into mm
00340
00341 if (pathL>0 && exp(-pathL/attenL) > ran)
00342 {
00343 //propagation time in scintillator
00344 G4double scinSwim=pathL/cvelScint;
00345 //scintillation timing
00346 //G4double scinTime=GenPhoton(partId);
00347 G4double scinTime=Scintillation(partId);
00348
00349 //PMT transit time
00350 G4double transitTime=TransitTime();
00351 //sum up all time components
00352 G4double endTime= timeFlight + ddT + scinSwim + scinTime + transitTime;
00353 //std::cout << "endtime: " << endTime << std::endl;
00354
00355 if (m_G4Svc->TofRootFlag())
00356 {
00357 //m_forb = forb;
00358 //m_timeFlight = timeFlight;
00359 //m_ddT = ddT;
00360 m_scinSwim = scinSwim;
00361 //m_scinTime = scinTime;
00362 //m_transitTime = transitTime;
00363 m_endTime = endTime;
00364 m_tupleTof3->write();
00365 }
00366 //store timing into binned buffer
00367 AccuSignal(endTime, forb);
00368
00369 //update 1st and last timings here
00370 if (m_t1st[forb]>endTime) m_t1st[forb] = endTime;
00371 if (m_tLast[forb]<endTime) m_tLast[forb]= endTime;
00372 }
00373 }
00374 }
00375 }
00376
|
|
|
|
|
|
00209 {
00210 Initialize();
00211
00212 m_t1st[0]=100;
00213 m_t1st[1]=100;
00214 m_tLast[0]=0.;
00215 m_tLast[1]=0;
00216 m_totalPhot[0]=0;
00217 m_totalPhot[1]=0;
00218 for (G4int i=0;i<2;i++)
00219 for (G4int j=0;j<m_profBinN;j++)
00220 m_nPhot[j][i]=0;
00221
00222 if (m_G4Svc->TofRootFlag())
00223 {
00224 m_partId = -9;
00225 m_scinNb = -9;
00226 m_edep = 0;
00227 m_nHits = 0;
00228 m_time1st0 = 100;
00229 m_time1st1 = 100;
00230 m_timelast0 = 0;
00231 m_timelast1 = 0;
00232 m_totalPhot0 = 0;
00233 m_totalPhot1 = 0;
00234 m_NphAllSteps = 0;
00235 m_max0 = 0;
00236 m_max1 = 0;
00237 m_tdc0 = -999;
00238 m_adc0 = -999;
00239 m_tdc1 = -999;
00240 m_adc1 = -999;
00241 }
00242 }
|
|
|
|
|
|
00756 {
00757 //std::cout << "BesTofDigitizerBrV2::TofPmtRspns()" << std::endl;
00758 BesTofGeoParameter* tofPara = BesTofGeoParameter::GetInstance();
00759 //to generate PMT signal shape for single photoelectron.
00760 //only one time for a job.
00761 G4double snpe[m_snpeBinN][2];
00762
00763 //Model: f(t)=Gain*mv_1pe* t**2 * exp-(t/tau)**2/normal-const
00764 //normalization const =sqrt(pi)*tau*tau*tau/4.0
00765 //G4double tau = m_riseTime;
00766 G4double norma_const;
00767 G4double echarge=1.6e-7; //in unit of PC
00768
00769 //time profile of pmt signals for Back and Forward PMT.
00770 G4double profPmt[m_profBinN][2];
00771
00772 G4double t;
00773 G4double tau;
00774 G4int n1, n2, ii;
00775 G4int phtn;
00776 // forb = 0, east
00777 for (G4int i=0;i<m_snpeBinN;i++)
00778 {
00779 tau = tofPara->GetBrERiseTime(scinNb);
00780 norma_const = sqrt(CLHEP::pi)*tau*tau*tau/4.0; //in unit of ns**3
00781 t = (i+1)*m_timeBinSize;
00782 snpe[i][0] = m_Ce*t*t*exp(- (t/tau) * (t/tau) )/norma_const;// Pulse of one single photoelectron
00783 }
00784 // forb = 1, west
00785 for (G4int i=0;i<m_snpeBinN;i++)
00786 {
00787 tau = tofPara->GetBrWRiseTime(scinNb);
00788 norma_const = sqrt(CLHEP::pi)*tau*tau*tau/4.0; //in unit of ns**3
00789 t = (i+1)*m_timeBinSize;
00790 snpe[i][1] = m_Ce*t*t*exp(- (t/tau) * (t/tau) )/norma_const;
00791 }
00792 //for barrel and endcap tof: fb=2 or 1
00793 G4int fb=2;
00794 G4double Npoisson;
00795 G4double pmtGain0,pmtGain,relativeGain,smearPMTgain;
00796 G4double smearADC[2] = {0.,0.};
00797 G4int runId = m_RealizationSvc->getRunId();
00798 pmtGain0 = m_tofSimSvc->BarPMTGain();
00799
00800 // if(runId>=-80000 && runId<=-9484)
00801 // {
00802 // // After TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 5.E5
00803 // // High/Low Threshold for barrel: 500/125 mV
00804 // pmtGain0 = 6.E5;
00805 // }
00806 // else
00807 // {
00808 // // Before TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 2.5E5
00809 // // High/Low Threshold for barrel: 600/150 mV
00810 // pmtGain0 = 5.E5;
00811 // }
00812
00813 G4double timeSmear = G4RandGauss::shoot(0,0.045);
00814
00815 for (G4int j=0; j<fb; j++)
00816 {
00817 if (m_totalPhot[j] > 0)
00818 {
00819 n1=G4int(m_t1st[j]/m_timeBinSize);
00820 n2=G4int(m_tLast[j]/m_timeBinSize);
00821 //std::cout << "n1: " << n1 << std::endl;
00822
00823 for (G4int i=0;i<m_profBinN;i++)
00824 profPmt[i][j]=0.0;
00825
00826 //generate PMT pulse
00827 n2 = n2<m_profBinN ? n2:m_profBinN;
00828 for (G4int i=n1;i<n2;i++)
00829 {
00830 phtn=m_nPhot[i][j];
00831 if (phtn>0)
00832 {
00833 while(1) {
00834 Npoisson = G4Poisson(10.0);
00835 if(Npoisson>0.) break;
00836 }
00837 while(1) {
00838 //pmtGain = j ? tofPara->GetBrWPMTgain(scinNb) : tofPara->GetBrEPMTgain(scinNb);
00839 //relativeGain = j ? m_tofCaliSvc->BGainBackward(scinNb) : m_tofCaliSvc->BGainForward(scinNb);
00840 relativeGain = j ? m_tofSimSvc->BarGain2(scinNb) : m_tofSimSvc->BarGain1(scinNb);
00841 pmtGain = pmtGain0*relativeGain;
00842 smearPMTgain = G4RandGauss::shoot(pmtGain,pmtGain/sqrt(Npoisson));
00843 //smearPMTgain = pmtGain;
00844 if(smearPMTgain>0) break;
00845 }
00846 smearADC[j] += phtn*smearPMTgain;
00847
00848 for (G4int ihst=0; ihst<m_snpeBinN; ihst++)
00849 {
00850 ii=i+ihst;
00851 if (ii<m_profBinN)
00852 profPmt[ii][j] += smearPMTgain*phtn*snpe[ihst][j];
00853 else
00854 break;
00855 }
00856 }
00857 }
00858
00859 //add preamplifier and noise
00860 for (int i=0;i<m_profBinN;i++)
00861 {
00862 if (profPmt[i][j]>0)
00863 profPmt[i][j] = m_preGain*profPmt[i][j]+G4RandGauss::shoot(0,m_noiseSigma);
00864 }
00865
00866 //get pulse height
00867 G4double max=0;
00868 for (int i=n1;i<m_profBinN;i++)
00869 {
00870 if (profPmt[i][j]>max)
00871 max=profPmt[i][j];
00872 }
00873 if (m_G4Svc->TofRootFlag())
00874 {
00875 if (j==0) m_max0=max;
00876 else m_max1=max;
00877 }
00878
00879 G4double tmp_HLthresh, tmp_LLthresh, adcFactor;
00880 G4double ratio;
00881
00882 if(runId>=-80000 && runId<=-9484)
00883 {
00884 // After TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 5.E5
00885 // High/Low Threshold for barrel: 500/125 mV
00886 tmp_HLthresh = m_HLthresh;
00887 tmp_LLthresh = m_LLthresh;
00888 adcFactor = 5.89;
00889 }
00890 else
00891 {
00892 // Before TOF HV adjustment, PMT Gain of barrel TOF in MC is set to 2.5E5
00893 // High/Low Threshold for barrel: 600/150 mV
00894 tmp_HLthresh = 600;
00895 tmp_LLthresh = 150;
00896 adcFactor = 4.8;
00897 }
00898 // if(runId>=-80000 && runId<=-9484)
00899 // {
00900 // ratio=2.16*1923.8/2197.8;
00901 // }
00902 // else
00903 // {
00904 // ratio = 2.11*2437.0/3102.9;
00905 // }
00906 tmp_HLthresh = m_tofSimSvc->BarHighThres();
00907 tmp_LLthresh = m_tofSimSvc->BarLowThres();
00908 ratio = m_tofSimSvc->BarConstant();
00909
00910 //get final tdc and adc
00911 if (max>=tmp_HLthresh)
00912 {
00913 for (int i=0;i<m_profBinN;i++)
00914 {
00915 if ( profPmt[i][j] >= tmp_LLthresh )
00916 {
00917 m_TDC[j] = i*m_timeBinSize + timeSmear +G4RandGauss::shoot(0,0.025); // Adding Electronical Uncertainty of 25ps
00918
00919 if( m_G4Svc->TofSaturationFlag())
00920 {
00921 //get ADC[j] using tofQElecSvc
00922 double x = m_preGain*smearADC[j]*echarge*ratio;
00923 if (j==0)
00924 {
00925 m_ADC[j] = m_tofQElecSvc->BQChannel1(scinNb,x);
00926 }
00927 else
00928 {
00929 m_ADC[j] = m_tofQElecSvc->BQChannel2(scinNb,x);
00930 }
00931 }
00932 else
00933 m_ADC[j] = m_preGain*smearADC[j]*echarge*adcFactor;
00934
00935
00936 if (m_G4Svc->TofRootFlag())
00937 {
00938 if (j==0) {
00939 m_tdc0 = m_TDC[0];
00940 m_adc0 = m_ADC[0];
00941 }
00942 else {
00943 m_tdc1 = m_TDC[1];
00944 m_adc1 = m_ADC[1];
00945 }
00946 }
00947 break;
00948 }
00949 }
00950 }
00951 }
00952 }
00953 }
00954 }
|
|
|
|
|
|
00739 {
00740 //get time transit spread
00741 return G4RandGauss::shoot(m_ttsMean,m_ttsSigma);
00742 }
00743
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Reimplemented from BesTofDigitizerV. |
|
|
Reimplemented from BesTofDigitizerV. |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
1.3.9.1