diff --git a/geometry_source/recoil_tof/GeoArgParse.groovy b/geometry_source/recoil_tof/GeoArgParse.groovy new file mode 100644 index 00000000..0d48250d --- /dev/null +++ b/geometry_source/recoil_tof/GeoArgParse.groovy @@ -0,0 +1,18 @@ + +class GeoArgParse { + static String defaultVariation="default"; + static int defaultRunNumber=11; + static String getVariation(String[] args) { + for (int ii=0; ii + w<{$vname}; + $detector{"description"} = "CLAS12 recoil TOF, Sector $s"; + $detector{"pos"} = $positions->{$vname}; + $detector{"rotation"} = $rotations->{$vname}; + $detector{"type"} = $types->{$vname}; + $detector{"dimensions"} = $dimensions->{$vname}; + $detector{"color"} = "aa0000"; + $detector{"material"} = "G4_Galactic"; + $detector{"style"} = 0; + $detector{"visible"} = 0; + print_det(\%main::configuration, \%detector); + + for(my $row=1; $row<=$nRows; $row++) { + for(my $col=1; $col<=$nCols; $col++) { + make_bar($s, $row, $col); + } + } + } + +} + +sub make_bar +{ + my $iSector = $_[0]; + my $iRow = $_[1]; + my $iCol =$_[2]; + + my %detector = init_det(); + + my $vname = "bar_sector${iSector}_row${iRow}_column${iCol}"; + + $detector{"name"} = $vname; + $detector{"mother"} = $mothers->{$vname}; + $detector{"description"} = "Sector $iSector, Bar, row $iRow, column $iCol"; + $detector{"pos"} = $positions->{$vname}; + $detector{"rotation"} = $rotations->{$vname}; + $detector{"type"} = $types->{$vname}; + $detector{"dimensions"} = $dimensions->{$vname}; + $detector{"color"} = "aa0000"; + $detector{"material"} = "scintillator"; + $detector{"style"} = 0; + $detector{"visible"} = 1; + $detector{"sensitivity"} = "recoil_tof"; + $detector{"hit_type"} = "recoil_tof"; + # set the identifiers + $detector{"identifiers"} = "sector manual $iSector row manual $iRow column manual $iCol order manual 0"; + print_det(\%main::configuration, \%detector); + +} + + +sub makeRECOILTOF +{ + + ($mothers, $positions, $rotations, $types, $dimensions, $ids) = @main::volumes; + + make_sector(); + + +} + +1; diff --git a/geometry_source/recoil_tof/hit.pl b/geometry_source/recoil_tof/hit.pl new file mode 100644 index 00000000..18c9ad85 --- /dev/null +++ b/geometry_source/recoil_tof/hit.pl @@ -0,0 +1,23 @@ +use strict; +use warnings; + +our %configuration; + +sub define_hit +{ + # uploading the hit definition + my %hit = init_hit(); + $hit{"name"} = "recoil_tof"; + $hit{"description"} = "recoil tof hit definitions "; + $hit{"identifiers"} = "sector row column"; + $hit{"signalThreshold"} = "0.5*MeV"; + $hit{"timeWindow"} = "400*ns"; + $hit{"prodThreshold"} = "1*mm"; + $hit{"maxStep"} = "1*cm"; + $hit{"delay"} = "50*ns"; + $hit{"riseTime"} = "0.7*ns"; + $hit{"fallTime"} = "1.8*ns"; + $hit{"mvToMeV"} = 100; + $hit{"pedestal"} = -20; + print_hit(\%configuration, \%hit); +} diff --git a/geometry_source/recoil_tof/materials.pl b/geometry_source/recoil_tof/materials.pl new file mode 100755 index 00000000..fce03b19 --- /dev/null +++ b/geometry_source/recoil_tof/materials.pl @@ -0,0 +1,21 @@ +use strict; +use warnings; + +our %configuration; + + +sub materials +{ + # uploading the mat definition + + # Scintillator + my %mat = init_mat(); + $mat{"name"} = "scintillator"; + $mat{"description"} = "recoil tof scintillator material"; + $mat{"density"} = "1.032"; + $mat{"ncomponents"} = "2"; + $mat{"components"} = "C 9 H 10"; + print_mat(\%configuration, \%mat); + +} + diff --git a/geometry_source/recoil_tof/recoil_tof.gcard b/geometry_source/recoil_tof/recoil_tof.gcard new file mode 100644 index 00000000..64ccac45 --- /dev/null +++ b/geometry_source/recoil_tof/recoil_tof.gcard @@ -0,0 +1,29 @@ + + + diff --git a/geometry_source/recoil_tof/recoil_tof.pl b/geometry_source/recoil_tof/recoil_tof.pl new file mode 100755 index 00000000..715a3095 --- /dev/null +++ b/geometry_source/recoil_tof/recoil_tof.pl @@ -0,0 +1,91 @@ +#!/usr/bin/perl -w + +use strict; +use lib ("$ENV{GEMC}/api/perl"); +use utils; +use parameters; +use geometry; +use hit; +use bank; +use math; +use materials; + +use Math::Trig; + +# Help Message +sub help() +{ + print "\n Usage: \n"; + print " recoil_tof.pl \n"; + print " Will create the CLAS12 recoil tof geometry, materials, bank and hit definitions\n"; + exit; +} + +# Make sure the argument list is correct +if( scalar @ARGV != 1) +{ + help(); + exit; +} + +# Loading configuration file and paramters +our %configuration = load_configuration($ARGV[0]); +#$configuration{"variation"} = "default" ; + +# Global pars - these should be read by the load_parameters from file or DB + + +# materials +require "./materials.pl"; + +# banks definitions +require "./bank.pl"; + +# hits definitions +require "./hit.pl"; + +#system("which run-groovy"); + +# sensitive geometry +require "./geometry_java.pl"; + + +# all the scripts must be run for every configuration +# Right now run both configurations, later on just ccdb +#my @allConfs = ("ccdb", "cosmicR1", "ddvcs", "java"); +#my @allConfs = ("default"); + +my @allConfs = ("default"); + +# bank definitions common to all variations +define_banks(); + +#foreach my $conf ( @allConfs ){ +for(my $ii=0; $ii<=$#allConfs; $ii++){ +# $configuration{"variation"} = $conf ; + $configuration{"variation"} = $allConfs[$ii]; + print "$configuration{variation}\n"; + + system("/w/hallb-scshelf2102/clas12/users/nwickjlb/coatjava_edit/coatjava/coatjava/bin/run-groovy factory.groovy --variation $configuration{variation} --runnumber 11"); + #system("groovy -cp '../*:..' factory.groovy --variation $configuration{variation} --runnumber 11"); + + # %configuration = $configuration{"variation"}; + + # materials + materials(); + + # hits + define_hit(); + + + # sensitive geometry + # Global pars - these should be read by the load_parameters from file or DB + our @volumes = get_volumes(%configuration); + + + print "variation default\n"; + print "configuration{$configuration{variation})\n"; + coatjava::makeRECOILTOF(); + +} + diff --git a/source/SConstruct b/source/SConstruct index 12b75f5a..170c38e6 100644 --- a/source/SConstruct +++ b/source/SConstruct @@ -109,7 +109,7 @@ env.Library(source = field_sources, target = "lib/gfields") # Hit Processes -env.Append(CPPPATH = ['hitprocess', 'hitprocess/clas12', 'hitprocess/clas12/svt', 'hitprocess/clas12/micromegas', 'hitprocess/clas12/alert', 'hitprocess/clas12/recoil']) +env.Append(CPPPATH = ['hitprocess', 'hitprocess/clas12', 'hitprocess/clas12/svt', 'hitprocess/clas12/micromegas', 'hitprocess/clas12/alert', 'hitprocess/clas12/recoil', 'hitprocess/clas12/recoil_tof']) env.Append(CPPPATH = ['hitprocess/Aprime', 'hitprocess/GlueX', 'hitprocess/solid']) hitp_sources = Split(""" hitprocess/HitProcess_MapRegister.cc @@ -142,8 +142,9 @@ hitp_sources = Split(""" hitprocess/clas12/band_hitprocess.cc hitprocess/clas12/uRwell/uRwell_hitprocess.cc hitprocess/clas12/uRwell/uRwell_strip.cc - hitprocess/clas12/recoil/recoil_hitprocess.cc - hitprocess/clas12/recoil/recoil_strip.cc + hitprocess/clas12/recoil/recoil_hitprocess.cc + hitprocess/clas12/recoil/recoil_strip.cc + hitprocess/clas12/recoil_tof/recoil_tof_hitprocess.cc hitprocess/bdx/cormo_hitprocess.cc hitprocess/bdx/veto_hitprocess.cc hitprocess/bdx/crs_hitprocess.cc diff --git a/source/hitprocess/HitProcess_MapRegister.cc b/source/hitprocess/HitProcess_MapRegister.cc index 0ddacd45..6143d4d5 100644 --- a/source/hitprocess/HitProcess_MapRegister.cc +++ b/source/hitprocess/HitProcess_MapRegister.cc @@ -27,6 +27,7 @@ #include "clas12/rtpc_hitprocess.h" ///< Radial Time Projection Chamber (RTPC) #include "clas12/uRwell/uRwell_hitprocess.h" ///< uRwell detector #include "clas12/recoil/recoil_hitprocess.h" ///< RGH recoil detector +#include "clas12/recoil_tof/recoil_tof_hitprocess.h" ///< RGH recoil TOF detector // Beam Dump eXperiment #include "bdx/cormo_hitprocess.h" ///< Cormorino detector @@ -93,6 +94,7 @@ map HitProcess_Map(string experiments) hitMap["rtpc"] = &rtpc_HitProcess::createHitClass; hitMap["urwell"] = &uRwell_HitProcess::createHitClass; hitMap["recoil"] = &recoil_HitProcess::createHitClass; + hitMap["recoil_tof"] = &recoil_tof_HitProcess::createHitClass; } else if(EXP == "HPS") { // Aprime hitMap["SVT"] = &SVT_HitProcess::createHitClass; diff --git a/source/hitprocess/clas12/recoil_tof/recoil_tof_hitprocess.cc b/source/hitprocess/clas12/recoil_tof/recoil_tof_hitprocess.cc new file mode 100644 index 00000000..976d36db --- /dev/null +++ b/source/hitprocess/clas12/recoil_tof/recoil_tof_hitprocess.cc @@ -0,0 +1,298 @@ +// G4 headers +#include "G4Poisson.hh" +#include "Randomize.hh" + +#include +#include +#include +using namespace ccdb; + +// gemc headers +#include "recoil_tof_hitprocess.h" + +// CLHEP units +#include "CLHEP/Units/PhysicalConstants.h" +using namespace CLHEP; + +static recoil_tofConstants initializeRECOILTOFConstants(int runno, string digiVariation = "default") { + recoil_tofConstants atc; + + // 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 atc; + + atc.runNo = runno; + if (getenv("CCDB_CONNECTION") != nullptr) + atc.connection = (string) getenv("CCDB_CONNECTION"); + else + atc.connection = "mysql://clas12reader@clasdb.jlab.org/clas12"; + + return atc; +} + + +// this methos is for implementation of digitized outputs and the first one that needs to be implemented. +map recoil_tof_HitProcess::integrateDgt(MHit* aHit, int hitn) { + + // digitized output + map dgtz; + vector identity = aHit->GetId(); + rejectHitConditions = false; + writeHit = true; + + + int recoil_tof_sector = identity[0].id; + int recoil_tof_row = identity[1].id; + int recoil_tof_column = identity[2].id; + int recoil_tof_order = identity[3].id; + //cout << "sector " << recoil_tof_sector << " row "<< recoil_tof_row << " column " << recoil_tof_column << " order " << recoil_tof_order << endl; + double time_to_tdc = 1./0.015625; + + if(aHit->isBackgroundHit == 1) { + + double totEdep = aHit->GetEdep()[0]; + double stepTime = aHit->GetTime()[0]; + double tdc = stepTime * time_to_tdc; + + dgtz["hitn"] = hitn; + dgtz["sector"] = recoil_tof_sector; //Sector ranges from 1 to 2 + dgtz["row"] = recoil_tof_row; //Row ranges from 1 to 5 + dgtz["column"] = recoil_tof_column; //Column ranges from 1 to 63 + dgtz["TDC_order"] = recoil_tof_order; //order is 0/1 for front(upstream)/back(downstream) + dgtz["TDC_ToT"] = (int) totEdep; + dgtz["TDC_TDC"] = tdc; + return dgtz; + } + + trueInfos tInfos(aHit); + + double length = aHit->GetDetector().dimensions[1]; // half length of bar along y direction + double thickness = 2 * aHit->GetDetector().dimensions[2]; // thickness of bar + //cout << "half length = " << length << " mm " << "thickness = " << thickness << " mm" << endl; + + vector Edep = aHit->GetEdep(); + vector Lpos = aHit->GetLPos(); // local position at each step + vector times = aHit->GetTime(); + + double adc_CC_front, adc_CC_back, tdc_CC_front, tdc_CC_back; + + double LposX=0.0; + double LposY=0.0; + double LposZ=0.0; + + // Simple output not equal to real physics, just to feel the adc, time values + // Should be: double energy = tInfos.eTot*att; + + double totEdep=0.0; + + //Distance calculation from the hit to the front or back SIPM + double dFront = 0.0; + double dBack = 0.0; + double e_Front = 0.0; + double e_Back = 0.0; + double E_tot_Front = 0.0; + double E_tot_Back = 0.0; + + double attlength = 1600.0; // here in mm! because all lengths from the volume are in mm!! EJ-204 160 cm + double pmtPEYld = 10400.0; // EJ-204 10400 (photons / [1MeV*e-]) + double dEdxMIP = 1.956; // energy deposited by MIP per cm of scintillator material, to be adapted for SiPM case, it is a function of ? + + //Variables for tdc calculation (time) + double EtimesTime_Front=0.0; + double EtimesTime_Back=0.0; + double v_eff_Front = 200.0; // mm/ns! CND v_eff = 16 cm/ns + double v_eff_Back = 200.0; + + /* + double dist_h_SiPMFront =0.0; + double dist_h_SiPMBack =0.0; + */ + // cout << "First loop on steps begins" << endl; + + + // notice these calculations are done both for front and back + // this can be optimized to have just one calculation using order as discriminating value + for(unsigned int s=0; s 0.0) || (E_tot_Back > 0.0)) + { + double nphe_fr = G4Poisson(E_tot_Front*pmtPEYld); + double energy_fr = nphe_fr/pmtPEYld; + + double nphe_bck = G4Poisson(E_tot_Back*pmtPEYld); + double energy_bck = nphe_bck/pmtPEYld; + + adc_front = energy_fr *adc_CC_front *(1/(dEdxMIP * thickness * 0.1)); // thickness of bar (mm -> cm) + adc_back = energy_bck *adc_CC_back *(1/(dEdxMIP * thickness * 0.1)); + + + time_front = EtimesTime_Front/E_tot_Front; + time_back = EtimesTime_Back/E_tot_Back; + tdc_front = G4RandGauss::shoot(time_front, sigma_time) / tdc_CC_front; + tdc_back = G4RandGauss::shoot(time_back, sigma_time) / tdc_CC_back; + } + + double adc = 0; + double time = 0; + + if ( recoil_tof_order == 0 ) + { + adc = adc_front; + time = tdc_front; + } + else + { + adc = adc_back; + time = tdc_back; + } + + + dgtz["hitn"] = hitn; + dgtz["sector"] = recoil_tof_sector; //Sector ranges from 1 to 2 + dgtz["row"] = recoil_tof_row; //Row ranges from 1 to 5 + dgtz["column"] = recoil_tof_column; //Column ranges from 1 to 63 + dgtz["TDC_order"] = recoil_tof_order; //order is 0/1 for front(upstream)/back(downstream) + dgtz["TDC_ToT"] = (int)adc*100; + dgtz["TDC_TDC"] = time * time_to_tdc; + + // define conditions to reject hit + if (rejectHitConditions) { + writeHit = false; + } + + return dgtz; +} + +vector recoil_tof_HitProcess::processID(vector id, G4Step* aStep, detector Detector) { + + vector yid = id; + + yid[0].id_sharing = 1; // sector + yid[1].id_sharing = 1; // row + yid[2].id_sharing = 1; // column + yid[3].id_sharing = 1; // order + + if (yid[3].id != 0) { + cout << "*****WARNING***** in recoil_tof_HitProcess :: processID, order of the original hit should be 0 " << endl; + cout << "yid[3].id = " << yid[3].id << endl; + } + + // Now we want to have similar identifiers, but the only difference be id order to be 1, instead of 0 + identifier this_id = yid[0]; + yid.push_back(this_id); + this_id = yid[1]; + yid.push_back(this_id); + this_id = yid[2]; + yid.push_back(this_id); + this_id = yid[3]; + this_id.id = 1; + yid.push_back(this_id); + + return yid; +} + +// - electronicNoise: returns a vector of hits generated / by electronics. + +vector recoil_tof_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 > recoil_tof_HitProcess::multiDgt(MHit* aHit, int hitn) { + map< string, vector > MH; + + return MH; +} + +// - charge: returns charge/time digitized information / step + +map< int, vector > recoil_tof_HitProcess::chargeTime(MHit* aHit, int hitn) { + map< int, vector > CT; + + return CT; +} + +// - voltage: returns a voltage value for a given time. The inputs are: +// charge value (coming from chargeAtElectronics) +// time (coming from timeAtElectronics) + +double recoil_tof_HitProcess::voltage(double charge, double time, double forTime) { + return 0.0; +} + +void recoil_tof_HitProcess::initWithRunNumber(int runno) +{ + string digiVariation = gemcOpt.optMap["DIGITIZATION_VARIATION"].args; + + if (atc.runNo != runno) { + cout << " > Initializing " << HCname << " digitization for run number " << runno << endl; + atc = initializeRECOILTOFConstants(runno, digiVariation); + atc.runNo = runno; + } +} + +// this static function will be loaded first thing by the executable +recoil_tofConstants recoil_tof_HitProcess::atc = initializeRECOILTOFConstants(-1); + + + + diff --git a/source/hitprocess/clas12/recoil_tof/recoil_tof_hitprocess.h b/source/hitprocess/clas12/recoil_tof/recoil_tof_hitprocess.h new file mode 100644 index 00000000..1bb1797d --- /dev/null +++ b/source/hitprocess/clas12/recoil_tof/recoil_tof_hitprocess.h @@ -0,0 +1,68 @@ +#ifndef MYRECOIL_TOF_HITPROCESS_H +#define MYRECOIL_TOF_HITPROCESS_H 1 + +// gemc headers +#include "HitProcess.h" + + +class recoil_tofConstants +{ +public: + + // Database parameters + int runNo; + string date; + string connection; + char database[80]; + + // translation table + TranslationTable TT; +}; + + +// Class definition +/// \class recoil_tof_HitProcess +/// Recoil Time of Flight Hit Process Routine\n\n + +class recoil_tof_HitProcess : public HitProcess +{ +public: + + ~recoil_tof_HitProcess(){;} + + // constants initialized with initWithRunNumber + static recoil_tofConstants atc; + + 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 recoil_tof_HitProcess;} + + // - electronicNoise: returns a vector of hits generated / by electronics. + vector electronicNoise(); + +}; + + + +#endif + + + + diff --git a/source/output/hipoSchemas.cc b/source/output/hipoSchemas.cc index a45669e4..71963a74 100644 --- a/source/output/hipoSchemas.cc +++ b/source/output/hipoSchemas.cc @@ -67,6 +67,7 @@ HipoSchema::HipoSchema() { helONLINESchema = hipo::schema("HEL::online", 22000, 13); urwellADCSchema = hipo::schema("URWELL::adc", 22300, 11); recoilADCSchema = hipo::schema("RECOIL::adc", 22600, 11); + recoiltofTDCSchema = hipo::schema("RECOILTOF::tdc", 22700, 12); rawADCSchema = hipo::schema("RAW::adc", 20000, 11); rawTDCSchema = hipo::schema("RAW::tdc", 20000, 12); rawSCALERSchema = hipo::schema("RAW::scaler", 20000, 13); @@ -146,6 +147,7 @@ HipoSchema::HipoSchema() { urwellADCSchema.parse("sector/B, layer/B, component/S, order/B, ADC/I, time/F, ped/S"); recoilADCSchema.parse( "sector/B, layer/B, component/S, order/B, ADC/I, time/F, ped/S"); + recoiltofTDCSchema.parse("sector/B, row/B, column/B, order/B, TDC/I, ToT/I"); rawADCSchema.parse("crate/B, slot/B, channel/S, order/B, ADC/I, time/F, ped/S"); rawTDCSchema.parse("crate/B, slot/B, channel/S, order/B, TDC/I"); @@ -201,6 +203,7 @@ HipoSchema::HipoSchema() { schemasToLoad["RASTER::adc"] = rasterADCSchema; schemasToLoad["URWELL::adc"] = urwellADCSchema; schemasToLoad["RECOIL::adc"] = recoilADCSchema; + schemasToLoad["RECOILTOF::tdc"] = recoiltofTDCSchema; cout << " Done defining Hipo4 schemas." << endl; @@ -228,6 +231,8 @@ hipo::schema HipoSchema::getSchema(string schemaName, int type) { return fthodoADCSchema; } else if (schemaName == "ft_trk") { return ftrkTDCSchema; + } else if (schemaName == "recoil_tof") { + return recoiltofTDCSchema; } else { if (non_registered_detectors(schemaName, type)) { cout << " SCHEMA " << schemaName << " " << " not found for type " << type << " = " << schemaType << endl; @@ -241,7 +246,7 @@ bool HipoSchema::non_registered_detectors(string schemaName, int type) { if (type == 0) { // non adc detectors: - if (schemaName == "dc" || schemaName == "rich" || schemaName == "atof" || schemaName == "ahdc") { + if (schemaName == "dc" || schemaName == "rich" || schemaName == "atof" || schemaName == "ahdc" || schemaName == "recoil_tof") { return false; } } else if (type == 1) { // non tdc detectors @@ -252,7 +257,7 @@ bool HipoSchema::non_registered_detectors(string schemaName, int type) { if (schemaName == "atof" || schemaName == "band" || schemaName == "bmt" || schemaName == "fmt" || schemaName == "ftm" || schemaName == "dc" || schemaName == "bst" || schemaName == "cnd" || schemaName == "ctof" || schemaName == "ecal" || schemaName == "ftof" || schemaName == "ft_cal" || schemaName == "ft_hodo" || schemaName == "ft_trk" - || schemaName == "htcc" || schemaName == "ltcc" || schemaName == "rich" || schemaName == "rtpc" || schemaName == "urwell" || schemaName == "recoil" || schemaName == "flux") { + || schemaName == "htcc" || schemaName == "ltcc" || schemaName == "rich" || schemaName == "rtpc" || schemaName == "urwell" || schemaName == "recoil" || schemaName == "recoil_tof" || schemaName == "flux") { return false; } } diff --git a/source/output/hipoSchemas.h b/source/output/hipoSchemas.h index 7c7b7190..6a941847 100644 --- a/source/output/hipoSchemas.h +++ b/source/output/hipoSchemas.h @@ -56,6 +56,7 @@ class HipoSchema { hipo::schema helONLINESchema; hipo::schema urwellADCSchema; hipo::schema recoilADCSchema; + hipo::schema recoiltofTDCSchema; hipo::schema rasterADCSchema; hipo::schema rawADCSchema; hipo::schema rawTDCSchema; diff --git a/source/output/hipo_output.cc b/source/output/hipo_output.cc index cca913e0..2e4e1c15 100644 --- a/source/output/hipo_output.cc +++ b/source/output/hipo_output.cc @@ -741,7 +741,7 @@ void hipo_output::writeG4DgtIntegrated(outputContainer *output, vector