TreeFiller/src/FillerTracks.cc

// $Id: FillerTracks.cc,v 1.29 2009/03/15 11:20:41 loizides Exp $

#include "MitProd/TreeFiller/interface/FillerTracks.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/TrackReco/interface/HitPattern.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "DataFormats/RecoCandidate/interface/TrackAssociation.h"
#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "MitEdm/DataFormats/interface/Types.h"

using namespace std;
using namespace edm;
using namespace mithep;

//--------------------------------------------------------------------------------------------------
FillerTracks::FillerTracks(const ParameterSet &cfg, const char *name, bool active) :
  BaseFiller(cfg,name,active),
  ecalAssocActive_(Conf().getUntrackedParameter<bool>("ecalAssocActive",0)),
  edmName_(Conf().getUntrackedParameter<string>("edmName","generalTracks")),
  mitName_(Conf().getUntrackedParameter<string>("mitName",Names::gkTrackBrn)),
  edmSimAssocName_(Conf().getUntrackedParameter<string>("edmSimAssociationName","")),
  trackingMapName_(Conf().getUntrackedParameter<string>("trackingMapName","TrackingMap")),
  barrelSuperClusterIdMapName_(Conf().getUntrackedParameter<string>
                               ("superClusterIdMapName","barrelSuperClusterIdMap")),
  endcapSuperClusterIdMapName_(Conf().getUntrackedParameter<string>
                               ("endcapClusterIdMapName","endcapSuperClusterIdMap")),
  trackMapName_(Conf().getUntrackedParameter<string>("trackMapName",
                                                     Form("%sMapName",mitName_.c_str()))),
  trackingMap_(0),
  tracks_(new mithep::TrackArr(250)), 
  trackMap_(new mithep::TrackMap)
{
  // Constructor.
  
  if (ecalAssocActive_) // initialize track associator configuration if needed
    assocParams_.loadParameters(
      cfg.getUntrackedParameter<ParameterSet>("TrackAssociatorParameters"));
}

//--------------------------------------------------------------------------------------------------
FillerTracks::~FillerTracks()
{
  // Destructor.

  delete tracks_;
  delete trackMap_;
}

//--------------------------------------------------------------------------------------------------
void FillerTracks::BookDataBlock(TreeWriter &tws)
{
  // Add tracks branch to tree, publish and get our objects.

  tws.AddBranch(mitName_,&tracks_);
  OS()->add<TrackArr>(tracks_,mitName_);

  if (!trackingMapName_.empty()) {
    trackingMap_ = OS()->get<TrackingParticleMap>(trackingMapName_);
    if (trackingMap_)
      AddBranchDep(mitName_,trackingMap_->GetBrName());
  }
  if (ecalAssocActive_ && !barrelSuperClusterIdMapName_.empty()) {
    barrelSuperClusterIdMap_ = OS()->get<SuperClusterIdMap>(barrelSuperClusterIdMapName_);
    if (barrelSuperClusterIdMap_)
      AddBranchDep(mitName_,barrelSuperClusterIdMap_->GetBrName());
  }
  if (ecalAssocActive_ && !endcapSuperClusterIdMapName_.empty()) {
    endcapSuperClusterIdMap_ = OS()->get<SuperClusterIdMap>(endcapSuperClusterIdMapName_);
    if (endcapSuperClusterIdMap_)
      AddBranchDep(mitName_,endcapSuperClusterIdMap_->GetBrName());
  }
  if (!trackMapName_.empty()) {
    trackMap_->SetBrName(mitName_);
    OS()->add<TrackMap>(trackMap_,trackMapName_);
  }
}

//--------------------------------------------------------------------------------------------------
void FillerTracks::FillDataBlock(const edm::Event      &event, 
                                 const edm::EventSetup &setup)
{
  // Fill tracks from edm collection into our collection.

  tracks_  ->Delete();
  trackMap_->Reset();

  // initialize handle and get product,
  // this usage allows also to get collections of classes which inherit from reco::Track
  Handle<View<reco::Track> > hTrackProduct;
  GetProduct(edmName_, hTrackProduct, event);  
        
  trackMap_->SetEdmProductId(hTrackProduct.id().id());
  const View<reco::Track> inTracks = *(hTrackProduct.product());  
  
  // for MC SimParticle association (reco->sim mappings)
  reco::RecoToSimCollection simAssociation;
  if (trackingMap_ && !edmSimAssocName_.empty()) {
    Handle<reco::RecoToSimCollection> simAssociationProduct;
    GetProduct(edmSimAssocName_, simAssociationProduct, event);  
    simAssociation = *(simAssociationProduct.product());
    if (verbose_>1)
      printf("SimAssociation Map Size = %i\n",simAssociation.size());
  }
  
  // set up associator for Track-Ecal associations
  TrackDetectorAssociator trackAssociator;
  trackAssociator.useDefaultPropagator();
  edm::ESHandle<MagneticField> bField;
  if (ecalAssocActive_)
    setup.get<IdealMagneticFieldRecord>().get(bField);
  
  // loop through all tracks and fill the information
  for (View<reco::Track>::const_iterator it = inTracks.begin();
         it != inTracks.end(); ++it) {
         
    mithep::Track *outTrack = tracks_->Allocate();
    // create track and set the core parameters
    new (outTrack) mithep::Track(it->qoverp(),it->lambda(),
                                 it->phi(),it->dxy(),it->dsz());
    outTrack->SetErrors(it->qoverpError(),it->lambdaError(),
                        it->phiError(),it->dxyError(),it->dszError());
        
    // fill track quality information
    outTrack->SetChi2(it->chi2());
    outTrack->SetNdof(static_cast<Int_t>(it->ndof()));
    
    //fill hit layer map
    const reco::HitPattern &hits = it->hitPattern();
    for (Int_t i=0; i<hits.numberOfHits(); i++) {
      uint32_t hit = hits.getHitPattern(i);
      if (hits.validHitFilter(hit))
        if (hits.trackerHitFilter(hit))
          outTrack->SetHit(hitReader_.Layer(hit));
                
      if (verbose_>2) {
        if (hits.muonDTHitFilter(hit))
          printf("Muon DT Layer = %i\n", hits.getLayer(hit));
        if (hits.muonCSCHitFilter(hit))
          printf("Muon CSC Layer = %i\n", hits.getLayer(hit));        
        if (hits.muonRPCHitFilter(hit))
          printf("Muon RPC Layer = %i\n", hits.getLayer(hit));
      }
    }

    //make ecal associations
    if (ecalAssocActive_) {
      TrackDetMatchInfo matchInfo;
      //Extra check to allow propagation to work in AOD where no TrackExtra is available
      if (it->extra().isAvailable())
        matchInfo = trackAssociator.associate(event,setup,*it,assocParams_);
      else {
        TrajectoryStateTransform tsTransform;
        FreeTrajectoryState initialState = tsTransform.initialFreeState(*it,&*bField);
        matchInfo = trackAssociator.associate(event, setup, assocParams_, &initialState);
      }
        
      outTrack->SetEtaEcal(matchInfo.trkGlobPosAtEcal.eta());
      outTrack->SetPhiEcal(matchInfo.trkGlobPosAtEcal.phi());

      //fill supercluster link
      if (barrelSuperClusterIdMap_ || endcapSuperClusterIdMap_) {
        mithep::SuperCluster *cluster = 0;
        for (std::vector<const ::CaloTower*>::const_iterator iTower =
               matchInfo.crossedTowers.begin(); 
             iTower<matchInfo.crossedTowers.end() && !cluster; iTower++) {

          for (uint ihit=0; ihit<(*iTower)->constituentsSize() && !cluster; ihit++) {
            DetId hit = (*iTower)->constituent(ihit);
            if (barrelSuperClusterIdMap_ && barrelSuperClusterIdMap_->HasMit(hit))
              cluster = barrelSuperClusterIdMap_->GetMit(hit);
            else if (endcapSuperClusterIdMap_ && endcapSuperClusterIdMap_->HasMit(hit))
              cluster = endcapSuperClusterIdMap_->GetMit(hit);
          }
        }
        if (cluster)
          outTrack->SetSCluster(cluster);
      }
    }
    
    // add reference between mithep and edm object
    if (trackMap_) {
      mitedm::TrackPtr thePtr = inTracks.ptrAt(it - inTracks.begin());
      trackMap_->Add(thePtr, outTrack);
    }

    //do dim associations
    if (trackingMap_ && !edmSimAssocName_.empty()) {
      if (verbose_>1)
        printf("Trying Track-Sim association\n");
      reco::TrackBaseRef theBaseRef = inTracks.refAt(it - inTracks.begin());
      vector<pair<TrackingParticleRef, double> > simRefs;
      Bool_t noSimParticle = 0;
      try {
        simRefs = simAssociation[theBaseRef]; //try to get the sim references if existing
      }
      catch (edm::Exception &ex) {
        noSimParticle = 1;
      }

      if (!noSimParticle) { //loop through sim match candidates
        if (verbose_>1)
          printf("Applying track-sim association\n");
        for (vector<pair<TrackingParticleRef, double> >::const_iterator simRefPair=simRefs.begin(); 
             simRefPair != simRefs.end(); ++simRefPair) 
          if (simRefPair->second > 0.5) // require more than 50% shared hits between reco and sim
            outTrack->SetMCPart(trackingMap_->GetMit(simRefPair->first)); //add reco->sim reference
      }
    }
  }
  tracks_->Trim();
}
Revision:	1.30
Committed:	Thu Mar 19 22:19:16 2009 UTC (16 years, 1 month ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_009a, Mit_009, Mit_008
Changes since 1.29:	+8 -7 lines
Log Message:	Add check for ecalAssocActive.
#	User	Rev	Content
1	loizides	1.30	// $Id: FillerTracks.cc,v 1.29 2009/03/15 11:20:41 loizides Exp $
2	loizides	1.1
3	bendavid	1.11	#include "MitProd/TreeFiller/interface/FillerTracks.h"
4	loizides	1.1	#include "FWCore/MessageLogger/interface/MessageLogger.h"
5			#include "DataFormats/Common/interface/Handle.h"
6	bendavid	1.14	#include "DataFormats/TrackReco/interface/HitPattern.h"
7	loizides	1.1	#include "DataFormats/TrackReco/interface/Track.h"
8			#include "DataFormats/TrackReco/interface/TrackFwd.h"
9			#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
10			#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
11			#include "DataFormats/RecoCandidate/interface/TrackAssociation.h"
12	bendavid	1.23	#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
13	bendavid	1.28	#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
14			#include "MagneticField/Engine/interface/MagneticField.h"
15			#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
16	bendavid	1.26	#include "MitEdm/DataFormats/interface/Types.h"
17	loizides	1.1
18			using namespace std;
19			using namespace edm;
20			using namespace mithep;
21
22	loizides	1.4	//--------------------------------------------------------------------------------------------------
23	loizides	1.29	FillerTracks::FillerTracks(const ParameterSet &cfg, const char *name, bool active) :
24	loizides	1.7	BaseFiller(cfg,name,active),
25	loizides	1.29	ecalAssocActive_(Conf().getUntrackedParameter<bool>("ecalAssocActive",0)),
26			edmName_(Conf().getUntrackedParameter<string>("edmName","generalTracks")),
27			mitName_(Conf().getUntrackedParameter<string>("mitName",Names::gkTrackBrn)),
28	loizides	1.19	edmSimAssocName_(Conf().getUntrackedParameter<string>("edmSimAssociationName","")),
29	bendavid	1.20	trackingMapName_(Conf().getUntrackedParameter<string>("trackingMapName","TrackingMap")),
30	loizides	1.29	barrelSuperClusterIdMapName_(Conf().getUntrackedParameter<string>
31			("superClusterIdMapName","barrelSuperClusterIdMap")),
32			endcapSuperClusterIdMapName_(Conf().getUntrackedParameter<string>
33			("endcapClusterIdMapName","endcapSuperClusterIdMap")),
34	loizides	1.13	trackMapName_(Conf().getUntrackedParameter<string>("trackMapName",
35			Form("%sMapName",mitName_.c_str()))),
36	bendavid	1.17	trackingMap_(0),
37	loizides	1.7	tracks_(new mithep::TrackArr(250)),
38	loizides	1.1	trackMap_(new mithep::TrackMap)
39			{
40	loizides	1.13	// Constructor.
41	bendavid	1.23
42	loizides	1.29	if (ecalAssocActive_) // initialize track associator configuration if needed
43			assocParams_.loadParameters(
44			cfg.getUntrackedParameter<ParameterSet>("TrackAssociatorParameters"));
45	loizides	1.1	}
46
47	loizides	1.4	//--------------------------------------------------------------------------------------------------
48	loizides	1.1	FillerTracks::~FillerTracks()
49			{
50	loizides	1.13	// Destructor.
51
52	loizides	1.5	delete tracks_;
53	loizides	1.1	delete trackMap_;
54			}
55
56	loizides	1.4	//--------------------------------------------------------------------------------------------------
57	loizides	1.1	void FillerTracks::BookDataBlock(TreeWriter &tws)
58			{
59	loizides	1.13	// Add tracks branch to tree, publish and get our objects.
60
61	loizides	1.29	tws.AddBranch(mitName_,&tracks_);
62			OS()->add<TrackArr>(tracks_,mitName_);
63	loizides	1.12
64	loizides	1.29	if (!trackingMapName_.empty()) {
65			trackingMap_ = OS()->get<TrackingParticleMap>(trackingMapName_);
66			if (trackingMap_)
67			AddBranchDep(mitName_,trackingMap_->GetBrName());
68			}
69	loizides	1.30	if (ecalAssocActive_ && !barrelSuperClusterIdMapName_.empty()) {
70	loizides	1.29	barrelSuperClusterIdMap_ = OS()->get<SuperClusterIdMap>(barrelSuperClusterIdMapName_);
71			if (barrelSuperClusterIdMap_)
72			AddBranchDep(mitName_,barrelSuperClusterIdMap_->GetBrName());
73			}
74	loizides	1.30	if (ecalAssocActive_ && !endcapSuperClusterIdMapName_.empty()) {
75	loizides	1.29	endcapSuperClusterIdMap_ = OS()->get<SuperClusterIdMap>(endcapSuperClusterIdMapName_);
76			if (endcapSuperClusterIdMap_)
77			AddBranchDep(mitName_,endcapSuperClusterIdMap_->GetBrName());
78			}
79			if (!trackMapName_.empty()) {
80			trackMap_->SetBrName(mitName_);
81			OS()->add<TrackMap>(trackMap_,trackMapName_);
82			}
83	loizides	1.1	}
84
85	loizides	1.4	//--------------------------------------------------------------------------------------------------
86	loizides	1.1	void FillerTracks::FillDataBlock(const edm::Event &event,
87			const edm::EventSetup &setup)
88			{
89	loizides	1.13	// Fill tracks from edm collection into our collection.
90
91	bendavid	1.27	tracks_ ->Delete();
92	loizides	1.1	trackMap_->Reset();
93	loizides	1.6
94	loizides	1.29	// initialize handle and get product,
95			// this usage allows also to get collections of classes which inherit from reco::Track
96	bendavid	1.24	Handle<View<reco::Track> > hTrackProduct;
97	loizides	1.6	GetProduct(edmName_, hTrackProduct, event);
98	loizides	1.1
99	loizides	1.6	trackMap_->SetEdmProductId(hTrackProduct.id().id());
100	bendavid	1.24	const View<reco::Track> inTracks = *(hTrackProduct.product());
101	loizides	1.1
102	paus	1.10	// for MC SimParticle association (reco->sim mappings)
103	loizides	1.1	reco::RecoToSimCollection simAssociation;
104	loizides	1.19	if (trackingMap_ && !edmSimAssocName_.empty()) {
105	loizides	1.1	Handle<reco::RecoToSimCollection> simAssociationProduct;
106	loizides	1.19	GetProduct(edmSimAssocName_, simAssociationProduct, event);
107	loizides	1.1	simAssociation = *(simAssociationProduct.product());
108	loizides	1.29	if (verbose_>1)
109			printf("SimAssociation Map Size = %i\n",simAssociation.size());
110	loizides	1.1	}
111
112	loizides	1.29	// set up associator for Track-Ecal associations
113	bendavid	1.23	TrackDetectorAssociator trackAssociator;
114			trackAssociator.useDefaultPropagator();
115	bendavid	1.28	edm::ESHandle<MagneticField> bField;
116			if (ecalAssocActive_)
117			setup.get<IdealMagneticFieldRecord>().get(bField);
118
119	paus	1.10	// loop through all tracks and fill the information
120	bendavid	1.24	for (View<reco::Track>::const_iterator it = inTracks.begin();
121			it != inTracks.end(); ++it) {
122
123	loizides	1.5	mithep::Track *outTrack = tracks_->Allocate();
124	paus	1.10	// create track and set the core parameters
125	loizides	1.19	new (outTrack) mithep::Track(it->qoverp(),it->lambda(),
126			it->phi(),it->dxy(),it->dsz());
127			outTrack->SetErrors(it->qoverpError(),it->lambdaError(),
128			it->phiError(),it->dxyError(),it->dszError());
129	loizides	1.1
130	loizides	1.19	// fill track quality information
131	bendavid	1.14	outTrack->SetChi2(it->chi2());
132	loizides	1.16	outTrack->SetNdof(static_cast<Int_t>(it->ndof()));
133	bendavid	1.14
134	loizides	1.19	//fill hit layer map
135	bendavid	1.14	const reco::HitPattern &hits = it->hitPattern();
136			for (Int_t i=0; i<hits.numberOfHits(); i++) {
137			uint32_t hit = hits.getHitPattern(i);
138	loizides	1.19	if (hits.validHitFilter(hit))
139			if (hits.trackerHitFilter(hit))
140	bendavid	1.22	outTrack->SetHit(hitReader_.Layer(hit));
141	bendavid	1.14
142	loizides	1.29	if (verbose_>2) {
143	loizides	1.19	if (hits.muonDTHitFilter(hit))
144			printf("Muon DT Layer = %i\n", hits.getLayer(hit));
145			if (hits.muonCSCHitFilter(hit))
146			printf("Muon CSC Layer = %i\n", hits.getLayer(hit));
147			if (hits.muonRPCHitFilter(hit))
148			printf("Muon RPC Layer = %i\n", hits.getLayer(hit));
149			}
150	bendavid	1.14	}
151	bendavid	1.23
152			//make ecal associations
153			if (ecalAssocActive_) {
154	bendavid	1.28	TrackDetMatchInfo matchInfo;
155			//Extra check to allow propagation to work in AOD where no TrackExtra is available
156			if (it->extra().isAvailable())
157			matchInfo = trackAssociator.associate(event,setup,*it,assocParams_);
158			else {
159			TrajectoryStateTransform tsTransform;
160			FreeTrajectoryState initialState = tsTransform.initialFreeState(it,&bField);
161			matchInfo = trackAssociator.associate(event, setup, assocParams_, &initialState);
162			}
163
164	bendavid	1.23	outTrack->SetEtaEcal(matchInfo.trkGlobPosAtEcal.eta());
165			outTrack->SetPhiEcal(matchInfo.trkGlobPosAtEcal.phi());
166
167			//fill supercluster link
168			if (barrelSuperClusterIdMap_ \|\| endcapSuperClusterIdMap_) {
169			mithep::SuperCluster *cluster = 0;
170	loizides	1.29	for (std::vector<const ::CaloTower*>::const_iterator iTower =
171			matchInfo.crossedTowers.begin();
172			iTower<matchInfo.crossedTowers.end() && !cluster; iTower++) {
173	bendavid	1.23
174			for (uint ihit=0; ihit<(*iTower)->constituentsSize() && !cluster; ihit++) {
175			DetId hit = (*iTower)->constituent(ihit);
176			if (barrelSuperClusterIdMap_ && barrelSuperClusterIdMap_->HasMit(hit))
177			cluster = barrelSuperClusterIdMap_->GetMit(hit);
178			else if (endcapSuperClusterIdMap_ && endcapSuperClusterIdMap_->HasMit(hit))
179			cluster = endcapSuperClusterIdMap_->GetMit(hit);
180			}
181			}
182			if (cluster)
183			outTrack->SetSCluster(cluster);
184			}
185			}
186	bendavid	1.14
187	loizides	1.1	// add reference between mithep and edm object
188	loizides	1.30	if (trackMap_) {
189			mitedm::TrackPtr thePtr = inTracks.ptrAt(it - inTracks.begin());
190			trackMap_->Add(thePtr, outTrack);
191			}
192
193	bendavid	1.23	//do dim associations
194	loizides	1.19	if (trackingMap_ && !edmSimAssocName_.empty()) {
195	loizides	1.29	if (verbose_>1)
196			printf("Trying Track-Sim association\n");
197	bendavid	1.24	reco::TrackBaseRef theBaseRef = inTracks.refAt(it - inTracks.begin());
198	loizides	1.1	vector<pair<TrackingParticleRef, double> > simRefs;
199	paus	1.10	Bool_t noSimParticle = 0;
200	loizides	1.1	try {
201			simRefs = simAssociation[theBaseRef]; //try to get the sim references if existing
202			}
203	loizides	1.5	catch (edm::Exception &ex) {
204	paus	1.10	noSimParticle = 1;
205	loizides	1.1	}
206	paus	1.10
207	loizides	1.1	if (!noSimParticle) { //loop through sim match candidates
208	loizides	1.29	if (verbose_>1)
209			printf("Applying track-sim association\n");
210	loizides	1.1	for (vector<pair<TrackingParticleRef, double> >::const_iterator simRefPair=simRefs.begin();
211			simRefPair != simRefs.end(); ++simRefPair)
212	paus	1.10	if (simRefPair->second > 0.5) // require more than 50% shared hits between reco and sim
213	bendavid	1.17	outTrack->SetMCPart(trackingMap_->GetMit(simRefPair->first)); //add reco->sim reference
214	loizides	1.1	}
215			}
216			}
217			tracks_->Trim();
218			}