ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerTracks.cc

Comparing UserCode/MitProd/TreeFiller/src/FillerTracks.cc (file contents):
Revision 1.2 by loizides, Tue Jul 1 21:11:47 2008 UTC vs.
Revision 1.41 by bendavid, Fri Mar 11 04:03:55 2011 UTC

#	Line 1 | Line 1
1		// $Id$
2
3		#include "MitProd/TreeFiller/interface/FillerTracks.h"
4	<	#include "FWCore/MessageLogger/interface/MessageLogger.h"
5	<	#include "DataFormats/Common/interface/Handle.h"
4	>	#include "DataFormats/RecoCandidate/interface/TrackAssociation.h"
5	>	#include "DataFormats/TrackReco/interface/HitPattern.h"
6		#include "DataFormats/TrackReco/interface/Track.h"
7		#include "DataFormats/TrackReco/interface/TrackFwd.h"
8	+	#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
9		#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
10		#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
11	<	#include "DataFormats/RecoCandidate/interface/TrackAssociation.h"
12	<
13	<	#include "MitAna/DataTree/interface/Track.h"
11	>	#include "MagneticField/Engine/interface/MagneticField.h"
12	>	#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
13	>	#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
14	>	#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
15		#include "MitAna/DataTree/interface/Names.h"
16	+	#include "MitAna/DataTree/interface/TrackCol.h"
17	+	#include "MitProd/ObjectService/interface/ObjectService.h"
18	+	#include "MitEdm/DataFormats/interface/Types.h"
19
20		using namespace std;
21		using namespace edm;
22		using namespace mithep;
23
24	<	//-------------------------------------------------------------------------------------------------
25	<	FillerTracks::FillerTracks(const ParameterSet &cfg, const char *name,
26	<	bool active, const SimParticleMap *sm) :
27	<	BaseFiller(cfg, name, active),
28	<	edmName_(Conf().getUntrackedParameter<string>("edmName","")),
29	<	edmDataName_(Conf().getUntrackedParameter<string>("edmDataName","")),
30	<	mitName_(Conf().getUntrackedParameter<string>("mitName","")),
31	<	edmSimAssociationName_(Conf().getUntrackedParameter<string>("edmSimAssociationName","")),
32	<	simMap_(sm),
33	<	tracks_(new mithep::TrackArr),
24	>	//--------------------------------------------------------------------------------------------------
25	>	FillerTracks::FillerTracks(const ParameterSet &cfg, const char *name, bool active) :
26	>	BaseFiller(cfg,name,active),
27	>	ecalAssocActive_(Conf().getUntrackedParameter<bool>("ecalAssocActive",0)),
28	>	edmName_(Conf().getUntrackedParameter<string>("edmName","generalTracks")),
29	>	mitName_(Conf().getUntrackedParameter<string>("mitName",Names::gkTrackBrn)),
30	>	edmSimAssocName_(Conf().getUntrackedParameter<string>("edmSimAssociationName","")),
31	>	trackingMapName_(Conf().getUntrackedParameter<string>("trackingMapName","TrackingMap")),
32	>	barrelSuperClusterIdMapName_(Conf().getUntrackedParameter<string>
33	>	("superClusterIdMapName","barrelSuperClusterIdMap")),
34	>	endcapSuperClusterIdMapName_(Conf().getUntrackedParameter<string>
35	>	("endcapClusterIdMapName","endcapSuperClusterIdMap")),
36	>	trackMapName_(Conf().getUntrackedParameter<string>("trackMapName",
37	>	Form("%sMapName",mitName_.c_str()))),
38	>	trackingMap_(0),
39	>	tracks_(new mithep::TrackArr(250)),
40		trackMap_(new mithep::TrackMap)
41		{
42		// Constructor.
43	+
44	+	if (ecalAssocActive_) // initialize track associator configuration if needed
45	+	assocParams_.loadParameters(
46	+	cfg.getUntrackedParameter<ParameterSet>("TrackAssociatorParameters"));
47		}
48
49	<	//-------------------------------------------------------------------------------------------------
49	>	//--------------------------------------------------------------------------------------------------
50		FillerTracks::~FillerTracks()
51		{
52		// Destructor.
53
54	+	delete tracks_;
55		delete trackMap_;
56		}
57
58	<	//-------------------------------------------------------------------------------------------------
58	>	//--------------------------------------------------------------------------------------------------
59		void FillerTracks::BookDataBlock(TreeWriter &tws)
60		{
61	<	// Add tracks branch to tree.
61	>	// Add tracks branch to tree, publish and get our objects.
62
63	<	tws.AddBranch(mitName_.c_str(),&tracks_);
63	>	tws.AddBranch(mitName_,&tracks_);
64	>	OS()->add<TrackArr>(tracks_,mitName_);
65	>
66	>	if (!trackingMapName_.empty()) {
67	>	trackingMap_ = OS()->get<TrackingParticleMap>(trackingMapName_);
68	>	if (trackingMap_)
69	>	AddBranchDep(mitName_,trackingMap_->GetBrName());
70	>	}
71	>	if (ecalAssocActive_ && !barrelSuperClusterIdMapName_.empty()) {
72	>	barrelSuperClusterIdMap_ = OS()->get<SuperClusterIdMap>(barrelSuperClusterIdMapName_);
73	>	if (barrelSuperClusterIdMap_)
74	>	AddBranchDep(mitName_,barrelSuperClusterIdMap_->GetBrName());
75	>	}
76	>	if (ecalAssocActive_ && !endcapSuperClusterIdMapName_.empty()) {
77	>	endcapSuperClusterIdMap_ = OS()->get<SuperClusterIdMap>(endcapSuperClusterIdMapName_);
78	>	if (endcapSuperClusterIdMap_)
79	>	AddBranchDep(mitName_,endcapSuperClusterIdMap_->GetBrName());
80	>	}
81	>	if (!trackMapName_.empty()) {
82	>	trackMap_->SetBrName(mitName_);
83	>	OS()->add<TrackMap>(trackMap_,trackMapName_);
84	>	}
85		}
86
87	<	//-------------------------------------------------------------------------------------------------
87	>	//--------------------------------------------------------------------------------------------------
88		void FillerTracks::FillDataBlock(const edm::Event      &event,
89		const edm::EventSetup &setup)
90		{
91		// Fill tracks from edm collection into our collection.
92
93	<	tracks_->Reset();
93	>	tracks_  ->Delete();
94		trackMap_->Reset();
95	<
96	<	// get the tracks collection
97	<	try {
98	<	if (edmDataName_ == "")
99	<	event.getByLabel(edmName_,trackProduct_);
63	<	else
64	<	event.getByLabel(edmName_,edmDataName_,trackProduct_);
65	<	} catch (cms::Exception& ex) {
66	<	edm::LogError("FillerTracks") << "Error! Cannot get collection with label "
67	<	<< edmName_ << endl;
68	<	throw edm::Exception(edm::errors::Configuration, "FillerTracks:FillDataBlock()\n")
69	<	<< "Error! Cannot get collection with label " << edmName_ << endl;
70	<	}
95	>
96	>	// initialize handle and get product,
97	>	// this usage allows also to get collections of classes which inherit from reco::Track
98	>	Handle<View<reco::Track> > hTrackProduct;
99	>	GetProduct(edmName_, hTrackProduct, event);
100
101	<	trackMap_->SetEdmProductId(trackProduct_.id().id());
102	<	const reco::TrackCollection inTracks = *(trackProduct_.product());
101	>	trackMap_->SetEdmProductId(hTrackProduct.id().id());
102	>	const View<reco::Track> inTracks = *(hTrackProduct.product());
103	>
104	>	if (verbose_>1)
105	>	printf("Track Collection: %s, id=%i\n",edmName_.c_str(),hTrackProduct.id().id());
106
107	<	// if we have a Sim Particle association (for monte carlo), initialize the reco->sim mappings
107	>	// for MC SimParticle association (reco->sim mappings)
108		reco::RecoToSimCollection simAssociation;
109	<	if (simMap_ && !edmSimAssociationName_.empty()) {
109	>	if (trackingMap_ && !edmSimAssocName_.empty()) {
110		Handle<reco::RecoToSimCollection> simAssociationProduct;
111	<	try {
80	<	event.getByLabel(edmSimAssociationName_, simAssociationProduct);
81	<	}
82	<	catch (cms::Exception& ex) {
83	<	edm::LogError("FillerTracks") << "Error! Cannot get collection with label "
84	<	<< edmSimAssociationName_ << endl;
85	<	throw edm::Exception(edm::errors::Configuration, "FillerTracks:FillDataBlock()\n")
86	<	<< "Error! Cannot get collection with label " << edmSimAssociationName_ << endl;
87	<	}
111	>	GetProduct(edmSimAssocName_, simAssociationProduct, event);
112		simAssociation = *(simAssociationProduct.product());
113	+	if (verbose_>1)
114	+	printf("SimAssociation Map Size = %i\n",(int)simAssociation.size());
115		}
116
117	<	// loop through all tracks
118	<	for (reco::TrackCollection::const_iterator inTrack = inTracks.begin();
119	<	inTrack != inTracks.end(); ++inTrack) {
120	<
121	<	mithep::Track* outTrack = tracks_->Allocate();
122	<	new (outTrack) mithep::Track(inTrack->phi(),inTrack->d0(),inTrack->pt(),inTrack->dz(),inTrack->theta());
123	<
117	>	// set up associator for Track-Ecal associations
118	>	TrackDetectorAssociator trackAssociator;
119	>	trackAssociator.useDefaultPropagator();
120	>	edm::ESHandle<MagneticField> bField;
121	>	if (ecalAssocActive_) {
122	>	setup.get<IdealMagneticFieldRecord>().get(bField);
123	>	}
124	>
125	>
126	>	// loop through all tracks and fill the information
127	>	for (View<reco::Track>::const_iterator it = inTracks.begin();
128	>	it != inTracks.end(); ++it) {
129	>
130	>	mithep::Track *outTrack = tracks_->Allocate();
131	>	// create track and set the core parameters
132	>	new (outTrack) mithep::Track(it->qoverp(),it->lambda(),
133	>	it->phi(),it->dxy(),it->dsz());
134	>	outTrack->SetErrors(it->qoverpError(),it->lambdaError(),
135	>	it->phiError(),it->dxyError(),it->dszError());
136	>
137	>	outTrack->SetPtErr(it->ptError());
138
139	<	outTrack->SetErrors(inTrack->phiError(),
140	<	inTrack->d0Error(),
141	<	inTrack->ptError(),
102	<	inTrack->dzError(),
103	<	inTrack->thetaError());
139	>	// fill track quality information
140	>	outTrack->SetChi2(it->chi2());
141	>	outTrack->SetNdof(static_cast<Int_t>(it->ndof()));
142
143	<	outTrack->SetCharge(inTrack->charge());
144	<
143	>	//fill algo enum
144	>	outTrack->SetAlgo(static_cast<mithep::Track::ETrackAlgorithm>(it->algo()));
145	>
146	>	//fill quality bitmask
147	>	outTrack->Quality().SetQualityMask(BitMask8(UChar_t(it->qualityMask())));
148	>
149	>	if (verbose_>2) {
150	>	printf("reco::Track   high purity = %i,  ", it->quality(reco::TrackBase::highPurity));
151	>	printf("mithep::Track high purity = %i\n",outTrack->Quality().Quality(mithep::TrackQuality::highPurity));
152	>	}
153	>
154	>	//fill gsf flag, some type gymastics needed...
155	>	if (typeid(*it)==typeid(reco::GsfTrack))
156	>	outTrack->SetIsGsf(kTRUE);
157	>	else
158	>	outTrack->SetIsGsf(kFALSE);
159	>
160	>	if (verbose_>2) {
161	>	printf("Track pt = %5f, eta = %5f, phi = %5f, nhits = %i, numberofvalidhits = %i, numberofmissinghits = %i, expectedhitsinner = %i, expectedhitsouter = %i\n",it->pt(),it->eta(),it->phi(),it->hitPattern().numberOfHits(),it->numberOfValidHits(),it->hitPattern().numberOfLostHits(),it->trackerExpectedHitsInner().numberOfHits(),it->trackerExpectedHitsOuter().numberOfHits());
162	>	}
163	>
164	>	//fill hit layer map and missing hits
165	>	const reco::HitPattern &hits = it->hitPattern();
166	>	BitMask48 hitLayers;
167	>	BitMask48 missingHitLayers;
168	>	for (Int_t i=0; i<hits.numberOfHits(); i++) {
169	>	uint32_t hit = hits.getHitPattern(i);
170	>	if (hits.validHitFilter(hit)) {
171	>	if (hits.trackerHitFilter(hit)) {
172	>	hitLayers.SetBit(hitReader_.Layer(hit));
173	>	if (verbose_>2)
174	>	printf("Found valid hit with layer %i\n",hitReader_.Layer(hit));
175	>	}
176	>	}
177	>
178	>	if (hits.getHitType(hit)==1) {
179	>	if (hits.trackerHitFilter(hit)) {
180	>	missingHitLayers.SetBit(hitReader_.Layer(hit));
181	>	}
182	>	}
183	>
184	>	if (verbose_>2) {
185	>	if (hits.muonDTHitFilter(hit))
186	>	printf("Muon DT Layer = %i\n", hits.getLayer(hit));
187	>	if (hits.muonCSCHitFilter(hit))
188	>	printf("Muon CSC Layer = %i\n", hits.getLayer(hit));
189	>	if (hits.muonRPCHitFilter(hit))
190	>	printf("Muon RPC Layer = %i\n", hits.getLayer(hit));
191	>	}
192	>	}
193	>
194	>	outTrack->SetHits(hitLayers);
195	>	outTrack->SetMissingHits(missingHitLayers);
196	>
197	>
198	>	//set expected inner hits
199	>	const reco::HitPattern &expectedInnerHitPattern = it->trackerExpectedHitsInner();
200	>	BitMask48 expectedHitsInner;
201	>	// search for all good crossed layers (with or without hit)
202	>	for (Int_t hi=0; hi<expectedInnerHitPattern.numberOfHits(); hi++) {
203	>	uint32_t hit = expectedInnerHitPattern.getHitPattern(hi);
204	>	if (expectedInnerHitPattern.getHitType(hit)<=1) {
205	>	if (expectedInnerHitPattern.trackerHitFilter(hit)) {
206	>	expectedHitsInner.SetBit(hitReader_.Layer(hit));
207	>	}
208	>	}
209	>	}
210	>
211	>	outTrack->SetExpectedHitsInner(expectedHitsInner);
212	>
213	>	//set expected outer hits
214	>	const reco::HitPattern &expectedOuterHitPattern = it->trackerExpectedHitsOuter();
215	>	BitMask48 expectedHitsOuter;
216	>	// search for all good crossed layers (with or without hit)
217	>	for (Int_t hi=0; hi<expectedOuterHitPattern.numberOfHits(); hi++) {
218	>	uint32_t hit = expectedOuterHitPattern.getHitPattern(hi);
219	>	if (expectedOuterHitPattern.getHitType(hit)<=1) {
220	>	if (expectedOuterHitPattern.trackerHitFilter(hit)) {
221	>	expectedHitsOuter.SetBit(hitReader_.Layer(hit));
222	>	}
223	>	}
224	>	}
225	>
226	>	outTrack->SetExpectedHitsOuter(expectedHitsOuter);
227	>
228	>	//fill actual hit counts
229	>	outTrack->SetNHits(it->numberOfValidHits());
230	>	outTrack->SetNPixelHits(it->hitPattern().numberOfValidPixelHits());
231	>	outTrack->SetNMissingHits(it->hitPattern().numberOfLostHits());
232	>	outTrack->SetNExpectedHitsInner(it->trackerExpectedHitsInner().numberOfHits());
233	>	outTrack->SetNExpectedHitsOuter(it->trackerExpectedHitsOuter().numberOfHits());
234	>
235	>
236	>	if (verbose_>2)
237	>	printf("OutTrack NHits = %i, NMissingHits = %i, NExpectedHitsInner = %i, NExpectedHitsOuter = %i\n",outTrack->NHits(),outTrack->NMissingHits(),outTrack->NExpectedHitsInner(),outTrack->NExpectedHitsOuter());
238	>
239	>
240	>	//make ecal associations
241	>	if (ecalAssocActive_) {
242	>	TrackDetMatchInfo matchInfo;
243	>	//Extra check to allow propagation to work in AOD where no TrackExtra is available
244	>	if (it->extra().isAvailable())
245	>	matchInfo = trackAssociator.associate(event,setup,*it,assocParams_);
246	>	else {
247	>	TrajectoryStateTransform tsTransform;
248	>	FreeTrajectoryState initialState = tsTransform.initialFreeState(*it,&*bField);
249	>	matchInfo = trackAssociator.associate(event, setup, assocParams_, &initialState);
250	>	}
251	>
252	>	if (matchInfo.isGoodEcal) {
253	>	outTrack->SetEtaEcal(matchInfo.trkGlobPosAtEcal.eta());
254	>	outTrack->SetPhiEcal(matchInfo.trkGlobPosAtEcal.phi());
255	>	}
256	>
257	>	//fill supercluster link
258	>	if (barrelSuperClusterIdMap_ || endcapSuperClusterIdMap_) {
259	>	mithep::SuperCluster *cluster = 0;
260	>	for (std::vector<const ::CaloTower*>::const_iterator iTower =
261	>	matchInfo.crossedTowers.begin();
262	>	iTower<matchInfo.crossedTowers.end() && !cluster; iTower++) {
263	>
264	>	for (uint ihit=0; ihit<(*iTower)->constituentsSize() && !cluster; ihit++) {
265	>	DetId hit = (*iTower)->constituent(ihit);
266	>	if (barrelSuperClusterIdMap_ && barrelSuperClusterIdMap_->HasMit(hit))
267	>	cluster = barrelSuperClusterIdMap_->GetMit(hit);
268	>	else if (endcapSuperClusterIdMap_ && endcapSuperClusterIdMap_->HasMit(hit))
269	>	cluster = endcapSuperClusterIdMap_->GetMit(hit);
270	>	}
271	>	}
272	>	if (cluster)
273	>	outTrack->SetSCluster(cluster);
274	>	}
275	>	}
276	>
277		// add reference between mithep and edm object
278	<	reco::TrackRef theRef(trackProduct_, inTrack-inTracks.begin());
279	<	trackMap_->Add(theRef, outTrack);
280	<
281	<	if (simMap_ && !edmSimAssociationName_.empty()) {
282	<	reco::TrackBaseRef theBaseRef(theRef);
278	>	if (trackMap_) {
279	>	mitedm::TrackPtr thePtr = inTracks.ptrAt(it - inTracks.begin());
280	>	trackMap_->Add(thePtr, outTrack);
281	>	}
282	>
283	>	//do sim associations
284	>	if (trackingMap_ && !edmSimAssocName_.empty()) {
285	>	if (verbose_>1)
286	>	printf("Trying Track-Sim association\n");
287	>	reco::TrackBaseRef theBaseRef = inTracks.refAt(it - inTracks.begin());
288		vector<pair<TrackingParticleRef, double> > simRefs;
289	<	Bool_t noSimParticle=0;
289	>	Bool_t noSimParticle = 0;
290		try {
291		simRefs = simAssociation[theBaseRef]; //try to get the sim references if existing
292		}
293	<	catch (edm::Exception& ex) {
294	<	noSimParticle=1;
293	>	catch (edm::Exception &ex) {
294	>	noSimParticle = 1;
295		}
296	+
297		if (!noSimParticle) { //loop through sim match candidates
298	+	if (verbose_>1)
299	+	printf("Applying track-sim association\n");
300		for (vector<pair<TrackingParticleRef, double> >::const_iterator simRefPair=simRefs.begin();
301		simRefPair != simRefs.end(); ++simRefPair)
302	<
303	<	if ( simRefPair->second > 0.5 ) // require more than 50% shared hits between reco and sim
126	<	outTrack->SetSimParticle(simMap_->GetMit(simRefPair->first)); //add reco->sim reference
302	>	if (simRefPair->second > 0.5) // require more than 50% shared hits between reco and sim
303	>	outTrack->SetMCPart(trackingMap_->GetMit(simRefPair->first)); //add reco->sim reference
304		}
305		}
306		}
130	–
307		tracks_->Trim();
308		}

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines