CMS 3D CMS Logo

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