00001
00002
00003
00004
00005
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
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
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
00099
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
00118 for(int i=0;i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;i++){
00119
00120 for(int j=0;j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;j++){
00121
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
00133 }else{
00134 Prob_ChargePath->SetBinContent (i, j, k, 0);
00135 }
00136 }
00137 }
00138 }
00139 }else if(strcmp(ProbabilityMode.c_str(),"Integral")==0){
00140
00141 for(int i=0;i<=Prob_ChargePath->GetXaxis()->GetNbins()+1;i++){
00142
00143 for(int j=0;j<=Prob_ChargePath->GetYaxis()->GetNbins()+1;j++){
00144
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
00155 }else{
00156 Prob_ChargePath->SetBinContent (i, j, k, 0);
00157 }
00158 }
00159 }
00160 }
00161 }else{
00162
00163
00164 exit(0);
00165 }
00166
00167
00168
00169
00170
00171
00172
00173
00174
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
00226
00227
00228
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
00299
00300 Error = NClusterSaturating;
00301 dEdxDiscrims[track_index] = DeDxData(estimator, Error, size );
00302
00303
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
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
00335 stModInfo* MOD = MODsColl[detId];
00336
00337
00338
00339 if( ampls.size()>MaxNrStrips) {return -1;}
00340
00341
00342
00343
00344
00345
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
00452 DEFINE_FWK_MODULE(DeDxDiscriminatorProducer);