CMS 3D CMS Logo

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

Go to the documentation of this file.
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.2 2012/01/17 13:46:51 innocent Exp $
00019 //
00020 //
00021 
00022 
00023 // system include files
00024 #include <memory>
00025 #include "DataFormats/Common/interface/ValueMap.h"
00026 #include "DataFormats/TrackReco/interface/DeDxData.h"
00027 
00028 #include "RecoTracker/DeDx/plugins/DeDxDiscriminatorProducer.h"
00029 #include "DataFormats/TrackReco/interface/Track.h"
00030 #include "TrackingTools/PatternTools/interface/TrajTrackAssociation.h"
00031 
00032 #include "FWCore/Framework/interface/ESHandle.h"
00033 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00034 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00035 
00036 #include "CondFormats/DataRecord/interface/SiStripDeDxMip_3D_Rcd.h"
00037 #include "CondFormats/DataRecord/interface/SiStripDeDxElectron_3D_Rcd.h"
00038 #include "CondFormats/DataRecord/interface/SiStripDeDxProton_3D_Rcd.h"
00039 #include "CondFormats/DataRecord/interface/SiStripDeDxPion_3D_Rcd.h"
00040 #include "CondFormats/DataRecord/interface/SiStripDeDxKaon_3D_Rcd.h"
00041 
00042 
00043 #include "TFile.h"
00044 
00045 using namespace reco;
00046 using namespace std;
00047 using namespace edm;
00048 
00049 DeDxDiscriminatorProducer::DeDxDiscriminatorProducer(const edm::ParameterSet& iConfig)
00050 {
00051 
00052    produces<ValueMap<DeDxData> >();
00053 
00054    m_tracksTag = iConfig.getParameter<edm::InputTag>("tracks");
00055    m_trajTrackAssociationTag   = iConfig.getParameter<edm::InputTag>("trajectoryTrackAssociation");
00056 
00057    usePixel = iConfig.getParameter<bool>("UsePixel"); 
00058    useStrip = iConfig.getParameter<bool>("UseStrip");
00059    if(!usePixel && !useStrip)
00060    edm::LogWarning("DeDxHitsProducer") << "Pixel Hits AND Strip Hits will not be used to estimate dEdx --> BUG, Please Update the config file";
00061 
00062    Formula             = iConfig.getUntrackedParameter<unsigned>("Formula"            ,  0);
00063    Reccord             = iConfig.getUntrackedParameter<string>  ("Reccord"            , "SiStripDeDxMip_3D_Rcd");
00064    ProbabilityMode     = iConfig.getUntrackedParameter<string>  ("ProbabilityMode"    , "Accumulation");
00065 
00066 
00067    MinTrackMomentum    = iConfig.getUntrackedParameter<double>  ("minTrackMomentum"   ,  0.0);
00068    MaxTrackMomentum    = iConfig.getUntrackedParameter<double>  ("maxTrackMomentum"   ,  99999.0); 
00069    MinTrackEta         = iConfig.getUntrackedParameter<double>  ("minTrackEta"        , -5.0);
00070    MaxTrackEta         = iConfig.getUntrackedParameter<double>  ("maxTrackEta"        ,  5.0);
00071    MaxNrStrips         = iConfig.getUntrackedParameter<unsigned>("maxNrStrips"        ,  255);
00072    MinTrackHits        = iConfig.getUntrackedParameter<unsigned>("MinTrackHits"       ,  3);
00073 
00074    shapetest           = iConfig.getParameter<bool>("ShapeTest");
00075    useCalibration      = iConfig.getParameter<bool>("UseCalibration");
00076    m_calibrationPath   = iConfig.getParameter<string>("calibrationPath");
00077 
00078    Prob_ChargePath = NULL;
00079 }
00080 
00081 
00082 DeDxDiscriminatorProducer::~DeDxDiscriminatorProducer(){}
00083 
00084 // ------------ method called once each job just before starting event loop  ------------
00085 void  DeDxDiscriminatorProducer::beginRun(edm::Run & run, const edm::EventSetup& iSetup)
00086 {
00087    edm::ESHandle<PhysicsTools::Calibration::HistogramD3D> DeDxMapHandle_;    
00088    if(      strcmp(Reccord.c_str(),"SiStripDeDxMip_3D_Rcd")==0){
00089       iSetup.get<SiStripDeDxMip_3D_Rcd>() .get(DeDxMapHandle_);
00090    }else if(strcmp(Reccord.c_str(),"SiStripDeDxPion_3D_Rcd")==0){
00091       iSetup.get<SiStripDeDxPion_3D_Rcd>().get(DeDxMapHandle_);
00092    }else if(strcmp(Reccord.c_str(),"SiStripDeDxKaon_3D_Rcd")==0){
00093       iSetup.get<SiStripDeDxKaon_3D_Rcd>().get(DeDxMapHandle_);
00094    }else if(strcmp(Reccord.c_str(),"SiStripDeDxProton_3D_Rcd")==0){
00095       iSetup.get<SiStripDeDxProton_3D_Rcd>().get(DeDxMapHandle_);
00096    }else if(strcmp(Reccord.c_str(),"SiStripDeDxElectron_3D_Rcd")==0){
00097       iSetup.get<SiStripDeDxElectron_3D_Rcd>().get(DeDxMapHandle_);
00098    }else{
00099 //      printf("The reccord %s is not known by the DeDxDiscriminatorProducer\n", Reccord.c_str());
00100 //      printf("Program will exit now\n");
00101       exit(0);
00102    }
00103    DeDxMap_ = *DeDxMapHandle_.product();
00104 
00105    double xmin = DeDxMap_.rangeX().min;
00106    double xmax = DeDxMap_.rangeX().max;
00107    double ymin = DeDxMap_.rangeY().min;
00108    double ymax = DeDxMap_.rangeY().max;
00109    double zmin = DeDxMap_.rangeZ().min;
00110    double zmax = DeDxMap_.rangeZ().max;
00111 
00112    if(Prob_ChargePath)delete Prob_ChargePath;
00113    Prob_ChargePath  = new TH3D ("Prob_ChargePath"     , "Prob_ChargePath" , DeDxMap_.numberOfBinsX(), xmin, xmax, DeDxMap_.numberOfBinsY() , ymin, ymax, DeDxMap_.numberOfBinsZ(), zmin, zmax);
00114 
00115    
00116 
00117    if(strcmp(ProbabilityMode.c_str(),"Accumulation")==0){
00118 //      printf("LOOOP ON P\n");
00119       for(int i=0;i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;i++){
00120 //         printf("LOOOP ON PATH\n");
00121          for(int j=0;j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;j++){
00122 //            printf("LOOOP ON CHARGE\n");
00123 
00124             double Ni = 0;
00125             for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){ Ni+=DeDxMap_.binContent(i,j,k);} 
00126 
00127             for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){
00128                double tmp = 0;
00129                for(int l=0;l<=k;l++){ tmp+=DeDxMap_.binContent(i,j,l);}
00130 
00131                if(Ni>0){
00132                   Prob_ChargePath->SetBinContent (i, j, k, tmp/Ni);
00133 //                printf("P=%6.2f Path=%6.2f Charge%8.2f --> Prob=%8.3f\n",Prob_ChargePath->GetXaxis()->GetBinCenter(i), Prob_ChargePath->GetYaxis()->GetBinCenter(j), Prob_ChargePath->GetZaxis()->GetBinCenter(k), tmp/Ni);
00134                }else{
00135                   Prob_ChargePath->SetBinContent (i, j, k, 0);
00136                }
00137             }
00138          }
00139       }
00140    }else if(strcmp(ProbabilityMode.c_str(),"Integral")==0){
00141 //      printf("LOOOP ON P\n");
00142       for(int i=0;i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;i++){
00143 //         printf("LOOOP ON PATH\n");
00144          for(int j=0;j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;j++){
00145 //            printf("LOOOP ON CHARGE\n");
00146 
00147             double Ni = 0;
00148             for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){ Ni+=DeDxMap_.binContent(i,j,k);}
00149 
00150             for(int k=0;k<=Prob_ChargePath->GetZaxis()->GetNbins()+1;k++){
00151                double tmp = DeDxMap_.binContent(i,j,k);
00152 
00153                if(Ni>0){
00154                   Prob_ChargePath->SetBinContent (i, j, k, tmp/Ni);
00155 //                  printf("P=%6.2f Path=%6.2f Charge%8.2f --> Prob=%8.3f\n",Prob_ChargePath->GetXaxis()->GetBinCenter(i), Prob_ChargePath->GetYaxis()->GetBinCenter(j), Prob_ChargePath->GetZaxis()->GetBinCenter(k), tmp/Ni);
00156                }else{
00157                   Prob_ChargePath->SetBinContent (i, j, k, 0);
00158                }
00159             }
00160          }
00161       }   
00162    }else{
00163 //      printf("The ProbabilityMode: %s is unknown\n",ProbabilityMode.c_str());
00164 //      printf("The program will stop now\n");
00165       exit(0);
00166    }
00167 
00168 
00169 
00170 /*
00171    for(int i=0;i<Prob_ChargePath->GetXaxis()->GetNbins();i++){
00172       for(int j=0;j<Prob_ChargePath->GetYaxis()->GetNbins();j++){
00173          double tmp = DeDxMap_.binContent(i,j);
00174          Prob_ChargePath->SetBinContent (i, j, tmp);
00175          printf("%i %i --> %f\n",i,j,tmp);
00176       }
00177    }
00178 */
00179 
00180 
00181 
00182    if(MODsColl.size()==0){
00183       edm::ESHandle<TrackerGeometry> tkGeom;
00184       iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
00185       m_tracker = tkGeom.product();
00186 
00187       vector<GeomDet*> Det = tkGeom->dets();
00188       for(unsigned int i=0;i<Det.size();i++){
00189          DetId  Detid  = Det[i]->geographicalId();
00190          int    SubDet = Detid.subdetId();
00191 
00192          if( SubDet == StripSubdetector::TIB ||  SubDet == StripSubdetector::TID ||
00193              SubDet == StripSubdetector::TOB ||  SubDet == StripSubdetector::TEC  ){
00194 
00195              StripGeomDetUnit* DetUnit     = dynamic_cast<StripGeomDetUnit*> (Det[i]);
00196              if(!DetUnit)continue;
00197 
00198              const StripTopology& Topo     = DetUnit->specificTopology();
00199              unsigned int         NAPV     = Topo.nstrips()/128;
00200 
00201              double Eta     = DetUnit->position().basicVector().eta();
00202              double R       = DetUnit->position().basicVector().transverse();
00203              double Thick   = DetUnit->surface().bounds().thickness();
00204 
00205              stModInfo* MOD = new stModInfo;
00206              MOD->DetId     = Detid.rawId();
00207              MOD->SubDet    = SubDet;
00208              MOD->Eta       = Eta;
00209              MOD->R         = R;
00210              MOD->Thickness = Thick;
00211              MOD->NAPV      = NAPV;
00212              MODsColl[MOD->DetId] = MOD;
00213          }
00214       }
00215  
00216       MakeCalibrationMap();
00217    }
00218 
00219 }
00220 
00221 void  DeDxDiscriminatorProducer::endJob()
00222 {
00223    MODsColl.clear();
00224 }
00225 
00226 
00227 
00228 void DeDxDiscriminatorProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup)
00229 {
00230   auto_ptr<ValueMap<DeDxData> > trackDeDxDiscrimAssociation(new ValueMap<DeDxData> );  
00231   ValueMap<DeDxData>::Filler filler(*trackDeDxDiscrimAssociation);
00232 
00233   Handle<TrajTrackAssociationCollection> trajTrackAssociationHandle;
00234   iEvent.getByLabel(m_trajTrackAssociationTag, trajTrackAssociationHandle);
00235   const TrajTrackAssociationCollection TrajToTrackMap = *trajTrackAssociationHandle.product();
00236 
00237   edm::Handle<reco::TrackCollection> trackCollectionHandle;
00238   iEvent.getByLabel(m_tracksTag,trackCollectionHandle);
00239 
00240    edm::ESHandle<TrackerGeometry> tkGeom;
00241    iSetup.get<TrackerDigiGeometryRecord>().get( tkGeom );
00242    m_tracker = tkGeom.product();
00243  
00244    std::vector<DeDxData> dEdxDiscrims( TrajToTrackMap.size() );
00245 
00246    unsigned track_index = 0;
00247    for(TrajTrackAssociationCollection::const_iterator it = TrajToTrackMap.begin(); it!=TrajToTrackMap.end(); ++it, track_index++) {
00248       dEdxDiscrims[track_index] = DeDxData(-1, -2, 0 );
00249 
00250       const Track      track = *it->val;
00251       const Trajectory traj  = *it->key;
00252 
00253       if(track.eta()  <MinTrackEta      || track.eta()>MaxTrackEta     ){continue;}
00254       if(track.p()    <MinTrackMomentum || track.p()  >MaxTrackMomentum){continue;}
00255       if(track.found()<MinTrackHits                                    ){continue;}
00256 
00257       std::vector<double> vect_probs;
00258       vector<TrajectoryMeasurement> measurements = traj.measurements();
00259 
00260       unsigned int NClusterSaturating = 0;
00261       for(vector<TrajectoryMeasurement>::const_iterator measurement_it = measurements.begin(); measurement_it!=measurements.end(); measurement_it++){
00262 
00263          TrajectoryStateOnSurface trajState = measurement_it->updatedState();
00264          if( !trajState.isValid() ) continue;
00265 
00266          const TrackingRecHit*         hit                 = (*measurement_it->recHit()).hit();
00267          const SiStripRecHit2D*        sistripsimplehit    = dynamic_cast<const SiStripRecHit2D*>(hit);
00268          const SiStripMatchedRecHit2D* sistripmatchedhit   = dynamic_cast<const SiStripMatchedRecHit2D*>(hit);
00269          const SiStripRecHit1D*        sistripsimple1dhit  = dynamic_cast<const SiStripRecHit1D*>(hit);
00270          
00271          double Prob;
00272          if(sistripsimplehit){
00273                 Prob = GetProbability(DeDxTools::GetCluster(sistripsimplehit), trajState,sistripsimplehit->geographicalId());       
00274                 if(shapetest && !(DeDxTools::shapeSelection(DeDxTools::GetCluster(sistripsimplehit)->amplitudes()))) Prob=-1.0;
00275                 if(Prob>=0) vect_probs.push_back(Prob);             
00276             
00277                 if(ClusterSaturatingStrip(DeDxTools::GetCluster(sistripsimplehit),sistripsimplehit->geographicalId())>0)NClusterSaturating++;
00278 
00279          }else if(sistripmatchedhit){
00280                 Prob = GetProbability(DeDxTools::GetCluster(sistripmatchedhit->monoHit()), trajState, sistripmatchedhit->monoId());
00281                 if(shapetest && !(DeDxTools::shapeSelection(DeDxTools::GetCluster(sistripmatchedhit->monoHit())->amplitudes()))) Prob=-1.0;
00282                 if(Prob>=0) vect_probs.push_back(Prob);
00283            
00284                 Prob = GetProbability(DeDxTools::GetCluster(sistripmatchedhit->stereoHit()), trajState,sistripmatchedhit->stereoId());
00285                 if(Prob>=0) vect_probs.push_back(Prob);
00286             
00287                 if(ClusterSaturatingStrip(DeDxTools::GetCluster(sistripmatchedhit->monoHit()),sistripmatchedhit->monoId())  >0)NClusterSaturating++;
00288                 if(ClusterSaturatingStrip(DeDxTools::GetCluster(sistripmatchedhit->stereoHit()),sistripmatchedhit->stereoId())>0)NClusterSaturating++;
00289          }else if(sistripsimple1dhit){ 
00290                 Prob = GetProbability(DeDxTools::GetCluster(sistripsimple1dhit), trajState, sistripsimple1dhit->geographicalId());
00291                 if(shapetest && !(DeDxTools::shapeSelection(DeDxTools::GetCluster(sistripsimple1dhit)->amplitudes()))) Prob=-1.0;
00292                 if(Prob>=0) vect_probs.push_back(Prob);
00293                 if(ClusterSaturatingStrip(DeDxTools::GetCluster(sistripsimple1dhit),sistripsimple1dhit->geographicalId())>0)NClusterSaturating++;
00294          }else{
00295          }
00296       }
00297 
00298       double estimator          = ComputeDiscriminator(vect_probs);
00299       int    size               = vect_probs.size();
00300       float  Error              = -1;
00301 
00302       //WARNING: Since the dEdX Error is not properly computed for the moment
00303       //It was decided to store the number of saturating cluster in that dataformat
00304       Error = NClusterSaturating;
00305       dEdxDiscrims[track_index] = DeDxData(estimator, Error, size );
00306 
00307 //      printf("%i --> %g\n",size,estimator);
00308    }
00309 
00310   filler.insert(trackCollectionHandle, dEdxDiscrims.begin(), dEdxDiscrims.end());
00311   filler.fill();
00312   iEvent.put(trackDeDxDiscrimAssociation);
00313 }
00314 
00315 
00316 int DeDxDiscriminatorProducer::ClusterSaturatingStrip(const SiStripCluster*   cluster,
00317                                                       const uint32_t &               detId){
00318    const vector<uint8_t>&  ampls          = cluster->amplitudes();
00319 
00320    stModInfo* MOD                         = NULL;
00321    if(useCalibration)MOD                  = MODsColl[detId];
00322 
00323    int SaturatingStrip = 0;
00324    for(unsigned int i=0;i<ampls.size();i++){
00325       int StripCharge = ampls[i];
00326       if(MOD){StripCharge = (int)(StripCharge / MOD->Gain);}
00327       if(StripCharge>=254)SaturatingStrip++;
00328    }
00329    return SaturatingStrip;
00330 }
00331 
00332 double DeDxDiscriminatorProducer::GetProbability(const SiStripCluster*   cluster, TrajectoryStateOnSurface trajState,
00333                                                  const uint32_t &               detId)
00334 {
00335    // Get All needed variables
00336    LocalVector             trackDirection = trajState.localDirection();
00337    double                  cosine         = trackDirection.z()/trackDirection.mag();
00338    const vector<uint8_t>&  ampls          = cluster->amplitudes();
00339    //   uint32_t                detId          = cluster->geographicalId();
00340    //   int                     firstStrip     = cluster->firstStrip();
00341    stModInfo* MOD                         = MODsColl[detId];
00342 
00343 
00344    // Sanity Checks
00345    if( ampls.size()>MaxNrStrips)                                                                      {return -1;}
00346 // if( DeDxDiscriminatorTools::IsSaturatingStrip  (ampls))                                            {return -1;}
00347 // if( DeDxDiscriminatorTools::IsSpanningOver2APV (firstStrip, ampls.size()))                         {return -1;}
00348 // if(!DeDxDiscriminatorTools::IsFarFromBorder    (trajState, m_tracker->idToDetUnit(DetId(detId)) )) {return -1;}
00349 
00350 
00351    // Find Probability for this given Charge and Path
00352    double charge = 0;
00353    if(useCalibration){
00354       for(unsigned int i=0;i<ampls.size();i++){
00355          int CalibratedCharge = ampls[i];
00356          CalibratedCharge = (int)(CalibratedCharge / MOD->Gain);
00357          if(CalibratedCharge>=1024){
00358             CalibratedCharge = 255;
00359          }else if(CalibratedCharge>254){
00360             CalibratedCharge = 254;
00361          }
00362          charge+=CalibratedCharge;
00363       }
00364    }else{
00365       charge = DeDxDiscriminatorTools::charge(ampls);
00366    }
00367    double path   = DeDxDiscriminatorTools::path(cosine,MOD->Thickness);
00368 
00369    int    BinX   = Prob_ChargePath->GetXaxis()->FindBin(trajState.localMomentum().mag());
00370    int    BinY   = Prob_ChargePath->GetYaxis()->FindBin(path);
00371    int    BinZ   = Prob_ChargePath->GetZaxis()->FindBin(charge/path);
00372    return Prob_ChargePath->GetBinContent(BinX,BinY,BinZ);
00373 }
00374 
00375 
00376 double DeDxDiscriminatorProducer::ComputeDiscriminator(std::vector<double>& vect_probs)
00377 {
00378    double estimator = -1;
00379    int    size      = vect_probs.size();
00380    if(size<=0)return estimator;
00381 
00382    if(Formula==0){
00383       double P = 1;
00384       for(int i=0;i<size;i++){
00385          if(vect_probs[i]<=0.0001){P *= pow(0.0001       , 1.0/size);}
00386          else                     {P *= pow(vect_probs[i], 1.0/size);}
00387       }
00388       estimator = P;
00389    }else if(Formula==1){
00390       std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
00391       for(int i=0;i<size;i++){if(vect_probs[i]<=0.0001)vect_probs[i] = 0.0001;    }
00392 
00393       double SumJet = 0.;
00394       for(int i=0;i<size;i++){ SumJet+= log(vect_probs[i]); }
00395 
00396       double Loginvlog=log(-SumJet);
00397       double Prob =1.;
00398       double lfact=1.;
00399 
00400       for(int l=1; l!=size; l++){
00401          lfact*=l;
00402          Prob+=exp(l*Loginvlog-log(1.*lfact));
00403       }
00404 
00405       double LogProb=log(Prob);
00406       double ProbJet=std::min(exp(std::max(LogProb+SumJet,-30.)),1.);
00407       estimator = -log10(ProbJet)/4.;
00408       estimator = 1-estimator;
00409    }else if(Formula==2){
00410       std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
00411       double P = 1.0/(12*size);
00412       for(int i=1;i<=size;i++){
00413          P += pow(vect_probs[i-1] - ((2.0*i-1.0)/(2.0*size)),2);
00414       }
00415       P *= (3.0/size);
00416       estimator = P;
00417    }else{
00418       std::sort(vect_probs.begin(), vect_probs.end(), std::less<double>() );
00419       double P = 1.0/(12*size);
00420       for(int i=1;i<=size;i++){
00421          P += vect_probs[i-1] * pow(vect_probs[i-1] - ((2.0*i-1.0)/(2.0*size)),2);
00422       }
00423       P *= (3.0/size);
00424       estimator = P;
00425    }
00426 
00427    return estimator;
00428 }
00429 
00430 
00431 void DeDxDiscriminatorProducer::MakeCalibrationMap(){
00432    if(!useCalibration)return;
00433 
00434 
00435    TChain* t1 = new TChain("SiStripCalib/APVGain");
00436    t1->Add(m_calibrationPath.c_str());
00437 
00438    unsigned int  tree_DetId;
00439    unsigned char tree_APVId;
00440    double        tree_Gain;
00441 
00442    t1->SetBranchAddress("DetId"             ,&tree_DetId      );
00443    t1->SetBranchAddress("APVId"             ,&tree_APVId      );
00444    t1->SetBranchAddress("Gain"              ,&tree_Gain       );
00445 
00446    for (unsigned int ientry = 0; ientry < t1->GetEntries(); ientry++) {
00447        t1->GetEntry(ientry);
00448        stModInfo* MOD  = MODsColl[tree_DetId];
00449        MOD->Gain = tree_Gain;
00450    }
00451 
00452    delete t1;
00453 
00454 }
00455 
00456 //define this as a plug-in
00457 DEFINE_FWK_MODULE(DeDxDiscriminatorProducer);