Go to the documentation of this file.00001 #include "CalibTracker/SiStripCommon/interface/ShallowGainCalibration.h"
00002 #include "CalibFormats/SiStripObjects/interface/SiStripGain.h"
00003 #include "CalibTracker/Records/interface/SiStripGainRcd.h"
00004
00005 using namespace edm;
00006 using namespace reco;
00007 using namespace std;
00008
00009
00010
00011
00012 ShallowGainCalibration::ShallowGainCalibration(const edm::ParameterSet& iConfig)
00013 : theTracksLabel( iConfig.getParameter<edm::InputTag>("Tracks") ),
00014 Suffix ( iConfig.getParameter<std::string>("Suffix") ),
00015 Prefix ( iConfig.getParameter<std::string>("Prefix") )
00016 {
00017 produces <std::vector<int> > ( Prefix + "trackindex" + Suffix );
00018 produces <std::vector<unsigned int> > ( Prefix + "rawid" + Suffix );
00019 produces <std::vector<float> > ( Prefix + "localdirx" + Suffix );
00020 produces <std::vector<float> > ( Prefix + "localdiry" + Suffix );
00021 produces <std::vector<float> > ( Prefix + "localdirz" + Suffix );
00022 produces <std::vector<unsigned short> > ( Prefix + "firststrip" + Suffix );
00023 produces <std::vector<unsigned short> > ( Prefix + "nstrips" + Suffix );
00024 produces <std::vector<bool> > ( Prefix + "saturation" + Suffix );
00025 produces <std::vector<bool> > ( Prefix + "overlapping" + Suffix );
00026 produces <std::vector<bool> > ( Prefix + "farfromedge" + Suffix );
00027 produces <std::vector<unsigned int> > ( Prefix + "charge" + Suffix );
00028 produces <std::vector<float> > ( Prefix + "path" + Suffix );
00029 produces <std::vector<float> > ( Prefix + "chargeoverpath" + Suffix );
00030 produces <std::vector<unsigned char> > ( Prefix + "amplitude" + Suffix );
00031 produces <std::vector<double> > ( Prefix + "gainused" + Suffix );
00032 }
00033
00034 void ShallowGainCalibration::
00035 produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
00036 std::auto_ptr<std::vector<int> > trackindex ( new std::vector<int> );
00037 std::auto_ptr<std::vector<unsigned int> > rawid ( new std::vector<unsigned int> );
00038 std::auto_ptr<std::vector<float> > localdirx ( new std::vector<float> );
00039 std::auto_ptr<std::vector<float> > localdiry ( new std::vector<float> );
00040 std::auto_ptr<std::vector<float> > localdirz ( new std::vector<float> );
00041 std::auto_ptr<std::vector<unsigned short> > firststrip ( new std::vector<unsigned short> );
00042 std::auto_ptr<std::vector<unsigned short> > nstrips ( new std::vector<unsigned short> );
00043 std::auto_ptr<std::vector<bool> > saturation ( new std::vector<bool> );
00044 std::auto_ptr<std::vector<bool> > overlapping ( new std::vector<bool> );
00045 std::auto_ptr<std::vector<bool> > farfromedge ( new std::vector<bool> );
00046 std::auto_ptr<std::vector<unsigned int> > charge ( new std::vector<unsigned int> );
00047 std::auto_ptr<std::vector<float> > path ( new std::vector<float> );
00048 std::auto_ptr<std::vector<float> > chargeoverpath( new std::vector<float> );
00049 std::auto_ptr<std::vector<unsigned char> > amplitude ( new std::vector<unsigned char> );
00050 std::auto_ptr<std::vector<double> > gainused ( new std::vector<double> );
00051
00052 edm::ESHandle<TrackerGeometry> theTrackerGeometry; iSetup.get<TrackerDigiGeometryRecord>().get( theTrackerGeometry );
00053 m_tracker=&(* theTrackerGeometry );
00054 edm::ESHandle<SiStripGain> gainHandle; iSetup.get<SiStripGainRcd>().get(gainHandle);
00055 edm::Handle<edm::View<reco::Track> > tracks; iEvent.getByLabel(theTracksLabel, tracks);
00056 edm::Handle<TrajTrackAssociationCollection> associations; iEvent.getByLabel(theTracksLabel, associations);
00057
00058 for( TrajTrackAssociationCollection::const_iterator association = associations->begin(); association != associations->end(); association++) {
00059 const Trajectory* traj = association->key.get();
00060 const reco::Track* track = association->val.get();
00061
00062 vector<TrajectoryMeasurement> measurements = traj->measurements();
00063 for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
00064 TrajectoryStateOnSurface trajState = measurement_it->updatedState();
00065 if( !trajState.isValid() ) continue;
00066
00067 const TrackingRecHit* hit = (*measurement_it->recHit()).hit();
00068 const SiStripRecHit1D* sistripsimple1dhit = dynamic_cast<const SiStripRecHit1D*>(hit);
00069 const SiStripRecHit2D* sistripsimplehit = dynamic_cast<const SiStripRecHit2D*>(hit);
00070 const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
00071
00072 const SiStripCluster* Cluster = NULL;
00073
00074
00075 for(unsigned int h=0;h<2;h++){
00076 if(!sistripmatchedhit && h==1){
00077 continue;
00078 }else if(sistripmatchedhit && h==0){
00079 Cluster = (sistripmatchedhit->monoHit() ->cluster()).get();
00080 }else if(sistripmatchedhit && h==1){
00081 Cluster = (sistripmatchedhit->stereoHit()->cluster()).get();
00082 }else if(sistripsimplehit){
00083 Cluster = (sistripsimplehit->cluster()).get();
00084 }else if(sistripsimple1dhit){
00085 Cluster = (sistripsimple1dhit->cluster()).get();
00086 }else{
00087 continue;
00088 }
00089
00090 LocalVector trackDirection = trajState.localDirection();
00091 double cosine = trackDirection.z()/trackDirection.mag();
00092 const vector<uint8_t>& Ampls = Cluster->amplitudes();
00093 uint32_t DetId = Cluster->geographicalId();
00094 int FirstStrip = Cluster->firstStrip();
00095 int APVId = FirstStrip/128;
00096 bool Saturation = false;
00097 bool Overlapping = false;
00098 unsigned int Charge = 0;
00099 double Path = (10.0*thickness(DetId))/fabs(cosine);
00100 double PrevGain = -1;
00101
00102 if(gainHandle.isValid()){
00103 SiStripApvGain::Range detGainRange = gainHandle->getRange(DetId);
00104 PrevGain = *(detGainRange.first + APVId);
00105 }
00106
00107 for(unsigned int a=0;a<Ampls.size();a++){
00108 Charge+=Ampls[a];
00109 if(Ampls[a] >=254)Saturation =true;
00110 amplitude->push_back( Ampls[a] );
00111 }
00112 double ChargeOverPath = (double)Charge / Path ;
00113
00114 if(FirstStrip==0 )Overlapping=true;
00115 if(FirstStrip==128 )Overlapping=true;
00116 if(FirstStrip==256 )Overlapping=true;
00117 if(FirstStrip==384 )Overlapping=true;
00118 if(FirstStrip==512 )Overlapping=true;
00119 if(FirstStrip==640 )Overlapping=true;
00120
00121 if(FirstStrip<=127 && FirstStrip+Ampls.size()>127)Overlapping=true;
00122 if(FirstStrip<=255 && FirstStrip+Ampls.size()>255)Overlapping=true;
00123 if(FirstStrip<=383 && FirstStrip+Ampls.size()>383)Overlapping=true;
00124 if(FirstStrip<=511 && FirstStrip+Ampls.size()>511)Overlapping=true;
00125 if(FirstStrip<=639 && FirstStrip+Ampls.size()>639)Overlapping=true;
00126
00127 if(FirstStrip+Ampls.size()==127 )Overlapping=true;
00128 if(FirstStrip+Ampls.size()==255 )Overlapping=true;
00129 if(FirstStrip+Ampls.size()==383 )Overlapping=true;
00130 if(FirstStrip+Ampls.size()==511 )Overlapping=true;
00131 if(FirstStrip+Ampls.size()==639 )Overlapping=true;
00132 if(FirstStrip+Ampls.size()==767 )Overlapping=true;
00133
00134 trackindex ->push_back( shallow::findTrackIndex(tracks, track) );
00135 rawid ->push_back( DetId );
00136 localdirx ->push_back( trackDirection.x() );
00137 localdiry ->push_back( trackDirection.y() );
00138 localdirz ->push_back( trackDirection.z() );
00139 firststrip ->push_back( FirstStrip );
00140 nstrips ->push_back( Ampls.size() );
00141 saturation ->push_back( Saturation );
00142 overlapping ->push_back( Overlapping );
00143 farfromedge ->push_back( IsFarFromBorder(&trajState,DetId, &iSetup) );
00144 charge ->push_back( Charge );
00145 path ->push_back( Path );
00146 chargeoverpath->push_back( ChargeOverPath );
00147 gainused ->push_back( PrevGain );
00148 }
00149 }
00150 }
00151
00152 iEvent.put(trackindex, Prefix + "trackindex" + Suffix );
00153 iEvent.put(rawid , Prefix + "rawid" + Suffix );
00154 iEvent.put(localdirx , Prefix + "localdirx" + Suffix );
00155 iEvent.put(localdiry , Prefix + "localdiry" + Suffix );
00156 iEvent.put(localdirz , Prefix + "localdirz" + Suffix );
00157 iEvent.put(firststrip, Prefix + "firststrip" + Suffix );
00158 iEvent.put(nstrips, Prefix + "nstrips" + Suffix );
00159 iEvent.put(saturation, Prefix + "saturation" + Suffix );
00160 iEvent.put(overlapping, Prefix + "overlapping" + Suffix );
00161 iEvent.put(farfromedge, Prefix + "farfromedge" + Suffix );
00162 iEvent.put(charge, Prefix + "charge" + Suffix );
00163 iEvent.put(path, Prefix + "path" + Suffix );
00164 iEvent.put(chargeoverpath,Prefix + "chargeoverpath"+ Suffix );
00165 iEvent.put(amplitude, Prefix + "amplitude" + Suffix );
00166 iEvent.put(gainused, Prefix + "gainused" + Suffix );
00167 }
00168
00169
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224 bool ShallowGainCalibration::IsFarFromBorder(TrajectoryStateOnSurface* trajState, const uint32_t detid, const edm::EventSetup* iSetup)
00225 {
00226 edm::ESHandle<TrackerGeometry> tkGeom; iSetup->get<TrackerDigiGeometryRecord>().get( tkGeom );
00227
00228 LocalPoint HitLocalPos = trajState->localPosition();
00229 LocalError HitLocalError = trajState->localError().positionError() ;
00230
00231 const GeomDetUnit* it = tkGeom->idToDetUnit(DetId(detid));
00232 if (dynamic_cast<const StripGeomDetUnit*>(it)==0 && dynamic_cast<const PixelGeomDetUnit*>(it)==0) {
00233 std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
00234 return false;
00235 }
00236
00237 const BoundPlane plane = it->surface();
00238 const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
00239 const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
00240
00241 double DistFromBorder = 1.0;
00242 double HalfWidth = it->surface().bounds().width() /2.0;
00243 double HalfLength = it->surface().bounds().length() /2.0;
00244
00245 if(trapezoidalBounds)
00246 {
00247 std::vector<float> const & parameters = (*trapezoidalBounds).parameters();
00248 HalfLength = parameters[3];
00249 double t = (HalfLength + HitLocalPos.y()) / (2*HalfLength) ;
00250 HalfWidth = parameters[0] + (parameters[1]-parameters[0]) * t;
00251 }else if(rectangularBounds){
00252 HalfWidth = it->surface().bounds().width() /2.0;
00253 HalfLength = it->surface().bounds().length() /2.0;
00254 }else{return false;}
00255
00256
00257 if (fabs(HitLocalPos.y())+HitLocalError.yy() >= (HalfLength - DistFromBorder) ) return false;
00258
00259 return true;
00260 }
00261
00262
00263 double ShallowGainCalibration::thickness(DetId id)
00264 {
00265 map<DetId,double>::iterator th=m_thicknessMap.find(id);
00266 if(th!=m_thicknessMap.end())
00267 return (*th).second;
00268 else {
00269 double detThickness=1.;
00270
00271 const GeomDetUnit* it = m_tracker->idToDetUnit(DetId(id));
00272 bool isPixel = dynamic_cast<const PixelGeomDetUnit*>(it)!=0;
00273 bool isStrip = dynamic_cast<const StripGeomDetUnit*>(it)!=0;
00274 if (!isPixel && ! isStrip) {
00275
00276 edm::LogWarning("DeDxHitsProducer") << "\t\t this detID doesn't seem to belong to the Tracker";
00277 detThickness = 1.;
00278 }else{
00279 detThickness = it->surface().bounds().thickness();
00280 }
00281
00282 m_thicknessMap[id]=detThickness;
00283 return detThickness;
00284 }
00285
00286 }
00287