Merge pull request cms-sw#20 from mverzett/master

Ready for new ntuples

interface/ElectronNtuple.h

@@ -11,6 +11,7 @@
 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
 #include "DataFormats/Math/interface/Point3D.h"
 #include "FWCore/Framework/interface/Event.h"
+#include <vector>
 class TTree;
 /*namespace reco {
 	class GsfTrackRef;
@@ -60,7 +61,8 @@ class ElectronNtuple {
 	void fill_gen(const reco::GenParticleRef genp);
 	void fill_gsf_trk(const reco::GsfTrackRef trk, const reco::BeamSpot &spot);
 	void fill_preid(const reco::PreId &preid, const reco::BeamSpot &spot, const int num_gsf);
-	void fill_ele(const reco::GsfElectronRef ele, float mvaid_v1=-2, float mvaid_v2=-2, float ele_conv_vtx_fit_prob = -1.);
+	void fill_ele(const reco::GsfElectronRef ele, float mvaid_v1=-2, float mvaid_v2=-2, 
+								float ele_conv_vtx_fit_prob = -1., const std::vector<float>& iso_rings={0., 0., 0., 0.});
 	void fill_supercluster(const reco::GsfElectronRef ele);
 	void fill_ktf_trk(const reco::TrackRef trk, const reco::BeamSpot &spot);
 	void fill_GSF_ECAL_cluster_info(
@@ -238,6 +240,10 @@ class ElectronNtuple {
 	float ele_mvaIdV1_ = -2.;
 	float ele_mvaIdV2_ = -2.;
 	float ele_conv_vtx_fit_prob_ = -1.;
+	float ele_iso01_ = 0.;
+	float ele_iso02_ = 0.;
+	float ele_iso03_ = 0.;
+	float ele_iso04_ = 0.;
 
 	// Bottom up approach	
 	float gsf_ecal_cluster_e_ = -1;
@@ -271,6 +277,8 @@ class ElectronNtuple {
 	float ktf_ecal_cluster_covEtaPhi_ = -42;
 	float ktf_ecal_cluster_covPhiPhi_ = -42;
 	float ktf_ecal_cluster_ematrix_[ECAL_CLUSTER_SIZE][ECAL_CLUSTER_SIZE] = {{0}};
+	float ktf_ecal_cluster_r9_ = -0.1;
+	float ktf_ecal_cluster_circularity_ = -0.1;
 
 	float ktf_hcal_cluster_e_ = -1;
 	float ktf_hcal_cluster_eta_ = -1;

macros/datasets.py

@@ -8,6 +8,7 @@
 import socket
 path = ""
 if "cern.ch" in socket.gethostname() : path = '/eos/cms/store/cmst3/group/bpark/electron_training'
+elif 'cmg-gpu1080' in socket.gethostname() : path = '/eos/cms/store/cmst3/group/bpark/electron_training'
 elif "hep.ph.ic.ac.uk" in socket.gethostname() : path = '/vols/cms/bainbrid/BParking/electron_training'
 print socket.gethostname()
 

macros/features.py

@@ -212,7 +212,7 @@ def get_features(ftype):
       features = mva_id_inputs
       additional = id_additional
    elif ftype == 'combined_id':
-      features = list(set(mva_id_inputs+id_features)-to_drop-useless)
+      features = list(set(mva_id_inputs+id_features))#-to_drop-useless)
       additional = id_additional
    else:
       raise ValueError('%s is not among the possible feature collection' % ftype)

plugins/TrackerElectronsFeatures.cc

@@ -115,6 +115,7 @@ class TrackerElectronsFeatures: public edm::EDAnalyzer {
   std::pair<uint,uint> indices( float logpt, float eta );
   bool empty_weights();
   void print_weights();
+	vector<float> get_iso_rings(const GsfElectronRef& ele, const reco::TrackRef& ele_trk, const vector<reco::TrackRef> &tracks);
 
   std::pair<float,float> printPfBlock( const reco::GenParticleRef gen,
 				       const reco::PreIdRef preid,
@@ -210,6 +211,19 @@ TrackerElectronsFeatures::TrackerElectronsFeatures(const ParameterSet& cfg):
 //
 // member functions
 //
+vector<float> 
+TrackerElectronsFeatures::get_iso_rings(const GsfElectronRef& ele, const reco::TrackRef& ele_trk, const vector<reco::TrackRef> &tracks) {
+	vector<float> ret = {0., 0., 0., 0.};
+	for(const auto& trk : tracks) {
+		if(ele_trk == trk) continue;
+	 	double dr = deltaR(*ele, *trk);
+	 	size_t idx = std::floor(dr/0.1);
+	 	if(idx < ret.size())
+			ret[idx] += trk->pt(); //this is the conflicting line, for some reason
+	}
+	return ret;
+}
+
 PFClusterRef 
 TrackerElectronsFeatures::closest_cluster(const reco::PFTrajectoryPoint& point, const edm::Handle<reco::PFClusterCollection>& clusters) {
 	PFClusterRef best_ref;
@@ -296,8 +310,12 @@ TrackerElectronsFeatures::analyze(const Event& iEvent, const EventSetup& iSetup)
 	iEvent.getByToken(convVtxFitProb_, convVtxFitProb);
 
 	int ntrks = 0;
+	std::vector<reco::TrackRef> tracks;
 	for ( unsigned int ii = 0; ii < preids->size(); ++ii ) {
-	  if ( preids.product()->at(ii).trackRef().isNonnull() ) { ++ntrks; }
+	  if ( preids.product()->at(ii).trackRef().isNonnull() ) { 
+			tracks.push_back(preids.product()->at(ii).trackRef());
+			++ntrks; 
+		}
 	}
 
 	//Match gen to seed
@@ -652,7 +670,8 @@ TrackerElectronsFeatures::analyze(const Event& iEvent, const EventSetup& iSetup)
           float id1 = (*mvaid_v1)[ele_match->second];
           float id2 = (*mvaid_v2)[ele_match->second];
           float conv_vtx_fit_prob = (*convVtxFitProb)[ele_match->second];
-          ntuple_.fill_ele(ele_match->second, id1, id2, conv_vtx_fit_prob);
+					vector<float> iso_rings = get_iso_rings(ele_match->second, ktf, tracks);
+          ntuple_.fill_ele(ele_match->second, id1, id2, conv_vtx_fit_prob, iso_rings);
 					ele_ref = &(ele_match->second);
         } //matched to GED Electron
 
@@ -753,7 +772,8 @@ TrackerElectronsFeatures::analyze(const Event& iEvent, const EventSetup& iSetup)
 				float id1 = (*mvaid_v1)[ele_match->second];
 				float id2 = (*mvaid_v2)[ele_match->second];
 				float conv_vtx_fit_prob = (*convVtxFitProb)[ele_match->second];
-				ntuple_.fill_ele(ele_match->second, id1, id2, conv_vtx_fit_prob);
+				vector<float> iso_rings = get_iso_rings(ele_match->second, ktf, tracks);
+				ntuple_.fill_ele(ele_match->second, id1, id2, conv_vtx_fit_prob, iso_rings);
 			} //matched to GED Electron
 		}//matched to GSF Track
 
@@ -860,7 +880,8 @@ TrackerElectronsFeatures::analyze(const Event& iEvent, const EventSetup& iSetup)
 				float id1 = (*mvaid_v1)[ele_match->second];
 				float id2 = (*mvaid_v2)[ele_match->second];
 				float conv_vtx_fit_prob = (*convVtxFitProb)[ele_match->second];
-				ntuple_.fill_ele(ele_match->second, id1, id2, conv_vtx_fit_prob);
+				vector<float> iso_rings = get_iso_rings(ele_match->second, ktf, tracks);
+				ntuple_.fill_ele(ele_match->second, id1, id2, conv_vtx_fit_prob, iso_rings);
 			} //matched to GED Electron
 		}//matched to GSF Track
 

production/mc/crab_config.py

@@ -40,7 +40,7 @@
 request = '%s_%s' % (date, name)
 idx = 2
 while True:
-   if not os.path.isdir(request):
+   if not os.path.isdir('crab_%s' % request):
       break
    request = '%sv%d_%s' % (date, idx, name)
    idx += 1

src/ElectronNtuple.cc

@@ -128,6 +128,10 @@ void ElectronNtuple::link_tree(TTree *tree) {
 	tree->Branch("ele_mvaIdV1",		 	  &ele_mvaIdV1_     		 , "ele_mvaIdV1/f");
 	tree->Branch("ele_mvaIdV2",		 	  &ele_mvaIdV2_     		 , "ele_mvaIdV2/f");
 	tree->Branch("ele_conv_vtx_fit_prob", &ele_conv_vtx_fit_prob_, "ele_conv_vtx_fit_prob/f");
+	tree->Branch("ele_iso01", &ele_iso01_, "ele_iso01/f");
+	tree->Branch("ele_iso02", &ele_iso02_, "ele_iso02/f");
+	tree->Branch("ele_iso03", &ele_iso03_, "ele_iso03/f");
+	tree->Branch("ele_iso04", &ele_iso04_, "ele_iso04/f");
 
 	//Bottom up approach
 	tree->Branch("gsf_ecal_cluster_e", &gsf_ecal_cluster_e_, "gsf_ecal_cluster_e/f");
@@ -163,6 +167,8 @@ void ElectronNtuple::link_tree(TTree *tree) {
 	tree->Branch("ktf_ecal_cluster_covEtaPhi", &ktf_ecal_cluster_covEtaPhi_, "ktf_ecal_cluster_covEtaPhi/f");
 	tree->Branch("ktf_ecal_cluster_covPhiPhi", &ktf_ecal_cluster_covPhiPhi_, "ktf_ecal_cluster_covPhiPhi/f");
 	tree->Branch("ktf_ecal_cluster_ematrix", &ktf_ecal_cluster_ematrix_, ("ktf_ecal_cluster_ematrix"+buffer.str()).c_str());
+	tree->Branch("ktf_ecal_cluster_r9", &ktf_ecal_cluster_r9_, "ktf_ecal_cluster_r9/f");
+	tree->Branch("ktf_ecal_cluster_circularity_", &ktf_ecal_cluster_circularity_, "ktf_ecal_cluster_circularity/f");
 
 	tree->Branch("ktf_hcal_cluster_e",    &ktf_hcal_cluster_e_,    "ktf_hcal_cluster_e/f");
 	tree->Branch("ktf_hcal_cluster_eta",  &ktf_hcal_cluster_eta_,  "ktf_hcal_cluster_eta/f");
@@ -380,14 +386,18 @@ void ElectronNtuple::fill_preid( const PreId &preid, const reco::BeamSpot &spot,
 	preid_mva_pass_ = preid.mvaSelected();
 }
 
-void ElectronNtuple::fill_ele(const reco::GsfElectronRef ele, float mvaid_v1, float mvaid_v2, float ele_conv_vtx_fit_prob ) {
+void ElectronNtuple::fill_ele(const reco::GsfElectronRef ele, float mvaid_v1, float mvaid_v2, float ele_conv_vtx_fit_prob, const std::vector<float>& iso_rings) {
 	ele_p_			 = ele->p();
 	ele_pt_			 = ele->pt();
 	ele_eta_		 = ele->eta();
 	ele_phi_     = ele->phi();
 	ele_mvaIdV1_ = mvaid_v1;
 	ele_mvaIdV2_ = mvaid_v2;
 	ele_conv_vtx_fit_prob_ = ele_conv_vtx_fit_prob;
+	ele_iso01_ = iso_rings.at(0);
+	ele_iso02_ = iso_rings.at(1);
+	ele_iso03_ = iso_rings.at(2);
+	ele_iso04_ = iso_rings.at(3);
 	fill_supercluster(ele);
 }
 
@@ -655,6 +665,10 @@ void ElectronNtuple::fill_KTF_ECAL_cluster_info(
 	ktf_ecal_cluster_covEtaPhi_ = covs[1];
 	ktf_ecal_cluster_covPhiPhi_ = covs[2];
 
+	ktf_ecal_cluster_r9_ = ktf_ecal_cluster_e3x3_/ktf_ecal_cluster_e_;
+	float e1x5 = tools.e1x5(*cluster);
+	ktf_ecal_cluster_circularity_ = (ktf_ecal_cluster_e5x5_ > 0) ? 1 - e1x5/ktf_ecal_cluster_e5x5_ : -0.1;
+
 	int cluster_window = (ECAL_CLUSTER_SIZE-1)/2;
 	DetId seedid = cluster->hitsAndFractions().front().first;
 	auto energies = tools.fullMatrixEnergy(