7 #include "CLHEP/Units/defs.h"
8 #include "CLHEP/Units/PhysicalConstants.h"
18 wmanager_(iPSet,consumesCollector())
19 ,genparticleCollection_(iPSet.getParameter<edm::
InputTag>(
"genparticleCollection"))
20 ,hepmcCollection_(iPSet.getParameter<edm::
InputTag>(
"hepmcCollection"))
41 nTaus = i.
book1D(
"nTaus",
"n analyzed Taus", 1, 0., 1.);
131 TauFSRPhotonsN=i.
book1D(
"TauFSRPhotonsN",
"FSR Photons radiating from/with tau (Gauge Boson)", 5 ,-0.5,4.5);
133 TauFSRPhotonsPt=i.
book1D(
"TauFSRPhotonsPt",
"Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
135 TauFSRPhotonsPtSum=i.
book1D(
"TauFSRPhotonsPtSum",
"Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
147 MODEInvMass.push_back(std::vector<MonitorElement *>());
172 for (reco::GenParticleCollection::const_iterator
iter = genParticles->begin();
iter != genParticles->end(); ++
iter) {
189 unsigned int jak_id, TauBitMask;
198 TLorentzVector LVQ(0,0,0,0);
199 TLorentzVector LVS12(0,0,0,0);
200 TLorentzVector LVS13(0,0,0,0);
201 TLorentzVector LVS23(0,0,0,0);
205 for(
unsigned int i=0;
i<part.size();
i++){
210 SV=TVector3(part.at(
i)->vx(),part.at(
i)->vy(),part.at(
i)->vz());
213 double c(2.99792458E8),Ltau(DL.Mag()/100),
beta(
iter->p()/
iter->mass());
219 TLorentzVector
LV(part.at(
i)->px(),part.at(
i)->py(),part.at(
i)->pz(),part.at(
i)->energy());
273 std::vector<const reco::GenParticle*> mothers;
277 mothers.push_back(mother);
294 TauList.insert(TauList.begin(),
tau);
304 std::vector<const reco::GenParticle*> &ListofBrem){
313 if(
abs((dau)->
pdgId()) ==
abs(photo_ID) && !doBrem){ListofFSR.push_back(dau);}
314 if(
abs((dau)->
pdgId()) ==
abs(photo_ID) && doBrem){ListofBrem.push_back(dau);}
324 int mother_pid=m->
pdgId();
329 ListofFSR.push_back(dau);
333 if(
abs(mother_pid) == 24) BosonScale=1.0;
334 if(
abs(mother_pid) == 23) BosonScale=2.0;
335 if(
abs(mother_pid) == 22) BosonScale=2.0;
336 if(
abs(mother_pid) == 25) BosonScale=2.0;
337 if(
abs(mother_pid) == 35) BosonScale=2.0;
338 if(
abs(mother_pid) == 36) BosonScale=2.0;
339 if(
abs(mother_pid) == 37) BosonScale=1.0;
343 if(
abs(tau->
pdgId()) != 15 )
return -3;
345 if(mother_pid == -2)
return -2;
347 if(
abs(mother_pid) == 24) label =
W;
348 if(
abs(mother_pid) == 23) label =
Z;
349 if(
abs(mother_pid) == 22) label =
gamma;
350 if(
abs(mother_pid) == 25) label =
HSM;
351 if(
abs(mother_pid) == 35) label =
H0;
352 if(
abs(mother_pid) == 36) label =
A0;
353 if(
abs(mother_pid) == 37) label =
Hpm;
354 int mother_shortpid=(
abs(mother_pid)%10000);
355 if(mother_shortpid>500 && mother_shortpid<600 )label =
B;
356 if(mother_shortpid>400 && mother_shortpid<500)label =
D;
358 if(label==
B || label ==
D || label ==
other)
return -1;
375 countParticles(tau,allCount,eCount,muCount,pi0Count,piCount,rhoCount,a1Count,KCount,KstarCount);
378 if(KCount >= 1) channel =
K;
379 if(KstarCount >= 1) channel =
Kstar;
380 if(a1Count >= 1) channel =
a1;
381 if(rhoCount >= 1) channel =
rho;
385 if(piCount == 1 && pi0Count == 0) channel =
pi;
386 if(piCount == 1 && pi0Count == 1) channel =
pi1pi0;
387 if(piCount == 1 && pi0Count > 1) channel =
pinpi0;
388 if(piCount == 3 && pi0Count == 0) channel =
tripi;
389 if(piCount == 3 && pi0Count > 0) channel =
tripinpi0;
391 if(muCount == 1) channel =
muon;
397 int &pi0Count,
int &piCount,
int &rhoCount,
int &a1Count,
int &KCount,
int &KstarCount){
402 if(
abs(pid) == 11) eCount++;
403 if(
abs(pid) == 13) muCount++;
404 if(
abs(pid) == 111) pi0Count++;
405 if(
abs(pid) == 211) piCount++;
406 if(
abs(pid) == 213) rhoCount++;
407 if(
abs(pid) == 20213) a1Count++;
408 if(
abs(pid) == 321) KCount++;
409 if(
abs(pid) == 323) KstarCount++;
410 countParticles(dau,allCount,eCount,muCount,pi0Count,piCount,rhoCount,a1Count,KCount,KstarCount);
418 TLorentzVector momP4 =
motherP4(tau);
420 pionP4.Boost(-1*momP4.BoostVector());
421 double energy = pionP4.E()/(momP4.M()/2);
436 TLorentzVector
rho(0,0,0,0),
pi(0,0,0,0);
437 for(
unsigned int i=0;
i<part.size();
i++){
438 TLorentzVector
LV(part.at(
i)->px(),part.at(
i)->py(),part.at(
i)->pz(),part.at(
i)->energy());
446 TLorentzVector
a1(0,0,0,0),pi_p(0,0,0,0),pi_m(0,0,0,0);
447 int nplus(0),nminus(0);
448 for(
unsigned int i=0;
i<part.size();
i++){
449 TLorentzVector
LV(part.at(
i)->px(),part.at(
i)->py(),part.at(
i)->pz(),part.at(
i)->energy());
454 if(nplus+nminus==3 && nplus==1) gamma=2*pi_p.P()/
a1.P()-1;
455 if(nplus+nminus==3 && nminus==1) gamma=2*pi_m.P()/
a1.P()-1;
457 pi_p+=pi_m; gamma=2*pi_p.P()/
a1.P()-1;
468 if(ntau==1 && dau->
pdgId() == 15)
return;
473 TLorentzVector tautau(0,0,0,0);
474 TLorentzVector pipi(0,0,0,0);
475 TLorentzVector taum(0,0,0,0);
476 TLorentzVector taup(0,0,0,0);
477 TLorentzVector rho_plus,rho_minus,pi_rhominus,pi0_rhominus,pi_rhoplus,pi0_rhoplus,pi_plus,pi_minus;
478 bool hasrho_minus(
false),hasrho_plus(
false),haspi_minus(
false),haspi_plus(
false);
479 int nSinglePionDecays(0),nSingleMuonDecays(0),nSingleElectronDecays(0);
481 TLorentzVector Zboson(boson->
px(),boson->
py(),boson->
pz(),boson->
energy());
487 unsigned int jak_id, TauBitMask;
488 if(TD.
AnalyzeTau(dau,jak_id,TauBitMask,
false,
false)){
494 TLorentzVector LVtau(dau->
px(),dau->
py(),dau->
pz(),dau->
energy());
499 int charge = (int) pd->charge();
500 LVtau.Boost(-1*Zboson.BoostVector());
501 LVpi.Boost(-1*Zboson.BoostVector());
527 if(charge<0){x1=LVpi.P()/LVtau.E(); taum=LVtau;}
528 else{ x2=LVpi.P()/LVtau.E();}
530 TLorentzVector LVtau(dau->
px(),dau->
py(),dau->
pz(),dau->
energy());
531 if(pid == 15)taum=LVtau;
532 if(pid ==-15)taup=LVtau;
534 for(
unsigned int i=0;
i<part.size();
i++){
535 int pid_d = part.at(
i)->pdgId();
536 if(
abs(pid_d)==211 ||
abs(pid_d)==111){
537 TLorentzVector
LV(part.at(
i)->px(),part.at(
i)->py(),part.at(
i)->pz(),part.at(
i)->energy());
540 if(pid_d==-211 ){ pi_rhominus=
LV;}
541 if(
abs(pid_d)==111 ){ pi0_rhominus=
LV;}
545 if(pid_d==211 ){pi_rhoplus=
LV;}
546 if(
abs(pid_d)==111 ){pi0_rhoplus=
LV;}
552 for(
unsigned int i=0;
i<part.size();
i++){
553 int pid_d = part.at(
i)->pdgId();
554 if(
abs(pid_d)==211 ){
555 TLorentzVector
LV(part.at(
i)->px(),part.at(
i)->py(),part.at(
i)->pz(),part.at(
i)->energy());
558 if(pid_d==-211 ){ pi_minus=
LV;}
562 if(pid_d==211 ){pi_plus=
LV;}
570 if(hasrho_minus && hasrho_plus){
572 rho_minus=pi_rhominus;
573 rho_minus+=pi0_rhominus;
575 rho_plus+=pi0_rhoplus;
576 TLorentzVector rhorho=rho_minus;rhorho+=rho_plus;
579 TLorentzVector pi_rhoplusb=pi_rhoplus; pi_rhoplusb.Boost(-1*rhorho.BoostVector());
580 TLorentzVector pi0_rhoplusb=pi0_rhoplus; pi0_rhoplusb.Boost(-1*rhorho.BoostVector());
581 TLorentzVector pi_rhominusb=pi_rhominus; pi_rhominusb.Boost(-1*rhorho.BoostVector());
582 TLorentzVector pi0_rhominusb=pi0_rhominus; pi0_rhominusb.Boost(-1*rhorho.BoostVector());
585 TVector3 n_plus=pi_rhoplusb.Vect().Cross(pi0_rhoplusb.Vect());
586 TVector3 n_minus=pi_rhominusb.Vect().Cross(pi0_rhominusb.Vect());
589 double Acoplanarity=acos(n_plus.Dot(n_minus)/(n_plus.Mag()*n_minus.Mag()));
590 if(pi_rhominusb.Vect().Dot(n_plus)>0){Acoplanarity*=-1;Acoplanarity+=2*
TMath::Pi();}
593 pi_rhoplus.Boost(-1*taup.BoostVector());
594 pi0_rhoplus.Boost(-1*taup.BoostVector());
595 pi_rhominus.Boost(-1*taum.BoostVector());
596 pi0_rhominus.Boost(-1*taum.BoostVector());
599 double y1=(pi_rhoplus.E()-pi0_rhoplus.E())/(pi_rhoplus.E()+pi0_rhoplus.E());
600 double y2=(pi_rhominus.E()-pi0_rhominus.E())/(pi_rhominus.E()+pi0_rhominus.E());
606 if(haspi_minus && haspi_plus){
607 TLorentzVector tauporig=taup;
608 TLorentzVector taumorig=taum;
611 pi_plus.Boost(-1*Zboson.BoostVector());
612 pi_minus.Boost(-1*Zboson.BoostVector());
614 taup.Boost(-1*Zboson.BoostVector());
615 taum.Boost(-1*Zboson.BoostVector());
622 double proj_m=taum.Vect().Dot(pi_minus.Vect())/(taum.P()*taum.P());
623 double proj_p=taup.Vect().Dot(pi_plus.Vect())/(taup.P()*taup.P());
624 TVector3 Tau_m=taum.Vect();
625 TVector3 Tau_p=taup.Vect();
628 TVector3 Pit_m=pi_minus.Vect()-Tau_m;
629 TVector3 Pit_p=pi_plus.Vect()-Tau_p;
631 double Acoplanarity=acos(Pit_m.Dot(Pit_p)/(Pit_p.Mag()*Pit_m.Mag()));
632 TVector3
n=Pit_p.Cross(Pit_m);
633 if(n.Dot(Tau_m)/Tau_m.Mag()>0){Acoplanarity*=-1; Acoplanarity+=2*
TMath::Pi();}
639 if(nSinglePionDecays == 2 && tautau.M()!= 0) {
644 if(x2-x1>alow && x2-x1<aup){
645 double zs=(zsup+zslow)/2;
655 const std::vector<const reco::GenParticle*>
m=
GetMothers(boson);
657 TLorentzVector
Z(0,0,0,0);
658 for(
unsigned int i=0;
i<m.size();
i++){
663 if(taum.Vect().Dot(Zboson.Vect())/(Zboson.P()*taum.P())>0){
678 double a=1-
sqrt(fabs(1.0-2*fabs(zs)));
691 TLorentzVector
p4(0,0,0,0);
696 if(!(
abs(pid)==211 ||
abs(pid)==13 ||
abs(pid)==11))
continue;
697 if(dau->
p() > p4.P()) p4 = TLorentzVector(dau->
px(),dau->
py(),dau->
pz(),dau->
energy());
704 return TLorentzVector(m->
px(),m->
py(),m->
pz(),m->
energy());
713 if(
abs(pid) == 12 ||
abs(pid) == 14 ||
abs(pid) == 16) {
722 std::vector<const reco::GenParticle*> TauList;
727 std::vector<const reco::GenParticle*> ListofFSR; ListofFSR.clear();
728 std::vector<const reco::GenParticle*> ListofBrem; ListofBrem.clear();
729 std::vector<const reco::GenParticle*> FSR_photos; FSR_photos.clear();
730 double BosonScale(1);
731 if(TauList.size()>0){
737 double photonPtSum=0;
738 for(
unsigned int i=0;
i<ListofBrem.size();
i++){
739 photonPtSum+=ListofBrem.at(
i)->pt();
748 for(
unsigned int i=0;
i<ListofFSR.size();
i++){
749 photonPtSum+=ListofFSR.at(
i)->pt();
752 double FSR_photosSum(0);
753 for(
unsigned int i=0;
i<FSR_photos.size();
i++){
754 FSR_photosSum+=FSR_photos.at(
i)->pt();
MonitorElement * TauFSRPhotonsN
virtual int pdgId() const
PDG identifier.
virtual double p() const
magnitude of momentum vector
int findMother(const reco::GenParticle *)
MonitorElement * TauSpinEffectsH_rhorhoAcoplanarityplus
MonitorElement * TauSpinEffectsW_X
MonitorElement * TauSpinEffectsHpm_UpsilonA1
MonitorElement * TauSpinEffectsZ_Xb
MonitorElement * DecayLength
std::vector< std::vector< MonitorElement * > > MODEInvMass
bool getByToken(EDGetToken token, Handle< PROD > &result) const
virtual int status() const
status word
TauValidation(const edm::ParameterSet &)
MonitorElement * TauSpinEffectsH_MVis
MonitorElement * TauMothers
MonitorElement * TauSpinEffectsHpm_eX
MonitorElement * TauSpinEffectsZ_X100to120
void FindPhotosFSR(const reco::GenParticle *p, std::vector< const reco::GenParticle * > &ListofFSR, double &BosonScale)
MonitorElement * TauSpinEffectsW_muX
MonitorElement * TauBremPhotonsN
void spinEffectsWHpm(const reco::GenParticle *, int, int, std::vector< const reco::GenParticle * > &part, double weight)
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
bool isLastTauinChain(const reco::GenParticle *tau)
edm::EDGetTokenT< edm::HepMCProduct > hepmcCollectionToken_
MonitorElement * TauBremPhotonsPtSum
MonitorElement * TauSpinEffectsH_pipiAcoplanarity
MonitorElement * TauSpinEffectsHpm_UpsilonRho
double leadingPionMomentum(const reco::GenParticle *, double weight)
void getData(T &iHolder) const
MonitorElement * TauSpinEffectsZ_eX
MonitorElement * TauProngs
TLorentzVector motherP4(const reco::GenParticle *)
bool isTauFinalStateParticle(int pdgid)
virtual void dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) override
MonitorElement * TauSpinEffectsH_pipiAcollinearityzoom
math::XYZTLorentzVectorD LV
void photons(const reco::GenParticle *, double weight)
virtual double energy() const
energy
MonitorElement * TauSpinEffectsW_UpsilonRho
const reco::GenParticle * GetMother(const reco::GenParticle *tau)
MonitorElement * TauSpinEffectsZ_X88to100
MonitorElement * TauSpinEffectsH_rhorhoAcoplanarityminus
virtual void analyze(edm::Event const &, edm::EventSetup const &) override
virtual size_t numberOfMothers() const
number of mothers
MonitorElement * TauSpinEffectsZ_X50to75
virtual size_t numberOfDaughters() const
number of daughters
virtual const Candidate * daughter(size_type) const
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
MonitorElement * TauSpinEffectsH_Xf
MonitorElement * book1D(Args &&...args)
MonitorElement * TauFSRPhotonsPt
Abs< T >::type abs(const T &t)
MonitorElement * TauSpinEffectsZ_muX
MonitorElement * TauDecayChannels
edm::InputTag hepmcCollection_
MonitorElement * TauSpinEffectsH_Zs
MonitorElement * TauFSRPhotonsPtSum
HepPDT::ParticleData ParticleData
int tauMother(const reco::GenParticle *, double weight)
virtual void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override
std::vector< const reco::GenParticle * > Get_TauDecayProducts()
const std::vector< const reco::GenParticle * > GetMothers(const reco::GenParticle *boson)
MonitorElement * TauSpinEffectsZ_MVis
virtual int pdgId() const =0
PDG identifier.
MonitorElement * TauSpinEffectsZ_X75to88
double visibleTauEnergy(const reco::GenParticle *)
static std::string DecayMode(unsigned int &MODE_ID)
int tauDecayChannel(const reco::GenParticle *tau, int jak_id, unsigned int TauBitMask, double weight)
void setCurrentFolder(const std::string &fullpath)
MonitorElement * TauSpinEffectsH_pipiAcollinearity
virtual double px() const
x coordinate of momentum vector
MonitorElement * TauSpinEffectsZ_X
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
PDT table.
virtual double pz() const
z coordinate of momentum vector
MonitorElement * TauSpinEffectsH_eX
MonitorElement * TauSpinEffectsW_eX
std::vector< std::vector< double > > tmp
edm::EDGetTokenT< reco::GenParticleCollection > genparticleCollectionToken_
void findTauList(const reco::GenParticle *tau, std::vector< const reco::GenParticle * > &TauList)
void findFSRandBrem(const reco::GenParticle *p, bool doBrem, std::vector< const reco::GenParticle * > &ListofFSR, std::vector< const reco::GenParticle * > &ListofBrem)
MonitorElement * LifeTime
MonitorElement * TauSpinEffectsH_muX
MonitorElement * TauSpinEffectsW_UpsilonA1
void countParticles(const reco::GenParticle *p, int &allCount, int &eCount, int &muCount, int &pi0Count, int &piCount, int &rhoCount, int &a1Count, int &KCount, int &KstarCount)
edm::InputTag genparticleCollection_
MonitorElement * TauSpinEffectsZ_Zs
MonitorElement * TauSpinEffectsHpm_X
MonitorElement * TauSpinEffectsHpm_muX
double weight(const edm::Event &)
bool AnalyzeTau(const reco::GenParticle *Tau, unsigned int &MODE_ID, unsigned int &TauBitMask, bool dores, bool dopi0)
MonitorElement * TauSpinEffectsZ_X120UP
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
MonitorElement * TauSpinEffectsH_Xb
virtual const Candidate * mother(size_type=0) const
return mother at a given position, i = 0, ... numberOfMothers() - 1 (read only mode) ...
MonitorElement * nPrimeTaus
MonitorElement * TauBremPhotonsPt
virtual double py() const
y coordinate of momentum vector
void spinEffectsZH(const reco::GenParticle *boson, double weight)
unsigned int nProng(unsigned int &TauBitMask)
MonitorElement * TauSpinEffectsH_X
TLorentzVector leadingPionP4(const reco::GenParticle *)