TreeFiller/src/FillerTracks.cc

// $Id: FillerTracks.cc,v 1.24 2008/11/03 18:11:10 bendavid 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"

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

//--------------------------------------------------------------------------------------------------
FillerTracks::FillerTracks(const ParameterSet &cfg, const char *name, bool active, bool ecalActive,
                           const char *edmName, const char *mitName) :
  BaseFiller(cfg,name,active),
  ecalAssocActive_(Conf().getUntrackedParameter<bool>("ecalAssocActive",ecalActive)),
  edmName_(Conf().getUntrackedParameter<string>("edmName",edmName)),
  mitName_(Conf().getUntrackedParameter<string>("mitName",mitName)),
  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.
  
  //initialize track associator configuration if needed
  if (ecalAssocActive_)
    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_.c_str(),&tracks_);

  trackingMap_ = OS()->get<TrackingParticleMap>(trackingMapName_.c_str());
  barrelSuperClusterIdMap_ = OS()->get<SuperClusterIdMap>(barrelSuperClusterIdMapName_.c_str());
  endcapSuperClusterIdMap_ = OS()->get<SuperClusterIdMap>(endcapSuperClusterIdMapName_.c_str());
  OS()->add<TrackMap>(trackMap_,trackMapName_.c_str());
  OS()->add<TrackArr>(tracks_,mitName_.c_str());
}

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

  tracks_  ->Reset();
  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());
    //printf("SimAssociation Map Size = %i\n",simAssociation.size());
  }
  
  //set up associator for Track-Ecal associations
  TrackDetectorAssociator trackAssociator;
  trackAssociator.useDefaultPropagator();

  // 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 (0) {
        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 = trackAssociator.associate(event,setup,
                                                             *it,
                                                              assocParams_);
      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
    edm::Ptr<reco::Track> thePtr = inTracks.ptrAt(it - inTracks.begin());
    trackMap_->Add(thePtr, outTrack);
        
    //do dim associations
    if (trackingMap_ && !edmSimAssocName_.empty()) {
      //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
        //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.25
Committed:	Mon Nov 3 18:16:23 2008 UTC (16 years, 6 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
Changes since 1.24:	+3 -1 lines
Log Message:	comments
#	User	Rev	Content
1	bendavid	1.25	// $Id: FillerTracks.cc,v 1.24 2008/11/03 18:11:10 bendavid 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	loizides	1.1
14			using namespace std;
15			using namespace edm;
16			using namespace mithep;
17
18	loizides	1.4	//--------------------------------------------------------------------------------------------------
19	bendavid	1.23	FillerTracks::FillerTracks(const ParameterSet &cfg, const char *name, bool active, bool ecalActive,
20	loizides	1.16	const char edmName, const char mitName) :
21	loizides	1.7	BaseFiller(cfg,name,active),
22	bendavid	1.23	ecalAssocActive_(Conf().getUntrackedParameter<bool>("ecalAssocActive",ecalActive)),
23	loizides	1.16	edmName_(Conf().getUntrackedParameter<string>("edmName",edmName)),
24			mitName_(Conf().getUntrackedParameter<string>("mitName",mitName)),
25	loizides	1.19	edmSimAssocName_(Conf().getUntrackedParameter<string>("edmSimAssociationName","")),
26	bendavid	1.20	trackingMapName_(Conf().getUntrackedParameter<string>("trackingMapName","TrackingMap")),
27	bendavid	1.23	barrelSuperClusterIdMapName_(Conf().getUntrackedParameter<string>("superClusterIdMapName","barrelSuperClusterIdMap")),
28			endcapSuperClusterIdMapName_(Conf().getUntrackedParameter<string>("endcapClusterIdMapName","endcapSuperClusterIdMap")),
29	loizides	1.13	trackMapName_(Conf().getUntrackedParameter<string>("trackMapName",
30			Form("%sMapName",mitName_.c_str()))),
31	bendavid	1.17	trackingMap_(0),
32	loizides	1.7	tracks_(new mithep::TrackArr(250)),
33	loizides	1.1	trackMap_(new mithep::TrackMap)
34			{
35	loizides	1.13	// Constructor.
36	bendavid	1.23
37			//initialize track associator configuration if needed
38			if (ecalAssocActive_)
39			assocParams_.loadParameters(cfg.getUntrackedParameter<ParameterSet>("TrackAssociatorParameters"));
40
41	loizides	1.1	}
42
43	loizides	1.4	//--------------------------------------------------------------------------------------------------
44	loizides	1.1	FillerTracks::~FillerTracks()
45			{
46	loizides	1.13	// Destructor.
47
48	loizides	1.5	delete tracks_;
49	loizides	1.1	delete trackMap_;
50			}
51
52	loizides	1.4	//--------------------------------------------------------------------------------------------------
53	loizides	1.1	void FillerTracks::BookDataBlock(TreeWriter &tws)
54			{
55	loizides	1.13	// Add tracks branch to tree, publish and get our objects.
56
57	loizides	1.1	tws.AddBranch(mitName_.c_str(),&tracks_);
58	loizides	1.12
59	bendavid	1.17	trackingMap_ = OS()->get<TrackingParticleMap>(trackingMapName_.c_str());
60	bendavid	1.23	barrelSuperClusterIdMap_ = OS()->get<SuperClusterIdMap>(barrelSuperClusterIdMapName_.c_str());
61			endcapSuperClusterIdMap_ = OS()->get<SuperClusterIdMap>(endcapSuperClusterIdMapName_.c_str());
62	loizides	1.13	OS()->add<TrackMap>(trackMap_,trackMapName_.c_str());
63			OS()->add<TrackArr>(tracks_,mitName_.c_str());
64	loizides	1.1	}
65
66	loizides	1.4	//--------------------------------------------------------------------------------------------------
67	loizides	1.1	void FillerTracks::FillDataBlock(const edm::Event &event,
68			const edm::EventSetup &setup)
69			{
70	loizides	1.13	// Fill tracks from edm collection into our collection.
71
72	paus	1.10	tracks_ ->Reset();
73	loizides	1.1	trackMap_->Reset();
74	loizides	1.6
75	bendavid	1.25	//initialize Handle and get product, this usage allows also to get collections of classes which
76			//inherit from reco::Track
77	bendavid	1.24	Handle<View<reco::Track> > hTrackProduct;
78	loizides	1.6	GetProduct(edmName_, hTrackProduct, event);
79	loizides	1.1
80	loizides	1.6	trackMap_->SetEdmProductId(hTrackProduct.id().id());
81	bendavid	1.24	const View<reco::Track> inTracks = *(hTrackProduct.product());
82	loizides	1.1
83	paus	1.10	// for MC SimParticle association (reco->sim mappings)
84	loizides	1.1	reco::RecoToSimCollection simAssociation;
85	loizides	1.19	if (trackingMap_ && !edmSimAssocName_.empty()) {
86	loizides	1.1	Handle<reco::RecoToSimCollection> simAssociationProduct;
87	loizides	1.19	GetProduct(edmSimAssocName_, simAssociationProduct, event);
88	loizides	1.1	simAssociation = *(simAssociationProduct.product());
89	bendavid	1.17	//printf("SimAssociation Map Size = %i\n",simAssociation.size());
90	loizides	1.1	}
91
92	bendavid	1.23	//set up associator for Track-Ecal associations
93			TrackDetectorAssociator trackAssociator;
94			trackAssociator.useDefaultPropagator();
95
96	paus	1.10	// loop through all tracks and fill the information
97	bendavid	1.24	for (View<reco::Track>::const_iterator it = inTracks.begin();
98			it != inTracks.end(); ++it) {
99
100	loizides	1.5	mithep::Track *outTrack = tracks_->Allocate();
101	paus	1.10	// create track and set the core parameters
102	loizides	1.19	new (outTrack) mithep::Track(it->qoverp(),it->lambda(),
103			it->phi(),it->dxy(),it->dsz());
104			outTrack->SetErrors(it->qoverpError(),it->lambdaError(),
105			it->phiError(),it->dxyError(),it->dszError());
106	loizides	1.1
107	loizides	1.19	// fill track quality information
108	bendavid	1.14	outTrack->SetChi2(it->chi2());
109	loizides	1.16	outTrack->SetNdof(static_cast<Int_t>(it->ndof()));
110	bendavid	1.14
111	loizides	1.19	//fill hit layer map
112	bendavid	1.14	const reco::HitPattern &hits = it->hitPattern();
113			for (Int_t i=0; i<hits.numberOfHits(); i++) {
114			uint32_t hit = hits.getHitPattern(i);
115	loizides	1.19	if (hits.validHitFilter(hit))
116			if (hits.trackerHitFilter(hit))
117	bendavid	1.22	outTrack->SetHit(hitReader_.Layer(hit));
118	bendavid	1.14
119	loizides	1.19	if (0) {
120			if (hits.muonDTHitFilter(hit))
121			printf("Muon DT Layer = %i\n", hits.getLayer(hit));
122			if (hits.muonCSCHitFilter(hit))
123			printf("Muon CSC Layer = %i\n", hits.getLayer(hit));
124			if (hits.muonRPCHitFilter(hit))
125			printf("Muon RPC Layer = %i\n", hits.getLayer(hit));
126			}
127	bendavid	1.14	}
128	bendavid	1.23
129			//make ecal associations
130			if (ecalAssocActive_) {
131			TrackDetMatchInfo matchInfo = trackAssociator.associate(event,setup,
132			*it,
133			assocParams_);
134			outTrack->SetEtaEcal(matchInfo.trkGlobPosAtEcal.eta());
135			outTrack->SetPhiEcal(matchInfo.trkGlobPosAtEcal.phi());
136
137			//fill supercluster link
138			if (barrelSuperClusterIdMap_ \|\| endcapSuperClusterIdMap_) {
139			mithep::SuperCluster *cluster = 0;
140			for (std::vector<const ::CaloTower*>::const_iterator iTower = matchInfo.crossedTowers.begin();
141			iTower<matchInfo.crossedTowers.end() && !cluster; iTower++) {
142
143			for (uint ihit=0; ihit<(*iTower)->constituentsSize() && !cluster; ihit++) {
144			DetId hit = (*iTower)->constituent(ihit);
145			if (barrelSuperClusterIdMap_ && barrelSuperClusterIdMap_->HasMit(hit))
146			cluster = barrelSuperClusterIdMap_->GetMit(hit);
147			else if (endcapSuperClusterIdMap_ && endcapSuperClusterIdMap_->HasMit(hit))
148			cluster = endcapSuperClusterIdMap_->GetMit(hit);
149			}
150			}
151			if (cluster)
152			outTrack->SetSCluster(cluster);
153
154			}
155			}
156	bendavid	1.14
157	loizides	1.1	// add reference between mithep and edm object
158	bendavid	1.24	edm::Ptr<reco::Track> thePtr = inTracks.ptrAt(it - inTracks.begin());
159			trackMap_->Add(thePtr, outTrack);
160	loizides	1.1
161	bendavid	1.23	//do dim associations
162	loizides	1.19	if (trackingMap_ && !edmSimAssocName_.empty()) {
163	bendavid	1.17	//printf("Trying Track-Sim association\n");
164	bendavid	1.24	reco::TrackBaseRef theBaseRef = inTracks.refAt(it - inTracks.begin());
165	loizides	1.1	vector<pair<TrackingParticleRef, double> > simRefs;
166	paus	1.10	Bool_t noSimParticle = 0;
167	loizides	1.1	try {
168			simRefs = simAssociation[theBaseRef]; //try to get the sim references if existing
169			}
170	loizides	1.5	catch (edm::Exception &ex) {
171	paus	1.10	noSimParticle = 1;
172	loizides	1.1	}
173	paus	1.10
174	loizides	1.1	if (!noSimParticle) { //loop through sim match candidates
175	bendavid	1.17	//printf("Applying track-sim association\n");
176	loizides	1.1	for (vector<pair<TrackingParticleRef, double> >::const_iterator simRefPair=simRefs.begin();
177			simRefPair != simRefs.end(); ++simRefPair)
178
179	paus	1.10	if (simRefPair->second > 0.5) // require more than 50% shared hits between reco and sim
180	bendavid	1.17	outTrack->SetMCPart(trackingMap_->GetMit(simRefPair->first)); //add reco->sim reference
181	loizides	1.1	}
182			}
183			}
184			tracks_->Trim();
185			}
186	bendavid	1.14