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