9 #include "CLHEP/Units/defs.h"
10 #include "CLHEP/Units/PhysicalConstants.h"
22 ,hepmcCollection_(iPSet.getParameter<edm::
InputTag>(
"hepmcCollection"))
23 ,tauEtCut(iPSet.getParameter<double>(
"tauEtCutForRtau"))
116 TauFSRPhotonsPt=
dbe->
book1D(
"TauFSRPhotonsPt",
"Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
130 JAKInvMass.push_back(std::vector<MonitorElement *>());
165 HepMC::GenEvent *myGenEvent =
new HepMC::GenEvent(*(evt->GetEvent()));
180 for(HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin(); iter != myGenEvent->particles_end(); iter++) {
181 if(
abs((*iter)->pdg_id())==23){
184 if(
abs((*iter)->pdg_id())==15){
195 rtau(*iter,mother,decaychannel,weight);
202 unsigned int jak_id, TauBitMask;
203 if(TD.
AnalyzeTau((*iter),jak_id,TauBitMask,
false,
false)){
206 int tcharge=(*iter)->pdg_id()/
abs((*iter)->pdg_id());
209 TLorentzVector LVQ(0,0,0,0);
210 TLorentzVector LVS12(0,0,0,0);
211 TLorentzVector LVS13(0,0,0,0);
212 TLorentzVector LVS23(0,0,0,0);
214 for(
unsigned int i=0;
i<part.size();
i++){
219 TLorentzVector
LV(part.at(
i)->momentum().px(),part.at(
i)->momentum().py(),part.at(
i)->momentum().pz(),part.at(
i)->momentum().e());
264 if ( tau->production_vertex() ) {
265 HepMC::GenVertex::particle_iterator mother;
266 for (mother = tau->production_vertex()->particles_begin(
HepMC::parents); mother!= tau->production_vertex()->particles_end(
HepMC::parents); mother++ ) {
267 if((*mother)->pdg_id() == tau->pdg_id())
return GetMother(*mother);
276 std::vector<HepMC::GenParticle*> mothers;
277 if ( boson->production_vertex() ) {
278 HepMC::GenVertex::particle_iterator mother;
279 for (mother = boson->production_vertex()->particles_begin(
HepMC::parents); mother!= boson->production_vertex()->particles_end(
HepMC::parents); mother++ ) {
280 if((*mother)->pdg_id() == boson->pdg_id())
return GetMothers(*mother);
281 mothers.push_back(*mother);
290 if ( tau->production_vertex() ) {
291 HepMC::GenVertex::particle_iterator mother;
292 for (mother = tau->production_vertex()->particles_begin(
HepMC::parents); mother!= tau->production_vertex()->particles_end(
HepMC::parents); mother++ ) {
293 mother_pid = (*mother)->pdg_id();
294 if(mother_pid == tau->pdg_id())
return findMother(*mother);
302 if ( tau->end_vertex() ) {
303 HepMC::GenVertex::particle_iterator dau;
304 for (dau = tau->end_vertex()->particles_begin(HepMC::children); dau!= tau->end_vertex()->particles_end(HepMC::children); dau++ ) {
305 int dau_pid = (*dau)->pdg_id();
306 if(dau_pid == tau->pdg_id())
return false;
314 TauList.insert(TauList.begin(),
tau);
315 if ( tau->production_vertex() ) {
316 HepMC::GenVertex::particle_iterator mother;
317 for (mother = tau->production_vertex()->particles_begin(
HepMC::parents); mother!= tau->production_vertex()->particles_end(
HepMC::parents);mother++) {
318 if((*mother)->pdg_id() == tau->pdg_id()){
326 std::vector<const HepMC::GenParticle*> &ListofBrem){
328 if(
abs(p->pdg_id())==15){
333 if ( p->end_vertex() ) {
334 HepMC::GenVertex::particle_iterator dau;
335 for (dau = p->end_vertex()->particles_begin(HepMC::children); dau!= p->end_vertex()->particles_end(HepMC::children); dau++ ) {
338 if(
abs((*dau)->pdg_id()) ==
abs(photo_ID) && !doBrem){ListofFSR.push_back(*dau);}
339 if(
abs((*dau)->pdg_id()) ==
abs(photo_ID) && doBrem){ListofBrem.push_back(*dau);}
340 if((*dau)->end_vertex() && (*dau)->end_vertex()->particles_out_size()>0 &&
abs((*dau)->pdg_id()) != 111 &&
abs((*dau)->pdg_id()) != 221){
352 double mother_pid=m->pdg_id();
353 if(m->end_vertex() && mother_pid!=p->pdg_id()){
354 HepMC::GenVertex::particle_iterator dau;
355 for (dau = m->end_vertex()->particles_begin(HepMC::children); dau!= m->end_vertex()->particles_end(HepMC::children); dau++ ) {
356 int dau_pid = (*dau)->pdg_id();
357 if(fabs(dau_pid) == 22) {
358 ListofFSR.push_back(*dau);
362 if(
abs(mother_pid) == 24) BosonScale=1.0;
363 if(
abs(mother_pid) == 23) BosonScale=2.0;
364 if(
abs(mother_pid) == 22) BosonScale=2.0;
365 if(
abs(mother_pid) == 25) BosonScale=2.0;
366 if(
abs(mother_pid) == 35) BosonScale=2.0;
367 if(
abs(mother_pid) == 36) BosonScale=2.0;
368 if(
abs(mother_pid) == 37) BosonScale=1.0;
374 if(
abs(tau->pdg_id()) != 15 )
return -3;
377 if(mother_pid == -2)
return -2;
380 if(
abs(mother_pid) == 24) label =
W;
381 if(
abs(mother_pid) == 23) label =
Z;
382 if(
abs(mother_pid) == 22) label =
gamma;
383 if(
abs(mother_pid) == 25) label =
HSM;
384 if(
abs(mother_pid) == 35) label =
H0;
385 if(
abs(mother_pid) == 36) label =
A0;
386 if(
abs(mother_pid) == 37) label =
Hpm;
388 int mother_shortpid=(
abs(mother_pid)%10000);
389 if(mother_shortpid>500 && mother_shortpid<600 )label =
B;
390 if(mother_shortpid>400 && mother_shortpid<500)label =
D;
392 if(label==
B || label ==
D || label ==
other)
return -1;
399 if ( tau->end_vertex() ) {
400 HepMC::GenVertex::particle_iterator des;
401 for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
402 des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
403 int pid = (*des)->pdg_id();
405 if((*des)->status() != 1)
continue;
408 int charge = (int) pd->charge();
409 if(charge == 0)
continue;
420 if(tau->status() == 1) channel =
stable;
431 if ( tau->end_vertex() ) {
432 HepMC::GenVertex::particle_iterator des;
433 for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
434 des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
435 int pid = (*des)->pdg_id();
439 if(
abs(pid) == 11) eCount++;
440 if(
abs(pid) == 13) muCount++;
441 if(
abs(pid) == 111) pi0Count++;
442 if(
abs(pid) == 211) piCount++;
443 if(
abs(pid) == 213) rhoCount++;
444 if(
abs(pid) == 20213) a1Count++;
445 if(
abs(pid) == 321) KCount++;
446 if(
abs(pid) == 323) KstarCount++;
451 if(KCount >= 1) channel =
K;
452 if(KstarCount >= 1) channel =
Kstar;
453 if(a1Count >= 1) channel =
a1;
454 if(rhoCount >= 1) channel =
rho;
458 if(piCount == 1 && pi0Count == 0) channel =
pi;
459 if(piCount == 1 && pi0Count == 1) channel =
pi1pi0;
460 if(piCount == 1 && pi0Count > 1) channel =
pinpi0;
461 if(piCount == 3 && pi0Count == 0) channel =
tripi;
462 if(piCount == 3 && pi0Count > 0) channel =
tripinpi0;
464 if(muCount == 1) channel =
muon;
472 if(decay !=
pi1pi0)
return;
474 if(tau->momentum().perp() <
tauEtCut)
return;
480 if(visibleTauE != 0) rTau = ltrack/visibleTauE;
488 TLorentzVector momP4 =
motherP4(tau);
490 pionP4.Boost(-1*momP4.BoostVector());
491 double energy = pionP4.E()/(momP4.M()/2);
507 TLorentzVector
rho(0,0,0,0),
pi(0,0,0,0);
508 for(
unsigned int i=0;
i<part.size();
i++){
509 TLorentzVector
LV(part.at(
i)->momentum().px(),part.at(
i)->momentum().py(),part.at(
i)->momentum().pz(),part.at(
i)->momentum().e());
517 TLorentzVector
a1(0,0,0,0),pi_p(0,0,0,0),pi_m(0,0,0,0);
518 int nplus(0),nminus(0);
519 for(
unsigned int i=0;
i<part.size();
i++){
520 TLorentzVector
LV(part.at(
i)->momentum().px(),part.at(
i)->momentum().py(),part.at(
i)->momentum().pz(),part.at(
i)->momentum().e());
525 if(nplus+nminus==3 && nplus==1) gamma=2*pi_p.Pt()/
a1.Pt()-1;
526 if(nplus+nminus==3 && nminus==1) gamma=2*pi_m.Pt()/
a1.Pt()-1;
528 pi_p+=pi_m; gamma=2*pi_p.Pt()/
a1.Pt()-1;
537 TLorentzVector tautau(0,0,0,0);
538 TLorentzVector pipi(0,0,0,0);
539 TLorentzVector taum(0,0,0,0);
540 int nSinglePionDecays(0),nSingleMuonDecays(0),nSingleElectronDecays(0);
542 TLorentzVector Zboson(boson->momentum().px(),boson->momentum().py(),boson->momentum().pz(),boson->momentum().e());
543 if ( boson->end_vertex() ) {
544 HepMC::GenVertex::particle_iterator des;
545 for(des = boson->end_vertex()->particles_begin(HepMC::children); des!= boson->end_vertex()->particles_end(HepMC::children);++des ) {
546 int pid = (*des)->pdg_id();
551 TLorentzVector LVtau((*des)->momentum().px(),(*des)->momentum().py(),(*des)->momentum().pz(),(*des)->momentum().e());
556 int charge = (int) pd->charge();
557 LVtau.Boost(-1*Zboson.BoostVector());
558 LVpi.Boost(-1*Zboson.BoostVector());
559 if(charge<0){x1=LVpi.P()/LVtau.E(); taum=LVtau;}
560 else{ x2=LVpi.P()/LVtau.E();}
564 if(nSingleMuonDecays==2){
568 if(nSingleElectronDecays==2){
572 if(nSinglePionDecays == 2 && tautau.M()!= 0) {
577 if(x2-x1>alow && x2-x1<aup){
578 double zs=(zsup+zslow)/2;
588 const std::vector<HepMC::GenParticle*>
m=
GetMothers(boson);
590 TLorentzVector
Z(0,0,0,0);
591 for(
unsigned int i=0;
i<m.size();
i++){
596 if(taum.Vect().Dot(Zboson.Vect())/(Zboson.P()*taum.P())>0){
610 double a=1-
sqrt(fabs(1.0-2*fabs(zs)));
624 TLorentzVector
p4(0,0,0,0);
626 if ( tau->end_vertex() ) {
627 HepMC::GenVertex::particle_iterator des;
628 for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
629 des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
630 int pid = (*des)->pdg_id();
634 if(!(
abs(pid)==211 ||
abs(pid)==13 ||
abs(pid)==11))
continue;
636 if((*des)->momentum().rho() > p4.P()) {
637 p4 = TLorentzVector((*des)->momentum().px(),
638 (*des)->momentum().py(),
639 (*des)->momentum().pz(),
640 (*des)->momentum().e());
650 return TLorentzVector(m->momentum().px(),m->momentum().py(),m->momentum().pz(),m->momentum().e());
654 TLorentzVector
p4(tau->momentum().px(),
655 tau->momentum().py(),
656 tau->momentum().pz(),
657 tau->momentum().e());
659 if ( tau->end_vertex() ) {
660 HepMC::GenVertex::particle_iterator des;
661 for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
662 des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
663 int pid = (*des)->pdg_id();
667 if(
abs(pid) == 12 ||
abs(pid) == 14 ||
abs(pid) == 16) {
668 p4 -= TLorentzVector((*des)->momentum().px(),
669 (*des)->momentum().py(),
670 (*des)->momentum().pz(),
671 (*des)->momentum().e());
681 std::vector<const HepMC::GenParticle*> TauList;
686 std::vector<const HepMC::GenParticle*> ListofFSR; ListofFSR.clear();
687 std::vector<const HepMC::GenParticle*> ListofBrem; ListofBrem.clear();
688 std::vector<const HepMC::GenParticle*> FSR_photos; FSR_photos.clear();
689 double BosonScale(1);
690 if(TauList.size()>0){
696 double photonPtSum=0;
697 for(
unsigned int i=0;
i<ListofBrem.size();
i++){
698 photonPtSum+=ListofBrem.at(
i)->momentum().perp();
707 for(
unsigned int i=0;
i<ListofFSR.size();
i++){
708 photonPtSum+=ListofFSR.at(
i)->momentum().perp();
711 double FSR_photosSum(0);
712 for(
unsigned int i=0;
i<FSR_photos.size();
i++){
713 FSR_photosSum+=FSR_photos.at(
i)->momentum().perp();
void findTauList(const HepMC::GenParticle *tau, std::vector< const HepMC::GenParticle * > &TauList)
MonitorElement * TauFSRPhotonsN
MonitorElement * TauSpinEffectsW_X
MonitorElement * TauSpinEffectsHpm_UpsilonA1
std::vector< std::vector< MonitorElement * > > JAKInvMass
MonitorElement * TauSpinEffectsZ_Xb
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
MonitorElement * TauRtauW
TauValidation(const edm::ParameterSet &)
void spinEffects(const HepMC::GenParticle *, int, int, std::vector< HepMC::GenParticle * > &part, double weight)
MonitorElement * TauSpinEffectsH_MVis
MonitorElement * TauMothers
MonitorElement * TauSpinEffectsHpm_eX
MonitorElement * TauSpinEffectsW_muX
MonitorElement * TauBremPhotonsN
void setBinLabel(int bin, const std::string &label, int axis=1)
set bin label for x, y or z axis (axis=1, 2, 3 respectively)
MonitorElement * TauSpinEffectsZ_Xf
int tauProngs(const HepMC::GenParticle *, double weight)
MonitorElement * TauBremPhotonsPtSum
void spinEffectsZ(const HepMC::GenParticle *boson, double weight)
MonitorElement * TauSpinEffectsHpm_UpsilonRho
virtual void analyze(const edm::Event &, const edm::EventSetup &)
void getData(T &iHolder) const
MonitorElement * TauSpinEffectsZ_eX
int findMother(const HepMC::GenParticle *)
MonitorElement * TauProngs
int tauMother(const HepMC::GenParticle *, double weight)
int tauDecayChannel(const HepMC::GenParticle *, double weight=0.0)
bool isTauFinalStateParticle(int pdgid)
const std::vector< HepMC::GenParticle * > GetMothers(const HepMC::GenParticle *boson)
math::XYZTLorentzVectorD LV
MonitorElement * TauSpinEffectsW_UpsilonRho
void photons(const HepMC::GenParticle *, double weight)
double leadingPionMomentum(const HepMC::GenParticle *, double weight)
MonitorElement * TauSpinEffectsH_Xf
MonitorElement * TauFSRPhotonsPt
MonitorElement * TauSpinEffectsZ_muX
MonitorElement * TauDecayChannels
edm::InputTag hepmcCollection_
MonitorElement * TauSpinEffectsH_Zs
MonitorElement * TauFSRPhotonsPtSum
HepPDT::ParticleData ParticleData
double visibleTauEnergy(const HepMC::GenParticle *)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
TLorentzVector leadingPionP4(const HepMC::GenParticle *)
bool isLastTauinChain(const HepMC::GenParticle *tau)
MonitorElement * TauSpinEffectsZ_MVis
TLorentzVector motherP4(const HepMC::GenParticle *)
void rtau(const HepMC::GenParticle *, int, int, double weight)
void FindPhotosFSR(const HepMC::GenParticle *p, std::vector< const HepMC::GenParticle * > &ListofFSR, double &BosonScale)
virtual void endRun(const edm::Run &, const edm::EventSetup &)
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
PDT table.
MonitorElement * TauSpinEffectsH_eX
MonitorElement * TauSpinEffectsW_eX
std::vector< std::vector< double > > tmp
const HepMC::GenParticle * GetMother(const HepMC::GenParticle *tau)
MonitorElement * TauSpinEffectsH_muX
MonitorElement * TauSpinEffectsW_UpsilonA1
DQMStore * dbe
ME's "container".
void findFSRandBrem(const HepMC::GenParticle *p, bool doBrem, std::vector< const HepMC::GenParticle * > &ListofFSR, std::vector< const HepMC::GenParticle * > &ListofBrem)
MonitorElement * TauSpinEffectsZ_Zs
MonitorElement * TauSpinEffectsHpm_X
virtual void beginRun(const edm::Run &, const edm::EventSetup &)
MonitorElement * TauSpinEffectsHpm_muX
double weight(const edm::Event &)
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
MonitorElement * TauSpinEffectsH_Xb
std::vector< HepMC::GenParticle * > Get_TauDecayProducts()
bool AnalyzeTau(HepMC::GenParticle *Tau, unsigned int &JAK_ID, unsigned int &TauBitMask, bool dores=true, bool dopi0=true)
MonitorElement * nPrimeTaus
MonitorElement * TauBremPhotonsPt
MonitorElement * TauRtauHpm
void setCurrentFolder(const std::string &fullpath)
unsigned int nProng(unsigned int &TauBitMask)