CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/SimTracker/SiPixelDigitizer/src/SiPixelDigitizerAlgorithm.cc

Go to the documentation of this file.
00001 //class SiPixelDigitizerAlgorithm SimTracker/SiPixelDigitizer/src/SiPixelDigitizerAlgoithm.cc 
00002 
00003 // Original Author Danek Kotlinski
00004 // Ported in CMSSW by  Michele Pioppi-INFN perugia
00005 // Added DB capabilities by F.Blekman, Cornell University
00006 //         Created:  Mon Sep 26 11:08:32 CEST 2005
00007 // Add tof, change AddNoise to tracked. 4/06
00008 // Change drift direction. 6/06 d.k.
00009 // Add the statuis (non-rate dependent) inefficiency.
00010 //     -1 - no ineffciency
00011 //      0 - static inefficency only
00012 //    1,2 - low-lumi rate dependent inefficency added
00013 //     10 - high-lumi inefficiency added
00014 // Adopt the correct drift sign convetion from Morris Swartz. d.k. 8/06
00015 // Add more complex misscalinbration, change kev/e to 3.61, diff=3.7,d.k.9/06
00016 // Add the readout channel electronic noise. d.k. 3/07
00017 // Lower the pixel noise from 500 to 175elec.
00018 // Change the input threshold from noise units to electrons.
00019 // Lower the amount of static dead pixels from 0.01 to 0.001.
00020 // Modify to the new random number services. d.k. 5/07
00021 // Protect against sigma=0 (delta tracks on the surface). d.k.5/07
00022 // Change the TOF cut to lower and upper limit. d.k. 7/07
00023 //
00024 // July 2008: Split Lorentz Angle configuration in BPix/FPix (V. Cuplov)
00025 // tanLorentzAngleperTesla_FPix=0.0912 and tanLorentzAngleperTesla_BPix=0.106
00026 // Sept. 2008: Disable Pixel modules which are declared dead in the configuration python file. (V. Cuplov)
00027 // Oct. 2008: Accessing/Reading the Lorentz angle from the DataBase instead of the cfg file. (V. Cuplov)
00028 // Accessing dead modules from the DB. Implementation done and tested on a test.db
00029 // Do not use this option for now. The PixelQuality Objects are not in the official DB yet.
00030 // Feb. 2009: Split Fpix and Bpix threshold and use official numbers (V. Cuplov)
00031 // ThresholdInElectrons_FPix = 2870 and ThresholdInElectrons_BPix = 3700
00032 // update the electron to VCAL conversion using: VCAL_electrons = VCAL * 65.5 - 414
00033 // Feb. 2009: Threshold gaussian smearing (V. Cuplov)
00034 // March, 2009: changed DB access to *SimRcd objects (to de-couple the DB objects from reco chain) (F. Blekman) 
00035 // May, 2009: Pixel charge VCAL smearing. (V. Cuplov)
00036 // November, 2009: new parameterization of the pixel response. (V. Cuplov)
00037 // December, 2009: Fix issue with different compilers.
00038 // October, 2010: Improvement: Removing single dead ROC (V. Cuplov)
00039 // November, 2010: Bug fix in removing TBMB/A half-modules (V. Cuplov)
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 //#include "SimTracker/SiPixelDigitizer/interface/PixelIndices.h"
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 // Accessing dead pixel modules from the DB:
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 // Geometry
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_){// load gain calibration service fromdb...
00098     theSiPixelGainCalibrationService_= new SiPixelGainCalibrationOfflineSimService(conf_);
00099     theSiPixelGainCalibrationService_->setESObjects( es );
00100   }
00101 
00102   fillDeadModules(es); // gets the dead module from config file or DB.
00103   fillLorentzAngle(es); // gets the Lorentz angle from the config file or DB.
00104   fillMapandGeom(es);  //gets the map and geometry from the DB (to kill ROCs)  
00105 }
00106 
00107 //=========================================================================
00108 
00109 void SiPixelDigitizerAlgorithm::fillDeadModules(const edm::EventSetup& es){
00110   if(!use_deadmodule_DB_){
00111     DeadModules = conf_.getParameter<Parameters>("DeadModules"); // get dead module from cfg file
00112   }
00113   else{  // Get dead module from DB record 
00114     // ESHandle was defined in the header file   edm::ESHandle<SiPixelQuality> SiPixelBadModule_;
00115     es.get<SiPixelQualityRcd>().get(SiPixelBadModule_);
00116   }
00117 }
00118 
00119 void SiPixelDigitizerAlgorithm::fillLorentzAngle(const edm::EventSetup& es){
00120   if(!use_LorentzAngle_DB_){
00121     // Get the Lorentz angle from the cfg file:
00122     tanLorentzAnglePerTesla_FPix=conf_.getParameter<double>("TanLorentzAnglePerTesla_FPix");
00123     tanLorentzAnglePerTesla_BPix=conf_.getParameter<double>("TanLorentzAnglePerTesla_BPix");
00124   } 
00125   else {
00126     // Get Lorentz angle from DB record 
00127     // ESHandle was defined in the header file edm::ESHandle<SiPixelLorentzAngle> SiPixelLorentzAngle_;
00128     es.get<SiPixelLorentzAngleSimRcd>().get(SiPixelLorentzAngle_);
00129   }
00130 }
00131 
00132 // For ROC killing:
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")), // boolean to kill or not modules
00144   use_deadmodule_DB_(conf_.getParameter<bool>("DeadModules_DB")), // boolean to access dead modules from DB
00145   use_LorentzAngle_DB_(conf_.getParameter<bool>("LorentzAngle_DB")), // boolean to access Lorentz angle from DB 
00146   theSiPixelGainCalibrationService_(0),rndEngine(eng)
00147 {
00148   using std::cout;
00149   using std::endl;
00150 
00151   // Common pixel parameters
00152   // These are parameters which are not likely to be changed
00153   NumberOfSegments = 20; // Default number of track segment divisions
00154   ClusterWidth = 3.;     // Charge integration spread on the collection plane
00155   GeVperElectron = 3.61E-09; // 1 electron =3.61eV, 1keV=277e, mod 9/06 d.k.
00156   Sigma0 = 0.00037;           // Charge diffusion constant 7->3.7 
00157   Dist300 = 0.0300;          //   normalized to 300micron Silicon
00158 
00159   alpha2Order = conf_.getParameter<bool>("Alpha2Order");   // switch on/off of E.B effect   
00160 
00161   // get external parameters:
00162   // ADC calibration 1adc count = 135e.
00163   // Corresponds to 2adc/kev, 270[e/kev]/135[e/adc]=2[adc/kev]
00164   // Be carefull, this parameter is also used in SiPixelDet.cc to 
00165   // calculate the noise in adc counts from noise in electrons.
00166   // Both defaults should be the same.
00167   theElectronPerADC=conf_.getParameter<double>("ElectronPerAdc");
00168 
00169   // ADC saturation value, 255=8bit adc.
00170   //theAdcFullScale=conf_.getUntrackedParameter<int>("AdcFullScale",255);
00171   theAdcFullScale=conf_.getParameter<int>("AdcFullScale");
00172 
00173   // Pixel threshold in units of noise:
00174   // thePixelThreshold=conf_.getParameter<double>("ThresholdInNoiseUnits");
00175   // Pixel threshold in electron units.
00176   theThresholdInE_FPix=conf_.getParameter<double>("ThresholdInElectrons_FPix");
00177   theThresholdInE_BPix=conf_.getParameter<double>("ThresholdInElectrons_BPix");
00178 
00179   // Add threshold gaussian smearing:
00180   addThresholdSmearing = conf_.getParameter<bool>("AddThresholdSmearing");
00181   theThresholdSmearing_FPix = conf_.getParameter<double>("ThresholdSmearing_FPix");
00182   theThresholdSmearing_BPix = conf_.getParameter<double>("ThresholdSmearing_BPix");
00183 
00184   // signal response new parameterization: split Fpix and BPix 
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   // electrons to VCAL conversion needed in misscalibrate()
00196   electronsPerVCAL = conf_.getParameter<double>("ElectronsPerVcal");
00197   electronsPerVCAL_Offset = conf_.getParameter<double>("ElectronsPerVcal_Offset");
00198 
00199   // Add noise   
00200   addNoise=conf_.getParameter<bool>("AddNoise");
00201 
00202   // Smear the pixel charge with a gaussian which RMS is a function of the
00203   // pixel charge (Danek's study)
00204   addChargeVCALSmearing=conf_.getParameter<bool>("ChargeVCALSmearing");
00205 
00206   // Add noisy pixels 
00207   addNoisyPixels=conf_.getParameter<bool>("AddNoisyPixels");
00208   // Noise in electrons:
00209   // Pixel cell noise, relevant for generating noisy pixels 
00210   theNoiseInElectrons=conf_.getParameter<double>("NoiseInElectrons");
00211   // Fill readout noise, including all readout chain, relevant for smearing
00212   //theReadoutNoise=conf_.getUntrackedParameter<double>("ReadoutNoiseInElec",500.);
00213   theReadoutNoise=conf_.getParameter<double>("ReadoutNoiseInElec");
00214 
00215   //theTofCut 12.5, cut in particle TOD +/- 12.5ns
00216   //theTofCut=conf_.getUntrackedParameter<double>("TofCut",12.5);
00217   theTofLowerCut=conf_.getParameter<double>("TofLowerCut");
00218   theTofUpperCut=conf_.getParameter<double>("TofUpperCut");
00219   
00220   // Fluctuate charge in track subsegments
00221   fluctuateCharge=conf_.getUntrackedParameter<bool>("FluctuateCharge",true);
00222 
00223   // delta cutoff in MeV, has to be same as in OSCAR=0.030/cmsim=1.0 MeV
00224   //tMax = 0.030; // In MeV.  
00225   //tMax =conf_.getUntrackedParameter<double>("DeltaProductionCut",0.030);  
00226   tMax =conf_.getParameter<double>("DeltaProductionCut");  
00227  
00228   // Control the pixel inefficiency
00229   thePixelLuminosity=conf_.getParameter<int>("AddPixelInefficiency");
00230 
00231   // Get the constants for the miss-calibration studies
00232   doMissCalibrate=conf_.getParameter<bool>("MissCalibrate"); // Enable miss-calibration
00233   theGainSmearing=conf_.getParameter<double>("GainSmearing"); // sigma of the gain smearing
00234   theOffsetSmearing=conf_.getParameter<double>("OffsetSmearing"); //sigma of the offset smearing
00235 
00236   //pixel inefficiency
00237   // the first 3 settings [0],[1],[2] are for the barrel pixels
00238   // the next  3 settings [3],[4],[5] are for the endcaps (undecided how)  
00239 
00240   if (thePixelLuminosity==-1) {  // No inefficiency, all 100% efficient
00241     pixelInefficiency=false;
00242     for (int i=0; i<6;i++) {
00243       thePixelEfficiency[i]     = 1.;  // pixels = 100%
00244       thePixelColEfficiency[i]  = 1.;  // columns = 100%
00245       thePixelChipEfficiency[i] = 1.; // chips = 100%
00246     }
00247     
00248     // include only the static (non rate depedent) efficiency 
00249     // Usefull for very low rates (luminosity)
00250   } else if (thePixelLuminosity==0) { // static effciency
00251     pixelInefficiency=true;
00252     // Default efficiencies 
00253     for (int i=0; i<6;i++) {
00254       if(i<3) {  // For the barrel
00255         // Assume 1% inefficiency for single pixels, 
00256         // this is given by faulty bump-bonding and seus.  
00257         thePixelEfficiency[i]     = 1.-0.001;  // pixels = 99.9%
00258         // For columns make 0.1% default.
00259         thePixelColEfficiency[i]  = 1.-0.001;  // columns = 99.9%
00260         // A flat 0.1% inefficiency due to lost rocs
00261         thePixelChipEfficiency[i] = 1.-0.001; // chips = 99.9%
00262       } else { // For the endcaps
00263         // Assume 1% inefficiency for single pixels, 
00264         // this is given by faulty bump-bonding and seus.  
00265         thePixelEfficiency[i]     = 1.-0.001;  // pixels = 99.9%
00266         // For columns make 0.1% default.
00267         thePixelColEfficiency[i]  = 1.-0.001;  // columns = 99.9%
00268         // A flat 0.1% inefficiency due to lost rocs
00269         thePixelChipEfficiency[i] = 1.-0.001; // chips = 99.9%
00270       }
00271     }
00272     
00273     // Include also luminosity rate dependent inefficieny
00274   } else if (thePixelLuminosity>0) { // Include effciency
00275     pixelInefficiency=true;
00276     // Default efficiencies 
00277     for (int i=0; i<6;i++) {
00278       if(i<3) { // For the barrel
00279         // Assume 1% inefficiency for single pixels, 
00280         // this is given by faulty bump-bonding and seus.  
00281         thePixelEfficiency[i]     = 1.-0.01;  // pixels = 99%
00282         // For columns make 1% default.
00283         thePixelColEfficiency[i]  = 1.-0.01;  // columns = 99%
00284         // A flat 0.25% inefficiency due to lost data packets from TBM
00285         thePixelChipEfficiency[i] = 1.-0.0025; // chips = 99.75%
00286       } else { // For the endcaps
00287         // Assume 1% inefficiency for single pixels, 
00288         // this is given by faulty bump-bonding and seus.  
00289         thePixelEfficiency[i]     = 1.-0.01;  // pixels = 99%
00290         // For columns make 1% default.
00291         thePixelColEfficiency[i]  = 1.-0.01;  // columns = 99%
00292         // A flat 0.25% inefficiency due to lost data packets from TBM
00293         thePixelChipEfficiency[i] = 1.-0.0025; // chips = 99.75%
00294       }
00295     }
00296    
00297     // Special cases ( High-lumi for 4cm layer) where the readout losses are higher
00298     if(thePixelLuminosity==10) { // For high luminosity, bar layer 1
00299       thePixelColEfficiency[0] = 1.-0.034; // 3.4% for r=4 only
00300       thePixelEfficiency[0]    = 1.-0.015; // 1.5% for r=4
00301     }
00302     
00303   } // end the pixel inefficiency part
00304 
00305   // Init the random number services  
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     // Threshold smearing with gaussian distribution:
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     } //end Init the random number services
00325 
00326 
00327   // Prepare for the analog amplitude miss-calibration
00328   if(doMissCalibrate) {
00329     LogDebug ("PixelDigitizer ") 
00330       << " miss-calibrate the pixel amplitude "; 
00331     
00332     const bool ReadCalParameters = false;
00333     if(ReadCalParameters) {   // Read the calibration files from file
00334       // read the calibration constants from a file (testing only)
00335       ifstream in_file;  // data file pointer
00336       char filename[80] = "phCalibrationFit_C0.dat";
00337       
00338       in_file.open(filename, ios::in ); // in C++
00339       if (in_file.bad()) {
00340         cout << " File not found " << endl;
00341         return; // signal error
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       // Read MC tracks
00357       for(int i=0;i<(52*80);i++)  { // loop over tracks    
00358         in_file >> par0 >> par1 >> par2 >> par3 >> name >> colid
00359                 >> rowid;
00360         if (in_file.bad()) { // check for errors
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         //cout << " line " << i << " " <<par0<<" "<<par1<<" "<<par2<<" "<<par3<<" "
00372         //   <<colid<<" "<<rowid<<endl;
00373         
00374         CalParameters onePix;
00375         onePix.p0=par0;
00376         onePix.p1=par1;
00377         onePix.p2=par2;
00378         onePix.p3=par3;
00379         
00380         // Convert ROC pixel index to channel 
00381         int chan = PixelIndices::pixelToChannelROC(rowid,colid);
00382         calmap.insert(pair<int,CalParameters>(chan,onePix));
00383         
00384         // Testing the index conversion, can be skipped 
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       } // pixel loop in a ROC
00390  
00391       cout << " map size  " << calmap.size() <<" max "<<calmap.max_size() << " "
00392            <<calmap.empty()<< endl;
00393       
00394 //     cout << " map size  " << calmap.size()  << endl;
00395 //     map<int,CalParameters,less<int> >::iterator ix,it;
00396 //     map<int,CalParameters,less<int> >::const_iterator ip;
00397 //     for (ix = calmap.begin(); ix != calmap.end(); ++ix) {
00398 //       int i = (*ix).first;
00399 //       pair<int,int> p = channelToPixelROC(i);
00400 //       it  = calmap.find(i);
00401 //       CalParameters y  = (*it).second;
00402 //       CalParameters z = (*ix).second;
00403 //       cout << i <<" "<<p.first<<" "<<p.second<<" "<<y.p0<<" "<<z.p0<<" "<<calmap[i].p0<<endl; 
00404       
00405 //       //int dummy=0;
00406 //       //cin>>dummy;
00407 //     }
00408 
00409     } // end if readparameters
00410   } // end if missCalibration 
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   // Destructor
00430   delete gaussDistribution_;
00431   delete gaussDistributionVCALNoise_;
00432   delete flatDistribution_;
00433 
00434   delete theSiPixelGainCalibrationService_;
00435 
00436   // Threshold gaussian smearing:
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; //cache the PixGeomDetUnit
00455   _PixelHits=input; //cache the SimHit
00456   _bfield=bfield; //cache the drift direction
00457 
00458   // Pixel Efficiency moved from the constructor to the method run because
00459   // the information of the det are not available in the constructor
00460   // Effciency parameters. 0 - no inefficiency, 1-low lumi, 10-high lumi
00461   
00462   detID= _detp->geographicalId().rawId();
00463   
00464   _signal.clear();
00465   
00466   // initalization  of pixeldigisimlinks
00467   link_coll.clear();
00468   
00469   //Digitization of the SimHits of a given pixdet
00470   vector<PixelDigi> collector =digitize(pixdet);
00471 
00472   // edm::DetSet<PixelDigi> collector;
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(); // cache topology
00487     numColumns = topol->ncolumns();  // det module number of cols&rows
00488     numRows = topol->nrows();
00489     
00490     // full detector thickness
00491     moduleThickness = det->specificSurface().bounds().thickness(); 
00492     
00493     // The index converter is only needed when inefficiencies or misscalibration
00494     // are simulated.
00495     if((pixelInefficiency>0) || doMissCalibrate ) {  // Init pixel indices
00496       pIndexConverter = new PixelIndices(numColumns,numRows);
00497     }
00498     
00499     // Noise already defined in electrons
00500     //thePixelThresholdInE = thePixelThreshold * theNoiseInElectrons ; 
00501     // Find the threshold in noise units, needed for the noiser.
00502 
00503   unsigned int Sub_detid=DetId(detID).subdetId();
00504 
00505   if(theNoiseInElectrons>0.){
00506     if(Sub_detid == PixelSubdetector::PixelBarrel){ // Barrel modules
00507       if(addThresholdSmearing) { 
00508         thePixelThresholdInE = smearedThreshold_BPix_->fire(); // gaussian smearing 
00509       } else {
00510         thePixelThresholdInE = theThresholdInE_BPix; // no smearing
00511       }
00512       
00513       thePixelThreshold = thePixelThresholdInE/theNoiseInElectrons; 
00514       
00515     } else { // Forward disks modules 
00516       if(addThresholdSmearing) {
00517         thePixelThresholdInE = smearedThreshold_FPix_->fire(); // gaussian smearing
00518       } else {
00519         thePixelThresholdInE = theThresholdInE_FPix; // no smearing
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     // produce SignalPoint's for all SimHit's in detector
00537     // Loop over hits
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();  // Clear the container
00552       // fill _collection_points for this SimHit, indpendent of topology
00553       // Check the TOF cut
00554       //if (std::abs( (*ssbegin).tof() )<theTofCut){ // old cut
00555       if ( ((*ssbegin).tof() >= theTofLowerCut) && ((*ssbegin).tof() <= theTofUpperCut) ) {
00556         primary_ionization(*ssbegin); // fills _ionization_points       
00557         drift(*ssbegin);  // transforms _ionization_points to _collection_points        
00558         // compute induced signal on readout elements and add to _signal
00559         induce_signal(*ssbegin); // *ihit needed only for SimHit<-->Digi link
00560       } //  end if 
00561     } // end for 
00562   
00563     if(addNoise) add_noise();  // generate noise
00564     // Do only if needed 
00565 
00566     if((pixelInefficiency>0) && (_signal.size()>0)) 
00567       pixel_inefficiency(); // Kill some pixels
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_) // remove dead modules using DB
00575       module_killing_DB();
00576 
00577     if(use_module_killing_ && !use_deadmodule_DB_) // remove dead modules using the list in cfg file
00578       module_killing_conf();
00579 
00580   }
00581 
00582   make_digis();
00583 
00584   return internal_coll;
00585 }
00586 
00587 //***********************************************************************/
00588 // Generate primary ionization along the track segment. 
00589 // Divide the track into small sub-segments  
00590 void SiPixelDigitizerAlgorithm::primary_ionization(const PSimHit& hit) {
00591 
00592 // Straight line approximation for trajectory inside active media
00593 
00594   const float SegmentLength = 0.0010; //10microns in cm
00595   float energy;
00596   
00597   // Get the 3D segment direction vector 
00598   LocalVector direction = hit.exitPoint() - hit.entryPoint(); 
00599   
00600   float eLoss = hit.energyLoss();  // Eloss in GeV
00601   float length = direction.mag();  // Track length in Silicon
00602   
00603   NumberOfSegments = int ( length / SegmentLength); // Number of segments
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];  // Eloss vector
00617 
00618   if( fluctuateCharge ) {
00619     //MP DA RIMUOVERE ASSOLUTAMENTE
00620     int pid = hit.particleType();
00621     //int pid=211;  // assume it is a pion
00622   
00623     float momentum = hit.pabs();
00624     // Generate fluctuated charge points
00625     fluctuateEloss(pid, momentum, eLoss, length, NumberOfSegments, 
00626                    elossVector);
00627   }
00628   
00629   _ionization_points.resize( NumberOfSegments); // set size
00630 
00631   // loop over segments
00632   for ( int i = 0; i != NumberOfSegments; i++) {
00633     // Divide the segment into equal length subsegments 
00634     Local3DPoint point = hit.entryPoint() + 
00635       float((i+0.5)/NumberOfSegments) * direction;
00636 
00637     if( fluctuateCharge ) 
00638       energy = elossVector[i]/GeVperElectron; // Convert charge to elec.
00639     else
00640       energy = hit.energyLoss()/GeVperElectron/float(NumberOfSegments);
00641     
00642     EnergyDepositUnit edu( energy, point); //define position,energy point
00643     _ionization_points[i] = edu; // save
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   }  // end for loop 
00654 
00655   delete[] elossVector;
00656 
00657 }
00658 //****************************************************************************** 
00659 
00660 // Fluctuate the charge comming from a small (10um) track segment.
00661 // Use the G4 routine. For mip pions for the moment.  
00662 void SiPixelDigitizerAlgorithm::fluctuateEloss(int pid, float particleMomentum, 
00663                                       float eloss, float length, 
00664                                       int NumberOfSegs,float elossVector[]) {
00665 
00666   // Get dedx for this track
00667   float dedx;
00668   if( length > 0.) dedx = eloss/length;
00669   else dedx = eloss;
00670 
00671   double particleMass = 139.6; // Mass in MeV, Assume pion
00672   pid = abs(pid);
00673   if(pid!=211) {       // Mass in MeV
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   // What is the track segment length.
00680   float segmentLength = length/NumberOfSegs;
00681 
00682   // Generate charge fluctuations.
00683   float de=0.;
00684   float sum=0.;
00685   double segmentEloss = (1000.*eloss)/NumberOfSegs; //eloss in MeV
00686   for (int i=0;i<NumberOfSegs;i++) {
00687     //       material,*,   momentum,energy,*, *,  mass
00688     //myglandz_(14.,segmentLength,2.,2.,dedx,de,0.14);
00689     // The G4 routine needs momentum in MeV, mass in Mev, delta-cut in MeV,
00690     // track segment length in mm, segment eloss in MeV 
00691     // Returns fluctuated eloss in MeV
00692     double deltaCutoff = tMax; // the cutoff is sometimes redefined inside, so fix it.
00693     de = fluctuate->SampleFluctuations(double(particleMomentum*1000.),
00694                                       particleMass, deltaCutoff, 
00695                                       double(segmentLength*10.),
00696                                       segmentEloss )/1000.; //convert to GeV
00697    
00698     elossVector[i]=de;
00699     sum +=de;
00700   }
00701   
00702   if(sum>0.) {  // If fluctuations give eloss>0.
00703     // Rescale to the same total eloss
00704     float ratio = eloss/sum;
00705    
00706     for (int ii=0;ii<NumberOfSegs;ii++) elossVector[ii]= ratio*elossVector[ii];
00707   } else {  // If fluctuations gives 0 eloss
00708     float averageEloss = eloss/NumberOfSegs;
00709     for (int ii=0;ii<NumberOfSegs;ii++) elossVector[ii]= averageEloss; 
00710   }
00711   return;
00712 }
00713 
00714 //*******************************************************************************
00715 // Drift the charge segments to the sensor surface (collection plane)
00716 // Include the effect of E-field and B-field 
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()); // set size
00725   
00726   LocalVector driftDir=DriftDirection();  // get the charge drift direction
00727   if(driftDir.z() ==0.) {
00728     LogWarning("Magnetic field") << " pxlx: drift in z is zero ";
00729     return;
00730   }  
00731 
00732   // tangent of Lorentz angle
00733   //float TanLorenzAngleX = driftDir.x()/driftDir.z(); 
00734   //float TanLorenzAngleY = 0.; // force to 0, driftDir.y()/driftDir.z();
00735 
00736   float TanLorenzAngleX, TanLorenzAngleY,dir_z, CosLorenzAngleX,
00737     CosLorenzAngleY;
00738   if ( alpha2Order) {
00739     
00740     TanLorenzAngleX = driftDir.x(); // tangen of Lorentz angle
00741     TanLorenzAngleY = driftDir.y();
00742     dir_z = driftDir.z(); // The z drift direction
00743     CosLorenzAngleX = 1./sqrt(1.+TanLorenzAngleX*TanLorenzAngleX); //cosine
00744     CosLorenzAngleY = 1./sqrt(1.+TanLorenzAngleY*TanLorenzAngleY); //cosine;
00745     
00746   } else{
00747     
00748     TanLorenzAngleX = driftDir.x();
00749     TanLorenzAngleY = 0.; // force to 0, driftDir.y()/driftDir.z();
00750     dir_z = driftDir.z(); // The z drift direction
00751     CosLorenzAngleX = 1./sqrt(1.+TanLorenzAngleX*TanLorenzAngleX); //cosine to estimate the path length
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.;  // Charge spread 
00764   float Sigma_y = 1.;
00765   float DriftDistance; // Distance between charge generation and collection 
00766   float DriftLength;   // Actual Drift Lentgh
00767   float Sigma;
00768   
00769   
00770   for (unsigned int i = 0; i != _ionization_points.size(); i++) {
00771     
00772     float SegX, SegY, SegZ; // position
00773     SegX = _ionization_points[i].x();
00774     SegY = _ionization_points[i].y();
00775     SegZ = _ionization_points[i].z();
00776     
00777     // Distance from the collection plane
00778     //DriftDistance = (moduleThickness/2. + SegZ); // Drift to -z 
00779     // Include explixitely the E drift direction (for CMS dir_z=-1)
00780     DriftDistance = moduleThickness/2. - (dir_z * SegZ); // Drift to -z 
00781     
00782     //if( DriftDistance <= 0.) 
00783     //cout<<" <=0 "<<DriftDistance<<" "<<i<<" "<<SegZ<<" "<<dir_z<<" "
00784     //  <<SegX<<" "<<SegY<<" "<<(moduleThickness/2)<<" "
00785     //  <<_ionization_points[i].energy()<<" "
00786     //  <<hit.particleType()<<" "<<hit.pabs()<<" "<<hit.energyLoss()<<" "
00787     //  <<hit.entryPoint()<<" "<<hit.exitPoint()
00788     //  <<endl;
00789 
00790     if( DriftDistance < 0.) {
00791       DriftDistance = 0.;
00792     } else if ( DriftDistance > moduleThickness )
00793       DriftDistance = moduleThickness;
00794     
00795     // Assume full depletion now, partial depletion will come later.
00796     float XDriftDueToMagField = DriftDistance * TanLorenzAngleX;
00797     float YDriftDueToMagField = DriftDistance * TanLorenzAngleY;
00798     
00799     // Shift cloud center
00800     float CloudCenterX = SegX + XDriftDueToMagField;
00801     float CloudCenterY = SegY + YDriftDueToMagField;
00802 
00803     // Calculate how long is the charge drift path
00804     DriftLength = sqrt( DriftDistance*DriftDistance + 
00805                         XDriftDueToMagField*XDriftDueToMagField +
00806                         YDriftDueToMagField*YDriftDueToMagField );
00807 
00808     // What is the charge diffusion after this path
00809     Sigma = sqrt(DriftLength/Dist300) * Sigma0;
00810 
00811     // Project the diffusion sigma on the collection plane
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     // Load the Charge distribution parameters
00819     _collection_points[i] = (sp);
00820 
00821   } // loop over ionization points, i.
00822  
00823 } // end drift
00824 
00825 //*************************************************************************
00826 // Induce the signal on the collection plane of the active sensor area.
00827 void SiPixelDigitizerAlgorithm::induce_signal( const PSimHit& hit) {
00828 
00829   // X  - Rows, Left-Right, 160, (1.6cm)   for barrel
00830   // Y  - Columns, Down-Up, 416, (6.4cm)
00831 
00832 #ifdef TP_DEBUG
00833     LogDebug ("Pixel Digitizer") 
00834       << " enter induce_signal, " 
00835       << topol->pitch().first << " " << topol->pitch().second; //OK
00836 #endif
00837 
00838    // local map to store pixels hit by 1 Hit.      
00839    typedef map< int, float, less<int> > hit_map_type;
00840    hit_map_type hit_signal;
00841 
00842    // map to store pixel integrals in the x and in the y directions
00843    map<int, float, less<int> > x,y; 
00844    
00845    // Assign signals to readout channels and store sorted by channel number
00846    
00847    // Iterate over collection points on the collection plane
00848    for ( vector<SignalPoint>::const_iterator i=_collection_points.begin();
00849          i != _collection_points.end(); i++) {
00850      
00851      float CloudCenterX = i->position().x(); // Charge position in x
00852      float CloudCenterY = i->position().y(); //                 in y
00853      float SigmaX = i->sigma_x();            // Charge spread in x
00854      float SigmaY = i->sigma_y();            //               in y
00855      float Charge = i->amplitude();          // Charge amplitude
00856 
00857 
00858      //if(SigmaX==0 || SigmaY==0) {
00859      //cout<<SigmaX<<" "<<SigmaY
00860      //   << " cloud " << i->position().x() << " " << i->position().y() << " " 
00861      //   << i->sigma_x() << " " << i->sigma_y() << " " << i->amplitude()<<endl;
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      // Find the maximum cloud spread in 2D plane , assume 3*sigma
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      // Define 2D cloud limit points
00877      LocalPoint PointRightUp  = LocalPoint(CloudRight,CloudUp);
00878      LocalPoint PointLeftDown = LocalPoint(CloudLeft,CloudDown);
00879      
00880      // This points can be located outside the sensor area.
00881      // The conversion to measurement point does not check for that
00882      // so the returned pixel index might be wrong (outside range).
00883      // We rely on the limits check below to fix this.
00884      // But remember whatever we do here THE CHARGE OUTSIDE THE ACTIVE
00885      // PIXEL ARE IS LOST, it should not be collected.
00886 
00887      // Convert the 2D points to pixel indices
00888      MeasurementPoint mp = topol->measurementPosition(PointRightUp ); //OK
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 ); //OK
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      // Check detector limits to correct for pixels outside range.
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(); // clear temporary integration array
00917      y.clear();
00918 
00919      // First integrate cahrge strips in x
00920      int ix; // TT for compatibility
00921      for (ix=IPixLeftDownX; ix<=IPixRightUpX; ix++) {  // loop over x index
00922        float xUB, xLB, UpperBound, LowerBound;
00923       
00924        // Why is set to 0 if ix=0, does it meen that we accept charge 
00925        // outside the sensor? CHeck How it was done in ORCA? 
00926        //if (ix == 0) LowerBound = 0.;
00927        if (ix == 0 || SigmaX==0. )  // skip for surface segemnts 
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; // get strip
00952        x[ix] = TotalIntegrationRange; // save strip integral 
00953        //if(SigmaX==0 || SigmaY==0) 
00954        //cout<<TotalIntegrationRange<<" "<<ix<<endl;
00955 
00956      }
00957 
00958     // Now integarte strips in y
00959     int iy; // TT for compatibility
00960     for (iy=IPixLeftDownY; iy<=IPixRightUpY; iy++) { //loope over y ind  
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; // save strip integral
00989       //if(SigmaX==0 || SigmaY==0) 
00990       //cout<<TotalIntegrationRange<<" "<<iy<<endl;
00991     }       
00992 
00993     // Get the 2D charge integrals by folding x and y strips
00994     int chan;
00995     for (ix=IPixLeftDownX; ix<=IPixRightUpX; ix++) {  // loop over x index
00996       for (iy=IPixLeftDownY; iy<=IPixRightUpY; iy++) { //loope over y ind  
00997         
00998         float ChargeFraction = Charge*x[ix]*y[iy];
00999 
01000         if( ChargeFraction > 0. ) {
01001           chan = PixelDigi::pixelToChannel( ix, iy);  // Get index 
01002           // Load the amplitude                                          
01003           hit_signal[chan] += ChargeFraction;
01004         } // endif
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() << " "  // givex edge position
01017           << chan; // edge belongs to previous ?
01018 #endif
01019 
01020       } // endfor iy
01021     } //endfor ix
01022    
01023 
01024     // Test conversions (THIS IS FOR TESTING ONLY) comment-out.
01025     //     mp = topol->measurementPosition( i->position() ); //OK
01026     //     LocalPoint lp = topol->localPosition(mp);     //OK
01027     //     pair<float,float> p = topol->pixel( i->position() );  //OK
01028     //     chan = PixelDigi::pixelToChannel( int(p.first), int(p.second));
01029     //     pair<int,int> ip = PixelDigi::channelToPixel(chan);
01030     //     MeasurementPoint mp1 = MeasurementPoint( float(ip.first),
01031     //                                       float(ip.second) );
01032     //     LogDebug ("Pixel Digitizer") << " Test "<< mp.x() << " " << mp.y() 
01033     //                           << " "<< lp.x() << " " << lp.y() << " "<<" "
01034     //                           <<p.first <<" "<<p.second<<" "<<chan<< " "
01035     //                           <<" " << ip.first << " " << ip.second << " "
01036     //                           << mp1.x() << " " << mp1.y() << " " //OK
01037     //                           << topol->localPosition(mp1).x() << " "  //OK
01038     //                           << topol->localPosition(mp1).y() << " "
01039     //                           << topol->channel( i->position() ); //OK
01040     
01041     
01042   } // loop over charge distributions
01043    
01044   // Fill the global map with all hit pixels from this event
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 } // end induce_signal
01061 
01062 /***********************************************************************/
01063 
01064 // Build pixels, check threshold, add misscalibration, ... 
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   // Loop over hit pixels
01076 
01077   for ( signal_map_iterator i = _signal.begin(); i != _signal.end(); i++) {
01078     
01079         float signalInElectrons = (*i).second ;   // signal in electrons
01080 
01081     // Do the miss calibration for calibration studies only.
01082     //if(doMissCalibrate) signalInElectrons = missCalibrate(signalInElectrons)
01083 
01084     // Do only for pixels above threshold
01085 
01086        if ( signalInElectrons >= thePixelThresholdInE) { // check threshold
01087              
01088       int chan =  (*i).first;  // channel number
01089       pair<int,int> ip = PixelDigi::channelToPixel(chan);
01090       int adc=0;  // ADC count as integer
01091 
01092       // Do the miss calibration for calibration studies only.
01093       if(doMissCalibrate) {
01094         int row = ip.first;  // X in row
01095         int col = ip.second; // Y is in col
01096         adc = int(missCalibrate(col,row,signalInElectrons)); //full misscalib.
01097       } else { // Just do a simple electron->adc conversion
01098         adc = int( signalInElectrons / theElectronPerADC ); // calibrate gain
01099       }
01100       adc = min(adc, theAdcFullScale); // Check maximum value
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       // Load digis
01109       internal_coll.push_back( PixelDigi( ip.first, ip.second, adc));
01110   
01111       //digilink     
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         //sum the contribution of the same trackid 
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 //  Add electronic noise to pixel charge
01144 void SiPixelDigitizerAlgorithm::add_noise() {
01145 
01146 #ifdef TP_DEBUG
01147   LogDebug ("Pixel Digitizer") << " enter add_noise " << theNoiseInElectrons;
01148 #endif
01149  
01150   // First add noise to hit pixels
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         // Noise from Vcal smearing:
01166         float noise_ChargeVCALSmearing = theSmearedChargeRMS * gaussDistributionVCALNoise_->fire() ;
01167         // Noise from full readout:
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       } // End if addChargeVCalSmearing
01177          else
01178      {
01179         // Noise: ONLY full READOUT Noise.
01180         // Use here the FULL readout noise, including TBM,ALT,AOH,OPT-REC.
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      } // end if only Noise from full readout
01189     
01190   }
01191 
01192   if(!addNoisyPixels)  // Option to skip noise in non-hit pixels
01193     return;
01194 
01195   // Add noise on non-hit pixels
01196   // Use here the pixel noise 
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, //thr. in un. of nois
01203                       theNoiseInElectrons, // noise in elec. 
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   // Add noisy pixels
01215   for (mapI = otherPixels.begin(); mapI!= otherPixels.end(); mapI++) {
01216     int iy = ((*mapI).first) / numRows;
01217     int ix = ((*mapI).first) - (iy*numRows);
01218 
01219     // Keep for a while for testing.
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       //      float noise = float( (*mapI).second );
01235       int noise=int( (*mapI).second );
01236       _signal[chan] = Amplitude (noise, 0,-1.);
01237     }
01238   }
01239 
01240 }
01241 
01242 /***********************************************************************/
01243 
01244 // Simulate the readout inefficiencies. 
01245 // Delete a selected number of single pixels, dcols and rocs.
01246 void SiPixelDigitizerAlgorithm::pixel_inefficiency() {
01247 
01248 
01249   // Predefined efficiencies
01250   float pixelEfficiency  = 1.0;
01251   float columnEfficiency = 1.0;
01252   float chipEfficiency   = 1.0;
01253 
01254   // setup the chip indices conversion
01255   unsigned int Subid=DetId(detID).subdetId();
01256   if    (Subid==  PixelSubdetector::PixelBarrel){// barrel layers
01257     int layerIndex=PXBDetId(detID).layer();
01258     pixelEfficiency  = thePixelEfficiency[layerIndex-1];
01259     columnEfficiency = thePixelColEfficiency[layerIndex-1];
01260     chipEfficiency   = thePixelChipEfficiency[layerIndex-1];
01261  
01262     // This should never happen
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 {                // forward disks
01267    
01268     // For endcaps take same for each endcap
01269     pixelEfficiency  = thePixelEfficiency[3];
01270     columnEfficiency = thePixelColEfficiency[3];
01271     chipEfficiency   = thePixelChipEfficiency[3];
01272  
01273     // Sometimes the forward pixels have wrong size, 
01274     // this crashes the index conversion, so exit.
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   } // if barrel/forward
01281 
01282 #ifdef TP_DEBUG
01283   LogDebug ("Pixel Digitizer") << " enter pixel_inefficiency " << pixelEfficiency << " " 
01284                                << columnEfficiency << " " << chipEfficiency;
01285 #endif
01286   
01287   // Initilize the index converter
01288   //PixelIndices indexConverter(numColumns,numRows);
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   // Find out the number of columns and rocs hits
01296   // Loop over hit pixels, amplitude in electrons, channel = coded row,col
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;  // X in row
01302     int col = ip.second; // Y is in col
01303     //transform to ROC index coordinates   
01304     pIndexConverter->transformToROC(col,row,chipIndex,colROC,rowROC);
01305     int dColInChip = pIndexConverter->DColumn(colROC); // get ROC dcol from ROC col 
01306     //dcol in mod
01307     int dColInDet = pIndexConverter->DColumnInModule(dColInChip,chipIndex); 
01308       
01309     chips[chipIndex]++;
01310     columns[dColInDet]++;
01311   }
01312   
01313   // Delete some ROC hits.
01314   for ( iter = chips.begin(); iter != chips.end() ; iter++ ) {
01315     //float rand  = RandFlat::shoot();
01316     float rand  = flatDistribution_->fire();
01317     if( rand > chipEfficiency ) chips[iter->first]=0;
01318   }
01319 
01320   // Delete some Dcol hits.
01321   for ( iter = columns.begin(); iter != columns.end() ; iter++ ) {
01322     //float rand  = RandFlat::shoot();
01323     float rand  = flatDistribution_->fire();
01324     if( rand > columnEfficiency ) columns[iter->first]=0;
01325   }
01326   
01327   // Now loop again over pixels to kill some of them.
01328   // Loop over hit pixels, amplitude in electrons, channel = coded row,col
01329   for(signal_map_iterator i = _signal.begin();i != _signal.end(); i++) {    
01330 
01331     //    int chan = i->first;
01332     pair<int,int> ip = PixelDigi::channelToPixel(i->first);//get pixel pos
01333     int row = ip.first;  // X in row
01334     int col = ip.second; // Y is in col
01335     //transform to ROC index coordinates   
01336     pIndexConverter->transformToROC(col,row,chipIndex,colROC,rowROC);
01337     int dColInChip = pIndexConverter->DColumn(colROC); //get ROC dcol from ROC col 
01338     //dcol in mod
01339     int dColInDet = pIndexConverter->DColumnInModule(dColInChip,chipIndex); 
01340 
01341 
01342     //float rand  = RandFlat::shoot();
01343     float rand  = flatDistribution_->fire();
01344     if( chips[chipIndex]==0 || columns[dColInDet]==0 
01345         || rand>pixelEfficiency ) {
01346       // make pixel amplitude =0, pixel will be lost at clusterization    
01347       i->second.set(0.); // reset amplitude, 
01348     } // end if
01349 
01350   } // end pixel loop
01351   
01352 } // end pixel_indefficiency
01353 
01354 //***********************************************************************
01355 
01356 // Fluctuate the gain and offset for the amplitude calibration
01357 // Use gaussian smearing.
01358 //float SiPixelDigitizerAlgorithm::missCalibrate(const float amp) const {
01359   //float gain  = RandGaussQ::shoot(1.,theGainSmearing);
01360   //float offset  = RandGaussQ::shoot(0.,theOffsetSmearing);
01361   //float newAmp = amp * gain + offset;
01362   // More complex misscalibration 
01363 float SiPixelDigitizerAlgorithm::missCalibrate(int col,int row,
01364                                  const float signalInElectrons) const {
01365 
01366   // Central values
01367   //const float p0=0.00352, p1=0.868, p2=112., p3=113.; // pix(0,0,0)
01368   //  const float p0=0.00382, p1=0.886, p2=112.7, p3=113.0; // average roc=0
01369   //const float p0=0.00492, p1=1.998, p2=90.6, p3=134.1; // average roc=6
01370   // Smeared (rms)
01371   //const float s0=0.00020, s1=0.051, s2=5.4, s3=4.4; // average roc=0
01372   //const float s0=0.00015, s1=0.043, s2=3.2, s3=3.1; // col average roc=0
01373 
01374   // Make 2 sets of parameters for Fpix and BPIx:
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){// barrel layers
01384       p0 = BPix_p0;
01385       p1 = BPix_p1;
01386       p2 = BPix_p2;
01387       p3 = BPix_p3;
01388     } else {// forward disks
01389       p0 = FPix_p0;
01390       p1 = FPix_p1;
01391       p2 = FPix_p2;
01392       p3 = FPix_p3;
01393     }
01394 
01395   //  const float electronsPerVCAL = 65.5; // our present VCAL calibration (feb 2009)
01396   //  const float electronsPerVCAL_Offset = -414.0; // our present VCAL calibration (feb 2009)
01397   float newAmp = 0.; //Modified signal
01398 
01399   // Convert electrons to VCAL units
01400   float signal = (signalInElectrons-electronsPerVCAL_Offset)/electronsPerVCAL;
01401 
01402   // Simulate the analog response with fixed parametrization
01403   newAmp = p3 + p2 * tanh(p0*signal - p1);
01404 
01405 
01406   // Use the pixel-by-pixel calibrations
01407   //transform to ROC index coordinates
01408   //int chipIndex=0, colROC=0, rowROC=0;
01409   //pIndexConverter->transformToROC(col,row,chipIndex,colROC,rowROC);
01410 
01411   // Use calibration from a file
01412   //int chanROC = PixelIndices::pixelToChannelROC(rowROC,colROC); // use ROC coordinates
01413   //float pp0=0, pp1=0,pp2=0,pp3=0;
01414   //map<int,CalParameters,less<int> >::const_iterator it=calmap.find(chanROC);
01415   //CalParameters y  = (*it).second;
01416   //pp0 = y.p0;
01417   //pp1 = y.p1;
01418   //pp2 = y.p2;
01419   //pp3 = y.p3;
01420 
01421   //
01422   // Use random smearing 
01423   // Randomize the pixel response
01424   //float pp0  = RandGaussQ::shoot(p0,s0);
01425   //float pp1  = RandGaussQ::shoot(p1,s1);
01426   //float pp2  = RandGaussQ::shoot(p2,s2);
01427   //float pp3  = RandGaussQ::shoot(p3,s3);
01428 
01429   //newAmp = pp3 + pp2 * tanh(pp0*signal - pp1); // Final signal
01430 
01431   //cout<<" misscalibrate "<<col<<" "<<row<<" "<<chipIndex<<" "<<colROC<<" "
01432   //  <<rowROC<<" "<<signalInElectrons<<" "<<signal<<" "<<newAmp<<" "
01433   //  <<(signalInElectrons/theElectronPerADC)<<endl;
01434 
01435   return newAmp;
01436 }  
01437 //******************************************************************************
01438 
01439 // Set the drift direction accoring to the Bfield in local det-unit frame
01440 // Works for both barrel and forward pixels.
01441 // Replace the sign convention to fit M.Swartz's formulaes.
01442 // Configurations for barrel and foward pixels possess different tanLorentzAngleperTesla 
01443 // parameter value
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   //float dir_x = -tanLorentzAnglePerTesla * Bfield.y();
01454   //float dir_y = +tanLorentzAnglePerTesla * Bfield.x();
01455   //float dir_z = -1.; // E field always in z direction, so electrons go to -z
01456   // The dir_z has to be +/- 1. !
01457   // LocalVector theDriftDirection = LocalVector(dir_x,dir_y,dir_z);
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   // Read Lorentz angle from cfg file:**************************************************************
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){// barrel layers
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 {// forward disks
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     } // end: Read LA from cfg file.
01491   
01492   //Read Lorentz angle from DB:********************************************************************
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         //      std::cout << "detID is: " << it->first << "The LA per tesla is: " << it->second << std::endl;
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   }// end: Read LA from DataBase.
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   // Loop over hit pixels, amplitude in electrons, channel = coded row,col
01533   for(signal_map_iterator i = _signal.begin();i != _signal.end(); i++) {    
01534     
01535     //    int chan = i->first;
01536     pair<int,int> ip = PixelDigi::channelToPixel(i->first);//get pixel pos
01537     int row = ip.first;  // X in row
01538     int col = ip.second; // Y is in col
01539     uint32_t detid = detID;
01540     //transform to ROC index coordinates   
01541     if(theSiPixelGainCalibrationService_->isDead(detid, col, row)){
01542       //      std::cout << "now in isdead check, row " << detid << " " << col << "," << row << std::endl;
01543       // make pixel amplitude =0, pixel will be lost at clusterization    
01544       i->second.set(0.); // reset amplitude, 
01545     } // end if
01546   } // end pixel loop
01547 } // end pixel_indefficiency
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.); // reset amplitude
01577     }
01578   }
01579   
01580   for(signal_map_iterator i = _signal.begin();i != _signal.end(); i++) {    
01581     pair<int,int> ip = PixelDigi::channelToPixel(i->first);//get pixel pos
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){ // this is a whole dead module.
01621 
01622     for(signal_map_iterator i = _signal.begin();i != _signal.end(); i++) {
01623            i->second.set(0.); // reset amplitude
01624     }
01625   }
01626   else { // all other module types: half-modules and single ROCs.
01627            // Get Bad ROC position:
01628            //follow the example of getBadRocPositions in CondFormats/SiPixelObjects/src/SiPixelQuality.cc 
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};   //corresponding to center of ROC row, col
01640               GlobalPixel global = myroc->toGlobal( LocalPixel(local) );
01641               badrocpositions.push_back(global);
01642          break;
01643        }
01644       }
01645      }
01646    }// end of getBadRocPositions
01647 
01648   for(signal_map_iterator i = _signal.begin();i != _signal.end(); i++) {    
01649     pair<int,int> ip = PixelDigi::channelToPixel(i->first);//get pixel pos
01650 
01651   for(std::vector<GlobalPixel>::const_iterator it = badrocpositions.begin(); it != badrocpositions.end(); ++it){
01652     //  std::cout<<"getBadRocPositions, detidnumber: " << detid << "row is "<<it->row<<" col is "<<it->col<<std::endl;
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