| 1 |
#include "../interface/MCAssociator.h"
|
| 2 |
|
| 3 |
|
| 4 |
MCAssociator::MCAssociator(): verbosity_(0), nMC_(0), mcParticles_(0), genParticles_(), mcParticlesMap_()
|
| 5 |
{
|
| 6 |
}
|
| 7 |
|
| 8 |
|
| 9 |
MCAssociator::MCAssociator(const edm::ParameterSet& producersNames, int verbosity) : verbosity_(verbosity), nMC_(0), mcParticles_(0), genParticles_(), mcParticlesMap_()
|
| 10 |
{
|
| 11 |
genParticlesProducer_ = producersNames.getParameter<edm::InputTag>("genParticlesProducer");
|
| 12 |
}
|
| 13 |
|
| 14 |
|
| 15 |
|
| 16 |
void MCAssociator::init(const edm::Event& iEvent, TClonesArray* mcParticles)
|
| 17 |
{
|
| 18 |
// FIXME - Protect against no genParticles collection in PoolSource
|
| 19 |
iEvent.getByLabel( genParticlesProducer_, genParticles_ );
|
| 20 |
|
| 21 |
// fill map<igen,imc> where igen=index in genParticle collection and imc=index in mcParticles TClonesArray
|
| 22 |
if(verbosity_>1) cout << endl << "Matching recoParticles to mcParticles... " << endl;
|
| 23 |
int igen;
|
| 24 |
mcParticles_ = mcParticles;
|
| 25 |
nMC_ = mcParticles_->GetEntriesFast();
|
| 26 |
for (int imc=0; imc<nMC_; imc++)
|
| 27 |
{
|
| 28 |
// TODO - remove indexInList in TRootMCParticle
|
| 29 |
//igen = ( (TRootParticle*)mcParticles_->At(imc))->genParticleIndex_();
|
| 30 |
igen = ( (TRootMCParticle*)mcParticles_->At(imc))->indexInList();
|
| 31 |
mcParticlesMap_[igen]=imc;
|
| 32 |
}
|
| 33 |
}
|
| 34 |
|
| 35 |
|
| 36 |
void MCAssociator::matchGenParticlesTo(TClonesArray* recoParticles)
|
| 37 |
{
|
| 38 |
int igen;
|
| 39 |
std::map<int,int>::iterator it;
|
| 40 |
|
| 41 |
for (int ipart=0; ipart<recoParticles->GetEntriesFast(); ipart++)
|
| 42 |
{
|
| 43 |
TRootParticle* recoParticle = (TRootParticle*)recoParticles->At(ipart);
|
| 44 |
igen = recoParticle->genParticleIndex();
|
| 45 |
if(igen<0)
|
| 46 |
{
|
| 47 |
if(verbosity_>2) cout <<" "<< recoParticle->typeName() << "[" << ipart << "] not matched (at CMSSW level)" << endl;
|
| 48 |
continue;
|
| 49 |
}
|
| 50 |
|
| 51 |
it=mcParticlesMap_.find(igen);
|
| 52 |
if(it==mcParticlesMap_.end())
|
| 53 |
{
|
| 54 |
if(verbosity_>2) cout <<" "<< recoParticle->typeName() << "[" << ipart << "] not matched (at TotoAna level)... add new TRootMCParticle..." << endl;
|
| 55 |
// if igen not found in mcParticles[], add genParticle in mcParticles[]
|
| 56 |
const reco::GenParticle & p = (*genParticles_)[igen];
|
| 57 |
//find the mother ID
|
| 58 |
Int_t motherID = 0; Int_t grannyID = 0;
|
| 59 |
if (p.numberOfMothers() > 0 )
|
| 60 |
{
|
| 61 |
//sanity check
|
| 62 |
const reco::Candidate * mom = p.mother();
|
| 63 |
motherID = mom->pdgId();
|
| 64 |
if (mom->numberOfMothers() > 0)
|
| 65 |
{
|
| 66 |
const reco::Candidate * granny = mom->mother();
|
| 67 |
grannyID = granny->pdgId();
|
| 68 |
if ( motherID == p.pdgId() )
|
| 69 |
{
|
| 70 |
//check if the particle is "daugther of itself"
|
| 71 |
motherID = granny->pdgId();
|
| 72 |
if (granny->numberOfMothers() > 0) grannyID = (granny->mother())->pdgId();
|
| 73 |
}
|
| 74 |
}
|
| 75 |
}
|
| 76 |
|
| 77 |
TRootMCParticle localMCParticle( p.px(), p.py(), p.pz(), p.energy(), p.vx(), p.vy(), p.vz(), p.pdgId(), p.charge(), p.status(), p.numberOfDaughters(), motherID, grannyID, 0, 0, 0, 0, igen );
|
| 78 |
new( (*mcParticles_)[nMC_] ) TRootMCParticle(localMCParticle);
|
| 79 |
recoParticle->setGenParticle( (TObject*)mcParticles_->At(nMC_) );
|
| 80 |
if(verbosity_>2) cout <<" ===> now matched to mcParticle["<< nMC_<<"] "<< localMCParticle << endl;
|
| 81 |
nMC_++;
|
| 82 |
}
|
| 83 |
else
|
| 84 |
{
|
| 85 |
if(verbosity_>2) cout <<" "<< ( (TRootParticle*)recoParticles->At(ipart))->typeName() << "[" << ipart << "] matched to mcParticles[" << it->second << "]" << endl;
|
| 86 |
recoParticle->setGenParticle( (TObject*)mcParticles_->At(it->second) );
|
| 87 |
}
|
| 88 |
}
|
| 89 |
}
|
| 90 |
|
| 91 |
|
| 92 |
void MCAssociator::printParticleAssociation(TClonesArray* recoParticles)
|
| 93 |
{
|
| 94 |
std::map<int,int>::iterator it;
|
| 95 |
for (int ipart=0; ipart<recoParticles->GetEntriesFast(); ipart++)
|
| 96 |
{
|
| 97 |
TRootParticle* localParticle = (TRootParticle*)recoParticles->At(ipart);
|
| 98 |
it=mcParticlesMap_.find(localParticle->genParticleIndex());
|
| 99 |
int imc = it->second;
|
| 100 |
cout << endl <<" "<< localParticle->typeName() <<"["<< setw(3) << ipart << "] - Charge=" << setw(2) << localParticle->charge()
|
| 101 |
<< " (Et,eta,phi)=("<< setw(9) << localParticle->Et() <<","<< setw(9) << localParticle->Eta() <<","<< setw(9) << localParticle->Phi() << ")"
|
| 102 |
<< " vertex(x,y,z)=("<< setw(11) << localParticle->vx() <<","<< setw(11) << localParticle->vy() <<","<< setw(11) << localParticle->vz() << ")" << endl;
|
| 103 |
if(it==mcParticlesMap_.end())
|
| 104 |
{
|
| 105 |
cout <<" ===> not matched" << endl;
|
| 106 |
} else {
|
| 107 |
cout << " ===> matched to mcParticles[" << setw(3) << imc << "] - Charge=" << setw(2) << localParticle->genParticle()->charge()
|
| 108 |
<< " (Et,eta,phi)=("<< setw(9) << localParticle->genParticle()->Et() <<","<< setw(9) << localParticle->genParticle()->Eta() <<","<< setw(9) << localParticle->genParticle()->Phi() << ")"
|
| 109 |
<< " vertex(x,y,z)=("<< setw(11) << localParticle->genParticle()->vx() <<","<< setw(11) << localParticle->genParticle()->vy() <<","<< setw(11) << localParticle->genParticle()->vz() << ")" << endl;
|
| 110 |
}
|
| 111 |
}
|
| 112 |
}
|
| 113 |
|
| 114 |
|
| 115 |
void MCAssociator::processGenJets(TClonesArray* genJets, TClonesArray* recoJets)
|
| 116 |
{
|
| 117 |
|
| 118 |
TRootJet* localJet;
|
| 119 |
TRootParticle* localGenJet;
|
| 120 |
if(verbosity_>1) cout << endl << "Matching genJets to recoJets... " << endl;
|
| 121 |
|
| 122 |
Double_t energy;
|
| 123 |
map<Double_t,Int_t> genMap;
|
| 124 |
map<Double_t,Int_t>::iterator it;
|
| 125 |
|
| 126 |
for (int igen=0; igen<genJets->GetEntriesFast(); igen++)
|
| 127 |
{
|
| 128 |
energy= ((TRootParticle*) genJets->At(igen))->Energy();
|
| 129 |
genMap.insert ( pair<Double_t,Int_t>(energy,igen) );
|
| 130 |
}
|
| 131 |
|
| 132 |
for (int irec=0; irec<recoJets->GetEntriesFast(); irec++)
|
| 133 |
{
|
| 134 |
localJet = (TRootJet*) recoJets->At(irec);
|
| 135 |
energy = localJet->genJetEnergy();
|
| 136 |
if(energy<0)
|
| 137 |
{
|
| 138 |
if(verbosity_>2) cout <<" "<< "TRootJet[" << irec << "] not matched to GenJet (at CMSSW level)" << endl;
|
| 139 |
continue;
|
| 140 |
}
|
| 141 |
it=genMap.find(energy);
|
| 142 |
if(it==genMap.end())
|
| 143 |
{
|
| 144 |
if(verbosity_>2) cout <<" "<< "TRootJet[" << irec << "] not matched (at TotoAna level)" << endl;
|
| 145 |
// TODO - add missing GenJet
|
| 146 |
continue;
|
| 147 |
}
|
| 148 |
else
|
| 149 |
{
|
| 150 |
if(verbosity_>2) cout <<" "<< "TRootJet[" << irec << "] matched to GenJet[" << it->second << "]" << endl;
|
| 151 |
localGenJet = (TRootParticle*) genJets->At((*it).second);
|
| 152 |
localJet->setGenJetIndex( (*it).second );
|
| 153 |
localJet->setGenJet( localGenJet );
|
| 154 |
//TODO - Add TRef in GenJet pointing to TRootJet ?
|
| 155 |
// localGenJet->setRecoJet( localJet );
|
| 156 |
}
|
| 157 |
}
|
| 158 |
|
| 159 |
genMap.clear();
|
| 160 |
}
|
| 161 |
|
| 162 |
|
| 163 |
void MCAssociator::printJetAssociation(TClonesArray* recoJets)
|
| 164 |
{
|
| 165 |
std::map<int,int>::iterator it;
|
| 166 |
for(int ipart=0; ipart<recoJets->GetEntriesFast(); ipart++)
|
| 167 |
{
|
| 168 |
|
| 169 |
TRootJet* localJet = (TRootJet*)recoJets->At(ipart);
|
| 170 |
it=mcParticlesMap_.find(localJet->genParticleIndex());
|
| 171 |
int imc = it->second;
|
| 172 |
cout << endl << " TRootJet[" << setw(3) << ipart << "] - Charge=" << setw(2) << localJet->charge()
|
| 173 |
<< " (Et,eta,phi)=("<< setw(9) << localJet->Et() <<","<< setw(9) << localJet->Eta() <<","<< setw(9) << localJet->Phi() << ")"
|
| 174 |
<< " vertex(x,y,z)=("<< setw(11) << localJet->vx() <<","<< setw(11) << localJet->vy() <<","<< setw(11) << localJet->vz() << ")" << endl;
|
| 175 |
|
| 176 |
if( localJet->genJet()!=0 )
|
| 177 |
{
|
| 178 |
cout << " ===> matched to genJet[" << setw(3) << localJet->genJetIndex() << "] - Charge=" << setw(2) << localJet->genJet()->charge()
|
| 179 |
<< " (Et,eta,phi)=("<< setw(9) << localJet->genJet()->Et() <<","<< setw(9) << localJet->genJet()->Eta() <<","<< setw(9) << localJet->genJet()->Phi() << ")"
|
| 180 |
<< " vertex(x,y,z)=("<< setw(11) << localJet->genJet()->vx() <<","<< setw(11) << localJet->genJet()->vy() <<","<< setw(11) << localJet->genJet()->vz() << ")" << endl;
|
| 181 |
} else {
|
| 182 |
cout <<" ===> not matched to genJet" << endl;
|
| 183 |
}
|
| 184 |
|
| 185 |
if(it==mcParticlesMap_.end())
|
| 186 |
{
|
| 187 |
cout <<" ===> not matched to mcParticles" << endl;
|
| 188 |
} else {
|
| 189 |
cout << " ===> matched to mcParticles[" << setw(3) << imc << "] - Charge=" << setw(2) << localJet->genParticle()->charge()
|
| 190 |
<< " (Et,eta,phi)=("<< setw(9) << localJet->genParticle()->Et() <<","<< setw(9) << localJet->genParticle()->Eta() <<","<< setw(9) << localJet->genParticle()->Phi() << ")"
|
| 191 |
<< " vertex(x,y,z)=("<< setw(11) << localJet->genParticle()->vx() <<","<< setw(11) << localJet->genParticle()->vy() <<","<< setw(11) << localJet->genParticle()->vz() << ")" << endl;
|
| 192 |
}
|
| 193 |
|
| 194 |
}
|
| 195 |
}
|
