00001 #include "Validation/RecoTrack/interface/MultiTrackValidator.h"
00002 #include "DQMServices/ClientConfig/interface/FitSlicesYTool.h"
00003
00004 #include "FWCore/Framework/interface/MakerMacros.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006
00007 #include "DataFormats/TrackReco/interface/Track.h"
00008 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00009 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
00010 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
00011 #include "SimDataFormats/Track/interface/SimTrackContainer.h"
00012 #include "SimDataFormats/Vertex/interface/SimVertexContainer.h"
00013 #include "SimTracker/TrackAssociation/interface/TrackAssociatorByChi2.h"
00014 #include "SimTracker/TrackAssociation/interface/TrackAssociatorByHits.h"
00015 #include "SimTracker/TrackerHitAssociation/interface/TrackerHitAssociator.h"
00016 #include "SimTracker/Records/interface/TrackAssociatorRecord.h"
00017 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
00018 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
00019 #include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
00020 #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
00021 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
00022 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
00023 #include "TrackingTools/PatternTools/interface/TSCBLBuilderNoMaterial.h"
00024 #include "SimTracker/TrackAssociation/plugins/ParametersDefinerForTPESProducer.h"
00025 #include "SimTracker/TrackAssociation/plugins/CosmicParametersDefinerForTPESProducer.h"
00026 #include "Validation/RecoTrack/interface/MTVHistoProducerAlgoFactory.h"
00027
00028 #include "DataFormats/TrackReco/interface/DeDxData.h"
00029 #include "DataFormats/Common/interface/ValueMap.h"
00030 #include "DataFormats/Common/interface/Ref.h"
00031
00032 #include "TMath.h"
00033 #include <TF1.h>
00034
00035
00036
00037 using namespace std;
00038 using namespace edm;
00039
00040 typedef edm::Ref<edm::HepMCProduct, HepMC::GenParticle > GenParticleRef;
00041
00042 MultiTrackValidator::MultiTrackValidator(const edm::ParameterSet& pset):MultiTrackValidatorBase(pset){
00043
00044
00045 ParameterSet psetForHistoProducerAlgo = pset.getParameter<ParameterSet>("histoProducerAlgoBlock");
00046 string histoProducerAlgoName = psetForHistoProducerAlgo.getParameter<string>("ComponentName");
00047 histoProducerAlgo_ = MTVHistoProducerAlgoFactory::get()->create(histoProducerAlgoName ,psetForHistoProducerAlgo);
00048 histoProducerAlgo_->setDQMStore(dbe_);
00049
00050 dirName_ = pset.getParameter<std::string>("dirName");
00051 associatormap = pset.getParameter< edm::InputTag >("associatormap");
00052 UseAssociators = pset.getParameter< bool >("UseAssociators");
00053
00054 m_dEdx1Tag = pset.getParameter< edm::InputTag >("dEdx1Tag");
00055 m_dEdx2Tag = pset.getParameter< edm::InputTag >("dEdx2Tag");
00056
00057 tpSelector = TrackingParticleSelector(pset.getParameter<double>("ptMinTP"),
00058 pset.getParameter<double>("minRapidityTP"),
00059 pset.getParameter<double>("maxRapidityTP"),
00060 pset.getParameter<double>("tipTP"),
00061 pset.getParameter<double>("lipTP"),
00062 pset.getParameter<int>("minHitTP"),
00063 pset.getParameter<bool>("signalOnlyTP"),
00064 pset.getParameter<bool>("chargedOnlyTP"),
00065 pset.getParameter<bool>("stableOnlyTP"),
00066 pset.getParameter<std::vector<int> >("pdgIdTP"));
00067
00068 cosmictpSelector = CosmicTrackingParticleSelector(pset.getParameter<double>("ptMinTP"),
00069 pset.getParameter<double>("minRapidityTP"),
00070 pset.getParameter<double>("maxRapidityTP"),
00071 pset.getParameter<double>("tipTP"),
00072 pset.getParameter<double>("lipTP"),
00073 pset.getParameter<int>("minHitTP"),
00074 pset.getParameter<bool>("chargedOnlyTP"),
00075 pset.getParameter<std::vector<int> >("pdgIdTP"));
00076
00077 useGsf = pset.getParameter<bool>("useGsf");
00078 runStandalone = pset.getParameter<bool>("runStandalone");
00079
00080
00081
00082 if (!UseAssociators) {
00083 associators.clear();
00084 associators.push_back(associatormap.label());
00085 }
00086
00087 }
00088
00089
00090 MultiTrackValidator::~MultiTrackValidator(){delete histoProducerAlgo_;}
00091
00092 void MultiTrackValidator::beginRun(Run const&, EventSetup const& setup) {
00093
00094
00095
00096 for (unsigned int ww=0;ww<associators.size();ww++){
00097 for (unsigned int www=0;www<label.size();www++){
00098 dbe_->cd();
00099 InputTag algo = label[www];
00100 string dirName=dirName_;
00101 if (algo.process()!="")
00102 dirName+=algo.process()+"_";
00103 if(algo.label()!="")
00104 dirName+=algo.label()+"_";
00105 if(algo.instance()!="")
00106 dirName+=algo.instance()+"_";
00107 if (dirName.find("Tracks")<dirName.length()){
00108 dirName.replace(dirName.find("Tracks"),6,"");
00109 }
00110 string assoc= associators[ww];
00111 if (assoc.find("Track")<assoc.length()){
00112 assoc.replace(assoc.find("Track"),5,"");
00113 }
00114 dirName+=assoc;
00115 std::replace(dirName.begin(), dirName.end(), ':', '_');
00116
00117 dbe_->setCurrentFolder(dirName.c_str());
00118
00119
00120 histoProducerAlgo_->initialize();
00121
00122 dbe_->goUp();
00123 string subDirName = dirName + "/simulation";
00124 dbe_->setCurrentFolder(subDirName.c_str());
00125
00126
00127 histoProducerAlgo_->bookSimHistos();
00128
00129 dbe_->cd();
00130 dbe_->setCurrentFolder(dirName.c_str());
00131
00132
00133 histoProducerAlgo_->bookRecoHistos();
00134 if (runStandalone) histoProducerAlgo_->bookRecoHistosForStandaloneRunning();
00135
00136 if (UseAssociators) {
00137 edm::ESHandle<TrackAssociatorBase> theAssociator;
00138 for (unsigned int w=0;w<associators.size();w++) {
00139 setup.get<TrackAssociatorRecord>().get(associators[w],theAssociator);
00140 associator.push_back( theAssociator.product() );
00141 }
00142 }
00143 }
00144 }
00145 }
00146
00147 void MultiTrackValidator::analyze(const edm::Event& event, const edm::EventSetup& setup){
00148 using namespace reco;
00149
00150 edm::LogInfo("TrackValidator") << "\n====================================================" << "\n"
00151 << "Analyzing new event" << "\n"
00152 << "====================================================\n" << "\n";
00153 edm::ESHandle<ParametersDefinerForTP> parametersDefinerTP;
00154 setup.get<TrackAssociatorRecord>().get(parametersDefiner,parametersDefinerTP);
00155
00156 edm::Handle<TrackingParticleCollection> TPCollectionHeff ;
00157 event.getByLabel(label_tp_effic,TPCollectionHeff);
00158 const TrackingParticleCollection tPCeff = *(TPCollectionHeff.product());
00159
00160 edm::Handle<TrackingParticleCollection> TPCollectionHfake ;
00161 event.getByLabel(label_tp_fake,TPCollectionHfake);
00162 const TrackingParticleCollection tPCfake = *(TPCollectionHfake.product());
00163
00164
00165
00166
00167
00168
00169 edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
00170 event.getByLabel(bsSrc,recoBeamSpotHandle);
00171 reco::BeamSpot bs = *recoBeamSpotHandle;
00172
00173 edm::Handle< vector<PileupSummaryInfo> > puinfoH;
00174 event.getByLabel(label_pileupinfo,puinfoH);
00175 PileupSummaryInfo puinfo;
00176
00177 for (unsigned int puinfo_ite=0;puinfo_ite<(*puinfoH).size();++puinfo_ite){
00178 if ((*puinfoH)[puinfo_ite].getBunchCrossing()==0){
00179 puinfo=(*puinfoH)[puinfo_ite];
00180 break;
00181 }
00182 }
00183
00184 edm::Handle<TrackingVertexCollection> tvH;
00185 event.getByLabel(label_tv,tvH);
00186 TrackingVertexCollection tv = *tvH;
00187
00188 int w=0;
00189 for (unsigned int ww=0;ww<associators.size();ww++){
00190 for (unsigned int www=0;www<label.size();www++){
00191
00192
00193
00194 edm::Handle<View<Track> > trackCollection;
00195 if(!event.getByLabel(label[www], trackCollection)&&ignoremissingtkcollection_)continue;
00196
00197
00198
00199
00200 reco::RecoToSimCollection recSimColl;
00201 reco::SimToRecoCollection simRecColl;
00202
00203 if(UseAssociators){
00204 edm::LogVerbatim("TrackValidator") << "Analyzing "
00205 << label[www].process()<<":"
00206 << label[www].label()<<":"
00207 << label[www].instance()<<" with "
00208 << associators[ww].c_str() <<"\n";
00209
00210 LogTrace("TrackValidator") << "Calling associateRecoToSim method" << "\n";
00211 recSimColl=associator[ww]->associateRecoToSim(trackCollection,
00212 TPCollectionHfake,
00213 &event);
00214 LogTrace("TrackValidator") << "Calling associateSimToReco method" << "\n";
00215 simRecColl=associator[ww]->associateSimToReco(trackCollection,
00216 TPCollectionHeff,
00217 &event);
00218 }
00219 else{
00220 edm::LogVerbatim("TrackValidator") << "Analyzing "
00221 << label[www].process()<<":"
00222 << label[www].label()<<":"
00223 << label[www].instance()<<" with "
00224 << associatormap.process()<<":"
00225 << associatormap.label()<<":"
00226 << associatormap.instance()<<"\n";
00227
00228 Handle<reco::SimToRecoCollection > simtorecoCollectionH;
00229 event.getByLabel(associatormap,simtorecoCollectionH);
00230 simRecColl= *(simtorecoCollectionH.product());
00231
00232 Handle<reco::RecoToSimCollection > recotosimCollectionH;
00233 event.getByLabel(associatormap,recotosimCollectionH);
00234 recSimColl= *(recotosimCollectionH.product());
00235 }
00236
00237
00238
00239
00240
00241
00242
00243
00244
00245 edm::LogVerbatim("TrackValidator") << "\n# of TrackingParticles: " << tPCeff.size() << "\n";
00246 int ats(0);
00247 int st(0);
00248 unsigned sts(0);
00249 unsigned asts(0);
00250 for (TrackingParticleCollection::size_type i=0; i<tPCeff.size(); i++){
00251 TrackingParticleRef tpr(TPCollectionHeff, i);
00252 TrackingParticle* tp=const_cast<TrackingParticle*>(tpr.get());
00253 ParticleBase::Vector momentumTP;
00254 ParticleBase::Point vertexTP;
00255 double dxySim(0);
00256 double dzSim(0);
00257
00258
00259
00260 if(parametersDefiner=="LhcParametersDefinerForTP")
00261 {
00262 if(! tpSelector(*tp)) continue;
00263 momentumTP = tp->momentum();
00264 vertexTP = tp->vertex();
00265
00266 ParticleBase::Vector momentum = parametersDefinerTP->momentum(event,setup,*tp);
00267 ParticleBase::Point vertex = parametersDefinerTP->vertex(event,setup,*tp);
00268 dxySim = (-vertex.x()*sin(momentum.phi())+vertex.y()*cos(momentum.phi()));
00269 dzSim = vertex.z() - (vertex.x()*momentum.x()+vertex.y()*momentum.y())/sqrt(momentum.perp2())
00270 * momentum.z()/sqrt(momentum.perp2());
00271 }
00272
00273 if(parametersDefiner=="CosmicParametersDefinerForTP")
00274 {
00275 if(! cosmictpSelector(*tp,&bs,event,setup)) continue;
00276 momentumTP = parametersDefinerTP->momentum(event,setup,*tp);
00277 vertexTP = parametersDefinerTP->vertex(event,setup,*tp);
00278 dxySim = (-vertexTP.x()*sin(momentumTP.phi())+vertexTP.y()*cos(momentumTP.phi()));
00279 dzSim = vertexTP.z() - (vertexTP.x()*momentumTP.x()+vertexTP.y()*momentumTP.y())/sqrt(momentumTP.perp2())
00280 * momentumTP.z()/sqrt(momentumTP.perp2());
00281 }
00282
00283
00284 st++;
00285
00286
00287
00288
00289
00290
00291
00292 histoProducerAlgo_->fill_generic_simTrack_histos(w,momentumTP,vertexTP, tp->eventId().bunchCrossing());
00293
00294
00295
00296
00297
00298
00299 const reco::Track* matchedTrackPointer=0;
00300 std::vector<std::pair<RefToBase<Track>, double> > rt;
00301 if(simRecColl.find(tpr) != simRecColl.end()){
00302 rt = (std::vector<std::pair<RefToBase<Track>, double> >) simRecColl[tpr];
00303 if (rt.size()!=0) {
00304 ats++;
00305
00306 matchedTrackPointer = rt.begin()->first.get();
00307 edm::LogVerbatim("TrackValidator") << "TrackingParticle #" << st
00308 << " with pt=" << sqrt(momentumTP.perp2())
00309 << " associated with quality:" << rt.begin()->second <<"\n";
00310 }
00311 }else{
00312 edm::LogVerbatim("TrackValidator")
00313 << "TrackingParticle #" << st
00314 << " with pt,eta,phi: "
00315 << sqrt(momentumTP.perp2()) << " , "
00316 << momentumTP.eta() << " , "
00317 << momentumTP.phi() << " , "
00318 << " NOT associated to any reco::Track" << "\n";
00319 }
00320
00321
00322
00323
00324 std::vector<PSimHit> simhits=tp->trackPSimHit(DetId::Tracker);
00325 int nSimHits = simhits.end()-simhits.begin();
00326
00327 double vtx_z_PU = vertexTP.z();
00328 for (size_t j = 0; j < tv.size(); j++) {
00329 if (tp->eventId().event() == tv[j].eventId().event()) {
00330 vtx_z_PU = tv[j].position().z();
00331 break;
00332 }
00333 }
00334
00335 histoProducerAlgo_->fill_recoAssociated_simTrack_histos(w,*tp,momentumTP,vertexTP,dxySim,dzSim,nSimHits,matchedTrackPointer,puinfo.getPU_NumInteractions(), vtx_z_PU);
00336 sts++;
00337 if (matchedTrackPointer) asts++;
00338
00339
00340
00341
00342 }
00343
00344
00345
00346
00347
00348
00349
00350 edm::LogVerbatim("TrackValidator") << "\n# of reco::Tracks with "
00351 << label[www].process()<<":"
00352 << label[www].label()<<":"
00353 << label[www].instance()
00354 << ": " << trackCollection->size() << "\n";
00355
00356 int sat(0);
00357 int at(0);
00358 int rT(0);
00359
00360
00361
00362
00363
00364 edm::Handle<edm::ValueMap<reco::DeDxData> > dEdx1Handle;
00365 edm::Handle<edm::ValueMap<reco::DeDxData> > dEdx2Handle;
00366 std::vector<edm::ValueMap<reco::DeDxData> > v_dEdx;
00367 v_dEdx.clear();
00368
00369 if (label[www].label()=="generalTracks") {
00370 try {
00371 event.getByLabel(m_dEdx1Tag, dEdx1Handle);
00372 const edm::ValueMap<reco::DeDxData> dEdx1 = *dEdx1Handle.product();
00373 event.getByLabel(m_dEdx2Tag, dEdx2Handle);
00374 const edm::ValueMap<reco::DeDxData> dEdx2 = *dEdx2Handle.product();
00375 v_dEdx.push_back(dEdx1);
00376 v_dEdx.push_back(dEdx2);
00377 } catch (cms::Exception e){
00378 LogTrace("TrackValidator") << "exception found: " << e.what() << "\n";
00379 }
00380 }
00381
00382
00383 for(View<Track>::size_type i=0; i<trackCollection->size(); ++i){
00384
00385 RefToBase<Track> track(trackCollection, i);
00386 rT++;
00387
00388 bool isSigSimMatched(false);
00389 bool isSimMatched(false);
00390 bool isChargeMatched(true);
00391 int numAssocRecoTracks = 0;
00392 int tpbx = 0;
00393 int nSimHits = 0;
00394 double sharedFraction = 0.;
00395 std::vector<std::pair<TrackingParticleRef, double> > tp;
00396 if(recSimColl.find(track) != recSimColl.end()){
00397 tp = recSimColl[track];
00398 if (tp.size()!=0) {
00399 std::vector<PSimHit> simhits=tp[0].first->trackPSimHit(DetId::Tracker);
00400 nSimHits = simhits.end()-simhits.begin();
00401 sharedFraction = tp[0].second;
00402 isSimMatched = true;
00403 if (tp[0].first->charge() != track->charge()) isChargeMatched = false;
00404 if(simRecColl.find(tp[0].first) != simRecColl.end()) numAssocRecoTracks = simRecColl[tp[0].first].size();
00405
00406 tpbx = tp[0].first->eventId().bunchCrossing();
00407 at++;
00408 for (unsigned int tp_ite=0;tp_ite<tp.size();++tp_ite){
00409 TrackingParticle trackpart = *(tp[tp_ite].first);
00410 if ((trackpart.eventId().event() == 0) && (trackpart.eventId().bunchCrossing() == 0)){
00411 isSigSimMatched = true;
00412 sat++;
00413 break;
00414 }
00415 }
00416 edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
00417 << " associated with quality:" << tp.begin()->second <<"\n";
00418 }
00419 } else {
00420 edm::LogVerbatim("TrackValidator") << "reco::Track #" << rT << " with pt=" << track->pt()
00421 << " NOT associated to any TrackingParticle" << "\n";
00422 }
00423
00424
00425 histoProducerAlgo_->fill_generic_recoTrack_histos(w,*track,bs.position(),isSimMatched,isSigSimMatched, isChargeMatched, numAssocRecoTracks, puinfo.getPU_NumInteractions(), tpbx, nSimHits, sharedFraction);
00426
00427
00428
00429 if (v_dEdx.size() > 0) histoProducerAlgo_->fill_dedx_recoTrack_histos(w,track, v_dEdx);
00430
00431
00432
00433
00434
00435
00436 if (tp.size()==0) continue;
00437
00438 histoProducerAlgo_->fill_simAssociated_recoTrack_histos(w,*track);
00439
00440 TrackingParticleRef tpr = tp.begin()->first;
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453
00454
00455
00456
00457
00458 ParticleBase::Vector momentumTP = parametersDefinerTP->momentum(event,setup,*(tpr.get()));
00459 ParticleBase::Point vertexTP = parametersDefinerTP->vertex(event,setup,*(tpr.get()));
00460 int chargeTP = tpr->charge();
00461
00462 histoProducerAlgo_->fill_ResoAndPull_recoTrack_histos(w,momentumTP,vertexTP,chargeTP,
00463 *track,bs.position());
00464
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477 }
00478
00479 histoProducerAlgo_->fill_trackBased_histos(w,at,rT,st);
00480
00481 edm::LogVerbatim("TrackValidator") << "Total Simulated: " << st << "\n"
00482 << "Total Associated (simToReco): " << ats << "\n"
00483 << "Total Reconstructed: " << rT << "\n"
00484 << "Total Associated (recoToSim): " << at << "\n"
00485 << "Total Fakes: " << rT-at << "\n";
00486
00487 w++;
00488 }
00489 }
00490
00491 }
00492
00493 void MultiTrackValidator::endRun(Run const&, EventSetup const&) {
00494 int w=0;
00495 for (unsigned int ww=0;ww<associators.size();ww++){
00496 for (unsigned int www=0;www<label.size();www++){
00497 if(!skipHistoFit && runStandalone) histoProducerAlgo_->finalHistoFits(w);
00498 if (runStandalone) histoProducerAlgo_->fillProfileHistosFromVectors(w);
00499 histoProducerAlgo_->fillHistosFromVectors(w);
00500 w++;
00501 }
00502 }
00503 if ( out.size() != 0 && dbe_ ) dbe_->save(out);
00504 }
00505
00506
00507