TreeFiller/src/FillerTracks.cc

// $Id: FillerTracks.cc,v 1.41 2011/03/11 04:03:55 bendavid Exp $

#include "MitProd/TreeFiller/interface/FillerTracks.h"
#include "DataFormats/RecoCandidate/interface/TrackAssociation.h"
#include "DataFormats/TrackReco/interface/HitPattern.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "MagneticField/Engine/interface/MagneticField.h"
#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
#include "TrackingTools/TrackAssociator/interface/TrackDetectorAssociator.h"
#include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
#include "MitAna/DataTree/interface/Names.h"
#include "MitAna/DataTree/interface/TrackCol.h"
#include "MitProd/ObjectService/interface/ObjectService.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);  
        
  //printf("edmName = %s, product id = %i\n",edmName_.c_str(),hTrackProduct.id().id());
  
  trackMap_->SetEdmProductId(hTrackProduct.id().id());
  const View<reco::Track> inTracks = *(hTrackProduct.product());  
  
  if (verbose_>1)
    printf("Track Collection: %s, id=%i\n",edmName_.c_str(),hTrackProduct.id().id());
  
  // 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",(int)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());
                        
    outTrack->SetPtErr(it->ptError());
        
    // fill track quality information
    outTrack->SetChi2(it->chi2());
    outTrack->SetNdof(static_cast<Int_t>(it->ndof()));

    //fill algo enum
    outTrack->SetAlgo(static_cast<mithep::Track::ETrackAlgorithm>(it->algo()));

    //fill quality bitmask
    outTrack->Quality().SetQualityMask(BitMask8(UChar_t(it->qualityMask())));
    
    if (verbose_>2) {
      printf("reco::Track   high purity = %i,  ", it->quality(reco::TrackBase::highPurity));
      printf("mithep::Track high purity = %i\n",outTrack->Quality().Quality(mithep::TrackQuality::highPurity));
    }
    
    //fill gsf flag, some type gymastics needed...
    if (typeid(*it)==typeid(reco::GsfTrack))
      outTrack->SetIsGsf(kTRUE);
    else
      outTrack->SetIsGsf(kFALSE);
    
    if (verbose_>2) {
      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());
    }

    //fill hit layer map and missing hits
    const reco::HitPattern &hits = it->hitPattern();
    BitMask48 hitLayers;
    BitMask48 missingHitLayers;
    for (Int_t i=0; i<hits.numberOfHits(); i++) {
      uint32_t hit = hits.getHitPattern(i);
      if (hits.validHitFilter(hit)) {
        if (hits.trackerHitFilter(hit)) {
          hitLayers.SetBit(hitReader_.Layer(hit));
          if (verbose_>2)
            printf("Found valid hit with layer %i\n",hitReader_.Layer(hit));
        }
      }

      if (hits.getHitType(hit)==1) {
        if (hits.trackerHitFilter(hit)) {
          missingHitLayers.SetBit(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));
      }
    }

    outTrack->SetHits(hitLayers);
    outTrack->SetMissingHits(missingHitLayers);


    //set expected inner hits
    const reco::HitPattern &expectedInnerHitPattern = it->trackerExpectedHitsInner();
    BitMask48 expectedHitsInner;
    // search for all good crossed layers (with or without hit)
    for (Int_t hi=0; hi<expectedInnerHitPattern.numberOfHits(); hi++) {
      uint32_t hit = expectedInnerHitPattern.getHitPattern(hi);
      if (expectedInnerHitPattern.getHitType(hit)<=1) {
        if (expectedInnerHitPattern.trackerHitFilter(hit)) {
          expectedHitsInner.SetBit(hitReader_.Layer(hit));
        }
      }
    }

    outTrack->SetExpectedHitsInner(expectedHitsInner);

    //set expected outer hits
    const reco::HitPattern &expectedOuterHitPattern = it->trackerExpectedHitsOuter();
    BitMask48 expectedHitsOuter;
    // search for all good crossed layers (with or without hit)
    for (Int_t hi=0; hi<expectedOuterHitPattern.numberOfHits(); hi++) {
      uint32_t hit = expectedOuterHitPattern.getHitPattern(hi);
      if (expectedOuterHitPattern.getHitType(hit)<=1) {
        if (expectedOuterHitPattern.trackerHitFilter(hit)) {
          expectedHitsOuter.SetBit(hitReader_.Layer(hit));
        }
      }
    }

    outTrack->SetExpectedHitsOuter(expectedHitsOuter);

    //fill actual hit counts
    outTrack->SetNHits(it->numberOfValidHits());
    outTrack->SetNPixelHits(it->hitPattern().numberOfValidPixelHits());
    outTrack->SetNMissingHits(it->hitPattern().numberOfLostHits());
    outTrack->SetNExpectedHitsInner(it->trackerExpectedHitsInner().numberOfHits());
    outTrack->SetNExpectedHitsOuter(it->trackerExpectedHitsOuter().numberOfHits());


    if (verbose_>2)
      printf("OutTrack NHits = %i, NMissingHits = %i, NExpectedHitsInner = %i, NExpectedHitsOuter = %i\n",outTrack->NHits(),outTrack->NMissingHits(),outTrack->NExpectedHitsInner(),outTrack->NExpectedHitsOuter()); 


    //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);
      }
       
      if (matchInfo.isGoodEcal) {
        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 sim 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.42
Committed:	Wed Sep 28 16:50:07 2011 UTC (13 years, 7 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_025c_branch2, Mit_025c_branch1, Mit_025c_branch0, Mit_025e, Mit_025d, Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2
Branch point for:	Mit_025c_branch
Changes since 1.41:	+3 -1 lines
Log Message:	changes for 440 event content
#	User	Rev	Content
1	bendavid	1.42	// $Id: FillerTracks.cc,v 1.41 2011/03/11 04:03:55 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	bendavid	1.42	//printf("edmName = %s, product id = %i\n",edmName_.c_str(),hTrackProduct.id().id());
102
103	loizides	1.6	trackMap_->SetEdmProductId(hTrackProduct.id().id());
104	bendavid	1.24	const View<reco::Track> inTracks = *(hTrackProduct.product());
105	loizides	1.1
106	bendavid	1.32	if (verbose_>1)
107			printf("Track Collection: %s, id=%i\n",edmName_.c_str(),hTrackProduct.id().id());
108
109	paus	1.10	// for MC SimParticle association (reco->sim mappings)
110	loizides	1.1	reco::RecoToSimCollection simAssociation;
111	loizides	1.19	if (trackingMap_ && !edmSimAssocName_.empty()) {
112	loizides	1.1	Handle<reco::RecoToSimCollection> simAssociationProduct;
113	loizides	1.19	GetProduct(edmSimAssocName_, simAssociationProduct, event);
114	loizides	1.1	simAssociation = *(simAssociationProduct.product());
115	loizides	1.29	if (verbose_>1)
116	bendavid	1.41	printf("SimAssociation Map Size = %i\n",(int)simAssociation.size());
117	loizides	1.1	}
118
119	loizides	1.29	// set up associator for Track-Ecal associations
120	bendavid	1.23	TrackDetectorAssociator trackAssociator;
121			trackAssociator.useDefaultPropagator();
122	bendavid	1.28	edm::ESHandle<MagneticField> bField;
123	bendavid	1.36	if (ecalAssocActive_) {
124	bendavid	1.28	setup.get<IdealMagneticFieldRecord>().get(bField);
125	bendavid	1.36	}
126
127
128	paus	1.10	// loop through all tracks and fill the information
129	bendavid	1.24	for (View<reco::Track>::const_iterator it = inTracks.begin();
130			it != inTracks.end(); ++it) {
131
132	loizides	1.5	mithep::Track *outTrack = tracks_->Allocate();
133	paus	1.10	// create track and set the core parameters
134	loizides	1.19	new (outTrack) mithep::Track(it->qoverp(),it->lambda(),
135			it->phi(),it->dxy(),it->dsz());
136			outTrack->SetErrors(it->qoverpError(),it->lambdaError(),
137			it->phiError(),it->dxyError(),it->dszError());
138	bendavid	1.39
139			outTrack->SetPtErr(it->ptError());
140	loizides	1.1
141	loizides	1.19	// fill track quality information
142	bendavid	1.14	outTrack->SetChi2(it->chi2());
143	loizides	1.16	outTrack->SetNdof(static_cast<Int_t>(it->ndof()));
144	bendavid	1.33
145			//fill algo enum
146			outTrack->SetAlgo(static_cast<mithep::Track::ETrackAlgorithm>(it->algo()));
147
148	bendavid	1.35	//fill quality bitmask
149			outTrack->Quality().SetQualityMask(BitMask8(UChar_t(it->qualityMask())));
150
151			if (verbose_>2) {
152			printf("reco::Track high purity = %i, ", it->quality(reco::TrackBase::highPurity));
153			printf("mithep::Track high purity = %i\n",outTrack->Quality().Quality(mithep::TrackQuality::highPurity));
154			}
155
156	bendavid	1.33	//fill gsf flag, some type gymastics needed...
157			if (typeid(*it)==typeid(reco::GsfTrack))
158			outTrack->SetIsGsf(kTRUE);
159			else
160			outTrack->SetIsGsf(kFALSE);
161	bendavid	1.14
162	bendavid	1.40	if (verbose_>2) {
163			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());
164			}
165
166	bendavid	1.38	//fill hit layer map and missing hits
167	bendavid	1.14	const reco::HitPattern &hits = it->hitPattern();
168	bendavid	1.38	BitMask48 hitLayers;
169			BitMask48 missingHitLayers;
170	bendavid	1.14	for (Int_t i=0; i<hits.numberOfHits(); i++) {
171			uint32_t hit = hits.getHitPattern(i);
172	bendavid	1.38	if (hits.validHitFilter(hit)) {
173			if (hits.trackerHitFilter(hit)) {
174			hitLayers.SetBit(hitReader_.Layer(hit));
175	bendavid	1.40	if (verbose_>2)
176			printf("Found valid hit with layer %i\n",hitReader_.Layer(hit));
177	bendavid	1.38	}
178			}
179
180			if (hits.getHitType(hit)==1) {
181			if (hits.trackerHitFilter(hit)) {
182			missingHitLayers.SetBit(hitReader_.Layer(hit));
183			}
184			}
185
186	loizides	1.29	if (verbose_>2) {
187	loizides	1.19	if (hits.muonDTHitFilter(hit))
188			printf("Muon DT Layer = %i\n", hits.getLayer(hit));
189			if (hits.muonCSCHitFilter(hit))
190			printf("Muon CSC Layer = %i\n", hits.getLayer(hit));
191			if (hits.muonRPCHitFilter(hit))
192			printf("Muon RPC Layer = %i\n", hits.getLayer(hit));
193			}
194	bendavid	1.14	}
195	bendavid	1.23
196	bendavid	1.38	outTrack->SetHits(hitLayers);
197			outTrack->SetMissingHits(missingHitLayers);
198
199	bendavid	1.40
200	bendavid	1.38	//set expected inner hits
201			const reco::HitPattern &expectedInnerHitPattern = it->trackerExpectedHitsInner();
202			BitMask48 expectedHitsInner;
203			// search for all good crossed layers (with or without hit)
204			for (Int_t hi=0; hi<expectedInnerHitPattern.numberOfHits(); hi++) {
205			uint32_t hit = expectedInnerHitPattern.getHitPattern(hi);
206			if (expectedInnerHitPattern.getHitType(hit)<=1) {
207			if (expectedInnerHitPattern.trackerHitFilter(hit)) {
208			expectedHitsInner.SetBit(hitReader_.Layer(hit));
209			}
210			}
211			}
212
213			outTrack->SetExpectedHitsInner(expectedHitsInner);
214
215			//set expected outer hits
216			const reco::HitPattern &expectedOuterHitPattern = it->trackerExpectedHitsOuter();
217			BitMask48 expectedHitsOuter;
218			// search for all good crossed layers (with or without hit)
219			for (Int_t hi=0; hi<expectedOuterHitPattern.numberOfHits(); hi++) {
220			uint32_t hit = expectedOuterHitPattern.getHitPattern(hi);
221			if (expectedOuterHitPattern.getHitType(hit)<=1) {
222			if (expectedOuterHitPattern.trackerHitFilter(hit)) {
223			expectedHitsOuter.SetBit(hitReader_.Layer(hit));
224			}
225			}
226			}
227
228			outTrack->SetExpectedHitsOuter(expectedHitsOuter);
229
230	bendavid	1.40	//fill actual hit counts
231			outTrack->SetNHits(it->numberOfValidHits());
232			outTrack->SetNPixelHits(it->hitPattern().numberOfValidPixelHits());
233			outTrack->SetNMissingHits(it->hitPattern().numberOfLostHits());
234			outTrack->SetNExpectedHitsInner(it->trackerExpectedHitsInner().numberOfHits());
235			outTrack->SetNExpectedHitsOuter(it->trackerExpectedHitsOuter().numberOfHits());
236
237
238			if (verbose_>2)
239			printf("OutTrack NHits = %i, NMissingHits = %i, NExpectedHitsInner = %i, NExpectedHitsOuter = %i\n",outTrack->NHits(),outTrack->NMissingHits(),outTrack->NExpectedHitsInner(),outTrack->NExpectedHitsOuter());
240
241
242	bendavid	1.23	//make ecal associations
243			if (ecalAssocActive_) {
244	bendavid	1.28	TrackDetMatchInfo matchInfo;
245			//Extra check to allow propagation to work in AOD where no TrackExtra is available
246			if (it->extra().isAvailable())
247			matchInfo = trackAssociator.associate(event,setup,*it,assocParams_);
248			else {
249			TrajectoryStateTransform tsTransform;
250			FreeTrajectoryState initialState = tsTransform.initialFreeState(it,&bField);
251			matchInfo = trackAssociator.associate(event, setup, assocParams_, &initialState);
252			}
253	bendavid	1.36
254			if (matchInfo.isGoodEcal) {
255			outTrack->SetEtaEcal(matchInfo.trkGlobPosAtEcal.eta());
256			outTrack->SetPhiEcal(matchInfo.trkGlobPosAtEcal.phi());
257			}
258	bendavid	1.23
259			//fill supercluster link
260			if (barrelSuperClusterIdMap_ \|\| endcapSuperClusterIdMap_) {
261			mithep::SuperCluster *cluster = 0;
262	loizides	1.29	for (std::vector<const ::CaloTower*>::const_iterator iTower =
263			matchInfo.crossedTowers.begin();
264			iTower<matchInfo.crossedTowers.end() && !cluster; iTower++) {
265	bendavid	1.23
266			for (uint ihit=0; ihit<(*iTower)->constituentsSize() && !cluster; ihit++) {
267			DetId hit = (*iTower)->constituent(ihit);
268			if (barrelSuperClusterIdMap_ && barrelSuperClusterIdMap_->HasMit(hit))
269			cluster = barrelSuperClusterIdMap_->GetMit(hit);
270			else if (endcapSuperClusterIdMap_ && endcapSuperClusterIdMap_->HasMit(hit))
271			cluster = endcapSuperClusterIdMap_->GetMit(hit);
272			}
273			}
274			if (cluster)
275			outTrack->SetSCluster(cluster);
276			}
277			}
278	bendavid	1.14
279	loizides	1.1	// add reference between mithep and edm object
280	loizides	1.30	if (trackMap_) {
281			mitedm::TrackPtr thePtr = inTracks.ptrAt(it - inTracks.begin());
282			trackMap_->Add(thePtr, outTrack);
283			}
284
285	bendavid	1.35	//do sim associations
286	loizides	1.19	if (trackingMap_ && !edmSimAssocName_.empty()) {
287	loizides	1.29	if (verbose_>1)
288			printf("Trying Track-Sim association\n");
289	bendavid	1.24	reco::TrackBaseRef theBaseRef = inTracks.refAt(it - inTracks.begin());
290	loizides	1.1	vector<pair<TrackingParticleRef, double> > simRefs;
291	paus	1.10	Bool_t noSimParticle = 0;
292	loizides	1.1	try {
293			simRefs = simAssociation[theBaseRef]; //try to get the sim references if existing
294			}
295	loizides	1.5	catch (edm::Exception &ex) {
296	paus	1.10	noSimParticle = 1;
297	loizides	1.1	}
298	paus	1.10
299	loizides	1.1	if (!noSimParticle) { //loop through sim match candidates
300	loizides	1.29	if (verbose_>1)
301			printf("Applying track-sim association\n");
302	loizides	1.1	for (vector<pair<TrackingParticleRef, double> >::const_iterator simRefPair=simRefs.begin();
303			simRefPair != simRefs.end(); ++simRefPair)
304	paus	1.10	if (simRefPair->second > 0.5) // require more than 50% shared hits between reco and sim
305	bendavid	1.17	outTrack->SetMCPart(trackingMap_->GetMit(simRefPair->first)); //add reco->sim reference
306	loizides	1.1	}
307			}
308			}
309			tracks_->Trim();
310			}