TreeFiller/src/FillerDecayParts.cc

// $Id: FillerDecayParts.cc,v 1.19 2009/12/15 23:27:54 bendavid Exp $

#include "MitProd/TreeFiller/interface/FillerDecayParts.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 "DataFormats/RecoCandidate/interface/TrackAssociation.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
#include "MitAna/DataTree/interface/StableDataCol.h"
#include "MitAna/DataTree/interface/DecayDataCol.h"
#include "MitAna/DataTree/interface/DecayParticleCol.h"
#include "MitAna/DataTree/interface/Names.h"
#include "MitEdm/DataFormats/interface/Collections.h"
#include "MitProd/ObjectService/interface/ObjectService.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()));
    }
      
    //fill shared layer bitmask
    const reco::HitPattern &sharedHitPattern = p.sharedHits();
    BitMask48 sharedLayers;
    UInt_t numShared=0;
    // search for all good crossed layers (with or without hit)
    for (Int_t hi=0; hi<sharedHitPattern.numberOfHits(); hi++) {
      uint32_t hit = sharedHitPattern.getHitPattern(hi);
      if (sharedHitPattern.getHitType(hit)<=1) {
        if (sharedHitPattern.trackerHitFilter(hit)) {
          numShared++;
          sharedLayers.SetBit(hitReader_.Layer(hit));
        }
      }
    }

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