49 #include <gsl/gsl_sf_erf.h>
51 #include "CLHEP/Random/RandGaussQ.h"
52 #include "CLHEP/Random/RandFlat.h"
93 using namespace sipixelobjects;
95 #define TP_DEBUG // protect all LogDebug with ifdef. Takes too much CPU
98 if(use_ineff_from_db_){
100 theSiPixelGainCalibrationService_->setESObjects( es );
104 fillLorentzAngle(es);
111 if(!use_deadmodule_DB_){
112 DeadModules = conf_.getParameter<
Parameters>(
"DeadModules");
121 if(!use_LorentzAngle_DB_){
123 tanLorentzAnglePerTesla_FPix=conf_.getParameter<
double>(
"TanLorentzAnglePerTesla_FPix");
124 tanLorentzAnglePerTesla_BPix=conf_.getParameter<
double>(
"TanLorentzAnglePerTesla_BPix");
142 conf_(conf) , fluctuate(0), theNoiser(0), pIndexConverter(0),
143 use_ineff_from_db_(conf_.getParameter<bool>(
"useDB")),
144 use_module_killing_(conf_.getParameter<bool>(
"killModules")),
145 use_deadmodule_DB_(conf_.getParameter<bool>(
"DeadModules_DB")),
146 use_LorentzAngle_DB_(conf_.getParameter<bool>(
"LorentzAngle_DB")),
147 theSiPixelGainCalibrationService_(0),rndEngine(eng)
238 if (thePixelLuminosity==-20){
241 thePixelColEfficiency[2] =
conf_.
getParameter<
double>(
"thePixelColEfficiency_BPix3");
242 thePixelColEfficiency[3] =
conf_.
getParameter<
double>(
"thePixelColEfficiency_FPix1");
243 thePixelColEfficiency[4] =
conf_.
getParameter<
double>(
"thePixelColEfficiency_FPix2");
244 cout<<
"\nReading in custom Pixel efficiencies "<<thePixelColEfficiency[0]<<
" , "<<thePixelColEfficiency[1]<<
" , "
245 <<thePixelColEfficiency[2]<<
" , "<<thePixelColEfficiency[3]<<
" , "<<thePixelColEfficiency[4]<<
"\n";
246 if (thePixelColEfficiency[0]<=0.5) {
cout <<
"\n\nDid you mean to set the Pixel efficiency at "<<thePixelColEfficiency[0]
247 <<
", or did you mean for this to be the inefficiency?\n\n\n";}
252 if(thePixelLuminosity==-1) {
254 for (
int i=0;
i<6;
i++) {
262 }
else if(thePixelLuminosity==0) {
265 for (
int i=0;
i<6;
i++) {
286 }
else if(thePixelLuminosity>0) {
289 for (
int i=0;
i<6;
i++) {
310 if(thePixelLuminosity==10) {
318 if(addNoise || thePixelLuminosity || fluctuateCharge || addThresholdSmearing ) {
327 if(fluctuateCharge) {
332 if(addThresholdSmearing) {
340 if(doMissCalibrate) {
342 <<
" miss-calibrate the pixel amplitude ";
344 const bool ReadCalParameters =
false;
345 if(ReadCalParameters) {
348 char filename[80] =
"phCalibrationFit_C0.dat";
350 in_file.open(filename,
ios::in );
352 cout <<
" File not found " << endl;
355 cout <<
" file opened : " << filename << endl;
358 for (
int i = 0;
i < 3;
i++) {
359 in_file.getline(line, 500,
'\n');
363 cout <<
" test map" << endl;
365 float par0,par1,par2,par3;
369 for(
int i=0;
i<(52*80);
i++) {
370 in_file >> par0 >> par1 >> par2 >> par3 >> name >> colid
373 cerr <<
"Cannot read data file" << endl;
376 if( in_file.eof() != 0 ) {
377 cerr << in_file.eof() <<
" " << in_file.gcount() <<
" "
378 << in_file.fail() <<
" " << in_file.good() <<
" end of file "
394 calmap.insert(pair<int,CalParameters>(chan,onePix));
398 if(rowid!=p.first)
cout<<
" wrong channel row "<<rowid<<
" "<<p.first<<endl;
399 if(colid!=p.second)
cout<<
" wrong channel col "<<colid<<
" "<<p.second<<endl;
424 LogInfo (
"PixelDigitizer ") <<
"SiPixelDigitizerAlgorithm constructed"
425 <<
"Configuration parameters:"
426 <<
"Threshold/Gain = "
427 <<
"threshold in electron FPix = "
428 << theThresholdInE_FPix
429 <<
"threshold in electron BPix = "
430 << theThresholdInE_BPix
431 <<
" " << theElectronPerADC <<
" " << theAdcFullScale
432 <<
" The delta cut-off is set to " << tMax
439 LogDebug (
"PixelDigitizer")<<
"SiPixelDigitizerAlgorithm deleted";
462 const std::vector<PSimHit> &
input,
482 vector<PixelDigi> collector =
digitize(pixdet);
487 LogDebug (
"PixelDigitizer") <<
"[SiPixelDigitizerAlgorithm] converted " << collector.size() <<
" PixelDigis in DetUnit" <<
detID;
544 <<
" PixelDigitizer "
551 vector<PSimHit>::const_iterator ssbegin;
556 << (*ssbegin).particleType() <<
" " << (*ssbegin).pabs() <<
" "
557 << (*ssbegin).energyLoss() <<
" " << (*ssbegin).tof() <<
" "
558 << (*ssbegin).trackId() <<
" " << (*ssbegin).processType() <<
" "
559 << (*ssbegin).detUnitId()
560 << (*ssbegin).entryPoint() <<
" " << (*ssbegin).exitPoint() ;
606 const float SegmentLength = 0.0010;
613 float length = direction.
mag();
635 float momentum = hit.
pabs();
667 delete[] elossVector;
675 float eloss,
float length,
676 int NumberOfSegs,
float elossVector[]) {
683 double particleMass = 139.6;
686 if(pid==11) particleMass = 0.511;
687 else if(pid==13) particleMass = 105.7;
688 else if(pid==321) particleMass = 493.7;
689 else if(pid==2212) particleMass = 938.3;
692 float segmentLength = length/NumberOfSegs;
697 double segmentEloss = (1000.*eloss)/NumberOfSegs;
698 for (
int i=0;
i<NumberOfSegs;
i++) {
704 double deltaCutoff =
tMax;
706 particleMass, deltaCutoff,
707 double(segmentLength*10.),
708 segmentEloss )/1000.;
716 float ratio = eloss/sum;
718 for (
int ii=0;ii<NumberOfSegs;ii++) elossVector[ii]= ratio*elossVector[ii];
720 float averageEloss = eloss/NumberOfSegs;
721 for (
int ii=0;ii<NumberOfSegs;ii++) elossVector[ii]= averageEloss;
733 LogDebug (
"Pixel Digitizer") <<
" enter drift " ;
745 if(driftDir.
z() ==0.) {
746 LogWarning(
"Magnetic field") <<
" pxlx: drift in z is zero ";
754 float TanLorenzAngleX, TanLorenzAngleY,dir_z, CosLorenzAngleX,
758 TanLorenzAngleX = driftDir.
x();
759 TanLorenzAngleY = driftDir.
y();
760 dir_z = driftDir.
z();
761 CosLorenzAngleX = 1./
sqrt(1.+TanLorenzAngleX*TanLorenzAngleX);
762 CosLorenzAngleY = 1./
sqrt(1.+TanLorenzAngleY*TanLorenzAngleY);
766 TanLorenzAngleX = driftDir.
x();
767 TanLorenzAngleY = 0.;
768 dir_z = driftDir.
z();
769 CosLorenzAngleX = 1./
sqrt(1.+TanLorenzAngleX*TanLorenzAngleX);
770 CosLorenzAngleY = 1.;
776 <<
" Lorentz Tan " << TanLorenzAngleX <<
" " << TanLorenzAngleY <<
" "
777 << CosLorenzAngleX <<
" " << CosLorenzAngleY <<
" "
790 float SegX, SegY, SegZ;
808 if( DriftDistance < 0.) {
814 float XDriftDueToMagField = DriftDistance * TanLorenzAngleX;
815 float YDriftDueToMagField = DriftDistance * TanLorenzAngleY;
818 float CloudCenterX = SegX + XDriftDueToMagField;
819 float CloudCenterY = SegY + YDriftDueToMagField;
822 DriftLength =
sqrt( DriftDistance*DriftDistance +
823 XDriftDueToMagField*XDriftDueToMagField +
824 YDriftDueToMagField*YDriftDueToMagField );
830 Sigma_x = Sigma / CosLorenzAngleX ;
831 Sigma_y = Sigma / CosLorenzAngleY ;
852 <<
" enter induce_signal, "
857 typedef map< int, float, less<int> > hit_map_type;
858 hit_map_type hit_signal;
861 map<int, float, less<int> >
x,
y;
869 float CloudCenterX =
i->position().x();
870 float CloudCenterY =
i->position().y();
871 float SigmaX =
i->sigma_x();
872 float SigmaY =
i->sigma_y();
873 float Charge =
i->amplitude();
884 <<
" cloud " <<
i->position().x() <<
" " <<
i->position().y() <<
" "
885 <<
i->sigma_x() <<
" " <<
i->sigma_y() <<
" " <<
i->amplitude();
908 int IPixRightUpX = int( floor( mp.
x()));
909 int IPixRightUpY = int( floor( mp.
y()));
912 LogDebug (
"Pixel Digitizer") <<
" right-up " << PointRightUp <<
" "
913 << mp.
x() <<
" " << mp.
y() <<
" "
914 << IPixRightUpX <<
" " << IPixRightUpY ;
919 int IPixLeftDownX = int( floor( mp.
x()));
920 int IPixLeftDownY = int( floor( mp.
y()));
923 LogDebug (
"Pixel Digitizer") <<
" left-down " << PointLeftDown <<
" "
924 << mp.
x() <<
" " << mp.
y() <<
" "
925 << IPixLeftDownX <<
" " << IPixLeftDownY ;
931 IPixLeftDownX = 0<IPixLeftDownX ? IPixLeftDownX : 0 ;
932 IPixLeftDownY = 0<IPixLeftDownY ? IPixLeftDownY : 0 ;
939 for (ix=IPixLeftDownX; ix<=IPixRightUpX; ix++) {
940 float xUB, xLB, UpperBound, LowerBound;
945 if(ix == 0 || SigmaX==0. )
951 int status = gsl_sf_erf_Q_e( (xLB-CloudCenterX)/SigmaX, &result);
953 LogWarning (
"Integration")<<
"could not compute gaussian probability";
954 LowerBound = 1-result.val;
957 if(ix ==
numRows-1 || SigmaX==0. )
963 int status = gsl_sf_erf_Q_e( (xUB-CloudCenterX)/SigmaX, &result);
965 LogWarning (
"Integration")<<
"could not compute gaussian probability";
966 UpperBound = 1. - result.val;
969 float TotalIntegrationRange = UpperBound - LowerBound;
970 x[ix] = TotalIntegrationRange;
978 for (iy=IPixLeftDownY; iy<=IPixRightUpY; iy++) {
979 float yUB, yLB, UpperBound, LowerBound;
981 if(iy == 0 || SigmaY==0.)
987 int status = gsl_sf_erf_Q_e( (yLB-CloudCenterY)/SigmaY, &result);
989 LogWarning (
"Integration")<<
"could not compute gaussian probability";
990 LowerBound = 1. - result.val;
999 int status = gsl_sf_erf_Q_e( (yUB-CloudCenterY)/SigmaY, &result);
1001 LogWarning (
"Integration")<<
"could not compute gaussian probability";
1002 UpperBound = 1. - result.val;
1005 float TotalIntegrationRange = UpperBound - LowerBound;
1006 y[iy] = TotalIntegrationRange;
1013 for (ix=IPixLeftDownX; ix<=IPixRightUpX; ix++) {
1014 for (iy=IPixLeftDownY; iy<=IPixRightUpY; iy++) {
1016 float ChargeFraction = Charge*x[ix]*y[iy];
1018 if( ChargeFraction > 0. ) {
1021 hit_signal[chan] += ChargeFraction;
1031 <<
" pixel " << ix <<
" " << iy <<
" - "<<
" "
1032 << chan <<
" " << ChargeFraction<<
" "
1033 << mp.
x() <<
" " << mp.
y() <<
" "
1034 << lp.
x() <<
" " << lp.
y() <<
" "
1064 for ( hit_map_type::const_iterator im = hit_signal.begin();
1065 im != hit_signal.end(); im++) {
1068 int chan = (*im).first;
1073 <<
" pixel " << ip.first <<
" " << ip.second <<
" "
1087 LogDebug (
"Pixel Digitizer") <<
" make digis "<<
" "
1090 <<
" List pixels passing threshold ";
1097 float signalInElectrons = (*i).second ;
1106 int chan = (*i).first;
1113 int col = ip.second;
1122 << (*i).first <<
" " << (*i).second <<
" " << signalInElectrons
1123 <<
" " << adc << ip.first <<
" " << ip.second ;
1130 if((*i).second.hits().size()>0){
1133 for( vector<const PSimHit*>::const_iterator ihit = (*i).second.hits().begin();
1134 ihit != (*i).second.hits().end(); ihit++) {
1135 simi[(**ihit).trackId()].push_back((*i).second.individualampl()[il]);
1140 for( simlink_map::iterator simiiter=
simi.begin();
1141 simiiter!=
simi.end();
1144 float sum_samechannel=0;
1145 for (
unsigned int iii=0;iii<(*simiiter).second.size();iii++){
1146 sum_samechannel+=(*simiiter).second[iii];
1148 float fraction=sum_samechannel/(*i).second;
1149 if(fraction>1.) fraction=1.;
1174 if((*i).second < 3000)
1177 }
else if((*i).second < 6000){
1188 if(((*i).second +
Amplitude(noise+noise_ChargeVCALSmearing,0,-1.)) < 0. ) {
1189 (*i).second.set(0);}
1191 (*i).second +=
Amplitude(noise+noise_ChargeVCALSmearing,0,-1.);
1201 if(((*i).second +
Amplitude(noise,0,-1.)) < 0. ) {
1202 (*i).second.set(0);}
1216 map<int,float, less<int> > otherPixels;
1217 map<int,float, less<int> >::iterator mapI;
1226 <<
" Add noisy pixels " <<
numRows <<
" "
1229 << otherPixels.size() ;
1233 for (mapI = otherPixels.begin(); mapI!= otherPixels.end(); mapI++) {
1234 int iy = ((*mapI).first) /
numRows;
1235 int ix = ((*mapI).first) - (iy*
numRows);
1239 LogWarning (
"Pixel Geometry") <<
" error in iy " << iy ;
1240 if( ix < 0 || ix > (numRows-1) )
1241 LogWarning (
"Pixel Geometry") <<
" error in ix " << ix ;
1247 <<
" Storing noise = " << (*mapI).first <<
" " << (*mapI).second
1248 <<
" " << ix <<
" " << iy <<
" " << chan ;
1253 int noise=int( (*mapI).second );
1268 float pixelEfficiency = 1.0;
1269 float columnEfficiency = 1.0;
1270 float chipEfficiency = 1.0;
1301 LogDebug (
"Pixel Digitizer") <<
" enter pixel_inefficiency " << pixelEfficiency <<
" "
1302 << columnEfficiency <<
" " << chipEfficiency;
1310 map<int, int, less<int> >chips, columns;
1311 map<int, int, less<int> >::iterator iter;
1317 int chan =
i->first;
1320 int col = ip.second;
1328 columns[dColInDet]++;
1332 for ( iter = chips.begin(); iter != chips.end() ; iter++ ) {
1335 if( rand > chipEfficiency ) chips[iter->first]=0;
1339 for ( iter = columns.begin(); iter != columns.end() ; iter++ ) {
1342 if( rand > columnEfficiency ) columns[iter->first]=0;
1352 int col = ip.second;
1362 if( chips[chipIndex]==0 || columns[dColInDet]==0
1363 || rand>pixelEfficiency ) {
1382 const float signalInElectrons)
const {
1421 newAmp = p3 + p2 * tanh(p0*signal - p1);
1500 dir_z = -(1 + alpha2_BPix* Bfield.z()*Bfield.z() );
1501 scale = (1 + alpha2_BPix* Bfield.z()*Bfield.z() );
1506 dir_z = -(1 + alpha2_FPix* Bfield.z()*Bfield.z() );
1507 scale = (1 + alpha2_FPix* Bfield.z()*Bfield.z() );
1515 dir_x = -( lorentzAngle * Bfield.y() + alpha2 * Bfield.z()* Bfield.x() );
1516 dir_y = +( lorentzAngle * Bfield.x() - alpha2 * Bfield.z()* Bfield.y() );
1517 dir_z = -(1 + alpha2 * Bfield.z()*Bfield.z() );
1518 scale = (1 + alpha2 * Bfield.z()*Bfield.z() );
1524 LogDebug (
"Pixel Digitizer") <<
" The drift direction in local coordinate is "
1525 << theDriftDirection ;
1528 return theDriftDirection;
1543 int col = ip.second;
1564 Parameters::iterator itDeadModules=
DeadModules.begin();
1566 for(; itDeadModules !=
DeadModules.end(); ++itDeadModules){
1567 int Dead_detID = itDeadModules->getParameter<
int>(
"Dead_detID");
1568 if(detid==Dead_detID){
1577 std::string
Module = itDeadModules->getParameter<std::string>(
"Module");
1579 if(Module==
"whole"){
1588 if(Module==
"tbmA" && ip.first>=80 && ip.first<=159){
1592 if( Module==
"tbmB" && ip.first<=79){
1606 std::vector<SiPixelQuality::disabledModuleType>disabledModules =
SiPixelBadModule_->getBadComponentList();
1610 for (
size_t id=0;
id<disabledModules.size();
id++)
1612 if(detid==disabledModules[
id].DetID){
1614 badmodule = disabledModules[
id];
1632 std::vector<GlobalPixel> badrocpositions (0);
1633 for(
unsigned int j = 0;
j < 16;
j++){
1636 std::vector<CablingPathToDetUnit>
path =
map_.
product()->pathToDetUnit(detid);
1637 typedef std::vector<CablingPathToDetUnit>::const_iterator
IT;
1638 for (IT it = path.begin(); it != path.end(); ++it) {
1643 badrocpositions.push_back(global);
1654 for(std::vector<GlobalPixel>::const_iterator it = badrocpositions.begin(); it != badrocpositions.end(); ++it){
1655 if(it->row >= 80 && ip.first >= 80 ){
1656 if((fabs(ip.second - it->col) < 26) ) {
i->second.set(0.);}
1657 else if(it->row==120 && ip.second-it->col==26){
i->second.set(0.);}
1658 else if(it->row==119 && it->col-ip.second==26){
i->second.set(0.);}
1660 else if(it->row < 80 && ip.first < 80 ){
1661 if((fabs(ip.second - it->col) < 26) ){
i->second.set(0.);}
1662 else if(it->row==40 && ip.second-it->col==26){
i->second.set(0.);}
1663 else if(it->row==39 && it->col-ip.second==26){
i->second.set(0.);}
int adc(sample_type sample)
get the ADC sample (12 bits)
void init(const edm::EventSetup &es)
std::vector< T > collection_type
T getParameter(std::string const &) const
float thePixelThresholdInE
T getUntrackedParameter(std::string const &, T const &) const
bool use_LorentzAngle_DB_
Local3DVector LocalVector
float tof() const
deprecated name for timeOfFlight()
edm::ESHandle< SiPixelFedCablingMap > map_
GaussianTailNoiseGenerator * theNoiser
float tanLorentzAnglePerTesla_BPix
GlobalPixel toGlobal(const LocalPixel &loc) const
static int DColumn(const int colROC)
virtual LocalPoint localPosition(const MeasurementPoint &) const =0
virtual int ncolumns() const =0
const PixelGeomDetUnit * _detp
SiG4UniversalFluctuation * fluctuate
float thePixelColEfficiency[6]
float thePixelChipEfficiency[6]
virtual int nrows() const =0
int transformToROC(const int col, const int row, int &rocId, int &colROC, int &rowROC) const
std::vector< PixelDigi > digitize(PixelGeomDetUnit *det)
~SiPixelDigitizerAlgorithm()
void fluctuateEloss(int particleId, float momentum, float eloss, float length, int NumberOfSegments, float elossVector[])
identify pixel inside single ROC
unsigned int layer() const
layer id
const BoundPlane & specificSurface() const
Same as surface(), kept for backward compatibility.
SiPixelGainCalibrationOfflineSimService * theSiPixelGainCalibrationService_
static int pixelToChannel(int row, int col)
global coordinates (row and column in DetUnit, as in PixelDigi)
uint32_t rawId() const
get the raw id
CLHEP::RandGaussQ * gaussDistributionVCALNoise_
virtual float thickness() const =0
double SampleFluctuations(const double momentum, const double mass, double &tmax, const double length, const double meanLoss)
double theThresholdSmearing_BPix
Measurement2DPoint MeasurementPoint
Measurement points are two-dimensional by default.
Local3DPoint exitPoint() const
Exit point in the local Det frame.
float theThresholdInE_BPix
void fillDeadModules(const edm::EventSetup &es)
float thePixelEfficiency[6]
static int pixelToChannelROC(const int rowROC, const int colROC)
std::vector< PixelDigiSimLink > link_coll
edm::DetSet< PixelDigi >::collection_type run(const std::vector< PSimHit > &input, PixelGeomDetUnit *pixdet, GlobalVector)
static std::pair< int, int > channelToPixelROC(const int chan)
unsigned int idInDetUnit() const
id of this ROC in DetUnit etermined by token path
float tanLorentzAnglePerTesla_FPix
PixelIndices * pIndexConverter
virtual int channel(const LocalPoint &p) const =0
DetId geographicalId() const
The label of this GeomDet.
virtual MeasurementPoint measurementPosition(const LocalPoint &) const =0
void fillLorentzAngle(const edm::EventSetup &es)
std::vector< LinkConnSpec >::const_iterator IT
float theThresholdInE_FPix
float pabs() const
fast and more accurate access to momentumAtEntry().mag()
signal_map_type::iterator signal_map_iterator
int subdetId() const
get the contents of the subdetector field (not cast into any detector's numbering enum) ...
std::vector< PSimHit > _PixelHits
float theNoiseInElectrons
void module_killing_conf()
edm::ESHandle< SiPixelLorentzAngle > SiPixelLorentzAngle_
std::map< int, CalParameters, std::less< int > > calmap
double electronsPerVCAL_Offset
CLHEP::RandGaussQ * gaussDistribution_
Double_t Sigma(double mHreq, TString proc)
virtual std::pair< float, float > pitch() const =0
const PixelTopology * topol
const Bounds & bounds() const
std::vector< PixelDigi > internal_coll
edm::ESHandle< SiPixelQuality > SiPixelBadModule_
void fillMapandGeom(const edm::EventSetup &es)
static std::pair< int, int > channelToPixel(int ch)
void pixel_inefficiency()
virtual const PixelTopology & specificTopology() const
Returns a reference to the pixel proxy topology.
T const * product() const
bool addChargeVCALSmearing
row and collumn in ROC representation
CLHEP::RandGaussQ * smearedThreshold_FPix_
double theThresholdSmearing_FPix
std::vector< EnergyDepositUnit > _ionization_points
std::vector< edm::ParameterSet > Parameters
CLHEP::RandGaussQ * smearedThreshold_BPix_
void generate(int NumberOfchannels, float threshold, float noiseRMS, std::map< int, float > &theMap)
float energyLoss() const
The energy deposit in the PSimHit, in ???.
const BoundPlane & surface() const
The nominal surface of the GeomDet.
void drift(const PSimHit &hit)
float theSmearedChargeRMS
bool addThresholdSmearing
const RotationType & rotation() const
void primary_ionization(const PSimHit &hit)
CLHEP::RandFlat * flatDistribution_
SiPixelDigitizerAlgorithm(const edm::ParameterSet &conf, CLHEP::HepRandomEngine &)
std::vector< SignalPoint > _collection_points
float missCalibrate(int col, int row, float amp) const
static int DColumnInModule(const int dcol, const int chipIndex)
void induce_signal(const PSimHit &hit)
const PositionType & position() const
Local3DPoint entryPoint() const
Entry point in the local Det frame.
bool isDead(const uint32_t &detID, const int &col, const int &row)
void pixel_inefficiency_db()
*vegas h *****************************************************used in the default bin number in original ***version of VEGAS is ***a higher bin number might help to derive a more precise ***grade subtle point
unsigned int detUnitId() const
LocalVector DriftDirection()
CLHEP::HepRandomEngine & rndEngine