PingTan/src/EdmAnalysisTools.cc

#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/Common/interface/RefToBase.h" 

#include "DataFormats/Math/interface/Vector3D.h"
#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"

#include "DataFormats/CaloTowers/interface/CaloTower.h"
#include "DataFormats/BTauReco/interface/JetTag.h"
#include "DataFormats/JetReco/interface/Jet.h"


#include <DataFormats/VertexReco/interface/Vertex.h>
#include <DataFormats/VertexReco/interface/VertexFwd.h>
#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"

// Math
#include "Math/GenVector/VectorUtil.h"
#include "Math/GenVector/PxPyPzE4D.h"
#include "TVector3.h"

#include "EdmAnalysisTools.h"


#define TRACK_CHI2 5

using namespace std;
using namespace edm;
using namespace reco;


// jet tagging
float btaggingAssociation(Jet jet, const reco::JetTagCollection *btags, float matching_deltaR) {

  // float taginfo = -9999;

  for (unsigned int i = 0; i != (*btags).size(); ++i) {

    float  delta = ROOT::Math::VectorUtil::DeltaR(  (*btags)[i].first->momentum(), jet.momentum() );

    if (delta < matching_deltaR) return (*btags)[i].second;
  }

  return -9999;
}

// jet MC flavor
float mcflavorAssociation(Jet jet, edm::Handle<reco::JetFlavourMatchingCollection> tagList, float matching_deltaR) {


  for ( JetFlavourMatchingCollection::const_iterator j  = tagList->begin();
        j != tagList->end();
        j ++ ) {


    float  delta = ROOT::Math::VectorUtil::DeltaR(  (*j).first->momentum(), jet.momentum() );
    
    if (delta < matching_deltaR) return (*j).second.getFlavour()
;
 }

  return -9999;
}


// still use the global muons only at this moment
void calPtRel(reco::CaloJetCollection::const_iterator jet, const reco::MuonCollection *recoMuons, _jet_ *myJet, float coneSize, float scale) {

  double delta = 0, muonpt = 0;

  Int_t muonIndex = -1;


  for (reco::MuonCollection::const_iterator muon = (*recoMuons).begin(); muon != (*recoMuons).end(); muon ++) {

    muonIndex ++;
    if (!muon->isGlobalMuon()) {continue;}

    delta = ROOT::Math::VectorUtil::DeltaR( (*(muon->track())).momentum(), jet->momentum() );
    if (delta < coneSize) {

      //      myJet->muonNum ++;
      if ((*(muon->track())).pt() < muonpt) continue;

      muonpt = (*(muon->track())).pt();

      //      myJet->muonIndex = muonIndex;
      // myJet->muonType  = muon->type();
      //myJet->muonPt    = muonpt;
      //myJet->deltaR    = delta;
      //myJet->muonChi2  = (*(muon->track())).normalizedChi2();


      // calculate the ptrel
      TVector3 tmp_muon( (*(muon->track())).momentum().X(),
                         (*(muon->track())).momentum().Y(),
                         (*(muon->track())).momentum().Z()) ;
      TVector3 tt( scale*jet->p4().Vect().X(),
                   scale*jet->p4().Vect().Y(),
                   scale*jet->p4().Vect().Z());


      tt+= tmp_muon;
      //  myJet->ptRel =    tmp_muon.Perp(tt);
    }
  }
}


double calIsolation( const CaloTowerCollection   *caloTowers, reco::MuonCollection::const_iterator muon, double threshold) {

  double isolation = 0;
  if (!caloTowers) return isolation;
  for (CaloTowerCollection::const_iterator atower = (*caloTowers).begin(); atower != (*caloTowers).end(); atower ++) {
  
    double delta  = ROOT::Math::VectorUtil::DeltaR(muon->momentum(), atower->momentum() );

    if (delta < threshold) {

      isolation += (atower->emEt() +  atower->hadEt());
    }
  }

  //  isolation -= (muon->calEnergy().em + muon->calEnergy().had) * sin( muon->theta() );
  return isolation;
}


double calIsolation( const CaloTowerCollection   *caloTowers, reco::GsfElectronCollection::const_iterator electron, double threshold) {

  double isolation = 0;
  if (!caloTowers) return isolation;
  for (CaloTowerCollection::const_iterator atower = (*caloTowers).begin(); atower != (*caloTowers).end(); atower ++) {
  
    double delta  = ROOT::Math::VectorUtil::DeltaR(electron->gsfTrack()->momentum(), atower->momentum() );

    if (delta < threshold) {

      isolation += (atower->emEt() +  atower->hadEt());
    }
  }

  return isolation;
}


void calMet( const CaloTowerCollection   *caloTowers,   const reco::MuonCollection *muons, double &px, double &py) {


  if (!caloTowers || !muons) return;
  px = 0;
  py = 0;

  for (CaloTowerCollection::const_iterator atower = (*caloTowers).begin(); atower != (*caloTowers).end(); atower ++) {
 
    px += atower->momentum().x();
    py += atower->momentum().y();
  }


  for (reco::MuonCollection::const_iterator amuon = (*muons).begin(); amuon != (*muons).end(); amuon ++) {

    //    emenergy = amuon->getCalEnergy().em + amuon->getCalEnergy().had + amuon->getCalEnergy().ho;
    if (!amuon->isGlobalMuon()) continue;
    if ((*(amuon->track())).normalizedChi2() > TRACK_CHI2 ) continue;
    
    px += amuon->px();
    py += amuon->py();
    
  }

  px = -px;
  py = -py;
}


// track isolation in trackers
// only tracks with chi2 <5 & nhits >=7 are used in the calculation
float trackIsolation(const reco::Track & theTrack,  const reco::TrackCollection &tracks, 
                     double threshold, double larger_th, 
                     int *nTracks, 
                     float *larger_iso, int *larger_nTrks){

  double isolation = 0;
  (*nTracks) = 0;
  (*larger_iso) =0;
  (*larger_nTrks) = 0;


  for (reco::TrackCollection::const_iterator track = tracks.begin(); track != tracks.end(); track ++) {


    if (track->chi2() >=TRACK_CHI2 ) continue;

    double delta  = ROOT::Math::VectorUtil::DeltaR(theTrack.momentum(), track->momentum() );

    if (delta < larger_th && delta >= 1e-4) {

      (*larger_iso)   += track->pt();
      (*larger_nTrks) ++;

      if (delta < threshold) {
        isolation += track->pt();
        (*nTracks) ++;
      }
    }
  }

  return isolation;
}


// track isolation in tracker with Gsf tracks
// only tracks with chi2 < TRACK_CHI2are used in the calculation
float trackIsolation(const reco::GsfTrack & theTrack,  const reco::TrackCollection &tracks, 
                     double threshold, double larger_th, 
                     int *nTracks, 
                     float *larger_iso, int *larger_nTrks){

  double isolation = 0;
  (*nTracks) = 0;
  (*larger_iso) =0;
  (*larger_nTrks) = 0;


  for (reco::TrackCollection::const_iterator track = tracks.begin(); track != tracks.end(); track ++) {


    if (track->chi2() >=TRACK_CHI2) continue;

    double delta  = ROOT::Math::VectorUtil::DeltaR(theTrack.momentum(), track->momentum() );
    
    if (delta < larger_th && delta >= 1e-3) {
      
      (*larger_iso)   += track->pt();
      (*larger_nTrks) ++;
      
      if (delta < threshold) {
        isolation += track->pt();
        (*nTracks) ++;
      }
    }
  }

  return isolation;
}


float jetIsolation(const reco::Jet & jet, const reco::CaloJetCollection & jets, double threshold) {

  double isolation = 10;
  for (CaloJetCollection::const_iterator aJet = jets.begin(); aJet != jets.end(); aJet ++) {


    if (aJet->pt() < threshold) continue;

    double delta  = ROOT::Math::VectorUtil::DeltaR(jet.p4().Vect(), aJet->p4().Vect() );

    if (delta < 0.001) continue; // skip the identical jet. 
    if (delta < isolation) isolation = delta;
  }

  return isolation;
}

// use with caution: ...
// requires the sim track has minimum 5 GeV pt value
const SimTrack *closestSimTrack(const reco::Track &atrack, const edm::SimTrackContainer *simTrkColl) {


  double predelta = +10;
  const SimTrack *matchedOne = 0;
  if (!simTrkColl) return matchedOne;

  for (SimTrackContainer::const_iterator gentrk = (*simTrkColl).begin(); gentrk != (*simTrkColl).end(); gentrk++) {

    //    double delta  = ROOT::Math::VectorUtil::DeltaR( TVector3((*gentrk).momentum().x(),(*gentrk).momentum().y(),(*gentrk).momentum().z())
    double delta  = ROOT::Math::VectorUtil::DeltaR((*gentrk).momentum().Vect()
                                                   ,atrack.momentum() );
    if (delta < predelta &&(*gentrk).momentum().Pt()>5 ) {
      predelta = delta;
      matchedOne =&(*gentrk);
    }
  }  

  return matchedOne;
}


const reco::Vertex  *vertexAssociation(const TrackBaseRef &track, Handle<reco::VertexCollection> recVtxs) {


  for (reco::VertexCollection::const_iterator v=recVtxs->begin(); 
       v!=recVtxs->end(); 
       ++v){


    //    if ( v->refittedTracks().empty() ) continue;

    std::vector<TrackBaseRef >::const_iterator it = find(v->tracks_begin(), v->tracks_end(), track);
    if (it!=v->tracks_end()) return &(*v);
  }

  return 0;
}


// !!!!!!!!!!! WARNING: filter is hard-coded here !!!!!!!!!!!!!
int  whichVertex(const TrackBaseRef &track, Handle<reco::VertexCollection> recVtxs, float & distance) {


  int index = -1, associated = -1;
  distance = 99999;
  for (reco::VertexCollection::const_iterator v=recVtxs->begin(); 
       v!=recVtxs->end(); 
       ++v){


    index ++;

    if (  fabs( v->position().z() ) > 15 || 
          v->position().rho() >2 ||
          v->isFake() ||
          v->ndof() <=4 ) continue;
    //    if ( v->refittedTracks().empty() ) continue;

    if (distance > track->dz(v->position()) ) {

      distance = track->dz(v->position());
      associated = index;

    }

    //    if (fabs( v->position().z() - track->vz() ) < distance) {
    //      distance = fabs( v->position().z() - track->vz() );
    // associated = index;
    //}
    //    std::vector<TrackBaseRef >::const_iterator it = find(v->tracks_begin(), v->tracks_end(), track);
    //  if (it!=v->tracks_end()) return index;
  }

  return associated;
}
Revision:	1.1
Committed:	Tue Oct 5 17:09:41 2010 UTC (14 years, 6 months ago) by ptan
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Log Message:	* empty log message *
#	User	Rev	Content
1	ptan	1.1	#include "DataFormats/TrackReco/interface/TrackFwd.h"
2			#include "DataFormats/TrackReco/interface/Track.h"
3			#include "DataFormats/Common/interface/RefToBase.h"
4
5			#include "DataFormats/Math/interface/Vector3D.h"
6			#include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
7
8			#include "DataFormats/CaloTowers/interface/CaloTower.h"
9			#include "DataFormats/BTauReco/interface/JetTag.h"
10			#include "DataFormats/JetReco/interface/Jet.h"
11
12
13
14
15			#include <DataFormats/VertexReco/interface/Vertex.h>
16			#include <DataFormats/VertexReco/interface/VertexFwd.h>
17			#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h"
18
19			// Math
20			#include "Math/GenVector/VectorUtil.h"
21			#include "Math/GenVector/PxPyPzE4D.h"
22			#include "TVector3.h"
23
24			#include "EdmAnalysisTools.h"
25
26
27			#define TRACK_CHI2 5
28
29			using namespace std;
30			using namespace edm;
31			using namespace reco;
32
33
34			// jet tagging
35			float btaggingAssociation(Jet jet, const reco::JetTagCollection *btags, float matching_deltaR) {
36
37			// float taginfo = -9999;
38
39			for (unsigned int i = 0; i != (*btags).size(); ++i) {
40
41			float delta = ROOT::Math::VectorUtil::DeltaR( (*btags)[i].first->momentum(), jet.momentum() );
42
43			if (delta < matching_deltaR) return (*btags)[i].second;
44			}
45
46			return -9999;
47			}
48
49			// jet MC flavor
50			float mcflavorAssociation(Jet jet, edm::Handle<reco::JetFlavourMatchingCollection> tagList, float matching_deltaR) {
51
52
53			for ( JetFlavourMatchingCollection::const_iterator j = tagList->begin();
54			j != tagList->end();
55			j ++ ) {
56
57
58			float delta = ROOT::Math::VectorUtil::DeltaR( (*j).first->momentum(), jet.momentum() );
59
60			if (delta < matching_deltaR) return (*j).second.getFlavour()
61			;
62			}
63
64			return -9999;
65			}
66
67
68
69			// still use the global muons only at this moment
70			void calPtRel(reco::CaloJetCollection::const_iterator jet, const reco::MuonCollection recoMuons, _jet_ myJet, float coneSize, float scale) {
71
72			double delta = 0, muonpt = 0;
73
74			Int_t muonIndex = -1;
75
76
77			for (reco::MuonCollection::const_iterator muon = (recoMuons).begin(); muon != (recoMuons).end(); muon ++) {
78
79			muonIndex ++;
80			if (!muon->isGlobalMuon()) {continue;}
81
82			delta = ROOT::Math::VectorUtil::DeltaR( (*(muon->track())).momentum(), jet->momentum() );
83			if (delta < coneSize) {
84
85			// myJet->muonNum ++;
86			if ((*(muon->track())).pt() < muonpt) continue;
87
88			muonpt = (*(muon->track())).pt();
89
90			// myJet->muonIndex = muonIndex;
91			// myJet->muonType = muon->type();
92			//myJet->muonPt = muonpt;
93			//myJet->deltaR = delta;
94			//myJet->muonChi2 = (*(muon->track())).normalizedChi2();
95
96
97			// calculate the ptrel
98			TVector3 tmp_muon( (*(muon->track())).momentum().X(),
99			(*(muon->track())).momentum().Y(),
100			(*(muon->track())).momentum().Z()) ;
101			TVector3 tt( scale*jet->p4().Vect().X(),
102			scale*jet->p4().Vect().Y(),
103			scale*jet->p4().Vect().Z());
104
105
106			tt+= tmp_muon;
107			// myJet->ptRel = tmp_muon.Perp(tt);
108			}
109			}
110			}
111
112
113			double calIsolation( const CaloTowerCollection *caloTowers, reco::MuonCollection::const_iterator muon, double threshold) {
114
115			double isolation = 0;
116			if (!caloTowers) return isolation;
117			for (CaloTowerCollection::const_iterator atower = (caloTowers).begin(); atower != (caloTowers).end(); atower ++) {
118
119			double delta = ROOT::Math::VectorUtil::DeltaR(muon->momentum(), atower->momentum() );
120
121			if (delta < threshold) {
122
123			isolation += (atower->emEt() + atower->hadEt());
124			}
125			}
126
127			// isolation -= (muon->calEnergy().em + muon->calEnergy().had) * sin( muon->theta() );
128			return isolation;
129			}
130
131
132			double calIsolation( const CaloTowerCollection *caloTowers, reco::GsfElectronCollection::const_iterator electron, double threshold) {
133
134			double isolation = 0;
135			if (!caloTowers) return isolation;
136			for (CaloTowerCollection::const_iterator atower = (caloTowers).begin(); atower != (caloTowers).end(); atower ++) {
137
138			double delta = ROOT::Math::VectorUtil::DeltaR(electron->gsfTrack()->momentum(), atower->momentum() );
139
140			if (delta < threshold) {
141
142			isolation += (atower->emEt() + atower->hadEt());
143			}
144			}
145
146			return isolation;
147			}
148
149
150
151			void calMet( const CaloTowerCollection caloTowers, const reco::MuonCollection muons, double &px, double &py) {
152
153
154			if (!caloTowers \|\| !muons) return;
155			px = 0;
156			py = 0;
157
158			for (CaloTowerCollection::const_iterator atower = (caloTowers).begin(); atower != (caloTowers).end(); atower ++) {
159
160			px += atower->momentum().x();
161			py += atower->momentum().y();
162			}
163
164
165			for (reco::MuonCollection::const_iterator amuon = (muons).begin(); amuon != (muons).end(); amuon ++) {
166
167			// emenergy = amuon->getCalEnergy().em + amuon->getCalEnergy().had + amuon->getCalEnergy().ho;
168			if (!amuon->isGlobalMuon()) continue;
169			if ((*(amuon->track())).normalizedChi2() > TRACK_CHI2 ) continue;
170
171			px += amuon->px();
172			py += amuon->py();
173
174			}
175
176			px = -px;
177			py = -py;
178			}
179
180
181
182
183			// track isolation in trackers
184			// only tracks with chi2 <5 & nhits >=7 are used in the calculation
185			float trackIsolation(const reco::Track & theTrack, const reco::TrackCollection &tracks,
186			double threshold, double larger_th,
187			int *nTracks,
188			float larger_iso, int larger_nTrks){
189
190			double isolation = 0;
191			(*nTracks) = 0;
192			(*larger_iso) =0;
193			(*larger_nTrks) = 0;
194
195
196			for (reco::TrackCollection::const_iterator track = tracks.begin(); track != tracks.end(); track ++) {
197
198
199			if (track->chi2() >=TRACK_CHI2 ) continue;
200
201			double delta = ROOT::Math::VectorUtil::DeltaR(theTrack.momentum(), track->momentum() );
202
203			if (delta < larger_th && delta >= 1e-4) {
204
205			(*larger_iso) += track->pt();
206			(*larger_nTrks) ++;
207
208			if (delta < threshold) {
209			isolation += track->pt();
210			(*nTracks) ++;
211			}
212			}
213			}
214
215			return isolation;
216			}
217
218
219			// track isolation in tracker with Gsf tracks
220			// only tracks with chi2 < TRACK_CHI2are used in the calculation
221			float trackIsolation(const reco::GsfTrack & theTrack, const reco::TrackCollection &tracks,
222			double threshold, double larger_th,
223			int *nTracks,
224			float larger_iso, int larger_nTrks){
225
226			double isolation = 0;
227			(*nTracks) = 0;
228			(*larger_iso) =0;
229			(*larger_nTrks) = 0;
230
231
232			for (reco::TrackCollection::const_iterator track = tracks.begin(); track != tracks.end(); track ++) {
233
234
235			if (track->chi2() >=TRACK_CHI2) continue;
236
237			double delta = ROOT::Math::VectorUtil::DeltaR(theTrack.momentum(), track->momentum() );
238
239			if (delta < larger_th && delta >= 1e-3) {
240
241			(*larger_iso) += track->pt();
242			(*larger_nTrks) ++;
243
244			if (delta < threshold) {
245			isolation += track->pt();
246			(*nTracks) ++;
247			}
248			}
249			}
250
251			return isolation;
252			}
253
254
255			float jetIsolation(const reco::Jet & jet, const reco::CaloJetCollection & jets, double threshold) {
256
257			double isolation = 10;
258			for (CaloJetCollection::const_iterator aJet = jets.begin(); aJet != jets.end(); aJet ++) {
259
260
261			if (aJet->pt() < threshold) continue;
262
263			double delta = ROOT::Math::VectorUtil::DeltaR(jet.p4().Vect(), aJet->p4().Vect() );
264
265			if (delta < 0.001) continue; // skip the identical jet.
266			if (delta < isolation) isolation = delta;
267			}
268
269			return isolation;
270			}
271
272			// use with caution: ...
273			// requires the sim track has minimum 5 GeV pt value
274			const SimTrack closestSimTrack(const reco::Track &atrack, const edm::SimTrackContainer simTrkColl) {
275
276
277			double predelta = +10;
278			const SimTrack *matchedOne = 0;
279			if (!simTrkColl) return matchedOne;
280
281			for (SimTrackContainer::const_iterator gentrk = (simTrkColl).begin(); gentrk != (simTrkColl).end(); gentrk++) {
282
283			// double delta = ROOT::Math::VectorUtil::DeltaR( TVector3((gentrk).momentum().x(),(gentrk).momentum().y(),(*gentrk).momentum().z())
284			double delta = ROOT::Math::VectorUtil::DeltaR((*gentrk).momentum().Vect()
285			,atrack.momentum() );
286			if (delta < predelta &&(*gentrk).momentum().Pt()>5 ) {
287			predelta = delta;
288			matchedOne =&(*gentrk);
289			}
290			}
291
292			return matchedOne;
293			}
294
295
296
297			const reco::Vertex *vertexAssociation(const TrackBaseRef &track, Handle<reco::VertexCollection> recVtxs) {
298
299
300
301			for (reco::VertexCollection::const_iterator v=recVtxs->begin();
302			v!=recVtxs->end();
303			++v){
304
305
306			// if ( v->refittedTracks().empty() ) continue;
307
308			std::vector<TrackBaseRef >::const_iterator it = find(v->tracks_begin(), v->tracks_end(), track);
309			if (it!=v->tracks_end()) return &(*v);
310			}
311
312			return 0;
313			}
314
315
316			// !!!!!!!!!!! WARNING: filter is hard-coded here !!!!!!!!!!!!!
317			int whichVertex(const TrackBaseRef &track, Handle<reco::VertexCollection> recVtxs, float & distance) {
318
319
320			int index = -1, associated = -1;
321			distance = 99999;
322			for (reco::VertexCollection::const_iterator v=recVtxs->begin();
323			v!=recVtxs->end();
324			++v){
325
326
327			index ++;
328
329			if ( fabs( v->position().z() ) > 15 \|\|
330			v->position().rho() >2 \|\|
331			v->isFake() \|\|
332			v->ndof() <=4 ) continue;
333			// if ( v->refittedTracks().empty() ) continue;
334
335			if (distance > track->dz(v->position()) ) {
336
337			distance = track->dz(v->position());
338			associated = index;
339
340			}
341
342			// if (fabs( v->position().z() - track->vz() ) < distance) {
343			// distance = fabs( v->position().z() - track->vz() );
344			// associated = index;
345			//}
346			// std::vector<TrackBaseRef >::const_iterator it = find(v->tracks_begin(), v->tracks_end(), track);
347			// if (it!=v->tracks_end()) return index;
348			}
349
350			return associated;
351			}