CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/CalibTracker/SiStripCommon/plugins/ShallowGainCalibration.cc

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 void ShallowGainCalibration::beginJob(const edm::EventSetup& iSetup)
00171 {
00172    printf("Befin JOB\n");
00173 
00174    edm::ESHandle<TrackerGeometry> tkGeom;
00175    iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
00176    vector<GeomDet*> Det = tkGeom->dets();
00177 
00178    edm::ESHandle<SiStripGain> gainHandle;
00179    iSetup.get<SiStripGainRcd>().get(gainHandle);
00180    if(!gainHandle.isValid()){printf("\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");exit(0);}
00181 
00182    for(unsigned int i=0;i<Det.size();i++){
00183       DetId  Detid  = Det[i]->geographicalId();
00184       int    SubDet = Detid.subdetId();
00185 
00186       if( SubDet == StripSubdetector::TIB ||  SubDet == StripSubdetector::TID ||
00187           SubDet == StripSubdetector::TOB ||  SubDet == StripSubdetector::TEC  ){
00188 
00189           StripGeomDetUnit* DetUnit     = dynamic_cast<StripGeomDetUnit*> (Det[i]);
00190           if(!DetUnit)continue;
00191 
00192           const StripTopology& Topo     = DetUnit->specificTopology();
00193           unsigned int         NAPV     = Topo.nstrips()/128;
00194    
00195           for(unsigned int j=0;j<NAPV;j++){
00196                 stAPVGain* APV = new stAPVGain;
00197                 APV->DetId         = Detid.rawId();
00198                 APV->APVId         = j;
00199                 APV->PreviousGain  = 1;
00200 
00201                 APVsCollOrdered.push_back(APV);
00202                 APVsColl[(APV->DetId<<3) | APV->APVId] = APV;
00203           }
00204       }
00205    }
00206 }
00207 
00208 
00209 void ShallowGainCalibration::beginRun(edm::Run &, const edm::EventSetup &iSetup){
00210     printf("BEFIN RUN\n");
00211 
00212     edm::ESHandle<SiStripGain> gainHandle;
00213     iSetup.get<SiStripGainRcd>().get(gainHandle);
00214     if(!gainHandle.isValid()){printf("\n#####################\n\nERROR --> gainHandle is not valid\n\n#####################\n\n");exit(0);}
00215 
00216     for(std::vector<stAPVGain*>::iterator it = APVsCollOrdered.begin();it!=APVsCollOrdered.end();it++){
00217        stAPVGain* APV = *it;
00218        SiStripApvGain::Range detGainRange = gainHandle->getRange(APV->DetId);
00219        APV->PreviousGain = *(detGainRange.first + APV->APVId);
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 //  if (fabs(HitLocalPos.x())+HitLocalError.xx() >= (HalfWidth  - DistFromBorder) ) return false;//Don't think is really necessary
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    //compute thickness normalization
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    //FIXME throw exception
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;//computed value
00283    return detThickness;
00284  }
00285 
00286 }
00287