7 #include "CLHEP/Units/defs.h"
8 #include "CLHEP/Units/PhysicalConstants.h"
19 wmanager_(iPSet,consumesCollector())
20 ,genparticleCollection_(iPSet.getParameter<edm::
InputTag>(
"genparticleCollection"))
21 ,hepmcCollection_(iPSet.getParameter<edm::
InputTag>(
"hepmcCollection"))
42 nTaus = i.
book1D(
"nTaus",
"n analyzed Taus", 1, 0., 1.);
132 TauFSRPhotonsN=i.
book1D(
"TauFSRPhotonsN",
"FSR Photons radiating from/with tau (Gauge Boson)", 5 ,-0.5,4.5);
134 TauFSRPhotonsPt=i.
book1D(
"TauFSRPhotonsPt",
"Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
136 TauFSRPhotonsPtSum=i.
book1D(
"TauFSRPhotonsPtSum",
"Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
148 JAKInvMass.push_back(std::vector<MonitorElement *>());
155 JAKInvMass.at(
j).push_back(i.
book1D(
"M13"+tmp,
"M_{13,"+tmp+
"} (GeV)", 80 ,0,2.0));
156 JAKInvMass.at(
j).push_back(i.
book1D(
"M23"+tmp,
"M_{23,"+tmp+
"} (GeV)", 80 ,0,2.0));
157 JAKInvMass.at(
j).push_back(i.
book1D(
"M12"+tmp,
"M_{12,"+tmp+
"} (GeV)", 80 ,0,2.0));
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++){
209 SV=TVector3(part.at(
i)->vx(),part.at(
i)->vy(),part.at(
i)->vz());
212 double c(2.99792458E8),Ltau(DL.Mag()/100),
beta(
iter->p()/
iter->mass());
218 TLorentzVector
LV(part.at(
i)->px(),part.at(
i)->py(),part.at(
i)->pz(),part.at(
i)->energy());
272 std::vector<const reco::GenParticle*> mothers;
276 mothers.push_back(mother);
293 TauList.insert(TauList.begin(),
tau);
303 std::vector<const reco::GenParticle*> &ListofBrem){
312 if(
abs((dau)->
pdgId()) ==
abs(photo_ID) && !doBrem){ListofFSR.push_back(dau);}
313 if(
abs((dau)->
pdgId()) ==
abs(photo_ID) && doBrem){ListofBrem.push_back(dau);}
323 int mother_pid=m->
pdgId();
328 ListofFSR.push_back(dau);
332 if(
abs(mother_pid) == 24) BosonScale=1.0;
333 if(
abs(mother_pid) == 23) BosonScale=2.0;
334 if(
abs(mother_pid) == 22) BosonScale=2.0;
335 if(
abs(mother_pid) == 25) BosonScale=2.0;
336 if(
abs(mother_pid) == 35) BosonScale=2.0;
337 if(
abs(mother_pid) == 36) BosonScale=2.0;
338 if(
abs(mother_pid) == 37) BosonScale=1.0;
342 if(
abs(tau->
pdgId()) != 15 )
return -3;
344 if(mother_pid == -2)
return -2;
346 if(
abs(mother_pid) == 24) label =
W;
347 if(
abs(mother_pid) == 23) label =
Z;
348 if(
abs(mother_pid) == 22) label =
gamma;
349 if(
abs(mother_pid) == 25) label =
HSM;
350 if(
abs(mother_pid) == 35) label =
H0;
351 if(
abs(mother_pid) == 36) label =
A0;
352 if(
abs(mother_pid) == 37) label =
Hpm;
353 int mother_shortpid=(
abs(mother_pid)%10000);
354 if(mother_shortpid>500 && mother_shortpid<600 )label =
B;
355 if(mother_shortpid>400 && mother_shortpid<500)label =
D;
357 if(label==
B || label ==
D || label ==
other)
return -1;
374 countParticles(tau,allCount,eCount,muCount,pi0Count,piCount,rhoCount,a1Count,KCount,KstarCount);
377 if(KCount >= 1) channel =
K;
378 if(KstarCount >= 1) channel =
Kstar;
379 if(a1Count >= 1) channel =
a1;
380 if(rhoCount >= 1) channel =
rho;
384 if(piCount == 1 && pi0Count == 0) channel =
pi;
385 if(piCount == 1 && pi0Count == 1) channel =
pi1pi0;
386 if(piCount == 1 && pi0Count > 1) channel =
pinpi0;
387 if(piCount == 3 && pi0Count == 0) channel =
tripi;
388 if(piCount == 3 && pi0Count > 0) channel =
tripinpi0;
390 if(muCount == 1) channel =
muon;
396 int &pi0Count,
int &piCount,
int &rhoCount,
int &a1Count,
int &KCount,
int &KstarCount){
401 if(
abs(pid) == 11) eCount++;
402 if(
abs(pid) == 13) muCount++;
403 if(
abs(pid) == 111) pi0Count++;
404 if(
abs(pid) == 211) piCount++;
405 if(
abs(pid) == 213) rhoCount++;
406 if(
abs(pid) == 20213) a1Count++;
407 if(
abs(pid) == 321) KCount++;
408 if(
abs(pid) == 323) KstarCount++;
409 countParticles(dau,allCount,eCount,muCount,pi0Count,piCount,rhoCount,a1Count,KCount,KstarCount);
417 TLorentzVector momP4 =
motherP4(tau);
419 pionP4.Boost(-1*momP4.BoostVector());
420 double energy = pionP4.E()/(momP4.M()/2);
435 TLorentzVector
rho(0,0,0,0),
pi(0,0,0,0);
436 for(
unsigned int i=0;
i<part.size();
i++){
437 TLorentzVector
LV(part.at(
i)->px(),part.at(
i)->py(),part.at(
i)->pz(),part.at(
i)->energy());
445 TLorentzVector
a1(0,0,0,0),pi_p(0,0,0,0),pi_m(0,0,0,0);
446 int nplus(0),nminus(0);
447 for(
unsigned int i=0;
i<part.size();
i++){
448 TLorentzVector
LV(part.at(
i)->px(),part.at(
i)->py(),part.at(
i)->pz(),part.at(
i)->energy());
453 if(nplus+nminus==3 && nplus==1) gamma=2*pi_p.P()/
a1.P()-1;
454 if(nplus+nminus==3 && nminus==1) gamma=2*pi_m.P()/
a1.P()-1;
456 pi_p+=pi_m; gamma=2*pi_p.P()/
a1.P()-1;
467 if(ntau==1 && dau->
pdgId() == 15)
return;
472 TLorentzVector tautau(0,0,0,0);
473 TLorentzVector pipi(0,0,0,0);
474 TLorentzVector taum(0,0,0,0);
475 TLorentzVector taup(0,0,0,0);
476 TLorentzVector rho_plus,rho_minus,pi_rhominus,pi0_rhominus,pi_rhoplus,pi0_rhoplus,pi_plus,pi_minus;
477 bool hasrho_minus(
false),hasrho_plus(
false),haspi_minus(
false),haspi_plus(
false);
478 int nSinglePionDecays(0),nSingleMuonDecays(0),nSingleElectronDecays(0);
480 TLorentzVector Zboson(boson->
px(),boson->
py(),boson->
pz(),boson->
energy());
486 unsigned int jak_id, TauBitMask;
487 if(TD.
AnalyzeTau(dau,jak_id,TauBitMask,
false,
false)){
493 TLorentzVector LVtau(dau->
px(),dau->
py(),dau->
pz(),dau->
energy());
498 int charge = (int) pd->charge();
499 LVtau.Boost(-1*Zboson.BoostVector());
500 LVpi.Boost(-1*Zboson.BoostVector());
512 if(charge<0){x1=LVpi.P()/LVtau.E(); taum=LVtau;}
513 else{ x2=LVpi.P()/LVtau.E();}
515 TLorentzVector LVtau(dau->
px(),dau->
py(),dau->
pz(),dau->
energy());
516 if(pid == 15)taum=LVtau;
517 if(pid ==-15)taup=LVtau;
519 for(
unsigned int i=0;
i<part.size();
i++){
520 int pid_d = part.at(
i)->pdgId();
521 if(
abs(pid_d)==211 ||
abs(pid_d)==111){
522 TLorentzVector
LV(part.at(
i)->px(),part.at(
i)->py(),part.at(
i)->pz(),part.at(
i)->energy());
525 if(pid_d==-211 ){ pi_rhominus=
LV;}
526 if(
abs(pid_d)==111 ){ pi0_rhominus=
LV;}
530 if(pid_d==211 ){pi_rhoplus=
LV;}
531 if(
abs(pid_d)==111 ){pi0_rhoplus=
LV;}
537 for(
unsigned int i=0;
i<part.size();
i++){
538 int pid_d = part.at(
i)->pdgId();
539 if(
abs(pid_d)==211 ){
540 TLorentzVector
LV(part.at(
i)->px(),part.at(
i)->py(),part.at(
i)->pz(),part.at(
i)->energy());
543 if(pid_d==-211 ){ pi_minus=
LV;}
547 if(pid_d==211 ){pi_plus=
LV;}
555 if(hasrho_minus && hasrho_plus){
557 rho_minus=pi_rhominus;
558 rho_minus+=pi0_rhominus;
560 rho_plus+=pi0_rhoplus;
561 TLorentzVector rhorho=rho_minus;rhorho+=rho_plus;
564 TLorentzVector pi_rhoplusb=pi_rhoplus; pi_rhoplusb.Boost(-1*rhorho.BoostVector());
565 TLorentzVector pi0_rhoplusb=pi0_rhoplus; pi0_rhoplusb.Boost(-1*rhorho.BoostVector());
566 TLorentzVector pi_rhominusb=pi_rhominus; pi_rhominusb.Boost(-1*rhorho.BoostVector());
567 TLorentzVector pi0_rhominusb=pi0_rhominus; pi0_rhominusb.Boost(-1*rhorho.BoostVector());
570 TVector3 n_plus=pi_rhoplusb.Vect().Cross(pi0_rhoplusb.Vect());
571 TVector3 n_minus=pi_rhominusb.Vect().Cross(pi0_rhominusb.Vect());
574 double Acoplanarity=acos(n_plus.Dot(n_minus)/(n_plus.Mag()*n_minus.Mag()));
575 if(pi_rhominusb.Vect().Dot(n_plus)>0){Acoplanarity*=-1;Acoplanarity+=2*
TMath::Pi();}
578 pi_rhoplus.Boost(-1*taup.BoostVector());
579 pi0_rhoplus.Boost(-1*taup.BoostVector());
580 pi_rhominus.Boost(-1*taum.BoostVector());
581 pi0_rhominus.Boost(-1*taum.BoostVector());
584 double y1=(pi_rhoplus.E()-pi0_rhoplus.E())/(pi_rhoplus.E()+pi0_rhoplus.E());
585 double y2=(pi_rhominus.E()-pi0_rhominus.E())/(pi_rhominus.E()+pi0_rhominus.E());
591 if(haspi_minus && haspi_plus){
592 TLorentzVector tauporig=taup;
593 TLorentzVector taumorig=taum;
596 pi_plus.Boost(-1*Zboson.BoostVector());
597 pi_minus.Boost(-1*Zboson.BoostVector());
599 taup.Boost(-1*Zboson.BoostVector());
600 taum.Boost(-1*Zboson.BoostVector());
607 double proj_m=taum.Vect().Dot(pi_minus.Vect())/(taum.P()*taum.P());
608 double proj_p=taup.Vect().Dot(pi_plus.Vect())/(taup.P()*taup.P());
609 TVector3 Tau_m=taum.Vect();
610 TVector3 Tau_p=taup.Vect();
613 TVector3 Pit_m=pi_minus.Vect()-Tau_m;
614 TVector3 Pit_p=pi_plus.Vect()-Tau_p;
616 double Acoplanarity=acos(Pit_m.Dot(Pit_p)/(Pit_p.Mag()*Pit_m.Mag()));
617 TVector3
n=Pit_p.Cross(Pit_m);
618 if(n.Dot(Tau_m)/Tau_m.Mag()>0){Acoplanarity*=-1; Acoplanarity+=2*
TMath::Pi();}
625 if(nSingleMuonDecays==2){
629 if(nSingleElectronDecays==2){
633 if(nSinglePionDecays == 2 && tautau.M()!= 0) {
638 if(x2-x1>alow && x2-x1<aup){
639 double zs=(zsup+zslow)/2;
649 const std::vector<const reco::GenParticle*>
m=
GetMothers(boson);
651 TLorentzVector
Z(0,0,0,0);
652 for(
unsigned int i=0;
i<m.size();
i++){
657 if(taum.Vect().Dot(Zboson.Vect())/(Zboson.P()*taum.P())>0){
672 double a=1-
sqrt(fabs(1.0-2*fabs(zs)));
685 TLorentzVector
p4(0,0,0,0);
690 if(!(
abs(pid)==211 ||
abs(pid)==13 ||
abs(pid)==11))
continue;
691 if(dau->
p() > p4.P()) p4 = TLorentzVector(dau->
px(),dau->
py(),dau->
pz(),dau->
energy());
698 return TLorentzVector(m->
px(),m->
py(),m->
pz(),m->
energy());
707 if(
abs(pid) == 12 ||
abs(pid) == 14 ||
abs(pid) == 16) {
716 std::vector<const reco::GenParticle*> TauList;
721 std::vector<const reco::GenParticle*> ListofFSR; ListofFSR.clear();
722 std::vector<const reco::GenParticle*> ListofBrem; ListofBrem.clear();
723 std::vector<const reco::GenParticle*> FSR_photos; FSR_photos.clear();
724 double BosonScale(1);
725 if(TauList.size()>0){
731 double photonPtSum=0;
732 for(
unsigned int i=0;
i<ListofBrem.size();
i++){
733 photonPtSum+=ListofBrem.at(
i)->pt();
742 for(
unsigned int i=0;
i<ListofFSR.size();
i++){
743 photonPtSum+=ListofFSR.at(
i)->pt();
746 double FSR_photosSum(0);
747 for(
unsigned int i=0;
i<FSR_photos.size();
i++){
748 FSR_photosSum+=FSR_photos.at(
i)->pt();
MonitorElement * TauFSRPhotonsN
virtual double energy() const GCC11_FINAL
energy
int findMother(const reco::GenParticle *)
MonitorElement * TauSpinEffectsH_rhorhoAcoplanarityplus
MonitorElement * TauSpinEffectsW_X
MonitorElement * TauSpinEffectsHpm_UpsilonA1
std::vector< std::vector< MonitorElement * > > JAKInvMass
MonitorElement * TauSpinEffectsZ_Xb
virtual double p() const GCC11_FINAL
magnitude of momentum vector
MonitorElement * DecayLength
bool getByToken(EDGetToken token, Handle< PROD > &result) const
virtual int pdgId() const GCC11_FINAL
PDG identifier.
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
virtual double pz() const GCC11_FINAL
z coordinate of momentum vector
virtual double py() const GCC11_FINAL
y coordinate of momentum vector
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)
MonitorElement * TauSpinEffectsW_UpsilonRho
const reco::GenParticle * GetMother(const reco::GenParticle *tau)
MonitorElement * TauSpinEffectsZ_X88to100
MonitorElement * TauSpinEffectsH_rhorhoAcoplanarityminus
virtual int status() const GCC11_FINAL
status word
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
virtual double px() const GCC11_FINAL
x coordinate of momentum vector
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 *)
int tauDecayChannel(const reco::GenParticle *tau, int jak_id, unsigned int TauBitMask, double weight)
void setCurrentFolder(const std::string &fullpath)
MonitorElement * TauSpinEffectsH_pipiAcollinearity
MonitorElement * TauSpinEffectsZ_X
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
PDT table.
bool AnalyzeTau(const reco::GenParticle *Tau, unsigned int &JAK_ID, unsigned int &TauBitMask, bool dores, bool dopi0)
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 &)
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
void spinEffectsZH(const reco::GenParticle *boson, double weight)
unsigned int nProng(unsigned int &TauBitMask)
MonitorElement * TauSpinEffectsH_X
TLorentzVector leadingPionP4(const reco::GenParticle *)