TreeFiller/src/FillerDecayParts.cc

// $Id: FillerDecayParts.cc,v 1.20 2010/03/18 20:21:00 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()));
    }

    if(0) {
      printf("In  Particle position: x = %5f, y = %5f, z = %5f\n",p.position().x(),p.position().y(),p.position().z());
      printf("Out Particle position: x = %5f, y = %5f, z = %5f\n",d->Position().x(),d->Position().y(),d->Position().z());
    }      

    //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.21
Committed:	Wed Jun 23 09:02:45 2010 UTC (14 years, 10 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, Mit_014d, Mit_014c, Mit_014b
Changes since 1.20:	+8 -2 lines
Log Message:	some extra conversion variables
#	User	Rev	Content
1	bendavid	1.21	// $Id: FillerDecayParts.cc,v 1.20 2010/03/18 20:21:00 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	bendavid	1.21
117			if(0) {
118			printf("In Particle position: x = %5f, y = %5f, z = %5f\n",p.position().x(),p.position().y(),p.position().z());
119			printf("Out Particle position: x = %5f, y = %5f, z = %5f\n",d->Position().x(),d->Position().y(),d->Position().z());
120			}
121
122	bendavid	1.19	//fill shared layer bitmask
123			const reco::HitPattern &sharedHitPattern = p.sharedHits();
124			BitMask48 sharedLayers;
125			UInt_t numShared=0;
126			// search for all good crossed layers (with or without hit)
127			for (Int_t hi=0; hi<sharedHitPattern.numberOfHits(); hi++) {
128			uint32_t hit = sharedHitPattern.getHitPattern(hi);
129			if (sharedHitPattern.getHitType(hit)<=1) {
130			if (sharedHitPattern.trackerHitFilter(hit)) {
131			numShared++;
132			sharedLayers.SetBit(hitReader_.Layer(hit));
133			}
134			}
135			}
136
137			d->SetSharedLayers(sharedLayers);
138
139	bendavid	1.6	if (basePartMaps_.size()) {
140	loizides	1.16	// loop through and add stable daughters
141	bendavid	1.9	for (Int_t j=0; j<p.nStableChild(); ++j) {
142			const mitedm::StableData &stable = p.getStableData(j);
143			mitedm::BasePartPtr thePtr = stable.originalPtr();
144			mithep::Particle *daughter = getMitParticle(thePtr);
145
146			mithep::StableData *outStable = stableData_->Allocate();
147	loizides	1.16	new (outStable) mithep::StableData(daughter,
148			stable.p3().x(),stable.p3().y(),stable.p3().z());
149	bendavid	1.6
150	loizides	1.16	// fill bad layer mask information using corrected HitPattern
151	bendavid	1.11	const ChargedParticle cDaughter = dynamic_cast<const ChargedParticle>(daughter);
152	bendavid	1.14	if (stable.HitsFilled() && cDaughter) {
153	bendavid	1.11	const reco::HitPattern &hitPattern = stable.Hits();
154	bendavid	1.13	BitMask48 crossedLayers;
155	bendavid	1.11	UInt_t numCrossed=0;
156	loizides	1.16	// search for all good crossed layers (with or without hit)
157	bendavid	1.11	for (Int_t hi=0; hi<hitPattern.numberOfHits(); hi++) {
158			uint32_t hit = hitPattern.getHitPattern(hi);
159	loizides	1.16	if (hitPattern.getHitType(hit)<=1) {
160	bendavid	1.11	if (hitPattern.trackerHitFilter(hit)) {
161			numCrossed++;
162			crossedLayers.SetBit(hitReader_.Layer(hit));
163			}
164	loizides	1.16	}
165	bendavid	1.11	}
166	loizides	1.16
167	bendavid	1.13	BitMask48 badLayers = crossedLayers ^ cDaughter->Trk()->Hits();
168	bendavid	1.11	outStable->SetBadLayers(badLayers);
169
170	loizides	1.16	if(verbose_>1) {
171	loizides	1.12	if (outStable->NWrongHits()) {
172			printf("numCrossed: %i\n",numCrossed);
173			printf("Hit Pattern Size: %i\n",hitPattern.numberOfHits());
174			printf("# of hits: %i\n",cDaughter->Trk()->NHits());
175			printf("# of crossed layers: %i\n",crossedLayers.NBitsSet());
176			printf("# of wrong/missing layers: %i\n",badLayers.NBitsSet());
177			printf("# of missed hits: %i\n",outStable->NMissedHits());
178			printf("# of wrong hits: %i\n",outStable->NWrongHits());
179			printf("Hits : ");
180			for (Int_t biti=63; biti >= 0; --biti) {
181			printf("%i",cDaughter->Trk()->Hits().TestBit(biti));
182			}
183			printf("\n");
184			printf("Crossed Layers : ");
185			for (Int_t biti=63; biti >= 0; --biti) {
186			printf("%i",crossedLayers.TestBit(biti));
187			}
188			printf("\n");
189			printf("WrongOrMissingHits: ");
190			for (Int_t biti=63; biti>=0; --biti) {
191			printf("%i",badLayers.TestBit(biti));
192			}
193			printf("\n");
194			}
195			}
196	bendavid	1.10	}
197	loizides	1.12
198	bendavid	1.9	d->AddDaughterData(outStable);
199	bendavid	1.21
200	bendavid	1.9	}
201	loizides	1.12
202	loizides	1.16	// loop through and add decay daughters
203	bendavid	1.9	for (Int_t j=0; j<p.nDecayChild(); ++j) {
204			const mitedm::DecayData &decay = p.getDecayData(j);
205			mitedm::BasePartPtr thePtr = decay.originalPtr();
206			mithep::Particle *daughter = getMitParticle(thePtr);
207
208			mithep::DecayData *outDecay = decayData_->Allocate();
209			new (outDecay) mithep::DecayData(daughter,decay.p4());
210			outDecay->SetMassError(decay.massErr());
211			outDecay->SetLxy(decay.lxy());
212			outDecay->SetLxyError(decay.lxyErr());
213			outDecay->SetLz(decay.lz());
214			outDecay->SetLzError(decay.lzErr());
215			outDecay->SetDxy(decay.dxy());
216			outDecay->SetDxyError(decay.dxyErr());
217
218			d->AddDaughterData(outDecay);
219	bendavid	1.2	}
220			}
221	paus	1.1	}
222			decays_->Trim();
223	bendavid	1.9	stableData_->Trim();
224			decayData_->Trim();
225			}
226
227			//--------------------------------------------------------------------------------------------------
228			mithep::Particle *FillerDecayParts::getMitParticle(mitedm::BasePartPtr ptr) const
229			{
230	loizides	1.12	// Return our particle referenced by the edm pointer.
231
232	bendavid	1.9	mithep::Particle *mitPart = 0;
233			for (std::vector<const mithep::BasePartMap*>::const_iterator bmap = basePartMaps_.begin();
234			bmap!=basePartMaps_.end(); ++bmap) {
235			const mithep::BasePartMap theMap = bmap;
236	loizides	1.16	if (theMap->HasMit(ptr)) {
237	bendavid	1.9	mitPart = theMap->GetMit(ptr);
238	loizides	1.16	return mitPart;
239			}
240	bendavid	1.9	}
241
242			if (!mitPart)
243			throw edm::Exception(edm::errors::Configuration, "FillerDecayParts::FillDataBlock()\n")
244			<< "Error! MITHEP Object "
245			<< "not found in AssociationMaps (" << typeid(*this).name() << ")." << std::endl;
246
247			return mitPart;
248	paus	1.1	}