00001 #include "RecoParticleFlow/PFSimProducer/plugins/PFSimParticleProducer.h"
00002
00003 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
00004 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyResolution.h"
00005
00006 #include "DataFormats/ParticleFlowReco/interface/PFLayer.h"
00007 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
00008
00009 #include "DataFormats/ParticleFlowReco/interface/PFSimParticle.h"
00010 #include "DataFormats/ParticleFlowReco/interface/PFSimParticleFwd.h"
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "FWCore/Framework/interface/ESHandle.h"
00022 #include "MagneticField/Engine/interface/MagneticField.h"
00023 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
00024 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
00025 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
00026 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00027 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00028 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
00029 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
00030 #include "DataFormats/GeometrySurface/interface/TkRotation.h"
00031 #include "DataFormats/GeometrySurface/interface/SimpleCylinderBounds.h"
00032 #include "DataFormats/GeometrySurface/interface/SimpleDiskBounds.h"
00033 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
00034 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00035 #include "TrackingTools/PatternTools/interface/Trajectory.h"
00036 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
00037
00038 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
00039
00040 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00041 #include "FWCore/Utilities/interface/Exception.h"
00042 #include "FWCore/Framework/interface/EventSetup.h"
00043
00044 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
00045 #include "FastSimulation/Event/interface/FSimEvent.h"
00046 #include "FastSimulation/Event/interface/FSimTrack.h"
00047 #include "FastSimulation/Event/interface/FSimVertex.h"
00048 #include "FastSimulation/Particle/interface/ParticleTable.h"
00049
00050 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
00051 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
00052 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00053 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
00054 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00055 #include "DataFormats/ParticleFlowReco/interface/PFRecTrack.h"
00056 #include "DataFormats/TrackReco/interface/Track.h"
00057 #include "DataFormats/TrackReco/interface/TrackFwd.h"
00058 #include "DataFormats/TrackerRecHit2D/interface/SiTrackerGSMatchedRecHit2D.h"
00059 #include "DataFormats/TrackerRecHit2D/interface/SiTrackerGSMatchedRecHit2DCollection.h"
00060
00061 #include <set>
00062 #include <sstream>
00063
00064 using namespace std;
00065 using namespace edm;
00066
00067 PFSimParticleProducer::PFSimParticleProducer(const edm::ParameterSet& iConfig)
00068 {
00069
00070
00071 processParticles_ =
00072 iConfig.getUntrackedParameter<bool>("process_Particles",true);
00073
00074
00075 inputTagSim_
00076 = iConfig.getParameter<InputTag>("sim");
00077
00078
00079
00080
00081 inputTagFamosSimHits_
00082 = iConfig.getUntrackedParameter<InputTag>("famosSimHits");
00083 mctruthMatchingInfo_ =
00084 iConfig.getUntrackedParameter<bool>("MCTruthMatchingInfo",false);
00085
00086
00087 inputTagRecTracks_
00088 = iConfig.getParameter<InputTag>("RecTracks");
00089 inputTagEcalRecHitsEB_
00090 = iConfig.getParameter<InputTag>("ecalRecHitsEB");
00091 inputTagEcalRecHitsEE_
00092 = iConfig.getParameter<InputTag>("ecalRecHitsEE");
00093
00094 verbose_ =
00095 iConfig.getUntrackedParameter<bool>("verbose",false);
00096
00097
00098
00099 produces<reco::PFSimParticleCollection>();
00100
00101 particleFilter_ = iConfig.getParameter<ParameterSet>
00102 ( "ParticleFilter" );
00103
00104 mySimEvent = new FSimEvent( particleFilter_ );
00105 }
00106
00107
00108
00109 PFSimParticleProducer::~PFSimParticleProducer()
00110 {
00111 delete mySimEvent;
00112 }
00113
00114
00115 void
00116 PFSimParticleProducer::beginRun(edm::Run& run,
00117 const edm::EventSetup & es)
00118 {
00119
00120
00121 edm::ESHandle < HepPDT::ParticleDataTable > pdt;
00122
00123 es.getData(pdt);
00124 if ( !ParticleTable::instance() ) ParticleTable::instance(&(*pdt));
00125 mySimEvent->initializePdt(&(*pdt));
00126 }
00127
00128
00129 void PFSimParticleProducer::produce(Event& iEvent,
00130 const EventSetup& iSetup)
00131 {
00132
00133 LogDebug("PFSimParticleProducer")<<"START event: "<<iEvent.id().event()
00134 <<" in run "<<iEvent.id().run()<<endl;
00135
00136
00137
00138
00139
00140 std::vector<unsigned> recTrackSimID;
00141
00142
00143
00144
00145
00146 typedef std::pair<double, unsigned> hitSimID;
00147 typedef std::list< std::pair<double, unsigned> >::iterator ITM;
00148 std::vector< std::list <hitSimID> > caloHitsEBID(62000);
00149 std::vector<double> caloHitsEBTotE(62000,0.0);
00150
00151 if(mctruthMatchingInfo_){
00152
00153
00154 edm::Handle<edm::PCaloHitContainer> pcalohits;
00155
00156
00157
00158
00159 bool found_phit
00160 = iEvent.getByLabel(inputTagFamosSimHits_,pcalohits);
00161
00162
00163 if(!found_phit) {
00164 ostringstream err;
00165 err<<"could not find pcaloHit "<<"famosSimHits:EcalHitsEB";
00166 LogError("PFSimParticleProducer")<<err.str()<<endl;
00167
00168 throw cms::Exception( "MissingProduct", err.str());
00169 }
00170 else {
00171 assert( pcalohits.isValid() );
00172
00173
00174
00175
00176 edm::PCaloHitContainer::const_iterator it
00177 = pcalohits.product()->begin();
00178 edm::PCaloHitContainer::const_iterator itend
00179 = pcalohits.product()->end();
00180
00181
00182 for(;it!=itend;++it)
00183 {
00184 EBDetId detid(it->id());
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194 if(it->energy() > 0.0) {
00195 std::pair<double, unsigned> phitsimid
00196 = make_pair(it->energy(),it->geantTrackId());
00197 caloHitsEBID[detid.hashedIndex()].push_back(phitsimid);
00198 caloHitsEBTotE[detid.hashedIndex()]
00199 += it->energy();
00200 }
00201
00202 }
00203 }
00204
00205
00206
00207 Handle< reco::PFRecTrackCollection > recTracks;
00208 try{
00209 LogDebug("PFSimParticleProducer")<<"getting PFRecTracks"<<endl;
00210 iEvent.getByLabel(inputTagRecTracks_, recTracks);
00211
00212 } catch (cms::Exception& err) {
00213 LogError("PFSimParticleProducer")<<err
00214 <<" cannot get collection "
00215 <<"particleFlowBlock"<<":"
00216 <<""
00217 <<endl;
00218 }
00219
00220
00221
00222 getSimIDs( recTracks, recTrackSimID );
00223
00224 }
00225
00226
00227 if( processParticles_) {
00228
00229 auto_ptr< reco::PFSimParticleCollection >
00230 pOutputPFSimParticleCollection(new reco::PFSimParticleCollection );
00231
00232 Handle<vector<SimTrack> > simTracks;
00233 bool found = iEvent.getByLabel(inputTagSim_,simTracks);
00234 if(!found) {
00235
00236 ostringstream err;
00237 err<<"cannot find sim tracks: "<<inputTagSim_;
00238 LogError("PFSimParticleProducer")<<err.str()<<endl;
00239
00240 throw cms::Exception( "MissingProduct", err.str());
00241 }
00242
00243
00244
00245 Handle<vector<SimVertex> > simVertices;
00246 found = iEvent.getByLabel(inputTagSim_,simVertices);
00247 if(!found) {
00248 LogError("PFSimParticleProducer")
00249 <<"cannot find sim vertices: "<<inputTagSim_<<endl;
00250 return;
00251 }
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261 mySimEvent->fill( *simTracks, *simVertices );
00262
00263 if(verbose_)
00264 mySimEvent->print();
00265
00266
00267 for(unsigned i=0; i<mySimEvent->nTracks(); i++) {
00268
00269 const FSimTrack& fst = mySimEvent->track(i);
00270
00271 int motherId = -1;
00272 if( ! fst.noMother() )
00273 motherId = fst.mother().id();
00274
00275
00276
00277
00278
00279
00280
00281
00282 unsigned recTrackID = 99999;
00283 vector<unsigned> recHitContrib;
00284 vector<double> recHitContribFrac;
00285
00286 if(mctruthMatchingInfo_){
00287
00288 for(unsigned lo=0; lo<recTrackSimID.size();
00289 lo++) {
00290 if( i == recTrackSimID[lo] ) {
00291
00292
00293 recTrackID = lo;
00294 }
00295 }
00296
00297
00298
00299
00300
00301 edm::Handle<EcalRecHitCollection> rhcHandle;
00302 bool found = iEvent.getByLabel(inputTagEcalRecHitsEB_,
00303 rhcHandle);
00304 if(!found) {
00305 ostringstream err;
00306 err<<"could not find rechits "<< inputTagEcalRecHitsEB_;
00307 LogError("PFSimParticleProducer")<<err.str()<<endl;
00308
00309 throw cms::Exception( "MissingProduct", err.str());
00310 }
00311 else {
00312 assert( rhcHandle.isValid() );
00313
00314
00315
00316 EBRecHitCollection::const_iterator it_rh
00317 = rhcHandle.product()->begin();
00318 EBRecHitCollection::const_iterator itend_rh
00319 = rhcHandle.product()->end();
00320
00321 for(;it_rh!=itend_rh;++it_rh)
00322 {
00323 unsigned rhit_hi
00324 = EBDetId(it_rh->id()).hashedIndex();
00325 EBDetId detid(it_rh->id());
00326
00327
00328
00329
00330 ITM it_phit = caloHitsEBID[rhit_hi].begin();
00331 ITM itend_phit = caloHitsEBID[rhit_hi].end();
00332 for(;it_phit!=itend_phit;++it_phit)
00333 {
00334 if(i == it_phit->second)
00335 {
00336
00337
00338 bool alreadyin = false;
00339 for( unsigned ihit = 0; ihit < recHitContrib.size();
00340 ++ihit )
00341 if(detid.rawId() == recHitContrib[ihit])
00342 alreadyin = true;
00343
00344 if(!alreadyin){
00345 double pcalofraction = 0.0;
00346 if(caloHitsEBTotE[rhit_hi] != 0.0)
00347 pcalofraction
00348 = (it_phit->first/caloHitsEBTotE[rhit_hi])*100.0;
00349
00350
00351 recHitContrib.push_back(it_rh->id());
00352 recHitContribFrac.push_back(pcalofraction);
00353 }
00354 }
00355 }
00356
00357 }
00358
00359 }
00360
00361 }
00362
00363
00364
00365
00366
00367
00368
00369
00370 reco::PFSimParticle particle( fst.charge(),
00371 fst.type(),
00372 fst.id(),
00373 motherId,
00374 fst.daughters(),
00375 recTrackID,
00376 recHitContrib,
00377 recHitContribFrac);
00378
00379
00380 const FSimVertex& originVtx = fst.vertex();
00381
00382 math::XYZPoint posOrig( originVtx.position().x(),
00383 originVtx.position().y(),
00384 originVtx.position().z() );
00385
00386 math::XYZTLorentzVector momOrig( fst.momentum().px(),
00387 fst.momentum().py(),
00388 fst.momentum().pz(),
00389 fst.momentum().e() );
00390 reco::PFTrajectoryPoint
00391 pointOrig(-1,
00392 reco::PFTrajectoryPoint::ClosestApproach,
00393 posOrig, momOrig);
00394
00395
00396 particle.addPoint(pointOrig);
00397
00398
00399 if( ! fst.noEndVertex() ) {
00400 const FSimVertex& endVtx = fst.endVertex();
00401
00402 math::XYZPoint posEnd( endVtx.position().x(),
00403 endVtx.position().y(),
00404 endVtx.position().z() );
00405
00406
00407
00408
00409 math::XYZTLorentzVector momEnd;
00410
00411 reco::PFTrajectoryPoint
00412 pointEnd( -1,
00413 reco::PFTrajectoryPoint::BeamPipeOrEndVertex,
00414 posEnd, momEnd);
00415
00416 particle.addPoint(pointEnd);
00417 }
00418 else {
00419 reco::PFTrajectoryPoint dummy;
00420 particle.addPoint(dummy);
00421 }
00422
00423
00424 if( fst.onLayer1() ) {
00425 const RawParticle& rp = fst.layer1Entrance();
00426
00427 math::XYZPoint posLayer1( rp.x(), rp.y(), rp.z() );
00428 math::XYZTLorentzVector momLayer1( rp.px(), rp.py(), rp.pz(), rp.e() );
00429 reco::PFTrajectoryPoint layer1Pt(-1, reco::PFTrajectoryPoint::PS1,
00430 posLayer1, momLayer1);
00431
00432 particle.addPoint( layer1Pt );
00433
00434
00435 }
00436 else {
00437 reco::PFTrajectoryPoint dummy;
00438 particle.addPoint(dummy);
00439 }
00440
00441 if( fst.onLayer2() ) {
00442 const RawParticle& rp = fst.layer2Entrance();
00443
00444 math::XYZPoint posLayer2( rp.x(), rp.y(), rp.z() );
00445 math::XYZTLorentzVector momLayer2( rp.px(), rp.py(), rp.pz(), rp.e() );
00446 reco::PFTrajectoryPoint layer2Pt(-1, reco::PFTrajectoryPoint::PS2,
00447 posLayer2, momLayer2);
00448
00449 particle.addPoint( layer2Pt );
00450
00451
00452 }
00453 else {
00454 reco::PFTrajectoryPoint dummy;
00455 particle.addPoint(dummy);
00456 }
00457
00458 if( fst.onEcal() ) {
00459 const RawParticle& rp = fst.ecalEntrance();
00460
00461 math::XYZPoint posECAL( rp.x(), rp.y(), rp.z() );
00462 math::XYZTLorentzVector momECAL( rp.px(), rp.py(), rp.pz(), rp.e() );
00463 reco::PFTrajectoryPoint ecalPt(-1,
00464 reco::PFTrajectoryPoint::ECALEntrance,
00465 posECAL, momECAL);
00466
00467 particle.addPoint( ecalPt );
00468
00469
00470 }
00471 else {
00472 reco::PFTrajectoryPoint dummy;
00473 particle.addPoint(dummy);
00474 }
00475
00476
00477 reco::PFTrajectoryPoint dummy;
00478 particle.addPoint(dummy);
00479
00480 if( fst.onHcal() ) {
00481
00482 const RawParticle& rpin = fst.hcalEntrance();
00483
00484 math::XYZPoint posHCALin( rpin.x(), rpin.y(), rpin.z() );
00485 math::XYZTLorentzVector momHCALin( rpin.px(), rpin.py(), rpin.pz(),
00486 rpin.e() );
00487 reco::PFTrajectoryPoint hcalPtin(-1,
00488 reco::PFTrajectoryPoint::HCALEntrance,
00489 posHCALin, momHCALin);
00490
00491 particle.addPoint( hcalPtin );
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504 }
00505 else {
00506 reco::PFTrajectoryPoint dummy;
00507 particle.addPoint(dummy);
00508 }
00509
00510 pOutputPFSimParticleCollection->push_back( particle );
00511 }
00512
00513 iEvent.put(pOutputPFSimParticleCollection);
00514 }
00515
00516
00517 LogDebug("PFSimParticleProducer")<<"STOP event: "<<iEvent.id().event()
00518 <<" in run "<<iEvent.id().run()<<endl;
00519 }
00520
00521 void PFSimParticleProducer::getSimIDs( const TrackHandle& trackh,
00522 std::vector<unsigned>& recTrackSimID )
00523 {
00524
00525 if( trackh.isValid() ) {
00526
00527
00528 for(unsigned i=0;i<trackh->size(); i++) {
00529
00530 reco::PFRecTrackRef ref( trackh,i );
00531 const reco::PFRecTrack& PFT = *ref;
00532 const reco::TrackRef trackref = PFT.trackRef();
00533
00534
00535
00536
00537
00538
00539 trackingRecHit_iterator rhitbeg
00540 = trackref->recHitsBegin();
00541 trackingRecHit_iterator rhitend
00542 = trackref->recHitsEnd();
00543 for (trackingRecHit_iterator it = rhitbeg;
00544 it != rhitend; it++){
00545
00546 if( it->get()->isValid() ){
00547
00548 const SiTrackerGSMatchedRecHit2D * rechit
00549 = (const SiTrackerGSMatchedRecHit2D*) (it->get());
00550
00551
00552
00553
00554
00555
00556 recTrackSimID.push_back( rechit->simtrackId() );
00557 break;
00558 }
00559 }
00560 }
00561 }
00562 }