ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitProd/TreeFiller/src/FillerTracks.cc
Revision: 1.38
Committed: Tue Mar 23 16:52:58 2010 UTC (15 years, 1 month ago) by bendavid
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_013d, Mit_013c, Mit_013b, Mit_013a, Mit_013
Changes since 1.37: +49 -6 lines
Log Message:
Fill expected and missing hit information for tracks

File Contents

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