CMSSW Reference Manual

00001 // -*- C++ -*-
00002 //
00003 // Package:    DeDxDiscriminatorProducer
00004 // Class:      DeDxDiscriminatorProducer
00005 // 
00013 //
00014 // Original Author:  andrea
00015 //         Created:  Thu May 31 14:09:02 CEST 2007
00016 //    Code Updates:  loic Quertenmont (querten)
00017 //         Created:  Thu May 10 14:09:02 CEST 2008
00018 // $Id: DeDxDiscriminatorProducer.cc,v 1.3 2008/09/15 09:43:12 querten Exp $
00019 //
00020 //
00021 
00022 
00023 // system include files
00024 #include <memory>
00025 #include "DataFormats/Common/interface/ValueMap.h"
00026 
00027 #include "RecoTracker/DeDx/interface/DeDxDiscriminatorProducer.h"
00028 #include "DataFormats/TrackReco/interface/DeDxData.h"
00029 #include "DataFormats/TrackReco/interface/TrackDeDxHits.h"
00030 #include "DataFormats/TrackReco/interface/DeDxHit.h"
00031 #include "DataFormats/TrackReco/interface/Track.h"
00032 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
00033 
00034 #include "RecoTracker/DeDx/interface/GenericAverageDeDxEstimator.h"
00035 #include "RecoTracker/DeDx/interface/TruncatedAverageDeDxEstimator.h"
00036 #include "RecoTracker/DeDx/interface/MedianDeDxEstimator.h"
00037 
00038 #include "FWCore/Framework/interface/ESHandle.h"
00039 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00040 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00041 
00042 using namespace reco;
00043 using namespace std;
00044 using namespace edm;
00045 
00046 DeDxDiscriminatorProducer::DeDxDiscriminatorProducer(const edm::ParameterSet& iConfig)
00047 {
00048 
00049    produces<ValueMap<DeDxData> >();
00050 
00051    m_tracksTag = iConfig.getParameter<edm::InputTag>("tracks");
00052    m_trajTrackAssociationTag   = iConfig.getParameter<edm::InputTag>("trajectoryTrackAssociation");
00053 
00054    usePixel = iConfig.getParameter<bool>("UsePixel"); 
00055    useStrip = iConfig.getParameter<bool>("UseStrip");
00056    if(!usePixel && !useStrip)
00057    edm::LogWarning("DeDxHitsProducer") << "Pixel Hits AND Strip Hits will not be used to estimate dEdx --> BUG, Please Update the config file";
00058 
00059    DiscriminatorMode   = iConfig.getUntrackedParameter<bool>("DiscriminatorMode", true);
00060    MapFileName         = iConfig.getParameter<std::string>("MapFile");
00061    Formula             = iConfig.getUntrackedParameter<unsigned>("Formula"            ,  0);
00062 
00063    MinTrackMomentum    = iConfig.getUntrackedParameter<double>  ("minTrackMomentum"   ,  3.0);
00064    MaxTrackMomentum    = iConfig.getUntrackedParameter<double>  ("maxTrackMomentum"   ,  99999.0); 
00065    MinTrackEta         = iConfig.getUntrackedParameter<double>  ("minTrackEta"        , -5.0);
00066    MaxTrackEta         = iConfig.getUntrackedParameter<double>  ("maxTrackEta"        ,  5.0);
00067    MaxNrStrips         = iConfig.getUntrackedParameter<unsigned>("maxNrStrips"        ,  2);
00068    MinTrackHits        = iConfig.getUntrackedParameter<unsigned>("MinTrackHits"       ,  8);
00069    AllowSaturation     = iConfig.getUntrackedParameter<bool>    ("AllowSaturation"    ,  false);
00070 }
00071 
00072 
00073 DeDxDiscriminatorProducer::~DeDxDiscriminatorProducer(){}
00074 
00075 // ------------ method called once each job just before starting event loop  ------------
00076 void  DeDxDiscriminatorProducer::beginJob(const edm::EventSetup& iSetup){
00077    TH1::AddDirectory(kTRUE);
00078 
00079    iSetup_                  = &iSetup;
00080 
00081    if(!DiscriminatorMode){
00082       MapFile                   = new TFile(MapFileName.c_str(), "RECREATE");
00083       Charge_Vs_Path_Barrel     = new TH2F ("Charge_Vs_Path_Barrel"     , "Charge_Vs_Path_Barrel" , 250, 0.2, 1.4, 1000, 0, 5000);
00084       Charge_Vs_Path_Endcap     = new TH2F ("Charge_Vs_Path_Endcap"     , "Charge_Vs_Path_Endcap" , 250, 0.2, 1.4, 1000, 0, 5000);
00085    }else{
00086       MapFile                   = new TFile(MapFileName.c_str());
00087 
00088       Charge_Vs_Path_Barrel     = (TH2F*) MapFile->FindObjectAny("Charge_Vs_Path_Barrel");
00089       Charge_Vs_Path_Endcap     = (TH2F*) MapFile->FindObjectAny("Charge_Vs_Path_Endcap");
00090 
00091       PCharge_Vs_Path_Barrel    = new TH2F ("PCharge_Vs_Path_Barrel"     , "PCharge_Vs_Path_Barrel" , 250, 0.2, 1.4, 1000, 0, 5000);
00092       PCharge_Vs_Path_Endcap    = new TH2F ("PCharge_Vs_Path_Endcap"     , "PCharge_Vs_Path_Endcap" , 250, 0.2, 1.4, 1000, 0, 5000);
00093 
00094       TH1D* Proj;
00095       for(int i=0;i<Charge_Vs_Path_Barrel->GetXaxis()->GetNbins();i++){
00096          Proj      = Charge_Vs_Path_Barrel->ProjectionY(" ",i,i,"e");
00097          for(int j=0;j<Proj->GetXaxis()->GetNbins();j++){
00098             float Prob;
00099             if(Proj->GetEntries()>0){ Prob = Proj->Integral(0,j) / Proj->GetEntries();
00100             }else{                    Prob = 0;}
00101             PCharge_Vs_Path_Barrel->SetBinContent (i, j, Prob);
00102          }
00103          delete Proj;
00104       }
00105 
00106       for(int i=0;i<Charge_Vs_Path_Endcap->GetXaxis()->GetNbins();i++){
00107          Proj      = Charge_Vs_Path_Endcap->ProjectionY(" ",i,i,"e");
00108          for(int j=0;j<Proj->GetXaxis()->GetNbins();j++){
00109             float Prob;
00110             if(Proj->GetEntries()>0){ Prob = Proj->Integral(0,j) / Proj->GetEntries();
00111             }else{                    Prob = 0;}
00112             PCharge_Vs_Path_Endcap->SetBinContent (i, j, Prob);
00113          }
00114          delete Proj;
00115       }
00116 
00117    }
00118 
00119 
00120    edm::ESHandle<TrackerGeometry> tkGeom;
00121    iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
00122    m_tracker=&(* tkGeom );
00123 
00124    vector<GeomDet*> Det = tkGeom->dets();
00125    for(unsigned int i=0;i<Det.size();i++){
00126       DetId  Detid  = Det[i]->geographicalId();
00127       int    SubDet = Detid.subdetId();
00128 
00129       if( SubDet == StripSubdetector::TIB ||  SubDet == StripSubdetector::TID ||
00130           SubDet == StripSubdetector::TOB ||  SubDet == StripSubdetector::TEC  ){
00131 
00132           StripGeomDetUnit* DetUnit     = dynamic_cast<StripGeomDetUnit*> (Det[i]);
00133           if(!DetUnit)continue;
00134 
00135           const StripTopology& Topo     = DetUnit->specificTopology();
00136           unsigned int         NAPV     = Topo.nstrips()/128;
00137 
00138           double Eta     = DetUnit->position().basicVector().eta();
00139           double R       = DetUnit->position().basicVector().transverse();
00140           double Thick   = DetUnit->surface().bounds().thickness();
00141 
00142           stModInfo* MOD = new stModInfo;
00143           MOD->DetId     = Detid.rawId();
00144           MOD->SubDet    = SubDet;
00145           MOD->Eta       = Eta;
00146           MOD->R         = R;
00147           MOD->Thickness = Thick;
00148           MOD->NAPV      = NAPV;
00149           MODsColl[MOD->DetId] = MOD;
00150       }
00151    }
00152 
00153 
00154 }
00155 
00156 // ------------ method called once each job just after ending the event loop  ------------
00157 void  DeDxDiscriminatorProducer::endJob(){
00158    if(!DiscriminatorMode) MapFile->Write();
00159    MapFile->Close();
00160 }
00161 
00162 
00163 
00164 void DeDxDiscriminatorProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00165 {
00166   iEvent_ = &iEvent;  
00167 
00168   auto_ptr<ValueMap<DeDxData> > trackDeDxDiscrimAssociation(new ValueMap<DeDxData> );  
00169   ValueMap<DeDxData>::Filler filler(*trackDeDxDiscrimAssociation);
00170 
00171   Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
00172   iEvent.getByLabel(m_trajTrackAssociationTag, trajTrackAssociationHandle);
00173   const TrajTrackAssociationCollection TrajToTrackMap = *trajTrackAssociationHandle.product();
00174 
00175   edm::Handle<reco::TrackCollection> trackCollectionHandle;
00176   iEvent.getByLabel(m_tracksTag,trackCollectionHandle);
00177  
00178    std::vector<DeDxData> dEdxDiscrims( TrajToTrackMap.size() );
00179 
00180 
00181    cout << "DDP  TEST1\n";
00182 
00183 
00184    unsigned track_index = 0;
00185    for(TrajTrackAssociationCollection::const_iterator it = TrajToTrackMap.begin(); it!=TrajToTrackMap.end(); ++it, track_index++) {
00186       dEdxDiscrims[track_index] = DeDxData(-1, -1, 0 );
00187 
00188       const Track      track = *it->val;
00189       const Trajectory traj  = *it->key;
00190 
00191 /*
00192       if(OnlyHSCP){
00193          bool GoodTrack = false;
00194          HepMC::GenEvent * myGenEvent = new  HepMC::GenEvent(*(HepMC->GetEvent()));
00195          for ( HepMC::GenEvent::particle_iterator p = myGenEvent->particles_begin(); p != myGenEvent->particles_end(); ++p ) {
00196             if ( abs((*p)->pdg_id()) > 100000 && (*p)->status() == 3 ){
00197                 double Dist = deltaR(track, (*p)->momentum());
00198                 if(Dist<0.3)GoodTrack = true;
00199 
00200             }
00201          }
00202         if(!GoodTrack){printf("HSCPCut\n");continue;}
00203       }
00204 */
00205 
00206    cout << "DDP  TEST2\n";
00207 
00208 
00209       if(track.eta()<MinTrackEta || track.eta()>MaxTrackEta){printf("Eta Cut\n");continue;}
00210       if(track.p()<MinTrackMomentum || track.p()>MaxTrackMomentum){printf("Pt Cut\n");continue;}
00211       if(track.found()<MinTrackHits){printf("Hits Cut\n");continue;}
00212 
00213 
00214       vector<TrajectoryMeasurement> measurements = traj.measurements();
00215       if(traj.foundHits()<(int)MinTrackHits)continue;
00216 
00217       MeasurementProbabilities.clear();
00218       for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
00219 
00220          TrajectoryStateOnSurface trajState = measurement_it->updatedState();
00221          if( !trajState.isValid() ) continue;
00222 
00223          const TrackingRecHit*         hit               = (*measurement_it->recHit()).hit();
00224          const SiStripRecHit2D*        sistripsimplehit  = dynamic_cast<const SiStripRecHit2D*>(hit);
00225          const SiStripMatchedRecHit2D* sistripmatchedhit = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
00226 
00227          if(sistripsimplehit)
00228          {
00229              ComputeChargeOverPath(sistripsimplehit, trajState, &iSetup, &track, traj.chiSquared());
00230          }else if(sistripmatchedhit){
00231              ComputeChargeOverPath(sistripmatchedhit->monoHit()  ,trajState, &iSetup, &track, traj.chiSquared());
00232              ComputeChargeOverPath(sistripmatchedhit->stereoHit(),trajState, &iSetup, &track, traj.chiSquared());
00233          }else{
00234          }
00235       }
00236 
00237    cout << "DDP  TEST3\n";
00238 
00239 
00240       if(DiscriminatorMode){
00241          int size = MeasurementProbabilities.size();
00242 
00243          double estimator = 1;
00244          if(Formula==0){
00245             double P = 1;
00246             for(int i=0;i<size;i++){
00247                if(MeasurementProbabilities[i]<=0.001){P *= pow(0.001f, 1.0f/size);}
00248                else                                  {P *= pow(MeasurementProbabilities[i], 1.0f/size);}
00249             }
00250             estimator = P;
00251          }else if(Formula==1){
00252 
00253             if(MeasurementProbabilities.size()>0){
00254 
00255               std::sort(MeasurementProbabilities.begin(), MeasurementProbabilities.end(), std::less<double>() );
00256               for(int i=0;i<size;i++){if(MeasurementProbabilities[i]<=0.001)MeasurementProbabilities[i] = 0.001f;    }
00257 
00258                double SumJet = 0.;
00259                for(int i=0;i<size;i++){ SumJet+= log(MeasurementProbabilities[i]); }
00260 
00261               double Loginvlog=log(-SumJet);
00262               double Prob =1.;
00263               double lfact=1.;
00264 
00265               for(int l=1; l!=size; l++){
00266                    lfact*=l;
00267                    Prob+=exp(l*Loginvlog-log(1.*lfact));
00268                }
00269 
00270                double LogProb=log(Prob);
00271                double ProbJet=std::min(exp(std::max(LogProb+SumJet,-30.)),1.);
00272                estimator = -log10(ProbJet)/4.;
00273                estimator = 1-estimator;
00274             }else{
00275                estimator = -1;
00276             }
00277          }else if(Formula==2){
00278            estimator = -2;
00279            if(size>0){
00280                std::sort(MeasurementProbabilities.begin(), MeasurementProbabilities.end(), std::less<double>() );
00281                double P = 1.0/(12*size);
00282                for(int i=1;i<=size;i++){
00283                   P += pow(MeasurementProbabilities[i-1] - ((2.0*i-1.0)/(2.0*size)),2);
00284                }
00285                P *= (1.0/size);
00286                estimator = P;
00287                if(estimator>=0.333)printf("BUG\n");
00288             }
00289          }else{
00290            estimator = -2;
00291            if(size>0){
00292                std::sort(MeasurementProbabilities.begin(), MeasurementProbabilities.end(), std::less<double>() );
00293                double P = 1.0/(12*size);
00294                for(int i=1;i<=size;i++){
00295                   P += MeasurementProbabilities[i-1] * pow(MeasurementProbabilities[i-1] - ((2.0*i-1.0)/(2.0*size)),2);
00296                }
00297                P *= (1.0/size);
00298                estimator = P;
00299                if(estimator>=0.333)printf("BUG\n");
00300            }
00301          }
00302 
00303          dEdxDiscrims[track_index] = DeDxData(estimator, -1, size );
00304       }
00305    }
00306 
00307   filler.insert(trackCollectionHandle, dEdxDiscrims.begin(), dEdxDiscrims.end());
00308   filler.fill();
00309   iEvent.put(trackDeDxDiscrimAssociation);
00310 }
00311 
00312 
00313 double
00314 DeDxDiscriminatorProducer::ComputeChargeOverPath(const SiStripRecHit2D* sistripsimplehit,TrajectoryStateOnSurface trajState, const edm::EventSetup* iSetup,  const Track* track, double trajChi2OverN)
00315 {
00316 
00317    LocalVector          trackDirection = trajState.localDirection();
00318    double                  cosine      = trackDirection.z()/trackDirection.mag();
00319    const SiStripCluster*   Cluster     = (sistripsimplehit->cluster()).get();
00320    const vector<uint8_t>&  Ampls       = Cluster->amplitudes();
00321    uint32_t                DetId       = Cluster->geographicalId();
00322    int                     FirstStrip  = Cluster->firstStrip();
00323    bool                    Saturation  = false;
00324    bool                    Overlaping  = false;
00325    int                     Charge      = 0;
00326    stModInfo* MOD                      = MODsColl[DetId];
00327 
00328 
00329    if(!IsFarFromBorder(trajState,DetId, iSetup)){/*printf("tooCloseFromBorder\n");*/return -1;}
00330 
00331 
00332    if(FirstStrip==0                                  )Overlaping=true;
00333    if(FirstStrip==128                                )Overlaping=true;
00334    if(FirstStrip==256                                )Overlaping=true;
00335    if(FirstStrip==384                                )Overlaping=true;
00336    if(FirstStrip==512                                )Overlaping=true;
00337    if(FirstStrip==640                                )Overlaping=true;
00338 
00339    if(FirstStrip<=127 && FirstStrip+Ampls.size()>127)Overlaping=true;
00340    if(FirstStrip<=255 && FirstStrip+Ampls.size()>255)Overlaping=true;
00341    if(FirstStrip<=383 && FirstStrip+Ampls.size()>383)Overlaping=true;
00342    if(FirstStrip<=511 && FirstStrip+Ampls.size()>511)Overlaping=true;
00343    if(FirstStrip<=639 && FirstStrip+Ampls.size()>639)Overlaping=true;
00344 
00345    if(FirstStrip+Ampls.size()==127                   )Overlaping=true;
00346    if(FirstStrip+Ampls.size()==255                   )Overlaping=true;
00347    if(FirstStrip+Ampls.size()==383                   )Overlaping=true;
00348    if(FirstStrip+Ampls.size()==511                   )Overlaping=true;
00349    if(FirstStrip+Ampls.size()==639                   )Overlaping=true;
00350    if(FirstStrip+Ampls.size()==767                   )Overlaping=true;
00351 //   if(!DiscriminatorMode && Overlaping){printf("Overlapping\n");return -1;}
00352 
00353 
00354    for(unsigned int a=0;a<Ampls.size();a++){Charge+=Ampls[a];if(Ampls[a]>=254)Saturation=true;}
00355    double path                    = (10.0*MOD->Thickness)/fabs(cosine);
00356    double ClusterChargeOverPath   = (double)Charge / path ;
00357 
00358    if(Ampls.size()>MaxNrStrips)      {/*printf("tooMuchStrips\n");*/return -1;}
00359 //   if(!DiscriminatorMode && Saturation && !AllowSaturation){printf("Saturation\n");return -1;}
00360 
00361    if(!DiscriminatorMode){
00362       if(MOD->SubDet == StripSubdetector::TIB || MOD->SubDet == StripSubdetector::TOB) Charge_Vs_Path_Barrel->Fill(path,ClusterChargeOverPath);
00363       if(MOD->SubDet == StripSubdetector::TID || MOD->SubDet == StripSubdetector::TEC) Charge_Vs_Path_Endcap->Fill(path,ClusterChargeOverPath);
00364    }else{
00365       TH2F* MapToUse = NULL;
00366       if(MOD->SubDet == StripSubdetector::TIB || MOD->SubDet == StripSubdetector::TOB) MapToUse = PCharge_Vs_Path_Barrel;
00367       if(MOD->SubDet == StripSubdetector::TID || MOD->SubDet == StripSubdetector::TEC) MapToUse = PCharge_Vs_Path_Endcap;
00368 
00369       int   BinX  = MapToUse->GetXaxis()->FindBin(path);
00370       int   BinY  = MapToUse->GetYaxis()->FindBin(ClusterChargeOverPath);
00371 
00372       float Prob = MapToUse->GetBinContent(BinX,BinY);
00373 
00374       if(Prob>=0)MeasurementProbabilities.push_back(Prob);
00375    }
00376 
00377    return ClusterChargeOverPath;
00378 }
00379 
00380 
00381 bool DeDxDiscriminatorProducer::IsFarFromBorder(TrajectoryStateOnSurface trajState, const uint32_t detid, const edm::EventSetup* iSetup)
00382 {
00383   edm::ESHandle<TrackerGeometry> tkGeom; iSetup->get<TrackerDigiGeometryRecord>().get( tkGeom );
00384 
00385   LocalPoint  HitLocalPos   = trajState.localPosition();
00386   LocalError  HitLocalError = trajState.localError().positionError() ;
00387 
00388   const GeomDetUnit* it = tkGeom->idToDetUnit(DetId(detid));
00389   if (dynamic_cast<const StripGeomDetUnit*>(it)==0 && dynamic_cast<const PixelGeomDetUnit*>(it)==0) {
00390      std::cout << "this detID doesn't seem to belong to the Tracker" << std::endl;
00391      return false;
00392   }
00393 
00394   const BoundPlane plane = it->surface();
00395   const TrapezoidalPlaneBounds* trapezoidalBounds( dynamic_cast<const TrapezoidalPlaneBounds*>(&(plane.bounds())));
00396   const RectangularPlaneBounds* rectangularBounds( dynamic_cast<const RectangularPlaneBounds*>(&(plane.bounds())));
00397 
00398   double DistFromBorder = 1.0;
00399   double HalfWidth      = it->surface().bounds().width()  /2.0;
00400   double HalfLength     = it->surface().bounds().length() /2.0;
00401 
00402   if(trapezoidalBounds)
00403   {
00404      std::vector<float> const & parameters = (*trapezoidalBounds).parameters();
00405      HalfLength     = parameters[3];
00406      double t       = (HalfLength + HitLocalPos.y()) / (2*HalfLength) ;
00407      HalfWidth      = parameters[0] + (parameters[1]-parameters[0]) * t;
00408   }else if(rectangularBounds){
00409      HalfWidth      = it->surface().bounds().width()  /2.0;
00410      HalfLength     = it->surface().bounds().length() /2.0;
00411   }else{return false;}
00412 
00413 //  if (fabs(HitLocalPos.x())+HitLocalError.xx() >= (HalfWidth  - DistFromBorder) ) return false;//Don't think is really necessary
00414   if (fabs(HitLocalPos.y())+HitLocalError.yy() >= (HalfLength - DistFromBorder) ) return false;
00415 
00416   return true;
00417 }
00418 
00419 
00420 
00421 //define this as a plug-in
00422 DEFINE_FWK_MODULE(DeDxDiscriminatorProducer);