00001 //---------------------------------------------------------------------------// 00002 // BOOST --- BESIII Object_Oriented Simulation Tool // 00003 //---------------------------------------------------------------------------// 00004 //Description: 00005 //Author: Hemiao 00006 //Created: Oct 25, 2004 00007 //Modified: 00008 //Comment: 00009 //---------------------------------------------------------------------------// 00010 // $Id: BesEmcWaveform.cc 00011 00012 #include "BesEmcWaveform.hh" 00013 #include "Randomize.hh" 00014 #include "CLHEP/Random/RanecuEngine.h" 00015 #include "CLHEP/Random/RandGauss.h" 00016 #include "BesEmcParameter.hh" 00017 #include "BesEmcDigi.hh" 00018 using namespace CLHEP; 00019 // Constructors 00020 00021 BesEmcWaveform::BesEmcWaveform() 00022 { 00023 BesEmcParameter& emcPara=BesEmcParameter::GetInstance(); 00024 //BesEmcParameter emcPara; 00025 //emcPara.ReadData(); 00026 array_size = emcPara.GetArraySize(); 00027 m_tau = emcPara.GetTau(); 00028 m_highRange = emcPara.GetHighRange(); 00029 m_midRange = emcPara.GetMidRange(); 00030 m_lowRange = emcPara.GetLowRange(); 00031 m_sampleTime = emcPara.GetSampleTime(); 00032 m_peakTime = emcPara.GetPeakTime(); 00033 m_timeOffset = emcPara.GetTimeOffset(); 00034 m_bitNb = emcPara.GetADCbit(); 00035 m_photonsPerMeV = emcPara.GetPhotonsPerMeV(); 00036 m_nonuniformity = emcPara.GetNonuniformity(); 00037 m_flag = -1; 00038 m_highPrecision = m_highRange/((G4double)(1<<m_bitNb)); 00039 m_midPrecision = m_midRange/((G4double)(1<<m_bitNb)); 00040 m_lowPrecision = m_lowRange/((G4double)(1<<m_bitNb)); 00041 emcWave=new G4double[array_size]; 00042 00043 for (G4long i=0; i<array_size;i++) 00044 emcWave[i]=0; 00045 } 00046 00047 BesEmcWaveform::BesEmcWaveform(G4long size, G4double tau, G4double sampleTime) 00048 :m_tau(tau),m_sampleTime(sampleTime) 00049 { 00050 if(size>0){ 00051 array_size=size; 00052 emcWave=new G4double[array_size]; 00053 G4double *init=emcWave+array_size; 00054 while(init!=emcWave) *--init=0.0; 00055 } 00056 else{ 00057 G4Exception("Invalid size"); 00058 } 00059 } 00060 00061 // Destructors 00062 BesEmcWaveform::~BesEmcWaveform(){ 00063 delete []emcWave; 00064 emcWave=0; 00065 } 00066 00067 // Operators 00068 BesEmcWaveform &BesEmcWaveform::operator*=(const G4double &scale) 00069 { 00070 for (G4long i=0; i<array_size;i++) emcWave[i]*=scale; 00071 return *this; 00072 } 00073 00074 BesEmcWaveform &BesEmcWaveform::operator/=(const G4double &scale) 00075 { 00076 for (G4long i=0; i<array_size;i++) emcWave[i]/=scale; 00077 return *this; 00078 } 00079 00080 BesEmcWaveform &BesEmcWaveform::operator+=(const BesEmcWaveform &assign) 00081 { 00082 for (G4long i=0; i<array_size;i++) emcWave[i]+=assign[i]; 00083 return *this; 00084 } 00085 00086 BesEmcWaveform &BesEmcWaveform::operator=(const BesEmcWaveform &assign) 00087 { 00088 if (this != &assign) { 00089 if (emcWave!=0) delete []emcWave; 00090 emcWave=new G4double[assign.array_size]; 00091 array_size=assign.array_size; 00092 for (G4long i=0;i<array_size;i++) emcWave[i]=assign[i]; 00093 } 00094 return *this; 00095 } 00096 00097 G4double BesEmcWaveform::max(G4long &binOfMax) const 00098 { 00099 G4double maxi=emcWave[0]; 00100 binOfMax = 0; 00101 for (G4long i=1;i<array_size;i++) 00102 { 00103 if (emcWave[i]>maxi) 00104 { 00105 maxi=emcWave[i]; 00106 binOfMax = i; 00107 } 00108 } 00109 return maxi; 00110 } 00111 00112 void BesEmcWaveform::updateWaveform(BesEmcHit* hit) 00113 { 00114 G4double energy = hit->GetEdepCrystal(); 00115 G4double time = hit->GetTimeCrystal(); 00116 00117 makeWaveform(energy, time); 00118 } 00119 00120 void BesEmcWaveform::updateWaveform(BesEmcDigi* digi) 00121 { 00122 G4double energy = digi->GetEnergy(); 00123 G4double time = digi->GetTime()*m_sampleTime+m_timeOffset-m_peakTime; 00124 00125 makeWaveform(energy, time); 00126 } 00127 00128 void BesEmcWaveform::makeWaveform(G4double energy, G4double time) 00129 { 00130 G4double amplitude; 00131 if(m_photonsPerMeV==0) 00132 { 00133 amplitude = energy; 00134 } 00135 else 00136 { 00137 G4double photons = energy*m_photonsPerMeV; 00138 if(photons>0) { 00139 G4double photonStatFactor = RandGauss::shoot(photons, sqrt(photons))/photons; 00140 amplitude = energy*photonStatFactor; 00141 } else { 00142 amplitude = 0; 00143 } 00144 } 00145 00146 G4double tempTime; 00147 tempTime = 0; 00148 00149 G4double peak = m_peakTime*m_peakTime*m_peakTime*m_peakTime*exp(-m_peakTime/m_tau)/24; 00150 00151 for (G4long i=0;i<array_size;i++) 00152 { 00153 tempTime = i*m_sampleTime + m_timeOffset - time; 00154 if(tempTime>0) 00155 emcWave[i] += amplitude*tempTime*tempTime*tempTime*tempTime*exp(-tempTime/m_tau)/(24*peak); 00156 } 00157 } 00158 00159 void BesEmcWaveform::digitize() 00160 { 00161 G4double oneBitResolution; 00162 oneBitResolution = m_midRange*2/((G4double)(1<<m_bitNb)); //not used now 00163 G4double energy; 00164 00165 for(G4long i=0;i<array_size;i++) 00166 { 00167 energy = emcWave[i]; 00168 if(energy > m_highRange) 00169 G4cout<<"---In BesEmcWaveform: Over measurement!--- energy="<<energy<<G4endl; 00170 else if(energy > m_midRange) 00171 { 00172 m_flag = 2; 00173 emcWave[i] = (G4double)((G4long)(emcWave[i]/m_highPrecision))*m_highPrecision; 00174 } 00175 else if(energy > m_lowRange) 00176 { 00177 m_flag = 1; 00178 emcWave[i] = (G4double)((G4long)(emcWave[i]/m_midPrecision))*m_midPrecision; 00179 } 00180 else 00181 { 00182 m_flag = 0; 00183 emcWave[i] = (G4double)((G4long)(emcWave[i]/m_lowPrecision))*m_lowPrecision; 00184 } 00185 } 00186 } 00187 00188 void BesEmcWaveform::print() 00189 { 00190 G4cout<<"New Wave!"<<G4endl; 00191 for(G4long i=0;i<array_size;i++) 00192 G4cout<<emcWave[i]<<"\t"; 00193 G4cout<<G4endl; 00194 } 00195 00196 void BesEmcWaveform::addElecNoise(G4double width, G4double coherentNoise) 00197 { 00198 for(G4int i=0;i<array_size;i++) 00199 { 00200 emcWave[i] += RandGauss::shoot()*width; 00201 emcWave[i] += coherentNoise; 00202 emcWave[i] = emcWave[i]>0?emcWave[i]:0; 00203 } 00204 }