diff --git a/CMakeLists.txt b/CMakeLists.txt index 30304f00..79b5d229 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,8 +1,8 @@ # The name of our project is "ShipRoot". CMakeLists files in this project can -# refer to the root source directory of the project as ${ShipRoot_SOURCE_DIR} -# or as ${CMAKE_SOURCE_DIR} and to the root binary directory of the project as +# refer to the root source directory of the project as ${ShipRoot_SOURCE_DIR} +# or as ${CMAKE_SOURCE_DIR} and to the root binary directory of the project as # ${ShipRoot_BINARY_DIR} or ${CMAKE_BINARY_DIR}. -# This difference is important for the base classes which are in FAIRROOT +# This difference is important for the base classes which are in FAIRROOT # and ShipRoot. # Check if cmake has the required version @@ -18,7 +18,7 @@ foreach(p endforeach() # Set name of our project to "ShipRoot". Has to be done -# after check of cmake version since this is a new feature +# after check of cmake version since this is a new feature project(ShipRoot) FIND_PATH(FAIRBASE NAMES FairRun.h PATHS @@ -38,7 +38,7 @@ Else (FAIRBASE) SET(FAIRROOTPATH $ENV{FAIRROOTPATH}) EndIf (FAIRBASE) -# where to look first for cmake modules, before ${CMAKE_ROOT}/Modules/ +# where to look first for cmake modules, before ${CMAKE_ROOT}/Modules/ # is checked set(CMAKE_MODULE_PATH "${CMAKE_SOURCE_DIR}/cmake/modules") set(CMAKE_MODULE_PATH "${FAIRROOTPATH}/share/fairbase/cmake/modules" ${CMAKE_MODULE_PATH}) @@ -106,14 +106,14 @@ CHECK_OUT_OF_SOURCE_BUILD() # message IF(NOT UNIX) MESSAGE(FATAL_ERROR "You're not on an UNIX system. The project was up to now only tested on UNIX systems, so we break here. If you want to go on please edit the CMakeLists.txt in the source directory.") -ENDIF(NOT UNIX) +ENDIF(NOT UNIX) # Check if the external packages are installed into a separate install # directory CHECK_EXTERNAL_PACKAGE_INSTALL_DIR() -# Set the build type. Possibilities are None, Debug, Release, -# RelWithDebInfo and MinSizeRel +# Set the build type. Possibilities are None, Debug, Release, +# RelWithDebInfo and MinSizeRel #SET(CMAKE_BUILD_TYPE Debug) # searches for needed packages @@ -124,6 +124,7 @@ CHECK_EXTERNAL_PACKAGE_INSTALL_DIR() find_package(ROOT 6.10.06 REQUIRED) find_package2(PUBLIC Pythia8 REQUIRED) find_package2(PUBLIC EvtGen REQUIRED) + # (GENERATORS REQUIRED) find_package2(PUBLIC GEANT3) find_package2(PUBLIC GEANT4) @@ -157,8 +158,8 @@ ENDIF(DEFINED $ENV{GSL_ROOT}) Message("-- Looking for Boost ...") # If an older version of boost is found both of the variables below are - # cached and in a second cmake run, a good boost version is found even - # if the version is to old. + # cached and in a second cmake run, a good boost version is found even + # if the version is to old. # To overcome this problem both variables are cleared before checking # for boost. Unset(Boost_INCLUDE_DIR CACHE) @@ -186,8 +187,8 @@ add_compile_definitions(PYTHIA8_V=${ver}) # compiler will not generate any warnings. This is usefull since one is only # interested about warnings from the own project. SYSTEM_INCLUDE_DIRECTORIES # is defined in FairMacros.cmake. In the moment the defined directories are -# the root and boost include directories. -Set(SYSTEM_INCLUDE_DIRECTORIES +# the root and boost include directories. +Set(SYSTEM_INCLUDE_DIRECTORIES ${SYSTEM_INCLUDE_DIRECTORIES} ${BASE_INCLUDE_DIRECTORIES} ${FAIRLOGGER_INCLUDE_DIR} @@ -223,7 +224,7 @@ SET(LD_LIBRARY_PATH ${CBMLIBDIR} ${LD_LIBRARY_PATH}) # Message("Set default install path ...") #ENDIF(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT) -install(DIRECTORY geometry DESTINATION pnd_install +install(DIRECTORY geometry DESTINATION pnd_install PATTERN ".svn" EXCLUDE ) diff --git a/shipLHC/AdvTargetHit.cxx b/shipLHC/AdvTargetHit.cxx index 15686ab5..dc4a6cde 100644 --- a/shipLHC/AdvTargetHit.cxx +++ b/shipLHC/AdvTargetHit.cxx @@ -1,10 +1,24 @@ #include "AdvTargetHit.h" - +#include "AdvTargetPoint.h" +#include "digitisation/AdvSignal.h" +#include "digitisation/AdvDigitisation.h" #include "FairLogger.h" #include "TGeoBBox.h" #include "TGeoManager.h" #include "TGeoNavigator.h" #include "TROOT.h" +#include "TRandom.h" +#include "TVector3.h" +#include "TStopwatch.h" + +#include +#include +#include +#include +#include +#include +using std::cout; +using std::endl; // ----- Default constructor ------------------------------------------- AdvTargetHit::AdvTargetHit() @@ -29,7 +43,10 @@ AdvTargetHit::AdvTargetHit(Int_t detID) AdvTargetHit::AdvTargetHit(Int_t detID, const std::vector& V) : SndlhcHit(detID) { + AdvDigitisation advdigi{}; + fDigitisedHit = advdigi.digirunoutput(detID, V); flag = true; + for (Int_t i = 0; i < 16; i++) { fMasked[i] = kFALSE; } diff --git a/shipLHC/AdvTargetHit.h b/shipLHC/AdvTargetHit.h index 39901a86..c3d01906 100644 --- a/shipLHC/AdvTargetHit.h +++ b/shipLHC/AdvTargetHit.h @@ -3,9 +3,15 @@ #include "SiSensor.h" #include "SndlhcHit.h" +#include "TArrayD.h" +#include "TVector3.h" +#include "digitisation/AdvSignal.h" +#include class AdvTargetPoint; +class TArrayD; + class AdvTargetHit : public SndlhcHit { public: @@ -24,6 +30,7 @@ class AdvTargetHit : public SndlhcHit bool isValid() const { return flag; } bool isMasked(Int_t i) const { return fMasked[i]; } void SetMasked(Int_t i) { fMasked[i] = kTRUE; } + std::map> GetHit() { return fDigitisedHit; } int constexpr GetStation() { return fDetectorID >> 17; } int constexpr GetPlane() { return (fDetectorID >> 16) % 2; } // 0 is X-plane, 1 is Y-pane int constexpr GetRow() { return (fDetectorID >> 13) % 8; } @@ -32,11 +39,12 @@ class AdvTargetHit : public SndlhcHit int constexpr GetStrip() { return (fDetectorID) % 1024; } int constexpr GetModule() { return advsnd::target::columns * GetRow() + 1 + GetColumn(); } bool constexpr isVertical() { return GetPlane() == 0; }; + private: bool flag; ///< flag bool fMasked[16]; /// masked signal - + std::map> fDigitisedHit; ClassDef(AdvTargetHit, 1); }; diff --git a/shipLHC/CMakeLists.txt b/shipLHC/CMakeLists.txt index fa3e0c4e..3c7912a4 100644 --- a/shipLHC/CMakeLists.txt +++ b/shipLHC/CMakeLists.txt @@ -6,6 +6,8 @@ set(INCLUDE_DIRECTORIES ${BASE_INCLUDE_DIRECTORIES} ${CMAKE_SOURCE_DIR}/shipdata ${CMAKE_SOURCE_DIR}/shipLHC +${CMAKE_SOURCE_DIR}/shipLHC/digitisation +${CMAKE_SOURCE_DIR}/shipLHC/external/clhep/include ${CMAKE_SOURCE_DIR}/veto ${CMAKE_SOURCE_DIR}/genfit/core/include ${CMAKE_SOURCE_DIR}/genfit/trackReps/include @@ -23,7 +25,7 @@ ${ROOT_LIBRARY_DIR} ${FAIRROOT_LIBRARY_DIR} ) - + link_directories( ${LINK_DIRECTORIES}) set(SRCS @@ -49,6 +51,13 @@ AdvMuFilter.cxx AdvMuFilterPoint.cxx SNDLHCEventHeader.cxx sndRecoTrack.cxx +digitisation/ChargeDivision.cxx +digitisation/ChargeDrift.cxx +digitisation/InducedCharge.cxx +digitisation/FrontendDriver.cxx +digitisation/StripNoise.cxx +digitisation/SiG4UniversalFluctuation.cxx +digitisation/AdvDigitisation.cxx ) Set(HEADERS ) diff --git a/shipLHC/digitisation/AdvDigitisation.cxx b/shipLHC/digitisation/AdvDigitisation.cxx new file mode 100644 index 00000000..718d6f54 --- /dev/null +++ b/shipLHC/digitisation/AdvDigitisation.cxx @@ -0,0 +1,66 @@ +#include "AdvDigitisation.h" +#include "AdvTargetPoint.h" +#include "ChargeDivision.h" +#include "ChargeDrift.h" +#include "InducedCharge.h" +#include "FrontendDriver.h" +#include "StripNoise.h" +#include "EnergyFluctUnit.h" +#include "SurfaceSignal.h" +#include "AdvSignal.h" +#include "SiDigiParameters.h" + +#include "TFile.h" +#include "TTree.h" +#include "TVector3.h" +#include "TStopwatch.h" + +#include +#include +#include +#include +#include +#include +using namespace std; + +// Running the digitisation + +/* To be included : + Saturation of FED dynamic range depending on mode + Add noise + Add FED modes */ + +AdvDigitisation::AdvDigitisation() {} + +std::map> AdvDigitisation::digirunoutput(Int_t detID, const std::vector &V) +{ + // Charge Division + ChargeDivision chargedivision{}; + std::vector EnergyLossVector = chargedivision.Divide(detID, V); + + //Charge Drift + ChargeDrift chargedrift{}; + std::vector DiffusionSignal = chargedrift.Drift(EnergyLossVector); + + //Induced Charge on strips + InducedCharge inducedcharge{}; + AdvSignal ResponseSignal = inducedcharge.IntegrateCharge(DiffusionSignal); + + //FED Response + FrontendDriver frontenddriver{}; + AdvSignal FEDResponseSignal = frontenddriver.FEDResponse(ResponseSignal); + + //Creating map of hit + std::map> DigitisedHit; + std::vector Charge = FEDResponseSignal.getIntegratedSignal(); + std::vector Strips = FEDResponseSignal.getStrips(); + std::vector ADC(Charge.size()); + std::transform(Charge.begin(), Charge.end(), ADC.begin(), [](Double_t x) { + if (x>0){return (int)x;} + else {return 0;} + }); + DigitisedHit["Strips"] = Strips; + DigitisedHit["ADC"] = ADC; + + return DigitisedHit; +} \ No newline at end of file diff --git a/shipLHC/digitisation/AdvDigitisation.h b/shipLHC/digitisation/AdvDigitisation.h new file mode 100644 index 00000000..69c334f6 --- /dev/null +++ b/shipLHC/digitisation/AdvDigitisation.h @@ -0,0 +1,52 @@ +#ifndef SHIPLHC_ADVDIGITISATION_H_ +#define SHIPLHC_ADVDIGITISATION_H_ + +#include "AdvTargetPoint.h" +#include "AdvSignal.h" +#include "TVector3.h" +#include "TGeoNavigator.h" + +#include "FairGeoBuilder.h" +#include "FairGeoInterface.h" +#include "FairGeoLoader.h" +#include "FairGeoMedia.h" +#include "FairGeoMedium.h" +#include "FairGeoNode.h" +#include "FairGeoTransform.h" +#include "FairGeoVolume.h" +#include "FairRootManager.h" +#include "FairVolume.h" +#include "ShipDetectorList.h" +#include "ShipStack.h" +#include "ShipUnit.h" +#include "SiSensor.h" +#include "TGeoArb8.h" +#include "TGeoBBox.h" +#include "TGeoCompositeShape.h" +#include "TGeoGlobalMagField.h" +#include "TGeoManager.h" +#include "TGeoMaterial.h" +#include "TGeoMedium.h" +#include "TGeoSphere.h" +#include "TGeoTrd1.h" +#include "TGeoTrd2.h" +#include "TGeoTube.h" +#include "TGeoUniformMagField.h" +#include "TParticle.h" +#include "TString.h" // for TString +#include "TVector3.h" +#include "TVirtualMC.h" + +#include +#include +using namespace std; + + + +class AdvDigitisation +{ + public: + AdvDigitisation(); + std::map> digirunoutput(Int_t detID, const std::vector& V); +}; +#endif diff --git a/shipLHC/digitisation/AdvSignal.h b/shipLHC/digitisation/AdvSignal.h new file mode 100644 index 00000000..166de92d --- /dev/null +++ b/shipLHC/digitisation/AdvSignal.h @@ -0,0 +1,33 @@ +#ifndef SHIPLHC_ADVSIGNAL_H_ +#define SHIPLHC_ADVSIGNAL_H_ + +#include "TVector3.h" + +#include + +class AdvSignal +{ + public: + AdvSignal() + : Strips_(), IntegratedSignal_() + { + } + + AdvSignal(std::vector Strips, std::vector IntegratedSignal) + : Strips_(Strips) + , IntegratedSignal_(IntegratedSignal) + { + } + + std::vector getStrips() const { return Strips_; } + std::vector getIntegratedSignal() const { return IntegratedSignal_; } + + void setStrips(std::vector value) { Strips_ = value; } + void setIntegratedSignal(std::vector value) { IntegratedSignal_ = value; } + + private: + std::vector Strips_; + std::vector IntegratedSignal_; +}; + +#endif // SHIPLHC_ADVSIGNAL_H_ diff --git a/shipLHC/digitisation/ChargeDivision.cxx b/shipLHC/digitisation/ChargeDivision.cxx new file mode 100644 index 00000000..4400b2c3 --- /dev/null +++ b/shipLHC/digitisation/ChargeDivision.cxx @@ -0,0 +1,200 @@ +#include "ChargeDivision.h" +#include "AdvTargetPoint.h" +#include "EnergyFluctUnit.h" +#include "SiG4UniversalFluctuation.h" +#include "TVector3.h" +#include "TGeoNavigator.h" +#include "SiDigiParameters.h" + +#include "FairGeoBuilder.h" +#include "FairGeoInterface.h" +#include "FairGeoLoader.h" +#include "FairGeoMedia.h" +#include "FairGeoMedium.h" +#include "FairGeoNode.h" +#include "FairGeoTransform.h" +#include "FairGeoVolume.h" +#include "FairRootManager.h" +#include "FairVolume.h" +#include "ShipDetectorList.h" +#include "ShipStack.h" +#include "ShipUnit.h" +#include "SiSensor.h" +#include "TGeoArb8.h" +#include "TGeoBBox.h" +#include "TGeoCompositeShape.h" +#include "TGeoGlobalMagField.h" +#include "TGeoManager.h" +#include "TGeoMaterial.h" +#include "TGeoMedium.h" +#include "TGeoSphere.h" +#include "TGeoTrd1.h" +#include "TGeoTrd2.h" +#include "TGeoTube.h" +#include "TGeoUniformMagField.h" +#include "TParticle.h" +#include "TString.h" // for TString +#include "TVector3.h" +#include "TVirtualMC.h" + + +#include +#include +#include +#include +#include +#include +#include +#include + +/*Class for division of energy deposited along particle track in each module. + +Not included : - Possibility of Lorentz angle for shift in track due to magnetic field i + + */ + +ChargeDivision::ChargeDivision() {} + +TVector3 ChargeDivision::getLocal(Int_t detID, TVector3 point) +{ + + /*Global to Local coordinate conversion*/ + + TVector3 local_point; + // Calculate the detector id as per the geofile, where strips are disrespected + // int strip = (detID) % 1024; // actual strip ID + int geofile_detID = detID; // the det id number needed to read the geometry + int station = geofile_detID >> 17; + int plane = (geofile_detID >> 16) % 2; + int row = (geofile_detID >> 13) % 8; + int column = (geofile_detID >> 11) % 4; + int sensor = geofile_detID; + int sensor_module = advsnd::target::columns * row + 1 + column; + + TString path = TString::Format("/cave_1/" + "Detector_0/" + "volAdvTarget_1/" + "TrackingStation_%d/" + "TrackerPlane_%d/" + "SensorModule_%d/" + "SensorVolumeTarget_%d", + station, + plane, + sensor_module, + sensor); + TGeoNavigator *nav = gGeoManager->GetCurrentNavigator(); + if (nav->CheckPath(path)) { + nav->cd(path); + } else { + LOG(FATAL) << path; + } + // Get the corresponding node, which is a sensor made of strips + TGeoNode *W = nav->GetCurrentNode(); + Double_t x = point.X(); + Double_t y = point.Y(); + Double_t z = point.Z(); + double global_pos[3] = {x, y, z}; + double local_pos[3]; + nav->MasterToLocal(global_pos, local_pos); + + local_point = local_pos; + return local_point; +} + + + +TVector3 ChargeDivision::DriftDir(TVector3 EntryPoint, TVector3 ExitPoint, float length) +{ + /*Getting the position of the segment by using the length and direction of the particle track*/ + + TVector3 DriftDirUnit = (EntryPoint - ExitPoint).Unit(); + TVector3 DriftMul = DriftDirUnit * length; + TVector3 DriftPos = EntryPoint - (DriftDirUnit * length); + return DriftPos; +} + +std::vector ChargeDivision::Divide(Int_t detID, const std::vector& V) +{ + std::vector ELossVector; + + for (int i = 0; i < V.size(); i++) { + + std::vector fluctEnergy; + std::vector driftPos; + std::vector glob_driftPos; + Int_t pdgcode = V[i]->PdgCode(); + + //Getting the mass and charge of particle, otherwise assigning pion mass and charge + + if (TDatabasePDG::Instance()->GetParticle(pdgcode)) { + ParticleMass = (TDatabasePDG::Instance()->GetParticle(V[i]->PdgCode())->Mass()) * 1000; // in MeV + ParticleCharge = TDatabasePDG::Instance()->GetParticle(V[i]->PdgCode())->Charge(); + } else { + std::cout << "Could not find particle " << pdgcode << " , assuming pion mass and charge." << std::endl; + ParticleMass = 139.57; + ParticleCharge = 1; + + }; + + //Conversion of global entry and exit point to local coordinates. + + TVector3 local_entry_point = getLocal(V[i] -> GetDetectorID(), V[i]->GetEntryPoint()); + TVector3 local_exit_point = getLocal(V[i] -> GetDetectorID(), V[i]->GetExitPoint()); + + //Calculating the number of segments in the track by dividing the tracklength by the number of divisions per strip. To get an approximate of the number of strips, divide delta x by the strip pitch. Number of segments is 1 if the particle is neutral or very low mass. + + double len = (local_entry_point - local_exit_point).Mag(); + + if (fabs(ParticleMass) < 1e-6 || ParticleCharge == 0) { + NumberofSegments = 1; + } else { + NumberofSegments = 1 + + (stripsensor::chargedivision::ChargeDivisionsperStrip + * abs((local_entry_point.X() - local_exit_point.X()) + / stripsensor::chargedivision::StripPitch)); + } + + segLen = (len / NumberofSegments) * 10; // in mm + + // Getting the energy fluctuations per segment along with the local and global segment position. + + SiG4UniversalFluctuation sig4fluct{}; + + Double_t Etotal = V[i]->GetEnergyLoss() * 1000; // in MeV + Double_t Emean = Etotal / NumberofSegments; + Double_t momentum = sqrt(pow(V[i]->GetPx(), 2) + pow(V[i]->GetPy(), 2) + pow(V[i]->GetPz(), 2)) * 1000; + + if (NumberofSegments > 1) { + for (Int_t j = 0; j < NumberofSegments; j++) { + fluctEnergy.push_back( + sig4fluct.SampleFluctuations(ParticleMass, ParticleCharge, Emean, momentum, segLen)); + driftPos.push_back(DriftDir(local_entry_point, local_exit_point, (segLen * j) / 10)); + glob_driftPos.push_back(DriftDir(V[i]->GetEntryPoint(), V[i]->GetExitPoint(), (segLen * j) / 10)); + } + } else { + fluctEnergy.push_back(V[i]->GetEnergyLoss() * 1000); + driftPos.push_back(DriftDir(local_entry_point, local_exit_point, (segLen) / 10)); + glob_driftPos.push_back(DriftDir(V[i]->GetEntryPoint(), V[i]->GetExitPoint(), (segLen) / 10)); + } + + // Scaling each fluctuation accordingly with the total energy loss + + double sume = 0; + for (int n = 0; n < size(fluctEnergy); n++) { + sume = sume + fluctEnergy[n]; + } + if (sume > 0.) { + float rescale_ratio = Etotal / sume; + for (int m = 0; m < size(fluctEnergy); m++) { + fluctEnergy[m] = (fluctEnergy[m] * rescale_ratio) / 1000; + } + } + EnergyFluctUnit EnergyFluctuations(fluctEnergy, segLen / 10, driftPos, glob_driftPos); + ELossVector.push_back(EnergyFluctuations); + } + + // Returns the vector with segments and their energy depositon (in GeV) and position (in cm) + + return ELossVector; + +} diff --git a/shipLHC/digitisation/ChargeDivision.h b/shipLHC/digitisation/ChargeDivision.h new file mode 100644 index 00000000..eefda0df --- /dev/null +++ b/shipLHC/digitisation/ChargeDivision.h @@ -0,0 +1,27 @@ +#ifndef SHIPLHC_CHARGEDIVISION_H_ +#define SHIPLHC_CHARGEDIVISION_H_ + +#include "AdvTargetPoint.h" +#include "EnergyFluctUnit.h" +#include "TVector3.h" + +#include +#include + +class ChargeDivision +{ + public: + ChargeDivision(); + std::vector Divide(Int_t detID, const std::vector& V); + TVector3 DriftDir(TVector3 EntryPoint, TVector3 ExitPoint, float length); + TVector3 getLocal(Int_t detID,TVector3 point); + + private: + std::vector PulseValues; + Double_t ParticleCharge; + Double_t ParticleMass; + Int_t NumberofSegments; + float segLen; +}; + +#endif // SHIPLHC_CHARGEDIVISION_H_ diff --git a/shipLHC/digitisation/ChargeDrift.cxx b/shipLHC/digitisation/ChargeDrift.cxx new file mode 100644 index 00000000..4a8f4466 --- /dev/null +++ b/shipLHC/digitisation/ChargeDrift.cxx @@ -0,0 +1,145 @@ +#include "ChargeDrift.h" +#include "SurfaceSignal.h" +#include "TVector3.h" +#include "TGraph.h" +#include "TCanvas.h" +#include "SiDigiParameters.h" + +#include +#include +#include + +// Class for drifting the charge to the surface of the detector and diffusing the charge + +ChargeDrift::ChargeDrift() {} + +std::vector ChargeDrift::Drift(std::vector EnergyLossVector) +{ + std::vector DiffusionSignal; + for (int k = 0; k < EnergyLossVector.size(); k++) + { + std::vector diffusionarea; + std::vector diffusionpos; + std::vector amplitude; + + DriftPos = EnergyLossVector[k].getDriftPos(); + EnergyFluct = EnergyLossVector[k].getEfluct(); + + for (int i = 0; i < DriftPos.size(); i ++) + { + + //Using CMS approximation for calculating the drift time + + DriftDistance = 0.025 - DriftPos[i].Z(); + DriftDistanceFraction = DriftDistance / stripsensor::drift::module_thickness; + DriftDistanceFraction = DriftDistanceFraction > 0. ? DriftDistanceFraction : 0. ; + DriftDistanceFraction = DriftDistanceFraction < 1. ? DriftDistanceFraction : 1. ; + + //DriftTime = GetDriftTime(DriftDistance); + + Double_t tn = (stripsensor::drift::module_thickness * stripsensor::drift::module_thickness) / (2 * stripsensor::drift::depletion_voltage * stripsensor::drift::charge_mobility); + DriftTime = -tn * log(1 - 2 * stripsensor::drift::depletion_voltage * DriftDistanceFraction / (stripsensor::drift::depletion_voltage + stripsensor::drift::applied_voltage)) + stripsensor::drift::chargedistributionRMS; + + //Position of drifted charge on the surface of the detector + + DriftPositiononSurface.SetXYZ(DriftPos[i].X(), DriftPos[i].Y(), DriftPos[i].Z()+DriftDistance); + + // Calculating the diffusion area of the signal on the surface + + DiffusionConstant = (1.38E-23 / 1.6E-19) * stripsensor::drift::charge_mobility * stripsensor::drift::temperature; + + DiffusionArea = sqrt(2* DiffusionConstant * DriftTime); + + // Calculating the number of electrons produced per energy deposit + + Amplitude = EnergyFluct[i]> 0. ? floor(EnergyFluct[i]*1e9 / stripsensor::drift::perGeV) : 0. ; + + diffusionarea.push_back(DiffusionArea); + diffusionpos.push_back(DriftPositiononSurface); + amplitude.push_back(Amplitude); + + } + // Returning the surface signal with the area of the diffusion, position at surface and the total number of electrons + + SurfaceSignal Diffused(diffusionarea, diffusionpos, amplitude); + + DiffusionSignal.push_back(Diffused); + } + return DiffusionSignal; +} + +Double_t ChargeDrift::GetDriftTime(Double_t distance) +{ + // Calculate drift time using textbook formulas for pn junction + + Double_t z_width = 400e-6; // in um + Double_t p = 120e-6; + Double_t w = 30e-6; + Double_t h = 5e-6; + + Double_t Na = 1e15; // in cm + Double_t Nd = 1e12; // in cm + + Double_t e = 1.6e-19; // in C + + Double_t k = 1.38e-23; + Double_t T = 300; + Double_t epsilon = 8.854e-14*11.9; + + Double_t ni = 1.45e10; // in cm + + Double_t V_a = 470; + Double_t V_0 = (((k*T)/e)*log((Na*Nd)/(pow(ni, 2)))) + V_a; + + Double_t W = sqrt((2*epsilon*(Nd + Na)*V_0)/(e*Nd*Na)); + + Double_t Wp = (Nd * W) / (Nd + Na); + Double_t Wn = (Na * W) / (Nd + Na); + + Int_t num = 80; + + Double_t x[num] ; + Double_t y[num] ; + Double_t E_strip[num] ; + + Double_t x_start = 0 ; // in cm + Double_t x_end = 500e-4; + + + for (int j = 0; j < num ; j ++ ) + { + x[j] = x_start + abs(j*((x_start - x_end)/num)); + + } + + for (int i = 0; i < num; i++) + { + if ((-Wp <= x[i]) && (x[i] < 0)) + { + E_strip[i] = -((e*Na*x[i])/epsilon) - ((e*Na*Wp)/epsilon); + + } + else if ((0 <= x[i]) && (x[i] < Wn)) + { + E_strip[i] = ((e*Nd*x[i])/epsilon) - ((e*Nd*Wn)/epsilon); + } + else + { + E_strip[i] = 0; + } + } + + Double_t dt[num]; + Double_t tfrac[num]; + + for(int l = 0; l < num ; l++) + { + //tfrac[l] = -x[l]/250e-6; + dt[l] = -x[l]/(480*E_strip[l]); + } + + auto gr = new TGraph (num, x, dt); + gr->Draw(); + Double_t drifttime = gr->Eval(distance); + return drifttime; +} \ No newline at end of file diff --git a/shipLHC/digitisation/ChargeDrift.h b/shipLHC/digitisation/ChargeDrift.h new file mode 100644 index 00000000..ab616a44 --- /dev/null +++ b/shipLHC/digitisation/ChargeDrift.h @@ -0,0 +1,30 @@ +#ifndef SHIPLHC_CHARGEDRIFT_H_ +#define SHIPLHC_CHARGEDRIFT_H_ + +#include "EnergyFluctUnit.h" +#include "SurfaceSignal.h" + +#include +#include + +class ChargeDrift +{ + public: + ChargeDrift(); + std::vector Drift(std::vector EnergyLossVector); + Double_t GetDriftTime(Double_t distance); + + protected: + std::vector DriftPos; + std::vector EnergyFluct; + Double_t DriftDistance; + Double_t DriftDistanceFraction; + Double_t DriftTime; + TVector3 DriftPositiononSurface; + Double_t DiffusionArea; + Double_t DiffusionConstant; + Double_t Amplitude; + +}; + +#endif // SHIPLHC_CHARGEDRIFT_H_ \ No newline at end of file diff --git a/shipLHC/digitisation/EnergyFluctUnit.h b/shipLHC/digitisation/EnergyFluctUnit.h new file mode 100644 index 00000000..73397a66 --- /dev/null +++ b/shipLHC/digitisation/EnergyFluctUnit.h @@ -0,0 +1,37 @@ +#ifndef SHIPLHC_ENERGYFLUCTUNIT_H_ +#define SHIPLHC_ENERGYFLUCTUNIT_H_ + +#include "TVector3.h" + +#include + +class EnergyFluctUnit +{ + public: + EnergyFluctUnit(std::vector Efluct, float segLen, std::vector DriftPos, std::vector glob_DriftPos) + : Efluct_(Efluct) + , segLen_(segLen) + , DriftPos_(DriftPos) + , glob_DriftPos_(glob_DriftPos) + { + } + + std::vector getEfluct() const { return Efluct_; } + float getsegLen() const { return segLen_; } + std::vector getDriftPos() const { return DriftPos_; } + std::vector getglobDriftPos() const { return glob_DriftPos_; } + + + void setEfluct(std::vector value) { Efluct_ = value; } + void setsegLen(float value) { segLen_ = value; } + void setDriftPos(std::vector value) { DriftPos_ = value; } + void setglobDriftPos(std::vector value) { glob_DriftPos_ = value; } + + private: + std::vector Efluct_; + float segLen_; + std::vector DriftPos_; + std::vector glob_DriftPos_; +}; + +#endif diff --git a/shipLHC/digitisation/FrontendDriver.cxx b/shipLHC/digitisation/FrontendDriver.cxx new file mode 100644 index 00000000..96319620 --- /dev/null +++ b/shipLHC/digitisation/FrontendDriver.cxx @@ -0,0 +1,245 @@ +#include "FrontendDriver.h" +#include "StripNoise.h" +#include "InducedCharge.h" +#include "SiDigiParameters.h" +#include "AdvSignal.h" + +#include +#include +#include +#include +#include +using namespace std; + + +FrontendDriver::FrontendDriver() {} + +AdvSignal FrontendDriver::FEDResponse(AdvSignal Signal) +{ + AdvSignal ADCResponse = ADCConversion(Signal); + + AdvSignal FEDResponseSignal ; + std::vector temp_FEDResponseSignal; + + StripNoise stripnoise{}; + InducedCharge inducedcharge{}; + + if (stripsensor::frontend::ZSModeOption && stripsensor::frontend::NoiseOption) + { + FEDResponseSignal = stripnoise.AddGaussianTailNoise(Signal); + temp_FEDResponseSignal.push_back(FEDResponseSignal); + FEDResponseSignal = inducedcharge.Combine(temp_FEDResponseSignal); + FEDResponseSignal = ZeroSuppressionAlgorithms(FEDResponseSignal); + } + if (!stripsensor::frontend::ZSModeOption) + { + FEDResponseSignal = stripnoise.AddGaussianNoise(Signal); + FEDResponseSignal = stripnoise.AddCMNoise(FEDResponseSignal); + temp_FEDResponseSignal.push_back(FEDResponseSignal); + FEDResponseSignal = inducedcharge.Combine(temp_FEDResponseSignal); + } + FEDResponseSignal = SaturateRange(FEDResponseSignal); + return FEDResponseSignal; +} + +AdvSignal FrontendDriver::ADCConversion(AdvSignal ResponseSignal) +{ + std::vector NumberofElectrons = ResponseSignal.getIntegratedSignal(); + for (int i = 0; i < NumberofElectrons.size(); i++) + { + ADCcount.push_back(stripsensor::frontend::StripNoise + std::ceil(NumberofElectrons[i]/stripsensor::frontend::ElectronperADC)); + } + + AdvSignal ADCResponse(ResponseSignal.getStrips(), ADCcount); + + return ADCResponse; +} + +AdvSignal FrontendDriver::SaturateRange(AdvSignal Signal) +{ + std::vector Charge = Signal.getIntegratedSignal(); + for (int i = 0; i < Charge.size(); i++) + { + if (stripsensor::frontend::ZSModeOption) + { + if (Charge[i] > 1022) + { + Charge[i] = 255; + } + if (Charge[i] > 253) + { + Charge[i] = 254; + } + if (Charge[i] < 0) + { + Charge[i] = 0; + } + + } else { + if (Charge[i] > 1023) + { + Charge[i] = 1023; + } + if (Charge[i] < 0) + { + Charge[i] = 0; + } + } + } + AdvSignal SaturatedSignal(Signal.getStrips(), Charge); + return SaturatedSignal; +} + +AdvSignal FrontendDriver::ZeroSuppressionAlgorithms(AdvSignal Signal) +{ + std::vector Amplitude = Signal.getIntegratedSignal(); + std::vector Strips = Signal.getStrips(); + + std::vector ClusterStrips; + std::vector ClusterAmplitudes; + + std::vector temp_ClusterStrips; + std::vector temp_ClusterAmplitudes; + + Int_t mode = stripsensor::frontend::ZeroSuppressionMode; + + switch (mode) { + case 1: + for (int i = 0; i < Amplitude.size(); i++) + { + if (Amplitude[i] > stripsensor::frontend::ZeroSuppressionMode1T*stripsensor::frontend::StripNoise) + { + ClusterStrips.push_back(Strips[i]); + ClusterAmplitudes.push_back(Amplitude[i]); + } + } + break; + case 2: + for (int i = 0; i < Amplitude.size(); i++) + { + if (Amplitude[i] > 2*stripsensor::frontend::StripNoise) + { + temp_ClusterStrips.push_back(Strips[i]); + temp_ClusterAmplitudes.push_back(Amplitude[i]); + } + } + + for (int j = 0; j < temp_ClusterStrips.size(); j++) + { + if((std::find(temp_ClusterStrips.begin(), temp_ClusterStrips.end(), temp_ClusterStrips[j]+1)==temp_ClusterStrips.end()) && (std::find(temp_ClusterStrips.begin(), temp_ClusterStrips.end(), temp_ClusterStrips[j]-1)==temp_ClusterStrips.end())) + { + if (temp_ClusterAmplitudes[j] > 5*stripsensor::frontend::StripNoise) + { + ClusterStrips.push_back(temp_ClusterStrips[j]); ClusterAmplitudes.push_back(temp_ClusterAmplitudes[j]); + } + } else { + ClusterStrips.push_back(temp_ClusterStrips[j]); ClusterAmplitudes.push_back(temp_ClusterAmplitudes[j]); + } + } + break; + case 3: + for (int i = 0; i < Amplitude.size(); i++) + { + if (Amplitude[i] > 3*stripsensor::frontend::StripNoise) + { + temp_ClusterStrips.push_back(Strips[i]); + temp_ClusterAmplitudes.push_back(Amplitude[i]); + } + } + + ClusterStrips = temp_ClusterStrips; + ClusterAmplitudes = temp_ClusterAmplitudes; + + for (int j = 0; j < temp_ClusterStrips.size(); j++) + { + if((std::find(temp_ClusterStrips.begin(), temp_ClusterStrips.end(), temp_ClusterStrips[j]+1)!=temp_ClusterStrips.end()) && (std::find(ClusterStrips.begin(), ClusterStrips.end(), temp_ClusterStrips[j]+1)==ClusterStrips.end())) + { + auto it = std::find(temp_ClusterStrips.begin(), temp_ClusterStrips.end(), temp_ClusterStrips[j]+1 ); + if ((temp_ClusterAmplitudes[it - temp_ClusterStrips.begin()])>0) + { + ClusterStrips.push_back(temp_ClusterStrips[it - temp_ClusterStrips.begin()]); + ClusterAmplitudes.push_back(temp_ClusterAmplitudes[it - temp_ClusterStrips.begin()]); + } + + } + if((std::find(temp_ClusterStrips.begin(), temp_ClusterStrips.end(), temp_ClusterStrips[j]-1)!=temp_ClusterStrips.end()) && (std::find(ClusterStrips.begin(), ClusterStrips.end(), temp_ClusterStrips[j]-1)==ClusterStrips.end())) + { + auto itlower = std::find(temp_ClusterStrips.begin(), temp_ClusterStrips.end(), temp_ClusterStrips[j]-1); + if ((temp_ClusterAmplitudes[itlower - temp_ClusterStrips.begin()])>0) + { + ClusterStrips.push_back(temp_ClusterStrips[itlower - temp_ClusterStrips.begin()]); + ClusterAmplitudes.push_back(temp_ClusterAmplitudes[itlower - temp_ClusterStrips.begin()]); + } + + } + + } + break; + case 4: + Double_t SumCharge = 0; + Double_t SumNoise = 0; + + for (int i = 0; i < Amplitude.size(); i++) + { + if (Amplitude[i] > 3*stripsensor::frontend::StripNoise) + { + temp_ClusterStrips.push_back(Strips[i]); + temp_ClusterAmplitudes.push_back(Amplitude[i]); + } + } + Int_t neighbhour = 0 ; + Int_t neighbhourlower = 0; + Int_t j = 0; + std::vector Clusters ; + while (j < temp_ClusterStrips.size()) + { + neighbhour = 1; + neighbhourlower = 1; + while((std::find(temp_ClusterStrips.begin(), temp_ClusterStrips.end(), temp_ClusterStrips[j]+neighbhour)!=temp_ClusterStrips.end())) + { + if((std::find(ClusterStrips.begin(), ClusterStrips.end(), temp_ClusterStrips[j])==ClusterStrips.end())) + { + ClusterStrips.push_back(temp_ClusterStrips[j]); + ClusterAmplitudes.push_back(temp_ClusterAmplitudes[j]); + } + auto it = std::find(temp_ClusterStrips.begin(), temp_ClusterStrips.end(), temp_ClusterStrips[j]+neighbhour ); + if ((temp_ClusterAmplitudes[it - temp_ClusterStrips.begin()])>2*stripsensor::frontend::StripNoise) + { + ClusterStrips.push_back(temp_ClusterStrips[it - temp_ClusterStrips.begin()]); + ClusterAmplitudes.push_back(temp_ClusterAmplitudes[it - temp_ClusterStrips.begin()]); + neighbhour += 1; + } + + } + if (neighbhour > 1) + { + Clusters.push_back(j); + Clusters.push_back(j+neighbhour); + } + j = j + neighbhour; + } + + for (int l = 0; l < Clusters.size(); l+=2) + { + for (int m = Clusters[l]; m < Clusters[l+1]; m++) + { + SumCharge += temp_ClusterAmplitudes[m]; + SumNoise += pow(stripsensor::frontend::StripNoise, 2); + } + if (SumCharge < 5*SumNoise) + { + + for (int m = Clusters[l]; m < Clusters[l+1]; m++) + { + ClusterStrips.erase(find(ClusterStrips.begin(), ClusterStrips.end(), temp_ClusterStrips[m])); + ClusterAmplitudes.erase(find(ClusterAmplitudes.begin(), ClusterAmplitudes.end(), temp_ClusterAmplitudes[m])); + } + } + SumCharge = 0 ; + SumNoise = 0 ; + } + break; + } + AdvSignal ClusterSignal(ClusterStrips, ClusterAmplitudes); + return ClusterSignal; +} \ No newline at end of file diff --git a/shipLHC/digitisation/FrontendDriver.h b/shipLHC/digitisation/FrontendDriver.h new file mode 100644 index 00000000..9cba67f8 --- /dev/null +++ b/shipLHC/digitisation/FrontendDriver.h @@ -0,0 +1,26 @@ +#ifndef SHIPLHC_FRONTENDDRIVER_H_ +#define SHIPLHC_FRONTENDDRIVER_H_ + +#include "AdvSignal.h" +#include "StripNoise.h" + +#include +#include + +class FrontendDriver +{ + public: + FrontendDriver(); + AdvSignal FEDResponse(AdvSignal Signal); + AdvSignal ADCConversion(AdvSignal ResponseSignal); + AdvSignal SaturateRange(AdvSignal Signal); + AdvSignal ZeroSuppressionAlgorithms(AdvSignal Signal); + + void TestingAlgorithm(); + + private: + std::vector ADCcount; + +}; + +#endif \ No newline at end of file diff --git a/shipLHC/digitisation/InducedCharge.cxx b/shipLHC/digitisation/InducedCharge.cxx new file mode 100644 index 00000000..6e942bdd --- /dev/null +++ b/shipLHC/digitisation/InducedCharge.cxx @@ -0,0 +1,269 @@ +#include "InducedCharge.h" +#include "SurfaceSignal.h" +#include "TVector3.h" +#include "TGraph.h" +#include "TRandom.h" +#include "SiSensor.h" +#include "AdvSignal.h" +#include "SiDigiParameters.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +// Class for calculating the induced charge on the strips + +InducedCharge::InducedCharge() {} + +AdvSignal InducedCharge::IntegrateCharge(std::vector DiffusionSignal) +{ + std::vector ResponseSignal; + for (int k = 0; k < DiffusionSignal.size(); k++) + { + std::vector amplitude = DiffusionSignal[k].getAmplitude(); + std::vector surfacepos = DiffusionSignal[k].getSurfacePos(); + std::vector diffusionarea = DiffusionSignal[k].getDiffusionArea(); + + std::vector AffectedStrips; + std::vector temp_AffectedStrips; + std::vector UniqueAffectedStrips; + std::vector ChargeDeposited; + std::vector TotalChargeDeposited; + + Double_t x_start; + Double_t x_end; + Double_t z_start; + Double_t z_end; + Double_t integratedcharge; + + Double_t e = 1.6e-19; + + for (int i = 0; i < surfacepos.size(); i++) + { + //Getting the strips that would see the charge for each energy segment + + AffectedStrips = GetStrips(surfacepos[i], diffusionarea[i]); + + // Calculating the charge deposited by using the error function with the z position of the concerned strip + + for(int j = 0; j < AffectedStrips.size(); j++) + { + x_start = (AffectedStrips[j] - (advsnd::strips / 2))*(advsnd::sensor_width / advsnd::strips) - (stripsensor::inducedcharge::strip_pitch / 2); // check calculation + x_end = (AffectedStrips[j] - (advsnd::strips / 2))*(advsnd::sensor_width / advsnd::strips) + (stripsensor::inducedcharge::strip_pitch / 2); // check calculation + z_start = (x_start - surfacepos[i].X()) / diffusionarea[j]; + z_end = abs(x_end - surfacepos[i].X()) / diffusionarea[j]; + integratedcharge = (erf((z_end) / TMath::Sqrt2()) / 2) - (erf((z_start) / TMath::Sqrt2()) / 2); + temp_AffectedStrips.push_back(AffectedStrips[j]); + ChargeDeposited.push_back(integratedcharge*amplitude[j]); + } + } + + Double_t z = accumulate(amplitude.begin(), amplitude.end(), 0); + + // Getting the unique strips + + std::vector temp_Signalvector; + AdvSignal temp_signal(temp_AffectedStrips, ChargeDeposited); + temp_Signalvector.push_back(temp_signal); + AdvSignal TotalSignal = Combine(temp_Signalvector); + + + TotalChargeDeposited = TotalSignal.getIntegratedSignal(); + UniqueAffectedStrips = TotalSignal.getStrips(); + + Double_t r = accumulate(TotalChargeDeposited.begin(), TotalChargeDeposited.end(), 0); + Double_t rescale_ratio = r/z; + for (int m = 0; m < UniqueAffectedStrips.size(); m++) + { + //the total number of electrons + TotalChargeDeposited[m] = std::ceil((TotalChargeDeposited[m]/stripsensor::frontend::ElectronperADC) * rescale_ratio) ; + } + // Add coupling to neighbour strips + + if (stripsensor::inducedcharge::Coupling) + { + AdvSignal CoupledChargeDeposit = Coupling(TotalChargeDeposited, UniqueAffectedStrips); + ResponseSignal.push_back(CoupledChargeDeposit); + + } else { + AdvSignal PulseSignal(UniqueAffectedStrips, TotalChargeDeposited); + ResponseSignal.push_back(PulseSignal); + } + } + AdvSignal TotalSignal = Combine(ResponseSignal); + return TotalSignal; +} + +AdvSignal InducedCharge::Combine(std::vector Signal) +{ + std::vector CombinedStrips; + std::vector CombinedCharge; + + for (int n = 0; n < Signal.size(); n++) + { + for (int p = 0; p < (Signal[n].getStrips()).size(); p++) + { + CombinedStrips.push_back(Signal[n].getStrips()[p]); + CombinedCharge.push_back(Signal[n].getIntegratedSignal()[p]); + } + } + + std::vector UniqueAffectedStrips = CombinedStrips; + std::vector SummedCharge ; + sort(UniqueAffectedStrips.begin(), UniqueAffectedStrips.end()); + std::vector::iterator itstrip; + itstrip = unique(UniqueAffectedStrips.begin(), UniqueAffectedStrips.end()); + UniqueAffectedStrips.resize(distance(UniqueAffectedStrips.begin(),itstrip)); + + for (int l = 0; l < UniqueAffectedStrips.size(); l++) + { + Double_t temp_totalcharge = 0.0; + auto itfind = find(CombinedStrips.begin(), CombinedStrips.end(), UniqueAffectedStrips[l]); + while (itfind != CombinedStrips.end()) { + temp_totalcharge = temp_totalcharge + CombinedCharge[itfind - CombinedStrips.begin()]; + itfind = find(itfind + 1, CombinedStrips.end(), UniqueAffectedStrips[l]); + } + SummedCharge.push_back(temp_totalcharge); + } + + AdvSignal CombinedSignal(UniqueAffectedStrips, SummedCharge); + + return CombinedSignal; + +} + +std::vector InducedCharge::GetStrips(TVector3 point, Double_t area) +{ + std::vector affectedstrips; + + // Calculating the strips that see a charge in the Nsigma diffusion area + + int fromstrip = floor(((point.X()-(stripsensor::inducedcharge::NSigma*area)) / (advsnd::sensor_width / advsnd::strips)) + (advsnd::strips / 2)); // check calculation + fromstrip = std::max(0, fromstrip); + fromstrip = std::min(advsnd::strips - 1, fromstrip); + + int tostrip = floor(((point.X()+(stripsensor::inducedcharge::NSigma*area)) / (advsnd::sensor_width / advsnd::strips)) + (advsnd::strips / 2)); + tostrip = std::max(0, tostrip); + tostrip = std::min(advsnd::strips - 1, tostrip); + + Int_t N; + N = tostrip - fromstrip; + + for (int i = 0 ; i <= N ; i++) + { + affectedstrips.push_back(fromstrip + i); + } + + return affectedstrips; +} + +std::vector> InducedCharge::GetPulseShape(std::string PulseFileName, std::vector ChargeDeposited) +{ + // For pulse response for time correction, not used currently + + std::vector PulseValues; + + std::ifstream inputFile(PulseFileName); + + if (!inputFile.is_open()) { + std::cout << "Error opening the file!" << std::endl; + } + std::string line; + std::string res_find = "resolution:"; + float res; + std::string s; + + while (getline(inputFile, line)) { + if ((!line.empty()) && (line.substr(0, 1) != "#")) { + std::stringstream ss(line); + if (line.find(res_find) != std::string::npos) { + res = stof(line.substr(line.find(res_find) + res_find.length())); // implement check for + // resolution + } else { + std::string value; + while (getline(ss, value, ' ')) { + PulseValues.push_back(stod(value)); + } + } + } + } + const auto max_value = max_element(PulseValues.begin(), PulseValues.end()); + if (abs(*max_value - 1) > std::numeric_limits::epsilon()) { + throw std::invalid_argument("Maximum value of pulse shape not 1."); + } + + std::vector> PulseResponse; + std::vector temp_response ; + Double_t response_value; + + Int_t sampling_step = stripsensor::sampling / res; + + for(int i = 0; i < ChargeDeposited.size(); i++) + { + Double_t amplitude_max = ChargeDeposited[i]; + + std::default_random_engine generator; + std::normal_distribution dist(stripsensor::noise_mean, stripsensor::noise_std_dev); + response_value = ((stripsensor::baseline + (amplitude_max * stripsensor::amplificaton_factor)) > stripsensor::rail ? stripsensor::rail : (stripsensor::baseline + (amplitude_max * stripsensor::amplificaton_factor))); + temp_response.push_back(response_value + dist(generator)); + PulseResponse.push_back(temp_response); + // if (stripsensor::peakmode == 1) + // { + // response_value = ((stripsensor::baseline + (amplitude_max * stripsensor::amplificaton_factor)) > stripsensor::rail ? stripsensor::rail : (stripsensor::baseline + (amplitude_max * stripsensor::amplificaton_factor))); + // temp_response.push_back(response_value + dist(generator)); + // PulseResponse.push_back(temp_response); + // } + // else + // { + // response_value = ((stripsensor::baseline + (amplitude_max * stripsensor::amplificaton_factor)) > stripsensor::rail ? stripsensor::rail : (stripsensor::baseline + (amplitude_max * stripsensor::amplificaton_factor))); + // temp_response.push_back(response_value + dist(generator)); + // PulseResponse.push_back(temp_response); + // } + } + + return PulseResponse; + // get the vector of beginning to max of the pulse + // time response not included! +} + +AdvSignal InducedCharge::Coupling(std::vector TotalCharge, std::vector AffectedStrips) +{ + + // Using the coupling constants, charge sharing between strips and their neighbours + + Int_t couplingsize = stripsensor::inducedcharge::CouplingConstants.size(); + std::vector CoupledCharge; + for (int k = 1; k < couplingsize; k++){ + AffectedStrips.emplace(AffectedStrips.begin(), AffectedStrips[k-1]- k); + AffectedStrips.push_back(AffectedStrips[AffectedStrips.size()-k] + k); + TotalCharge.emplace(TotalCharge.begin(), 0); + TotalCharge.push_back(0); + + } + + CoupledCharge = TotalCharge; + for(int n = 0; n < TotalCharge.size(); n++) + { + CoupledCharge[n] = TotalCharge[n]*stripsensor::inducedcharge::CouplingConstants[0]; + } + + for (int i = 0; i < TotalCharge.size(); i++) + { + if (TotalCharge[i] == 0) + continue; + for (int j = 1; j < couplingsize; j++) + { + CoupledCharge[i-j] += TotalCharge[i]*stripsensor::inducedcharge::CouplingConstants[j]; + CoupledCharge[i+j] += TotalCharge[i]*stripsensor::inducedcharge::CouplingConstants[j]; + } + + + } + AdvSignal Coupledresult(AffectedStrips, CoupledCharge); + return Coupledresult; +} \ No newline at end of file diff --git a/shipLHC/digitisation/InducedCharge.h b/shipLHC/digitisation/InducedCharge.h new file mode 100644 index 00000000..f3461722 --- /dev/null +++ b/shipLHC/digitisation/InducedCharge.h @@ -0,0 +1,22 @@ +#ifndef SHIPLHC_INDUCEDCHARGE_H_ +#define SHIPLHC_INDUCEDCHARGE_H_ + +#include "EnergyFluctUnit.h" +#include "SurfaceSignal.h" +#include "AdvSignal.h" + +#include +#include + +class InducedCharge +{ + public: + InducedCharge(); + AdvSignal IntegrateCharge(std::vector DiffusionSignal); + AdvSignal Combine(std::vector Signal); + std::vector GetStrips(TVector3 point, Double_t area); + std::vector> GetPulseShape(std::string PulseFileName, std::vector ChargeDeposited); + AdvSignal Coupling(std::vector TotalCharge, std::vector AffectedStrips); +}; + +#endif \ No newline at end of file diff --git a/shipLHC/digitisation/SiDigiParameters.h b/shipLHC/digitisation/SiDigiParameters.h new file mode 100644 index 00000000..eb3d69c7 --- /dev/null +++ b/shipLHC/digitisation/SiDigiParameters.h @@ -0,0 +1,72 @@ +#ifndef SHIPLHC_SIDIGIPARAMETERS_H_ +#define SHIPLHC_SIDIGIPARAMETERS_H_ + +#include "ShipUnit.h" + +/* Parameters for strip digitization */ + +namespace stripsensor { + +const bool peakmode = 0; +const Double_t sampling = 25; +const std::string APVpeakpulse = "advsndsw/shipLHC/data/APVShapePeak_default.txt"; +const std::string APVdecopulse = "advsndsw/shipLHC/data/APVShapeDeco_default.txt"; + +const Double_t value_per_mip = 1; // in mA +const Double_t amplificaton_factor = value_per_mip/(30000 * 1.6 * 10e-19); // output in mA +const Double_t rail = 4; // in mA +const Double_t baseline = 2; // in mA +const Double_t threshold = 2.5; // in mA +const Double_t mode = 20 ; // in MHz - the frequency that output is cycled + +const Double_t noise_mean = 0; +const Double_t noise_std_dev = 0.01; + +namespace chargedivision{ + const Int_t ChargeDivisionsperStrip = 10; // Number of divisions of the track length per strip + const float StripPitch = 120e-4; // Strip pitch of the detector +} + +namespace drift{ + const bool cms_approximation = 0; // + const Double_t module_thickness = 0.05; //check value + const int depletion_voltage = 170; + const int applied_voltage = 300; + const int charge_mobility = 480; + const Double_t chargedistributionRMS = 6.5e-10; + const int temperature = 300; + const Double_t perGeV = 3.61; + +} + +namespace inducedcharge{ + const Int_t NSigma = 3; + const Double_t strip_width = 30e-4; //in cm + const Double_t strip_pitch = 120e-4; //in cm + const bool Coupling = 0; + const std::vector CouplingConstants = {0.964, 0.018}; +} + +namespace frontend{ + const Int_t muxgain_register = 2; // 5 possible values with middle sized register giving gain of 1 mA/mip with +/- 20% + const Double_t reference_current = 0.128; // in mA, above which is the gain + // const bool peakmode = 0; // 0 for peak mode operation of the APV, 1 for deconvolution mode + const Double_t gain = 0.1; // in uA. Need to check which value to use + const bool write_digi_to_text = 1; + const Int_t ElectronperADC = 250; + const bool NoiseOption = 1; + const bool CMNoiseOption = 0; + const bool ZSModeOption = 0; + const Int_t StripNoise = 1; + const Double_t NoiseRMS = 2; + const Double_t NoiseSigmaThreshold = 2; + const Double_t CMNoise = 2; + const Int_t NumberofStrips = 768; + const Int_t ZeroSuppressionMode = 4; + const Int_t ZeroSuppressionMode1T = 2; +} +} // namespace advsnd + +#endif // SHIPLHC_SIDIGIPARAMETERS_H_ + + diff --git a/shipLHC/digitisation/SiG4UniversalFluctuation.cxx b/shipLHC/digitisation/SiG4UniversalFluctuation.cxx new file mode 100644 index 00000000..25c5cdf8 --- /dev/null +++ b/shipLHC/digitisation/SiG4UniversalFluctuation.cxx @@ -0,0 +1,235 @@ +//Taken from GEANT4 class G4UniversalFluctuation with some modifications + + +// ******************************************************************** +// * License and Disclaimer * +// * * +// * The Geant4 software is copyright of the Copyright Holders of * +// * the Geant4 Collaboration. It is provided under the terms and * +// * conditions of the Geant4 Software License, included in the file * +// * LICENSE and available at http://cern.ch/geant4/license . These * +// * include a list of copyright holders. * +// * * +// * Neither the authors of this software system, nor their employing * +// * institutes,nor the agencies providing financial support for this * +// * work make any representation or warranty, express or implied, * +// * regarding this software system or assume any liability for its * +// * use. Please see the license in the file LICENSE and URL above * +// * for the full disclaimer and the limitation of liability. * +// * * +// * This code implementation is the result of the scientific and * +// * technical work of the GEANT4 collaboration. * +// * By using, copying, modifying or distributing the software (or * +// * any work based on the software) you agree to acknowledge its * +// * use in resulting scientific publications, and indicate your * +// * acceptance of all terms of the Geant4 Software license. * +// ******************************************************************** +// +// +// ------------------------------------------------------------------- +// +// GEANT4 Class file +// +// +// File name: G4UniversalFluctuation +// +// Author: V. Ivanchenko for Laszlo Urban +// +// Creation date: 03.01.2002 +// +// Modifications: +// +// + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +// #include +// #include +// #include + +#include "ShipUnit.h" +#include "TF1.h" +#include "TMath.h" +#include "TRandom.h" +#include + +#include +#include +#include +using namespace std; + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... +SiG4UniversalFluctuation::SiG4UniversalFluctuation() { rndmarray = new Double_t[sizearray]; } +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +SiG4UniversalFluctuation::~SiG4UniversalFluctuation() { delete[] rndmarray; } + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +Double_t SiG4UniversalFluctuation::SampleFluctuations(Double_t particleMass, + Double_t q, + Double_t averageLoss, + Double_t mom, + Double_t length) +{ + // Calculate actual loss from the mean loss. + // The model used to get the fluctuations is essentially the same + // as in Glandz in Geant3 (Cern program library W5013, phys332). + // L. Urban et al. NIM A362, p.416 (1995) and Geant4 Physics Reference Manual + + // shortcut for very small loss or from a step nearly equal to the range + // (out of validity of the model) + + m_Inv_particleMass = 1.0 / particleMass; + m_massrate = ShipUnit::electron_mass_c2 * m_Inv_particleMass; + chargeSquare = q * q; + + Double_t energyLoss = 0; + + minLoss = 10. * ShipUnit::eV; + if (averageLoss < minLoss) { + energyLoss = averageLoss; + return energyLoss; + } + + meanLoss = averageLoss; + + Double_t tkin = (mom * mom) / (particleMass * particleMass) + 1.0; + const Double_t gam = tkin * m_Inv_particleMass + 1.0; + const Double_t gam2 = gam * gam; + const Double_t beta = 1.0 - 1.0 / gam2; + const Double_t beta2 = beta * beta; + + Double_t loss(0.), siga(0.); + + // Gaussian regime + // for heavy particles only and conditions + // for Gauusian fluct. has been changed + + tmax = 2. * ShipUnit::electron_mass_c2 * beta2 * gam2 / (1. + m_massrate * (2. * gam + m_massrate)); + + if (particleMass > ShipUnit::electron_mass_c2 && meanLoss >= minNumberInteractionsBohr * tcut + && tmax <= 2. * tcut) { + + siga = + std::sqrt((tmax / beta2 - 0.5 * tcut) * ShipUnit::twopi_mc2_rcl2 * length * chargeSquare * electronDensity); + + const Double_t sn = meanLoss / siga; + + // thick target case + if (sn >= 2.0) { + + const Double_t twomeanLoss = meanLoss + meanLoss; + do { + + gRandom->SetSeed(0); + TRandom* rngSaved = static_cast(gRandom->Clone()); + energyLoss = gRandom->Gaus(meanLoss, siga); + // // Loop checking, 03-Aug-2015, Vladimir Ivanchenko + } while (0.0 > loss || twomeanLoss < loss); + } else { + const Double_t neff = sn * sn; + TF1* f1 = new TF1("Gamma(x)", "ROOT::Math::tgamma(x)", neff, 1.0); + double r = f1->GetRandom(); + gRandom->SetSeed(0); + TRandom* rngSaved = static_cast(gRandom->Clone()); + energyLoss = meanLoss * r / neff; + } + return energyLoss; + + } + if (tcut <= e0) { + return meanLoss; + } + + const Double_t scaling = std::min(1. + 0.5 * ShipUnit::keV / tcut, 1.50); + meanLoss /= scaling; + + w2 = (tcut > ipotFluct) ? TMath::Log(2. * ShipUnit::electron_mass_c2 * beta2 * gam2) - beta2 : 0.0; + + return SampleGlandz() * scaling; +} + +Double_t SiG4UniversalFluctuation::SampleGlandz() +{ + Double_t a1(0.0), a3(0.0); + Double_t loss = 0.0; + Double_t e1 = ipotFluct; + Double_t rate = rateFluct; + + if (tcut > e1) { + a1 = meanLoss * (1. - rate) / e1; + if (a1 < a0) { + const Double_t fwnow = 0.1 + (fw - 0.1) * std::sqrt(a1 / a0); + a1 /= fwnow; + e1 *= fwnow; + } else { + a1 /= fw; + e1 *= fw; + } + } + + const Double_t w1 = tcut / e0; + a3 = rate * meanLoss * (tcut - e0) / (e0 * tcut * TMath::Log(w1)); + if (a1 <= 0.) { + a3 /= rate; + } + + //'nearly' Gaussian fluctuation if a1>nmaxCont&&a2>nmaxCont&&a3>nmaxCont + Double_t emean = 0.; + Double_t sig2e = 0.; + + // AddExcitation(a1, e1, emean, loss, sig2e); + + // excitation of type 1 + if (a1 > 0.0) { + AddExcitation(a1, e1, emean, loss, sig2e); + } + + if (sig2e > 0.0) { + SampleGauss(emean, sig2e, loss); + } + + // ionisation + if (a3 > 0.) { + emean = 0.; + sig2e = 0.; + Double_t p3 = a3; + Double_t alfa = 1.; + if (a3 > nmaxCont) { + alfa = w1 * (nmaxCont + a3) / (w1 * nmaxCont + a3); + const Double_t alfa1 = alfa * TMath::Log(alfa) / (alfa - 1.); + const Double_t namean = a3 * w1 * (alfa - 1.) / ((w1 - 1.) * alfa); + emean += namean * e0 * alfa1; + sig2e += e0 * e0 * namean * (alfa - alfa1 * alfa1); + p3 = a3 - namean; + } + + const Double_t w3 = alfa * e0; + if (tcut > w3) { + gRandom->SetSeed(0); + TRandom* rndm = static_cast(gRandom->Clone()); + const Double_t w = (tcut - w3) / tcut; + const Int_t nnb = (Int_t)gRandom->Poisson(p3); + if (nnb > 0) { + if (nnb > sizearray) { + sizearray = nnb; + delete[] rndmarray; + rndmarray = new Double_t[nnb]; + } + for (Int_t k = 0; k < nnb; ++k) { + // check if this loop is needed to handle the NaN values + if(rndmarray[k]!=rndmarray[k]){loss += w3 / (1. - w * 1e-310);} + else{loss += w3 / (1. - w * rndmarray[k]);} + } + } + + } + + if (sig2e > 0.0) { + SampleGauss(emean, sig2e, loss); + // if(loss != loss) {cout << sig2e << endl; } + } + } + return loss; +} diff --git a/shipLHC/digitisation/SiG4UniversalFluctuation.h b/shipLHC/digitisation/SiG4UniversalFluctuation.h new file mode 100644 index 00000000..c218d0b3 --- /dev/null +++ b/shipLHC/digitisation/SiG4UniversalFluctuation.h @@ -0,0 +1,113 @@ +// +// GEANT4 Class header file +// +// +// File name: SiG4UniversalFluctuation +// +// Author: Vladimir Ivanchenko make a class for Laszlo Urban model +// +// Modified for standalone use in CMSSW. Danek K. 02/2006 +// +// Class Description: +// +// Implementation of energy loss fluctuations in Silicon + +// ------------------------------------------------------------------- +// + +#ifndef SiG4UniversalFluctuation_h +#define SiG4UniversalFluctuation_h + +#include "TRandom.h" +#include +using namespace std; + +namespace CLHEP { +class HepRandomEngine; +} + +class SiG4UniversalFluctuation +{ + public: + SiG4UniversalFluctuation(); + ~SiG4UniversalFluctuation(); + Double_t SampleFluctuations(Double_t ParticleMass, + Double_t ParticleCharge, + Double_t averageLoss, + Double_t mom, + Double_t length); + Double_t SampleGlandz(); + Int_t test_temp = 0; + + protected: + Double_t m_Inv_particleMass; + Double_t m_massrate; + Double_t chargeSquare; + + Double_t minLoss; + Double_t meanLoss; + Double_t tmax; + Double_t massrate; + + Double_t minNumberInteractionsBohr = 10.0; + Double_t tcut = 0.120425; // in MeV, value taken from CMS + Double_t electronDensity = 6.797E+20; // value taken from CMS + Double_t e0 = 1.E-5; // material->GetIonisation()->GetEnergy0fluct(); + Double_t ipotFluct = 0.0001736; // material->GetIonisation()->GetMeanExcitationEnergy(); + Double_t ipotLogFluct = -8.659; // material->GetIonisation()->GetLogMeanExcEnergy(); + Double_t rateFluct = 0.4; // material->GetIonisation()->GetRateionexcfluct(); + + Double_t a0 = 42.0; + Double_t fw = 4.0; + Double_t nmaxCont = 8.0; + Double_t w2 = 0.0; + Int_t sizearray = 30; + Double_t* rndmarray = nullptr; + + + + // virtual Double_t SampleGlandz(); + inline void AddExcitation(const Double_t ax, const Double_t ex, Double_t& eav, Double_t& eloss, Double_t& esig2); + + inline void SampleGauss(const Double_t eav, const Double_t esig2, Double_t& eloss); +}; + +inline void SiG4UniversalFluctuation::AddExcitation(const Double_t ax, + const Double_t ex, + Double_t& eav, + Double_t& eloss, + Double_t& esig2) +{ + if (ax > nmaxCont) { + eav += ax * ex; + esig2 += ax * ex * ex; + } else { + gRandom->SetSeed(0); + TRandom* rngSaved = static_cast(gRandom->Clone()); + const Int_t p = (Int_t)gRandom->Poisson(ax); + if (p > 0) { + eloss += ((p + 1) - 2. * gRandom->Uniform()) * ex; + } + } +} + +inline void SiG4UniversalFluctuation::SampleGauss(const Double_t eav, const Double_t esig2, Double_t& eloss) +{ + Double_t x = eav; + const Double_t sig = std::sqrt(esig2); + gRandom->SetSeed(0); + TRandom* rndm = static_cast(gRandom->Clone()); + if (eav < 0.25 * sig) { + x += (2. * rndm->Uniform() - 1.) * eav; + } else { + do { + x = rndm->Gaus(eav, sig); + } while (x < 0.0 || x > 2 * eav); + // Loop checking, 23-Feb-2016, Vladimir Ivanchenko + } + eloss += x; +} + +//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... + +#endif diff --git a/shipLHC/digitisation/StripNoise.cxx b/shipLHC/digitisation/StripNoise.cxx new file mode 100644 index 00000000..f9a90c72 --- /dev/null +++ b/shipLHC/digitisation/StripNoise.cxx @@ -0,0 +1,132 @@ +#include "StripNoise.h" +#include "SiDigiParameters.h" +#include "AdvSignal.h" +#include "TRandom.h" + + +#include +#include +#include +#include +#include +using namespace std; + +/* To Do : +Add Baseline Shift? +Add Pedestals for each channel after calibration */ + +StripNoise::StripNoise() {} + +AdvSignal StripNoise::AddGaussianNoise(AdvSignal Signal) +{ + std::vector Strips = Signal.getStrips(); + std::vector Amplitude = Signal.getIntegratedSignal(); + + TRandom* rndm = gRandom; + for (int i = 0; i < stripsensor::frontend::NumberofStrips; i++) + { + Double_t x = rndm->Gaus(0, stripsensor::frontend::NoiseRMS); + Strips.push_back(i); + Amplitude.push_back(x); + } + AdvSignal NoiseSignal(Strips, Amplitude); + return NoiseSignal; +} + +AdvSignal StripNoise::AddGaussianTailNoise(AdvSignal Signal) +{ + std::vector Strips = Signal.getStrips(); + Int_t StripsSize = Strips.size(); + std::vector Amplitude = Signal.getIntegratedSignal(); + + std::vector NoisyStrips; + std::vector NoisyAmplitudes; + + std::vector NoiseAddedStrips = Strips; + std::vector NoiseAddedAmplitudes = Amplitude; + + TRandom* rndm = gRandom; + for (int j = 0; j < StripsSize; j++) + { + Amplitude[j] += rndm->Gaus(0, stripsensor::frontend::NoiseRMS); + } + + double probabilityLeft = 0.5 * std::erfc(stripsensor::frontend::NoiseSigmaThreshold / std::sqrt(2.0)); + + Double_t MeanNumberofNoisyChannels = probabilityLeft * stripsensor::frontend::NumberofStrips; + + Double_t NumberofNoisyChannels = int(rndm->Poisson(MeanNumberofNoisyChannels)); + for (int i = 0; i < NumberofNoisyChannels; i++) + { + Strips.push_back(rndm->Integer(stripsensor::frontend::NumberofStrips)); + Amplitude.push_back(generate_gaussian_tail(stripsensor::frontend::NoiseSigmaThreshold*stripsensor::frontend::NoiseRMS, stripsensor::frontend::NoiseRMS)); + } + AdvSignal NoiseSignal(Strips, Amplitude); + + return NoiseSignal; + +} + +AdvSignal StripNoise::AddCMNoise(AdvSignal Signal) +{ + Double_t CMnoise = stripsensor::frontend::CMNoise; + std::vector Strips = Signal.getStrips(); + std::vector Charge = Signal.getIntegratedSignal(); + Int_t nAPVs = int(stripsensor::frontend::NumberofStrips / 128); + std::vector CMVector ; + TRandom* rndm = gRandom; + for (int i = 0 ; i < stripsensor::frontend::NumberofStrips; i += 128) + { + CMVector.push_back(rndm->Gaus(0, CMnoise)); + } + + for (int j = 0; j < Strips.size(); j ++) + { + Charge[j] += CMVector[int(Strips[j]/128)]; + } + + AdvSignal NoiseSignal(Strips, Charge); + return NoiseSignal; +} + +double StripNoise::generate_gaussian_tail(const double a,const double sigma) +{ + // taken from CMS + + /* Returns a gaussian random variable larger than a + * This implementation does one-sided upper-tailed deviates. + */ + TRandom* rndm = gRandom; + double s = a / sigma; + + if (s < 1) { + /* + For small s, use a direct rejection method. The limit s < 1 + can be adjusted to optimise the overall efficiency + */ + double x; + + do { + x = rndm->Gaus(0., 1.0); + } while (x < s); + return x * sigma; + + } else { + /* Use the "supertail" deviates from the last two steps + * of Marsaglia's rectangle-wedge-tail method, as described + * in Knuth, v2, 3rd ed, pp 123-128. (See also exercise 11, p139, + * and the solution, p586.) + */ + + double u, v, x; + + do { + u = rndm->Rndm(); + do { + v = rndm->Rndm(); + } while (v == 0.0); + x = sqrt(s * s - 2 * log(v)); + } while (x * u > s); + return x * sigma; + } +} diff --git a/shipLHC/digitisation/StripNoise.h b/shipLHC/digitisation/StripNoise.h new file mode 100644 index 00000000..f03c0b21 --- /dev/null +++ b/shipLHC/digitisation/StripNoise.h @@ -0,0 +1,23 @@ +#ifndef SHIPLHC_STRIPNOISE_H_ +#define SHIPLHC_STRIPNOISE_H_ + +#include "StripNoise.h" +#include "AdvSignal.h" + +#include +#include + +class StripNoise +{ + public: + StripNoise(); + AdvSignal AddGaussianNoise(AdvSignal Signal); + AdvSignal AddGaussianTailNoise(AdvSignal Signal); + AdvSignal AddCMNoise(AdvSignal Signal); + double generate_gaussian_tail(const double a, const double sigma); + + void TestingNoise(); + +}; + +#endif \ No newline at end of file diff --git a/shipLHC/digitisation/SurfaceSignal.h b/shipLHC/digitisation/SurfaceSignal.h new file mode 100644 index 00000000..34bc0a35 --- /dev/null +++ b/shipLHC/digitisation/SurfaceSignal.h @@ -0,0 +1,37 @@ +#ifndef SHIPLHC_SURFACESIGNAL_H_ +#define SHIPLHC_SURFACESIGNAL_H_ + +#include "TVector3.h" + +#include + +class SurfaceSignal +{ + public: + SurfaceSignal() + : DiffusionArea_(), SurfacePos_(), Amplitude_() + { + } + + SurfaceSignal(std::vector DiffusionArea, std::vector SurfacePos, std::vector Amplitude) + : DiffusionArea_(DiffusionArea) + , SurfacePos_(SurfacePos) + , Amplitude_(Amplitude) + { + } + + std::vector getDiffusionArea() const { return DiffusionArea_; } + std::vector getSurfacePos() const { return SurfacePos_; } + std::vector getAmplitude() const { return Amplitude_; } + + void setDiffusionArea(std::vector value) { DiffusionArea_ = value; } + void setSurfacePos(std::vector value) { SurfacePos_ = value; } + void setAmplitude(std::vector value) { Amplitude_ = value; } + + private: + std::vector DiffusionArea_; + std::vector SurfacePos_; + std::vector Amplitude_; +}; + +#endif diff --git a/shipLHC/digitisation/data/APVShapePeak_default.txt b/shipLHC/digitisation/data/APVShapePeak_default.txt new file mode 100644 index 00000000..34cb3643 --- /dev/null +++ b/shipLHC/digitisation/data/APVShapePeak_default.txt @@ -0,0 +1,9 @@ +# Automatically generated using parametrizePulse::generateCode(low=-30, high=35, step=0.1) +# That pulse shapes correspond to the ones described in CMS NOTE 2007/027 +# with tau=50ns and delta=20ns +# It is fairly similar to the previous analytical forms, except in the tails +# The max value has to be 1 +# resolution is in ns +resolution: 0.5 + +0.000598718 0.00179789 0.00380417 0.00622833 0.00934726 0.0133483 0.0172689 0.0225072 0.0282688 0.0339136 0.040572 0.0482572 0.0557854 0.063159 0.0727605 0.082258 0.0915594 0.101262 0.112735 0.123969 0.134968 0.146844 0.159929 0.172739 0.185278 0.199281 0.21373 0.227871 0.242102 0.257926 0.273412 0.288564 0.305358 0.321845 0.338351 0.355808 0.372906 0.390895 0.408861 0.427103 0.445108 0.462713 0.479924 0.496746 0.513186 0.52925 0.544943 0.560272 0.575241 0.589857 0.604126 0.618052 0.631641 0.644898 0.657829 0.670439 0.682733 0.694716 0.706393 0.717768 0.728848 0.739635 0.750136 0.760355 0.770296 0.779964 0.789363 0.798498 0.807373 0.815993 0.824361 0.832481 0.840359 0.847997 0.8554 0.862572 0.869516 0.876237 0.882738 0.889022 0.895095 0.900958 0.906616 0.912073 0.917331 0.922394 0.927265 0.931949 0.936447 0.940763 0.944901 0.948864 0.952653 0.956274 0.959728 0.963019 0.966149 0.969122 0.97194 0.974606 0.977122 0.979492 0.981718 0.983803 0.985749 0.987559 0.989235 0.990781 0.992197 0.993487 0.994654 0.995699 0.996624 0.997433 0.998126 0.998708 0.999178 0.999541 0.999797 0.99995 1 0.99995 0.999803 0.999559 0.999221 0.998791 0.99827 0.997661 0.996966 0.996185 0.995321 0.994376 0.993351 0.992248 0.991068 0.989814 0.988487 0.987088 0.985619 0.984081 0.982477 0.980807 0.979073 0.977276 0.975419 0.973501 0.971525 0.969492 0.967403 0.96526 0.963064 0.960816 0.958517 0.956169 0.953772 0.951329 0.94884 0.946306 0.943729 0.941109 0.938448 0.935747 0.933006 0.930228 0.927412 0.924561 0.921674 0.918753 0.915799 0.912813 0.909796 0.906748 0.903671 0.900566 0.897432 0.894272 0.891086 0.887875 0.884639 0.88138 0.878099 0.874795 0.87147 0.868125 0.86476 0.861376 0.857973 0.854553 0.851117 0.847664 0.844195 0.840712 0.837214 0.833703 0.830178 0.826641 0.823093 0.819533 0.815963 0.812382 0.808792 0.805193 0.801586 0.79797 0.794347 0.790718 0.787081 0.783439 0.779792 0.776139 0.772482 0.768821 0.765157 0.761489 0.757818 0.754145 0.75047 0.746794 0.743116 0.739438 0.735759 0.73208 0.728402 0.724724 0.721048 0.717372 0.713699 0.710028 0.706359 0.702692 0.699029 0.695369 0.691713 0.688061 0.684413 0.680769 0.67713 0.673496 0.669868 0.666245 0.662627 0.659016 0.655411 0.651812 0.648221 0.644636 0.641058 0.637488 0.633925 0.63037 0.626823 0.623284 0.619754 0.616232 0.612719 0.609215 0.605719 0.602233 0.598757 0.59529 0.591833 0.588385 0.584948 0.581521 0.578104 0.574697 0.571301 0.567916 0.564541 0.561178 0.557825 0.554484 0.551154 0.547835 0.544528 0.541232 0.537948 0.534676 0.531416 0.528167 0.524931 0.521707 0.518495 0.515295 0.512107 0.508932 0.50577 0.50262 0.499482 0.496358 0.493246 0.490146 0.48706 0.483986 0.480926 0.477878 0.474844 0.471822 0.468814 0.465819 0.462837 0.459868 0.456913 0.45397 0.451041 0.448126 0.445224 0.442335 0.43946 0.436598 0.433749 0.430914 0.428092 0.425284 0.42249 0.419709 0.416941 0.414187 0.411446 0.408719 0.406006 0.403306 0.40062 0.397947 0.395287 0.392641 0.390009 0.38739 0.384785 0.382193 0.379615 0.37705 0.374499 0.371961 0.369436 0.366925 0.364427 0.361943 0.359472 0.357014 0.35457 0.352139 0.349721 0.347317 0.344926 0.342547 0.340183 0.337831 0.335492 0.333167 0.330854 0.328555 0.326268 0.323995 0.321734 0.319487 0.317252 0.31503 0.312821 0.310625 0.308441 0.30627 0.304112 0.301966 0.299833 0.297713 0.295605 0.293509 0.291426 0.289356 0.287297 0.285252 0.283218 0.281196 0.279187 0.27719 0.275205 0.273232 0.271271 0.269322 0.267385 0.26546 0.263546 0.261645 0.259755 0.257877 0.25601 0.254156 0.252312 0.250481 0.24866 0.246852 0.245054 0.243268 0.241493 0.239729 0.237977 0.236236 0.234506 0.232786 0.231078 0.229381 0.227695 0.226019 0.224355 0.222701 0.221057 0.219425 0.217803 0.216191 0.214591 0.213 0.21142 0.20985 0.208291 0.206742 0.205203 0.203674 0.202155 0.200647 0.199148 0.19766 0.196181 0.194712 0.193253 0.191804 0.190364 0.188934 0.187514 0.186103 0.184702 0.18331 0.181928 0.180555 0.179191 0.177836 0.176491 0.175155 0.173828 0.17251 0.171201 0.169901 0.16861 0.167328 0.166055 0.16479 0.163535 0.162287 0.161049 0.159819 0.158598 0.157385 0.15618 0.154984 0.153796 0.152617 0.151446 0.150283 0.149128 0.147981 0.146842 0.145712 0.144589 0.143474 0.142367 0.141268 0.140177 0.139093 0.138017 0.136949 0.135888 0.134835 0.133789 0.132751 0.131721 0.130697 0.129681 0.128672 0.127671 0.126676 0.125689 0.124709 0.123736 0.12277 0.121811 0.120859 0.119914 0.118975 0.118044 0.117119 0.116201 0.115289 0.114384 0.113486 0.112594 0.111709 0.110831 0.109958 0.109092 0.108233 0.10738 0.106533 0.105692 0.104857 0.104029 0.103207 0.10239 0.10158 0.100776 0.0999777 0.0991854 0.0983989 0.0976182 0.0968433 0.0960742 0.0953108 0.094553 0.093801 0.0930545 0.0923136 0.0915782 0.0908483 0.0901239 0.0894049 0.0886913 0.087983 0.08728 0.0865824 0.0858899 0.0852027 0.0845206 0.0838437 0.0831719 0.0825052 0.0818435 0.0811867 0.080535 0.0798882 0.0792463 0.0786092 0.077977 0.0773496 0.0767269 0.076109 0.0754958 0.0748872 0.0742833 0.073684 0.0730893 0.0724991 0.0719134 0.0713321 0.0707554 0.070183 0.0696151 0.0690514 0.0684922 0.0679372 0.0673864 0.0668399 0.0662976 0.0657595 0.0652255 0.0646957 0.0641699 0.0636482 0.0631305 0.0626168 0.0621071 0.0616013 0.0610995 0.0606015 0.0601074 0.0596171 0.0591307 0.058648 0.058169 0.0576938 0.0572223 0.0567545 0.0562903 0.0558297 0.0553727 0.0549193 0.0544694 0.0540231 0.0535802 0.0531408 0.0527048 0.0522722 0.0518431 0.0514173 0.0509948 0.0505757 0.0501598 0.0497472 0.0493379 0.0489318 0.0485289 0.0481291 0.0477325 0.0473391 0.0469487 0.0465614 0.0461772 0.045796 0.0454178 0.0450427 0.0446704 0.0443012 0.0439348 0.0435714 0.0432109 0.0428532 0.0424983 0.0421463 0.041797 0.0414506 0.0411069 0.0407659 0.0404277 0.0400921 0.0397593 0.039429 0.0391014 0.0387765 0.0384541 0.0381343 0.037817 0.0375023 0.0371902 0.0368805 0.0365733 0.0362685 0.0359662 0.0356663 0.0353689 0.0350738 0.0347811 0.0344907 0.0342027 0.033917 0.0336336 0.0333525 0.0330736 0.032797 0.0325226 0.0322504 0.0319805 0.0317127 0.031447 0.0311836 0.0309222 0.030663 0.0304058 0.0301508 0.0298978 0.0296468 0.0293979 0.029151 0.0289061 0.0286632 0.0284223 0.0281833 0.0279463 0.0277112 0.027478 0.0272467 0.0270173 0.0267897 0.026564 0.0263402 0.0261181 0.0258979 0.0256795 0.0254628 0.025248 0.0250349 0.0248235 0.0246138 0.0244059 0.0241997 0.0239951 0.0237923 0.023591 0.0233915 0.0231936 0.0229973 0.0228026 0.0226095 0.022418 0.022228 0.0220397 0.0218528 0.0216675 0.0214838 0.0213015 0.0211208 0.0209415 0.0207637 0.0205874 0.0204125 0.0202391 0.0200671 0.0198965 0.0197273 0.0195595 0.0193931 0.0192281 0.0190645 0.0189022 0.0187412 0.0185816 0.0184233 0.0182663 0.0181106 0.0179562 0.0178031 0.0176512 0.0175006 0.0173513 0.0172032 0.0170563 0.0169106 0.0167662 0.0166229 0.0164809 0.01634 0.0162003 0.0160617 0.0159244 0.0157881 0.015653 0.015519 0.0153862 0.0152544 0.0151237 0.0149942 0.0148657 0.0147383 0.0146119 0.0144866 0.0143624 0.0142391 0.014117 0.0139958 0.0138756 0.0137565 0.0136384 0.0135212 0.013405 0.0132898 0.0131756 0.0130623 0.01295 0.0128386 0.0127281 0.0126186 0.01251 0.0124023 0.0122955 0.0121896 0.0120846 0.0119805 0.0118773 0.0117749 0.0116734 0.0115727 0.0114729 0.0113739 0.0112758 0.0111785 0.011082 0.0109863 0.0108914 0.0107974 0.0107041 0.0106116 0.0105199 0.010429 0.0103388 0.0102494 0.0101607 0.0100728 0.00998569 0.00989927 0.00981358 0.00972862 0.00964437 0.00956084 0.00947802 0.0093959 0.00931447 0.00923374 0.00915369 0.00907432 0.00899562 0.00891759 0.00884022 0.00876351 0.00868745 0.00861203 0.00853726 0.00846312 0.00838961 0.00831673 0.00824446 0.00817281 0.00810178 0.00803134 0.0079615 0.00789226 0.00782361 0.00775555 0.00768806 0.00762115 0.00755481 0.00748903 0.00742382 0.00735916 0.00729506 0.0072315 0.00716848 0.007106 0.00704406 0.00698264 0.00692175 0.00686138 0.00680153 0.00674218 0.00668335 0.00662502 0.00656719 0.00650985 0.006453 0.00639664 0.00634077 0.00628537 0.00623044 0.00617599 0.006122 0.00606848 0.00601542 0.00596281 0.00591065 0.00585894 0.00580767 0.00575685 0.00570646 0.0056565 0.00560697 0.00555787 0.00550919 0.00546093 0.00541308 0.00536565 0.00531862 0.005272 0.00522578 0.00517995 0.00513452 0.00508949 0.00504484 0.00500057 0.00495669 0.00491318 0.00487005 0.00482729 0.0047849 0.00474288 0.00470122 0.00465992 0.00461897 0.00457838 0.00453814 0.00449824 0.00445869 0.00441948 0.00438062 0.00434208 0.00430388 0.00426601 0.00422847 0.00419125 0.00415436 0.00411778 0.00408152 0.00404558 0.00400994 0.00397462 0.0039396 0.00390488 0.00387046 0.00383635 0.00380253 0.003769 0.00373576 0.00370281 0.00367015 0.00363777 0.00360567 0.00357385 0.0035423 0.00351103 0.00348003 0.0034493 0.00341884 0.00338864 0.00335871 0.00332903 0.00329962 0.00327046 0.00324155 0.00321289 0.00318449 0.00315633 0.00312841 0.00310074 0.00307331 0.00304612 0.00301916 diff --git a/shipLHC/shipLHCLinkDef.h b/shipLHC/shipLHCLinkDef.h old mode 100755 new mode 100644 index 0bedf3ef..7498ef91 --- a/shipLHC/shipLHCLinkDef.h +++ b/shipLHC/shipLHCLinkDef.h @@ -4,26 +4,35 @@ #pragma link off all classes; #pragma link off all functions; -#pragma link C++ class Floor+; -#pragma link C++ class boxTarget+; -#pragma link C++ class EmulsionDet+; -#pragma link C++ class EmulsionDetPoint+; -#pragma link C++ class EmulsionDetContFact+; -#pragma link C++ class Scifi+; -#pragma link C++ class ScifiPoint+; -#pragma link C++ class MuFilter+; -#pragma link C++ class MuFilterPoint+; -#pragma link C++ class MuFilterHit+; -#pragma link C++ class AdvTargetHit+; -#pragma link C++ class AdvMuFilterHit+; -#pragma link C++ class sndScifiHit+; +#pragma link C++ class Floor + ; +#pragma link C++ class boxTarget + ; +#pragma link C++ class EmulsionDet + ; +#pragma link C++ class EmulsionDetPoint + ; +#pragma link C++ class EmulsionDetContFact + ; +#pragma link C++ class Scifi + ; +#pragma link C++ class ScifiPoint + ; +#pragma link C++ class MuFilter + ; +#pragma link C++ class MuFilterPoint + ; +#pragma link C++ class MuFilterHit + ; +#pragma link C++ class AdvTargetHit + ; +#pragma link C++ class AdvMuFilterHit + ; +#pragma link C++ class sndScifiHit + ; #pragma link C++ class sndCluster; -#pragma link C++ class SNDLHCEventHeader+; -#pragma link C++ class sndRecoTrack+; -#pragma link C++ class Magnet+; -#pragma link C++ class MagnetPoint+; -#pragma link C++ class AdvTarget+; -#pragma link C++ class AdvTargetPoint+; -#pragma link C++ class AdvMuFilter+; -#pragma link C++ class AdvMuFilterPoint+; +#pragma link C++ class SNDLHCEventHeader + ; +#pragma link C++ class sndRecoTrack + ; +#pragma link C++ class Magnet + ; +#pragma link C++ class MagnetPoint + ; +#pragma link C++ class AdvTarget + ; +#pragma link C++ class AdvTargetPoint + ; +#pragma link C++ class AdvMuFilter + ; +#pragma link C++ class AdvMuFilterPoint + ; +#pragma link C++ class digitisation/ChargeDivision + ; +#pragma link C++ class digitisation/ChargeDrift +; +#pragma link C++ class digitisation/InducedCharge +; +#pragma link C++ class digitisation/FrontendDriver +; +#pragma link C++ class digitisation/StripNoise +; +#pragma link C++ class digitisation/SiG4UniversalFluctuation + ; +#pragma link C++ class digitisation/EnergyFluctUnit + ; +#pragma link C++ class digitisation/SurfaceSignal + ; +#pragma link C++ class digitisation/AdvDigitisation + ; #endif