CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoTracker/DeDx/plugins/DeDxDiscriminatorLearner.cc

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    DeDxDiscriminatorLearner
00004 // Class:      DeDxDiscriminatorLearner
00005 // 
00013 //
00014 // Original Author:  Loic Quertenmont(querten)
00015 //         Created:  Mon October 20 10:09:02 CEST 2008
00016 //
00017 
00018 
00019 // system include files
00020 #include <memory>
00021 
00022 #include "RecoTracker/DeDx/plugins/DeDxDiscriminatorLearner.h"
00023 
00024 //#include "DataFormats/TrackReco/interface/DeDxData.h"
00025 #include "DataFormats/TrackReco/interface/Track.h"
00026 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
00027 
00028 #include "FWCore/Framework/interface/ESHandle.h"
00029 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00030 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00031 
00032 using namespace reco;
00033 using namespace std;
00034 using namespace edm;
00035 
00036 DeDxDiscriminatorLearner::DeDxDiscriminatorLearner(const edm::ParameterSet& iConfig) : ConditionDBWriter<PhysicsTools::Calibration::HistogramD3D>(iConfig)
00037 {
00038    m_tracksTag                 = iConfig.getParameter<edm::InputTag>("tracks");
00039    m_trajTrackAssociationTag   = iConfig.getParameter<edm::InputTag>("trajectoryTrackAssociation");
00040 
00041    usePixel = iConfig.getParameter<bool>("UsePixel"); 
00042    useStrip = iConfig.getParameter<bool>("UseStrip");
00043    if(!usePixel && !useStrip)
00044    edm::LogWarning("DeDxDiscriminatorLearner") << "Pixel Hits AND Strip Hits will not be used to estimate dEdx --> BUG, Please Update the config file";
00045 
00046    P_Min               = iConfig.getParameter<double>  ("P_Min"  );
00047    P_Max               = iConfig.getParameter<double>  ("P_Max"  );
00048    P_NBins             = iConfig.getParameter<int>     ("P_NBins");
00049    Path_Min            = iConfig.getParameter<double>  ("Path_Min"  );
00050    Path_Max            = iConfig.getParameter<double>  ("Path_Max"  );
00051    Path_NBins          = iConfig.getParameter<int>     ("Path_NBins");
00052    Charge_Min          = iConfig.getParameter<double>  ("Charge_Min"  );
00053    Charge_Max          = iConfig.getParameter<double>  ("Charge_Max"  );
00054    Charge_NBins        = iConfig.getParameter<int>     ("Charge_NBins");
00055 
00056    MinTrackMomentum    = iConfig.getUntrackedParameter<double>  ("minTrackMomentum"   ,  3.0);
00057    MaxTrackMomentum    = iConfig.getUntrackedParameter<double>  ("maxTrackMomentum"   ,  99999.0); 
00058    MinTrackEta         = iConfig.getUntrackedParameter<double>  ("minTrackEta"        , -5.0);
00059    MaxTrackEta         = iConfig.getUntrackedParameter<double>  ("maxTrackEta"        ,  5.0);
00060    MaxNrStrips         = iConfig.getUntrackedParameter<unsigned>("maxNrStrips"        ,  255);
00061    MinTrackHits        = iConfig.getUntrackedParameter<unsigned>("MinTrackHits"       ,  4);
00062 
00063    algoMode            = iConfig.getUntrackedParameter<string>  ("AlgoMode"           ,  "SingleJob");
00064    HistoFile           = iConfig.getUntrackedParameter<string>  ("HistoFile"        ,  "out.root");
00065 }
00066 
00067 
00068 DeDxDiscriminatorLearner::~DeDxDiscriminatorLearner(){}
00069 
00070 // ------------ method called once each job just before starting event loop  ------------
00071 
00072 void  DeDxDiscriminatorLearner::algoBeginJob(const edm::EventSetup& iSetup)
00073 {
00074 //   Charge_Vs_Path = new TH2F ("Charge_Vs_Path"     , "Charge_Vs_Path" , 24, 0.2, 1.4, 250, 0, 5000);
00075    Charge_Vs_Path = new TH3F ("Charge_Vs_Path"     , "Charge_Vs_Path" , P_NBins, P_Min, P_Max, Path_NBins, Path_Min, Path_Max, Charge_NBins, Charge_Min, Charge_Max);
00076 
00077 
00078    edm::ESHandle<TrackerGeometry> tkGeom;
00079    iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
00080    m_tracker = tkGeom.product();
00081 
00082    vector<GeomDet*> Det = tkGeom->dets();
00083    for(unsigned int i=0;i<Det.size();i++){
00084       DetId  Detid  = Det[i]->geographicalId();
00085       int    SubDet = Detid.subdetId();
00086 
00087       if( SubDet == StripSubdetector::TIB ||  SubDet == StripSubdetector::TID ||
00088           SubDet == StripSubdetector::TOB ||  SubDet == StripSubdetector::TEC  ){
00089 
00090           StripGeomDetUnit* DetUnit     = dynamic_cast<StripGeomDetUnit*> (Det[i]);
00091           if(!DetUnit)continue;
00092 
00093           const StripTopology& Topo     = DetUnit->specificTopology();
00094           unsigned int         NAPV     = Topo.nstrips()/128;
00095 
00096           double Eta           = DetUnit->position().basicVector().eta();
00097           double R             = DetUnit->position().basicVector().transverse();
00098           double Thick         = DetUnit->surface().bounds().thickness();
00099 
00100           stModInfo* MOD       = new stModInfo;
00101           MOD->DetId           = Detid.rawId();
00102           MOD->SubDet          = SubDet;
00103           MOD->Eta             = Eta;
00104           MOD->R               = R;
00105           MOD->Thickness       = Thick;
00106           MOD->NAPV            = NAPV;
00107           MODsColl[MOD->DetId] = MOD;
00108       }
00109    }
00110 
00111 }
00112 
00113 // ------------ method called once each job just after ending the event loop  ------------
00114 
00115 void DeDxDiscriminatorLearner::algoEndJob()
00116 {
00117 
00118    if( strcmp(algoMode.c_str(),"MultiJob")==0){
00119         TFile* Output = new TFile(HistoFile.c_str(), "RECREATE");
00120         Charge_Vs_Path->Write();
00121         Output->Write();
00122         Output->Close();
00123    }else if( strcmp(algoMode.c_str(),"WriteOnDB")==0){
00124         TFile* Input = new TFile(HistoFile.c_str() );
00125         Charge_Vs_Path = (TH3F*)(Input->FindObjectAny("Charge_Vs_Path"))->Clone();  
00126         Input->Close();
00127    }
00128 
00129 }
00130 
00131 void DeDxDiscriminatorLearner::algoAnalyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
00132 {
00133    edm::ESHandle<TrackerGeometry> tkGeom;
00134    iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
00135    m_tracker = tkGeom.product();
00136 
00137 
00138 
00139   Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
00140   iEvent.getByLabel(m_trajTrackAssociationTag, trajTrackAssociationHandle);
00141   const TrajTrackAssociationCollection TrajToTrackMap = *trajTrackAssociationHandle.product();
00142 
00143   edm::Handle<reco::TrackCollection> trackCollectionHandle;
00144   iEvent.getByLabel(m_tracksTag,trackCollectionHandle);
00145 
00146   unsigned track_index = 0;
00147   for(TrajTrackAssociationCollection::const_iterator it = TrajToTrackMap.begin(); it!=TrajToTrackMap.end(); ++it, track_index++) {
00148 
00149       const Track      track = *it->val;
00150       const Trajectory traj  = *it->key;
00151 
00152       if(track.eta()  <MinTrackEta      || track.eta()>MaxTrackEta     ){continue;}
00153       if(track.p()    <MinTrackMomentum || track.p()  >MaxTrackMomentum){continue;}
00154       if(track.found()<MinTrackHits                                    ){continue;}
00155 
00156       vector<TrajectoryMeasurement> measurements = traj.measurements();
00157       for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++)
00158       {
00159          TrajectoryStateOnSurface trajState = measurement_it->updatedState();
00160          if( !trajState.isValid() ) continue;
00161 
00162          const TrackingRecHit*         hit               = (*measurement_it->recHit()).hit();
00163          const SiStripRecHit2D*        sistripsimplehit  = dynamic_cast<const SiStripRecHit2D*>(hit);
00164          const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
00165          const SiStripRecHit1D*        sistripsimple1dhit  = dynamic_cast<const SiStripRecHit1D*>(hit);
00166 
00167          if(sistripsimplehit){
00168              Learn((sistripsimplehit->cluster()).get(), trajState);
00169          }else if(sistripmatchedhit){
00170              Learn(&sistripmatchedhit->monoCluster(),trajState);
00171              Learn(&sistripmatchedhit->stereoCluster(),trajState);
00172          }else if(sistripsimple1dhit){
00173              Learn((sistripsimple1dhit->cluster()).get(), trajState);
00174          }else{
00175          }
00176 
00177       }
00178    }
00179 
00180 }
00181 
00182 
00183 void DeDxDiscriminatorLearner::Learn(const SiStripCluster*   cluster,TrajectoryStateOnSurface trajState)
00184 {
00185    // Get All needed variables
00186    LocalVector             trackDirection = trajState.localDirection();
00187    double                  cosine         = trackDirection.z()/trackDirection.mag();
00188    const vector<uint8_t>&  ampls          = cluster->amplitudes();
00189    uint32_t                detId          = cluster->geographicalId();
00190    int                     firstStrip     = cluster->firstStrip();
00191    stModInfo* MOD                         = MODsColl[detId];
00192    // Sanity Checks
00193    if( ampls.size()>MaxNrStrips)                                                                      { return; }
00194 // if( DeDxDiscriminatorTools::IsSaturatingStrip  (ampls))                                            { return; }
00195    if( DeDxDiscriminatorTools::IsSpanningOver2APV (firstStrip, ampls.size()))                         { return; }
00196    if(!DeDxDiscriminatorTools::IsFarFromBorder    (trajState, m_tracker->idToDetUnit(DetId(detId)) )) { return; }
00197    // Fill Histo Path Vs Charge/Path
00198    double charge = DeDxDiscriminatorTools::charge(ampls);
00199    double path   = DeDxDiscriminatorTools::path(cosine,MOD->Thickness);
00200    Charge_Vs_Path->Fill(trajState.localMomentum().mag(),path,charge/path);
00201 }
00202 
00203 
00204 PhysicsTools::Calibration::HistogramD3D* DeDxDiscriminatorLearner::getNewObject()
00205 {
00206 //   if( strcmp(algoMode.c_str(),"MultiJob")==0)return NULL;
00207 
00208    PhysicsTools::Calibration::HistogramD3D* obj;
00209    obj = new PhysicsTools::Calibration::HistogramD3D(
00210                 Charge_Vs_Path->GetNbinsX(), Charge_Vs_Path->GetXaxis()->GetXmin(),  Charge_Vs_Path->GetXaxis()->GetXmax(),
00211                 Charge_Vs_Path->GetNbinsY(), Charge_Vs_Path->GetYaxis()->GetXmin(),  Charge_Vs_Path->GetYaxis()->GetXmax(),
00212                 Charge_Vs_Path->GetNbinsZ(), Charge_Vs_Path->GetZaxis()->GetXmin(),  Charge_Vs_Path->GetZaxis()->GetXmax());
00213 
00214    for(int ix=0; ix<=Charge_Vs_Path->GetNbinsX()+1; ix++){
00215       for(int iy=0; iy<=Charge_Vs_Path->GetNbinsY()+1; iy++){
00216          for(int iz=0; iz<=Charge_Vs_Path->GetNbinsZ()+1; iz++){
00217             obj->setBinContent(ix, iy, iz, Charge_Vs_Path->GetBinContent(ix,iy, iz) );       
00218 //          if(Charge_Vs_Path->GetBinContent(ix,iy)!=0)printf("%i %i %i --> %f\n",ix,iy, iz, Charge_Vs_Path->GetBinContent(ix,iy,iz)); 
00219          }
00220       }
00221    }
00222 
00223    return obj;
00224 }
00225 
00226 
00227 
00228 //define this as a plug-in
00229 DEFINE_FWK_MODULE(DeDxDiscriminatorLearner);