TreeFiller/src/FillerTracks.cc

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