Merge pull request cms-sw#9765 from inugent/7_5_X_ValidationDoubleCou…

…ntingFix_v2

adding David Lange's comments built on rebase of ValidationDoubleCountin...

Validation/EventGenerator/interface/BPhysicsSpectrum.h

@@ -0,0 +1,55 @@
+#ifndef BPhysicsSpectrum_H
+#define BPhysicsSpectrum_H
+
+/*class BPhysicsSpectrum
+ *  
+ *  Class to fill Event Generator dqm monitor elements; works on HepMCProduct
+ *
+ *
+ */
+#include <iostream>
+#include "TMath.h"
+// framework & common header files
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/Run.h"
+
+#include "DataFormats/Common/interface/Handle.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+
+//DQM services
+#include "DQMServices/Core/interface/DQMStore.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "DQMServices/Core/interface/MonitorElement.h"
+#include "DQMServices/Core/interface/DQMEDAnalyzer.h"
+
+#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
+
+#include "Validation/EventGenerator/interface/DQMHelper.h"
+
+
+
+class BPhysicsSpectrum : public DQMEDAnalyzer {
+ public:
+  explicit BPhysicsSpectrum(const edm::ParameterSet&);
+  virtual ~BPhysicsSpectrum();
+
+  virtual void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override;
+  virtual void dqmBeginRun(const edm::Run& r, const edm::EventSetup& c) override;
+  virtual void analyze(edm::Event const&, edm::EventSetup const&) override;
+
+ private:
+  MonitorElement *mass, *Nobj; 
+  edm::InputTag genparticleCollection_;
+  edm::EDGetTokenT<reco::GenParticleCollection> genparticleCollectionToken_;
+  std::string name;
+  double mass_min, mass_max;
+  std::vector<int> Particles;
+};
+
+#endif

Validation/EventGenerator/interface/BPhysicsValidation.h

@@ -0,0 +1,91 @@
+#ifndef BPhysicsValidation_H
+#define BPhysicsValidation_H
+
+/*class BPhysicsValidation
+ *  
+ *  Class to fill Event Generator dqm monitor elements; works on HepMCProduct
+ *
+ *
+ */
+#include <iostream>
+#include "TMath.h"
+// framework & common header files
+
+#include "FWCore/Framework/interface/Event.h"
+#include "FWCore/Framework/interface/EventSetup.h"
+#include "FWCore/Framework/interface/Run.h"
+
+#include "DataFormats/Common/interface/Handle.h"
+#include "FWCore/Framework/interface/ESHandle.h"
+#include "FWCore/ParameterSet/interface/ParameterSet.h"
+#include "FWCore/Utilities/interface/InputTag.h"
+
+//DQM services
+#include "DQMServices/Core/interface/DQMStore.h"
+#include "FWCore/ServiceRegistry/interface/Service.h"
+#include "DQMServices/Core/interface/MonitorElement.h"
+#include "DQMServices/Core/interface/DQMEDAnalyzer.h"
+
+#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticle.h"
+#include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
+
+#include "Validation/EventGenerator/interface/DQMHelper.h"
+
+
+
+class BPhysicsValidation : public DQMEDAnalyzer {
+ public:
+  explicit BPhysicsValidation(const edm::ParameterSet&);
+  virtual ~BPhysicsValidation();
+
+  virtual void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override;
+  virtual void dqmBeginRun(const edm::Run& r, const edm::EventSetup& c) override;
+  virtual void analyze(edm::Event const&, edm::EventSetup const&) override;
+
+ private:
+
+  class ParticleMonitor{
+  public:
+  ParticleMonitor(std::string name_,const edm::ParameterSet &p_):p(p_),name(name_),pdgid(p.getParameter<int>("pdgid")){};
+    ~ParticleMonitor(){};
+
+    void Configure(DQMStore::IBooker &i){
+      std::string pname=p.getParameter<std::string>("pname");
+      double mass_min=p.getParameter<double>("massmin");
+      double mass_max=p.getParameter<double>("massmax");
+      DQMHelper dqm(&i); i.setCurrentFolder("Generator/BPhysics");
+      // Number of analyzed events
+      pt   = dqm.book1dHisto(name+"PT", "P_{t} of the "+pname+"s", 100, 0., 100,"P_{t} (GeV)","Number of Events");
+      eta  = dqm.book1dHisto(name+"ETA", "#eta of the "+pname+"s", 100, -5., 5.,"#eta","Number of Events");
+      phi  = dqm.book1dHisto(name+"PHI", "#phi of the "+pname+"s", 100, 0, 2*TMath::Pi(),"#phi","Number of Events");
+      mass = dqm.book1dHisto(name+"MASS", "Mass of the "+pname+"s", 100, mass_min, mass_max,"Mass (GeV)","Number of Events");
+    }
+
+    void Fill(const reco::GenParticle* p, double weight){
+      if(abs(p->pdgId())==abs(pdgid)){
+	pt->Fill(p->pt(),weight);
+	eta->Fill(p->eta(),weight);
+	phi->Fill(p->phi(),weight);
+	mass->Fill(p->mass(),weight);
+      }
+    }
+    int PDGID(){return pdgid;}
+
+  private:
+    const edm::ParameterSet p;
+    std::string name;
+    int pdgid;
+    MonitorElement *pt, *eta, *phi, *mass; 
+  };
+
+  void FillDaughters(const reco::GenParticle* p);
+  edm::InputTag genparticleCollection_;
+  edm::EDGetTokenT<reco::GenParticleCollection> genparticleCollectionToken_;
+  std::string name;
+  ParticleMonitor particle;
+  std::vector<ParticleMonitor> daughters;
+  MonitorElement *Nobj;
+};
+
+#endif

Validation/EventGenerator/interface/BasicHepMCValidation.h

@@ -23,12 +23,14 @@
 #include "FWCore/ServiceRegistry/interface/Service.h"
 #include "DQMServices/Core/interface/MonitorElement.h"
 #include "DQMServices/Core/interface/DQMEDAnalyzer.h"
+#include "Validation/EventGenerator/interface/DQMHelper.h"
 
 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
 
 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
 
 #include "Validation/EventGenerator/interface/WeightManager.h"
+#include "TVector3.h"
 
 class BasicHepMCValidation : public DQMEDAnalyzer{
     public:
@@ -43,77 +45,148 @@ class BasicHepMCValidation : public DQMEDAnalyzer{
 	WeightManager wmanager_;
     edm::InputTag hepmcCollection_;
 
-	/// PDT table
-	edm::ESHandle<HepPDT::ParticleDataTable> fPDGTable ;
+    /// PDT table
+    edm::ESHandle<HepPDT::ParticleDataTable> fPDGTable ;
+
+
+    class ParticleMonitor{
+    public:
+    ParticleMonitor(std::string name_,int pdgid_, DQMStore::IBooker &i,bool nlog_=false):name(name_),pdgid(pdgid_),count(0),nlog(nlog_){
+	DQMHelper dqm(&i);
+	// Number of analyzed events
+	if(!nlog){
+	numberPerEvent= dqm.book1dHisto(name+"Number", "Number of  "+name+"'s per event", 
+					20, 0, 20,"No. of "+name,"Number of Events");
+	}
+	else{
+	  numberPerEvent= dqm.book1dHisto(name+"Number", "Number of  "+name+"'s per event",
+					  20, 0, 20,"log_{10}(No. of "+name+")","Number of Events");
+	}
+	p_init  = dqm.book1dHisto(name+"Momentum", "log_{10}(P) of the "+name+"s", 
+				  60, -2, 4,"log_{10}(P) (log_{10}(GeV))","Number of "+name );
 
+	eta_init  = dqm.book1dHisto(name+"Eta", "#eta of the "+name+"s", 
+				   100, -5., 5.,"#eta","Number of "+name);
 
-    MonitorElement* nEvt;
-
-	///multiplicity ME's
-	MonitorElement *uNumber, *dNumber, *sNumber, *cNumber, *bNumber, *tNumber;
-	MonitorElement *ubarNumber, *dbarNumber, *sbarNumber, *cbarNumber, *bbarNumber, *tbarNumber;
-	//
-	MonitorElement *eminusNumber, *nueNumber, *muminusNumber, *numuNumber, *tauminusNumber, *nutauNumber;
-	MonitorElement *eplusNumber, *nuebarNumber, *muplusNumber, *numubarNumber, *tauplusNumber, *nutaubarNumber;
-	//
-	MonitorElement *gluNumber, *WplusNumber,*WminusNumber, *ZNumber, *gammaNumber;
-	//
-	MonitorElement *piplusNumber, *piminusNumber, *pizeroNumber, *KplusNumber, *KminusNumber, *KlzeroNumber, *KszeroNumber;
-	MonitorElement *pNumber, *pbarNumber, *nNumber, *nbarNumber, *l0Number, *l0barNumber;
-	//
-	MonitorElement *DplusNumber, *DminusNumber, *DzeroNumber, *BplusNumber, *BminusNumber, *BzeroNumber, *BszeroNumber;
-	//
-	MonitorElement *otherPtclNumber;
+	lifetime_init  = dqm.book1dHisto(name+"LifeTime", "#phi of the "+name+"s", 
+					 100, -15, -5,"Log_{10}(life-time^{final}) (log_{10}(s))","Number of "+name);
 
-	///Momentum ME's
-	MonitorElement *uMomentum, *dMomentum, *sMomentum, *cMomentum, *bMomentum, *tMomentum;
-	MonitorElement *ubarMomentum, *dbarMomentum, *sbarMomentum, *cbarMomentum, *bbarMomentum, *tbarMomentum;
-	//
-	MonitorElement *eminusMomentum, *nueMomentum, *muminusMomentum, *numuMomentum, *tauminusMomentum, *nutauMomentum;
-	MonitorElement *eplusMomentum, *nuebarMomentum, *muplusMomentum, *numubarMomentum, *tauplusMomentum, *nutaubarMomentum;
-	//
-	MonitorElement *gluMomentum, *WplusMomentum,*WminusMomentum, *ZMomentum, *gammaMomentum;
-	//
-	MonitorElement *piplusMomentum, *piminusMomentum, *pizeroMomentum, *KplusMomentum, *KminusMomentum, *KlzeroMomentum,  *KszeroMomentum;
-	//
-	MonitorElement *pMomentum, *pbarMomentum, *nMomentum, *nbarMomentum, *l0Momentum, *l0barMomentum;
-	//
-	MonitorElement *DplusMomentum, *DminusMomentum, *DzeroMomentum,  *BplusMomentum, *BminusMomentum, *BzeroMomentum, *BszeroMomentum;
-	//
-	MonitorElement *otherPtclMomentum;
-
-	///other ME's
-	MonitorElement *genPtclNumber; 
-	MonitorElement *genVrtxNumber;
-	MonitorElement *unknownPDTNumber;
-	MonitorElement *outVrtxPtclNumber;
+	p_final = dqm.book1dHisto(name+"MomentumFinal", "log_{10}(P^{final}) of "+name+"s at end of decay chain", 
+				  60, -2, 4,"log_{10}(P^{final}) (log_{10}(GeV))","Number of "+name);
+
+	lifetime_final=dqm.book1dHisto(name+"LifeTimeFinal", "Log_{10}(life-time^{final}) of "+name+"s at end of decay chain",
+					100,-15,-5,"Log_{10}(life-time^{final}) (log_{10}(s))","Number of "+name);
+      }
+
+      ~ParticleMonitor(){};
+
+      bool Fill(const HepMC::GenParticle* p, double weight){
+	if(p->pdg_id()==pdgid){
+	  if(isFirst(p)){
+	    p_init->Fill(log10(p->momentum().rho()),weight);
+	    eta_init->Fill(p->momentum().eta(),weight);
+            const HepMC::GenParticle* pf=GetFinal(p); // inlcude mixing
+            p_final->Fill(log10(pf->momentum().rho()),weight);
+	    // compute lifetime...
+	    if(p->production_vertex() && p->end_vertex()){
+	      TVector3 PV(p->production_vertex()->point3d().x(),p->production_vertex()->point3d().y(),p->production_vertex()->point3d().z()); 
+	      TVector3 SV(p->end_vertex()->point3d().x(),p->end_vertex()->point3d().y(),p->end_vertex()->point3d().z()); 
+	      TVector3 DL=SV-PV; 
+	      double c(2.99792458E8),Ltau(DL.Mag()/100)/*cm->m*/,beta(p->momentum().rho()/p->momentum().m()); 
+	      double lt=Ltau/(c*beta);
+	      if(lt>1E-16)lifetime_init->Fill(log10(lt),weight);
+	      if(pf->end_vertex()){
+		TVector3 SVf(pf->end_vertex()->point3d().x(),pf->end_vertex()->point3d().y(),pf->end_vertex()->point3d().z());
+		DL=SVf-PV;
+		Ltau=DL.Mag()/100;
+		lt=Ltau/(c*beta);
+		if(lt>1E-16)lifetime_final->Fill(log10(lt),weight);
+	      }
+	    }
+	    count++;
+	  }
+	  return true;
+	}
+	return false;
+      }
+
+      void FillCount(double weight){
+	if(nlog) numberPerEvent->Fill(log10(count),weight);
+	else numberPerEvent->Fill(count,weight);
+	count=0;
+      }
+
+      int PDGID(){return pdgid;}
+
+    private:
+      bool isFirst(const HepMC::GenParticle* p){
+	if(p->production_vertex()){
+	  for(HepMC::GenVertex::particles_in_const_iterator m=p->production_vertex()->particles_in_const_begin(); m!=p->production_vertex()->particles_in_const_end();m++){
+	    if(abs((*m)->pdg_id())==abs(p->pdg_id())) return false;
+	  }
+	}
+	return true;
+      }
+
+      const HepMC::GenParticle* GetFinal(const HepMC::GenParticle* p){ // includes mixing
+	if(p->end_vertex()){ 
+	  if(p->end_vertex()->particles_out_size()!=0){ 
+	    for(HepMC::GenVertex::particles_out_const_iterator d=p->end_vertex()->particles_out_const_begin(); d!=p->end_vertex()->particles_out_const_end();d++){ 
+	      if(abs((*d)->pdg_id())==abs(p->pdg_id())){ 
+		return GetFinal(*d); 
+	      } 
+	    } 
+	  } 
+	} 
+	return p; 
+      } 
+
+      std::string name;
+      int pdgid;
+      unsigned int count;
+      bool nlog;
+      MonitorElement *p_init, *p_final, *eta_init, *lifetime_init, *lifetime_final, *numberPerEvent;
+    };
+
+
+    MonitorElement* nEvt;
+    std::vector<ParticleMonitor> particles;
+
+    ///other ME's
+    MonitorElement *otherPtclNumber;
+    MonitorElement *otherPtclMomentum;
+    MonitorElement *genPtclNumber; 
+    MonitorElement *genVrtxNumber;
+    MonitorElement *unknownPDTNumber;
+    MonitorElement *outVrtxPtclNumber;
     MonitorElement *genPtclStatus;
     //
-	MonitorElement *stablePtclNumber;
-	MonitorElement *stableChaNumber;
-	MonitorElement *stablePtclPhi;
-	MonitorElement *stablePtclEta;
-	MonitorElement *stablePtclCharge;
-	MonitorElement *stablePtclp;
-	MonitorElement *stablePtclpT;
-	MonitorElement *partonNumber;
-	MonitorElement *partonpT;
-	MonitorElement *outVrtxStablePtclNumber;
-	//
-	MonitorElement *vrtxZ;
-	MonitorElement *vrtxRadius;
-	//
-	MonitorElement *Bjorken_x;
-
+    MonitorElement *stablePtclNumber;
+    MonitorElement *stableChaNumber;
+    MonitorElement *stablePtclPhi;
+    MonitorElement *stablePtclEta;
+    MonitorElement *stablePtclCharge;
+    MonitorElement *stablePtclp;
+    MonitorElement *stablePtclpT;
+    MonitorElement *partonNumber;
+    MonitorElement *partonpT;
+    MonitorElement *outVrtxStablePtclNumber;
+    //
+    MonitorElement *vrtxZ;
+    MonitorElement *vrtxRadius;
+    //
+    MonitorElement *Bjorken_x;
+    
     MonitorElement *status1ShortLived;
-
+
+    MonitorElement *log10DeltaEcms;
     MonitorElement *DeltaEcms;
     MonitorElement *DeltaPx;
     MonitorElement *DeltaPy;
     MonitorElement *DeltaPz;
-
+    
     edm::EDGetTokenT<edm::HepMCProduct> hepmcCollectionToken_;
-
+    
 };
 
 #endif

Validation/EventGenerator/interface/DQMHelper.h

@@ -17,6 +17,9 @@ class DQMHelper{
  public: 
   DQMHelper(DQMStore::IBooker *i);
   virtual ~DQMHelper();
+
+  MonitorElement* book1dHisto(std::string name,std::string title,int n,double xmin,double xmax,std::string xaxis, std::string yaxis);
+  MonitorElement* book2dHisto(std::string name,std::string title,int nx,double xmin,double xmax,int ny,double ymin,double ymax,std::string xaxis, std::string yaxis);
 
   MonitorElement* book1dHisto(const std::string &name,const std::string &title,int n,double xmin,double xmax);
   MonitorElement* book2dHisto(const std::string &name,const std::string &title,int nx,double xmin,double xmax,int ny,double ymin,double ymax);

Validation/EventGenerator/interface/TauValidation.h

@@ -23,8 +23,6 @@
 #include "SimGeneral/HepPDTRecord/interface/ParticleDataTable.h"
 #include "TLorentzVector.h"
 
-//#include "Validation/EventGenerator/interface/WeightManager.h"
-
 class TauValidation : public DQMEDAnalyzer
 {
     public: