| 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" |
| 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, bool active) : |
| 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 |
< |
edmSimAssociationName_(Conf().getUntrackedParameter<string>("edmSimAssociationName","")), |
| 31 |
< |
simMapName_(Conf().getUntrackedParameter<string>("simMapName","SimMap")), |
| 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 |
< |
simMap_(0), |
| 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 |
|
//-------------------------------------------------------------------------------------------------- |
| 60 |
|
{ |
| 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 |
< |
simMap_ = OS()->get<SimParticleMap>(simMapName_.c_str()); |
| 67 |
< |
OS()->add<TrackMap>(trackMap_,trackMapName_.c_str()); |
| 68 |
< |
OS()->add<TrackArr>(tracks_,mitName_.c_str()); |
| 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 |
|
//-------------------------------------------------------------------------------------------------- |
| 90 |
|
{ |
| 91 |
|
// Fill tracks from edm collection into our collection. |
| 92 |
|
|
| 93 |
< |
tracks_ ->Reset(); |
| 93 |
> |
tracks_ ->Delete(); |
| 94 |
|
trackMap_->Reset(); |
| 95 |
|
|
| 96 |
< |
Handle<reco::TrackCollection> hTrackProduct; |
| 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(hTrackProduct.id().id()); |
| 102 |
< |
const reco::TrackCollection inTracks = *(hTrackProduct.product()); |
| 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 |
|
// 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 |
< |
GetProduct(edmSimAssociationName_, simAssociationProduct, event); |
| 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 |
+ |
// 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 (reco::TrackCollection::const_iterator it = inTracks.begin(); |
| 128 |
< |
it != inTracks.end(); ++it) { |
| 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->phi(),it->d0(),it->pt(),it->dz(),it->theta()); |
| 133 |
< |
outTrack->SetErrors(it->phiError(),it->d0Error(),it->ptError(),it->dzError(),it->thetaError()); |
| 134 |
< |
outTrack->SetCharge(it->charge()); |
| 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 |
+ |
// fill track quality information |
| 140 |
+ |
outTrack->SetChi2(it->chi2()); |
| 141 |
+ |
outTrack->SetNdof(static_cast<Int_t>(it->ndof())); |
| 142 |
+ |
|
| 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(hTrackProduct, it - 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; |
| 290 |
|
try { |
| 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) |
| 104 |
– |
|
| 302 |
|
if (simRefPair->second > 0.5) // require more than 50% shared hits between reco and sim |
| 303 |
< |
outTrack->SetMCPart(simMap_->GetMit(simRefPair->first)); //add reco->sim reference |
| 303 |
> |
outTrack->SetMCPart(trackingMap_->GetMit(simRefPair->first)); //add reco->sim reference |
| 304 |
|
} |
| 305 |
|
} |
| 306 |
|
} |
