CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_4_1_8_patch12/src/RecoParticleFlow/PFSimProducer/plugins/PFSimParticleProducer.cc

Go to the documentation of this file.
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 // include files used for reconstructed tracks
00013 // #include "DataFormats/TrackReco/interface/Track.h"
00014 // #include "DataFormats/TrackReco/interface/TrackFwd.h"
00015 // #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
00016 // #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
00017 // #include "DataFormats/TrackCandidate/interface/TrackCandidate.h"
00018 // #include "DataFormats/TrackCandidate/interface/TrackCandidateCollection.h"
00019 // #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
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 // #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
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   //retrieving collections for MC Truth Matching
00079 
00080   //modif-beg
00081   inputTagFamosSimHits_ 
00082     = iConfig.getUntrackedParameter<InputTag>("famosSimHits");
00083   mctruthMatchingInfo_ = 
00084     iConfig.getUntrackedParameter<bool>("MCTruthMatchingInfo",false);
00085   //modif-end
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   // register products
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   // init Particle data table (from Pythia)
00121   edm::ESHandle < HepPDT::ParticleDataTable > pdt;
00122   //  edm::ESHandle < DefaultConfig::ParticleDataTable > pdt;
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   //MC Truth Matching only with Famos and UnFoldedMode option to true!!
00137 
00138   //vector to store the trackIDs of rectracks corresponding
00139   //to the simulated particle.
00140   std::vector<unsigned> recTrackSimID;
00141   
00142   //In order to know which simparticule contribute to 
00143   //a given Ecal RecHit energy, we need to access 
00144   //the PCAloHit from FastSim. 
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     //getting the PCAloHit
00154     edm::Handle<edm::PCaloHitContainer> pcalohits;
00155     //   bool found_phit 
00156     //     = iEvent.getByLabel("famosSimHits","EcalHitsEB",
00157     //                  pcalohits);  
00158     //modif-beg
00159     bool found_phit 
00160       = iEvent.getByLabel(inputTagFamosSimHits_,pcalohits);
00161     //modif-end
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       //     cout << "PFSimParticleProducer: number of pcalohits="
00174       //         << pcalohits->size() << endl;
00175       
00176       edm::PCaloHitContainer::const_iterator it    
00177         = pcalohits.product()->begin();
00178       edm::PCaloHitContainer::const_iterator itend 
00179         = pcalohits.product()->end();
00180       
00181       //loop on the PCaloHit from FastSim Calorimetry
00182       for(;it!=itend;++it)
00183         {
00184           EBDetId detid(it->id());
00185           
00186           //    cout << detid << " " << detid.rawId()
00187           //         << " " <<  detid.hashedIndex() 
00188           //         << " " << it->energy() 
00189           //         << " " << it->id() 
00190           //         << " trackId=" 
00191           //         << it->geantTrackId() 
00192           //         << endl;
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(); //summing pcalhit energy   
00200           }//energy > 0
00201           
00202         }//loop PcaloHits    
00203     }//pcalohit handle access
00204     
00205     //Retrieving the PFRecTrack collection for
00206     //Monte Carlo Truth Matching tool
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     }//pfrectrack handle access
00219     
00220     //getting the simID corresponding to 
00221     //each of the PFRecTracks
00222     getSimIDs( recTracks, recTrackSimID );
00223     
00224   }//mctruthMatchingInfo_ //modif
00225   
00226   // deal with true particles 
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     //     for(unsigned it = 0; it<simTracks->size(); it++ ) {
00256     //       cout<<"\t track "<< (*simTracks)[it]<<" "
00257     //    <<(*simTracks)[it].momentum().vect().perp()<<" "
00258     //    <<(*simTracks)[it].momentum().e()<<endl;
00259     //     }
00260 
00261     mySimEvent->fill( *simTracks, *simVertices );
00262       
00263     if(verbose_) 
00264       mySimEvent->print();
00265 
00266     // const std::vector<FSimTrack>& fsimTracks = *(mySimEvent->tracks() );
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() ) // a mother exist
00273         motherId = fst.mother().id();
00274 
00275       //This is finding out the simID corresponding 
00276       //to the recTrack
00277 //       cout << "Particle " << i 
00278 //         << " " << fst.genpartIndex() 
00279 //         << " -------------------------------------" << endl;
00280 
00281       //GETTING THE TRACK ID
00282       unsigned         recTrackID = 99999;
00283       vector<unsigned> recHitContrib; //modif
00284       vector<double>   recHitContribFrac; //modif
00285 
00286       if(mctruthMatchingInfo_){ //modif
00287 
00288         for(unsigned lo=0; lo<recTrackSimID.size(); 
00289             lo++) {
00290           if( i == recTrackSimID[lo] ) {
00291 //          cout << "Corresponding Rec Track " 
00292 //               << lo << endl;
00293             recTrackID = lo;
00294           }//match track
00295         }//loop rectrack
00296 //      if( recTrackID == 99999 ) 
00297 //        cout << "Sim Track not reconstructed pT=" <<  
00298 //          fst.momentum().pt() << endl;
00299         
00300         // get the ecalBarrel rechits for MC truth matching tool
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 //        cout << "PFSimParticleProducer: number of rechits="
00314 //             << rhcHandle->size() << endl;
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 //          cout << detid << " " << detid.rawId()
00327 //               << " " <<  detid.hashedIndex() 
00328 //               << " " << it_rh->energy() << endl;
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                       //Alex (08/10/08) TO BE REMOVED, eliminating
00337                       //duplicated rechits
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                         //store info
00351                         recHitContrib.push_back(it_rh->id());
00352                         recHitContribFrac.push_back(pcalofraction);
00353                       }//selected rechits           
00354                     }//matching
00355                 }//loop pcalohit
00356               
00357             }//loop rechits
00358           
00359         }//getting the rechits
00360 
00361       }//mctruthMatchingInfo_ //modif
00362 
00363 //       cout << "This particule has " << recHitContrib.size() 
00364 //         << " rechit contribution" << endl;
00365 //       for( unsigned ih = 0; ih < recHitContrib.size(); ++ih )
00366 //      cout << recHitContrib[ih] 
00367 //           << " f=" << recHitContribFrac[ih] << " ";
00368 //       cout << endl;
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       // point 0 is origin vertex
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         //       cout<<"end vertex : "
00406         //        <<endVtx.position().x()<<" "
00407         //        <<endVtx.position().y()<<endl;
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 { // add a dummy point
00419         reco::PFTrajectoryPoint dummy;
00420         particle.addPoint(dummy); 
00421       }
00422 
00423 
00424       if( fst.onLayer1() ) { // PS layer1
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         // extrapolate to cluster depth
00435       }
00436       else { // add a dummy point
00437         reco::PFTrajectoryPoint dummy;
00438         particle.addPoint(dummy); 
00439       }
00440 
00441       if( fst.onLayer2() ) { // PS layer2
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         // extrapolate to cluster depth
00452       }
00453       else { // add a dummy point
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         // extrapolate to cluster depth
00470       }
00471       else { // add a dummy point
00472         reco::PFTrajectoryPoint dummy;
00473         particle.addPoint(dummy); 
00474       }
00475         
00476       // add a dummy point for ECAL Shower max
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         //      const RawParticle& rpout = fst.hcalExit();
00495         
00496         //      math::XYZPoint posHCALout( rpout.x(), rpout.y(), rpout.z() );
00497         //      math::XYZTLorentzVector momHCALout( rpout.px(), rpout.py(), rpout.pz(),
00498         //                                          rpout.e() );
00499         //      reco::PFTrajectoryPoint 
00500         //        hcalPtout(0, reco::PFTrajectoryPoint::HCALEntrance, 
00501         //                  posHCAL, momHCAL);
00502         
00503         //      particle.addPoint( hcalPtout );         
00504       }
00505       else { // add a dummy point
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 //     cout << "Size=" << trackh->size() << endl;
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 //       double Pt  = trackref->pt(); 
00535 //       double DPt = trackref->ptError();
00536 //       cout << " PFBlockProducer: PFrecTrack->Track Pt= " 
00537 //         << Pt << " DPt = " << DPt << endl;
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 //        cout <<  "rechit"            
00552 //             << " corresponding simId " 
00553 //             << rechit->simtrackId() 
00554 //             << endl;
00555 
00556           recTrackSimID.push_back( rechit->simtrackId() );
00557           break;
00558         }
00559       }//loop track rechit
00560     }//loop recTracks
00561   }//track handle valid
00562 }