Reference Manual

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 //
00003 // Package:    MyPrimaryVertexAnalyzer4PU
00004 // Class:      MyPrimaryVertexAnalyzer4PU
00005 // 
00013 //
00014 // Original Author:  Wolfram Erdmann
00015 
00016 
00017 // system include files
00018 #include <memory>
00019 #include <string>
00020 #include <vector>
00021  
00022 
00023 // user include files
00024 #include "FWCore/Framework/interface/Frameworkfwd.h"
00025 #include "FWCore/Framework/interface/EDAnalyzer.h"
00026 #include "FWCore/Framework/interface/Event.h"
00027 #include "FWCore/Framework/interface/MakerMacros.h"
00028 
00029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00030 #include "DataFormats/TrackReco/interface/Track.h"
00031 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00032 #include "FWCore/Utilities/interface/InputTag.h"
00033 
00034 //generator level
00035 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
00036 #include "HepMC/GenEvent.h"
00037 #include "HepMC/GenVertex.h"
00038 #include "HepMC/GenParticle.h"
00039 
00040 // vertex stuff
00042 #include <DataFormats/VertexReco/interface/VertexFwd.h>
00043 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
00044 
00045 // simulated vertices,..., add <use name=SimDataFormats/Vertex> and <../Track>
00046 #include <SimDataFormats/Vertex/interface/SimVertex.h>
00047 #include <SimDataFormats/Vertex/interface/SimVertexContainer.h>
00048 #include <SimDataFormats/Track/interface/SimTrack.h>
00049 #include <SimDataFormats/Track/interface/SimTrackContainer.h>
00050 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
00051 //#include "DataFormats/Math/interface/LorentzVector.h"
00052 #include <SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h>
00053 #include <SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h>
00054 #include <SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h>
00055 #include "SimTracker/TrackAssociation/interface/TrackAssociatorBase.h"
00056 #include "SimDataFormats/EncodedEventId/interface/EncodedEventId.h"
00057 
00058 //Track et al
00059 #include "DataFormats/TrackReco/interface/Track.h"
00060 #include "DataFormats/VertexReco/interface/Vertex.h"
00061 #include "DataFormats/Math/interface/Point3D.h"
00062 #include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
00063 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
00064 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
00065 
00066 // Root
00067 #include <TH1.h>
00068 #include <TFile.h>
00069 
00070 #include "RecoVertex/PrimaryVertexProducer/interface/TrackFilterForPVFinding.h"
00071 
00072 
00073 
00074 
00075 // class declaration
00076 class PrimaryVertexAnalyzer4PU : public edm::EDAnalyzer {
00077 
00078  typedef math::XYZTLorentzVector LorentzVector;
00079 
00080 typedef reco::TrackBase::ParameterVector ParameterVector;
00081 struct SimPart{
00082   ParameterVector par;
00083   int type; // 0 = primary
00084   double zdcap;   // z@dca' (closest approach to the beam
00085   double ddcap;
00086   double zvtx;    // z of the vertex
00087   double xvtx;    // x of the vertex
00088   double yvtx;    // y of the vertex
00089   int pdg;        // particle pdg id
00090   int rec;
00091 };
00092 
00093 // auxiliary class holding simulated primary vertices
00094 class simPrimaryVertex {
00095 public:
00096   
00097   simPrimaryVertex(double x1,double y1,double z1):x(x1),y(y1),z(z1),ptsq(0),nGenTrk(0){
00098     ptot.setPx(0);
00099     ptot.setPy(0);
00100     ptot.setPz(0);
00101     ptot.setE(0);
00102     cluster=-1;
00103     nclutrk=0;
00104     p4=LorentzVector(0,0,0,0) ;
00105     //    event=0;
00106   };
00107   double x,y,z;
00108   HepMC::FourVector ptot;
00109   LorentzVector p4;
00110   double ptsq;
00111   int nGenTrk;
00112   int nMatchedTracks;
00113   int cluster;
00114   //int event;
00115   EncodedEventId eventId;
00116   double nclutrk;
00117   std::vector<int> finalstateParticles;
00118   std::vector<int> simTrackIndex;
00119   std::vector<int> matchedRecTrackIndex;
00120   std::vector<int> genVertex;
00121   std::vector<reco::Track> reconstructedTracks;
00122   const reco::Vertex *recVtx;
00123 };
00124 
00125 
00126 // auxiliary class holding simulated events
00127 class SimEvent {
00128 public:
00129   
00130   SimEvent(){
00131     //event=-1;
00132     nrecTrack=0;
00133     z=-99;
00134     zfit=-99;
00135     sumpt2rec=0.;
00136     sumpt2=0;
00137     sumpt=0;
00138     Tc=-1;
00139     dzmax=0;
00140     dztrim=0;
00141     chisq=0;
00142   };
00143   double x,y,z;
00144   double xfit,yfit,zfit;
00145   int nrecTrack;
00146   //int event;
00147   EncodedEventId eventId;
00148   std::vector<const TrackingParticle*> tp;
00149   std::vector<reco::TransientTrack> tk;
00150   std::vector<reco::TransientTrack> tkprim;
00151   std::vector<reco::TransientTrack> tkprimsel;
00152   double sumpt2rec;
00153   double sumpt2,sumpt;
00154   double Tc,chisq,dzmax,dztrim,m4m2;
00155   // rec vertex matching
00156   double zmatch;
00157   int nmatch;
00158   std::map<double, int> ntInRecVz;  // number of tracks in recvtx at z
00159 };
00160 
00161 
00162 
00163 public:
00164   explicit PrimaryVertexAnalyzer4PU(const edm::ParameterSet&);
00165   ~PrimaryVertexAnalyzer4PU();
00166   
00167   virtual void analyze(const edm::Event&, const edm::EventSetup&);
00168   //virtual void beginJob(edm::EventSetup const&);
00169   virtual void beginJob();
00170   virtual void endJob();
00171 
00172 private:
00173   void printPVTrks(const edm::Handle<reco::TrackCollection> &recTrks, 
00174                    const edm::Handle<reco::VertexCollection> &recVtxs,  
00175                    std::vector<SimPart>& tsim,
00176                    std::vector<SimEvent>& simEvt,
00177                    const bool selectedOnly=true);
00178 
00179   int* supf(std::vector<SimPart>& simtrks, const reco::TrackCollection & trks);
00180   static bool match(const ParameterVector  &a, const ParameterVector &b);
00181   std::vector<SimPart> getSimTrkParameters( edm::Handle<edm::SimTrackContainer> & simTrks,
00182                                             edm::Handle<edm::SimVertexContainer> & simVtcs,
00183                                             double simUnit=1.0);
00184   void getTc(const std::vector<reco::TransientTrack>&,double &, double &, double &, double &, double&);
00185   void add(std::map<std::string, TH1*>& h, TH1* hist){  h[hist->GetName()]=hist; hist->StatOverflows(kTRUE);}
00186 
00187   void Fill(std::map<std::string, TH1*>& h, std::string s, double x){
00188     //    cout << "Fill1 " << s << endl;
00189     if(h.count(s)==0){
00190       std::cout << "Trying to fill non-exiting Histogram named " << s << std::endl;
00191       return;
00192     }
00193     h[s]->Fill(x);
00194   }
00195 
00196   void Fill(std::map<std::string, TH1*>& h, std::string s, double x, double y){
00197     //    cout << "Fill2 " << s << endl;
00198     if(h.count(s)==0){
00199       std::cout << "Trying to fill non-exiting Histogram named " << s << std::endl;
00200       return;
00201     }
00202     h[s]->Fill(x,y);
00203   }
00204 
00205   void Fill(std::map<std::string, TH1*>& h, std::string s, double x, bool signal){
00206     if(h.count(s)==0){
00207       std::cout << "Trying to fill non-exiting Histogram named " << s << std::endl;
00208       return;
00209     }
00210 
00211     h[s]->Fill(x);
00212     if(signal){
00213       if(h.count(s+"Signal")==0){
00214         std::cout << "Trying to fill non-exiting Histogram named " << s+"Signal" << std::endl;
00215         return;
00216       }
00217       h[s+"Signal"]->Fill(x);
00218     }else{
00219       if(h.count(s+"PU")==0){
00220         std::cout << "Trying to fill non-exiting Histogram named " << s+"PU" << std::endl;
00221         return;
00222       }
00223       h[s+"PU"]->Fill(x);
00224     }
00225   }
00226 
00227   void Fill(std::map<std::string, TH1*>& h, std::string s, bool yesno, bool signal){
00228     if (yesno){
00229       Fill(h, s,1.,signal);
00230     }else{
00231       Fill(h, s,0.,signal);
00232     }
00233   }
00234 
00235   void Cumulate(TH1* h){
00236     
00237     if((h->GetEntries()==0) || (h->Integral()<=0) ){
00238       std::cout << "DEBUG : Cumulate called with empty histogram " << h->GetTitle() << std::endl;
00239       return;
00240     }
00241     std::cout << "DEBUG : cumulating  " << h->GetTitle() << std::endl;
00242     try{
00243       h->ComputeIntegral();
00244       Double_t * integral=h->GetIntegral();
00245       h->SetContent(integral);
00246     }catch(...){
00247       std::cout << "DEBUG : an error occurred cumulating  " << h->GetTitle()  <<  std::endl;
00248     }
00249     std::cout << "DEBUG : cumulating  " << h->GetTitle() << "done " <<  std::endl;
00250   }
00251 
00252   std::map<std::string, TH1*> bookVertexHistograms();
00253 
00254   bool matchVertex(const simPrimaryVertex  &vsim, 
00255                    const reco::Vertex       &vrec);
00256   bool isResonance(const HepMC::GenParticle * p);
00257   bool isFinalstateParticle(const HepMC::GenParticle * p);
00258   bool isCharged(const HepMC::GenParticle * p);
00259   void fillTrackHistos(std::map<std::string, TH1*> & h, const std::string & ttype, const reco::Track & t, const reco::Vertex *v = NULL);
00260   void dumpHitInfo(const reco::Track & t);
00261   void printRecTrks( const edm::Handle<reco::TrackCollection> & recTrks);
00262   void printRecVtxs(const edm::Handle<reco::VertexCollection> recVtxs,  std::string title="Reconstructed Vertices");
00263   void printSimVtxs(const edm::Handle<edm::SimVertexContainer> simVtxs);
00264   void printSimTrks(const edm::Handle<edm::SimTrackContainer> simVtrks);
00265   std::vector<simPrimaryVertex> getSimPVs(const edm::Handle<edm::HepMCProduct> evtMC);
00266   std::vector<simPrimaryVertex> getSimPVs(const edm::Handle<edm::HepMCProduct> evt, 
00267                                           const edm::Handle<edm::SimVertexContainer> simVtxs, 
00268                                           const edm::Handle<edm::SimTrackContainer> simTrks);
00269   std::vector<PrimaryVertexAnalyzer4PU::simPrimaryVertex> getSimPVs(const edm::Handle<TrackingVertexCollection>);
00270 
00271   bool truthMatchedTrack( edm::RefToBase<reco::Track>, TrackingParticleRef &  );
00272   std::vector< edm::RefToBase<reco::Track> >  getTruthMatchedVertexTracks(
00273                                        const reco::Vertex&
00274                                        );
00275 
00276   std::vector<PrimaryVertexAnalyzer4PU::SimEvent> getSimEvents(
00277                                                               edm::Handle<TrackingParticleCollection>, 
00278                                                               edm::Handle<TrackingVertexCollection>,
00279                                                               edm::Handle<edm::View<reco::Track> >
00280                                                               );
00281 
00282   void matchRecTracksToVertex(simPrimaryVertex & pv, 
00283                               const std::vector<SimPart > & tsim,
00284                               const edm::Handle<reco::TrackCollection> & recTrks);
00285 
00286   void analyzeVertexCollection(std::map<std::string, TH1*> & h,
00287                                const edm::Handle<reco::VertexCollection> recVtxs,
00288                                const edm::Handle<reco::TrackCollection> recTrks, 
00289                                std::vector<simPrimaryVertex> & simpv,
00290                                const std::string message="");
00291 
00292   void analyzeVertexCollectionTP(std::map<std::string, TH1*> & h,
00293                                const edm::Handle<reco::VertexCollection> recVtxs,
00294                                const edm::Handle<reco::TrackCollection> recTrks, 
00295                                std::vector<SimEvent> & simEvt,
00296                                  const std::string message="");
00297 
00298   void printEventSummary(std::map<std::string, TH1*> & h,
00299                          const edm::Handle<reco::VertexCollection> recVtxs,
00300                                const edm::Handle<reco::TrackCollection> recTrks, 
00301                          std::vector<SimEvent> & simEvt,
00302                          const std::string message);
00303 
00304   void history(const edm::Handle<edm::View<reco::Track> > & tracks,const size_t trackindex=10000);
00305   std::string particleString(int) const;
00306   std::string vertexString(
00307     TrackingParticleRefVector,
00308     TrackingParticleRefVector
00309   ) const;
00310   std::string vertexString(
00311     HepMC::GenVertex::particles_in_const_iterator,
00312     HepMC::GenVertex::particles_in_const_iterator,
00313     HepMC::GenVertex::particles_out_const_iterator,
00314     HepMC::GenVertex::particles_out_const_iterator
00315   ) const;
00316 
00317 
00318   // ----------member data ---------------------------
00319   std::string recoTrackProducer_;
00320   std::string outputFile_;       // output file
00321   std::vector<std::string> vtxSample_;        // make this a a vector to keep cfg compatibility with PrimaryVertexAnalyzer
00322   double fBfield_;
00323   TFile*  rootFile_;             
00324   bool verbose_;
00325   bool doMatching_;
00326   bool printXBS_;
00327   edm::InputTag simG4_;
00328   double simUnit_;     
00329   double zmatch_;
00330   edm::ESHandle < ParticleDataTable > pdt_;
00331   math::XYZPoint myBeamSpot;
00332   // local counters
00333   int eventcounter_;
00334   int dumpcounter_;
00335   int ndump_;
00336   bool dumpThisEvent_;
00337   bool dumpPUcandidates_;
00338 
00339   // from the event setup
00340   int run_;
00341   int luminosityBlock_;
00342   int event_;
00343   int bunchCrossing_;
00344   int orbitNumber_;
00345 
00346   bool DEBUG_;
00347   
00348 
00349   std::map<std::string, TH1*> hBS;
00350   std::map<std::string, TH1*> hnoBS;
00351   std::map<std::string, TH1*> hDA;
00352   std::map<std::string, TH1*> hPIX;
00353   std::map<std::string, TH1*> hMVF;
00354   std::map<std::string, TH1*> hsimPV;
00355 
00356   TrackAssociatorBase * associatorByHits_;
00357   reco::RecoToSimCollection r2s_;
00358   std::map<double, TrackingParticleRef> z2tp_;
00359 
00360   TrackFilterForPVFinding theTrackFilter;
00361   reco::BeamSpot vertexBeamSpot_;
00362   double wxy2_;
00363   edm::Handle<reco::BeamSpot> recoBeamSpotHandle_;
00364   edm::ESHandle<TransientTrackBuilder> theB_;
00365   edm::InputTag beamSpot_;
00366 };
00367
/data/refman/pasoursint/CMSSW_6_1_1/src/Validation/RecoVertex/interface/PrimaryVertexAnalyzer4PU.h
