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