|
|
|
#include <SimG4CMS/Forward/interface/CastorSD.h>
Public Member Functions | |
| CastorSD (G4String, const DDCompactView &, SensitiveDetectorCatalog &clg, edm::ParameterSet const &, const SimTrackManager *) | |
| virtual double | getEnergyDeposit (G4Step *) |
| virtual uint32_t | setDetUnitId (G4Step *step) |
| void | setNumberingScheme (CastorNumberingScheme *scheme) |
| virtual | ~CastorSD () |
Private Attributes | |
| CastorNumberingScheme * | numberingScheme |
Usage: Used in sensitive detector builder
Definition at line 29 of file CastorSD.h.
| CastorSD::CastorSD | ( | G4String | name, | |
| const DDCompactView & | cpv, | |||
| SensitiveDetectorCatalog & | clg, | |||
| edm::ParameterSet const & | p, | |||
| const SimTrackManager * | manager | |||
| ) |
Definition at line 28 of file CastorSD.cc.
References edm::ParameterSet::getParameter(), and setNumberingScheme().
00031 : 00032 CaloSD(name, cpv, clg, p, manager), numberingScheme(0) { 00033 00034 edm::ParameterSet m_CastorSD = p.getParameter<edm::ParameterSet>("CastorSD"); 00035 00036 setNumberingScheme(new CastorNumberingScheme()); 00037 edm::LogInfo("ForwardSim") 00038 << "***************************************************\n" 00039 << "* *\n" 00040 << "* Constructing a CastorSD with name " << GetName() << "\n" 00041 << "* *\n" 00042 << "***************************************************"; 00043 }
| CastorSD::~CastorSD | ( | ) | [virtual] |
| double CastorSD::getEnergyDeposit | ( | G4Step * | aStep | ) | [virtual] |
Reimplemented from CaloSD.
Definition at line 47 of file CastorSD.cc.
References a, DeDxTools::beta(), d, eta, Exception, edm::Service< T >::isAvailable(), funct::log(), LogDebug, max, min, SensitiveDetector::name, NULL, phi, pi, CaloSD::preStepPoint, r, funct::sqrt(), funct::tan(), theta, and CaloSD::theTrack.
00047 { 00048 00049 float NCherPhot = 0.; 00050 00051 if (aStep == NULL) { 00052 return 0; 00053 } else { 00054 // preStepPoint information ********************************************* 00055 00056 G4SteppingControl stepControlFlag = aStep->GetControlFlag(); 00057 G4StepPoint* preStepPoint = aStep->GetPreStepPoint(); 00058 G4VPhysicalVolume* currentPV = preStepPoint->GetPhysicalVolume(); 00059 G4String name = currentPV->GetName(); 00060 std::string nameVolume; 00061 nameVolume.assign(name,0,4); 00062 00063 00064 00065 G4ThreeVector hitPoint = preStepPoint->GetPosition(); 00066 G4ThreeVector hit_mom = preStepPoint->GetMomentumDirection(); 00067 G4double stepl = aStep->GetStepLength()/cm; 00068 G4double beta = preStepPoint->GetBeta(); 00069 G4double charge = preStepPoint->GetCharge(); 00070 // G4VProcess* curprocess = preStepPoint->GetProcessDefinedStep(); 00071 // G4String namePr = preStepPoint->GetProcessDefinedStep()->GetProcessName(); 00072 // std::string nameProcess; 00073 // nameProcess.assign(namePr,0,4); 00074 00075 // G4LogicalVolume* lv = currentPV->GetLogicalVolume(); 00076 // G4Material* mat = lv->GetMaterial(); 00077 // G4double rad = mat->GetRadlen(); 00078 00079 00080 // postStepPoint information ********************************************* 00081 G4StepPoint* postStepPoint= aStep->GetPostStepPoint(); 00082 G4VPhysicalVolume* postPV = postStepPoint->GetPhysicalVolume(); 00083 G4String postname = postPV->GetName(); 00084 std::string postnameVolume; 00085 postnameVolume.assign(postname,0,4); 00086 00087 00088 // theTrack information ************************************************* 00089 G4Track* theTrack = aStep->GetTrack(); 00090 G4double entot = theTrack->GetTotalEnergy(); 00091 G4ThreeVector vert_mom = theTrack->GetVertexMomentumDirection(); 00092 00093 G4ThreeVector localPoint = theTrack->GetTouchable()->GetHistory()-> 00094 GetTopTransform().TransformPoint(hitPoint); 00095 00096 00097 // calculations... ************************************************* 00098 00099 float costheta =vert_mom.z()/sqrt(vert_mom.x()*vert_mom.x()+ 00100 vert_mom.y()*vert_mom.y()+ 00101 vert_mom.z()*vert_mom.z()); 00102 float theta = acos(std::min(std::max(costheta,float(-1.)),float(1.))); 00103 float eta = -log(tan(theta/2)); 00104 float phi = -100.; 00105 if (vert_mom.x() != 0) phi = atan2(vert_mom.y(),vert_mom.x()); 00106 if (phi < 0.) phi += twopi; 00107 G4String particleType = theTrack->GetDefinition()->GetParticleName(); 00108 G4int primaryID = theTrack->GetTrackID(); 00109 // ************************************************* 00110 00111 00112 // ************************************************* 00113 double edep = aStep->GetTotalEnergyDeposit(); 00114 00115 // ************************************************* 00116 00117 00118 // ************************************************* 00119 // take into account light collection curve for plate 00120 // double weight = curve_Castor(nameVolume, preStepPoint); 00121 // double edep = aStep->GetTotalEnergyDeposit() * weight; 00122 // ************************************************* 00123 00124 00125 // ************************************************* 00126 /* comments for sensitive volumes: 00127 C001 ...-... CP06,CBU1 ...-...CALM --- > fibres and bundle 00128 for first release of CASTOR 00129 CASF --- > quartz plate for first and second releases of CASTOR 00130 GF2Q, GFNQ, GR2Q, GRNQ 00131 for tests with my own test geometry of HF (on ask of Gavrilov) 00132 C3TF, C4TF - for third release of CASTOR 00133 */ 00134 double meanNCherPhot; 00135 00136 if(nameVolume == "C3EF" || nameVolume == "C4EF" || nameVolume == "C3HF" || 00137 nameVolume == "C4HF") { 00138 00139 float bThreshold = 0.67; 00140 float nMedium = 1.4925; 00141 // float photEnSpectrDL = (1./400.nm-1./700.nm)*10000000.cm/nm; /* cm-1 */ 00142 // float photEnSpectrDL = 10714.285714; 00143 00144 float photEnSpectrDE = 1.24; /* see below */ 00145 /* E = 2pi*(1./137.)*(eV*cm/370.)/lambda = */ 00146 /* = 12.389184*(eV*cm)/lambda */ 00147 /* Emax = 12.389184*(eV*cm)/400nm*10-7cm/nm = 3.01 eV */ 00148 /* Emin = 12.389184*(eV*cm)/700nm*10-7cm/nm = 1.77 eV */ 00149 /* delE = Emax - Emin = 1.24 eV --> */ 00150 /* */ 00151 /* default for Castor nameVolume == "CASF" or (C3TF & C4TF) */ 00152 float effPMTandTransport = 0.15; 00153 /* for test HF geometry volumes: */ 00154 if(nameVolume == "GF2Q" || nameVolume == "GFNQ" || 00155 nameVolume == "GR2Q" || nameVolume == "GRNQ") 00156 effPMTandTransport = 0.15; 00157 00158 float thFullRefl = 23.; /* 23.dergee */ 00159 float thFullReflRad = thFullRefl*pi/180.; 00160 00161 /* default for Castor nameVolume == "CASF" or (C3TF & C4TF) */ 00162 float thFibDir = 45.; /* .dergee */ 00163 /* for test HF geometry volumes: */ 00164 if(nameVolume == "GF2Q" || nameVolume == "GFNQ" || 00165 nameVolume == "GR2Q" || nameVolume == "GRNQ") 00166 thFibDir = 0.0; /* .dergee */ 00167 /* */ 00168 float thFibDirRad = thFibDir*pi/180.; 00169 /* */ 00170 /* */ 00171 00172 // at which theta the point is located: 00173 // float th1 = hitPoint.theta(); 00174 00175 // theta of charged particle in LabRF(hit momentum direction): 00176 float costh =hit_mom.z()/sqrt(hit_mom.x()*hit_mom.x()+ 00177 hit_mom.y()*hit_mom.y()+ 00178 hit_mom.z()*hit_mom.z()); 00179 float th = acos(std::min(std::max(costh,float(-1.)),float(1.))); 00180 00181 // just in case (can do bot use): 00182 if (th < 0.) th += twopi; 00183 00184 float thgrad = th*180./pi; 00185 00186 00187 // theta of cone with Cherenkov photons w.r.t.direction of charged part.: 00188 float costhcher =1./(nMedium*beta); 00189 float thcher = acos(std::min(std::max(costhcher,float(-1.)),float(1.))); 00190 float thchergrad = thcher*180./pi; 00191 00192 // diff thetas of charged part. and quartz direction in LabRF: 00193 float DelFibPart = fabs(th - thFibDirRad); 00194 float DelFibPartgrad = DelFibPart*180./pi; 00195 00196 00197 00198 // define real distances: 00199 float d = fabs(tan(th)-tan(thFibDirRad)); 00200 00201 // float a = fabs(tan(thFibDirRad)-tan(thFibDirRad+thFullReflRad)); 00202 // float r = fabs(tan(th)-tan(th+thcher)); 00203 00204 float a = tan(thFibDirRad)+tan(fabs(thFibDirRad-thFullReflRad)); 00205 float r = tan(th)+tan(fabs(th-thcher)); 00206 00207 00208 // define losses d_qz in cone of full reflection inside quartz direction 00209 float d_qz; 00210 float variant; 00211 00212 if(DelFibPart > (thFullReflRad + thcher) ) {d_qz = 0.; variant=0.;} 00213 // if(d > (r+a) ) {d_qz = 0.; variant=0.;} 00214 else { 00215 if((th + thcher) < (thFibDirRad+thFullReflRad) && 00216 (th - thcher) > (thFibDirRad-thFullReflRad) 00217 ) {d_qz = 1.; variant=1.;} 00218 // if((DelFibPart + thcher) < thFullReflRad ) {d_qz = 1.; variant=1.;} 00219 // if((d+r) < a ) {d_qz = 1.; variant=1.;} 00220 else { 00221 if((thFibDirRad + thFullReflRad) < (th + thcher) && 00222 (thFibDirRad - thFullReflRad) > (th - thcher) ) { 00223 // if((thcher - DelFibPart ) > thFullReflRad ) 00224 // if((r-d ) > a ) 00225 d_qz = 0.; variant=2.; 00226 00227 } else { 00228 // if((thcher + DelFibPart ) > thFullReflRad && 00229 // thcher < (DelFibPart+thFullReflRad) ) 00230 // { 00231 d_qz = 0.; variant=3.; 00232 00233 00234 // use crossed length of circles(cone projection) 00235 // dC1/dC2 : 00236 float arg_arcos = 0.; 00237 float tan_arcos = 2.*a*d; 00238 if(tan_arcos != 0.) arg_arcos =(r*r-a*a-d*d)/tan_arcos; 00239 arg_arcos = fabs(arg_arcos); 00240 float th_arcos = acos(std::min(std::max(arg_arcos,float(-1.)),float(1.))); 00241 d_qz = th_arcos/pi/2.; 00242 d_qz = fabs(d_qz); 00243 00244 00245 00246 // } 00247 // else 00248 // { 00249 // d_qz = 0.; variant=4.; 00250 //#ifdef debug 00251 // std::cout <<" ===============>variant 4 information: <===== " <<std::endl; 00252 // std::cout <<" !!!!!!!!!!!!!!!!!!!!!! variant = " << variant <<std::endl; 00253 //#endif 00254 // 00255 // } 00256 } 00257 } 00258 } 00259 00260 00261 // !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 00262 00263 if(charge != 0. && beta > bThreshold && d_qz != 0. ) { 00264 00265 meanNCherPhot = 370.*charge*charge* 00266 ( 1. - 1./(nMedium*nMedium*beta*beta) )* 00267 photEnSpectrDE*stepl; 00268 00269 // dLamdX: 00270 // meanNCherPhot = (2.*pi/137.)*charge*charge* 00271 // ( 1. - 1./(nMedium*nMedium*beta*beta) )* 00272 // photEnSpectrDL*stepl; 00273 00274 00275 // NCherPhot = meanNCherPhot; 00276 // Poisson: 00277 edm::Service<edm::RandomNumberGenerator> rng; 00278 if ( ! rng.isAvailable()) { 00279 throw cms::Exception("Configuration") 00280 << "ZdcSD requires the RandomNumberGeneratorService\n" 00281 << "which is not present in the configuration file. " 00282 << "You must add the service\n" << "in the configuration file " 00283 << "or remove the modules that require it."; 00284 } 00285 CLHEP::RandPoissonQ randPoisson(rng->getEngine()); 00286 G4int poissNCherPhot = (G4int) randPoisson.fire(meanNCherPhot); 00287 00288 // G4int poissNCherPhot = (G4int) G4Poisson(meanNCherPhot); 00289 00290 if(poissNCherPhot < 0) poissNCherPhot = 0; 00291 00292 NCherPhot = poissNCherPhot * effPMTandTransport * d_qz; 00293 00294 00295 00296 #ifdef debug 00297 LogDebug("ForwardSim") << " ==============================> start all " 00298 << "information:<========= \n" << " =====> for " 00299 << "test:<=== \n" << " variant = " << variant 00300 << "\n thgrad = " << thgrad << "\n thchergrad " 00301 << "= " << thchergrad << "\n DelFibPartgrad = " 00302 << DelFibPartgrad << "\n d_qz = " << d_qz 00303 << "\n =====> Start Step Information <=== \n" 00304 << " ===> calo preStepPoint info <=== \n" 00305 << " hitPoint = " << hitPoint << "\n" 00306 << " hitMom = " << hit_mom << "\n" 00307 << " stepControlFlag = " << stepControlFlag 00308 // << "\n curprocess = " << curprocess << "\n" 00309 // << " nameProcess = " << nameProcess 00310 << "\n charge = " << charge << "\n" 00311 << " beta = " << beta << "\n" 00312 << " bThreshold = " << bThreshold << "\n" 00313 << " thgrad =" << thgrad << "\n" 00314 << " effPMTandTransport=" << effPMTandTransport 00315 // << "\n volume = " << name 00316 << "\n nameVolume = " << nameVolume << "\n" 00317 << " nMedium = " << nMedium << "\n" 00318 // << " rad length = " << rad << "\n" 00319 // << " material = " << mat << "\n" 00320 << " stepl = " << stepl << "\n" 00321 << " photEnSpectrDE = " << photEnSpectrDE <<"\n" 00322 << " edep = " << edep << "\n" 00323 << " ===> calo theTrack info <=== " << "\n" 00324 << " particleType = " << particleType << "\n" 00325 << " primaryID = " << primaryID << "\n" 00326 << " entot= " << entot << "\n" 00327 << " vert_eta= " << eta << "\n" 00328 << " vert_phi= " << phi << "\n" 00329 << " vert_mom= " << vert_mom << "\n" 00330 << " ===> calo hit preStepPointinfo <=== "<<"\n" 00331 << " local point = " << localPoint << "\n" 00332 << " ==============================> final info" 00333 << ": <=== \n" 00334 << " meanNCherPhot = " << meanNCherPhot << "\n" 00335 << " poissNCherPhot = " << poissNCherPhot <<"\n" 00336 << " NCherPhot = " << NCherPhot; 00337 #endif 00338 00339 } 00340 } 00341 00342 00343 #ifdef debug 00344 LogDebug("ForwardSim") << "CastorSD:: " << nameVolume 00345 // << " Light Collection Efficiency " << weight 00346 << " Weighted Energy Deposit " << edep/MeV 00347 << " MeV\n"; 00348 #endif 00349 return NCherPhot; 00350 } 00351 }
| uint32_t CastorSD::setDetUnitId | ( | G4Step * | step | ) | [virtual] |
Implements CaloSD.
Definition at line 353 of file CastorSD.cc.
References CastorNumberingScheme::getUnitID(), and numberingScheme.
00353 { 00354 return (numberingScheme == 0 ? 0 : numberingScheme->getUnitID(aStep)); 00355 }
| void CastorSD::setNumberingScheme | ( | CastorNumberingScheme * | scheme | ) |
Definition at line 357 of file CastorSD.cc.
References numberingScheme.
Referenced by CastorSD().
00357 { 00358 00359 if (scheme != 0) { 00360 edm::LogWarning("ForwardSim") << "CastorSD: updates numbering scheme for " 00361 << GetName(); 00362 if (numberingScheme) delete numberingScheme; 00363 numberingScheme = scheme; 00364 } 00365 }
CastorNumberingScheme* CastorSD::numberingScheme [private] |
1.5.4