TreeFiller/src/FillerDecayParts.cc

// $Id: FillerDecayParts.cc,v 1.15 2009/02/26 17:04:03 bendavid Exp $

#include "MitAna/DataTree/interface/DecayParticle.h"
#include "MitAna/DataTree/interface/DaughterData.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/TrackReco/interface/Track.h"
#include "DataFormats/TrackReco/interface/TrackFwd.h"
#include "DataFormats/TrackReco/interface/HitPattern.h"
#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "DataFormats/RecoCandidate/interface/TrackAssociation.h"
#include "MitEdm/DataFormats/interface/Collections.h"
#include "MitEdm/DataFormats/interface/DecayPart.h"
#include "MitEdm/DataFormats/interface/BasePart.h"
#include "MitEdm/DataFormats/interface/BasePartFwd.h"
#include "MitProd/TreeFiller/interface/FillerDecayParts.h"
#include "MitAna/DataTree/interface/Names.h"

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

//--------------------------------------------------------------------------------------------------
FillerDecayParts::FillerDecayParts(const ParameterSet &cfg, const char *name, bool active) :
  BaseFiller(cfg,name,active),
  edmName_(Conf().getUntrackedParameter<string>("edmName","")),
  mitName_(Conf().getUntrackedParameter<string>("mitName","")),
  stableDataName_(mitName_ + "_StableDatas"),
  decayDataName_(mitName_ + "_DecayDatas"),
  vertexMapName_(Conf().getUntrackedParameter<string>("vertexMap","")),
  basePartMapNames_(Conf().exists("basePartMaps") ? 
                    Conf().getUntrackedParameter<vector<string> >("basePartMaps") : 
                    vector<string>()),
  vertexMap_(0),
  decays_(new mithep::DecayParticleArr(250)),
  stableData_(new mithep::StableDataArr(500)),
  decayData_(new mithep::DecayDataArr(500))
{
  // Constructor.
}

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

  delete decays_;
  delete stableData_;
  delete decayData_;
}

//--------------------------------------------------------------------------------------------------
void FillerDecayParts::BookDataBlock(TreeWriter &tws)
{
  // Add our branches to tree and get our map.
  
  tws.AddBranch(mitName_,        &decays_);
  OS()->add<mithep::DecayParticleArr>(decays_,mitName_);
  tws.AddBranch(stableDataName_, &stableData_);
  OS()->add<mithep::StableDataArr>(stableData_,stableDataName_);
  tws.AddBranch(decayDataName_,  &decayData_);
  OS()->add<mithep::DecayDataArr>(decayData_,decayDataName_);

  if (!vertexMapName_.empty()) {
    vertexMap_ = OS()->get<VertexMap>(vertexMapName_);
    if (vertexMap_)
      AddBranchDep(mitName_,vertexMap_->GetBrName());
  }

  for (std::vector<std::string>::const_iterator bmapName = basePartMapNames_.begin();
        bmapName!=basePartMapNames_.end(); ++bmapName) {
    if (!bmapName->empty()) {
      const BasePartMap *map = OS()->get<BasePartMap>(*bmapName);
      if (map) {
        basePartMaps_.push_back(map);
        AddBranchDep(stableDataName_,map->GetBrName());
        AddBranchDep(decayDataName_,map->GetBrName());
      }
    }
  }
}

//--------------------------------------------------------------------------------------------------
void FillerDecayParts::FillDataBlock(const edm::Event      &evt, 
                                     const edm::EventSetup &setup)
{
  // Fill the EDM DecayPart collection into the MIT DecayParticle collection.

  decays_->Delete();
  stableData_->Delete();
  decayData_->Delete();

  Handle<mitedm::DecayPartCol> hParts;
  GetProduct(edmName_, hParts, evt);  
  const mitedm::DecayPartCol *iParts = hParts.product();
  
  // loop through all DecayParts and fill the information
  for (UInt_t i=0; i<iParts->size(); ++i) {
    const mitedm::DecayPart &p = iParts->at(i); // for convenience
    mithep::DecayParticle *d = decays_->Allocate();
    new (d) mithep::DecayParticle(p.pid());
    
    d->SetMassError(p.fittedMassError());
    d->SetLxy(p.lxyToPv());
    d->SetLxyError(p.lxyToPvError());
    d->SetDxy(p.dxyToPv());
    d->SetDxyError(p.dxyToPvError());
    d->SetLz(p.lzToPv());
    d->SetLzError(p.lzToPvError());
    d->SetMom(p.fourMomentum());
    d->SetChi2(p.chi2());
    d->SetNdof(p.ndof());
    
    if (p.primaryVertex().isNonnull()) {
      d->SetPriVertex(vertexMap_->GetMit(p.primaryVertex()));
    }
      
    if (basePartMaps_.size()) {
       // loop through and add stable daughters
      for (Int_t j=0; j<p.nStableChild(); ++j) {
        const mitedm::StableData &stable = p.getStableData(j);
        mitedm::BasePartPtr thePtr = stable.originalPtr();
        mithep::Particle *daughter = getMitParticle(thePtr);
        
        mithep::StableData *outStable = stableData_->Allocate();
        new (outStable) mithep::StableData(daughter,
                                           stable.p3().x(),stable.p3().y(),stable.p3().z());
        
        // fill bad layer mask information using corrected HitPattern
        const ChargedParticle *cDaughter = dynamic_cast<const ChargedParticle*>(daughter);
        if (stable.HitsFilled() && cDaughter) {
          const reco::HitPattern &hitPattern = stable.Hits();
          BitMask48 crossedLayers;
          UInt_t numCrossed=0;
          // search for all good crossed layers (with or without hit)
          for (Int_t hi=0; hi<hitPattern.numberOfHits(); hi++) {
            uint32_t hit = hitPattern.getHitPattern(hi);
            if (hitPattern.getHitType(hit)<=1) {
              if (hitPattern.trackerHitFilter(hit)) {
                numCrossed++;
                crossedLayers.SetBit(hitReader_.Layer(hit));
              }
            }
          }

          BitMask48 badLayers = crossedLayers ^ cDaughter->Trk()->Hits();
          outStable->SetBadLayers(badLayers);
          
          if(verbose_>1) {
            if (outStable->NWrongHits()) {          
              printf("numCrossed:                %i\n",numCrossed);
              printf("Hit Pattern Size:          %i\n",hitPattern.numberOfHits());
              printf("# of hits:                 %i\n",cDaughter->Trk()->NHits());
              printf("# of crossed layers:       %i\n",crossedLayers.NBitsSet());            
              printf("# of wrong/missing layers: %i\n",badLayers.NBitsSet());
              printf("# of missed hits:          %i\n",outStable->NMissedHits());
              printf("# of wrong hits:           %i\n",outStable->NWrongHits());
              printf("Hits              : ");
              for (Int_t biti=63; biti >= 0; --biti) {
                printf("%i",cDaughter->Trk()->Hits().TestBit(biti));
              }
              printf("\n");
              printf("Crossed Layers    : ");
              for (Int_t biti=63; biti >= 0; --biti) {
                printf("%i",crossedLayers.TestBit(biti));
              }
              printf("\n");
              printf("WrongOrMissingHits: ");
              for (Int_t biti=63; biti>=0; --biti) {
                printf("%i",badLayers.TestBit(biti));
              }
              printf("\n");    
            }        
          }
        }

        d->AddDaughterData(outStable);
      }

      // loop through and add decay daughters
      for (Int_t j=0; j<p.nDecayChild(); ++j) {
        const mitedm::DecayData &decay = p.getDecayData(j);
        mitedm::BasePartPtr thePtr = decay.originalPtr();
        mithep::Particle *daughter = getMitParticle(thePtr);
        
        mithep::DecayData *outDecay = decayData_->Allocate();
        new (outDecay) mithep::DecayData(daughter,decay.p4());
        outDecay->SetMassError(decay.massErr());
        outDecay->SetLxy(decay.lxy());
        outDecay->SetLxyError(decay.lxyErr());
        outDecay->SetLz(decay.lz());
        outDecay->SetLzError(decay.lzErr());
        outDecay->SetDxy(decay.dxy());
        outDecay->SetDxyError(decay.dxyErr());
        
        d->AddDaughterData(outDecay);
      }
    }
  }
  decays_->Trim();
  stableData_->Trim();
  decayData_->Trim();
}

//--------------------------------------------------------------------------------------------------
mithep::Particle *FillerDecayParts::getMitParticle(mitedm::BasePartPtr ptr) const
{
  // Return our particle referenced by the edm pointer.

  mithep::Particle *mitPart = 0;
  for (std::vector<const mithep::BasePartMap*>::const_iterator bmap = basePartMaps_.begin();
        bmap!=basePartMaps_.end(); ++bmap) {
    const mithep::BasePartMap *theMap = *bmap;
    if (theMap->HasMit(ptr)) {
      mitPart = theMap->GetMit(ptr);
      return mitPart;
    }
  }
  
  if (!mitPart)
    throw edm::Exception(edm::errors::Configuration, "FillerDecayParts::FillDataBlock()\n")
    << "Error! MITHEP Object " 
    << "not found in AssociationMaps (" << typeid(*this).name() << ")." << std::endl;
    
  return mitPart;
}
Revision:	1.16
Committed:	Sun Mar 15 11:20:41 2009 UTC (16 years, 1 month ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_009a, Mit_009, Mit_008, Mit_008pre2
Changes since 1.15:	+39 -23 lines
Log Message:	Introduced BranchTable plus general cleanup.
#	User	Rev	Content
1	loizides	1.16	// $Id: FillerDecayParts.cc,v 1.15 2009/02/26 17:04:03 bendavid Exp $
2	paus	1.1
3	loizides	1.4	#include "MitAna/DataTree/interface/DecayParticle.h"
4	bendavid	1.6	#include "MitAna/DataTree/interface/DaughterData.h"
5	paus	1.1	#include "FWCore/MessageLogger/interface/MessageLogger.h"
6	bendavid	1.10	#include "FWCore/Framework/interface/ESHandle.h"
7	paus	1.1	#include "DataFormats/Common/interface/Handle.h"
8			#include "DataFormats/TrackReco/interface/Track.h"
9			#include "DataFormats/TrackReco/interface/TrackFwd.h"
10	bendavid	1.10	#include "DataFormats/TrackReco/interface/HitPattern.h"
11			#include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
12	paus	1.1	#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
13			#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
14			#include "DataFormats/RecoCandidate/interface/TrackAssociation.h"
15	loizides	1.8	#include "MitEdm/DataFormats/interface/Collections.h"
16	bendavid	1.2	#include "MitEdm/DataFormats/interface/DecayPart.h"
17			#include "MitEdm/DataFormats/interface/BasePart.h"
18			#include "MitEdm/DataFormats/interface/BasePartFwd.h"
19	paus	1.1	#include "MitProd/TreeFiller/interface/FillerDecayParts.h"
20			#include "MitAna/DataTree/interface/Names.h"
21
22			using namespace std;
23			using namespace edm;
24			using namespace mithep;
25
26			//--------------------------------------------------------------------------------------------------
27	loizides	1.4	FillerDecayParts::FillerDecayParts(const ParameterSet &cfg, const char *name, bool active) :
28	paus	1.1	BaseFiller(cfg,name,active),
29	loizides	1.4	edmName_(Conf().getUntrackedParameter<string>("edmName","")),
30			mitName_(Conf().getUntrackedParameter<string>("mitName","")),
31	loizides	1.16	stableDataName_(mitName_ + "_StableDatas"),
32			decayDataName_(mitName_ + "_DecayDatas"),
33	bendavid	1.9	vertexMapName_(Conf().getUntrackedParameter<string>("vertexMap","")),
34	loizides	1.8	basePartMapNames_(Conf().exists("basePartMaps") ?
35			Conf().getUntrackedParameter<vector<string> >("basePartMaps") :
36			vector<string>()),
37	bendavid	1.9	vertexMap_(0),
38			decays_(new mithep::DecayParticleArr(250)),
39			stableData_(new mithep::StableDataArr(500)),
40			decayData_(new mithep::DecayDataArr(500))
41	paus	1.1	{
42	loizides	1.4	// Constructor.
43	paus	1.1	}
44
45			//--------------------------------------------------------------------------------------------------
46			FillerDecayParts::~FillerDecayParts()
47			{
48	loizides	1.4	// Destructor.
49
50	paus	1.1	delete decays_;
51	bendavid	1.9	delete stableData_;
52			delete decayData_;
53	paus	1.1	}
54
55			//--------------------------------------------------------------------------------------------------
56			void FillerDecayParts::BookDataBlock(TreeWriter &tws)
57			{
58	loizides	1.12	// Add our branches to tree and get our map.
59	bendavid	1.9
60	loizides	1.16	tws.AddBranch(mitName_, &decays_);
61			OS()->add<mithep::DecayParticleArr>(decays_,mitName_);
62			tws.AddBranch(stableDataName_, &stableData_);
63			OS()->add<mithep::StableDataArr>(stableData_,stableDataName_);
64			tws.AddBranch(decayDataName_, &decayData_);
65			OS()->add<mithep::DecayDataArr>(decayData_,decayDataName_);
66
67			if (!vertexMapName_.empty()) {
68			vertexMap_ = OS()->get<VertexMap>(vertexMapName_);
69			if (vertexMap_)
70			AddBranchDep(mitName_,vertexMap_->GetBrName());
71			}
72	loizides	1.4
73	bendavid	1.6	for (std::vector<std::string>::const_iterator bmapName = basePartMapNames_.begin();
74			bmapName!=basePartMapNames_.end(); ++bmapName) {
75	loizides	1.16	if (!bmapName->empty()) {
76			const BasePartMap map = OS()->get<BasePartMap>(bmapName);
77			if (map) {
78			basePartMaps_.push_back(map);
79			AddBranchDep(stableDataName_,map->GetBrName());
80			AddBranchDep(decayDataName_,map->GetBrName());
81			}
82			}
83	bendavid	1.6	}
84	paus	1.1	}
85
86			//--------------------------------------------------------------------------------------------------
87			void FillerDecayParts::FillDataBlock(const edm::Event &evt,
88			const edm::EventSetup &setup)
89			{
90	loizides	1.16	// Fill the EDM DecayPart collection into the MIT DecayParticle collection.
91	loizides	1.12
92	bendavid	1.15	decays_->Delete();
93			stableData_->Delete();
94			decayData_->Delete();
95	loizides	1.4
96	bendavid	1.2	Handle<mitedm::DecayPartCol> hParts;
97	paus	1.1	GetProduct(edmName_, hParts, evt);
98	bendavid	1.2	const mitedm::DecayPartCol *iParts = hParts.product();
99	paus	1.1
100	loizides	1.16	// loop through all DecayParts and fill the information
101	bendavid	1.2	for (UInt_t i=0; i<iParts->size(); ++i) {
102	loizides	1.16	const mitedm::DecayPart &p = iParts->at(i); // for convenience
103	bendavid	1.2	mithep::DecayParticle *d = decays_->Allocate();
104	bendavid	1.9	new (d) mithep::DecayParticle(p.pid());
105	bendavid	1.2
106	bendavid	1.9	d->SetMassError(p.fittedMassError());
107			d->SetLxy(p.lxyToPv());
108			d->SetLxyError(p.lxyToPvError());
109			d->SetDxy(p.dxyToPv());
110			d->SetDxyError(p.dxyToPvError());
111			d->SetLz(p.lzToPv());
112			d->SetLzError(p.lzToPvError());
113	bendavid	1.2	d->SetMom(p.fourMomentum());
114	mrudolph	1.5	d->SetChi2(p.chi2());
115			d->SetNdof(p.ndof());
116	bendavid	1.10
117			if (p.primaryVertex().isNonnull()) {
118	bendavid	1.9	d->SetPriVertex(vertexMap_->GetMit(p.primaryVertex()));
119	bendavid	1.10	}
120
121	bendavid	1.6	if (basePartMaps_.size()) {
122	loizides	1.16	// loop through and add stable daughters
123	bendavid	1.9	for (Int_t j=0; j<p.nStableChild(); ++j) {
124			const mitedm::StableData &stable = p.getStableData(j);
125			mitedm::BasePartPtr thePtr = stable.originalPtr();
126			mithep::Particle *daughter = getMitParticle(thePtr);
127
128			mithep::StableData *outStable = stableData_->Allocate();
129	loizides	1.16	new (outStable) mithep::StableData(daughter,
130			stable.p3().x(),stable.p3().y(),stable.p3().z());
131	bendavid	1.6
132	loizides	1.16	// fill bad layer mask information using corrected HitPattern
133	bendavid	1.11	const ChargedParticle cDaughter = dynamic_cast<const ChargedParticle>(daughter);
134	bendavid	1.14	if (stable.HitsFilled() && cDaughter) {
135	bendavid	1.11	const reco::HitPattern &hitPattern = stable.Hits();
136	bendavid	1.13	BitMask48 crossedLayers;
137	bendavid	1.11	UInt_t numCrossed=0;
138	loizides	1.16	// search for all good crossed layers (with or without hit)
139	bendavid	1.11	for (Int_t hi=0; hi<hitPattern.numberOfHits(); hi++) {
140			uint32_t hit = hitPattern.getHitPattern(hi);
141	loizides	1.16	if (hitPattern.getHitType(hit)<=1) {
142	bendavid	1.11	if (hitPattern.trackerHitFilter(hit)) {
143			numCrossed++;
144			crossedLayers.SetBit(hitReader_.Layer(hit));
145			}
146	loizides	1.16	}
147	bendavid	1.11	}
148	loizides	1.16
149	bendavid	1.13	BitMask48 badLayers = crossedLayers ^ cDaughter->Trk()->Hits();
150	bendavid	1.11	outStable->SetBadLayers(badLayers);
151
152	loizides	1.16	if(verbose_>1) {
153	loizides	1.12	if (outStable->NWrongHits()) {
154			printf("numCrossed: %i\n",numCrossed);
155			printf("Hit Pattern Size: %i\n",hitPattern.numberOfHits());
156			printf("# of hits: %i\n",cDaughter->Trk()->NHits());
157			printf("# of crossed layers: %i\n",crossedLayers.NBitsSet());
158			printf("# of wrong/missing layers: %i\n",badLayers.NBitsSet());
159			printf("# of missed hits: %i\n",outStable->NMissedHits());
160			printf("# of wrong hits: %i\n",outStable->NWrongHits());
161			printf("Hits : ");
162			for (Int_t biti=63; biti >= 0; --biti) {
163			printf("%i",cDaughter->Trk()->Hits().TestBit(biti));
164			}
165			printf("\n");
166			printf("Crossed Layers : ");
167			for (Int_t biti=63; biti >= 0; --biti) {
168			printf("%i",crossedLayers.TestBit(biti));
169			}
170			printf("\n");
171			printf("WrongOrMissingHits: ");
172			for (Int_t biti=63; biti>=0; --biti) {
173			printf("%i",badLayers.TestBit(biti));
174			}
175			printf("\n");
176			}
177			}
178	bendavid	1.10	}
179	loizides	1.12
180	bendavid	1.9	d->AddDaughterData(outStable);
181			}
182	loizides	1.12
183	loizides	1.16	// loop through and add decay daughters
184	bendavid	1.9	for (Int_t j=0; j<p.nDecayChild(); ++j) {
185			const mitedm::DecayData &decay = p.getDecayData(j);
186			mitedm::BasePartPtr thePtr = decay.originalPtr();
187			mithep::Particle *daughter = getMitParticle(thePtr);
188
189			mithep::DecayData *outDecay = decayData_->Allocate();
190			new (outDecay) mithep::DecayData(daughter,decay.p4());
191			outDecay->SetMassError(decay.massErr());
192			outDecay->SetLxy(decay.lxy());
193			outDecay->SetLxyError(decay.lxyErr());
194			outDecay->SetLz(decay.lz());
195			outDecay->SetLzError(decay.lzErr());
196			outDecay->SetDxy(decay.dxy());
197			outDecay->SetDxyError(decay.dxyErr());
198
199			d->AddDaughterData(outDecay);
200	bendavid	1.2	}
201			}
202	paus	1.1	}
203			decays_->Trim();
204	bendavid	1.9	stableData_->Trim();
205			decayData_->Trim();
206			}
207
208			//--------------------------------------------------------------------------------------------------
209			mithep::Particle *FillerDecayParts::getMitParticle(mitedm::BasePartPtr ptr) const
210			{
211	loizides	1.12	// Return our particle referenced by the edm pointer.
212
213	bendavid	1.9	mithep::Particle *mitPart = 0;
214			for (std::vector<const mithep::BasePartMap*>::const_iterator bmap = basePartMaps_.begin();
215			bmap!=basePartMaps_.end(); ++bmap) {
216			const mithep::BasePartMap theMap = bmap;
217	loizides	1.16	if (theMap->HasMit(ptr)) {
218	bendavid	1.9	mitPart = theMap->GetMit(ptr);
219	loizides	1.16	return mitPart;
220			}
221	bendavid	1.9	}
222
223			if (!mitPart)
224			throw edm::Exception(edm::errors::Configuration, "FillerDecayParts::FillDataBlock()\n")
225			<< "Error! MITHEP Object "
226			<< "not found in AssociationMaps (" << typeid(*this).name() << ")." << std::endl;
227
228			return mitPart;
229	paus	1.1	}