00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041 #include <vector>
00042 #include <iostream>
00043
00044 #include "DataFormats/SiPixelDigi/interface/PixelDigi.h"
00045 #include "SimDataFormats/TrackingHit/interface/PSimHitContainer.h"
00046 #include "SimTracker/SiPixelDigitizer/interface/SiPixelDigitizerAlgorithm.h"
00047
00048 #include <gsl/gsl_sf_erf.h>
00049 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
00050 #include "CLHEP/Random/RandGaussQ.h"
00051 #include "CLHEP/Random/RandFlat.h"
00052
00053
00054 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
00055 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00056
00057 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00058 #include "FWCore/ServiceRegistry/interface/Service.h"
00059 #include "FWCore/Utilities/interface/Exception.h"
00060 #include "CalibTracker/SiPixelESProducers/interface/SiPixelGainCalibrationOfflineSimService.h"
00061
00062
00063 #include "DataFormats/DetId/interface/DetId.h"
00064
00065 #include "CondFormats/SiPixelObjects/interface/GlobalPixel.h"
00066 #include "FWCore/Framework/interface/ESHandle.h"
00067 #include "CondFormats/DataRecord/interface/SiPixelFedCablingMapRcd.h"
00068 #include "CondFormats/SiPixelObjects/interface/SiPixelFedCablingMap.h"
00069 #include "CondFormats/SiPixelObjects/interface/SiPixelFedCablingTree.h"
00070 #include "CondFormats/SiPixelObjects/interface/SiPixelFedCabling.h"
00071 #include "CondFormats/SiPixelObjects/interface/PixelROC.h"
00072 #include "CondFormats/SiPixelObjects/interface/LocalPixel.h"
00073 #include "CondFormats/SiPixelObjects/interface/CablingPathToDetUnit.h"
00074
00075 #include "CondFormats/SiPixelObjects/interface/SiPixelFrameReverter.h"
00076 #include "CondFormats/SiPixelObjects/interface/PixelFEDCabling.h"
00077 #include "CondFormats/SiPixelObjects/interface/PixelFEDLink.h"
00078 #include "DataFormats/FEDRawData/interface/FEDNumbering.h"
00079 #include "DataFormats/SiPixelDetId/interface/PixelBarrelName.h"
00080 #include "DataFormats/SiPixelDetId/interface/PixelEndcapName.h"
00081
00082
00083 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00084 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00085 #include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h"
00086 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
00087
00088 #include "CondFormats/SiPixelObjects/interface/PixelROC.h"
00089
00090 using namespace std;
00091 using namespace edm;
00092 using namespace sipixelobjects;
00093
00094 #define TP_DEBUG // protect all LogDebug with ifdef. Takes too much CPU
00095
00096 void SiPixelDigitizerAlgorithm::init(const edm::EventSetup& es){
00097 if(use_ineff_from_db_){
00098 theSiPixelGainCalibrationService_= new SiPixelGainCalibrationOfflineSimService(conf_);
00099 theSiPixelGainCalibrationService_->setESObjects( es );
00100 }
00101
00102 fillDeadModules(es);
00103 fillLorentzAngle(es);
00104 fillMapandGeom(es);
00105 }
00106
00107
00108
00109 void SiPixelDigitizerAlgorithm::fillDeadModules(const edm::EventSetup& es){
00110 if(!use_deadmodule_DB_){
00111 DeadModules = conf_.getParameter<Parameters>("DeadModules");
00112 }
00113 else{
00114
00115 es.get<SiPixelQualityRcd>().get(SiPixelBadModule_);
00116 }
00117 }
00118
00119 void SiPixelDigitizerAlgorithm::fillLorentzAngle(const edm::EventSetup& es){
00120 if(!use_LorentzAngle_DB_){
00121
00122 tanLorentzAnglePerTesla_FPix=conf_.getParameter<double>("TanLorentzAnglePerTesla_FPix");
00123 tanLorentzAnglePerTesla_BPix=conf_.getParameter<double>("TanLorentzAnglePerTesla_BPix");
00124 }
00125 else {
00126
00127
00128 es.get<SiPixelLorentzAngleSimRcd>().get(SiPixelLorentzAngle_);
00129 }
00130 }
00131
00132
00133 void SiPixelDigitizerAlgorithm::fillMapandGeom(const edm::EventSetup& es){
00134 es.get<SiPixelFedCablingMapRcd>().get(map_);
00135 es.get<TrackerDigiGeometryRecord>().get(geom_);
00136 }
00137
00138
00139
00140 SiPixelDigitizerAlgorithm::SiPixelDigitizerAlgorithm(const edm::ParameterSet& conf, CLHEP::HepRandomEngine& eng) :
00141 conf_(conf) , fluctuate(0), theNoiser(0), pIndexConverter(0),
00142 use_ineff_from_db_(conf_.getParameter<bool>("useDB")),
00143 use_module_killing_(conf_.getParameter<bool>("killModules")),
00144 use_deadmodule_DB_(conf_.getParameter<bool>("DeadModules_DB")),
00145 use_LorentzAngle_DB_(conf_.getParameter<bool>("LorentzAngle_DB")),
00146 theSiPixelGainCalibrationService_(0),rndEngine(eng)
00147 {
00148 using std::cout;
00149 using std::endl;
00150
00151
00152
00153 NumberOfSegments = 20;
00154 ClusterWidth = 3.;
00155 GeVperElectron = 3.61E-09;
00156 Sigma0 = 0.00037;
00157 Dist300 = 0.0300;
00158
00159 alpha2Order = conf_.getParameter<bool>("Alpha2Order");
00160
00161
00162
00163
00164
00165
00166
00167 theElectronPerADC=conf_.getParameter<double>("ElectronPerAdc");
00168
00169
00170
00171 theAdcFullScale=conf_.getParameter<int>("AdcFullScale");
00172
00173
00174
00175
00176 theThresholdInE_FPix=conf_.getParameter<double>("ThresholdInElectrons_FPix");
00177 theThresholdInE_BPix=conf_.getParameter<double>("ThresholdInElectrons_BPix");
00178
00179
00180 addThresholdSmearing = conf_.getParameter<bool>("AddThresholdSmearing");
00181 theThresholdSmearing_FPix = conf_.getParameter<double>("ThresholdSmearing_FPix");
00182 theThresholdSmearing_BPix = conf_.getParameter<double>("ThresholdSmearing_BPix");
00183
00184
00185 FPix_p0 = conf_.getParameter<double>("FPix_SignalResponse_p0");
00186 FPix_p1 = conf_.getParameter<double>("FPix_SignalResponse_p1");
00187 FPix_p2 = conf_.getParameter<double>("FPix_SignalResponse_p2");
00188 FPix_p3 = conf_.getParameter<double>("FPix_SignalResponse_p3");
00189
00190 BPix_p0 = conf_.getParameter<double>("BPix_SignalResponse_p0");
00191 BPix_p1 = conf_.getParameter<double>("BPix_SignalResponse_p1");
00192 BPix_p2 = conf_.getParameter<double>("BPix_SignalResponse_p2");
00193 BPix_p3 = conf_.getParameter<double>("BPix_SignalResponse_p3");
00194
00195
00196 electronsPerVCAL = conf_.getParameter<double>("ElectronsPerVcal");
00197 electronsPerVCAL_Offset = conf_.getParameter<double>("ElectronsPerVcal_Offset");
00198
00199
00200 addNoise=conf_.getParameter<bool>("AddNoise");
00201
00202
00203
00204 addChargeVCALSmearing=conf_.getParameter<bool>("ChargeVCALSmearing");
00205
00206
00207 addNoisyPixels=conf_.getParameter<bool>("AddNoisyPixels");
00208
00209
00210 theNoiseInElectrons=conf_.getParameter<double>("NoiseInElectrons");
00211
00212
00213 theReadoutNoise=conf_.getParameter<double>("ReadoutNoiseInElec");
00214
00215
00216
00217 theTofLowerCut=conf_.getParameter<double>("TofLowerCut");
00218 theTofUpperCut=conf_.getParameter<double>("TofUpperCut");
00219
00220
00221 fluctuateCharge=conf_.getUntrackedParameter<bool>("FluctuateCharge",true);
00222
00223
00224
00225
00226 tMax =conf_.getParameter<double>("DeltaProductionCut");
00227
00228
00229 thePixelLuminosity=conf_.getParameter<int>("AddPixelInefficiency");
00230
00231
00232 doMissCalibrate=conf_.getParameter<bool>("MissCalibrate");
00233 theGainSmearing=conf_.getParameter<double>("GainSmearing");
00234 theOffsetSmearing=conf_.getParameter<double>("OffsetSmearing");
00235
00236
00237
00238
00239
00240 if (thePixelLuminosity==-1) {
00241 pixelInefficiency=false;
00242 for (int i=0; i<6;i++) {
00243 thePixelEfficiency[i] = 1.;
00244 thePixelColEfficiency[i] = 1.;
00245 thePixelChipEfficiency[i] = 1.;
00246 }
00247
00248
00249
00250 } else if (thePixelLuminosity==0) {
00251 pixelInefficiency=true;
00252
00253 for (int i=0; i<6;i++) {
00254 if(i<3) {
00255
00256
00257 thePixelEfficiency[i] = 1.-0.001;
00258
00259 thePixelColEfficiency[i] = 1.-0.001;
00260
00261 thePixelChipEfficiency[i] = 1.-0.001;
00262 } else {
00263
00264
00265 thePixelEfficiency[i] = 1.-0.001;
00266
00267 thePixelColEfficiency[i] = 1.-0.001;
00268
00269 thePixelChipEfficiency[i] = 1.-0.001;
00270 }
00271 }
00272
00273
00274 } else if (thePixelLuminosity>0) {
00275 pixelInefficiency=true;
00276
00277 for (int i=0; i<6;i++) {
00278 if(i<3) {
00279
00280
00281 thePixelEfficiency[i] = 1.-0.01;
00282
00283 thePixelColEfficiency[i] = 1.-0.01;
00284
00285 thePixelChipEfficiency[i] = 1.-0.0025;
00286 } else {
00287
00288
00289 thePixelEfficiency[i] = 1.-0.01;
00290
00291 thePixelColEfficiency[i] = 1.-0.01;
00292
00293 thePixelChipEfficiency[i] = 1.-0.0025;
00294 }
00295 }
00296
00297
00298 if(thePixelLuminosity==10) {
00299 thePixelColEfficiency[0] = 1.-0.034;
00300 thePixelEfficiency[0] = 1.-0.015;
00301 }
00302
00303 }
00304
00305
00306 if(addNoise || thePixelLuminosity || fluctuateCharge || addThresholdSmearing ) {
00307 gaussDistribution_ = new CLHEP::RandGaussQ(rndEngine, 0., theReadoutNoise);
00308 gaussDistributionVCALNoise_ = new CLHEP::RandGaussQ(rndEngine, 0., 1.);
00309 flatDistribution_ = new CLHEP::RandFlat(rndEngine, 0., 1.);
00310
00311 if(addNoise) {
00312 theNoiser = new GaussianTailNoiseGenerator(rndEngine);
00313 }
00314
00315 if(fluctuateCharge) {
00316 fluctuate = new SiG4UniversalFluctuation(rndEngine);
00317 }
00318
00319
00320 if(addThresholdSmearing) {
00321 smearedThreshold_FPix_ = new CLHEP::RandGaussQ(rndEngine, theThresholdInE_FPix , theThresholdSmearing_FPix);
00322 smearedThreshold_BPix_ = new CLHEP::RandGaussQ(rndEngine, theThresholdInE_BPix , theThresholdSmearing_BPix);
00323 }
00324 }
00325
00326
00327
00328 if(doMissCalibrate) {
00329 LogDebug ("PixelDigitizer ")
00330 << " miss-calibrate the pixel amplitude ";
00331
00332 const bool ReadCalParameters = false;
00333 if(ReadCalParameters) {
00334
00335 ifstream in_file;
00336 char filename[80] = "phCalibrationFit_C0.dat";
00337
00338 in_file.open(filename, ios::in );
00339 if (in_file.bad()) {
00340 cout << " File not found " << endl;
00341 return;
00342 }
00343 cout << " file opened : " << filename << endl;
00344
00345 char line[500];
00346 for (int i = 0; i < 3; i++) {
00347 in_file.getline(line, 500,'\n');
00348 cout<<line<<endl;
00349 }
00350
00351 cout << " test map" << endl;
00352
00353 float par0,par1,par2,par3;
00354 int colid,rowid;
00355 string name;
00356
00357 for(int i=0;i<(52*80);i++) {
00358 in_file >> par0 >> par1 >> par2 >> par3 >> name >> colid
00359 >> rowid;
00360 if (in_file.bad()) {
00361 cerr << "Cannot read data file" << endl;
00362 return;
00363 }
00364 if ( in_file.eof() != 0 ) {
00365 cerr << in_file.eof() << " " << in_file.gcount() << " "
00366 << in_file.fail() << " " << in_file.good() << " end of file "
00367 << endl;
00368 return;
00369 }
00370
00371
00372
00373
00374 CalParameters onePix;
00375 onePix.p0=par0;
00376 onePix.p1=par1;
00377 onePix.p2=par2;
00378 onePix.p3=par3;
00379
00380
00381 int chan = PixelIndices::pixelToChannelROC(rowid,colid);
00382 calmap.insert(pair<int,CalParameters>(chan,onePix));
00383
00384
00385 pair<int,int> p = PixelIndices::channelToPixelROC(chan);
00386 if(rowid!=p.first) cout<<" wrong channel row "<<rowid<<" "<<p.first<<endl;
00387 if(colid!=p.second) cout<<" wrong channel col "<<colid<<" "<<p.second<<endl;
00388
00389 }
00390
00391 cout << " map size " << calmap.size() <<" max "<<calmap.max_size() << " "
00392 <<calmap.empty()<< endl;
00393
00394
00395
00396
00397
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409 }
00410 }
00411
00412 LogInfo ("PixelDigitizer ") <<"SiPixelDigitizerAlgorithm constructed"
00413 <<"Configuration parameters:"
00414 << "Threshold/Gain = "
00415 << "threshold in electron FPix = "
00416 << theThresholdInE_FPix
00417 << "threshold in electron BPix = "
00418 << theThresholdInE_BPix
00419 <<" " << theElectronPerADC << " " << theAdcFullScale
00420 << " The delta cut-off is set to " << tMax
00421 << " pix-inefficiency "<<thePixelLuminosity;
00422
00423 }
00424
00425 SiPixelDigitizerAlgorithm::~SiPixelDigitizerAlgorithm() {
00426
00427 LogDebug ("PixelDigitizer")<<"SiPixelDigitizerAlgorithm deleted";
00428
00429
00430 delete gaussDistribution_;
00431 delete gaussDistributionVCALNoise_;
00432 delete flatDistribution_;
00433
00434 delete theSiPixelGainCalibrationService_;
00435
00436
00437 if(addThresholdSmearing) {
00438 delete smearedThreshold_FPix_;
00439 delete smearedThreshold_BPix_;
00440 }
00441
00442 if(addNoise) delete theNoiser;
00443 if(fluctuateCharge) delete fluctuate;
00444
00445 }
00446
00447
00448 edm::DetSet<PixelDigi>::collection_type
00449 SiPixelDigitizerAlgorithm::run(
00450 const std::vector<PSimHit> &input,
00451 PixelGeomDetUnit *pixdet,
00452 GlobalVector bfield) {
00453
00454 _detp = pixdet;
00455 _PixelHits=input;
00456 _bfield=bfield;
00457
00458
00459
00460
00461
00462 detID= _detp->geographicalId().rawId();
00463
00464 _signal.clear();
00465
00466
00467 link_coll.clear();
00468
00469
00470 vector<PixelDigi> collector =digitize(pixdet);
00471
00472
00473
00474 #ifdef TP_DEBUG
00475 LogDebug ("PixelDigitizer") << "[SiPixelDigitizerAlgorithm] converted " << collector.size() << " PixelDigis in DetUnit" << detID;
00476 #endif
00477
00478 return collector;
00479 }
00480
00481
00482 vector<PixelDigi> SiPixelDigitizerAlgorithm::digitize(PixelGeomDetUnit *det){
00483
00484 if( _PixelHits.size() > 0 || addNoisyPixels) {
00485
00486 topol=&det->specificTopology();
00487 numColumns = topol->ncolumns();
00488 numRows = topol->nrows();
00489
00490
00491 moduleThickness = det->specificSurface().bounds().thickness();
00492
00493
00494
00495 if((pixelInefficiency>0) || doMissCalibrate ) {
00496 pIndexConverter = new PixelIndices(numColumns,numRows);
00497 }
00498
00499
00500
00501
00502
00503 unsigned int Sub_detid=DetId(detID).subdetId();
00504
00505 if(theNoiseInElectrons>0.){
00506 if(Sub_detid == PixelSubdetector::PixelBarrel){
00507 if(addThresholdSmearing) {
00508 thePixelThresholdInE = smearedThreshold_BPix_->fire();
00509 } else {
00510 thePixelThresholdInE = theThresholdInE_BPix;
00511 }
00512
00513 thePixelThreshold = thePixelThresholdInE/theNoiseInElectrons;
00514
00515 } else {
00516 if(addThresholdSmearing) {
00517 thePixelThresholdInE = smearedThreshold_FPix_->fire();
00518 } else {
00519 thePixelThresholdInE = theThresholdInE_FPix;
00520 }
00521
00522 thePixelThreshold = thePixelThresholdInE/theNoiseInElectrons;
00523
00524 }
00525 } else {
00526 thePixelThreshold = 0.;
00527 }
00528
00529
00530 #ifdef TP_DEBUG
00531 LogDebug ("PixelDigitizer")
00532 << " PixelDigitizer "
00533 << numColumns << " " << numRows << " " << moduleThickness;
00534 #endif
00535
00536
00537
00538
00539 vector<PSimHit>::const_iterator ssbegin;
00540 for (ssbegin= _PixelHits.begin();ssbegin !=_PixelHits.end(); ++ssbegin) {
00541
00542 #ifdef TP_DEBUG
00543 LogDebug ("Pixel Digitizer")
00544 << (*ssbegin).particleType() << " " << (*ssbegin).pabs() << " "
00545 << (*ssbegin).energyLoss() << " " << (*ssbegin).tof() << " "
00546 << (*ssbegin).trackId() << " " << (*ssbegin).processType() << " "
00547 << (*ssbegin).detUnitId()
00548 << (*ssbegin).entryPoint() << " " << (*ssbegin).exitPoint() ;
00549 #endif
00550
00551 _collection_points.clear();
00552
00553
00554
00555 if ( ((*ssbegin).tof() >= theTofLowerCut) && ((*ssbegin).tof() <= theTofUpperCut) ) {
00556 primary_ionization(*ssbegin);
00557 drift(*ssbegin);
00558
00559 induce_signal(*ssbegin);
00560 }
00561 }
00562
00563 if(addNoise) add_noise();
00564
00565
00566 if((pixelInefficiency>0) && (_signal.size()>0))
00567 pixel_inefficiency();
00568
00569 if(use_ineff_from_db_ && (_signal.size()>0))
00570 pixel_inefficiency_db();
00571
00572 delete pIndexConverter;
00573
00574 if(use_module_killing_ && use_deadmodule_DB_)
00575 module_killing_DB();
00576
00577 if(use_module_killing_ && !use_deadmodule_DB_)
00578 module_killing_conf();
00579
00580 }
00581
00582 make_digis();
00583
00584 return internal_coll;
00585 }
00586
00587
00588
00589
00590 void SiPixelDigitizerAlgorithm::primary_ionization(const PSimHit& hit) {
00591
00592
00593
00594 const float SegmentLength = 0.0010;
00595 float energy;
00596
00597
00598 LocalVector direction = hit.exitPoint() - hit.entryPoint();
00599
00600 float eLoss = hit.energyLoss();
00601 float length = direction.mag();
00602
00603 NumberOfSegments = int ( length / SegmentLength);
00604 if(NumberOfSegments < 1) NumberOfSegments = 1;
00605
00606 #ifdef TP_DEBUG
00607 LogDebug ("Pixel Digitizer")
00608 << " enter primary_ionzation " << NumberOfSegments
00609 << " shift = "
00610 << (hit.exitPoint().x()-hit.entryPoint().x()) << " "
00611 << (hit.exitPoint().y()-hit.entryPoint().y()) << " "
00612 << (hit.exitPoint().z()-hit.entryPoint().z()) << " "
00613 << hit.particleType() <<" "<< hit.pabs() ;
00614 #endif
00615
00616 float* elossVector = new float[NumberOfSegments];
00617
00618 if( fluctuateCharge ) {
00619
00620 int pid = hit.particleType();
00621
00622
00623 float momentum = hit.pabs();
00624
00625 fluctuateEloss(pid, momentum, eLoss, length, NumberOfSegments,
00626 elossVector);
00627 }
00628
00629 _ionization_points.resize( NumberOfSegments);
00630
00631
00632 for ( int i = 0; i != NumberOfSegments; i++) {
00633
00634 Local3DPoint point = hit.entryPoint() +
00635 float((i+0.5)/NumberOfSegments) * direction;
00636
00637 if( fluctuateCharge )
00638 energy = elossVector[i]/GeVperElectron;
00639 else
00640 energy = hit.energyLoss()/GeVperElectron/float(NumberOfSegments);
00641
00642 EnergyDepositUnit edu( energy, point);
00643 _ionization_points[i] = edu;
00644
00645 #ifdef TP_DEBUG
00646 LogDebug ("Pixel Digitizer")
00647 << i << " " << _ionization_points[i].x() << " "
00648 << _ionization_points[i].y() << " "
00649 << _ionization_points[i].z() << " "
00650 << _ionization_points[i].energy();
00651 #endif
00652
00653 }
00654
00655 delete[] elossVector;
00656
00657 }
00658
00659
00660
00661
00662 void SiPixelDigitizerAlgorithm::fluctuateEloss(int pid, float particleMomentum,
00663 float eloss, float length,
00664 int NumberOfSegs,float elossVector[]) {
00665
00666
00667 float dedx;
00668 if( length > 0.) dedx = eloss/length;
00669 else dedx = eloss;
00670
00671 double particleMass = 139.6;
00672 pid = abs(pid);
00673 if(pid!=211) {
00674 if(pid==11) particleMass = 0.511;
00675 else if(pid==13) particleMass = 105.7;
00676 else if(pid==321) particleMass = 493.7;
00677 else if(pid==2212) particleMass = 938.3;
00678 }
00679
00680 float segmentLength = length/NumberOfSegs;
00681
00682
00683 float de=0.;
00684 float sum=0.;
00685 double segmentEloss = (1000.*eloss)/NumberOfSegs;
00686 for (int i=0;i<NumberOfSegs;i++) {
00687
00688
00689
00690
00691
00692 double deltaCutoff = tMax;
00693 de = fluctuate->SampleFluctuations(double(particleMomentum*1000.),
00694 particleMass, deltaCutoff,
00695 double(segmentLength*10.),
00696 segmentEloss )/1000.;
00697
00698 elossVector[i]=de;
00699 sum +=de;
00700 }
00701
00702 if(sum>0.) {
00703
00704 float ratio = eloss/sum;
00705
00706 for (int ii=0;ii<NumberOfSegs;ii++) elossVector[ii]= ratio*elossVector[ii];
00707 } else {
00708 float averageEloss = eloss/NumberOfSegs;
00709 for (int ii=0;ii<NumberOfSegs;ii++) elossVector[ii]= averageEloss;
00710 }
00711 return;
00712 }
00713
00714
00715
00716
00717 void SiPixelDigitizerAlgorithm::drift(const PSimHit& hit){
00718
00719
00720 #ifdef TP_DEBUG
00721 LogDebug ("Pixel Digitizer") << " enter drift " ;
00722 #endif
00723
00724 _collection_points.resize( _ionization_points.size());
00725
00726 LocalVector driftDir=DriftDirection();
00727 if(driftDir.z() ==0.) {
00728 LogWarning("Magnetic field") << " pxlx: drift in z is zero ";
00729 return;
00730 }
00731
00732
00733
00734
00735
00736 float TanLorenzAngleX, TanLorenzAngleY,dir_z, CosLorenzAngleX,
00737 CosLorenzAngleY;
00738 if ( alpha2Order) {
00739
00740 TanLorenzAngleX = driftDir.x();
00741 TanLorenzAngleY = driftDir.y();
00742 dir_z = driftDir.z();
00743 CosLorenzAngleX = 1./sqrt(1.+TanLorenzAngleX*TanLorenzAngleX);
00744 CosLorenzAngleY = 1./sqrt(1.+TanLorenzAngleY*TanLorenzAngleY);
00745
00746 } else{
00747
00748 TanLorenzAngleX = driftDir.x();
00749 TanLorenzAngleY = 0.;
00750 dir_z = driftDir.z();
00751 CosLorenzAngleX = 1./sqrt(1.+TanLorenzAngleX*TanLorenzAngleX);
00752 CosLorenzAngleY = 1.;
00753 }
00754
00755
00756 #ifdef TP_DEBUG
00757 LogDebug ("Pixel Digitizer")
00758 << " Lorentz Tan " << TanLorenzAngleX << " " << TanLorenzAngleY <<" "
00759 << CosLorenzAngleX << " " << CosLorenzAngleY << " "
00760 << moduleThickness*TanLorenzAngleX << " " << driftDir;
00761 #endif
00762
00763 float Sigma_x = 1.;
00764 float Sigma_y = 1.;
00765 float DriftDistance;
00766 float DriftLength;
00767 float Sigma;
00768
00769
00770 for (unsigned int i = 0; i != _ionization_points.size(); i++) {
00771
00772 float SegX, SegY, SegZ;
00773 SegX = _ionization_points[i].x();
00774 SegY = _ionization_points[i].y();
00775 SegZ = _ionization_points[i].z();
00776
00777
00778
00779
00780 DriftDistance = moduleThickness/2. - (dir_z * SegZ);
00781
00782
00783
00784
00785
00786
00787
00788
00789
00790 if( DriftDistance < 0.) {
00791 DriftDistance = 0.;
00792 } else if ( DriftDistance > moduleThickness )
00793 DriftDistance = moduleThickness;
00794
00795
00796 float XDriftDueToMagField = DriftDistance * TanLorenzAngleX;
00797 float YDriftDueToMagField = DriftDistance * TanLorenzAngleY;
00798
00799
00800 float CloudCenterX = SegX + XDriftDueToMagField;
00801 float CloudCenterY = SegY + YDriftDueToMagField;
00802
00803
00804 DriftLength = sqrt( DriftDistance*DriftDistance +
00805 XDriftDueToMagField*XDriftDueToMagField +
00806 YDriftDueToMagField*YDriftDueToMagField );
00807
00808
00809 Sigma = sqrt(DriftLength/Dist300) * Sigma0;
00810
00811
00812 Sigma_x = Sigma / CosLorenzAngleX ;
00813 Sigma_y = Sigma / CosLorenzAngleY ;
00814
00815 SignalPoint sp( CloudCenterX, CloudCenterY,
00816 Sigma_x, Sigma_y, hit.tof(), _ionization_points[i].energy() );
00817
00818
00819 _collection_points[i] = (sp);
00820
00821 }
00822
00823 }
00824
00825
00826
00827 void SiPixelDigitizerAlgorithm::induce_signal( const PSimHit& hit) {
00828
00829
00830
00831
00832 #ifdef TP_DEBUG
00833 LogDebug ("Pixel Digitizer")
00834 << " enter induce_signal, "
00835 << topol->pitch().first << " " << topol->pitch().second;
00836 #endif
00837
00838
00839 typedef map< int, float, less<int> > hit_map_type;
00840 hit_map_type hit_signal;
00841
00842
00843 map<int, float, less<int> > x,y;
00844
00845
00846
00847
00848 for ( vector<SignalPoint>::const_iterator i=_collection_points.begin();
00849 i != _collection_points.end(); i++) {
00850
00851 float CloudCenterX = i->position().x();
00852 float CloudCenterY = i->position().y();
00853 float SigmaX = i->sigma_x();
00854 float SigmaY = i->sigma_y();
00855 float Charge = i->amplitude();
00856
00857
00858
00859
00860
00861
00862
00863
00864 #ifdef TP_DEBUG
00865 LogDebug ("Pixel Digitizer")
00866 << " cloud " << i->position().x() << " " << i->position().y() << " "
00867 << i->sigma_x() << " " << i->sigma_y() << " " << i->amplitude();
00868 #endif
00869
00870
00871 float CloudRight = CloudCenterX + ClusterWidth*SigmaX;
00872 float CloudLeft = CloudCenterX - ClusterWidth*SigmaX;
00873 float CloudUp = CloudCenterY + ClusterWidth*SigmaY;
00874 float CloudDown = CloudCenterY - ClusterWidth*SigmaY;
00875
00876
00877 LocalPoint PointRightUp = LocalPoint(CloudRight,CloudUp);
00878 LocalPoint PointLeftDown = LocalPoint(CloudLeft,CloudDown);
00879
00880
00881
00882
00883
00884
00885
00886
00887
00888 MeasurementPoint mp = topol->measurementPosition(PointRightUp );
00889
00890 int IPixRightUpX = int( floor( mp.x()));
00891 int IPixRightUpY = int( floor( mp.y()));
00892
00893 #ifdef TP_DEBUG
00894 LogDebug ("Pixel Digitizer") << " right-up " << PointRightUp << " "
00895 << mp.x() << " " << mp.y() << " "
00896 << IPixRightUpX << " " << IPixRightUpY ;
00897 #endif
00898
00899 mp = topol->measurementPosition(PointLeftDown );
00900
00901 int IPixLeftDownX = int( floor( mp.x()));
00902 int IPixLeftDownY = int( floor( mp.y()));
00903
00904 #ifdef TP_DEBUG
00905 LogDebug ("Pixel Digitizer") << " left-down " << PointLeftDown << " "
00906 << mp.x() << " " << mp.y() << " "
00907 << IPixLeftDownX << " " << IPixLeftDownY ;
00908 #endif
00909
00910
00911 IPixRightUpX = numRows>IPixRightUpX ? IPixRightUpX : numRows-1 ;
00912 IPixRightUpY = numColumns>IPixRightUpY ? IPixRightUpY : numColumns-1 ;
00913 IPixLeftDownX = 0<IPixLeftDownX ? IPixLeftDownX : 0 ;
00914 IPixLeftDownY = 0<IPixLeftDownY ? IPixLeftDownY : 0 ;
00915
00916 x.clear();
00917 y.clear();
00918
00919
00920 int ix;
00921 for (ix=IPixLeftDownX; ix<=IPixRightUpX; ix++) {
00922 float xUB, xLB, UpperBound, LowerBound;
00923
00924
00925
00926
00927 if (ix == 0 || SigmaX==0. )
00928 LowerBound = 0.;
00929 else {
00930 mp = MeasurementPoint( float(ix), 0.0);
00931 xLB = topol->localPosition(mp).x();
00932 gsl_sf_result result;
00933 int status = gsl_sf_erf_Q_e( (xLB-CloudCenterX)/SigmaX, &result);
00934 if (status != 0)
00935 LogWarning ("Integration")<<"could not compute gaussian probability";
00936 LowerBound = 1-result.val;
00937 }
00938
00939 if (ix == numRows-1 || SigmaX==0. )
00940 UpperBound = 1.;
00941 else {
00942 mp = MeasurementPoint( float(ix+1), 0.0);
00943 xUB = topol->localPosition(mp).x();
00944 gsl_sf_result result;
00945 int status = gsl_sf_erf_Q_e( (xUB-CloudCenterX)/SigmaX, &result);
00946 if (status != 0)
00947 LogWarning ("Integration")<<"could not compute gaussian probability";
00948 UpperBound = 1. - result.val;
00949 }
00950
00951 float TotalIntegrationRange = UpperBound - LowerBound;
00952 x[ix] = TotalIntegrationRange;
00953
00954
00955
00956 }
00957
00958
00959 int iy;
00960 for (iy=IPixLeftDownY; iy<=IPixRightUpY; iy++) {
00961 float yUB, yLB, UpperBound, LowerBound;
00962
00963 if (iy == 0 || SigmaY==0.)
00964 LowerBound = 0.;
00965 else {
00966 mp = MeasurementPoint( 0.0, float(iy) );
00967 yLB = topol->localPosition(mp).y();
00968 gsl_sf_result result;
00969 int status = gsl_sf_erf_Q_e( (yLB-CloudCenterY)/SigmaY, &result);
00970 if (status != 0)
00971 LogWarning ("Integration")<<"could not compute gaussian probability";
00972 LowerBound = 1. - result.val;
00973 }
00974
00975 if (iy == numColumns-1 || SigmaY==0. )
00976 UpperBound = 1.;
00977 else {
00978 mp = MeasurementPoint( 0.0, float(iy+1) );
00979 yUB = topol->localPosition(mp).y();
00980 gsl_sf_result result;
00981 int status = gsl_sf_erf_Q_e( (yUB-CloudCenterY)/SigmaY, &result);
00982 if (status != 0)
00983 LogWarning ("Integration")<<"could not compute gaussian probability";
00984 UpperBound = 1. - result.val;
00985 }
00986
00987 float TotalIntegrationRange = UpperBound - LowerBound;
00988 y[iy] = TotalIntegrationRange;
00989
00990
00991 }
00992
00993
00994 int chan;
00995 for (ix=IPixLeftDownX; ix<=IPixRightUpX; ix++) {
00996 for (iy=IPixLeftDownY; iy<=IPixRightUpY; iy++) {
00997
00998 float ChargeFraction = Charge*x[ix]*y[iy];
00999
01000 if( ChargeFraction > 0. ) {
01001 chan = PixelDigi::pixelToChannel( ix, iy);
01002
01003 hit_signal[chan] += ChargeFraction;
01004 }
01005
01006
01007 mp = MeasurementPoint( float(ix), float(iy) );
01008 LocalPoint lp = topol->localPosition(mp);
01009 chan = topol->channel(lp);
01010
01011 #ifdef TP_DEBUG
01012 LogDebug ("Pixel Digitizer")
01013 << " pixel " << ix << " " << iy << " - "<<" "
01014 << chan << " " << ChargeFraction<<" "
01015 << mp.x() << " " << mp.y() <<" "
01016 << lp.x() << " " << lp.y() << " "
01017 << chan;
01018 #endif
01019
01020 }
01021 }
01022
01023
01024
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039
01040
01041
01042 }
01043
01044
01045
01046 for ( hit_map_type::const_iterator im = hit_signal.begin();
01047 im != hit_signal.end(); im++) {
01048 _signal[(*im).first] += Amplitude( (*im).second, &hit, (*im).second);
01049
01050 int chan = (*im).first;
01051 pair<int,int> ip = PixelDigi::channelToPixel(chan);
01052
01053 #ifdef TP_DEBUG
01054 LogDebug ("Pixel Digitizer")
01055 << " pixel " << ip.first << " " << ip.second << " "
01056 << _signal[(*im).first];
01057 #endif
01058 }
01059
01060 }
01061
01062
01063
01064
01065 void SiPixelDigitizerAlgorithm::make_digis() {
01066 internal_coll.reserve(50); internal_coll.clear();
01067
01068 #ifdef TP_DEBUG
01069 LogDebug ("Pixel Digitizer") << " make digis "<<" "
01070 << " pixel threshold FPix" << theThresholdInE_FPix << " "
01071 << " pixel threshold BPix" << theThresholdInE_BPix << " "
01072 << " List pixels passing threshold ";
01073 #endif
01074
01075
01076
01077 for ( signal_map_iterator i = _signal.begin(); i != _signal.end(); i++) {
01078
01079 float signalInElectrons = (*i).second ;
01080
01081
01082
01083
01084
01085
01086 if ( signalInElectrons >= thePixelThresholdInE) {
01087
01088 int chan = (*i).first;
01089 pair<int,int> ip = PixelDigi::channelToPixel(chan);
01090 int adc=0;
01091
01092
01093 if(doMissCalibrate) {
01094 int row = ip.first;
01095 int col = ip.second;
01096 adc = int(missCalibrate(col,row,signalInElectrons));
01097 } else {
01098 adc = int( signalInElectrons / theElectronPerADC );
01099 }
01100 adc = min(adc, theAdcFullScale);
01101
01102 #ifdef TP_DEBUG
01103 LogDebug ("Pixel Digitizer")
01104 << (*i).first << " " << (*i).second << " " << signalInElectrons
01105 << " " << adc << ip.first << " " << ip.second ;
01106 #endif
01107
01108
01109 internal_coll.push_back( PixelDigi( ip.first, ip.second, adc));
01110
01111
01112 if((*i).second.hits().size()>0){
01113 simi.clear();
01114 unsigned int il=0;
01115 for( vector<const PSimHit*>::const_iterator ihit = (*i).second.hits().begin();
01116 ihit != (*i).second.hits().end(); ihit++) {
01117 simi[(**ihit).trackId()].push_back((*i).second.individualampl()[il]);
01118 il++;
01119 }
01120
01121
01122 for( simlink_map::iterator simiiter=simi.begin();
01123 simiiter!=simi.end();
01124 simiiter++){
01125
01126 float sum_samechannel=0;
01127 for (unsigned int iii=0;iii<(*simiiter).second.size();iii++){
01128 sum_samechannel+=(*simiiter).second[iii];
01129 }
01130 float fraction=sum_samechannel/(*i).second;
01131 if (fraction>1.) fraction=1.;
01132 link_coll.push_back(PixelDigiSimLink((*i).first,(*simiiter).first,((*i).second.hits().front())->eventId(),fraction));
01133 }
01134
01135 }
01136 }
01137
01138 }
01139 }
01140
01141
01142
01143
01144 void SiPixelDigitizerAlgorithm::add_noise() {
01145
01146 #ifdef TP_DEBUG
01147 LogDebug ("Pixel Digitizer") << " enter add_noise " << theNoiseInElectrons;
01148 #endif
01149
01150
01151
01152 for ( signal_map_iterator i = _signal.begin(); i != _signal.end(); i++) {
01153
01154 if(addChargeVCALSmearing)
01155 {
01156 if((*i).second < 3000)
01157 {
01158 theSmearedChargeRMS = 543.6 - (*i).second * 0.093;
01159 } else if((*i).second < 6000){
01160 theSmearedChargeRMS = 307.6 - (*i).second * 0.01;
01161 } else{
01162 theSmearedChargeRMS = -432.4 +(*i).second * 0.123;
01163 }
01164
01165
01166 float noise_ChargeVCALSmearing = theSmearedChargeRMS * gaussDistributionVCALNoise_->fire() ;
01167
01168 float noise = gaussDistribution_->fire() ;
01169
01170 if(((*i).second + Amplitude(noise+noise_ChargeVCALSmearing,0,-1.)) < 0. ) {
01171 (*i).second.set(0);}
01172 else{
01173 (*i).second +=Amplitude(noise+noise_ChargeVCALSmearing,0,-1.);
01174 }
01175
01176 }
01177 else
01178 {
01179
01180
01181
01182 float noise = gaussDistribution_->fire() ;
01183 if(((*i).second + Amplitude(noise,0,-1.)) < 0. ) {
01184 (*i).second.set(0);}
01185 else{
01186 (*i).second +=Amplitude(noise,0,-1.);
01187 }
01188 }
01189
01190 }
01191
01192 if(!addNoisyPixels)
01193 return;
01194
01195
01196
01197 int numberOfPixels = (numRows * numColumns);
01198 map<int,float, less<int> > otherPixels;
01199 map<int,float, less<int> >::iterator mapI;
01200
01201 theNoiser->generate(numberOfPixels,
01202 thePixelThreshold,
01203 theNoiseInElectrons,
01204 otherPixels );
01205
01206 #ifdef TP_DEBUG
01207 LogDebug ("Pixel Digitizer")
01208 << " Add noisy pixels " << numRows << " "
01209 << numColumns << " " << theNoiseInElectrons << " "
01210 << theThresholdInE_FPix << theThresholdInE_BPix <<" "<< numberOfPixels<<" "
01211 << otherPixels.size() ;
01212 #endif
01213
01214
01215 for (mapI = otherPixels.begin(); mapI!= otherPixels.end(); mapI++) {
01216 int iy = ((*mapI).first) / numRows;
01217 int ix = ((*mapI).first) - (iy*numRows);
01218
01219
01220 if( iy < 0 || iy > (numColumns-1) )
01221 LogWarning ("Pixel Geometry") << " error in iy " << iy ;
01222 if( ix < 0 || ix > (numRows-1) )
01223 LogWarning ("Pixel Geometry") << " error in ix " << ix ;
01224
01225 int chan = PixelDigi::pixelToChannel(ix, iy);
01226
01227 #ifdef TP_DEBUG
01228 LogDebug ("Pixel Digitizer")
01229 <<" Storing noise = " << (*mapI).first << " " << (*mapI).second
01230 << " " << ix << " " << iy << " " << chan ;
01231 #endif
01232
01233 if(_signal[chan] == 0){
01234
01235 int noise=int( (*mapI).second );
01236 _signal[chan] = Amplitude (noise, 0,-1.);
01237 }
01238 }
01239
01240 }
01241
01242
01243
01244
01245
01246 void SiPixelDigitizerAlgorithm::pixel_inefficiency() {
01247
01248
01249
01250 float pixelEfficiency = 1.0;
01251 float columnEfficiency = 1.0;
01252 float chipEfficiency = 1.0;
01253
01254
01255 unsigned int Subid=DetId(detID).subdetId();
01256 if (Subid== PixelSubdetector::PixelBarrel){
01257 int layerIndex=PXBDetId(detID).layer();
01258 pixelEfficiency = thePixelEfficiency[layerIndex-1];
01259 columnEfficiency = thePixelColEfficiency[layerIndex-1];
01260 chipEfficiency = thePixelChipEfficiency[layerIndex-1];
01261
01262
01263 if(numColumns>416) LogWarning ("Pixel Geometry") <<" wrong columns in barrel "<<numColumns;
01264 if(numRows>160) LogWarning ("Pixel Geometry") <<" wrong rows in barrel "<<numRows;
01265
01266 } else {
01267
01268
01269 pixelEfficiency = thePixelEfficiency[3];
01270 columnEfficiency = thePixelColEfficiency[3];
01271 chipEfficiency = thePixelChipEfficiency[3];
01272
01273
01274
01275 if(numColumns>260 || numRows>160) {
01276 if(numColumns>260) LogWarning ("Pixel Geometry") <<" wrong columns in endcaps "<<numColumns;
01277 if(numRows>160) LogWarning ("Pixel Geometry") <<" wrong rows in endcaps "<<numRows;
01278 return;
01279 }
01280 }
01281
01282 #ifdef TP_DEBUG
01283 LogDebug ("Pixel Digitizer") << " enter pixel_inefficiency " << pixelEfficiency << " "
01284 << columnEfficiency << " " << chipEfficiency;
01285 #endif
01286
01287
01288
01289 int chipIndex = 0;
01290 int rowROC = 0;
01291 int colROC = 0;
01292 map<int, int, less<int> >chips, columns;
01293 map<int, int, less<int> >::iterator iter;
01294
01295
01296
01297 for (signal_map_iterator i = _signal.begin();i != _signal.end();i++) {
01298
01299 int chan = i->first;
01300 pair<int,int> ip = PixelDigi::channelToPixel(chan);
01301 int row = ip.first;
01302 int col = ip.second;
01303
01304 pIndexConverter->transformToROC(col,row,chipIndex,colROC,rowROC);
01305 int dColInChip = pIndexConverter->DColumn(colROC);
01306
01307 int dColInDet = pIndexConverter->DColumnInModule(dColInChip,chipIndex);
01308
01309 chips[chipIndex]++;
01310 columns[dColInDet]++;
01311 }
01312
01313
01314 for ( iter = chips.begin(); iter != chips.end() ; iter++ ) {
01315
01316 float rand = flatDistribution_->fire();
01317 if( rand > chipEfficiency ) chips[iter->first]=0;
01318 }
01319
01320
01321 for ( iter = columns.begin(); iter != columns.end() ; iter++ ) {
01322
01323 float rand = flatDistribution_->fire();
01324 if( rand > columnEfficiency ) columns[iter->first]=0;
01325 }
01326
01327
01328
01329 for(signal_map_iterator i = _signal.begin();i != _signal.end(); i++) {
01330
01331
01332 pair<int,int> ip = PixelDigi::channelToPixel(i->first);
01333 int row = ip.first;
01334 int col = ip.second;
01335
01336 pIndexConverter->transformToROC(col,row,chipIndex,colROC,rowROC);
01337 int dColInChip = pIndexConverter->DColumn(colROC);
01338
01339 int dColInDet = pIndexConverter->DColumnInModule(dColInChip,chipIndex);
01340
01341
01342
01343 float rand = flatDistribution_->fire();
01344 if( chips[chipIndex]==0 || columns[dColInDet]==0
01345 || rand>pixelEfficiency ) {
01346
01347 i->second.set(0.);
01348 }
01349
01350 }
01351
01352 }
01353
01354
01355
01356
01357
01358
01359
01360
01361
01362
01363 float SiPixelDigitizerAlgorithm::missCalibrate(int col,int row,
01364 const float signalInElectrons) const {
01365
01366
01367
01368
01369
01370
01371
01372
01373
01374
01375
01376 float p0=0.0;
01377 float p1=0.0;
01378 float p2=0.0;
01379 float p3=0.0;
01380
01381 unsigned int Sub_detid=DetId(detID).subdetId();
01382
01383 if (Sub_detid == PixelSubdetector::PixelBarrel){
01384 p0 = BPix_p0;
01385 p1 = BPix_p1;
01386 p2 = BPix_p2;
01387 p3 = BPix_p3;
01388 } else {
01389 p0 = FPix_p0;
01390 p1 = FPix_p1;
01391 p2 = FPix_p2;
01392 p3 = FPix_p3;
01393 }
01394
01395
01396
01397 float newAmp = 0.;
01398
01399
01400 float signal = (signalInElectrons-electronsPerVCAL_Offset)/electronsPerVCAL;
01401
01402
01403 newAmp = p3 + p2 * tanh(p0*signal - p1);
01404
01405
01406
01407
01408
01409
01410
01411
01412
01413
01414
01415
01416
01417
01418
01419
01420
01421
01422
01423
01424
01425
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435 return newAmp;
01436 }
01437
01438
01439
01440
01441
01442
01443
01444
01445 LocalVector SiPixelDigitizerAlgorithm::DriftDirection(){
01446 Frame detFrame(_detp->surface().position(),_detp->surface().rotation());
01447 LocalVector Bfield=detFrame.toLocal(_bfield);
01448
01449 float alpha2_FPix;
01450 float alpha2_BPix;
01451 float alpha2;
01452
01453
01454
01455
01456
01457
01458
01459 float dir_x = 0.0;
01460 float dir_y = 0.0;
01461 float dir_z = 0.0;
01462 float scale = 0.0;
01463
01464 unsigned int Sub_detid=DetId(detID).subdetId();
01465
01466
01467
01468 if(!use_LorentzAngle_DB_){
01469
01470 if ( alpha2Order) {
01471 alpha2_FPix = tanLorentzAnglePerTesla_FPix*tanLorentzAnglePerTesla_FPix;
01472 alpha2_BPix = tanLorentzAnglePerTesla_BPix*tanLorentzAnglePerTesla_BPix;
01473 }else {
01474 alpha2_FPix = 0.0;
01475 alpha2_BPix = 0.0;
01476 }
01477
01478 if (Sub_detid == PixelSubdetector::PixelBarrel){
01479 dir_x = -( tanLorentzAnglePerTesla_BPix * Bfield.y() + alpha2_BPix* Bfield.z()* Bfield.x() );
01480 dir_y = +( tanLorentzAnglePerTesla_BPix * Bfield.x() - alpha2_BPix* Bfield.z()* Bfield.y() );
01481 dir_z = -(1 + alpha2_BPix* Bfield.z()*Bfield.z() );
01482 scale = (1 + alpha2_BPix* Bfield.z()*Bfield.z() );
01483
01484 } else {
01485 dir_x = -( tanLorentzAnglePerTesla_FPix * Bfield.y() + alpha2_FPix* Bfield.z()* Bfield.x() );
01486 dir_y = +( tanLorentzAnglePerTesla_FPix * Bfield.x() - alpha2_FPix* Bfield.z()* Bfield.y() );
01487 dir_z = -(1 + alpha2_FPix* Bfield.z()*Bfield.z() );
01488 scale = (1 + alpha2_FPix* Bfield.z()*Bfield.z() );
01489 }
01490 }
01491
01492
01493 if(use_LorentzAngle_DB_){
01494 std::map<unsigned int,float> detid_la= SiPixelLorentzAngle_->getLorentzAngles();
01495 std::map<unsigned int,float>::const_iterator it;
01496
01497
01498 for (it=detid_la.begin();it!=detid_la.end();it++)
01499 {
01500 if (detID==it->first) {
01501 if (alpha2Order) {
01502 alpha2 = it->second * it->second;
01503 }
01504 else {
01505 alpha2 = 0.0;
01506 }
01507
01508 dir_x = -( it->second * Bfield.y() + alpha2 * Bfield.z()* Bfield.x() );
01509 dir_y = +( it->second * Bfield.x() - alpha2 * Bfield.z()* Bfield.y() );
01510 dir_z = -(1 + alpha2 * Bfield.z()*Bfield.z() );
01511 scale = (1 + alpha2 * Bfield.z()*Bfield.z() );
01512 }
01513 }
01514 }
01515
01516 LocalVector theDriftDirection = LocalVector(dir_x/scale, dir_y/scale, dir_z/scale );
01517
01518 #ifdef TP_DEBUG
01519 LogDebug ("Pixel Digitizer") << " The drift direction in local coordinate is "
01520 << theDriftDirection ;
01521 #endif
01522
01523 return theDriftDirection;
01524 }
01525
01526
01527
01528 void SiPixelDigitizerAlgorithm::pixel_inefficiency_db(void){
01529 if(!use_ineff_from_db_)
01530 return;
01531
01532
01533 for(signal_map_iterator i = _signal.begin();i != _signal.end(); i++) {
01534
01535
01536 pair<int,int> ip = PixelDigi::channelToPixel(i->first);
01537 int row = ip.first;
01538 int col = ip.second;
01539 uint32_t detid = detID;
01540
01541 if(theSiPixelGainCalibrationService_->isDead(detid, col, row)){
01542
01543
01544 i->second.set(0.);
01545 }
01546 }
01547 }
01548
01549
01550
01551
01552 void SiPixelDigitizerAlgorithm::module_killing_conf(void){
01553 if(!use_module_killing_)
01554 return;
01555
01556 bool isbad=false;
01557 int detid = detID;
01558
01559 Parameters::iterator itDeadModules=DeadModules.begin();
01560
01561 for(; itDeadModules != DeadModules.end(); ++itDeadModules){
01562 int Dead_detID = itDeadModules->getParameter<int>("Dead_detID");
01563 if(detid==Dead_detID){
01564 isbad=true;
01565 break;
01566 }
01567 }
01568
01569 if(!isbad)
01570 return;
01571
01572 std::string Module = itDeadModules->getParameter<std::string>("Module");
01573
01574 if(Module=="whole"){
01575 for(signal_map_iterator i = _signal.begin();i != _signal.end(); i++) {
01576 i->second.set(0.);
01577 }
01578 }
01579
01580 for(signal_map_iterator i = _signal.begin();i != _signal.end(); i++) {
01581 pair<int,int> ip = PixelDigi::channelToPixel(i->first);
01582
01583 if(Module=="tbmA" && ip.first>=80 && ip.first<=159){
01584 i->second.set(0.);
01585 }
01586
01587 if( Module=="tbmB" && ip.first<=79){
01588 i->second.set(0.);
01589
01590 }
01591
01592 }
01593 }
01594
01595
01596
01597 void SiPixelDigitizerAlgorithm::module_killing_DB(void){
01598 if(!use_module_killing_)
01599 return;
01600
01601 bool isbad=false;
01602 uint32_t detid = detID;
01603
01604 std::vector<SiPixelQuality::disabledModuleType>disabledModules = SiPixelBadModule_->getBadComponentList();
01605
01606 SiPixelQuality::disabledModuleType badmodule;
01607
01608 for (size_t id=0;id<disabledModules.size();id++)
01609 {
01610 if(detid==disabledModules[id].DetID){
01611 isbad=true;
01612 badmodule = disabledModules[id];
01613 break;
01614 }
01615 }
01616
01617 if(!isbad)
01618 return;
01619
01620 if(badmodule.errorType == 0){
01621
01622 for(signal_map_iterator i = _signal.begin();i != _signal.end(); i++) {
01623 i->second.set(0.);
01624 }
01625 }
01626 else {
01627
01628
01629 std::vector<GlobalPixel> badrocpositions (0);
01630 std::pair<uint8_t, uint8_t> coord(1,1);
01631 for(unsigned int j = 0; j < 16; j++){
01632 if (SiPixelBadModule_->IsRocBad(detid, j) == true){
01633
01634 std::vector<CablingPathToDetUnit> path = map_.product()->pathToDetUnit(detid);
01635 typedef std::vector<CablingPathToDetUnit>::const_iterator IT;
01636 for (IT it = path.begin(); it != path.end(); ++it) {
01637 const PixelROC* myroc = map_.product()->findItem(*it);
01638 if( myroc->idInDetUnit() == j) {
01639 LocalPixel::RocRowCol local = { 39, 25};
01640 GlobalPixel global = myroc->toGlobal( LocalPixel(local) );
01641 badrocpositions.push_back(global);
01642 break;
01643 }
01644 }
01645 }
01646 }
01647
01648 for(signal_map_iterator i = _signal.begin();i != _signal.end(); i++) {
01649 pair<int,int> ip = PixelDigi::channelToPixel(i->first);
01650
01651 for(std::vector<GlobalPixel>::const_iterator it = badrocpositions.begin(); it != badrocpositions.end(); ++it){
01652
01653 if(it->row==120){
01654 if((ip.first>=(it->row)-40 && ip.first<=(it->row)+39) && (ip.second>=(it->col)-25 && ip.second<=(it->col)+26)){
01655 i->second.set(0.);
01656 }
01657 }
01658 else {
01659 if((ip.first>=(it->row)-39 && ip.first<=(it->row)+40) && (ip.second>=(it->col)-25 && ip.second<=(it->col)+26)){
01660 i->second.set(0.);
01661 }
01662 }
01663 }
01664 }
01665 }
01666 }
01667