diff --git a/diff.sh b/diff.sh index d9b65dc..d448540 100755 --- a/diff.sh +++ b/diff.sh @@ -11,7 +11,19 @@ # note: the script is meant to run one directory up from clas12tags prompt="no" -next_release=$(ls release_notes | grep '.md' | grep -v all | sort -u | tail -1 | awk -F. '{print $1"."$2}') + +next_release=$(ls release_notes | grep -v dev | grep '.md' | awk -F. '{print $1"."$2}' | sort -V | tail -1) + +# if dev.md is different than $next_release.md, then copy dev.md to $next_release.md +if [[ -f release_notes/dev.md ]]; then + if [[ -f release_notes/$next_release.md ]]; then + if [[ $(diff release_notes/dev.md release_notes/$next_release.md) ]]; then + cp release_notes/dev.md release_notes/$next_release.md + fi + else + cp release_notes/dev.md release_notes/$next_release.md + fi +fi # if argument is given, set prompt to yes if [[ $# -gt 0 ]]; then @@ -57,29 +69,3 @@ sed -i 's/const char.*/'$new_string'/' source/gemc.cc printf "\n- Changing initializeBMTConstants and initializeFMTConstants to initialize before processID" sed -i s/'initializeBMTConstants(-1)'/'initializeBMTConstants(1)'/ source/hitprocess/clas12/micromegas/BMT_hitprocess.cc sed -i s/'initializeFMTConstants(-1)'/'initializeFMTConstants(1)'/ source/hitprocess/clas12/micromegas/FMT_hitprocess.cc - -printf "\n- Removing evio support for clas12tags\n\n" -sed -i s/'env = init_environment("qt5 geant4 clhep evio xercesc ccdb mlibrary cadmesh hipo c12bfields")'/'env = init_environment("qt5 geant4 clhep xercesc ccdb mlibrary cadmesh hipo c12bfields")'/ source/SConstruct -sed -i s/'output\/evio_output.cc'/''/ source/SConstruct - - -sed -i s/'\/\/ EVIO'/''/ source/output/outputFactory.h -sed -i s/'#pragma GCC diagnostic push'/''/ source/output/outputFactory.h -sed -i s/'#pragma GCC diagnostic ignored "-Wdeprecated-declarations" '/''/ source/output/outputFactory.h -sed -i s/'#pragma GCC diagnostic ignored "-Wdeprecated"'/''/ source/output/outputFactory.h -sed -i s/'#include "evioUtil.hxx"'/''/ source/output/outputFactory.h -sed -i s/'#include "evioFileChannel.hxx"'/''/ source/output/outputFactory.h -sed -i s/'#pragma GCC diagnostic pop'/''/ source/output/outputFactory.h -sed -i s/'using namespace evio;'/''/ source/output/outputFactory.h -sed -i s/'evioFileChannel \*pchan;'/''/ source/output/outputFactory.h - - -sed -i s/'#include "evio_output.h"'/''/ source/output/outputFactory.cc -sed -i s/'\/\/ EVIO Buffer size set to 30M words'/''/ source/output/outputFactory.cc -sed -i s/'int evio_buffer = EVIO_BUFFER;'/''/ source/output/outputFactory.cc -sed -i s/'if(outType == "evio") {'/'{'/ source/output/outputFactory.cc -sed -i s/'pchan = new evioFileChannel(trimSpacesFromString(outFile).c_str(), "w", evio_buffer);'/''/ source/output/outputFactory.cc -sed -i s/'pchan->open();'/''/ source/output/outputFactory.cc -sed -i s/'outputMap\["evio"\] = &evio_output::createOutput;'/''/ source/output/outputFactory.cc -sed -i s/'pchan->close();'/''/ source/output/outputFactory.cc -sed -i s/'delete pchan;'/''/ source/output/outputFactory.cc diff --git a/release_notes/5.11.md b/release_notes/5.11.md new file mode 100644 index 0000000..3bf7477 --- /dev/null +++ b/release_notes/5.11.md @@ -0,0 +1,35 @@ +# clas12_tags dev + +- RG-D RG-D Flag Assembly geometry and variations (Lamiaa) +- added MIE scattering to api and source code (Connor) +- RICH hit process: PMT quantum efficiencies (Connor) +- Alert hit process improvements (M. Paolone and ) +- RTPC hit process improvements (YuchunHung) + +``` + > gemc version: gemc dev + + > Environment: + + > FIELD_DIR: /opt/jlab_software/noarch/data/magfield + > GEMC_DATA_DIR: /opt/jlab_software/1.1/macosx14-clang15/clas12Tags/dev + > G4DATA_DIR: /opt/jlab_software/1.1/macosx14-clang15/geant4/10.7.4/data/Geant4-10.7.4/data + > G4_VERSION: 10.7.4 + > G4INSTALL: /opt/jlab_software/1.1/macosx14-clang15/geant4/10.7.4 + +``` + +
+ +### To load production tag 5.10 at JLab, use the following commands: + +``` +source /group/clas12/packages/setup.[c]sh +module load clas12 +module switch gemc/5.10 +``` + + + +
+ \ No newline at end of file diff --git a/release_notes/dev.md b/release_notes/dev.md index e055d9a..06bd97e 100644 --- a/release_notes/dev.md +++ b/release_notes/dev.md @@ -1,6 +1,10 @@ # clas12_tags dev - RG-D RG-D Flag Assembly geometry and variations (Lamiaa) +- added MIE scattering to api and source code (Connor) +- RICH hit process: PMT quantum efficiencies (Connor) +- Alert hit process improvements (M. Paolone and fizikci0147) +- RTPC hit process improvements (YuchunHung) ``` > gemc version: gemc dev diff --git a/source/gemc.cc b/source/gemc.cc index 75998a9..51bce1c 100644 --- a/source/gemc.cc +++ b/source/gemc.cc @@ -28,7 +28,7 @@ /// \author \n © Maurizio Ungaro /// \author e-mail: ungaro@jlab.org\n\n\n -const char *GEMC_VERSION = "gemc dev.md" ; +const char *GEMC_VERSION = "gemc 5.11" ; // G4 headers #include "G4RunManager.hh" @@ -415,4 +415,3 @@ int main( int argc, char **argv ) // introducing OPTICALPHOTONPID here to be semi-transparent to G4 changes // this pid changed from 0 to -22 with geant4 10.7 int MHit::OPTICALPHOTONPID = -22; -// int MHit::OPTICALPHOTONPID = 0; diff --git a/source/hitprocess/clas12/alert/ahdc_hitprocess.cc b/source/hitprocess/clas12/alert/ahdc_hitprocess.cc index 060f19c..4dc5dc0 100644 --- a/source/hitprocess/clas12/alert/ahdc_hitprocess.cc +++ b/source/hitprocess/clas12/alert/ahdc_hitprocess.cc @@ -256,6 +256,19 @@ map ahdc_HitProcess::integrateDgt(MHit* aHit, int hitn) { // time calculation // signal_t = stepTime[s] + (H_abh/driftVelocity); // cout << "signal_t: " << signal_t << ", stepTime: " << stepTime[s] << endl; + + // docasig is a fit to sigma vs distance plot. A second order pol used for the fit (p0+p1*x+p2*x*x). + // drift velocity as a function of distance. pol2 fitted to t vs x plot and drift velocity is derived from the fit (1/(dt/dx)). + // both sig vs dist and t vs dist plots are taken from Lucien Causse's PhD thesis ("Development of a stereo drift chamber for the Jefferson Laboratory ALERT Experiment."). + // plots were digitized and then fitted to a pol2. + + double driftP1=-16.17; + double driftP2=24.81; + double docasig = 337.3-210.3*doca+34.7*pow(doca,2); + std::default_random_engine dseed(time(0)); //seed + std::normal_distribution ddist(doca, docasig); //a resuolution affect is added to doca. + double doca_r =ddist(dseed); + driftVelocity = 1./(driftP1+2.*driftP2*doca_r); // mm/ns // drift velocity as a function of distance. pol2 fitted to t vs x plot and drift velocity is then extracted from dx/dt, d/dt(p0+p1*x+p2*x^2)=p1+2*p2*x. signal_tTimesEdep = signal_tTimesEdep + (stepTime[s] + H_abh/driftVelocity) * E_wire; // cout << "signal_tTimesEdep: " << signal_tTimesEdep << endl; @@ -271,7 +284,8 @@ map ahdc_HitProcess::integrateDgt(MHit* aHit, int hitn) { // Just to test, time is linear with doca double a = 5.0; double b = 5.0; - double time = a*doca+b + signal_tTimesEdep/E_tot_wire; + //double time = a*doca+b + signal_tTimesEdep/E_tot_wire; + double time = signal_tTimesEdep/E_tot_wire; // double signal = 0.0; // signal = signal_tTimesEdep/E_tot_wire; diff --git a/source/hitprocess/clas12/rich_hitprocess.cc b/source/hitprocess/clas12/rich_hitprocess.cc index b14d936..98821c1 100644 --- a/source/hitprocess/clas12/rich_hitprocess.cc +++ b/source/hitprocess/clas12/rich_hitprocess.cc @@ -177,15 +177,43 @@ map rich_HitProcess :: integrateDgt(MHit* aHit, int hitn) writeHit = true; rejectHitConditions = false; - + + int pmtType = 12700; + // sector 4: mix of H12700 and H8500 + if(idsector == 4){ + pmtType = richc.pmtType[idpmt-1]; + } double energy = aHit->GetEs()[0]/electronvolt; double qeff = 0; - for(int i = 0; i < richc.nQEbins; i++){ - if(energy < richc.Ene[i] && energy > richc.Ene[i+1]){ - qeff = richc.QE[i]; - break; + + if(pmtType == 8500){ + for(int i = 0; i < richc.nQEbinsH8500; i++){ + if(energy < richc.Ene_H8500[i] && energy > richc.Ene_H8500[i+1]){ + if( std::abs(energy - richc.Ene_H8500[i]) < std::abs(energy - richc.Ene_H8500[i+1])){ + qeff = richc.QE_H8500[i]; + } + else{ + qeff = richc.QE_H8500[i+1]; + } + break; + } } } + else if(pmtType == 12700){ + for(int i = 0; i < richc.nQEbinsH12700; i++){ + if(energy < richc.Ene_H12700[i] && energy > richc.Ene_H12700[i+1]){ + if( std::abs(energy - richc.Ene_H12700[i]) < std::abs(energy - richc.Ene_H12700[i+1])){ + qeff = richc.QE_H12700[i]; + } + else{ + qeff = richc.QE_H12700[i+1]; + } + break; + } + } + } + + // applying quantum efficiency from thrown random value set in integrateDgt if( identity[2].userInfos[3] > qeff && !aHit->isElectronicNoise) { @@ -223,8 +251,14 @@ vector rich_HitProcess :: processID(vector id, G4Step* a G4ThreeVector pixelCenterLocal = getPixelCenter(pixel); G4ThreeVector pixelCenterGlobal = prestep->GetTouchableHandle()->GetHistory()->GetTopTransform().Inverse().TransformPoint(pixelCenterLocal); + int idsector = yid[0].id; int pmt = yid[1].id; - RichPixel richPixel(richc.pmtType[pmt-1]); + int pmtType = 12700; + // sector 4: mix of H8500 and H12700 + if(idsector==4){ + pmtType = richc.pmtType[pmt-1]; + } + RichPixel richPixel(pmtType); richPixel.Clear(); int t1 = -1; diff --git a/source/hitprocess/clas12/rich_hitprocess.h b/source/hitprocess/clas12/rich_hitprocess.h index 2ef2940..cd7ee65 100644 --- a/source/hitprocess/clas12/rich_hitprocess.h +++ b/source/hitprocess/clas12/rich_hitprocess.h @@ -254,9 +254,9 @@ class richConstants 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 3, 1, 2, 1, 2, 3, 1, 2, 3}; - // quantum efficiency table - const static int nQEbins=87; - double Ene[nQEbins]={6.148, 5.387, + // quantum efficiency tables + const static int nQEbinsH8500=87; + double Ene_H8500[nQEbinsH8500]={6.148, 5.387, 4.624, 4.558, 4.510, 4.432, 4.372, 4.285, 4.201, 4.081, 3.980, 3.815, 3.726, 3.591, 3.503, 3.410, 3.314, 3.214, 3.128, 3.054, 2.983, 2.902, 2.845, 2.783, 2.724, 2.668, 2.624, 2.572, 2.531, 2.487, 2.454, 2.431, @@ -267,7 +267,7 @@ class richConstants 1.919, 1.914, 1.908, 1.902, 1.897, 1.894, 1.888, 1.883, 1.877, 1.872, 1.867, 1.861, 1.858, 1.853, 1.848}; - double QE[nQEbins]={0.000, 0.000, + double QE_H8500[nQEbinsH8500]={0.000, 0.000, 0.11241, 0.12738, 0.14272, 0.16358, 0.18120, 0.20301, 0.22234, 0.24351, 0.25775, 0.27282, 0.27594, 0.27594, 0.27594, 0.27282, 0.26974, 0.26669, 0.26070, 0.25484, 0.24629, 0.23534, 0.22745, 0.21489, 0.20301, 0.18963, 0.17713, 0.16358, 0.15279, 0.13951, 0.12738, 0.11500, @@ -277,7 +277,25 @@ class richConstants 0.00317, 0.00276, 0.00247, 0.00215, 0.00192, 0.00169, 0.00150, 0.00130, 0.00118, 0.00103, 0.00093, 0.00082, 0.00072, 0.00063, 0.00056, 0.00050, 0.00043, 0.00039, 0.00035, 0.00030, 0.00027, 0.00023, 0.00021, 0.00018, 0.00016}; - + + // QE from official Hamamatsu table + const static int nQEbinsH12700 = 46; + double Ene_H12700[nQEbinsH12700] = {4.59185185, 4.42785714, 4.27517241, 4.13266667, 3.99935484, + 3.874375 , 3.7569697 , 3.64647059, 3.54228571, 3.44388889, + 3.35081081, 3.26263158, 3.17897436, 3.0995 , 3.02390244, + 2.95190476, 2.88325581, 2.81772727, 2.75511111, 2.69521739, + 2.63787234, 2.58291667, 2.53020408, 2.4796 , 2.43098039, + 2.38423077, 2.33924528, 2.29592593, 2.25418182, 2.21392857, + 2.17508772, 2.13758621, 2.10135593, 2.06633333, 2.03245902, + 1.99967742, 1.96793651, 1.9371875 , 1.90738462, 1.87848485, + 1.85044776, 1.82323529, 1.79681159, 1.77114286, 1.74619718, + 1.72194444}; + double QE_H12700[nQEbinsH12700] = { 0.1351, 0.2071, 0.2698, 0.3092, 0.3333, 0.3458, 0.3512, 0.3512, 0.3491, 0.3431, 0.3420, 0.3434, + 0.3348, 0.3284, 0.3198, 0.3112, 0.30, 0.2858, 0.2666, 0.2484, 0.2331, 0.2210, 0.2078, + 0.1817, 0.1440, 0.1169, 0.1010, 0.0895, 0.0793, 0.0698, 0.0605, 0.0515, 0.0428, + 0.0347, 0.0271, 0.0204, 0.0147, 0.01, 0.0065, 0.004, 0.0023, 0.0013, 0.0007, 0.0003, 0.0002, 0.0001}; + + // two types of pmts used in sector 4 rich int pmtType[391] = {12700, 12700, 12700, 12700, 12700, 12700, 8500, 8500, 12700, 12700, 12700, 8500, 8500, 8500, 8500, 8500, 8500, 8500, 8500, 8500, 8500, 8500, 8500, 8500, 8500, 8500, 8500, 12700, 12700, 12700, 8500, 8500, 12700, 12700, 12700, 12700, @@ -347,7 +365,7 @@ class rich_HitProcess : public HitProcess G4ThreeVector getPixelCenter(int pixel); // testing ccdb time paramters vs PMT simulation class - bool ccdbTiming = true; + bool ccdbTiming = false; // just converting double tdc to int for 1ns tdc precision double tdc_precision = 1.; diff --git a/source/hitprocess/clas12/rtpc_hitprocess.cc b/source/hitprocess/clas12/rtpc_hitprocess.cc index 9cb88b0..d655554 100644 --- a/source/hitprocess/clas12/rtpc_hitprocess.cc +++ b/source/hitprocess/clas12/rtpc_hitprocess.cc @@ -1,4 +1,3 @@ -#include "rtpc_hitprocess.h" #include #include @@ -9,66 +8,152 @@ #include "CLHEP/Units/PhysicalConstants.h" using namespace CLHEP; +// ccdb +#include +#include +#include +using namespace ccdb; + +// gemc headers +#include "rtpc_hitprocess.h" + +// 2023.10.03 copied from rich_hitprocess.cc. Added in case for the future needed. (y.-c.) +static rtpcConstants initializeRTPCConstants(int runno, string digiVariation = "default", string digiSnapshotTime = "no", bool accountForHardwareStatus = false) +{ + // all these constants should be read from CCDB + rtpcConstants rtpcc; + + // do not initialize at the beginning, only after the end of the first event, + // with the proper run number coming from options or run table + if(runno == -1) return rtpcc; + string timestamp = ""; + if(digiSnapshotTime != "no") { + timestamp = ":"+digiSnapshotTime; + } + + cout << "Entering initializeRTPCConstants" << endl; + + // database + rtpcc.runNo = runno; + + if(getenv ("CCDB_CONNECTION") != nullptr) + rtpcc.connection = (string) getenv("CCDB_CONNECTION"); + else + rtpcc.connection = "mysql://clas12reader@clasdb.jlab.org/clas12"; + + unique_ptr calib(CalibrationGenerator::CreateCalibration(rtpcc.connection)); + cout << "Connecting to " << rtpcc.connection << "/calibration/rtpc" << endl; + + + // get time parameter from CCDB (instead of hard code in the hitprocess file.) + vector > data; + + cout << "RTPC:Getting drift parameters" << endl; + snprintf(rtpcc.database, sizeof(rtpcc.database), "/calibration/rtpc/time_parms:%d:%s%s", rtpcc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); + calib -> GetCalib(data, rtpcc.database); + int iComponent = 0; + for (unsigned row = 0; row < data.size(); row++) { + //iSector = data[row][0]; + //iLayer = data [row][1]; + iComponent = data[row][2]; + rtpcc.z0[iComponent - 1] = data[row][3]; + rtpcc.z2[iComponent - 1] = data[row][5]; + rtpcc.z4[iComponent - 1] = data[row][7]; + + //cout << row << ": " << data[row][0] << " " << data[row][1] << " " << data[row][2] << " " << data[row][3] << " " << data[row][4] << " " << data[row][5] << " " << data[row][6] << " " << data[row][7] << endl; + } // row++ + // the last row is used for tdiff (not used in old hitprocess y.-c. Oct30.2023) + +/* + cout << "RTPC:Getting gain balance" << endl; // keep for future or reference + snprintf(rtpcc.database, sizeof(rtpcc.database), "/calibration/rtpc/gain_balance:%d:%s%s", rtpcc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); + calib -> GetCalib(data, rtpcc.database); + iComponent = 0; + int iLayer = 0; + for (unsigned row = 0; row < data.size(); row++) { // data.size() = 180*96 = 17280 + //iSector = data[row][0]; + iLayer = data[row][1]; + iComponent = data[row][2]; + rtpcc.gain_balance[iComponent][iLayer] = data[row][3]; // loop layer first, and then component in gain_balance table + } // row++ +*/ + + //rtpcc.date = "2016-03-15"; // useful? + //rtpcc.variation = "default"; + + // These are not in the CCDB (needed? y.-c. Oct30.2023) + + return rtpcc; +} + + map rtpc_HitProcess :: integrateDgt(MHit* aHit, int hitn) { - static const double PI=3.1415926535; + rejectHitConditions = false; writeHit = true; - // Establish constants - // NEW CONSTANTS FOR UPDATED RECONSTRUCTION - PAD_W =2.79; - PAD_L = 4.0; - PAD_S = 79.0;//=80.0 old version after comment - RTPC_L = 384.0; - phi_per_pad=(2.0*PI)/180; - // Drift time from first GEM to readout pad - t_2GEM2 = 0.0; // 169.183; - sigma_t_2GEM2 = 8.72728; - t_2GEM3 = 0.0; //222.415; - sigma_t_2GEM3 = 5.62223; - t_2PAD = 0.0; //414.459; - sigma_t_2PAD = 7.58056; - - phi_2GEM2 = 0.0; // 0.0416925; - sigma_phi_2GEM2 = 0.00384579; - phi_2GEM3 = 0.0; // 0.0416574; - sigma_phi_2GEM3 = 0.00160235; - phi_2PAD = 0.0; // 0.057566; - sigma_phi_2PAD = 0.00238653; + float t_2GEM2 = 0.0; // 169.183; + float t_2GEM3 = 0.0; //222.415; + float t_2PAD = 0.0; //414.459; // parameters to change for gas mixture and potential // gas mixture = He:CO2 at 85.5:14.5 // potential = 4450V // --------- Updated 25 June 2020 B-field parameters ----------- // - double a_t0 = 1470.0; // -0.607056; - double a_t1 = 0.009; // -0.00138137; + //double a_t0 = 1470.0; // -0.607056; + //double a_t1 = 0.009; // -0.00138137; //double a_t2 = -9.0e05; // 1.6494e-05; //Gave crazy big -ve values to tdc or time - double a_t2 = -9.0e-05; // 1.6494e-05; - - double b_t0 = 3290.0; // 3183.62; - double b_t1 = -0.0138; // -0.0269385; - double b_t2 = -0.000332; // -0.00037645; + //double a_t2 = -9.0e-05; // 1.6494e-05; - double c_t0 = 0.0; // -0.0200932; - double c_t1 = 0.0; // 0.00319441; - double c_t2 = 0.0; // -8.75299e-06; + //double b_t0 = 3290.0; // 3183.62; + //double b_t1 = -0.0138; // -0.0269385; + //double b_t2 = -0.000332; // -0.00037645; - double a_phi0 = -0.113; // 0; - double a_phi1 = -1.05e-6; // 0; - double a_phi2 = 1.58e-8; // 0; + //double c_t0 = 0.0; // -0.0200932; + //double c_t1 = 0.0; // 0.00319441; + //double c_t2 = 0.0; // -8.75299e-06; - double b_phi0 = -0.9697; // 1.22249; - double b_phi1 = 9.48e-6; // -2.3708e-05; - double b_phi2 = 9.9e-8; // -1.11055e-07; + //double a_phi0 = -0.113; // 0; + //double a_phi1 = -1.05e-6; // 0; + //double a_phi2 = 1.58e-8; // 0; - double c_phi0 = 0.0; // -1.13197; - double c_phi1 = 0.0; // 0.000118575; - double c_phi2 = 0.0; // 2.50094e-08; + //double b_phi0 = -0.9697; // 1.22249; + //double b_phi1 = 9.48e-6; // -2.3708e-05; + //double b_phi2 = 9.9e-8; // -1.11055e-07; + //double c_phi0 = 0.0; // -1.13197; + //double c_phi1 = 0.0; // 0.000118575; + //double c_phi2 = 0.0; // 2.50094e-08; + + //cout << "In integrateDgt: a_t para = " << a_t0 << " " << a_t1 << " " << a_t2 << endl; + //for (int i = 0; i < 7; i++) { + // cout << "In integrateDgt: CCDB cons " << i <<" = " << rtpcc.z0[i] << " " << rtpcc.z2[i] << " " << rtpcc.z4[i] << endl; + //} + //cout << " ------------------------------------- " << endl; + + // Read drift parameters from CCDB --> checked! Same as the hard code valuse above. + double a_t0 = rtpcc.z0[0]; + double a_t1 = rtpcc.z2[0]; + double a_t2 = rtpcc.z4[0]; + + double b_t0 = rtpcc.z0[1]; + double b_t1 = rtpcc.z2[1]; + double b_t2 = rtpcc.z4[1]; + // phi offset (delta_phi) + double a_phi0 = rtpcc.z0[3]; + double a_phi1 = rtpcc.z2[3]; + double a_phi2 = rtpcc.z4[3]; + // tan_phi + double b_phi0 = rtpcc.z0[4]; + double b_phi1 = rtpcc.z2[4]; + double b_phi2 = rtpcc.z4[4]; + // Diffusion parameters diff_at = 388.7449859; @@ -77,63 +162,17 @@ map rtpc_HitProcess :: integrateDgt(MHit* aHit, int hitn) diff_aphi = 6.00E-06; diff_bphi = 2.00E-06; - a_z=0.035972097; - b_z =-0.000739386; + a_z = 0.035972097; + b_z = -0.000739386; - t_2END = t_2GEM2 + t_2GEM3 + t_2PAD; + float t_2END = t_2GEM2 + t_2GEM3 + t_2PAD; sigma_t_gap = sqrt(pow(sigma_t_2GEM2,2)+pow(sigma_t_2GEM3,2)+pow(sigma_t_2PAD,2)); phi_2END = phi_2GEM2 + phi_2GEM3 + phi_2PAD; sigma_phi_gap = sqrt(pow(sigma_phi_2GEM2,2)+pow(sigma_phi_2GEM3,2)+pow(sigma_phi_2PAD,2)); - // ----------------------------------------- // - - // ------------ B=0 T Parameters ------------// - /*double a_t1 = 0; - double a_t2 = 0; - double a_t3 = 0; - double a_t4 = 0; - double a_t5 = 6.96387e+02; - - - double b_t1 = 0; - double b_t2 = 0; - double b_t3 = 0; - double b_t4 = 0; - double b_t5 = -4.73759e+01; - - double a_phi1 = 0; - double a_phi2 = 0; - double a_phi3 = 0; - double a_phi4 = 0; - double a_phi5 = 0; - - double b_phi1 = 0; - double b_phi2 = 0; - double b_phi3 = 0; - double b_phi4 = 0; - double b_phi5 = 0; - - // Diffusion parameters - c_t=0.0; - d_t=0.0; - - c_phi=0.0; - d_phi=0.0; - - a_z=0.0; - b_z =0.0; - - t_2END = 500.0; - sigma_t_gap = 0; - phi_2END = 0.0; - sigma_phi_gap = 0.0;*/ - // ------------------------------------------// TPC_TZERO = 0.0; - // not used anymore? - //int chan=0; - map dgtz; vector identity = aHit->GetId(); @@ -151,7 +190,6 @@ map rtpc_HitProcess :: integrateDgt(MHit* aHit, int hitn) // so tInfos.eTot is the sum of all steps s of Edep[s] vector Edep = aHit->GetEdep(); - // -------------------------- TIME SHIFT for non-primary tracks --------------------------- //M. U.: map.find(Key) tries to find an entry with that key, if it doesn't find it, it will return map.end() // So, to ensure that hits of a given track gets the same timeShift (but random among different tracks) @@ -178,134 +216,122 @@ map rtpc_HitProcess :: integrateDgt(MHit* aHit, int hitn) // // Get the information x,y,z and Edep at each ionization point // - double LposX=0.; - double LposY=0.; - double LposZ=0.; - double DiffEdep=0.; - if(tInfos.eTot > 0) - - { - - int adc =0; - double tdc =0; + if(tInfos.eTot > 0) { + + // Idea: After collecting the geant4 steps as a hit in processID, re-evaluate the destination on the readout padboard of this hit. Therefore, the final col and row write out in RTPC::adc bank might not be the same as the original(in processID). Futhermore, there might be some different hits have same col and row in RTPC::adc bank. + // y.-c. comment in Dec. 23. 2023 + + // total energy deposit for a hit + double DiffEdep = tInfos.eTot; + int adc = ((int)100000.0*DiffEdep); //cludge for tiny ADC numbers + + //convert (x0,y0,z0) into (r0,phi0,z0) + z_cm = tInfos.lz/10.0; // mm -> cm + + double r0 = (sqrt(tInfos.lx*tInfos.lx + tInfos.ly*tInfos.ly))/10.0; //in cm + double phi0_rad = atan2(tInfos.ly, tInfos.lx); //return (-Pi, + Pi) + if (phi0_rad < 0.) phi0_rad += 2.0*PI; // convert to (0, 2PI) + //if (phi0_rad >= 2.0*PI) phi0_rad -= 2.0*PI; useless because -pi <= phi0_rad <= pi + + // calculate drift time of a hit + // all in cm + a_t = a_t0 + a_t1*(pow(z_cm,2)) + a_t2*(pow(z_cm,4)); // t_offset + b_t = b_t0 + b_t1*(pow(z_cm,2)) + b_t2*(pow(z_cm,4)); // T_max + + // --------------------- Addition of Diffusion ----------------- // + // calculate drift time [ns] to first GEM + double t_drift = a_t + b_t*(7.0*7.0 - r0*r0)/40.0; + // determine sigma of drift time [ns] + double t_diff = sqrt(diff_at*(7.0 - r0) + diff_bt*(7.0 - r0)*(7.0 - r0)); + + // find t_s2pad by gaussians + double t_s2pad = G4RandGauss::shoot(t_drift, t_diff); + double t_gap = G4RandGauss::shoot(t_2END, sigma_t_gap); + // time shift + shift_t = timeShift_map.find(otids[0])->second; //To ensure secondaries have the same shift_t as primaries or acencestors + //tdc = t_s2pad + t_gap + shift_t; //kp: + double tdc = t_s2pad + t_gap;//+shift_t; scattered electron and spectator proton do not have shift_t + + + + // calculate drift angle to first GEM at 7 cm [rad] + // all in cm + a_phi = a_phi0 + a_phi1*(pow(z_cm,2)) + a_phi2*(pow(z_cm,4)); //delta_phi + b_phi = b_phi0 + b_phi1*(pow(z_cm,2)) + b_phi2*(pow(z_cm,4)); //tan_theta + // phi_drift = delta_phi + tanTheta*ln(r_max/r); + double phi_drift = a_phi + b_phi*log(7.0/r0); + + // determine sigma of drift angle [rad] + double phi_diff = sqrt(diff_aphi*(7.0 - r0) + diff_bphi*(7.0 - r0)*(7.0 - r0)); + + // find delta_phi by gaussians + double delta_phi = G4RandGauss::shoot(phi_drift, phi_diff); + double phi_gap = G4RandGauss::shoot(phi_2END, sigma_phi_gap); + + double phi_rad= phi0_rad + delta_phi + phi_gap; //phi at pad pcb board + if (phi_rad < 0.0) phi_rad += 2.0*PI; // return to (0, 2PI) + if (phi_rad >= 2.0*PI) phi_rad -= 2.0*PI; // return to (0, 2PI) + + + // calculate drift in z [cm] + double z_drift = 0.0; + + // determine sigma in z [cm] + double z_diff = sqrt(a_z*(7.0 - r0)+b_z*(7.0 - r0)*(7.0 - r0)); + + // find delta_z by gaussians + double delta_z = G4RandGauss::shoot(z_drift, z_diff); + + //double z_pos = (z_cm + delta_z)*10.0; // mm + double z_pos = z_cm*10.0 + delta_z; // mm + + + + // convert physical position into row & col of each pad. + int col = -999; + int row = -999; + + row = ceil(phi_rad/phi_per_pad); + float z_shift = (row - 1)%4; + + float col_min = -RTPC_L/2.0 + z_shift; + float col_max = RTPC_L/2.0 + z_shift; + + if (z_pos < col_min || z_pos > col_max) {row = -999; col = -999;} + else { + row = ceil(phi_rad/phi_per_pad); + col = ceil((z_pos + RTPC_L/2.0 - z_shift)/PAD_L); + } // consistent with PadVector.java in coatjava + + +// row = identity[0].id; +// col = identity[1].id; - for(unsigned int s=0; s70.) continue; - // if(LposZ<-(RTPC_L/2.0) || LposZ>(RTPC_L/2.0)) continue; - - DiffEdep = Edep[s]; - - z_cm = LposZ/10.0; - - // all in cm - a_t= a_t0 + a_t1*(pow(z_cm,2)) + a_t2*(pow(z_cm,4)); - b_t= b_t0 + b_t1*(pow(z_cm,2)) + b_t2*(pow(z_cm,4)); - c_t= c_t0 + c_t1*(pow(z_cm,2)) + c_t2*(pow(z_cm,4)); - - a_phi=a_phi2*(pow(z_cm,4)) + a_phi1*(pow(z_cm,2)) + a_phi0; - b_phi=b_phi2*(pow(z_cm,4)) + b_phi1*(pow(z_cm,2)) + b_phi0; - c_phi=c_phi2*(pow(z_cm,4)) + c_phi1*(pow(z_cm,2)) + c_phi0; - - double r0,phi0_rad; - //convert (x0,y0,z0) into (r0,phi0,z0) - r0=(sqrt(LposX*LposX+LposY*LposY))/10.0; //in cm - - phi0_rad=atan2(LposY,LposX); //return (-Pi, + Pi) - if( phi0_rad<0.) phi0_rad+=2.0*PI; - if( phi0_rad>=2.0*PI ) phi0_rad-=2.0*PI; - - - // -------------------------------- Addition of Diffusion ----------------------------- - - // calculate drift time [ns] to first GEM - double t_drift = a_t + b_t*(7.0*7.0-r0*r0)/40.0; - - // determine sigma of drift time [ns] - double t_diff = sqrt(diff_at*(7.0-r0)+diff_bt*(7.0-r0)*(7.0-r0)); - - // calculate drift angle to first GEM at 7 cm [rad] - double phi_drift = a_phi + b_phi*log(7.0/r0) + c_phi*((1.0/(r0*r0)) - (1.0/(49.0))); - //double phi_drift = a_phi*(7.0-r0)+b_phi*(7.0-r0)*(7.0-r0); - - // determine sigma of drift angle [rad] - double phi_diff = sqrt(diff_aphi*(7.0-r0)+diff_bphi*(7.0-r0)*(7.0-r0)); - - // calculate drift in z [mm] - double z_drift = 0.0; - - // determine sigma in z [mm] - double z_diff = sqrt(a_z*(7.0-r0)+b_z*(7.0-r0)*(7.0-r0)); - - // find t_s2pad and delta_phi by gaussians - double t_s2pad = G4RandGauss::shoot(t_drift, t_diff); - double delta_phi = G4RandGauss::shoot(phi_drift, phi_diff); - double delta_z = G4RandGauss::shoot(z_drift, z_diff); - double t_gap = G4RandGauss::shoot(t_2END, sigma_t_gap); - double phi_gap = G4RandGauss::shoot(phi_2END, sigma_phi_gap); - - // ------------------------------------------------------------------------------------ - - double phi_rad= phi0_rad+delta_phi+phi_gap; //phi at pad pcb board - if( phi_rad<0.0 ) phi_rad+=2.0*PI; - if( phi_rad>=2.0*PI ) phi_rad-=2.0*PI; - - // time shift - //shift_t = timeShift_map.find(aHit->GetTId())->second; - shift_t = timeShift_map.find(otids[s])->second; //To ensure secondaries have the same shift_t as primaries or acencestors - - // NO time shift - //shift_t = 0.0; - - tdc=t_s2pad+t_gap+shift_t; //kp: - adc=((int)100000.0*DiffEdep); //cludge for tiny ADC numbers - //adc = 0.5 + adc; //kp: For debugging purpose to avoid lots of zeros (due to int typecasting below) - //adc = 100000*DiffEdep;//kp: For same reason as bove - - //double z_pos = (z_cm+delta_z)*10.0;// must be mm - double z_pos = z_cm*10.0 + delta_z;// must be mm - int col = -999; - int row = -999; - - row = ceil(phi_rad/phi_per_pad); - float z_shift = row%4; - - float col_min = -RTPC_L/2.0+z_shift; - float col_max = RTPC_L/2.0+z_shift; - - if( z_pos < col_min || z_pos > col_max ) {row = -999; col = -999;} - else { - - row = ceil(phi_rad/phi_per_pad); - col = ceil((z_pos+RTPC_L/2.0-z_shift)/PAD_L); - } - - //int Num_of_Col = (int) (RTPC_L/PAD_L); - //chan = row*Num_of_Col+col; - - //Mohammad: kill all the hits in the last three rows that correspond to the seam effect - if(row>177){row = -999; col = -999;} - - dgtz["Sector"] = 1; - dgtz["Layer"] = col; - dgtz["Component"] = row; - dgtz["Order"] = 0; - dgtz["Time"] = tdc; - dgtz["ADC"] = (int) adc; - dgtz["Ped"] = 0; - //dgtz["TimeShift"] = shift_t; - dgtz["hitn"] = (int) hitn; - - } // end step + //Mohammad: kill all the hits in the last three rows that correspond to the seam effect + if(row > 177){row = -999; col = -999;} + + + // write out to adc bank + dgtz["sector"] = 1; + dgtz["layer"] = col; + dgtz["component"] = row; + dgtz["ADC_order"] = 0; + dgtz["ADC_time"] = tdc; + dgtz["ADC_ADC"] = (int) adc; + dgtz["ADC_ped"] = 0; - } - + //dgtz["ADC_hitn"] = (int) hitn; // not define in CLARA see in data.json + + //reject hit write out to the adc bank, 30 was a number assigned for test. + // apply with the efficiency table, for instance, fot protential future impovement + //if (adc < 30) + // rejectHitConditions = true; + + } // tInfos.eTot > 0 + + // define conditions to reject hit if(rejectHitConditions) { writeHit = false; @@ -318,37 +344,26 @@ map rtpc_HitProcess :: integrateDgt(MHit* aHit, int hitn) vector rtpc_HitProcess :: processID(vector id, G4Step* aStep, detector Detector) { - static const double PI=3.1415926535; - - // Establish constants - PAD_W = 2.758; //2.79 !! - PAD_L = 4.0; - PAD_S = 79.0; - RTPC_L = 384.0; - phi_per_pad=PAD_W/PAD_S; - - // Drift time from first GEM to readout pad - phi_2GEM2 = 0.0; // 0.0416925; - sigma_phi_2GEM2 = 0.00384579; - phi_2GEM3 = 0.0; // 0.0416574; - sigma_phi_2GEM3 = 0.00160235; - phi_2PAD = 0.0; // 0.057566; - sigma_phi_2PAD = 0.00238653; + //cout << " In processID ***************" << endl; // parameters to change for gas mixture and potential // gas mixture = He:CO2 at 85.5:14.5 // potential = 4450V // --------- Updated 14 Dec 2020 B-field parameters - double a_phi0 = -0.113; // 0; - double a_phi1 = -1.05e-6; // 0; - double a_phi2 = 1.58e-8; // 0; - double b_phi0 = -0.9697; // 1.22249; - double b_phi1 = 9.48e-6; // -2.3708e-05; - double b_phi2 = 9.9e-8; // -1.11055e-07; - double c_phi0 = 0.0; // -1.13197; - double c_phi1 = 0.0; // 0.000118575; - double c_phi2 = 0.0; // 2.50094e-08; + + // Read drift angle parameters from CCDB --> checked! Same as the hard code valuse above. + // phi offset (delta_phi) + double a_phi0 = rtpcc.z0[3]; + double a_phi1 = rtpcc.z2[3]; + double a_phi2 = rtpcc.z4[3]; + // tan_phi + double b_phi0 = rtpcc.z0[4]; + double b_phi1 = rtpcc.z2[4]; + double b_phi2 = rtpcc.z4[4]; + + //cout << "row 4: " << a_phi0 << " " << a_phi1 << " " << a_phi2 << endl; + //cout << "row 5: " << b_phi0 << " " << b_phi1 << " " << b_phi2 << endl; // Diffusion parameters diff_aphi = 6.00E-06; @@ -357,56 +372,54 @@ vector rtpc_HitProcess :: processID(vector id, G4Step* a_z = 0.035972097; b_z = -0.000739386; - + phi_2END = phi_2GEM2 + phi_2GEM3 + phi_2PAD; sigma_phi_gap = sqrt(pow(sigma_phi_2GEM2,2)+pow(sigma_phi_2GEM3,2)+pow(sigma_phi_2PAD,2)); vector yid = id; - G4ThreeVector xyz = aStep->GetPostStepPoint()->GetPosition(); - G4ThreeVector Lxyz = aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory() - ->GetTopTransform().TransformPoint(xyz);///< Local Coordinates of interaction + // Get Coordinates of interaction from Geant4 (aStep) + G4ThreeVector xyz = aStep -> GetPostStepPoint() -> GetPosition(); // global + G4ThreeVector Lxyz = aStep -> GetPreStepPoint() -> GetTouchableHandle() -> GetHistory() -> GetTopTransform().TransformPoint(xyz); // local - int chan = 0; - double LposX = 0; - double LposY = 0; - double LposZ = 0; - LposX = Lxyz.x(); - LposY = Lxyz.y(); - LposZ = Lxyz.z(); + //int chan = 0; + double LposX = Lxyz.x(); + double LposY = Lxyz.y(); + double LposZ = Lxyz.z(); - z_cm = LposZ/10.0; + z_cm = LposZ/10.0; // convert mm to cm - // all in cm - a_phi=a_phi2*(pow(z_cm,4)) + a_phi1*(pow(z_cm,2)) + a_phi0; - b_phi=b_phi2*(pow(z_cm,4)) + b_phi1*(pow(z_cm,2)) + b_phi0; - c_phi=c_phi2*(pow(z_cm,4)) + c_phi1*(pow(z_cm,2)) + c_phi0; + // all in cm, because the drift parameters were extracted in cm. + a_phi = a_phi0 + a_phi1*(pow(z_cm,2)) + a_phi2*(pow(z_cm,4)); // delta_phi + b_phi = b_phi0 + b_phi1*(pow(z_cm,2)) + b_phi2*(pow(z_cm,4)); // tan_theta + //c_phi=c_phi2*(pow(z_cm,4)) + c_phi1*(pow(z_cm,2)) + c_phi0; - double r0,phi0_rad; + //convert (x0,y0,z0) into (r0,phi0,z0) - r0=(sqrt(LposX*LposX+LposY*LposY))/10.0; //in cm - + double r0 = (sqrt(LposX*LposX + LposY*LposY))/10.0; //in cm - phi0_rad=atan2(LposY,LposX); //return (-Pi, + Pi) - if( phi0_rad<0.) phi0_rad+=2.0*PI; - if( phi0_rad>=2.0*PI) phi0_rad-=2.0*PI; + double phi0_rad = atan2(LposY,LposX); //return (-Pi, + Pi) + if (phi0_rad < 0.) phi0_rad += 2.0*PI; // convert to (0, 2PI) + //if (phi0_rad >= 2.0*PI) phi0_rad -= 2.0*PI; useless because -pi <= phi0_rad <= pi + // -------------------------------- Addition of Diffusion ----------------------------- // calculate drift angle to first GEM at 7 cm [rad] - double phi_drift = a_phi + b_phi*log(7.0/r0) + c_phi*((1.0/(r0*r0)) - (1.0/(49.0))); - //double phi_drift = a_phi*(7.0-r0)+b_phi*(7.0-r0)*(7.0-r0); + // phi_drift = delta_phi + tanTheta*ln(r_max/r); + double phi_drift = a_phi + b_phi*log(7.0/r0); + //double phi_drift = a_phi + b_phi*log(7.0/r0) + c_phi*((1.0/(r0*r0)) - (1.0/(49.0))); - // determine sigma of drift angle [rad] - double phi_diff = sqrt(diff_aphi*(7.0-r0)+diff_bphi*(7.0-r0)*(7.0-r0)); + // determine sigma of drift angle [rad] + double phi_diff = sqrt(diff_aphi*(7.0 - r0) + diff_bphi*(7.0 - r0)*(7.0 - r0)); - // calculate drift in z [mm] + // calculate drift in z [cm] double z_drift = 0.0; - // determine sigma in z [mm] - double z_diff = sqrt(a_z*(7.0-r0)+b_z*(7.0-r0)*(7.0-r0)); + // determine sigma in z [cm] + double z_diff = sqrt(a_z*(7.0 - r0)+b_z*(7.0 - r0)*(7.0 - r0)); - // find t_s2pad and delta_phi by gaussians + // find delta_phi by gaussians double delta_phi = G4RandGauss::shoot(phi_drift, phi_diff); double delta_z = G4RandGauss::shoot(z_drift, z_diff); double phi_gap = G4RandGauss::shoot(phi_2END, sigma_phi_gap); @@ -414,31 +427,35 @@ vector rtpc_HitProcess :: processID(vector id, G4Step* // ------------------------------------------------------------------------------------ - double phi_rad= phi0_rad+delta_phi+phi_gap; //phi at pad pcb board - if( phi_rad<0.0 ) phi_rad+=2.0*PI; - if( phi_rad>=2.0*PI ) phi_rad-=2.0*PI; + double phi_rad= phi0_rad + delta_phi + phi_gap; //phi at pad pcb board + if (phi_rad < 0.0) phi_rad += 2.0*PI; // return to (0, 2PI) + if (phi_rad >= 2.0*PI) phi_rad -= 2.0*PI; // return to (0, 2PI) - double z_pos = (z_cm+delta_z)*10.0; //mm + //double z_pos = (z_cm + delta_z)*10.0; // mm + double z_pos = z_cm*10.0 + delta_z; // mm int col = -999; int row = -999; row = ceil(phi_rad/phi_per_pad); - float z_shift = row%4; + float z_shift = (row - 1)%4; - float col_min = -RTPC_L/2.0+z_shift; - float col_max = RTPC_L/2.0+z_shift; + float col_min = -RTPC_L/2.0 + z_shift; + float col_max = RTPC_L/2.0 + z_shift; - if( z_pos < col_min || z_pos > col_max ) {row = -999; col = -999;} + if (z_pos < col_min || z_pos > col_max) {row = -999; col = -999;} else { - row = ceil(phi_rad/phi_per_pad); - col = ceil((z_pos+RTPC_L/2.0-z_shift)/PAD_L); + col = ceil((z_pos + RTPC_L/2.0 - z_shift)/PAD_L); } - int Num_of_Col = (int) (RTPC_L/PAD_L); - chan = row*Num_of_Col+col; - - yid[0].id=chan; + yid[0].id = row; + yid[1].id = col; + +/* cout << " In processID ***************" << endl; + cout << " input id size: " << id.size() << " name: " << id[0].name << ", " << id[1].name << endl; + cout << "row: " << row << " | " << col << endl; + cout << "id: " << yid[0].id << " | " << yid[1].id << " || rule: " << yid[0].rule << " id_sharing: " << yid[0].id_sharing << " | " << yid[1].id_sharing << " time: " << yid[0].time << " | " << yid[1].time << " name: " << yid[0].name << ", " << yid[1].name << endl; +*/ return yid; } @@ -485,10 +502,22 @@ double rtpc_HitProcess :: voltage(double charge, double time, double forTime) +void rtpc_HitProcess::initWithRunNumber(int runno) +{ + //cout << "*** In initWithRunNumber: runno = " << " " << runno << endl; + string digiVariation = gemcOpt.optMap["DIGITIZATION_VARIATION"].args; + string digiSnapshotTime = gemcOpt.optMap["DIGITIZATION_TIMESTAMP"].args; + if(rtpcc.runNo != runno) { + cout << " > Initializing " << HCname << " digitization for run number " << runno << endl; + rtpcc = initializeRTPCConstants(runno, digiVariation, digiSnapshotTime, accountForHardwareStatus); + rtpcc.runNo = runno; + } +} +// this static function will be loaded first thing by the executable +// setup 11 instead of -1 (in the clas12Tags) because z0, z2, z4 are also used in processID, they need to be loaded first. Moreover, only run 11 in CCDB has appropricate parameters for simulation. +rtpcConstants rtpc_HitProcess::rtpcc = initializeRTPCConstants(-1); - - - - +// add class and member function which can read constants from ccdb in rtpc_hitprocess.X files for the future needed, but comment out them for now (2023.10.03) +// diff --git a/source/hitprocess/clas12/rtpc_hitprocess.h b/source/hitprocess/clas12/rtpc_hitprocess.h index 703976e..1cb31a9 100644 --- a/source/hitprocess/clas12/rtpc_hitprocess.h +++ b/source/hitprocess/clas12/rtpc_hitprocess.h @@ -4,6 +4,31 @@ // gemc headers #include "HitProcess.h" +static const double PI=3.1415926535; + +// constants to be used in the digitization routine +class rtpcConstants +{ +public: + + // database + int runNo; + string date; + string connection; + char database[80]; + + // translation table + TranslationTable TT; + + // add constants here + // drift parameters + double z0[7], z2[7], z4[7]; + // gain balance (example) keep for future or reference + //double gain_balance[96][180]; // 180 rows, 96 cols + +}; + + // Class definition /// \class rtpc_HitProcess @@ -13,6 +38,12 @@ class rtpc_HitProcess : public HitProcess ~rtpc_HitProcess(){;} + // constants initialized with initWithRunNumber + static rtpcConstants rtpcc; + + void initWithRunNumber(int runno); + + // - integrateDgt: returns digitized information integrated over the hit map integrateDgt(MHit*, int); @@ -35,10 +66,14 @@ class rtpc_HitProcess : public HitProcess // - electronicNoise: returns a vector of hits generated / by electronics. vector electronicNoise(); -public: - // RTPC geometry parameters - float PAD_W, PAD_L, PAD_S, RTPC_L; - float phi_per_pad; +private: + + // RTPC geometry parameters (constant) These should be consistant with PadVector.java in coatjava. + const float PAD_W = 2.79; + const float PAD_L = 4.0; + const float PAD_S = 79.0; // old value was 80.0 + const float RTPC_L = 384.0; + const float phi_per_pad = (2.0*PI)/180; // parameters for drift and diffustion equations for drift time, // drift angle, and diffusion in z @@ -47,15 +82,26 @@ class rtpc_HitProcess : public HitProcess float a_z, b_z; // variables for storing drift times and diffusion in time - float t_2GEM2, t_2GEM3, t_2PAD, t_2END; - float sigma_t_2GEM2, sigma_t_2GEM3, sigma_t_2PAD, sigma_t_gap; + //float sigma_t_2GEM2, sigma_t_2GEM3, sigma_t_2PAD, sigma_t_gap; + const float sigma_t_2GEM2 = 8.72728; + const float sigma_t_2GEM3 = 5.62223; + const float sigma_t_2PAD = 7.58056; + float sigma_t_gap; // variables for storing drift angle and diffusion in phi - float phi_2GEM2, phi_2GEM3, phi_2PAD, phi_2END; - float sigma_phi_2GEM2, sigma_phi_2GEM3, sigma_phi_2PAD, sigma_phi_gap; + float phi_2GEM2 = 0.0; // 0.0416925; + float phi_2GEM3 = 0.0; // 0.0416574; + float phi_2PAD = 0.0; // 0.057566; + float phi_2END; + //float sigma_phi_2GEM2, sigma_phi_2GEM3, sigma_phi_2PAD, sigma_phi_gap; + const float sigma_phi_2GEM2 = 0.00384579; + const float sigma_phi_2GEM3 = 0.00160235; + const float sigma_phi_2PAD = 0.00238653; + float sigma_phi_gap; + float z_cm; - float TPC_TZERO; + float TPC_TZERO; // What's this? map timeShift_map; double shift_t; diff --git a/source/materials/material_factory.cc b/source/materials/material_factory.cc index 3653f39..c8085c0 100644 --- a/source/materials/material_factory.cc +++ b/source/materials/material_factory.cc @@ -126,6 +126,18 @@ void material::opticalsFromString(string property, string what) { if (what == "rayleigh") rayleigh.push_back(get_number(trimmedC)); + if (what == "mie") + mie.push_back(get_number(trimmedC)); + + if (what == "mieforward") + mieforward = get_number(trimmedC); + + if (what == "miebackward") + miebackward = get_number(trimmedC); + + if (what == "mieratio") + mieratio = get_number(trimmedC); + if (what == "birkConstant") birkConstant = get_number(trimmedC); } @@ -141,6 +153,7 @@ void material::opticalsFromString(string property, string what) { if (fastcomponent.size() != photonEnergy.size()) fastcomponent.clear(); if (slowcomponent.size() != photonEnergy.size()) slowcomponent.clear(); if (rayleigh.size() != photonEnergy.size()) rayleigh.clear(); + if (mie.size() != photonEnergy.size()) mie.clear(); } } @@ -369,6 +382,15 @@ map materials::materialsFromMap(map mma optTable.back()->AddProperty("RAYLEIGH", penergy, ray, nopts); } + // mie scattering + if (it->second.mie.size() == nopts) { + double mie[nopts]; + for (unsigned i = 0; i < nopts; i++) + mie[i] = it->second.mie[i]; + + optTable.back()->AddProperty("MIEHG", penergy, mie, nopts); + } + // in the API -1 is the default @@ -397,6 +419,16 @@ map materials::materialsFromMap(map mma if (it->second.birkConstant != -1) mats[it->first]->GetIonisation()->SetBirksConstant(it->second.birkConstant); + // additional mie scattering constants + if (it->second.miebackward != -1) + optTable.back()->AddConstProperty("MIEHG_BACKWARD", it->second.miebackward); + + if (it->second.mieforward != -1) + optTable.back()->AddConstProperty("MIEHG_FORWARD", it->second.mieforward); + + if (it->second.mieratio != -1) + optTable.back()->AddConstProperty("MIEHG_FORWARD_RATIO", it->second.mieratio); + mats[it->first]->SetMaterialPropertiesTable(optTable.back()); } diff --git a/source/materials/material_factory.h b/source/materials/material_factory.h index 20f0eae..b61bedd 100644 --- a/source/materials/material_factory.h +++ b/source/materials/material_factory.h @@ -33,6 +33,12 @@ class material vector absorptionLength; vector reflectivity; vector efficiency; + + // mie scattering + vector mie; + double mieforward; + double miebackward; + double mieratio; // scintillation vector fastcomponent; diff --git a/source/materials/mysql_materials.cc b/source/materials/mysql_materials.cc index 8e80db8..96c01b7 100644 --- a/source/materials/mysql_materials.cc +++ b/source/materials/mysql_materials.cc @@ -52,7 +52,8 @@ map mysql_materials::initMaterials(runConditions rc, goptio string dbexecute = "select name, description, density, ncomponents, components, photonEnergy, indexOfRefraction, "; dbexecute += "absorptionLength, reflectivity, efficiency, fastcomponent, slowcomponent, "; - dbexecute += "scintillationyield, resolutionscale, fasttimeconstant, slowtimeconstant, yieldratio from " + tname; + dbexecute += "scintillationyield, resolutionscale, fasttimeconstant, slowtimeconstant, yieldratio, rayleigh, birkConstant, "; + dbexecute += "mie, mieforward, miebackward, mieratio from " + tname; dbexecute += " where variation ='" + variation + "'"; // executing query - will exit if not successfull. @@ -94,6 +95,12 @@ map mysql_materials::initMaterials(runConditions rc, goptio // Birk Constant thisMat.birkConstant = q.value(18).toDouble(); // yieldratio + // Mie scattering + thisMat.opticalsFromString(qv_tostring( q.value(19)), "mie"); + thisMat.mieforward = q.value(20).toDouble(); // mie forward scattering + thisMat.miebackward = q.value(21).toDouble(); // mie backward scattering + thisMat.mieratio = q.value(22).toDouble(); // mie forward/backward scattering ratio + mymats[thisMat.name] = thisMat; } diff --git a/source/materials/sqlite_materials.cc b/source/materials/sqlite_materials.cc index 335c13c..b10edc2 100644 --- a/source/materials/sqlite_materials.cc +++ b/source/materials/sqlite_materials.cc @@ -52,7 +52,8 @@ map sqlite_materials::initMaterials(runConditions rc, gopt string dbexecute = "select name, description, density, ncomponents, components, photonEnergy, indexOfRefraction, "; dbexecute += "absorptionLength, reflectivity, efficiency, fastcomponent, slowcomponent, "; - dbexecute += "scintillationyield, resolutionscale, fasttimeconstant, slowtimeconstant, yieldratio from materials"; + dbexecute += "scintillationyield, resolutionscale, fasttimeconstant, slowtimeconstant, yieldratio, rayleigh, birkconstant, "; // rayleigh and birk constant were not here? + dbexecute += "mie, mieforward, miebackward, mieratio from materials"; dbexecute += " where variation ='" + variation + "'"; dbexecute += " and run = " + stringify(run_number); dbexecute += " and system = '" + dname + "'"; @@ -70,7 +71,7 @@ map sqlite_materials::initMaterials(runConditions rc, gopt } while (q.next()) { - material thisMat(trimSpacesFromString(qv_tostring( q.value(0)))); // name + material thisMat(trimSpacesFromString(qv_tostring( q.value(0)))); // name thisMat.desc = qv_tostring(q.value(1)); // description thisMat.density = q.value(2).toDouble(); // density thisMat.ncomponents = q.value(3).toInt(); // number of components @@ -94,6 +95,12 @@ map sqlite_materials::initMaterials(runConditions rc, gopt // Birk Constant thisMat.birkConstant = q.value(18).toDouble(); + // Mie scattering + thisMat.opticalsFromString(qv_tostring(q.value(19)), "mie"); + thisMat.mieforward = q.value(20).toDouble(); + thisMat.miebackward = q.value(21).toDouble(); + thisMat.mieratio = q.value(22).toDouble(); + mymats[thisMat.name] = thisMat; } diff --git a/source/materials/text_materials.cc b/source/materials/text_materials.cc index b0aeba8..056e995 100644 --- a/source/materials/text_materials.cc +++ b/source/materials/text_materials.cc @@ -110,7 +110,16 @@ map text_materials::initMaterials(runConditions rc, goptio if (gt.data.size() > 18) { thisMat.birkConstant = get_number(gt.data[18]); } - + // this condition is for backward compatibility, + // Birk Constant were added with gemc 2.6 + // 23 exact quantities + if (gt.data.size() > 19) { + thisMat.opticalsFromString(gt.data[19], "mie"); + thisMat.mieforward = get_number(gt.data[20]); + thisMat.miebackward = get_number(gt.data[21]); + thisMat.mieratio = get_number(gt.data[22]); + } + mymats[thisMat.name] = thisMat; } diff --git a/source/output/hipoSchemas.cc b/source/output/hipoSchemas.cc index cbe0dfb..da105ae 100644 --- a/source/output/hipoSchemas.cc +++ b/source/output/hipoSchemas.cc @@ -189,6 +189,8 @@ HipoSchema :: HipoSchema() schemasToLoad["LTCC::adc"] = ltccADCSchema; schemasToLoad["LTCC::tdc"] = ltccTDCSchema; schemasToLoad["RICH::tdc"] = richTDCSchema; + schemasToLoad["RTPC::adc"] = rtpcADCSchema; + schemasToLoad["RTPC::pos"] = rtpcPOSSchema; schemasToLoad["HEL::flip"] = helFLIPSchema; schemasToLoad["RASTER::adc"] = rasterADCSchema; schemasToLoad["URWELL::adc"] = urwellADCSchema; diff --git a/source/physics/PhysicsList.cc b/source/physics/PhysicsList.cc index ebf4953..69a3fed 100644 --- a/source/physics/PhysicsList.cc +++ b/source/physics/PhysicsList.cc @@ -386,7 +386,10 @@ void PhysicsList::cookPhysics() // see G4OpticalPhysics.hh G4OpticalPhysics* opticalPhysics = new G4OpticalPhysics(); - + // enable Mie scattering (do not think it is enabled by default) + // make this optional? + opticalPhysics->Configure(kMieHG, true); + // this was deprecated? // Method G4OpticalPhysics::SetWLSTimeProfile is deprecated. // Use G4OpticalParameters::SetWLSTimeProfile(G4String) instead.