From 84def80237527854b759e68386f7848c4b909a28 Mon Sep 17 00:00:00 2001 From: raffaelladevita Date: Tue, 28 Jan 2025 14:12:28 -0500 Subject: [PATCH 01/12] added ftof and ctof fadc time offset --- source/hitprocess/clas12/ctof_hitprocess.cc | 18 +++++++++++++++--- source/hitprocess/clas12/ctof_hitprocess.h | 3 +++ source/hitprocess/clas12/ftof_hitprocess.cc | 16 ++++++++++++++-- source/hitprocess/clas12/ftof_hitprocess.h | 3 +++ 4 files changed, 35 insertions(+), 5 deletions(-) diff --git a/source/hitprocess/clas12/ctof_hitprocess.cc b/source/hitprocess/clas12/ctof_hitprocess.cc index b3e35e75..31b78bcf 100644 --- a/source/hitprocess/clas12/ctof_hitprocess.cc +++ b/source/hitprocess/clas12/ctof_hitprocess.cc @@ -168,7 +168,18 @@ static ctofConstants initializeCTOFConstants(int runno, string digiVariation = " ctc.tdcconv[isec - 1][ilay - 1][1].push_back(data[row][4]); } - + snprintf(ctc.database, sizeof(ctc.database), "/calibration/ctof/fadc_offset:%d:%s%s", ctc.runNo, digiVariation.c_str(), timestamp.c_str()); + cout << "CTOF:Getting adc_offsets" << endl; + data.clear(); + calib->GetCalib(data, ctc.database); + for (unsigned row = 0; row < data.size(); row++) { + isec = data[row][0]; + ilay = data[row][1]; + // istr = data[row][2]; + ctc.adcoffset[isec - 1][ilay - 1][0].push_back(data[row][3]); + ctc.adcoffset[isec - 1][ilay - 1][1].push_back(data[row][4]); + } + snprintf(ctc.database, sizeof(ctc.database), "/geometry/ctof/ctof:%d:%s%s", ctc.runNo, digiVariation.c_str(), timestamp.c_str()); cout << "CTOF:Getting geometry" << endl; data.clear(); @@ -264,6 +275,7 @@ map ctof_HitProcess::integrateDgt(MHit* aHit, int hitn) // TDC conversion factors double tdcconv = ctc.tdcconv[sector - 1][layer - 1][side][paddle - 1]; + double adcoffset = ctc.adcoffset[sector - 1][layer - 1][side][paddle - 1]; double time_in_ns = 0; if(aHit->isBackgroundHit == 1) { @@ -279,7 +291,7 @@ map ctof_HitProcess::integrateDgt(MHit* aHit, int hitn) dgtz["component"] = paddle; dgtz["ADC_order"] = side; dgtz["ADC_ADC"] = (int) adc; - dgtz["ADC_time"] = convert_to_precision(stepTime); + dgtz["ADC_time"] = convert_to_precision(stepTime - adcoffset); dgtz["ADC_ped"] = 0; dgtz["TDC_order"] = side + 2; @@ -388,7 +400,7 @@ map ctof_HitProcess::integrateDgt(MHit* aHit, int hitn) } // standardizing fadc time and tdc info - double fadc_time = convert_to_precision(time_in_ns); + double fadc_time = convert_to_precision(time_in_ns - adcoffset); dgtz["hitn"] = hitn; diff --git a/source/hitprocess/clas12/ctof_hitprocess.h b/source/hitprocess/clas12/ctof_hitprocess.h index 5ed5e9eb..6bd5fcf4 100644 --- a/source/hitprocess/clas12/ctof_hitprocess.h +++ b/source/hitprocess/clas12/ctof_hitprocess.h @@ -38,6 +38,9 @@ class ctofConstants // tdc_conc: tdc conversion factors vector tdcconv[1][1][2]; + // adc time offsets + vector adcoffset[1][1][2]; + // veff: effective velocity vector veff[1][1][2]; diff --git a/source/hitprocess/clas12/ftof_hitprocess.cc b/source/hitprocess/clas12/ftof_hitprocess.cc index c363d727..0be111cb 100644 --- a/source/hitprocess/clas12/ftof_hitprocess.cc +++ b/source/hitprocess/clas12/ftof_hitprocess.cc @@ -164,6 +164,17 @@ static ftofConstants initializeFTOFConstants(int runno, string digiVariation = " ftc.tdcconv[isec - 1][ilay - 1][1].push_back(data[row][4]); } + cout << "FTOF:Getting adc_offsets" << endl; + snprintf(ftc.database, sizeof(ftc.database), "/calibration/ftof/fadc_offset:%d:%s%s", ftc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); + calib->GetCalib(data, ftc.database); + for (unsigned row = 0; row < data.size(); row++) { + isec = data[row][0]; + ilay = data[row][1]; + ftc.adcoffset[isec - 1][ilay - 1][0].push_back(data[row][3]); + ftc.adcoffset[isec - 1][ilay - 1][1].push_back(data[row][4]); + } + cout << "FTOF:Getting resolutions" << endl; snprintf(ftc.database, sizeof(ftc.database), "/calibration/ftof/tres:%d:%s%s", ftc.runNo, digiVariation.c_str(), timestamp.c_str()); data.clear(); @@ -249,6 +260,7 @@ map ftof_HitProcess::integrateDgt(MHit* aHit, int hitn) { int pmt = identity[3].id; // 0=> Left PMT, 1=> Right PMT. A better name would be pmtSide // TDC conversion factors + double adcoffset = ftc.adcoffset[sector - 1][panel - 1][pmt][paddle - 1]; double tdcconv = ftc.tdcconv[sector - 1][panel - 1][pmt][paddle - 1]; double time_in_ns = 0; @@ -265,7 +277,7 @@ map ftof_HitProcess::integrateDgt(MHit* aHit, int hitn) { dgtz["component"] = paddle; dgtz["ADC_order"] = pmt; dgtz["ADC_ADC"] = (int) adc; - dgtz["ADC_time"] = convert_to_precision(stepTime); + dgtz["ADC_time"] = convert_to_precision(stepTime - adcoffset); dgtz["ADC_ped"] = 0; dgtz["TDC_order"] = pmt + 2; @@ -396,7 +408,7 @@ map ftof_HitProcess::integrateDgt(MHit* aHit, int hitn) { // cout << " > FTOF status: " << ftc.status[sector-1][panel-1][0][paddle-1] << " for sector " << sector << ", panel " << panel << ", paddle " << paddle << " left: " << adcl << endl; // cout << " > FTOF status: " << ftc.status[sector-1][panel-1][1][paddle-1] << " for sector " << sector << ", panel " << panel << ", paddle " << paddle << " right: " << adcr << endl; - double fadc_time = convert_to_precision(time_in_ns); + double fadc_time = convert_to_precision(time_in_ns - ftc.adcoffset[sector - 1][panel - 1][pmt][paddle - 1]); tdc = time_in_ns / tdcconv; dgtz["sector"] = sector; diff --git a/source/hitprocess/clas12/ftof_hitprocess.h b/source/hitprocess/clas12/ftof_hitprocess.h index 9dd8ed90..63a58666 100644 --- a/source/hitprocess/clas12/ftof_hitprocess.h +++ b/source/hitprocess/clas12/ftof_hitprocess.h @@ -36,6 +36,9 @@ class ftofConstants // tdc_conc: tdc conversion factors vector tdcconv[6][3][2]; + // adc time offsets + vector adcoffset[6][3][2]; + // veff: effective velocity vector veff[6][3][2]; From 457e81811685c03db3f050fbc7c891ce9772ceaf Mon Sep 17 00:00:00 2001 From: raffaelladevita Date: Tue, 28 Jan 2025 14:14:20 -0500 Subject: [PATCH 02/12] added htcc npe -> adc conversion factors from ccdb, removing the hardcoded gain=100 --- source/hitprocess/clas12/htcc_hitprocess.cc | 12 ++++++++++-- source/hitprocess/clas12/htcc_hitprocess.h | 6 ++++-- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/source/hitprocess/clas12/htcc_hitprocess.cc b/source/hitprocess/clas12/htcc_hitprocess.cc index 9118aa6f..3a84d310 100644 --- a/source/hitprocess/clas12/htcc_hitprocess.cc +++ b/source/hitprocess/clas12/htcc_hitprocess.cc @@ -52,8 +52,16 @@ static htccConstants initializeHTCCConstants(int runno, string digiVariation = " isec = data[row][0]; ilay = data[row][1]; htccc.status[isec-1][ilay-1].push_back(data[row][3]); } + } + + cout<<"HTCC:Getting gain"<GetCalib(data,htccc.database); + for(unsigned row = 0; row < data.size(); row++) { + isec = data[row][0]; ilay = data[row][1]; + htccc.gain[isec-1][ilay-1].push_back(data[row][3]); + } - } cout<<"HTCC:Getting mc_gain"<GetCalib(data,htccc.database); @@ -287,7 +295,7 @@ map htcc_HitProcess :: integrateDgt(MHit* aHit, int hitn) } } - double adc = 100 * G4RandGauss::shoot(ndetected*htccc.mc_gain[idsector-1][idhalf-1][idring-1], ndetected*htccc.mc_smear[idsector-1][idhalf-1][idring-1]); + double adc = htccc.gain[idsector-1][idhalf-1][idring-1] * G4RandGauss::shoot(ndetected*htccc.mc_gain[idsector-1][idhalf-1][idring-1], ndetected*htccc.mc_smear[idsector-1][idhalf-1][idring-1]); time_in_ns = tInfos.time + htccc.tshift[idsector-1][idhalf-1][idring-1]; double fadc_time = convert_to_precision(time_in_ns); diff --git a/source/hitprocess/clas12/htcc_hitprocess.h b/source/hitprocess/clas12/htcc_hitprocess.h index 4c39ab13..20d02377 100644 --- a/source/hitprocess/clas12/htcc_hitprocess.h +++ b/source/hitprocess/clas12/htcc_hitprocess.h @@ -35,9 +35,11 @@ class htccConstants // veff: time shift vector tshift[6][2]; + // nephe to ADC conversion factors + vector gain[6][2]; + // gain to match data yield to MC - vector mc_gain[6][2]; - + vector mc_gain[6][2]; // smearing factor to match data yield to MC vector mc_smear[6][2]; From 8b01e783fe14a155b257ec326ca7a50588557882 Mon Sep 17 00:00:00 2001 From: raffaelladevita Date: Fri, 7 Nov 2025 14:20:37 -0500 Subject: [PATCH 03/12] added TDC jitter, OF adc-tdc offsets, updated TW for FTOF and added hpos for CTOF --- source/hitprocess/clas12/band_hitprocess.cc | 14 ++++- source/hitprocess/clas12/band_hitprocess.h | 5 ++ source/hitprocess/clas12/cnd_hitprocess.cc | 22 ++++++-- source/hitprocess/clas12/cnd_hitprocess.h | 5 ++ source/hitprocess/clas12/ctof_hitprocess.cc | 57 ++++++++++++++------- source/hitprocess/clas12/ctof_hitprocess.h | 11 +++- source/hitprocess/clas12/dc_hitprocess.cc | 10 +++- source/hitprocess/clas12/dc_hitprocess.h | 5 ++ source/hitprocess/clas12/ecal_hitprocess.cc | 11 +++- source/hitprocess/clas12/ecal_hitprocess.h | 5 ++ source/hitprocess/clas12/ftof_hitprocess.cc | 52 ++++++++++++++++--- source/hitprocess/clas12/ftof_hitprocess.h | 13 ++++- source/hitprocess/clas12/htcc_hitprocess.cc | 12 ++++- source/hitprocess/clas12/htcc_hitprocess.h | 3 ++ 14 files changed, 186 insertions(+), 39 deletions(-) diff --git a/source/hitprocess/clas12/band_hitprocess.cc b/source/hitprocess/clas12/band_hitprocess.cc index d2008394..e34f1019 100644 --- a/source/hitprocess/clas12/band_hitprocess.cc +++ b/source/hitprocess/clas12/band_hitprocess.cc @@ -93,6 +93,14 @@ static bandHitConstants initializeBANDHitConstants(int runno, string digiVariati calib->GetCalib(data, bhc.database); bhc.tdcconv = data[0][3]; + snprintf(bhc.database, sizeof(bhc.database), "/calibration/band/time_jitter:%d:%s%s", bhc.runNo, digiVariation.c_str(), timestamp.c_str()); + cout << "BAND:Getting time_jitter" << endl; + data.clear(); + calib->GetCalib(data, bhc.database); + bhc.jitter_period = data[0][3]; + bhc.jitter_phase = data[0][4]; + bhc.jitter_cycles = data[0][5]; + // These are not in the CCDB // Fill with constant values //cout<<"BAND:Getting MeV->ADC conversions"< band_HitProcess :: integrateDgt(MHit* aHit, int hitn) int ADC = (int) ( adcFactor*(side == 0 ? eTotL : eTotR) ); int amplitude = (int) ( adcFactor*(side == 0 ? xHit : zHit) ); time_in_ns = (side == 0 ? tL_fadc : tR_fadc); - double fadc_time = convert_to_precision(time_in_ns); - int TDC = (int) ( adcFactor*(side == 0 ? tL_tdc : tR_tdc)/tdcconv ); + double fadc_time = convert_to_precision(time_in_ns); + + double tdc_jitter = bhc.jitter_period * ((0 + bhc.jitter_phase) % bhc.jitter_cycles); // assumes event timestamp is zero + int TDC = (int) ( adcFactor*((side == 0 ? tL_tdc : tR_tdc)+tdc_jitter)/tdcconv ); dgtz["hitn"] = (int) hitn; dgtz["sector"] = (int) sector; diff --git a/source/hitprocess/clas12/band_hitprocess.h b/source/hitprocess/clas12/band_hitprocess.h index 45a6632c..56c9810e 100644 --- a/source/hitprocess/clas12/band_hitprocess.h +++ b/source/hitprocess/clas12/band_hitprocess.h @@ -20,6 +20,11 @@ class bandHitConstants int nlayer; int ncomp; + // tdc jitter parameters + double jitter_period; + int jitter_phase; + int jitter_cycles; + double mev_adc[6][6][7]; double eff_vel_tdc[6][6][7]; double eff_vel_fadc[6][6][7]; diff --git a/source/hitprocess/clas12/cnd_hitprocess.cc b/source/hitprocess/clas12/cnd_hitprocess.cc index 10317bb5..5cb146f0 100644 --- a/source/hitprocess/clas12/cnd_hitprocess.cc +++ b/source/hitprocess/clas12/cnd_hitprocess.cc @@ -75,6 +75,14 @@ static cndConstants initializeCNDConstants(int runno, string digiVariation = "de cndc.slope_R[isec-1][ilay-1][istr-1]=data[row][5]; } + snprintf(cndc.database, sizeof(cndc.database), "/calibration/cnd/time_jitter:%d:%s%s", cndc.runNo, digiVariation.c_str(), timestamp.c_str()); + cout << "CND:Getting time_jitter" << endl; + data.clear(); + calib->GetCalib(data, cndc.database); + cndc.jitter_period = data[0][3]; + cndc.jitter_phase = data[0][4]; + cndc.jitter_cycles = data[0][5]; + cout<<"CND:Getting attenuation"<GetCalib(data,cndc.database); @@ -281,7 +289,10 @@ map cnd_HitProcess :: integrateDgt(MHit* aHit, int hitn) double threshold = 0; double time_in_ns = 0; - + + double tdc_jitter = cndc.jitter_period * ((0 + cndc.jitter_phase) % cndc.jitter_cycles); // assumes event timestamp is zero + + // LEFT PADDLE if ( side == 1 ){ @@ -455,14 +466,14 @@ map cnd_HitProcess :: integrateDgt(MHit* aHit, int hitn) // MARK: TO delete later if (etotUp > 0.) { time_in_ns = G4RandGauss::shoot(timeD, sigmaTD/sqrt(etotUp)); - TDCD = (int) ( (G4RandGauss::shoot(timeD, sigmaTD/sqrt(etotUp)) ) / slope_D); + TDCD = (int) ( (G4RandGauss::shoot(timeD+tdc_jitter, sigmaTD/sqrt(etotUp)) ) / slope_D); double npheD = G4Poisson(etotUp*pmtPEYldD); double eneD = npheD/pmtPEYldD; ADCD = (int) (eneD*adc_mip_D*2./(dEdxMIP*thickness)); } if (etotDown > 0.) { time_in_ns = G4RandGauss::shoot(timeN, sigmaTD/sqrt(etotDown)); - TDCN = (int) ( (G4RandGauss::shoot(timeN, sigmaTN/sqrt(etotDown)) ) / slope_N); + TDCN = (int) ( (G4RandGauss::shoot(timeN+tdc_jitter, sigmaTN/sqrt(etotDown)) ) / slope_N); double npheN = G4Poisson(etotDown*pmtPEYldN); double eneN = npheN/pmtPEYldN; ADCN = (int) (eneN*adc_mip_N*2./(dEdxMIP*thickness)); @@ -472,7 +483,7 @@ map cnd_HitProcess :: integrateDgt(MHit* aHit, int hitn) // of the shooting of random numbers if ( eTotal > 0 ) { time_in_ns = G4RandGauss::shoot(eTime, sigma/sqrt(eTotal)); - TDC = (int) ( (G4RandGauss::shoot(eTime, sigma/sqrt(eTotal)) ) / slope); + TDC = (int) ( (G4RandGauss::shoot(eTime+tdc_jitter, sigma/sqrt(eTotal)) ) / slope); double nphe = G4Poisson(eTotal*pmtPEYld); double ene = nphe/pmtPEYld; ADC = (int) (ene*adc_mip*2./(dEdxMIP*thickness)); @@ -583,12 +594,13 @@ map cnd_HitProcess :: integrateDgt(MHit* aHit, int hitn) // standardizing fadc time and tdc info double fadc_time = convert_to_precision(time_in_ns); - + dgtz["hitn"] = hitn; dgtz["sector"] = sector; dgtz["layer"] = layer; dgtz["component"] = 1; dgtz["ADC_order"] = adc_order ; // 0 = left 1 = right + dgtz["ADC_ADC"] = ADC; dgtz["ADC_time"] = fadc_time; dgtz["ADC_ped"] = 0; diff --git a/source/hitprocess/clas12/cnd_hitprocess.h b/source/hitprocess/clas12/cnd_hitprocess.h index e38661d5..7758bca7 100644 --- a/source/hitprocess/clas12/cnd_hitprocess.h +++ b/source/hitprocess/clas12/cnd_hitprocess.h @@ -41,6 +41,11 @@ class cndConstants double ecalN[24][3][2]; */ + // tdc jitter parameters + double jitter_period; + int jitter_phase; + int jitter_cycles; + int status_L[24][3][2]; int status_R[24][3][2]; double threshold_L[24][3][2]; diff --git a/source/hitprocess/clas12/ctof_hitprocess.cc b/source/hitprocess/clas12/ctof_hitprocess.cc index 31b78bcf..453a666b 100644 --- a/source/hitprocess/clas12/ctof_hitprocess.cc +++ b/source/hitprocess/clas12/ctof_hitprocess.cc @@ -76,7 +76,7 @@ static ctofConstants initializeCTOFConstants(int runno, string digiVariation = " } if(accountForHardwareStatus) { - snprintf(ctc.database, sizeof(ctc.database), "/calibration/ctof/:%d:%s%s", ctc.runNo, digiVariation.c_str(), timestamp.c_str()); + snprintf(ctc.database, sizeof(ctc.database), "/calibration/ctof/status:%d:%s%s", ctc.runNo, digiVariation.c_str(), timestamp.c_str()); cout << "CTOF:Getting status" << endl; data.clear(); calib->GetCalib(data, ctc.database); @@ -113,8 +113,6 @@ static ctofConstants initializeCTOFConstants(int runno, string digiVariation = " ctc.efficiency[isec - 1][ilay - 1][1].push_back(data[row][3]); } - - snprintf(ctc.database, sizeof(ctc.database), "/calibration/ctof/gain_balance:%d:%s%s", ctc.runNo, digiVariation.c_str(), timestamp.c_str()); cout << "CTOF:Getting gain_balance" << endl; data.clear(); @@ -127,23 +125,30 @@ static ctofConstants initializeCTOFConstants(int runno, string digiVariation = " ctc.countsForMIP[isec - 1][ilay - 1][1].push_back(data[row][4]); } - /* For future use in HitProcess - cout<<"Getting time_walk"<GetCalib(data,ctc.database); for(unsigned row = 0; row < data.size(); row++) { - isec = data[row][0]; ilay = data[row][1]; istr = data[row][2]; - ctc.twlk[isec-1][ilay-1][0].push_back(data[row][3]); - ctc.twlk[isec-1][ilay-1][1].push_back(data[row][4]); - ctc.twlk[isec-1][ilay-1][2].push_back(data[row][5]); - ctc.twlk[isec-1][ilay-1][3].push_back(data[row][6]); - ctc.twlk[isec-1][ilay-1][4].push_back(data[row][7]); - ctc.twlk[isec-1][ilay-1][5].push_back(data[row][8]); + isec = data[row][0]; + ilay = data[row][1]; + ctc.hpos[isec-1][ilay-1][0].push_back(data[row][3]); + ctc.hpos[isec-1][ilay-1][1].push_back(data[row][4]); } - */ - snprintf(ctc.database, sizeof(ctc.database), "/calibration/ctof/time_offsets:%d:%s%s", ctc.runNo, digiVariation.c_str(), timestamp.c_str()); + snprintf(ctc.database, sizeof(ctc.database), "/calibration/ctof/hposbin:%d:%s%s",ctc.runNo, digiVariation.c_str(), timestamp.c_str()); + cout<<"CTOF: Getting hpos parameters"<GetCalib(data,ctc.database); + for(unsigned row = 0; row < data.size(); row++) + { + isec = data[row][0]; + ilay = data[row][1]; + for(unsigned bin=0; bin<=ctc.hposbins; bin++) { + ctc.hpos[isec-1][ilay-1][bin].push_back(data[row][3+bin]); + } + } + + snprintf(ctc.database, sizeof(ctc.database), "/calibration/ctof/time_offsets:%d:%s%s", ctc.runNo, digiVariation.c_str(), timestamp.c_str()); cout << "CTOF:Getting time_offsets" << endl; data.clear(); calib->GetCalib(data, ctc.database); @@ -168,6 +173,14 @@ static ctofConstants initializeCTOFConstants(int runno, string digiVariation = " ctc.tdcconv[isec - 1][ilay - 1][1].push_back(data[row][4]); } + snprintf(ctc.database, sizeof(ctc.database), "/calibration/ctof/time_jitter:%d:%s%s", ctc.runNo, digiVariation.c_str(), timestamp.c_str()); + cout << "CTOF:Getting time_jitter" << endl; + data.clear(); + calib->GetCalib(data, ctc.database); + ctc.jitter_period = data[0][3]; + ctc.jitter_phase = data[0][4]; + ctc.jitter_cycles = data[0][5]; + snprintf(ctc.database, sizeof(ctc.database), "/calibration/ctof/fadc_offset:%d:%s%s", ctc.runNo, digiVariation.c_str(), timestamp.c_str()); cout << "CTOF:Getting adc_offsets" << endl; data.clear(); @@ -310,7 +323,8 @@ map ctof_HitProcess::integrateDgt(MHit* aHit, int hitn) // ctof paddle center is offsetby ctc.zoffset from the CLAS12 target position, // so need to z is also the local coordinate // side = 0 or 1, - double d = length + (1. - 2. * side)*(tInfos.z - offset - ctc.targetZPos); // The distance between the hit and PMT? + double y = tInfos.z - offset - ctc.targetZPos; + double d = length + (1. - 2. * side)*y; // The distance between the hit and PMT? // attenuation length double attlen = ctc.attlen[sector - 1][layer - 1][side][paddle - 1]; @@ -362,16 +376,21 @@ map ctof_HitProcess::integrateDgt(MHit* aHit, int hitn) //double B = ctc.twlk[sector-1][layer-1][1][paddle-1]; //double C = ctc.twlk[sector-1][layer-1][2][paddle-1]; //double timeWalkUp = A/(B+C*sqrt(adcu)); - double timeWalk = 0.; + double hpos = ctc.hpos[sector-1][layer-1][0][paddle-1] * exp(ctc.hpos[sector-1][layer-1][1][paddle-1]/cm * y); + if(ctc.hposbin[sector-1][layer-1][0][paddle-1]>0) { + int bin = (int) floor(y/cm) + ctc.hposbins/2; + hpos += ctc.hposbin[sector-1][layer-1][bin+1][paddle-1]; + } double tU = tInfos.time + d/ctc.veff[sector-1][layer-1][side][paddle-1]/cm + (1. - 2. * side)*ctc.toff_UD[sector-1][layer-1][paddle-1]/2. - ctc.toff_RFpad[sector-1][layer-1][paddle-1] - ctc.toff_P2P[sector-1][layer-1][paddle-1] - + timeWalk; + - hpos; time_in_ns = G4RandGauss::shoot(tU, sqrt(2) * ctc.tres[paddle - 1]); // tdcu = tU / tdcconv; - tdc = (int) ( time_in_ns / tdcconv ); + double tdc_jitter = ctc.jitter_period * ((0 + ctc.jitter_phase) % ctc.jitter_cycles); // assumes event timestamp is zero + tdc = (int) ((time_in_ns + tdc_jitter) / tdcconv); } if(accountForHardwareStatus) { diff --git a/source/hitprocess/clas12/ctof_hitprocess.h b/source/hitprocess/clas12/ctof_hitprocess.h index 6bd5fcf4..385e4d85 100644 --- a/source/hitprocess/clas12/ctof_hitprocess.h +++ b/source/hitprocess/clas12/ctof_hitprocess.h @@ -41,6 +41,11 @@ class ctofConstants // adc time offsets vector adcoffset[1][1][2]; + // tdc jitter parameters + double jitter_period; + int jitter_phase; + int jitter_cycles; + // veff: effective velocity vector veff[1][1][2]; @@ -50,8 +55,10 @@ class ctofConstants // countsForMIP: Desired ADC channel for MIP peak calibration vector countsForMIP[1][1][2]; - // twlk: Time walk correction, 3 constants each for L and R - vector twlk[1][1][6]; + // hpos parameters + vector hpos[1][1][2]; + vector hposbin[1][1][101]; + unsigned hposbins = 100; vector toff_UD[1][1]; vector toff_RFpad[1][1]; diff --git a/source/hitprocess/clas12/dc_hitprocess.cc b/source/hitprocess/clas12/dc_hitprocess.cc index f50f13f9..1c09f6a3 100644 --- a/source/hitprocess/clas12/dc_hitprocess.cc +++ b/source/hitprocess/clas12/dc_hitprocess.cc @@ -139,6 +139,12 @@ static dcConstants initializeDCConstants(int runno, string digiVariation = "defa } //******************************************** + snprintf(dcc.database, sizeof(dcc.database), "/calibration/dc/time_jitter:%d:%s%s", dcc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); + calib->GetCalib(data, dcc.database); + dcc.jitter_period = data[0][3]; + dcc.jitter_phase = data[0][4]; + dcc.jitter_cycles = data[0][5]; // reading DC core parameters @@ -250,6 +256,8 @@ map dc_HitProcess :: integrateDgt(MHit* aHit, int hitn) double signal_t = 0; double hit_signal_t = 0; double prop_t = 0; + double tdc_jitter = dcc.jitter_period * ((0 + dcc.jitter_phase) % dcc.jitter_cycles); // assumes event timestamp is zero + for(unsigned int s=0; s dc_HitProcess :: integrateDgt(MHit* aHit, int hitn) // Now calculate the smeared time: // adding the time of hit from the start of the event (signal_t), which also has the drift velocity into it - double smeared_time = unsmeared_time + dt_random + hit_signal_t + prop_t + dcc.get_T0(SECI, SLI, LAYI, nwire); + double smeared_time = unsmeared_time + dt_random + hit_signal_t + prop_t + dcc.get_T0(SECI, SLI, LAYI, nwire) + tdc_jitter; // cout << " DC TIME stime: " << smeared_time << " X: " << X << " doca: " << doca/cm << " dmax: " << dcc.dmaxsuperlayer[SLI] << " tmax: " << dcc.tmaxsuperlayer[SECI][SLI] << " alpha: " << alpha << " thisMgnf: " << thisMgnf << " SECI: " << SECI << " SLI: " << SLI << endl; diff --git a/source/hitprocess/clas12/dc_hitprocess.h b/source/hitprocess/clas12/dc_hitprocess.h index 592654a1..6348bed5 100644 --- a/source/hitprocess/clas12/dc_hitprocess.h +++ b/source/hitprocess/clas12/dc_hitprocess.h @@ -46,6 +46,11 @@ class dcConstants double vmid[6][6], R[6][6]; double dmaxsuperlayer[6]; + // tdc jitter parameters + double jitter_period; + int jitter_phase; + int jitter_cycles; + // sector, SL, slot, cable double T0Correction[6][6][7][6]; diff --git a/source/hitprocess/clas12/ecal_hitprocess.cc b/source/hitprocess/clas12/ecal_hitprocess.cc index 858d613f..204b1f02 100644 --- a/source/hitprocess/clas12/ecal_hitprocess.cc +++ b/source/hitprocess/clas12/ecal_hitprocess.cc @@ -122,6 +122,12 @@ static ecConstants initializeECConstants(int runno, string digiVariation = "defa ecc.global_time_walk[isec-1][ilay-1].push_back(data[row][3]); } + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/time_jitter:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); calib->GetCalib(data, ecc.database); + ecc.jitter_period = data[0][3]; + ecc.jitter_phase = data[0][4]; + ecc.jitter_cycles = data[0][5]; + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/fadc_offset:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); data.clear(); calib->GetCalib(data, ecc.database); for(unsigned row = 0; row < data.size(); row++) @@ -402,6 +408,8 @@ map ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) double dtres2 = ecc.dtres[sector-1][layer-1][2][0]; double dtres3 = ecc.dtres[sector-1][layer-1][3][0]; + double tdc_jitter = ecc.jitter_period * ((0 + ecc.jitter_phase) % ecc.jitter_cycles); // assumes event timestamp is zero + for(unsigned int s=0; s ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) double ftim = (fa4==0||fa6==0) ? 0 : fa0 + fa2+exp(-(radc-fa3)/fa4)+1-exp( (radc-fa5)/fa6); double dtim = (da4==0||da6==0||da7==0) ? 0 : da0 + gtw/radc +da2+exp(-(radc-da3)/da4)+1-exp(-(da5-radc)/da6)-exp(-(radc-da3*0.95)/da7)*pow(radc,da8); ftime_in_ns = FTIME_raw + ftim + tgo - FTOFFSET - tmf - fo; - dtime_in_ns = DTIME_raw + dtim + tgo; + dtime_in_ns = DTIME_raw + dtim + tgo + tdc_jitter; ftime_in_ns_res = G4RandGauss::shoot(ftime_in_ns,getTRES(ADC,ftres0,ftres2,ftres3,ftres1)); dtime_in_ns_res = G4RandGauss::shoot(dtime_in_ns,getTRES(ADC,dtres0,dtres2,dtres3,dtres1)); } @@ -478,6 +486,7 @@ map ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) // simulated at present. double fadc_time = convert_to_precision(ftime_in_ns); + int tdc = dtime_in_ns/da1; dgtz["hitn"] = hitn; diff --git a/source/hitprocess/clas12/ecal_hitprocess.h b/source/hitprocess/clas12/ecal_hitprocess.h index a43af38d..c986ec7e 100644 --- a/source/hitprocess/clas12/ecal_hitprocess.h +++ b/source/hitprocess/clas12/ecal_hitprocess.h @@ -29,6 +29,11 @@ class ecConstants //gain: pmt gain vector gain[nsect][nlayer]; + // tdc jitter parameters + double jitter_period; + int jitter_phase; + int jitter_cycles; + //timing: TDC calibration constants vector dtime[nsect][nlayer][9]; vector ftime[nsect][nlayer][7]; diff --git a/source/hitprocess/clas12/ftof_hitprocess.cc b/source/hitprocess/clas12/ftof_hitprocess.cc index 0be111cb..7f820d89 100644 --- a/source/hitprocess/clas12/ftof_hitprocess.cc +++ b/source/hitprocess/clas12/ftof_hitprocess.cc @@ -137,9 +137,32 @@ static ftofConstants initializeFTOFConstants(int runno, string digiVariation = " ftc.twlk[isec - 1][ilay - 1][4].push_back(data[row][7]); ftc.twlk[isec - 1][ilay - 1][5].push_back(data[row][8]); } - + + cout << "FTOF:Getting energy-dependent time_walk" << endl; + snprintf(ftc.database, sizeof(ftc.database), "/calibration/ftof/time_walk_exp:%d:%s%s", ftc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); + calib->GetCalib(data, ftc.database); + for (unsigned row = 0; row < data.size(); row++) { + isec = data[row][0]; + ilay = data[row][1]; + ftc.twlke[isec - 1][ilay - 1][0].push_back(data[row][3]); + ftc.twlke[isec - 1][ilay - 1][1].push_back(data[row][4]); + ftc.twlke[isec - 1][ilay - 1][2].push_back(data[row][5]); + ftc.twlke[isec - 1][ilay - 1][3].push_back(data[row][6]); + } + + cout << "FTOF:Getting position-dependent time_walk" << endl; + snprintf(ftc.database, sizeof(ftc.database), "/calibration/ftof/time_walk_pos:%d:%s%s", ftc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); + calib->GetCalib(data, ftc.database); + for (unsigned row = 0; row < data.size(); row++) { + isec = data[row][0]; + ilay = data[row][1]; + ftc.twlkp[isec - 1][ilay - 1][0].push_back(data[row][3]); + ftc.twlkp[isec - 1][ilay - 1][1].push_back(data[row][4]); + } + cout << "FTOF:Getting time_offset" << endl; - snprintf(ftc.database, sizeof(ftc.database), "/calibration/ftof/time_offsets:%d:%s%s", ftc.runNo, digiVariation.c_str(), timestamp.c_str()); data.clear(); calib->GetCalib(data, ftc.database); @@ -164,6 +187,14 @@ static ftofConstants initializeFTOFConstants(int runno, string digiVariation = " ftc.tdcconv[isec - 1][ilay - 1][1].push_back(data[row][4]); } + snprintf(ftc.database, sizeof(ftc.database), "/calibration/ftof/time_jitter:%d:%s%s", ftc.runNo, digiVariation.c_str(), timestamp.c_str()); + cout << "FTOF:Getting time_jitter" << endl; + data.clear(); + calib->GetCalib(data, ftc.database); + ftc.jitter_period = data[0][3]; + ftc.jitter_phase = data[0][4]; + ftc.jitter_cycles = data[0][5]; + cout << "FTOF:Getting adc_offsets" << endl; snprintf(ftc.database, sizeof(ftc.database), "/calibration/ftof/fadc_offset:%d:%s%s", ftc.runNo, digiVariation.c_str(), timestamp.c_str()); data.clear(); @@ -262,6 +293,7 @@ map ftof_HitProcess::integrateDgt(MHit* aHit, int hitn) { // TDC conversion factors double adcoffset = ftc.adcoffset[sector - 1][panel - 1][pmt][paddle - 1]; double tdcconv = ftc.tdcconv[sector - 1][panel - 1][pmt][paddle - 1]; + double tdc_jitter = ftc.jitter_period * ((0 + ftc.jitter_phase) % ftc.jitter_cycles); // assumes event timestamp is zero double time_in_ns = 0; if(aHit->isBackgroundHit == 1) { @@ -357,9 +389,17 @@ map ftof_HitProcess::integrateDgt(MHit* aHit, int hitn) { adc = ene * ftc.countsForMIP[sector - 1][panel - 1][pmt][paddle - 1] / ftc.dEMIP[panel - 1] / gain; double A = ftc.twlk[sector - 1][panel - 1][3 * pmt + 0][paddle - 1]; double B = ftc.twlk[sector - 1][panel - 1][3 * pmt + 1][paddle - 1]; + double E0 = ftc.twlke[sector - 1][panel - 1][0][paddle - 1]; + double E1 = ftc.twlke[sector - 1][panel - 1][1][paddle - 1]; + double E2 = ftc.twlke[sector - 1][panel - 1][2][paddle - 1]; + double E3 = ftc.twlke[sector - 1][panel - 1][3][paddle - 1]; + double P1 = ftc.twlkp[sector - 1][panel - 1][0][paddle - 1]; + double P2 = ftc.twlkp[sector - 1][panel - 1][1][paddle - 1]; //double C = ftc.twlk[sector-1][panel-1][2][paddle-1]; - double timeWalk = A / pow(adc, B); + double timeWalk = (E0-1) * A / pow(adc, B); + double timeWalkE = E0*E1*exp(E2*tInfos.eTot/MeV)+E3/(tInfos.eTot/MeV); + double timeWalkP = P1*pow(tInfos.lx/cm, 2)+P2*tInfos.lx/cm; // double timeWalkU = A / pow(adcu, B); // double tU = tInfos.time + d/ftc.veff[sector-1][panel-1][pmt][paddle-1]/cm + (1. - 2. * pmt)*ftc.toff_LR[sector-1][panel-1][paddle-1]/2. @@ -376,8 +416,9 @@ map ftof_HitProcess::integrateDgt(MHit* aHit, int hitn) { // cout << " FTOF Unsmeared Time after p2p subtraction: " << tU << endl; // tdcu = (tU + timeWalkU) / tdcconv; - time_in_ns = G4RandGauss::shoot(tU+ timeWalk, sqrt(2) * ftc.tres[sector - 1][panel - 1][paddle - 1]) ; - // tdc = time / tdcconv; + time_in_ns = G4RandGauss::shoot(tU+ timeWalk + timeWalkE - timeWalkP, sqrt(2) * ftc.tres[sector - 1][panel - 1][paddle - 1]) ; + + tdc = (time_in_ns + tdc_jitter) / tdcconv; } @@ -409,7 +450,6 @@ map ftof_HitProcess::integrateDgt(MHit* aHit, int hitn) { // cout << " > FTOF status: " << ftc.status[sector-1][panel-1][1][paddle-1] << " for sector " << sector << ", panel " << panel << ", paddle " << paddle << " right: " << adcr << endl; double fadc_time = convert_to_precision(time_in_ns - ftc.adcoffset[sector - 1][panel - 1][pmt][paddle - 1]); - tdc = time_in_ns / tdcconv; dgtz["sector"] = sector; dgtz["layer"] = panel; diff --git a/source/hitprocess/clas12/ftof_hitprocess.h b/source/hitprocess/clas12/ftof_hitprocess.h index 63a58666..096e40df 100644 --- a/source/hitprocess/clas12/ftof_hitprocess.h +++ b/source/hitprocess/clas12/ftof_hitprocess.h @@ -36,6 +36,11 @@ class ftofConstants // tdc_conc: tdc conversion factors vector tdcconv[6][3][2]; + // tdc jitter parameters + double jitter_period; + int jitter_phase; + int jitter_cycles; + // adc time offsets vector adcoffset[6][3][2]; @@ -48,8 +53,14 @@ class ftofConstants // countsForMIP: Desired ADC channel for MIP peak calibration vector countsForMIP[6][3][2]; - // twlk: Time walk correction, 3 constants each for L and R + // twlk: Time walk correction, 5 constants each for L and R vector twlk[6][3][6]; + + // twlk: Energy-dependent Time walk correction + vector twlke[6][3][4]; + + // twlk: Position-depemdent Time walk correction + vector twlkp[6][3][2]; // toff_LR and tof_P2P: time offsets for Left-Right and Paddle-to-Paddle vector toff_LR[6][3]; diff --git a/source/hitprocess/clas12/htcc_hitprocess.cc b/source/hitprocess/clas12/htcc_hitprocess.cc index 3a84d310..54a409ec 100644 --- a/source/hitprocess/clas12/htcc_hitprocess.cc +++ b/source/hitprocess/clas12/htcc_hitprocess.cc @@ -86,6 +86,14 @@ static htccConstants initializeHTCCConstants(int runno, string digiVariation = " htccc.tshift[isec-1][ilay-1].push_back(data[row][3]); } + + cout<<"HTCC:Getting ring_time"<GetCalib(data,htccc.database); + for(unsigned row = 0; row < data.size(); row++) { + htccc.rshift.push_back(data[row][3]); + } + cout<<"HTCC:Getting tdc_conv"<GetCalib(data,htccc.database); @@ -296,10 +304,10 @@ map htcc_HitProcess :: integrateDgt(MHit* aHit, int hitn) } double adc = htccc.gain[idsector-1][idhalf-1][idring-1] * G4RandGauss::shoot(ndetected*htccc.mc_gain[idsector-1][idhalf-1][idring-1], ndetected*htccc.mc_smear[idsector-1][idhalf-1][idring-1]); - time_in_ns = tInfos.time + htccc.tshift[idsector-1][idhalf-1][idring-1]; + time_in_ns = tInfos.time + htccc.tshift[idsector-1][idhalf-1][idring-1] + htccc.rshift[idring-1]; double fadc_time = convert_to_precision(time_in_ns); - int tdc = time_in_ns / tdcconv; + int tdc = time_in_ns / tdcconv; dgtz["hitn"] = hitn; diff --git a/source/hitprocess/clas12/htcc_hitprocess.h b/source/hitprocess/clas12/htcc_hitprocess.h index 20d02377..b25864c7 100644 --- a/source/hitprocess/clas12/htcc_hitprocess.h +++ b/source/hitprocess/clas12/htcc_hitprocess.h @@ -35,6 +35,9 @@ class htccConstants // veff: time shift vector tshift[6][2]; + // ring time shift + vector rshift; + // nephe to ADC conversion factors vector gain[6][2]; From 0b3a2e765537bf5374e596bade99f2ddd530cf9a Mon Sep 17 00:00:00 2001 From: raffaelladevita Date: Tue, 6 Jan 2026 19:55:09 -0500 Subject: [PATCH 04/12] fix bug in TDC jitter formula when ncycle is 0 --- source/hitprocess/clas12/band_hitprocess.cc | 3 ++- source/hitprocess/clas12/cnd_hitprocess.cc | 3 ++- source/hitprocess/clas12/ctof_hitprocess.cc | 3 ++- source/hitprocess/clas12/dc_hitprocess.cc | 7 ++++--- source/hitprocess/clas12/ecal_hitprocess.cc | 3 ++- source/hitprocess/clas12/ftof_hitprocess.cc | 3 ++- 6 files changed, 14 insertions(+), 8 deletions(-) diff --git a/source/hitprocess/clas12/band_hitprocess.cc b/source/hitprocess/clas12/band_hitprocess.cc index e34f1019..a3da4d36 100644 --- a/source/hitprocess/clas12/band_hitprocess.cc +++ b/source/hitprocess/clas12/band_hitprocess.cc @@ -305,7 +305,8 @@ map band_HitProcess :: integrateDgt(MHit* aHit, int hitn) double fadc_time = convert_to_precision(time_in_ns); - double tdc_jitter = bhc.jitter_period * ((0 + bhc.jitter_phase) % bhc.jitter_cycles); // assumes event timestamp is zero + double tdc_jitter = 0; + if(bhc.jitter_cycles != 0) tdc_jitter = bhc.jitter_period * ((0 + bhc.jitter_phase) % bhc.jitter_cycles); // assumes event timestamp is zero int TDC = (int) ( adcFactor*((side == 0 ? tL_tdc : tR_tdc)+tdc_jitter)/tdcconv ); dgtz["hitn"] = (int) hitn; diff --git a/source/hitprocess/clas12/cnd_hitprocess.cc b/source/hitprocess/clas12/cnd_hitprocess.cc index 5cb146f0..d3e7d05e 100644 --- a/source/hitprocess/clas12/cnd_hitprocess.cc +++ b/source/hitprocess/clas12/cnd_hitprocess.cc @@ -290,7 +290,8 @@ map cnd_HitProcess :: integrateDgt(MHit* aHit, int hitn) double time_in_ns = 0; - double tdc_jitter = cndc.jitter_period * ((0 + cndc.jitter_phase) % cndc.jitter_cycles); // assumes event timestamp is zero + double tdc_jitter = 0; + if(cndc.jitter_cycles != 0) tdc_jitter = cndc.jitter_period * ((0 + cndc.jitter_phase) % cndc.jitter_cycles); // assumes event timestamp is zero // LEFT PADDLE diff --git a/source/hitprocess/clas12/ctof_hitprocess.cc b/source/hitprocess/clas12/ctof_hitprocess.cc index 453a666b..24dd1097 100644 --- a/source/hitprocess/clas12/ctof_hitprocess.cc +++ b/source/hitprocess/clas12/ctof_hitprocess.cc @@ -389,7 +389,8 @@ map ctof_HitProcess::integrateDgt(MHit* aHit, int hitn) time_in_ns = G4RandGauss::shoot(tU, sqrt(2) * ctc.tres[paddle - 1]); // tdcu = tU / tdcconv; - double tdc_jitter = ctc.jitter_period * ((0 + ctc.jitter_phase) % ctc.jitter_cycles); // assumes event timestamp is zero + double tdc_jitter = 0; + if(ctc.jitter_cycles != 0) tdc_jitter = ctc.jitter_period * ((0 + ctc.jitter_phase) % ctc.jitter_cycles); // assumes event timestamp is zero tdc = (int) ((time_in_ns + tdc_jitter) / tdcconv); } diff --git a/source/hitprocess/clas12/dc_hitprocess.cc b/source/hitprocess/clas12/dc_hitprocess.cc index 1c09f6a3..40f4d55b 100644 --- a/source/hitprocess/clas12/dc_hitprocess.cc +++ b/source/hitprocess/clas12/dc_hitprocess.cc @@ -256,7 +256,8 @@ map dc_HitProcess :: integrateDgt(MHit* aHit, int hitn) double signal_t = 0; double hit_signal_t = 0; double prop_t = 0; - double tdc_jitter = dcc.jitter_period * ((0 + dcc.jitter_phase) % dcc.jitter_cycles); // assumes event timestamp is zero + double tdc_jitter = 0; + if(dcc.jitter_cycles != 0) tdc_jitter = dcc.jitter_period * ((0 + dcc.jitter_phase) % dcc.jitter_cycles); // assumes event timestamp is zero for(unsigned int s=0; s dc_HitProcess :: integrateDgt(MHit* aHit, int hitn) // new hit time // (w/o the drift time) // TODO: After coatjava real run numner -// hit_signal_t = stepTime[s]/ns; -// prop_t = tprop/ns; + hit_signal_t = stepTime[s]/ns; + prop_t = tprop/ns; if(Edep[s] >= dcc.dcThreshold*eV) { diff --git a/source/hitprocess/clas12/ecal_hitprocess.cc b/source/hitprocess/clas12/ecal_hitprocess.cc index 204b1f02..27002fc9 100644 --- a/source/hitprocess/clas12/ecal_hitprocess.cc +++ b/source/hitprocess/clas12/ecal_hitprocess.cc @@ -408,7 +408,8 @@ map ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) double dtres2 = ecc.dtres[sector-1][layer-1][2][0]; double dtres3 = ecc.dtres[sector-1][layer-1][3][0]; - double tdc_jitter = ecc.jitter_period * ((0 + ecc.jitter_phase) % ecc.jitter_cycles); // assumes event timestamp is zero + double tdc_jitter = 0; + if(ecc.jitter_cycles != 0) tdc_jitter = ecc.jitter_period * ((0 + ecc.jitter_phase) % ecc.jitter_cycles); // assumes event timestamp is zero for(unsigned int s=0; s ftof_HitProcess::integrateDgt(MHit* aHit, int hitn) { // TDC conversion factors double adcoffset = ftc.adcoffset[sector - 1][panel - 1][pmt][paddle - 1]; double tdcconv = ftc.tdcconv[sector - 1][panel - 1][pmt][paddle - 1]; - double tdc_jitter = ftc.jitter_period * ((0 + ftc.jitter_phase) % ftc.jitter_cycles); // assumes event timestamp is zero + double tdc_jitter = 0; + if(ftc.jitter_cycles != 0) tdc_jitter = ftc.jitter_period * ((0 + ftc.jitter_phase) % ftc.jitter_cycles); // assumes event timestamp is zero double time_in_ns = 0; if(aHit->isBackgroundHit == 1) { From ce3c3c8f5d5d3809c81bd4f66e2a7c675cbac739 Mon Sep 17 00:00:00 2001 From: raffaelladevita Date: Tue, 13 Jan 2026 12:21:36 -0500 Subject: [PATCH 05/12] switching DC hit process to v2 constants tables --- source/hitprocess/clas12/dc_hitprocess.cc | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/source/hitprocess/clas12/dc_hitprocess.cc b/source/hitprocess/clas12/dc_hitprocess.cc index 40f4d55b..10fd2ea1 100644 --- a/source/hitprocess/clas12/dc_hitprocess.cc +++ b/source/hitprocess/clas12/dc_hitprocess.cc @@ -78,7 +78,7 @@ static dcConstants initializeDCConstants(int runno, string digiVariation = "defa //******************************************** //reading reference and current-run pressure: - snprintf(dcc.database, sizeof(dcc.database), "/calibration/dc/time_to_distance/ref_pressure:%d:%s%s", dcc.runNo, digiVariation.c_str(), timestamp.c_str()); + snprintf(dcc.database, sizeof(dcc.database), "/calibration/dc/v2/ref_pressure:%d:%s%s", dcc.runNo, digiVariation.c_str(), timestamp.c_str()); data.clear(); calib->GetCalib(data, dcc.database); double ref_pressure = data[0][3]; @@ -89,7 +89,7 @@ static dcConstants initializeDCConstants(int runno, string digiVariation = "defa double dpressure = current_pressure - ref_pressure; //******************************************** //calculating distance to time: - snprintf(dcc.database, sizeof(dcc.database), "/calibration/dc/time_to_distance/t2d_pressure:%d:%s%s", dcc.runNo, digiVariation.c_str(), timestamp.c_str()); + snprintf(dcc.database, sizeof(dcc.database), "/calibration/dc/v2/t2d_pressure:%d:%s%s", dcc.runNo, digiVariation.c_str(), timestamp.c_str()); data.clear(); calib->GetCalib(data, dcc.database); @@ -126,7 +126,7 @@ static dcConstants initializeDCConstants(int runno, string digiVariation = "defa dcc.vprop = 29.97924580*0.7*cm/ns; // hardcoded in reconstruction too // T0 corrections: a delay to be introduced (plus sign) to the TDC timing - snprintf(dcc.database, sizeof(dcc.database), "/calibration/dc/time_corrections/T0Corrections:%d:%s%s", dcc.runNo, digiVariation.c_str(), timestamp.c_str()); + snprintf(dcc.database, sizeof(dcc.database), "/calibration/dc/v2/t0:%d:%s%s", dcc.runNo, digiVariation.c_str(), timestamp.c_str()); data.clear(); calib->GetCalib(data, dcc.database); From b984674b6ea35c30722872a51cb346b5f5f38fd6 Mon Sep 17 00:00:00 2001 From: raffaelladevita Date: Tue, 13 Jan 2026 14:51:40 -0500 Subject: [PATCH 06/12] fixed sign of old FTOF time walk correction --- source/hitprocess/clas12/ftof_hitprocess.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/source/hitprocess/clas12/ftof_hitprocess.cc b/source/hitprocess/clas12/ftof_hitprocess.cc index e5902791..0d48fe59 100644 --- a/source/hitprocess/clas12/ftof_hitprocess.cc +++ b/source/hitprocess/clas12/ftof_hitprocess.cc @@ -398,7 +398,7 @@ map ftof_HitProcess::integrateDgt(MHit* aHit, int hitn) { double P2 = ftc.twlkp[sector - 1][panel - 1][1][paddle - 1]; //double C = ftc.twlk[sector-1][panel-1][2][paddle-1]; - double timeWalk = (E0-1) * A / pow(adc, B); + double timeWalk = (1-E0) * A / pow(adc, B); double timeWalkE = E0*E1*exp(E2*tInfos.eTot/MeV)+E3/(tInfos.eTot/MeV); double timeWalkP = P1*pow(tInfos.lx/cm, 2)+P2*tInfos.lx/cm; // double timeWalkU = A / pow(adcu, B); From 4103b1724916e1d6fb167909b43857112f92fa5d Mon Sep 17 00:00:00 2001 From: raffaelladevita Date: Tue, 7 Apr 2026 19:06:25 -0400 Subject: [PATCH 07/12] Revert "reverted changes to pre real run number for ecal" This reverts commit 104c6a60f532376fd7ffde5c536915d10e9e4ddf. --- source/hitprocess/clas12/ecal_hitprocess.cc | 296 ++++++-- source/hitprocess/clas12/ecal_hitprocess.h | 89 ++- .../hitprocess/clas12/ecal_hitprocess_rr.cc | 695 ------------------ source/hitprocess/clas12/ecal_hitprocess_rr.h | 131 ---- 4 files changed, 286 insertions(+), 925 deletions(-) delete mode 100644 source/hitprocess/clas12/ecal_hitprocess_rr.cc delete mode 100644 source/hitprocess/clas12/ecal_hitprocess_rr.h diff --git a/source/hitprocess/clas12/ecal_hitprocess.cc b/source/hitprocess/clas12/ecal_hitprocess.cc index 1c0773a1..858d613f 100644 --- a/source/hitprocess/clas12/ecal_hitprocess.cc +++ b/source/hitprocess/clas12/ecal_hitprocess.cc @@ -33,7 +33,7 @@ static ecConstants initializeECConstants(int runno, string digiVariation = "defa ecc.connection = "mysql://clas12reader@clasdb.jlab.org/clas12"; } - ecc.ADC_GeV_to_evio = 1./10000.; // MIP based calibration is nominally 10 channels/MeV + ecc.ADC_GeV_to_evio = 1./10000.; // MIP calibration: 10(15) ch/MeV for ECAL(PCAL) 50 for ECAL Sector 5 ecc.pmtQE = 0.27 ; ecc.pmtDynodeGain = 4.0 ; @@ -43,9 +43,10 @@ static ecConstants initializeECConstants(int runno, string digiVariation = "defa ecc.pmtFactor = sqrt(1 + 1/(ecc.pmtDynodeGain-1)); - // The callibration data will be filled in this vector data + // The calibration data will be filled in this vector data vector > data; unique_ptr calib(CalibrationGenerator::CreateCalibration(ecc.connection)); + // ======== Initialization of EC gains =========== snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/gain:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); @@ -58,7 +59,7 @@ static ecConstants initializeECConstants(int runno, string digiVariation = "defa } // ========= Initializations of attenuation lengths ======== - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/attenuation:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/atten:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); data.clear(); calib->GetCalib(data,ecc.database); for(unsigned row = 0; row < data.size(); row++) @@ -67,36 +68,142 @@ static ecConstants initializeECConstants(int runno, string digiVariation = "defa ecc.attlen[isec-1][ilay-1][0].push_back(data[row][3]); ecc.attlen[isec-1][ilay-1][1].push_back(data[row][5]); ecc.attlen[isec-1][ilay-1][2].push_back(data[row][7]); + ecc.attlen[isec-1][ilay-1][3].push_back(data[row][9]); + ecc.attlen[isec-1][ilay-1][4].push_back(data[row][11]); } // ========== Initialization of timings =========== - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/timing:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/ftime:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); data.clear(); calib->GetCalib(data,ecc.database); + + for(unsigned row = 0; row < data.size(); row++) + { + isec = data[row][0]; ilay = data[row][1]; + ecc.ftime[isec-1][ilay-1][0].push_back(data[row][3]); + ecc.ftime[isec-1][ilay-1][1].push_back(data[row][4]); + ecc.ftime[isec-1][ilay-1][2].push_back(data[row][5]); + ecc.ftime[isec-1][ilay-1][3].push_back(data[row][6]); + ecc.ftime[isec-1][ilay-1][4].push_back(data[row][7]); + ecc.ftime[isec-1][ilay-1][5].push_back(data[row][8]); + ecc.ftime[isec-1][ilay-1][6].push_back(data[row][9]); + } + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/dtime:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); calib->GetCalib(data,ecc.database); + for(unsigned row = 0; row < data.size(); row++) { isec = data[row][0]; ilay = data[row][1]; - ecc.timing[isec-1][ilay-1][0].push_back(data[row][3]); - ecc.timing[isec-1][ilay-1][1].push_back(data[row][4]); - ecc.timing[isec-1][ilay-1][2].push_back(data[row][5]); - ecc.timing[isec-1][ilay-1][3].push_back(data[row][6]); - ecc.timing[isec-1][ilay-1][4].push_back(data[row][7]); + ecc.dtime[isec-1][ilay-1][0].push_back(data[row][3]); + ecc.dtime[isec-1][ilay-1][1].push_back(data[row][4]); + ecc.dtime[isec-1][ilay-1][2].push_back(data[row][5]); + ecc.dtime[isec-1][ilay-1][3].push_back(data[row][6]); + ecc.dtime[isec-1][ilay-1][4].push_back(data[row][7]); + ecc.dtime[isec-1][ilay-1][5].push_back(data[row][8]); + ecc.dtime[isec-1][ilay-1][6].push_back(data[row][9]); + ecc.dtime[isec-1][ilay-1][7].push_back(data[row][10]); + ecc.dtime[isec-1][ilay-1][8].push_back(data[row][11]); } - // ========== Initialization of timing offset =========== + //========== Initialization of timing offsets and global constants =========== snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/tdc_global_offset:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); data.clear(); calib->GetCalib(data, ecc.database); ecc.tdc_global_offset = data[0][3]; + + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/fadc_global_offset:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); calib->GetCalib(data, ecc.database); + ecc.fadc_global_offset = data[0][3]; + + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/global_time_walk:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); calib->GetCalib(data, ecc.database); + for(unsigned row = 0; row < data.size(); row++) + { + isec = data[row][0]; ilay = data[row][1]; + ecc.global_time_walk[isec-1][ilay-1].push_back(data[row][3]); + } + + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/fadc_offset:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); calib->GetCalib(data, ecc.database); + for(unsigned row = 0; row < data.size(); row++) + { + isec = data[row][0]; ilay = data[row][1]; + ecc.fadc_offset[isec-1][ilay-1].push_back(data[row][3]); + } + + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/tmf_offset:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); calib->GetCalib(data, ecc.database); + for(unsigned row = 0; row < data.size(); row++) + { + isec = data[row][0]; ilay = data[row][1]; + ecc.tmf_offset[isec-1][ilay-1].push_back(data[row][3]); + } // ======== Initialization of EC effective velocities =========== - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/effective_velocity:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/fveff:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); data.clear(); calib->GetCalib(data,ecc.database); + + for(unsigned row = 0; row < data.size(); row++) + { + isec = data[row][0]; ilay = data[row][1]; + ecc.fveff[isec-1][ilay-1].push_back(data[row][3]); + } + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/dveff:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); calib->GetCalib(data,ecc.database); + + for(unsigned row = 0; row < data.size(); row++) + { + isec = data[row][0]; ilay = data[row][1]; + ecc.dveff[isec-1][ilay-1].push_back(data[row][3]); + } + + // ======== Initialization of EC DSC/TDC efficiency vs FADC parameters =========== + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/deff:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); calib->GetCalib(data,ecc.database); + for(unsigned row = 0; row < data.size(); row++) { isec = data[row][0]; ilay = data[row][1]; - ecc.veff[isec-1][ilay-1].push_back(data[row][3]); + ecc.deff[isec-1][ilay-1][0].push_back(data[row][3]); + ecc.deff[isec-1][ilay-1][1].push_back(data[row][4]); + ecc.deff[isec-1][ilay-1][2].push_back(data[row][5]); } + + // ======== Initialization of FADC thresholds =========== + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/fthr:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); calib->GetCalib(data,ecc.database); + + for(unsigned row = 0; row < data.size(); row++) + { + isec = data[row][0]; ilay = data[row][1]; + ecc.fthr[isec-1][ilay-1].push_back(data[row][3]); + } + + // ======== Initialization of timing resolution vs FADC parameters =========== + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/ftres:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); calib->GetCalib(data,ecc.database); + + for(unsigned row = 0; row < data.size(); row++) + { + isec = data[row][0]; ilay = data[row][1]; + ecc.ftres[isec-1][ilay-1][0].push_back(data[row][3]); + ecc.ftres[isec-1][ilay-1][1].push_back(data[row][4]); + ecc.ftres[isec-1][ilay-1][2].push_back(data[row][5]); + ecc.ftres[isec-1][ilay-1][3].push_back(data[row][6]); + } + + snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/dtres:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); + data.clear(); calib->GetCalib(data,ecc.database); + + for(unsigned row = 0; row < data.size(); row++) + { + isec = data[row][0]; ilay = data[row][1]; + ecc.dtres[isec-1][ilay-1][0].push_back(data[row][3]); + ecc.dtres[isec-1][ilay-1][1].push_back(data[row][4]); + ecc.dtres[isec-1][ilay-1][2].push_back(data[row][5]); + ecc.dtres[isec-1][ilay-1][3].push_back(data[row][6]); + } // ======== Initialization of EC status =========== if(accountForHardwareStatus) { @@ -159,7 +266,6 @@ static ecConstants initializeECConstants(int runno, string digiVariation = "defa return ecc; } - // Process the ID and hit for the EC using EC scintillator slab geometry instead of individual strips. map ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) { @@ -176,9 +282,6 @@ map ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) bool isPCAL = layer < 4 ; - double time_in_ns = 0; - - // layer = 1, 2 stays the same // subtract 3 from ec inner // subtract 6 from ec outer @@ -195,11 +298,10 @@ map ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) // Different for EC and PCAL double pmtPEYld = 3.5 ; if (isPCAL) { - // pcal pmtPEYld = 11.5 ; } - double a1 = ecc.timing[sector-1][layer-1][1][strip-1]; // tdc conversion + double a1 = ecc.ftime[sector-1][layer-1][1][strip-1]; // tdc conversion if(aHit->isBackgroundHit == 1) { @@ -233,9 +335,17 @@ map ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) // Get scintillator volume x dimension (mm) double pDx2 = aHit->GetDetector().dimensions[5]; ///< G4Trap Semilength. + // initialize ADC and TDC + double ADC = 0; + double ftime_in_ns = 0; + double dtime_in_ns = 0; + double ftime_in_ns_res = 0; + double dtime_in_ns_res = 0; + // Get Total Energy deposited double Etota = 0; - double Ttota = 0; + double FTtota = 0; + double DTtota = 0; double latt = 0; double att; @@ -243,37 +353,79 @@ map ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) double A = ecc.attlen[sector-1][layer-1][0][strip-1]; double B = ecc.attlen[sector-1][layer-1][1][strip-1]*10.; double C = ecc.attlen[sector-1][layer-1][2][strip-1]; + double D = ecc.attlen[sector-1][layer-1][3][strip-1]; + double E = ecc.attlen[sector-1][layer-1][4][strip-1]*10.; + double G = ecc.gain[sector-1][layer-1][strip-1]; - double a0 = ecc.timing[sector-1][layer-1][0][strip-1]; - double a2 = ecc.timing[sector-1][layer-1][2][strip-1]; - double veff = ecc.veff[sector-1][layer-1][strip-1]*10; + double tmf = ecc.tmf_offset[sector-1][layer-1][strip-1]; + double fo = ecc.fadc_offset[sector-1][layer-1][0]; + double gtw = ecc.global_time_walk[sector-1][layer-1][0]; + + double FTOFFSET = ecc.fadc_global_offset; + double tgo = ecc.tdc_global_offset; + + double fa0 = ecc.ftime[sector-1][layer-1][0][strip-1]; + double fa2 = ecc.ftime[sector-1][layer-1][2][strip-1]; + double fa3 = ecc.ftime[sector-1][layer-1][3][strip-1]; + double fa4 = ecc.ftime[sector-1][layer-1][4][strip-1]; + double fa5 = ecc.ftime[sector-1][layer-1][5][strip-1]; + double fa6 = ecc.ftime[sector-1][layer-1][6][strip-1]; + + double da0 = ecc.dtime[sector-1][layer-1][0][strip-1]; + double da1 = ecc.dtime[sector-1][layer-1][1][strip-1]; + double da2 = ecc.dtime[sector-1][layer-1][2][strip-1]; + double da3 = ecc.dtime[sector-1][layer-1][3][strip-1]; + double da4 = ecc.dtime[sector-1][layer-1][4][strip-1]; + double da5 = ecc.dtime[sector-1][layer-1][5][strip-1]; + double da6 = ecc.dtime[sector-1][layer-1][6][strip-1]; + double da7 = ecc.dtime[sector-1][layer-1][7][strip-1]; + double da8 = ecc.dtime[sector-1][layer-1][8][strip-1]; + ; + double fveff = ecc.fveff[sector-1][layer-1][strip-1]*10; + double dveff = ecc.dveff[sector-1][layer-1][strip-1]*10; + + double fthr = ecc.fthr[sector-1][layer-1][strip-1]; + + double def0 = ecc.deff[sector-1][layer-1][0][strip-1]; + double def1 = ecc.deff[sector-1][layer-1][1][strip-1]; + double def2 = ecc.deff[sector-1][layer-1][2][strip-1]; + + + double ftres0 = ecc.ftres[sector-1][layer-1][0][0]; + double ftres1 = ecc.ftres[sector-1][layer-1][1][0]; + double ftres2 = ecc.ftres[sector-1][layer-1][2][0]; + double ftres3 = ecc.ftres[sector-1][layer-1][3][0]; + + double dtres0 = ecc.dtres[sector-1][layer-1][0][0]; + double dtres1 = ecc.dtres[sector-1][layer-1][1][0]; + double dtres2 = ecc.dtres[sector-1][layer-1][2][0]; + double dtres3 = ecc.dtres[sector-1][layer-1][3][0]; + for(unsigned int s=0; s0) { - double xlocal = Lpos[s].x(); - if(view==1) latt = pDx2 + xlocal; - if(view==2) latt = pDx2 + xlocal; - if(view==3) { - if(layer > 3) { - // for ec, it's a minus sign - latt = pDx2-xlocal; - } else { - // for pcal, it's a plus sign - latt = pDx2+xlocal; - } - } - att = A*exp(-latt/B)+C; - Etota = Etota + Edep[s]*att; - Ttota = Ttota + latt/veff; - } else { - Etota = Etota + Edep[s]; + double xlocal = Lpos[s].x(); + if(view==1) latt = pDx2 + xlocal; + if(view==2) latt = pDx2 + xlocal; + if(view==3) { + if(layer > 3) { + // for ec, it's a minus sign + latt = pDx2-xlocal; + } else { + // for pcal, it's a plus sign + latt = pDx2+xlocal; + } } + att = A*(exp(-latt/B) + D*exp(-latt/E)) + C; //pass2 parameterization + Etota = Etota + Edep[s]*att; //reported in MeV + FTtota = FTtota + latt/fveff; //FADC based timing + DTtota = DTtota + latt/dveff; //DSC/TDC based timing } - - // initialize ADC and TDC - double ADC = 0; - - // simulate the adc value. + + //Used if ecc.outputRAW > 0 (no digitization, effieiency or resolution smearing) + double ADC_raw = Etota/1000/ecc.ADC_GeV_to_evio/G; + double FTIME_raw = tInfos.time+FTtota/tInfos.nsteps; + double DTIME_raw = tInfos.time+DTtota/tInfos.nsteps; + if (Etota > 0) { double EC_npe = G4Poisson(Etota*pmtPEYld); //number of photoelectrons if (EC_npe>0) { @@ -281,24 +433,33 @@ map ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) double EC_GeV = G4RandGauss::shoot(EC_npe,sigma)/1000./ecc.ADC_GeV_to_evio/G/pmtPEYld; if (EC_GeV>0) { ADC = EC_GeV; - time_in_ns = (tInfos.time+Ttota/tInfos.nsteps) + a0 + a2/sqrt(ADC) + ecc.tdc_global_offset; + double radc = sqrt(ADC); + double ftim = (fa4==0||fa6==0) ? 0 : fa0 + fa2+exp(-(radc-fa3)/fa4)+1-exp( (radc-fa5)/fa6); + double dtim = (da4==0||da6==0||da7==0) ? 0 : da0 + gtw/radc +da2+exp(-(radc-da3)/da4)+1-exp(-(da5-radc)/da6)-exp(-(radc-da3*0.95)/da7)*pow(radc,da8); + ftime_in_ns = FTIME_raw + ftim + tgo - FTOFFSET - tmf - fo; + dtime_in_ns = DTIME_raw + dtim + tgo; + ftime_in_ns_res = G4RandGauss::shoot(ftime_in_ns,getTRES(ADC,ftres0,ftres2,ftres3,ftres1)); + dtime_in_ns_res = G4RandGauss::shoot(dtime_in_ns,getTRES(ADC,dtres0,dtres2,dtres3,dtres1)); } } } - + + dtime_in_ns = ecc.outputRAW>0 ? DTIME_raw : dtime_in_ns_res; + ftime_in_ns = ecc.outputRAW>0 ? FTIME_raw : ftime_in_ns_res; + // Status flags if(accountForHardwareStatus) { switch (ecc.status[sector-1][layer-1][strip-1]) { case 0: break; case 1: - ADC = 0; + ADC = ADC_raw = ftime_in_ns = 0; break; case 2: - time_in_ns = 0; + dtime_in_ns = 0; break; case 3: - ADC = time_in_ns = 0; + ADC = ADC_raw = ftime_in_ns = dtime_in_ns = 0; break; case 5: @@ -308,27 +469,28 @@ map ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) cout << " > Unknown EC status: " << ecc.status[sector-1][layer-1][strip-1] << " for sector " << sector << ", layer " << layer << ", strip " << strip << endl; } } - + + if (ecc.outputRAW==0 && def0>0 && dtime_in_ns > 0 && G4UniformRand() > 1/pow(1+exp(-def0*(ADC/10-def1)),def2)) dtime_in_ns = 0; // DSC/TDC threshold + if (ecc.outputRAW==0 && ADC/10 < fthr) rejectHitConditions = true; // FADC threshold + // EVIO banks record time with offset determined by position of data in capture window. On forward carriage this is currently // around 7.9 us. This offset is omitted in the simulation. Also EVIO TDC time is relative to the trigger time, which is not // simulated at present. - double fadc_time = convert_to_precision(time_in_ns); - int tdc = time_in_ns/a1; + double fadc_time = convert_to_precision(ftime_in_ns); + int tdc = dtime_in_ns/da1; dgtz["hitn"] = hitn; dgtz["sector"] = sector; dgtz["layer"] = layer; dgtz["component"] = strip; dgtz["ADC_order"] = 0; - dgtz["ADC_ADC"] = ADC; + dgtz["ADC_ADC"] = ecc.outputRAW==1 ? ADC_raw : ADC; dgtz["ADC_time"] = fadc_time; dgtz["ADC_ped"] = 0; dgtz["TDC_order"] = 2; dgtz["TDC_TDC"] = tdc; - // cout << "sector = " << sector << " layer = " << view << " strip = " << strip << " ADC = " << ADC << " TDC = " << TDC << endl; - // define conditions to reject hit if(rejectHitConditions) { writeHit = false; @@ -445,11 +607,14 @@ map< int, vector > ecal_HitProcess :: chargeTime(MHit* aHit, int hitn) vector Edep = aHit->GetEdep(); vector time = aHit->GetTime(); - double A = ecc.attlen[sector-1][layer-1][0][strip-1]; - double B = ecc.attlen[sector-1][layer-1][1][strip-1]*10.; - double C = ecc.attlen[sector-1][layer-1][2][strip-1]; - double G = ecc.gain[sector-1][layer-1][strip-1]; - double veff = ecc.veff[sector-1][layer-1][strip-1]*10; + double A = ecc.attlen[sector-1][layer-1][0][strip-1]; + double B = ecc.attlen[sector-1][layer-1][1][strip-1]*10.; + double C = ecc.attlen[sector-1][layer-1][2][strip-1]; + double D = ecc.attlen[sector-1][layer-1][3][strip-1]; + double E = ecc.attlen[sector-1][layer-1][4][strip-1]*10.; + + double G = ecc.gain[sector-1][layer-1][strip-1]; + double fveff = ecc.fveff[sector-1][layer-1][strip-1]*10; for(unsigned int s=0; s0) { @@ -463,14 +628,15 @@ map< int, vector > ecal_HitProcess :: chargeTime(MHit* aHit, int hitn) // for ecal, it's a minus sign latt = pDx2-xlocal; } else { - // for ecal, it's a plus sign + // for pcal, it's a plus sign latt = pDx2+xlocal; } } - double att = A*exp(-latt/B)+C; + + double att = A*(exp(-latt/B) + D*exp(-latt/E)) + C; //pass2 parameterization double stepE = Edep[s]*att; - double stepTime = time[s] + latt/veff; + double stepTime = time[s] + latt/fveff; // cout<<"time[s] = "< -#include -#include -using namespace ccdb; - -// gemc headers -#include "ecal_hitprocess.h" - -static ecConstants initializeECConstants(int runno, string digiVariation = "default", string digiSnapshotTime = "no", bool accountForHardwareStatus = false) -{ - ecConstants ecc; - - // 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 ecc; - string timestamp = ""; - if(digiSnapshotTime != "no") { - timestamp = ":"+digiSnapshotTime; - } - - int isec,ilay; - - // database - ecc.runNo = runno; - - if(getenv ("CCDB_CONNECTION") != nullptr) { - ecc.connection = (string) getenv("CCDB_CONNECTION"); - } else { - ecc.connection = "mysql://clas12reader@clasdb.jlab.org/clas12"; - } - - ecc.ADC_GeV_to_evio = 1./10000.; // MIP calibration: 10(15) ch/MeV for ECAL(PCAL) 50 for ECAL Sector 5 - ecc.pmtQE = 0.27 ; - ecc.pmtDynodeGain = 4.0 ; - - // Fluctuations in PMT gain distributed using Gaussian with - // sigma = sqrt(npe)/SNR where 1/SNR = sqrt[(1 + 1/(ecc.pmtDynodeGain-1)) npe=number of photoelectrons - // Adapted from G-112 (pg. 174) of RCA PMT Handbook. - - ecc.pmtFactor = sqrt(1 + 1/(ecc.pmtDynodeGain-1)); - - // The calibration data will be filled in this vector data - vector > data; - unique_ptr calib(CalibrationGenerator::CreateCalibration(ecc.connection)); - - - // ======== Initialization of EC gains =========== - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/gain:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data,ecc.database); - - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.gain[isec-1][ilay-1].push_back(data[row][3]); - } - - // ========= Initializations of attenuation lengths ======== - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/atten:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data,ecc.database); - - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.attlen[isec-1][ilay-1][0].push_back(data[row][3]); - ecc.attlen[isec-1][ilay-1][1].push_back(data[row][5]); - ecc.attlen[isec-1][ilay-1][2].push_back(data[row][7]); - ecc.attlen[isec-1][ilay-1][3].push_back(data[row][9]); - ecc.attlen[isec-1][ilay-1][4].push_back(data[row][11]); - } - - // ========== Initialization of timings =========== - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/ftime:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data,ecc.database); - - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.ftime[isec-1][ilay-1][0].push_back(data[row][3]); - ecc.ftime[isec-1][ilay-1][1].push_back(data[row][4]); - ecc.ftime[isec-1][ilay-1][2].push_back(data[row][5]); - ecc.ftime[isec-1][ilay-1][3].push_back(data[row][6]); - ecc.ftime[isec-1][ilay-1][4].push_back(data[row][7]); - ecc.ftime[isec-1][ilay-1][5].push_back(data[row][8]); - ecc.ftime[isec-1][ilay-1][6].push_back(data[row][9]); - } - - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/dtime:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data,ecc.database); - - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.dtime[isec-1][ilay-1][0].push_back(data[row][3]); - ecc.dtime[isec-1][ilay-1][1].push_back(data[row][4]); - ecc.dtime[isec-1][ilay-1][2].push_back(data[row][5]); - ecc.dtime[isec-1][ilay-1][3].push_back(data[row][6]); - ecc.dtime[isec-1][ilay-1][4].push_back(data[row][7]); - ecc.dtime[isec-1][ilay-1][5].push_back(data[row][8]); - ecc.dtime[isec-1][ilay-1][6].push_back(data[row][9]); - ecc.dtime[isec-1][ilay-1][7].push_back(data[row][10]); - ecc.dtime[isec-1][ilay-1][8].push_back(data[row][11]); - } - - //========== Initialization of timing offsets and global constants =========== - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/tdc_global_offset:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data, ecc.database); - ecc.tdc_global_offset = data[0][3]; - - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/fadc_global_offset:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data, ecc.database); - ecc.fadc_global_offset = data[0][3]; - - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/global_time_walk:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data, ecc.database); - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.global_time_walk[isec-1][ilay-1].push_back(data[row][3]); - } - - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/fadc_offset:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data, ecc.database); - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.fadc_offset[isec-1][ilay-1].push_back(data[row][3]); - } - - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/tmf_offset:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data, ecc.database); - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.tmf_offset[isec-1][ilay-1].push_back(data[row][3]); - } - - // ======== Initialization of EC effective velocities =========== - - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/fveff:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data,ecc.database); - - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.fveff[isec-1][ilay-1].push_back(data[row][3]); - } - - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/dveff:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data,ecc.database); - - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.dveff[isec-1][ilay-1].push_back(data[row][3]); - } - - // ======== Initialization of EC DSC/TDC efficiency vs FADC parameters =========== - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/deff:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data,ecc.database); - - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.deff[isec-1][ilay-1][0].push_back(data[row][3]); - ecc.deff[isec-1][ilay-1][1].push_back(data[row][4]); - ecc.deff[isec-1][ilay-1][2].push_back(data[row][5]); - } - - // ======== Initialization of FADC thresholds =========== - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/fthr:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data,ecc.database); - - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.fthr[isec-1][ilay-1].push_back(data[row][3]); - } - - // ======== Initialization of timing resolution vs FADC parameters =========== - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/ftres:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data,ecc.database); - - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.ftres[isec-1][ilay-1][0].push_back(data[row][3]); - ecc.ftres[isec-1][ilay-1][1].push_back(data[row][4]); - ecc.ftres[isec-1][ilay-1][2].push_back(data[row][5]); - ecc.ftres[isec-1][ilay-1][3].push_back(data[row][6]); - } - - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/dtres:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); calib->GetCalib(data,ecc.database); - - for(unsigned row = 0; row < data.size(); row++) - { - isec = data[row][0]; ilay = data[row][1]; - ecc.dtres[isec-1][ilay-1][0].push_back(data[row][3]); - ecc.dtres[isec-1][ilay-1][1].push_back(data[row][4]); - ecc.dtres[isec-1][ilay-1][2].push_back(data[row][5]); - ecc.dtres[isec-1][ilay-1][3].push_back(data[row][6]); - } - - // ======== Initialization of EC status =========== - if(accountForHardwareStatus) { - snprintf(ecc.database, sizeof(ecc.database), "/calibration/ec/status:%d:%s%s", ecc.runNo, digiVariation.c_str(), timestamp.c_str()); - data.clear(); - calib->GetCalib(data, ecc.database); - for (unsigned row = 0; row < data.size(); row++) { - isec = data[row][0]; ilay = data[row][1]; - ecc.status[isec-1][ilay-1].push_back(data[row][3]); - } - } - - // =========== Initialization of FADC250 related informations, pedestals, nsa, nsb ====================== - // FOR now we will initialize pedestals and sigmas to a random value, - // in the future they should come from CCDB - const double const_ped_value = 101; - const double const_ped_sigm_value = 2; - // commands below fill all the elements of ecc.pedestal and ecc.pedestal_sigm with their values (const_ped_value, and const_ped_sigm_value respectively) - std::fill(&ecc.pedestal[0][0][0], &ecc.pedestal[0][0][0] + sizeof(ecc.pedestal)/sizeof(ecc.pedestal[0][0][0]), const_ped_value); - std::fill(&ecc.pedestal_sigm[0][0][0], &ecc.pedestal_sigm[0][0][0] + sizeof(ecc.pedestal_sigm)/sizeof(ecc.pedestal_sigm[0][0][0]), const_ped_sigm_value); - - // setting voltage signal parameters - ecc.vpar[0] = 0.; // delay, ns - ecc.vpar[1] = 2.8; // rise time, ns - ecc.vpar[2] = 20; // fall time, ns - ecc.vpar[3] = 1; // amplifier - - - // loading translation table - ecc.TT = TranslationTable("ecTT"); - - // loads translation table from CLAS12 Database: - // Translation table for EC (ECAL+PCAL). - // Crate sector assignments: ECAL/FADC=1,7,13,19,25,31 ECAL/TDC=2,8,14,20,26,32 - // PCAL/FADC=3,9,15,21,27,33 PCAL/TDC=4,10,16,22,28,34. - // ORDER: 0=FADC 2=TDC. - - string database = "/daq/tt/ec:1"; - - data.clear(); calib->GetCalib(data, database); - cout << " > " << ecc.TT.getName() << " TT Data loaded from CCDB with " << data.size() << " columns." << endl; - - // filling translation table - for(unsigned row = 0; row < data.size(); row++) - { - int crate = data[row][0]; - int slot = data[row][1]; - int channel = data[row][2]; - - int sector = data[row][3]; - int layer = data[row][4]; - int pmt = data[row][5]; - int order = data[row][6]; - - // order is important as we could have duplicate entries w/o it - ecc.TT.addHardwareItem({sector, layer, pmt, order}, Hardware(crate, slot, channel)); - } - cout << " > Data loaded in translation table " << ecc.TT.getName() << endl; - - return ecc; -} - -// Process the ID and hit for the EC using EC scintillator slab geometry instead of individual strips. -map ecal_HitProcess :: integrateDgt(MHit* aHit, int hitn) -{ - map dgtz; - vector identity = aHit->GetId(); - rejectHitConditions = false; - writeHit = true; - - int sector = identity[0].id; - int layer = identity[1].id; // layer=1-3 (PCAL) 4-9 (ECAL). Layer = view for pcal, ecinner, ecouter - int strip = identity[2].id; - // pcal - int view = layer; - - bool isPCAL = layer < 4 ; - - // layer = 1, 2 stays the same - // subtract 3 from ec inner - // subtract 6 from ec outer - - if (layer > 3 && layer <7) { - // ec inner (stack 1) - view = layer - 3; - } else if (layer >= 7) { - // ec outer (stack 2) - view = layer - 6; - } - - // Number of p.e. divided by the energy deposited in MeV. See EC NIM paper table 1. - // Different for EC and PCAL - double pmtPEYld = 3.5 ; - if (isPCAL) { - pmtPEYld = 11.5 ; - } - - double a1 = ecc.ftime[sector-1][layer-1][1][strip-1]; // tdc conversion - - if(aHit->isBackgroundHit == 1) { - - // background hit has all the energy in the first step. Time is also first step - double totEdep = aHit->GetEdep()[0]; - double stepTime = aHit->GetTime()[0]; - double adc = totEdep / ecc.ADC_GeV_to_evio ; // no gain as that comes from data already - int tdc = (stepTime / a1) ; - - dgtz["hitn"] = hitn; - dgtz["sector"] = sector; - dgtz["layer"] = layer; - dgtz["component"] = strip; - dgtz["ADC_order"] = 0; - dgtz["ADC_ADC"] = (int) adc; - dgtz["ADC_time"] = convert_to_precision(stepTime); - dgtz["ADC_ped"] = 0; - - dgtz["TDC_order"] = layer < 4 ? 2 : 1; // 1 ECAL, 2 PCAL - dgtz["TDC_TDC"] = tdc; - - return dgtz; - } - - HCname = "ECAL Hit Process"; - trueInfos tInfos(aHit); - - vector Lpos = aHit->GetLPos(); - vector Edep = aHit->GetEdep(); - - // Get scintillator volume x dimension (mm) - double pDx2 = aHit->GetDetector().dimensions[5]; ///< G4Trap Semilength. - - // initialize ADC and TDC - double ADC = 0; - double ftime_in_ns = 0; - double dtime_in_ns = 0; - double ftime_in_ns_res = 0; - double dtime_in_ns_res = 0; - - // Get Total Energy deposited - double Etota = 0; - double FTtota = 0; - double DTtota = 0; - double latt = 0; - - double att; - - double A = ecc.attlen[sector-1][layer-1][0][strip-1]; - double B = ecc.attlen[sector-1][layer-1][1][strip-1]*10.; - double C = ecc.attlen[sector-1][layer-1][2][strip-1]; - double D = ecc.attlen[sector-1][layer-1][3][strip-1]; - double E = ecc.attlen[sector-1][layer-1][4][strip-1]*10.; - - double G = ecc.gain[sector-1][layer-1][strip-1]; - - double tmf = ecc.tmf_offset[sector-1][layer-1][strip-1]; - double fo = ecc.fadc_offset[sector-1][layer-1][0]; - double gtw = ecc.global_time_walk[sector-1][layer-1][0]; - - double FTOFFSET = ecc.fadc_global_offset; - double tgo = ecc.tdc_global_offset; - - double fa0 = ecc.ftime[sector-1][layer-1][0][strip-1]; - double fa2 = ecc.ftime[sector-1][layer-1][2][strip-1]; - double fa3 = ecc.ftime[sector-1][layer-1][3][strip-1]; - double fa4 = ecc.ftime[sector-1][layer-1][4][strip-1]; - double fa5 = ecc.ftime[sector-1][layer-1][5][strip-1]; - double fa6 = ecc.ftime[sector-1][layer-1][6][strip-1]; - - double da0 = ecc.dtime[sector-1][layer-1][0][strip-1]; - double da1 = ecc.dtime[sector-1][layer-1][1][strip-1]; - double da2 = ecc.dtime[sector-1][layer-1][2][strip-1]; - double da3 = ecc.dtime[sector-1][layer-1][3][strip-1]; - double da4 = ecc.dtime[sector-1][layer-1][4][strip-1]; - double da5 = ecc.dtime[sector-1][layer-1][5][strip-1]; - double da6 = ecc.dtime[sector-1][layer-1][6][strip-1]; - double da7 = ecc.dtime[sector-1][layer-1][7][strip-1]; - double da8 = ecc.dtime[sector-1][layer-1][8][strip-1]; - ; - double fveff = ecc.fveff[sector-1][layer-1][strip-1]*10; - double dveff = ecc.dveff[sector-1][layer-1][strip-1]*10; - - double fthr = ecc.fthr[sector-1][layer-1][strip-1]; - - double def0 = ecc.deff[sector-1][layer-1][0][strip-1]; - double def1 = ecc.deff[sector-1][layer-1][1][strip-1]; - double def2 = ecc.deff[sector-1][layer-1][2][strip-1]; - - - double ftres0 = ecc.ftres[sector-1][layer-1][0][0]; - double ftres1 = ecc.ftres[sector-1][layer-1][1][0]; - double ftres2 = ecc.ftres[sector-1][layer-1][2][0]; - double ftres3 = ecc.ftres[sector-1][layer-1][3][0]; - - double dtres0 = ecc.dtres[sector-1][layer-1][0][0]; - double dtres1 = ecc.dtres[sector-1][layer-1][1][0]; - double dtres2 = ecc.dtres[sector-1][layer-1][2][0]; - double dtres3 = ecc.dtres[sector-1][layer-1][3][0]; - - for(unsigned int s=0; s 3) { - // for ec, it's a minus sign - latt = pDx2-xlocal; - } else { - // for pcal, it's a plus sign - latt = pDx2+xlocal; - } - } - att = A*(exp(-latt/B) + D*exp(-latt/E)) + C; //pass2 parameterization - Etota = Etota + Edep[s]*att; //reported in MeV - FTtota = FTtota + latt/fveff; //FADC based timing - DTtota = DTtota + latt/dveff; //DSC/TDC based timing - } - - //Used if ecc.outputRAW > 0 (no digitization, effieiency or resolution smearing) - double ADC_raw = Etota/1000/ecc.ADC_GeV_to_evio/G; - double FTIME_raw = tInfos.time+FTtota/tInfos.nsteps; - double DTIME_raw = tInfos.time+DTtota/tInfos.nsteps; - - if (Etota > 0) { - double EC_npe = G4Poisson(Etota*pmtPEYld); //number of photoelectrons - if (EC_npe>0) { - double sigma = sqrt(EC_npe)*ecc.pmtFactor; - double EC_GeV = G4RandGauss::shoot(EC_npe,sigma)/1000./ecc.ADC_GeV_to_evio/G/pmtPEYld; - if (EC_GeV>0) { - ADC = EC_GeV; - double radc = sqrt(ADC); - double ftim = (fa4==0||fa6==0) ? 0 : fa0 + fa2+exp(-(radc-fa3)/fa4)+1-exp( (radc-fa5)/fa6); - double dtim = (da4==0||da6==0||da7==0) ? 0 : da0 + gtw/radc +da2+exp(-(radc-da3)/da4)+1-exp(-(da5-radc)/da6)-exp(-(radc-da3*0.95)/da7)*pow(radc,da8); - ftime_in_ns = FTIME_raw + ftim + tgo - FTOFFSET - tmf - fo; - dtime_in_ns = DTIME_raw + dtim + tgo; - ftime_in_ns_res = G4RandGauss::shoot(ftime_in_ns,getTRES(ADC,ftres0,ftres2,ftres3,ftres1)); - dtime_in_ns_res = G4RandGauss::shoot(dtime_in_ns,getTRES(ADC,dtres0,dtres2,dtres3,dtres1)); - } - } - } - - dtime_in_ns = ecc.outputRAW>0 ? DTIME_raw : dtime_in_ns_res; - ftime_in_ns = ecc.outputRAW>0 ? FTIME_raw : ftime_in_ns_res; - - // Status flags - if(accountForHardwareStatus) { - switch (ecc.status[sector-1][layer-1][strip-1]) { - case 0: - break; - case 1: - ADC = ADC_raw = ftime_in_ns = 0; - break; - case 2: - dtime_in_ns = 0; - break; - case 3: - ADC = ADC_raw = ftime_in_ns = dtime_in_ns = 0; - break; - - case 5: - break; - - default: - cout << " > Unknown EC status: " << ecc.status[sector-1][layer-1][strip-1] << " for sector " << sector << ", layer " << layer << ", strip " << strip << endl; - } - } - - if (ecc.outputRAW==0 && def0>0 && dtime_in_ns > 0 && G4UniformRand() > 1/pow(1+exp(-def0*(ADC/10-def1)),def2)) dtime_in_ns = 0; // DSC/TDC threshold - if (ecc.outputRAW==0 && ADC/10 < fthr) rejectHitConditions = true; // FADC threshold - - // EVIO banks record time with offset determined by position of data in capture window. On forward carriage this is currently - // around 7.9 us. This offset is omitted in the simulation. Also EVIO TDC time is relative to the trigger time, which is not - // simulated at present. - - double fadc_time = convert_to_precision(ftime_in_ns); - int tdc = dtime_in_ns/da1; - - dgtz["hitn"] = hitn; - dgtz["sector"] = sector; - dgtz["layer"] = layer; - dgtz["component"] = strip; - dgtz["ADC_order"] = 0; - dgtz["ADC_ADC"] = ecc.outputRAW==1 ? ADC_raw : ADC; - dgtz["ADC_time"] = fadc_time; - dgtz["ADC_ped"] = 0; - dgtz["TDC_order"] = 2; - dgtz["TDC_TDC"] = tdc; - - // define conditions to reject hit - if(rejectHitConditions) { - writeHit = false; - } - - return dgtz; -} - -vector ecal_HitProcess :: processID(vector id, G4Step* aStep, detector Detector) -{ - id[id.size()-1].id_sharing = 1; - return id; -} - -// - electronicNoise: returns a vector of hits generated / by electronics. -vector ecal_HitProcess :: electronicNoise() -{ - vector noiseHits; - - // first, identify the cells that would have electronic noise - // then instantiate hit with energy E, time T, identifier IDF: - // - // MHit* thisNoiseHit = new MHit(E, T, IDF, pid); - - // push to noiseHits collection: - // noiseHits.push_back(thisNoiseHit) - - return noiseHits; -} - -map< string, vector > ecal_HitProcess :: multiDgt(MHit* aHit, int hitn) -{ - map< string, vector > MH; - - return MH; -} - -// - charge: returns charge/time digitized information / step -// index 0: hit number -// index 1: step index -// index 2: charge -// index 3: time at electronics -// index 4: vector of identifiers - have to match the translation table -// index 5: hardware object: crate/slot/channel from translation table -map< int, vector > ecal_HitProcess :: chargeTime(MHit* aHit, int hitn) -{ - map< int, vector > CT; - - vector hitNumbers; - vector stepIndex; - vector chargeAtElectronics; - vector timeAtElectronics; - vector identifiers; - vector hardware; - hitNumbers.push_back(hitn); - - // getting identifiers - vector identity = aHit->GetId(); - - // get sector, stack (inner or outer), view (U, V, W), and strip. - // The stack/view information is lost in the identifier but it can be recovered by looking at the volume description - // For example: - // grep "sector manual 1" ec__geometry_rga_fall2018.txt | awk -F\| '{print $2" "$18}' | grep "strip manual 1 " - int sector = identity[0].id; - int layer = identity[1].id; // layer=1-3 (PCAL) 4-9 (ECAL). Layer = view for pcal, ecinner, ecouter - int strip = identity[2].id; - - // pcal - int view = layer; - - if (layer > 3 && layer < 7) { - // ec inner (stack 1) - view = layer - 3; - } else if (layer > 7) { - // ec inner (stack 2) - view = layer - 6; - } - - - // Number of p.e. divided by the energy deposited in MeV. See EC NIM paper table 1. - // Different for EC and PCAL - double pmtPEYld = 3.5 ; - if (layer < 4) { - // pcal - pmtPEYld = 11.5 ; - } - - identifiers.push_back(sector); // sector - identifiers.push_back(layer); // laylayer=1-3 (PCAL) 4-9 (ECAL)er - identifiers.push_back(strip); // component (pmt) - identifiers.push_back(0); // order - - // getting hardware - Hardware thisHardware = ecc.TT.getHardware({sector, layer, strip, 0}); - hardware.push_back(thisHardware.getCrate()); - hardware.push_back(thisHardware.getSlot()); - hardware.push_back(thisHardware.getChannel()); - - // Adding pedestal mean and sigma into the hardware as well - // All of these variables start from 1, therefore -1 is subtracted, e.g. sector-1 - hardware.push_back(ecc.pedestal[sector - 1][layer - 1][view - 1]); - hardware.push_back(ecc.pedestal_sigm[sector - 1][layer - 1][view - 1]); - - // getting charge and time - trueInfos tInfos(aHit); - - // Get scintillator mother volume dimensions (mm) - double pDx2 = aHit->GetDetector().dimensions[5]; ///< G4Trap Semilength. - - vector pos = aHit->GetPos(); - vector Lpos = aHit->GetLPos(); - - - vector Edep = aHit->GetEdep(); - vector time = aHit->GetTime(); - - double A = ecc.attlen[sector-1][layer-1][0][strip-1]; - double B = ecc.attlen[sector-1][layer-1][1][strip-1]*10.; - double C = ecc.attlen[sector-1][layer-1][2][strip-1]; - double D = ecc.attlen[sector-1][layer-1][3][strip-1]; - double E = ecc.attlen[sector-1][layer-1][4][strip-1]*10.; - - double G = ecc.gain[sector-1][layer-1][strip-1]; - double fveff = ecc.fveff[sector-1][layer-1][strip-1]*10; - - for(unsigned int s=0; s0) { - double xlocal = Lpos[s].x(); - double latt = 0; - - if(view==1) latt = pDx2+xlocal; - if(view==2) latt = pDx2+xlocal; - if(view==3) { - if(layer > 3) { - // for ecal, it's a minus sign - latt = pDx2-xlocal; - } else { - // for pcal, it's a plus sign - latt = pDx2+xlocal; - } - } - - double att = A*(exp(-latt/B) + D*exp(-latt/E)) + C; //pass2 parameterization - - double stepE = Edep[s]*att; - double stepTime = time[s] + latt/fveff; - - // cout<<"time[s] = "< Initializing " << HCname << " digitization for run number " << runno << endl; - ecc = initializeECConstants(runno, digiVariation, digiSnapshotTime, accountForHardwareStatus); - ecc.runNo = runno; - } -} - -// this static function will be loaded first thing by the executable -ecConstants ecal_HitProcess::ecc = initializeECConstants(-1); diff --git a/source/hitprocess/clas12/ecal_hitprocess_rr.h b/source/hitprocess/clas12/ecal_hitprocess_rr.h deleted file mode 100644 index a43af38d..00000000 --- a/source/hitprocess/clas12/ecal_hitprocess_rr.h +++ /dev/null @@ -1,131 +0,0 @@ -#ifndef ECAL_HITPROCESS_H -#define ECAL_HITPROCESS_H 1 - -// gemc headers -#include "HitProcess.h" - -// ec constants -// these are loaded with initWithRunNumber -class ecConstants -{ -public: - // runNo is mandatory variable to keep track of run number changes - int runNo; - string date; - string connection; - char database[80]; - int outputRAW=0; //bypass calibration, resolution and efficiency - - static const int nsect = 6; // Number of sectors - static const int nlayer = 9; // layer=1-3 (PCAL) 4-6 (ECinner) 7-9 (ECouter) - static const int nview = 3; // Number of views, U,V and W - - // For strip dependent constants read from CCDB - // Array [6][9][3] -> sector,layer,view sector=1-6 layer=1-3 (PCAL) 4-6 (ECinner) 7-9 (ECouter) view=1-3 (U,V,W) - - //attlen: attenuation length - vector attlen[nsect][nlayer][10]; - - //gain: pmt gain - vector gain[nsect][nlayer]; - - //timing: TDC calibration constants - vector dtime[nsect][nlayer][9]; - vector ftime[nsect][nlayer][7]; - - vector fadc_offset[nsect][nlayer]; - vector tmf_offset[nsect][nlayer]; - vector global_time_walk[nsect][nlayer]; - - double tdc_global_offset; - double fadc_global_offset; - - //veff: effective velocity (cm/ns) - vector dveff[nsect][nlayer]; - vector fveff[nsect][nlayer]; - - // DSC/TSC efficiency - vector deff[nsect][nlayer][3]; - - // FADC threshold - vector fthr[nsect][nlayer]; - - // Timiing resolution - vector dtres[nsect][nlayer][4]; - vector ftres[nsect][nlayer][4]; - - // status: - // 0 - fully functioning - // 1 - noADC - // 2 - noTDC - // 3 - noADC, noTDC (PMT is dead) - // 5 - any other reconstruction problem - vector status[nsect][nlayer]; - - // ======== FADC Pedestals and sigmas =========== - double pedestal[nsect][nlayer][nview] = {}; - double pedestal_sigm[nsect][nlayer][nview] = {}; - - // double TDC_time_to_evio; // Conversion from time (ns) to TDC format - double ADC_GeV_to_evio; // Conversion from energy (GeV) to FADC250 format - double pmtQE; // Quantum efficiency of PMT - double pmtDynodeGain; // PMT dynode gain - double pmtDynodeK; // PMT dynode secondary emission statistics factor: K=0 (Poisson) K=1 (exponential) - double pmtFactor; // Contribution to FWHM from PMT statistical fluctuations. - - // voltage signal parameters, using double gaussian + delay (function DGauss, need documentation for it) - double vpar[4]; - - // translation table - TranslationTable TT; -}; - - -// Class definition -class ecal_HitProcess : public HitProcess -{ -public: - - ~ecal_HitProcess(){;} - - // constants initialized with initWithRunNumber - static ecConstants ecc; - - void initWithRunNumber(int runno); - - // - integrateDgt: returns digitized information integrated over the hit - map integrateDgt(MHit*, int); - - // - multiDgt: returns multiple digitized information / hit - map< string, vector > multiDgt(MHit*, int); - - // - charge: returns charge/time digitized information / step - virtual map< int, vector > chargeTime(MHit*, int); - - // - voltage: returns a voltage value for a given time. The input are charge value, time - virtual double voltage(double, double, double); - - // The pure virtual method processID returns a (new) identifier - // containing hit sharing information - vector processID(vector, G4Step*, detector); - - // creates the HitProcess - static HitProcess *createHitClass() {return new ecal_HitProcess;} - - // - electronicNoise: returns a vector of hits generated / by electronics. - vector electronicNoise(); - -private: - - double fadc_precision = 0.0625; // 62 picoseconds resolution - double convert_to_precision(double time) { - return (int( time / fadc_precision ) * fadc_precision); - } - - double getTRES(double x, double p0, double p1, double p2, double p3) { - return (p0*exp(pow(x,p3)/120000) + p1/x + p2/pow(x,0.5)); - } - -}; - -#endif From f5c4fb8dbf64d4d42cecc74e13d154bae1d70797 Mon Sep 17 00:00:00 2001 From: raffaelladevita Date: Thu, 9 Apr 2026 13:47:24 -0400 Subject: [PATCH 08/12] changed sign of RF offsets, so that the numbers in CCDB are effectively subtracted (since they are summed in reconstruction) --- source/src/MEventAction.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/source/src/MEventAction.cc b/source/src/MEventAction.cc index 83782f93..b2c5cd90 100644 --- a/source/src/MEventAction.cc +++ b/source/src/MEventAction.cc @@ -1135,8 +1135,8 @@ void MEventAction::setup_clas12_RF(int runno) { cout << " Connecting to " << database << endl; data.clear(); calib->GetCalib(data, database); - double offset1 = data[0][3]; - double offset2 = data[1][3]; + double offset1 = -data[0][3]; + double offset2 = -data[1][3]; rfvalue_strings = {to_string(clock), to_string(prescale), to_string(offset1), "20", to_string(offset2)}; set_and_show_rf_setup(); From 4dd7d2d97fbcd1d1fd750b21ce4eb0b7af9adac4 Mon Sep 17 00:00:00 2001 From: raffaelladevita Date: Mon, 13 Apr 2026 16:01:24 -0400 Subject: [PATCH 09/12] Reading BMT drift HV from CCDB, as in reconstruction --- .../clas12/micromegas/BMT_hitprocess.cc | 23 ++++++++++++------- 1 file changed, 15 insertions(+), 8 deletions(-) diff --git a/source/hitprocess/clas12/micromegas/BMT_hitprocess.cc b/source/hitprocess/clas12/micromegas/BMT_hitprocess.cc index b0dd1375..90636c98 100644 --- a/source/hitprocess/clas12/micromegas/BMT_hitprocess.cc +++ b/source/hitprocess/clas12/micromegas/BMT_hitprocess.cc @@ -80,13 +80,19 @@ static bmtConstants initializeBMTConstants(int runno, string digiVariation = "de } } } - - for (int j = 0; j GetCalib(data,bmtc.database); + for(unsigned row = 0; row < data.size(); row++) + { + int sector = data[row][0]; + int layer = data[row][1]; + bmtc.HV_DRIFT[layer-1][sector-1] = data[row][3]; + bmtc.HV_STRIPS[layer-1][sector-1]=520; } // all dimensions are in mm @@ -234,7 +240,8 @@ vector BMT_HitProcess :: processID(vector id, G4Step* a } double depe = aStep->GetTotalEnergyDeposit(); - + if(bmtc.HV_DRIFT[layer-1][sector-1]==0) depe=0; + // resetting depe for geantinos if (aStep->GetTrack()->GetDefinition() == G4ChargedGeantino::ChargedGeantinoDefinition() ){ int np=G4Poisson( (aStep->GetStepLength()/cm) *1e4); // Warning... StepLength must be in cm... because it is 10 e- per cm for MIP From b3f56f1486552c6dc2cd44c4d50077bc3a07bcd4 Mon Sep 17 00:00:00 2001 From: Veronique Ziegler Date: Tue, 28 Apr 2026 16:46:04 -0400 Subject: [PATCH 10/12] added beta-dependent timewalk contribution to time-2-distance --- source/hitprocess/clas12/dc_hitprocess.cc | 25 +++++++++++++++++++++-- source/hitprocess/clas12/dc_hitprocess.h | 2 +- 2 files changed, 24 insertions(+), 3 deletions(-) diff --git a/source/hitprocess/clas12/dc_hitprocess.cc b/source/hitprocess/clas12/dc_hitprocess.cc index f50f13f9..aeb43b39 100644 --- a/source/hitprocess/clas12/dc_hitprocess.cc +++ b/source/hitprocess/clas12/dc_hitprocess.cc @@ -99,7 +99,8 @@ static dcConstants initializeDCConstants(int runno, string digiVariation = "defa dcc.v0[sec][sl] = data[row][3] + data[row][4]*dpressure + data[row][5]*dpressure*dpressure; dcc.vmid[sec][sl] = data[row][6] + data[row][7]*dpressure + data[row][8]*dpressure*dpressure; dcc.tmaxsuperlayer[sec][sl] = data[row][9] + data[row][10]*dpressure + data[row][11]*dpressure*dpressure; - // Row left out, corresponds to distbeta + // Row corresponding to distbeta + dcc.distbeta[sec][sl] = data[row][12] + data[row][13]*dpressure + data[row][14]*dpressure*dpressure; dcc.delta_bfield_coefficient[sec][sl] = data[row][15] + data[row][16]*dpressure + data[row][17]*dpressure*dpressure; dcc.deltatime_bfield_par1[sec][sl] = data[row][18] + data[row][19]*dpressure + data[row][20]*dpressure*dpressure; dcc.deltatime_bfield_par2[sec][sl] = data[row][21] + data[row][22]*dpressure + data[row][23]*dpressure*dpressure; @@ -373,6 +374,12 @@ map dc_HitProcess :: integrateDgt(MHit* aHit, int hitn) // unsmeared time, based on the dist-time-function and alpha; double unsmeared_time = calc_Time(doca/cm,dcc.dmaxsuperlayer[SLI],dcc.tmaxsuperlayer[SECI][SLI],alpha,thisMgnf,SECI,SLI); + // add beta-dependent time-walk: + double beta_timewalk = calc_TimeBeta(doca/cm, beta_particle, SECI, SLI); + + // add it to the unsmeared time: + unsmeared_time += beta_timewalk; + // Include time smearing calculated from doca resolution double dt_random_in = doca_smearing(X, beta_particle, SECI, SLI); //double dt_random = dt_random_in*CLHEP::RandLandau::shoot(); @@ -555,7 +562,7 @@ double dc_HitProcess :: calc_Time_exp(double x, double dmax, double tmax, double return rtime; } -// NEW Polynomial function: returns a time in ns give: +// NEW Polynomial function: returns a time in ns given: // x = distance from the wire, in cm // dmax = cell size in superlayer // tmax = t max in superlayer @@ -597,6 +604,20 @@ double dc_HitProcess :: calc_Time(double x, double dmax, double tmax, double alp return time; } + +// Beta-dependent timewalk: returns a time in ns given: +// x = distance from the wire, in cm +// beta = beta of the particle +// sector = sector +// superlayer = superlayer +double dc_HitProcess :: calc_TimeBeta(double x, double beta, int sector, int superlayer) +{ double v0 = dcc.v0[sector][superlayer]; + double distbeta = dcc.distbeta[sector][superlayer]; + + double time = (0.5*pow(beta*beta*distbeta,3)*x/(pow(beta*beta*distbeta,3)+x*x*x))/v0; + + return time; +} // Define DOCA smearing based on data parameterization // x: distance from the wire normalized to the cell size // beta: beta of the particle diff --git a/source/hitprocess/clas12/dc_hitprocess.h b/source/hitprocess/clas12/dc_hitprocess.h index 592654a1..342c2cdb 100644 --- a/source/hitprocess/clas12/dc_hitprocess.h +++ b/source/hitprocess/clas12/dc_hitprocess.h @@ -43,7 +43,7 @@ class dcConstants //parameters for time to distance: double deltanm[6][6], v0[6][6], delta_bfield_coefficient[6][6],tmaxsuperlayer[6][6]; double deltatime_bfield_par1[6][6], deltatime_bfield_par2[6][6], deltatime_bfield_par3[6][6], deltatime_bfield_par4[6][6]; - double vmid[6][6], R[6][6]; + double vmid[6][6], R[6][6], distbeta[6][6]; double dmaxsuperlayer[6]; // sector, SL, slot, cable From 4cd3108ab6c51f60fe0503a821bfacfb10ddd2b7 Mon Sep 17 00:00:00 2001 From: raffaelladevita Date: Thu, 30 Apr 2026 16:35:22 -0400 Subject: [PATCH 11/12] added missing function prototype --- source/hitprocess/clas12/dc_hitprocess.h | 3 +++ 1 file changed, 3 insertions(+) diff --git a/source/hitprocess/clas12/dc_hitprocess.h b/source/hitprocess/clas12/dc_hitprocess.h index 26dbc87a..1c443d23 100644 --- a/source/hitprocess/clas12/dc_hitprocess.h +++ b/source/hitprocess/clas12/dc_hitprocess.h @@ -120,6 +120,9 @@ class dc_HitProcess : public HitProcess // returns a time given a distance: neew polynomial function double calc_Time(double x, double dmax, double tmax, double alpha, double bfield, int sector, int superlayer); + // beta-dependent ttime walk + double calc_TimeBeta(double x, double beta, int sector, int superlayer); + // returns time walks according to ionisation process: double doca_smearing(double x, double beta, int sector, int superlayer); From ba3444f697fcef557ac40cc313293046bdd1fdfe Mon Sep 17 00:00:00 2001 From: raffaelladevita Date: Tue, 7 Jul 2026 18:21:49 -0400 Subject: [PATCH 12/12] fixed bug in CTOF hposbin --- source/hitprocess/clas12/ctof_hitprocess.cc | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/source/hitprocess/clas12/ctof_hitprocess.cc b/source/hitprocess/clas12/ctof_hitprocess.cc index 5c0d59cf..eb16eb40 100644 --- a/source/hitprocess/clas12/ctof_hitprocess.cc +++ b/source/hitprocess/clas12/ctof_hitprocess.cc @@ -144,7 +144,7 @@ static ctofConstants initializeCTOFConstants(int runno, string digiVariation = " isec = data[row][0]; ilay = data[row][1]; for(unsigned bin=0; bin<=ctc.hposbins; bin++) { - ctc.hpos[isec-1][ilay-1][bin].push_back(data[row][3+bin]); + ctc.hposbin[isec-1][ilay-1][bin].push_back(data[row][3+bin]); } } @@ -214,6 +214,7 @@ static ctofConstants initializeCTOFConstants(int runno, string digiVariation = " calib->GetCalib(data, ctc.database); for (unsigned row = 0; row < data.size(); row++) { double sigma = data[row][3]; + ctc.tres.push_back(sigma); } @@ -376,10 +377,11 @@ map ctof_HitProcess::integrateDgt(MHit* aHit, int hitn) //double B = ctc.twlk[sector-1][layer-1][1][paddle-1]; //double C = ctc.twlk[sector-1][layer-1][2][paddle-1]; //double timeWalkUp = A/(B+C*sqrt(adcu)); - double hpos = ctc.hpos[sector-1][layer-1][0][paddle-1] * exp(ctc.hpos[sector-1][layer-1][1][paddle-1]/cm * y); + double ylab = tInfos.z; + double hpos = ctc.hpos[sector-1][layer-1][0][paddle-1] * exp(ctc.hpos[sector-1][layer-1][1][paddle-1]/cm * ylab); if(ctc.hposbin[sector-1][layer-1][0][paddle-1]>0) { - int bin = (int) floor(y/cm) + ctc.hposbins/2; - hpos += ctc.hposbin[sector-1][layer-1][bin+1][paddle-1]; + int bin = (int) floor(ylab/cm) + ctc.hposbins/2; + hpos = ctc.hposbin[sector-1][layer-1][bin+1][paddle-1]; } double tU = tInfos.time + d/ctc.veff[sector-1][layer-1][side][paddle-1]/cm + (1. - 2. * side)*ctc.toff_UD[sector-1][layer-1][paddle-1]/2.