CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
TauValidation.cc
Go to the documentation of this file.
1 /*class TauValidation
2  *
3  * Class to fill dqm monitor elements from existing EDM file
4  *
5  */
6 #include "Validation/EventGenerator/interface/TauValidation.h"
7 #include "CLHEP/Units/defs.h"
8 #include "CLHEP/Units/PhysicalConstants.h"
9 #include "DataFormats/Math/interface/LorentzVector.h"
10 #include "Validation/EventGenerator/interface/TauDecay_GenParticle.h"
11 #include "Validation/EventGenerator/interface/PdtPdgMini.h"
12 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
13 #include "DataFormats/HepMCCandidate/interface/GenParticleFwd.h"
14 #include <iostream>
15 #include "Validation/EventGenerator/interface/DQMHelper.h"
16 using namespace edm;
17 
18 TauValidation::TauValidation(const edm::ParameterSet& iPSet):
19  // wmanager_(iPSet,consumesCollector())
20  genparticleCollection_(iPSet.getParameter<edm::InputTag>("genparticleCollection"))
21  ,NMODEID(TauDecay::NMODEID-1)// fortran to C++ index
22  ,zsbins(20)
23  ,zsmin(-0.5)
24  ,zsmax(0.5)
25 {
26  genparticleCollectionToken_=consumes<reco::GenParticleCollection>(genparticleCollection_);
27 }
28 
29 TauValidation::~TauValidation(){}
30 
31 void TauValidation::dqmBeginRun(const edm::Run& r, const edm::EventSetup& c) {
32  c.getData( fPDGTable );
33 }
34 
35 void TauValidation::bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &){
37  DQMHelper dqm(&i); i.setCurrentFolder("Generator/Tau");
38  // Number of analyzed events
39  nTaus = dqm.book1dHisto("nTaus", "n analyzed Taus", 1, 0., 1.,"bin","Number of #tau's found");
40  nPrimeTaus = dqm.book1dHisto("nPrimeTaus", "n analyzed prime Taus", 1, 0., 1.,"bin","Number of #tau's from Gauge Bosons");
41 
42  //Kinematics
43  TauPt = dqm.book1dHisto("TauPt","Tau pT", 100 ,0,100,"P_{T}^{#tau}","Number of #tau's from Gauge Bosons");
44  TauEta = dqm.book1dHisto("TauEta","Tau eta", 100 ,-2.5,2.5,"#eta^{#tau}","Number of #tau's from Gauge Bosons");
45  TauPhi = dqm.book1dHisto("TauPhi","Tau phi", 100 ,-3.14,3.14,"#phi^{#tau}","Number of #tau's from Gauge Bosons");
46  TauProngs = dqm.book1dHisto("TauProngs","Tau n prongs", 7 ,0,7,"N_{prongs}","Number of #tau's from Gauge Bosons");
47  TauDecayChannels = dqm.book1dHisto("TauDecayChannels","Tau decay channels", 13 ,0,13,"Tau POG Decay Mode","Number of #tau's from Gauge Bosons");
48  TauDecayChannels->setBinLabel(1+undetermined,"?");
49  TauDecayChannels->setBinLabel(1+electron,"e");
50  TauDecayChannels->setBinLabel(1+muon,"mu");
51  TauDecayChannels->setBinLabel(1+pi,"#pi^{#pm}");
52  TauDecayChannels->setBinLabel(1+rho,"#rho^{#pm}");
53  TauDecayChannels->setBinLabel(1+a1,"a_{1}^{#pm}");
54  TauDecayChannels->setBinLabel(1+pi1pi0,"#pi^{#pm}#pi^{0}");
55  TauDecayChannels->setBinLabel(1+pinpi0,"#pi^{#pm}n#pi^{0}");
56  TauDecayChannels->setBinLabel(1+tripi,"3#pi^{#pm}");
57  TauDecayChannels->setBinLabel(1+tripinpi0,"3#pi^{#pm}n#pi^{0}");
58  TauDecayChannels->setBinLabel(1+K,"K");
59  TauDecayChannels->setBinLabel(1+Kstar,"K^{*}");
60  TauDecayChannels->setBinLabel(1+stable,"Stable");
61 
62  TauMothers = dqm.book1dHisto("TauMothers","Tau mother particles", 10 ,0,10,"Mother of #tau","Number of #tau's");
63 
64  TauMothers->setBinLabel(1+other,"?");
65  TauMothers->setBinLabel(1+B,"B Decays");
66  TauMothers->setBinLabel(1+D,"D Decays");
67  TauMothers->setBinLabel(1+gamma,"#gamma");
68  TauMothers->setBinLabel(1+Z,"Z");
69  TauMothers->setBinLabel(1+W,"W");
70  TauMothers->setBinLabel(1+HSM,"H_{SM}/h^{0}");
71  TauMothers->setBinLabel(1+H0,"H^{0}");
72  TauMothers->setBinLabel(1+A0,"A^{0}");
73  TauMothers->setBinLabel(1+Hpm,"H^{#pm}");
74 
75  DecayLength = dqm.book1dHisto("DecayLength","#tau Decay Length", 100 ,-20,20,"L_{#tau} (cm)","Number of #tau's from Gauge Bosons");
76  LifeTime = dqm.book1dHisto("LifeTime","#tau LifeTime ", 500 ,0,10000E-15,"#tau_{#tau} (s)","Number of #tau's from Gauge Bosons");
77 
78  TauSpinEffectsW_X = dqm.book1dHisto("TauSpinEffectsWX","X for pion", 50 ,0,1,"X","Number of #tau#rightarrow#pi#nu from W^{#pm} Bosons");
79  TauSpinEffectsHpm_X = dqm.book1dHisto("TauSpinEffectsHpmX","X for pion", 50 ,0,1,"X","Number of #tau#rightarrow#pi#nu from H^{#pm} Bosons");
80 
81  TauSpinEffectsW_eX = dqm.book1dHisto("TauSpinEffectsWeX","X for e", 50 ,0,1,"X","Number of #tau#rightarrowe#nu#nu from W^{#pm} Bosons");
82  TauSpinEffectsHpm_eX = dqm.book1dHisto("TauSpinEffectsHpmeX","X for e", 50 ,0,1,"X","Number of #tau#rightarrowe#nu#nu from H^{#pm} Bosons");
83 
84  TauSpinEffectsW_muX = dqm.book1dHisto("TauSpinEffectsWmuX","X for mu", 50 ,0,1,"X","Number of #tau#rightarrow#mu#nu#nu from W^{#pm} Bosons");
85  TauSpinEffectsHpm_muX = dqm.book1dHisto("TauSpinEffectsHpmmuX","X for mue", 50 ,0,1,"X","Number of #tau#rightarrow#mu#nu#nu from H^{#pm} Bosons");
86 
87  TauSpinEffectsW_UpsilonRho = dqm.book1dHisto("TauSpinEffectsWUpsilonRho","#Upsilon for #rho", 50 ,-1,1,"#Upsilon","Number of #tau#rightarrow#rho#nu from Gauge Bosons");
88  TauSpinEffectsHpm_UpsilonRho = dqm.book1dHisto("TauSpinEffectsHpmUpsilonRho","#Upsilon for #rho", 50 ,-1,1,"#Upsilon","Number of #tau#rightarrow#rho#nu from Gauge Bosons");
89 
90  TauSpinEffectsW_UpsilonA1 = dqm.book1dHisto("TauSpinEffectsWUpsilonA1","#Upsilon for a1", 50 ,-1,1,"#Upsilon","Number of #tau#rightarrow#pi#pi#pi#nu from Gauge Bosons");
91  TauSpinEffectsHpm_UpsilonA1 = dqm.book1dHisto("TauSpinEffectsHpmUpsilonA1","#Upsilon for a1", 50 ,-1,1,"#Upsilon","Number of #tau#rightarrow#pi#pi#pi#nu from Gauge Bosons");
92 
93  TauSpinEffectsH_pipiAcoplanarity = dqm.book1dHisto("TauSpinEffectsH_pipiAcoplanarity","H Acoplanarity for #pi^{-}#pi^{+}", 50 ,0,2*TMath::Pi(),"Acoplanarity","Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
94 
95  TauSpinEffectsH_pipiAcollinearity = dqm.book1dHisto("TauSpinEffectsH_pipiAcollinearity","H Acollinearity for #pi^{-}#pi^{+}", 50 ,0,TMath::Pi(),"Acollinearity","Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
96  TauSpinEffectsH_pipiAcollinearityzoom = dqm.book1dHisto("TauSpinEffectsH_pipiAcollinearityzoom","H Acollinearity for #pi^{-}#pi^{+}", 50 ,3,TMath::Pi(),"Acollinearity","Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
97 
98  TauSpinEffectsZ_MVis = dqm.book1dHisto("TauSpinEffectsZMVis","Mass of pi+ pi-", 25 ,0,1.1,"M_{#pi^{+}#pi^{-}} (GeV)","Number of Z#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
99  TauSpinEffectsH_MVis = dqm.book1dHisto("TauSpinEffectsHMVis","Mass of pi+ pi-", 25 ,0,1.1,"M_{#pi^{+}#pi^{-}} (GeV)","Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu) Events");
100 
101  TauSpinEffectsZ_Zs = dqm.book1dHisto("TauSpinEffectsZZs","Z_{s}", zsbins ,zsmin,zsmax,"Z_{s}","Number of Z#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu Events");
102  TauSpinEffectsH_Zs = dqm.book1dHisto("TauSpinEffectsHZs","Z_{s}", zsbins ,zsmin,zsmax,"Z_{s}","Number of H#rightarrow#tau^{-}(#rightarrow#pi^{-}#nu)#tau^{+}(#rightarrow#pi^{+}#nu Events");
103 
104  TauSpinEffectsZ_X= dqm.book1dHisto("TauSpinEffectsZX","X for pion of #tau^{-}", 25 ,0,1.0,"X","Number of #tau#rightarrow#pi#nu from Z Bosons");
105  TauSpinEffectsZ_X50to75= dqm.book1dHisto("TauSpinEffectsZX50to75","X for pion of #tau^{-} (50GeV-75GeV)", 10 ,0,1.0,"X","Number of #tau#rightarrow#pi#nu from Z(50GeV<M<75GeV) Bosons");
106  TauSpinEffectsZ_X75to88= dqm.book1dHisto("TauSpinEffectsZX75to88","X for pion of #tau^{-} (75GeV-88GeV)", 10 ,0,1.0,"X","Number of #tau#rightarrow#pi#nu from Z(75GeV<M<88GeV) Bosons");
107  TauSpinEffectsZ_X88to100= dqm.book1dHisto("TauSpinEffectsZX88to100","X for pion of #tau^{-} (88GeV-100GeV)", 10 ,0,1.0,"X","Number of #tau#rightarrow#pi#nu from Z(88GeV<M<100GeV) Bosons");
108  TauSpinEffectsZ_X100to120= dqm.book1dHisto("TauSpinEffectsZX100to120","X for pion of #tau^{-} (100GeV-120GeV)", 10 ,0,1.0,"X","Number of #tau#rightarrow#pi#nu from Z(100GeV<M<120GeV) Bosons");
109  TauSpinEffectsZ_X120UP= dqm.book1dHisto("TauSpinEffectsZX120UP","X for pion of #tau^{-} (>120GeV)", 10 ,0,1.0,"X","Number of #tau#rightarrow#pi#nu from Z(120GeV<MGeV) Bosons");
110 
111 
112  TauSpinEffectsH_X= dqm.book1dHisto("TauSpinEffectsH_X","X for pion of #tau^{-}", 25 ,0,1.0,"X","Number of #tau#rightarrow#pi#nu from H Bosons");
113 
114  TauSpinEffectsZ_Xf = dqm.book1dHisto("TauSpinEffectsZXf","X for pion of forward emitted #tau^{-}", 25 ,0,1.0,"X_{f}","Number of #tau#rightarrow#pi#nu from Z Bosons");
115  TauSpinEffectsH_Xf = dqm.book1dHisto("TauSpinEffectsHXf","X for pion of forward emitted #tau^{-}", 25 ,0,1.0,"X_{f}","Number of #tau#rightarrow#pi#nu from H Bosons");
116 
117  TauSpinEffectsZ_Xb = dqm.book1dHisto("TauSpinEffectsZXb","X for pion of backward emitted #tau^{-}", 25 ,0,1.0,"X_{b}","Number of #tau#rightarrow#pi#nu from Z Bosons");
118  TauSpinEffectsH_Xb = dqm.book1dHisto("TauSpinEffectsHXb","X for pion of backward emitted #tau^{-}", 25 ,0,1.0,"X_{b}","Number of #tau#rightarrow#pi#nu from H Bosons");
119 
120  TauSpinEffectsZ_eX = dqm.book1dHisto("TauSpinEffectsZeX","X for e", 50 ,0,1,"X","Number of #tau#rightarrowe#nu#nu from Gauge Bosons");
121  TauSpinEffectsH_eX = dqm.book1dHisto("TauSpinEffectsHeX","X for e", 50 ,0,1,"X","Number of #tau#rightarrowe#nu#nu from Gauge Bosons");
122 
123  TauSpinEffectsZ_muX = dqm.book1dHisto("TauSpinEffectsZmuX","X for mu", 50 ,0,1,"X","Number of #tau#rightarrow#mu#nu#nu from Gauge Bosons");
124  TauSpinEffectsH_muX = dqm.book1dHisto("TauSpinEffectsHmuX","X for mu", 50 ,0,1,"X","Number of #tau#rightarrow#mu#nu#nu from Gauge Bosons");
125 
126  TauSpinEffectsH_rhorhoAcoplanarityminus = dqm.book1dHisto("TauSpinEffectsH_rhorhoAcoplanarityminus","#phi^{*-} (acoplanarity) for Higgs #rightarrow #rho-#rho (y_{1}*y_{2}<0)", 32 ,0,2*TMath::Pi(),"#phi^{*-} (Acoplanarity)","Number of H#rightarrow#tau^{-}(#rightarrow#rho^{-}#nu)#tau^{+}(#rightarrow#rho^{+}#nu) Events");
127  TauSpinEffectsH_rhorhoAcoplanarityplus = dqm.book1dHisto("TauSpinEffectsH_rhorhoAcoplanarityplus","#phi^{*+} (acoplanarity) for Higgs #rightarrow #rho-#rho (y_{1}*y_{2}>0)", 32 ,0,2*TMath::Pi(),"#phi^{*+} (Acoplanarity)","Number of H#rightarrow#tau^{-}(#rightarrow#rho^{-}#nu)#tau^{+}(#rightarrow#rho^{+}#nu) Events");
128 
129  TauFSRPhotonsN=dqm.book1dHisto("TauFSRPhotonsN","FSR Photons radiating from/with tau (Gauge Boson)", 5 ,-0.5,4.5,"N^{FSR Photons radiating from/with #tau}","Number of #tau's from Gauge Bosons");
130  TauFSRPhotonsPt=dqm.book1dHisto("TauFSRPhotonsPt","Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100,"P_{t}^{FSR Photons radiating from/with #tau [per #tau]} (GeV)","Number of #tau's from Gauge Bosons");
131  TauFSRPhotonsPtSum=dqm.book1dHisto("TauFSRPhotonsPtSum","Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100,"P_{t}^{FSR Photons radiating from/with #tau [per #tau]} (GeV)","Number of #tau's from Gauge Bosons");
132 
133  TauBremPhotonsN=dqm.book1dHisto("TauBremPhotonsN","Brem. Photons radiating in tau decay", 5 ,-0.5,4.5,"N FSR Photons radiating from/with tau","Number of #tau's from Gauge Bosons");
134  TauBremPhotonsPt=dqm.book1dHisto("TauBremPhotonsPt","Sum Brem Pt ", 100 ,0,100,"P_{t}^{Brem. Photons radiating in tau decay} (GeV)","Number of #tau's from Gauge Bosons");
135  TauBremPhotonsPtSum =dqm.book1dHisto("TauBremPhotonsPtSum","Sum of Brem Pt ", 100 ,0,100,"Sum P_{t}^{Brem. Photons radiating in tau decay} (GeV)","Number of #tau's from Gauge Bosons");
136 
137  MODEID =dqm.book1dHisto("JAKID","JAK ID",NMODEID+1,-0.5,NMODEID+0.5);
138  for(unsigned int j=0; j<NMODEID+1;j++){
139  MODEInvMass.push_back(std::vector<MonitorElement *>());
140  std::string tmp="JAKID";
141  tmp+=std::to_string(j);
142  MODEInvMass.at(j).push_back(dqm.book1dHisto("M"+tmp,"M_{"+TauDecay::DecayMode(j)+"} (GeV)", 80 ,0,2.0,"M_{"+TauDecay::DecayMode(j)+"} (GeV)","Number of #tau's from Gauge Bosons"));
143  MODEID->setBinLabel(1+j,TauDecay::DecayMode(j));
144  if(j==TauDecay::MODE_3PI || j==TauDecay::MODE_PI2PI0 ||
145  j==TauDecay::MODE_KPIK ||
146  j==TauDecay::MODE_KPIPI ){
147  MODEInvMass.at(j).push_back(dqm.book1dHisto("M13"+tmp,"M_{13,"+TauDecay::DecayMode(j)+"} (GeV)", 80 ,0,2.0,"M_{13,"+TauDecay::DecayMode(j)+"} (GeV)","Number of #tau's from Gauge Bosons"));
148  MODEInvMass.at(j).push_back(dqm.book1dHisto("M23"+tmp,"M_{23,"+TauDecay::DecayMode(j)+"} (GeV)", 80 ,0,2.0,"M_{23,"+TauDecay::DecayMode(j)+"} (GeV)","Number of #tau's from Gauge Bosons"));
149  MODEInvMass.at(j).push_back(dqm.book1dHisto("M12"+tmp,"M_{12,"+TauDecay::DecayMode(j)+"} (GeV)", 80 ,0,2.0,"M_{12,"+TauDecay::DecayMode(j)+"} (GeV)","Number of #tau's from Gauge Bosons"));
150  }
151  }
152  return;
153 }
154 
155 
156 void TauValidation::analyze(const edm::Event& iEvent,const edm::EventSetup& iSetup){
158  edm::Handle<reco::GenParticleCollection> genParticles;
159  iEvent.getByToken(genparticleCollectionToken_, genParticles );
160 
161  double weight = 1.0;//= wmanager_.weight(iEvent);
163  // find taus
164  for (reco::GenParticleCollection::const_iterator iter = genParticles->begin(); iter != genParticles->end(); ++iter) {
165  if(abs(iter->pdgId())==PdtPdgMini::Z0 || abs(iter->pdgId())==PdtPdgMini::Higgs0){
166  spinEffectsZH(&(*iter),weight);
167  }
168  if(abs(iter->pdgId())==15){
169  if(isLastTauinChain(&(*iter))){
170  nTaus->Fill(0.5,weight);
171  int mother = tauMother(&(*iter),weight);
172  if(mother>-1){ // exclude B, D and other non-signal decay modes
173  nPrimeTaus->Fill(0.5,weight);
174  TauPt->Fill(iter->pt(),weight);
175  TauEta->Fill(iter->eta(),weight);
176  TauPhi->Fill(iter->phi(),weight);
177  photons(&(*iter),weight);
179  // Adding MODEID and Mass information
180  TauDecay_GenParticle TD;
181  unsigned int jak_id, TauBitMask;
182  if(TD.AnalyzeTau(&(*iter),jak_id,TauBitMask,false,false)){
183  MODEID->Fill(jak_id,weight);
184  TauProngs->Fill(TD.nProng(TauBitMask),weight);
185  tauDecayChannel(&(*iter),jak_id,TauBitMask,weight);
186  if(jak_id<=NMODEID){
187  int tcharge=iter->pdgId()/abs(iter->pdgId());
188  std::vector<const reco::GenParticle*> part=TD.Get_TauDecayProducts();
189  spinEffectsWHpm(&(*iter),mother,jak_id,part,weight);
190  TLorentzVector LVQ(0,0,0,0);
191  TLorentzVector LVS12(0,0,0,0);
192  TLorentzVector LVS13(0,0,0,0);
193  TLorentzVector LVS23(0,0,0,0);
194  bool haspart1=false;
195  TVector3 PV,SV;
196  bool hasDL(false);
197  for(unsigned int i=0;i<part.size();i++){
198  if(abs(part.at(i)->pdgId())!=PdtPdgMini::nu_tau && TD.isTauFinalStateParticle(part.at(i)->pdgId()) && !hasDL){
199  hasDL=true;
200  TLorentzVector tlv(iter->px(),iter->py(),iter->pz(),iter->energy());
201  PV=TVector3(iter->vx(),iter->vy(),iter->vz());
202  SV=TVector3(part.at(i)->vx(),part.at(i)->vy(),part.at(i)->vz());
203  TVector3 DL=SV-PV;
204  DecayLength->Fill(DL.Dot(tlv.Vect())/tlv.P(),weight);
205  double c(2.99792458E8),Ltau(DL.Mag()/100)/*cm->m*/,beta(iter->p()/iter->mass());
206  LifeTime->Fill( Ltau/(c*beta),weight);
207  }
208 
209  if(TD.isTauFinalStateParticle(part.at(i)->pdgId()) && abs(part.at(i)->pdgId())!=PdtPdgMini::nu_e &&
210  abs(part.at(i)->pdgId())!=PdtPdgMini::nu_mu && abs(part.at(i)->pdgId())!=PdtPdgMini::nu_tau ){
211  TLorentzVector LV(part.at(i)->px(),part.at(i)->py(),part.at(i)->pz(),part.at(i)->energy());
212  LVQ+=LV;
213  if(jak_id==TauDecay::MODE_3PI || jak_id==TauDecay::MODE_PI2PI0 ||
214  jak_id==TauDecay::MODE_KPIK ||
215  jak_id==TauDecay::MODE_KPIPI
216  ){
217  if((tcharge==part.at(i)->pdgId()/abs(part.at(i)->pdgId()) && TD.nProng(TauBitMask)==3) || ((jak_id==TauDecay::MODE_3PI || jak_id==TauDecay::MODE_PI2PI0) && TD.nProng(TauBitMask)==1 && abs(part.at(i)->pdgId())==PdtPdgMini::pi_plus) ){
218  LVS13+=LV;
219  LVS23+=LV;
220  }
221  else{
222  LVS12+=LV;
223  if(!haspart1 && ((jak_id==TauDecay::MODE_3PI || jak_id==TauDecay::MODE_PI2PI0) || ((jak_id!=TauDecay::MODE_3PI || jak_id==TauDecay::MODE_PI2PI0) && abs(part.at(i)->pdgId())==PdtPdgMini::K_plus) )){
224  LVS13+=LV;
225  haspart1=true;
226  }
227  else{
228  LVS23+=LV;
229  }
230  }
231  }
232  }
233  }
234  part.clear();
235  MODEInvMass.at(jak_id).at(0)->Fill(LVQ.M(),weight);
236  if(jak_id==TauDecay::MODE_3PI || jak_id==TauDecay::MODE_PI2PI0 ||
237  jak_id==TauDecay::MODE_KPIK ||
238  jak_id==TauDecay::MODE_KPIPI
239  ){
240  MODEInvMass.at(jak_id).at(1)->Fill(LVS13.M(),weight);
241  MODEInvMass.at(jak_id).at(2)->Fill(LVS23.M(),weight);
242  MODEInvMass.at(jak_id).at(3)->Fill(LVS12.M(),weight);
243  }
244  }
245  }
246  else{
247  MODEID->Fill(jak_id,weight);
248  }
249  }
250  }
251  }
252  }
253 }//analyze
254 
255 const reco::GenParticle* TauValidation::GetMother(const reco::GenParticle* tau){
256  for (unsigned int i=0;i<tau->numberOfMothers();i++) {
257  const reco::GenParticle *mother=static_cast<const reco::GenParticle*>(tau->mother(i));
258  if(mother->pdgId() == tau->pdgId()) return GetMother(mother);
259  return mother;
260  }
261  return tau;
262 }
263 
264 const std::vector<const reco::GenParticle*> TauValidation::GetMothers(const reco::GenParticle* boson){
265  std::vector<const reco::GenParticle*> mothers;
266  for (unsigned int i=0;i<boson->numberOfMothers();i++) {
267  const reco::GenParticle *mother=static_cast<const reco::GenParticle*>(boson->mother(i));
268  if(mother->pdgId() == boson->pdgId()) return GetMothers(mother);
269  mothers.push_back(mother);
270  }
271  return mothers;
272 }
273 
274 int TauValidation::findMother(const reco::GenParticle* tau){
275  return TauValidation::GetMother(tau)->pdgId();
276 }
277 
278 bool TauValidation::isLastTauinChain(const reco::GenParticle* tau){
279  for(unsigned int i = 0; i <tau->numberOfDaughters(); i++){
280  if(tau->daughter(i)->pdgId() == tau->pdgId()) return false;
281  }
282  return true;
283 }
284 
285 void TauValidation::findTauList(const reco::GenParticle* tau,std::vector<const reco::GenParticle*> &TauList){
286  TauList.insert(TauList.begin(),tau);
287  for(unsigned int i=0;i<tau->numberOfMothers();i++) {
288  const reco::GenParticle *mother=static_cast<const reco::GenParticle*>(tau->mother(i));
289  if(mother->pdgId() == tau->pdgId()){
290  findTauList(mother,TauList);
291  }
292  }
293 }
294 
295 void TauValidation::findFSRandBrem(const reco::GenParticle* p, bool doBrem, std::vector<const reco::GenParticle*> &ListofFSR,
296  std::vector<const reco::GenParticle*> &ListofBrem){
297  // note this code split the FSR and Brem based one if the tau decays into a tau+photon or not with the Fortran Tauola Interface, this is not 100% correct because photos puts the tau with the regular tau decay products.
298  if(abs(p->pdgId())==15){
299  if(isLastTauinChain(p)){ doBrem=true;}
300  else{ doBrem=false;}
301  }
302  int photo_ID=22;
303  for(unsigned int i = 0; i <p->numberOfDaughters(); i++){
304  const reco::GenParticle *dau=static_cast<const reco::GenParticle*>(p->daughter(i));
305  if(abs((dau)->pdgId()) == abs(photo_ID) && !doBrem){ListofFSR.push_back(dau);}
306  if(abs((dau)->pdgId()) == abs(photo_ID) && doBrem){ListofBrem.push_back(dau);}
307  if(abs((dau)->pdgId()) != 111 && abs((dau)->pdgId()) != 221){ // remove pi0 and eta decays
308  findFSRandBrem(dau,doBrem,ListofFSR,ListofBrem);
309  }
310  }
311 }
312 
313 void TauValidation::FindPhotosFSR(const reco::GenParticle* p,std::vector<const reco::GenParticle*> &ListofFSR,double &BosonScale){
314  BosonScale=0.0;
315  const reco::GenParticle* m=GetMother(p);
316  int mother_pid=m->pdgId();
317  if(m->pdgId()!=p->pdgId()){
318  for(unsigned int i=0; i <m->numberOfDaughters(); i++){
319  const reco::GenParticle *dau=static_cast<const reco::GenParticle*>(m->daughter(i));
320  if(abs(dau->pdgId()) == 22) {
321  ListofFSR.push_back(dau);
322  }
323  }
324  }
325  if(abs(mother_pid) == 24) BosonScale=1.0; // W
326  if(abs(mother_pid) == 23) BosonScale=2.0; // Z;
327  if(abs(mother_pid) == 22) BosonScale=2.0; // gamma;
328  if(abs(mother_pid) == 25) BosonScale=2.0; // HSM;
329  if(abs(mother_pid) == 35) BosonScale=2.0; // H0;
330  if(abs(mother_pid) == 36) BosonScale=2.0; // A0;
331  if(abs(mother_pid) == 37) BosonScale=1.0; //Hpm;
332 }
333 
334 int TauValidation::tauMother(const reco::GenParticle* tau, double weight){
335  if(abs(tau->pdgId()) != 15 ) return -3;
336  int mother_pid = findMother(tau);
337  if(mother_pid == -2) return -2;
338  int label = other;
339  if(abs(mother_pid) == 24) label = W;
340  if(abs(mother_pid) == 23) label = Z;
341  if(abs(mother_pid) == 22) label = gamma;
342  if(abs(mother_pid) == 25) label = HSM;
343  if(abs(mother_pid) == 35) label = H0;
344  if(abs(mother_pid) == 36) label = A0;
345  if(abs(mother_pid) == 37) label = Hpm;
346  int mother_shortpid=(abs(mother_pid)%10000);
347  if(mother_shortpid>500 && mother_shortpid<600 )label = B;
348  if(mother_shortpid>400 && mother_shortpid<500)label = D;
349  TauMothers->Fill(label,weight);
350  if(label==B || label == D || label == other) return -1;
351  return mother_pid;
352 }
353 
354 int TauValidation::tauDecayChannel(const reco::GenParticle* tau,int jak_id, unsigned int TauBitMask, double weight){
355  int channel = undetermined;
356  if(tau->status() == 1) channel = stable;
357  int allCount = 0,
358  eCount = 0,
359  muCount = 0,
360  pi0Count = 0,
361  piCount = 0,
362  rhoCount = 0,
363  a1Count = 0,
364  KCount = 0,
365  KstarCount = 0;
366 
367  countParticles(tau,allCount,eCount,muCount,pi0Count,piCount,rhoCount,a1Count,KCount,KstarCount);
368 
369  // resonances
370  if(KCount >= 1) channel = K;
371  if(KstarCount >= 1) channel = Kstar;
372  if(a1Count >= 1) channel = a1;
373  if(rhoCount >= 1) channel = rho;
374  if(channel!=undetermined && weight!=0.0) TauDecayChannels->Fill(channel,weight);
375 
376  // final state products
377  if(piCount == 1 && pi0Count == 0) channel = pi;
378  if(piCount == 1 && pi0Count == 1) channel = pi1pi0;
379  if(piCount == 1 && pi0Count > 1) channel = pinpi0;
380  if(piCount == 3 && pi0Count == 0) channel = tripi;
381  if(piCount == 3 && pi0Count > 0) channel = tripinpi0;
382  if(eCount == 1) channel = electron;
383  if(muCount == 1) channel = muon;
384  if(weight!=0.0) TauDecayChannels->Fill(channel,weight);
385  return channel;
386 }
387 
388 void TauValidation::countParticles(const reco::GenParticle* p,int &allCount, int &eCount, int &muCount,
389  int &pi0Count,int &piCount,int &rhoCount,int &a1Count,int &KCount,int &KstarCount){
390  for(unsigned int i=0; i<p->numberOfDaughters(); i++){
391  const reco::GenParticle *dau=static_cast<const reco::GenParticle*>(p->daughter(i));
392  int pid = dau->pdgId();
393  allCount++;
394  if(abs(pid) == 11) eCount++;
395  if(abs(pid) == 13) muCount++;
396  if(abs(pid) == 111) pi0Count++;
397  if(abs(pid) == 211) piCount++;
398  if(abs(pid) == 213) rhoCount++;
399  if(abs(pid) == 20213) a1Count++;
400  if(abs(pid) == 321) KCount++;
401  if(abs(pid) == 323) KstarCount++;
402  countParticles(dau,allCount,eCount,muCount,pi0Count,piCount,rhoCount,a1Count,KCount,KstarCount);
403  }
404 }
405 
406 
407 
408 void TauValidation::spinEffectsWHpm(const reco::GenParticle* tau,int mother, int decay, std::vector<const reco::GenParticle*> &part,double weight){
409  if(decay == TauDecay::MODE_PION || decay == TauDecay::MODE_MUON || decay == TauDecay::MODE_ELECTRON){ // polarization only for 1-prong hadronic taus with no neutral pions
410  TLorentzVector momP4 = motherP4(tau);
411  TLorentzVector pionP4 = leadingPionP4(tau);
412  pionP4.Boost(-1*momP4.BoostVector());
413  double energy = pionP4.E()/(momP4.M()/2);
414  if(decay == TauDecay::MODE_PION){
415  if(abs(mother) == 24) TauSpinEffectsW_X->Fill(energy,weight);
416  if(abs(mother) == 37) TauSpinEffectsHpm_X->Fill(energy,weight);
417  }
418  if(decay == TauDecay::MODE_MUON){
419  if(abs(mother) == 24) TauSpinEffectsW_muX->Fill(energy,weight);
420  if(abs(mother) == 37) TauSpinEffectsHpm_muX->Fill(energy,weight);
421  }
422  if(decay == TauDecay::MODE_ELECTRON){
423  if(abs(mother) == 24) TauSpinEffectsW_eX->Fill(energy,weight);
424  if(abs(mother) == 37) TauSpinEffectsHpm_eX->Fill(energy,weight);
425  }
426  }
427  else if(decay==TauDecay::MODE_PIPI0){
428  TLorentzVector rho(0,0,0,0),pi(0,0,0,0);
429  for(unsigned int i=0;i<part.size();i++){
430  TLorentzVector LV(part.at(i)->px(),part.at(i)->py(),part.at(i)->pz(),part.at(i)->energy());
431  if(abs(part.at(i)->pdgId())==PdtPdgMini::pi_plus){pi+=LV; rho+=LV;}
432  if(abs(part.at(i)->pdgId())==PdtPdgMini::pi0){rho+=LV;}
433  }
434  if(abs(mother) == 24) TauSpinEffectsW_UpsilonRho->Fill(2*pi.P()/rho.P()-1,weight);
435  if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonRho->Fill(2*pi.P()/rho.P()-1,weight);
436  }
437  else if(decay==TauDecay::MODE_3PI || decay==TauDecay::MODE_PI2PI0){ // only for pi2pi0 for now
438  TLorentzVector a1(0,0,0,0),pi_p(0,0,0,0),pi_m(0,0,0,0);
439  int nplus(0),nminus(0);
440  for(unsigned int i=0;i<part.size();i++){
441  TLorentzVector LV(part.at(i)->px(),part.at(i)->py(),part.at(i)->pz(),part.at(i)->energy());
442  if(part.at(i)->pdgId()==PdtPdgMini::pi_plus){ pi_p+=LV; a1+=LV; nplus++;}
443  if(part.at(i)->pdgId()==PdtPdgMini::pi_minus){pi_m+=LV; a1+=LV; nminus++;}
444  }
445  double gamma=0;
446  if(nplus+nminus==3 && nplus==1) gamma=2*pi_p.P()/a1.P()-1;
447  if(nplus+nminus==3 && nminus==1) gamma=2*pi_m.P()/a1.P()-1;
448  else{
449  pi_p+=pi_m; gamma=2*pi_p.P()/a1.P()-1;
450  }
451  if(abs(mother) == 24) TauSpinEffectsW_UpsilonA1->Fill(gamma,weight);
452  if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonA1->Fill(gamma,weight);
453  }
454 }
455 
456 void TauValidation::spinEffectsZH(const reco::GenParticle* boson, double weight){
457  int ntau(0);
458  for(unsigned int i = 0; i <boson->numberOfDaughters(); i++){
459  const reco::GenParticle *dau=static_cast<const reco::GenParticle*>(boson->daughter(i));
460  if(ntau==1 && dau->pdgId() == 15)return;
461  if(boson->pdgId()!= 15 && abs(dau->pdgId()) == 15)ntau++;
462  }
463  if(ntau!=2) return;
464  if(abs(boson->pdgId())==PdtPdgMini::Z0 || abs(boson->pdgId())==PdtPdgMini::Higgs0){
465  TLorentzVector tautau(0,0,0,0);
466  TLorentzVector pipi(0,0,0,0);
467  TLorentzVector taum(0,0,0,0);
468  TLorentzVector taup(0,0,0,0);
469  TLorentzVector rho_plus,rho_minus,pi_rhominus,pi0_rhominus,pi_rhoplus,pi0_rhoplus,pi_plus,pi_minus;
470  bool hasrho_minus(false),hasrho_plus(false),haspi_minus(false),haspi_plus(false);
471  int nSinglePionDecays(0),nSingleMuonDecays(0),nSingleElectronDecays(0);
472  double x1(0),x2(0);
473  TLorentzVector Zboson(boson->px(),boson->py(),boson->pz(),boson->energy());
474  for(unsigned int i = 0; i <boson->numberOfDaughters(); i++){
475  const reco::GenParticle *dau=static_cast<const reco::GenParticle*>(boson->daughter(i));
476  int pid = dau->pdgId();
477  if(abs(findMother(dau)) != 15 && abs(pid) == 15){
478  TauDecay_GenParticle TD;
479  unsigned int jak_id, TauBitMask;
480  if(TD.AnalyzeTau(dau,jak_id,TauBitMask,false,false)){
481  std::vector<const reco::GenParticle*> part=TD.Get_TauDecayProducts();
482  if(jak_id==TauDecay::MODE_PION || jak_id==TauDecay::MODE_MUON || jak_id==TauDecay::MODE_ELECTRON){
483  if(jak_id==TauDecay::MODE_PION) nSinglePionDecays++;
484  if(jak_id==TauDecay::MODE_MUON) nSingleMuonDecays++;
485  if(jak_id==TauDecay::MODE_ELECTRON) nSingleElectronDecays++;
486  TLorentzVector LVtau(dau->px(),dau->py(),dau->pz(),dau->energy());
487  tautau += LVtau;
488  TLorentzVector LVpi=leadingPionP4(dau);
489  pipi+=LVpi;
490  const HepPDT::ParticleData* pd = fPDGTable->particle(dau->pdgId ());
491  int charge = (int) pd->charge();
492  LVtau.Boost(-1*Zboson.BoostVector());
493  LVpi.Boost(-1*Zboson.BoostVector());
494 
495 
496 
497  if(jak_id==TauDecay::MODE_MUON){
498  if(abs(boson->pdgId())==PdtPdgMini::Z0) TauSpinEffectsZ_muX->Fill(LVpi.P()/LVtau.E(),weight);
499  if(abs(boson->pdgId())==PdtPdgMini::Higgs0) TauSpinEffectsH_muX->Fill(LVpi.P()/LVtau.E(),weight);
500 
501  }
502  if(jak_id==TauDecay::MODE_ELECTRON){
503  if(abs(boson->pdgId())==PdtPdgMini::Z0) TauSpinEffectsZ_eX->Fill(LVpi.P()/LVtau.E(),weight);
504  if(abs(boson->pdgId())==PdtPdgMini::Higgs0) TauSpinEffectsH_eX->Fill(LVpi.P()/LVtau.E(),weight);
505  }
506 
507 
508  if(jak_id==TauDecay::MODE_PION){
509  if(abs(boson->pdgId())==PdtPdgMini::Z0){
510  TauSpinEffectsZ_X->Fill(LVpi.P()/LVtau.E(),weight);
511  if(50.0<Zboson.M() && Zboson.M()<75.0) TauSpinEffectsZ_X50to75->Fill(LVpi.P()/LVtau.E(),weight);
512  if(75.0<Zboson.M() && Zboson.M()<88.0) TauSpinEffectsZ_X75to88->Fill(LVpi.P()/LVtau.E(),weight);
513  if(88.0<Zboson.M() && Zboson.M()<100.0) TauSpinEffectsZ_X88to100->Fill(LVpi.P()/LVtau.E(),weight);
514  if(100.0<Zboson.M() && Zboson.M()<120.0) TauSpinEffectsZ_X100to120->Fill(LVpi.P()/LVtau.E(),weight);
515  if(120.0<Zboson.M()) TauSpinEffectsZ_X120UP->Fill(LVpi.P()/LVtau.E(),weight);
516  }
517  if(abs(boson->pdgId())==PdtPdgMini::Higgs0) TauSpinEffectsH_X->Fill(LVpi.P()/LVtau.E(),weight);
518  }
519  if(charge<0){x1=LVpi.P()/LVtau.E(); taum=LVtau;}
520  else{ x2=LVpi.P()/LVtau.E();}
521  }
522  TLorentzVector LVtau(dau->px(),dau->py(),dau->pz(),dau->energy());
523  if(pid == 15)taum=LVtau;
524  if(pid ==-15)taup=LVtau;
525  if(jak_id==TauDecay::MODE_PIPI0){
526  for(unsigned int i=0; i<part.size();i++){
527  int pid_d = part.at(i)->pdgId();
528  if(abs(pid_d)==211 || abs(pid_d)==111){
529  TLorentzVector LV(part.at(i)->px(),part.at(i)->py(),part.at(i)->pz(),part.at(i)->energy());
530  if(pid==15){
531  hasrho_minus=true;
532  if(pid_d==-211 ){ pi_rhominus=LV;}
533  if(abs(pid_d)==111 ){ pi0_rhominus=LV;}
534  }
535  if(pid==-15){
536  hasrho_plus=true;
537  if(pid_d==211 ){pi_rhoplus=LV;}
538  if(abs(pid_d)==111 ){pi0_rhoplus=LV;}
539  }
540  }
541  }
542  }
543  if(jak_id==TauDecay::MODE_PION){
544  for(unsigned int i=0; i<part.size();i++){
545  int pid_d = part.at(i)->pdgId();
546  if(abs(pid_d)==211 ){
547  TLorentzVector LV(part.at(i)->px(),part.at(i)->py(),part.at(i)->pz(),part.at(i)->energy());
548  if(pid==15){
549  haspi_minus=true;
550  if(pid_d==-211 ){ pi_minus=LV;}
551  }
552  if(pid==-15){
553  haspi_plus=true;
554  if(pid_d==211 ){pi_plus=LV;}
555  }
556  }
557  }
558  }
559  }
560  }
561  }
562  if(hasrho_minus && hasrho_plus){
563  //compute rhorho
564  rho_minus=pi_rhominus;
565  rho_minus+=pi0_rhominus;
566  rho_plus=pi_rhoplus;
567  rho_plus+=pi0_rhoplus;
568  TLorentzVector rhorho=rho_minus;rhorho+=rho_plus;
569 
570  // boost to rhorho cm
571  TLorentzVector pi_rhoplusb=pi_rhoplus; pi_rhoplusb.Boost(-1*rhorho.BoostVector());
572  TLorentzVector pi0_rhoplusb=pi0_rhoplus; pi0_rhoplusb.Boost(-1*rhorho.BoostVector());
573  TLorentzVector pi_rhominusb=pi_rhominus; pi_rhominusb.Boost(-1*rhorho.BoostVector());
574  TLorentzVector pi0_rhominusb=pi0_rhominus; pi0_rhominusb.Boost(-1*rhorho.BoostVector());
575 
576  // compute n+/-
577  TVector3 n_plus=pi_rhoplusb.Vect().Cross(pi0_rhoplusb.Vect());
578  TVector3 n_minus=pi_rhominusb.Vect().Cross(pi0_rhominusb.Vect());
579 
580  // compute the acoplanarity
581  double Acoplanarity=acos(n_plus.Dot(n_minus)/(n_plus.Mag()*n_minus.Mag()));
582  if(pi_rhominusb.Vect().Dot(n_plus)>0){Acoplanarity*=-1;Acoplanarity+=2*TMath::Pi();}
583 
584  // now boost to tau frame
585  pi_rhoplus.Boost(-1*taup.BoostVector());
586  pi0_rhoplus.Boost(-1*taup.BoostVector());
587  pi_rhominus.Boost(-1*taum.BoostVector());
588  pi0_rhominus.Boost(-1*taum.BoostVector());
589 
590  // compute y1 and y2
591  double y1=(pi_rhoplus.E()-pi0_rhoplus.E())/(pi_rhoplus.E()+pi0_rhoplus.E());
592  double y2=(pi_rhominus.E()-pi0_rhominus.E())/(pi_rhominus.E()+pi0_rhominus.E());
593 
594  // fill histograms
595  if(abs(boson->pdgId())==PdtPdgMini::Higgs0 && y1*y2<0) TauSpinEffectsH_rhorhoAcoplanarityminus->Fill(Acoplanarity,weight);
596  if(abs(boson->pdgId())==PdtPdgMini::Higgs0 && y1*y2>0) TauSpinEffectsH_rhorhoAcoplanarityplus->Fill(Acoplanarity,weight);
597  }
598  if(haspi_minus && haspi_plus){
599  TLorentzVector tauporig=taup;
600  TLorentzVector taumorig=taum;
601 
602  // now boost to Higgs frame
603  pi_plus.Boost(-1*Zboson.BoostVector());
604  pi_minus.Boost(-1*Zboson.BoostVector());
605 
606  taup.Boost(-1*Zboson.BoostVector());
607  taum.Boost(-1*Zboson.BoostVector());
608 
609  if(abs(boson->pdgId())==PdtPdgMini::Higgs0){
610  TauSpinEffectsH_pipiAcollinearity->Fill(acos(pi_plus.Vect().Dot(pi_minus.Vect())/(pi_plus.P()*pi_minus.P())));
611  TauSpinEffectsH_pipiAcollinearityzoom->Fill(acos(pi_plus.Vect().Dot(pi_minus.Vect())/(pi_plus.P()*pi_minus.P())));
612  }
613 
614  double proj_m=taum.Vect().Dot(pi_minus.Vect())/(taum.P()*taum.P());
615  double proj_p=taup.Vect().Dot(pi_plus.Vect())/(taup.P()*taup.P());
616  TVector3 Tau_m=taum.Vect();
617  TVector3 Tau_p=taup.Vect();
618  Tau_m*=proj_m;
619  Tau_p*=proj_p;
620  TVector3 Pit_m=pi_minus.Vect()-Tau_m;
621  TVector3 Pit_p=pi_plus.Vect()-Tau_p;
622 
623  double Acoplanarity=acos(Pit_m.Dot(Pit_p)/(Pit_p.Mag()*Pit_m.Mag()));
624  TVector3 n=Pit_p.Cross(Pit_m);
625  if(n.Dot(Tau_m)/Tau_m.Mag()>0){Acoplanarity*=-1; Acoplanarity+=2*TMath::Pi();}
626  // fill histograms
627  if(abs(boson->pdgId())==PdtPdgMini::Higgs0) TauSpinEffectsH_pipiAcoplanarity->Fill(Acoplanarity,weight);
628  taup=tauporig;
629  taum=taumorig;
630  }
631  if(nSinglePionDecays == 2 && tautau.M()!= 0) {
632  for(int i=0;i<zsbins;i++){
633  double zslow=((double)i)*(zsmax-zsmin)/((double)zsbins)+zsmin;
634  double zsup=((double)i+1)*(zsmax-zsmin)/((double)zsbins)+zsmin;
635  double aup=Zstoa(zsup), alow=Zstoa(zslow);
636  if(x2-x1>alow && x2-x1<aup){
637  double zs=(zsup+zslow)/2;
638  if(abs(boson->pdgId())==PdtPdgMini::Z0) TauSpinEffectsZ_Zs->Fill(zs,weight);
639  if(abs(boson->pdgId())==PdtPdgMini::Higgs0) TauSpinEffectsH_Zs->Fill(zs,weight);
640  break;
641  }
642  }
643  if(abs(boson->pdgId())==PdtPdgMini::Z0) TauSpinEffectsZ_MVis->Fill(pipi.M()/tautau.M(),weight);
644  if(abs(boson->pdgId())==PdtPdgMini::Higgs0) TauSpinEffectsH_MVis->Fill(pipi.M()/tautau.M(),weight);
645 
646  if(x1!=0){
647  const std::vector<const reco::GenParticle*> m=GetMothers(boson);
648  int q(0),qbar(0);
649  TLorentzVector Z(0,0,0,0);
650  for(unsigned int i=0;i<m.size();i++){
651  if(m.at(i)->pdgId()==PdtPdgMini::d || m.at(i)->pdgId()==PdtPdgMini::u ){q++;}
652  if(m.at(i)->pdgId()==PdtPdgMini::anti_d || m.at(i)->pdgId()==PdtPdgMini::anti_u ){qbar++;}
653  }
654  if(q==1 && qbar==1){// assume q has largest E (valence vs see quarks)
655  if(taum.Vect().Dot(Zboson.Vect())/(Zboson.P()*taum.P())>0){
656  if(abs(boson->pdgId())==PdtPdgMini::Z0) TauSpinEffectsZ_Xf->Fill(x1,weight);
657  if(abs(boson->pdgId())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xf->Fill(x1,weight);
658  }
659  else{
660  if(abs(boson->pdgId())==PdtPdgMini::Z0) TauSpinEffectsZ_Xb->Fill(x1,weight);
661  if(abs(boson->pdgId())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xb->Fill(x1,weight);
662  }
663  }
664  }
665  }
666  }
667 }
668 
669 double TauValidation::Zstoa(double zs){
670  double a=1-sqrt(fabs(1.0-2*fabs(zs)));
671  if(zs<0){
672  a*=-1.0;
673  }
674  return a;
675 }
676 
677 
678 double TauValidation::leadingPionMomentum(const reco::GenParticle* tau, double weight){
679  return leadingPionP4(tau).P();
680 }
681 
682 TLorentzVector TauValidation::leadingPionP4(const reco::GenParticle* tau){
683  TLorentzVector p4(0,0,0,0);
684  for(unsigned int i = 0; i <tau->numberOfDaughters(); i++){
685  const reco::GenParticle *dau=static_cast<const reco::GenParticle*>(tau->daughter(i));
686  int pid = dau->pdgId();
687  if(abs(pid) == 15) return leadingPionP4(dau);
688  if(!(abs(pid)==211 || abs(pid)==13 || abs(pid)==11)) continue;
689  if(dau->p() > p4.P()) p4 = TLorentzVector(dau->px(),dau->py(),dau->pz(),dau->energy());
690  }
691  return p4;
692 }
693 
694 TLorentzVector TauValidation::motherP4(const reco::GenParticle* tau){
695  const reco::GenParticle* m=GetMother(tau);
696  return TLorentzVector(m->px(),m->py(),m->pz(),m->energy());
697 }
698 
699 double TauValidation::visibleTauEnergy(const reco::GenParticle* tau){
700  TLorentzVector p4(tau->px(),tau->py(),tau->pz(),tau->energy());
701  for(unsigned int i = 0; i <tau->numberOfDaughters(); i++){
702  const reco::GenParticle *dau=static_cast<const reco::GenParticle*>(tau->daughter(i));
703  int pid = dau->pdgId();
704  if(abs(pid) == 15) return visibleTauEnergy(dau);
705  if(abs(pid) == 12 || abs(pid) == 14 || abs(pid) == 16) {
706  p4-=TLorentzVector(dau->px(),dau->py(),dau->pz(),dau->energy());
707  }
708  }
709  return p4.E();
710 }
711 
712 void TauValidation::photons(const reco::GenParticle* tau, double weight){
713  // Find First tau in chain
714  std::vector<const reco::GenParticle*> TauList;
715  findTauList(tau,TauList);
716 
717  // Get List of Gauge Boson to tau(s) FSR and Brem
718  bool passedW=false;
719  std::vector<const reco::GenParticle*> ListofFSR; ListofFSR.clear();
720  std::vector<const reco::GenParticle*> ListofBrem; ListofBrem.clear();
721  std::vector<const reco::GenParticle*> FSR_photos; FSR_photos.clear();
722  double BosonScale(1);
723  if(TauList.size()>0){
724  TauValidation::findFSRandBrem(TauList.at(0),passedW,ListofFSR,ListofBrem);
725  TauValidation::FindPhotosFSR(TauList.at(0),FSR_photos,BosonScale);
726 
727  // Add the Tau Brem. information
728  TauBremPhotonsN->Fill(ListofBrem.size(),weight);
729  double photonPtSum=0;
730  for(unsigned int i=0;i<ListofBrem.size();i++){
731  photonPtSum+=ListofBrem.at(i)->pt();
732  TauBremPhotonsPt->Fill(ListofBrem.at(i)->pt(),weight);
733  }
734  TauBremPhotonsPtSum->Fill(photonPtSum,weight);
735 
736  // Now add the Gauge Boson FSR information
737  if(BosonScale!=0){
738  TauFSRPhotonsN->Fill(ListofFSR.size(),weight);
739  photonPtSum=0;
740  for(unsigned int i=0;i<ListofFSR.size();i++){
741  photonPtSum+=ListofFSR.at(i)->pt();
742  TauFSRPhotonsPt->Fill(ListofFSR.at(i)->pt(),weight);
743  }
744  double FSR_photosSum(0);
745  for(unsigned int i=0;i<FSR_photos.size();i++){
746  FSR_photosSum+=FSR_photos.at(i)->pt();
747  TauFSRPhotonsPt->Fill(FSR_photos.at(i)->pt()/BosonScale,weight*BosonScale);
748  }
749  TauFSRPhotonsPtSum->Fill(photonPtSum+FSR_photosSum/BosonScale,weight);
750  }
751  }
752 }
753 
TauValidation::TauFSRPhotonsN
MonitorElement * TauFSRPhotonsN
Definition: TauValidation.h:92
beta
const double beta
Definition: CosmicMuonParameters.h:69
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
i
int i
Definition: DBlmapReader.cc:9
reco::LeafCandidate::pdgId
virtual int pdgId() const
PDG identifier.
Definition: LeafCandidate.h:236
reco::LeafCandidate::p
virtual double p() const
magnitude of momentum vector
Definition: LeafCandidate.h:157
TauValidation::findMother
int findMother(const reco::GenParticle *)
Definition: TauValidation.cc:274
TauValidation::~TauValidation
virtual ~TauValidation()
Definition: TauValidation.cc:29
TauDecay_GenParticle
Definition: TauDecay_GenParticle.h:23
TauValidation::TauSpinEffectsH_rhorhoAcoplanarityplus
MonitorElement * TauSpinEffectsH_rhorhoAcoplanarityplus
Definition: TauValidation.h:92
TauValidation::TauPhi
MonitorElement * TauPhi
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsW_X
MonitorElement * TauSpinEffectsW_X
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsHpm_UpsilonA1
MonitorElement * TauSpinEffectsHpm_UpsilonA1
Definition: TauValidation.h:92
TauValidation::MODEID
MonitorElement * MODEID
Definition: TauValidation.h:105
TauValidation::TauSpinEffectsZ_Xb
MonitorElement * TauSpinEffectsZ_Xb
Definition: TauValidation.h:92
TauValidation::nTaus
MonitorElement * nTaus
Definition: TauValidation.h:91
TauValidation::DecayLength
MonitorElement * DecayLength
Definition: TauValidation.h:92
TauValidation::MODEInvMass
std::vector< std::vector< MonitorElement * > > MODEInvMass
Definition: TauValidation.h:106
edm::Event::getByToken
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:434
reco::LeafCandidate::status
virtual int status() const
status word
Definition: LeafCandidate.h:240
TauValidation::TauValidation
TauValidation(const edm::ParameterSet &)
Definition: TauValidation.cc:18
TauValidation::TauPt
MonitorElement * TauPt
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsH_MVis
MonitorElement * TauSpinEffectsH_MVis
Definition: TauValidation.h:92
TauValidation::TauMothers
MonitorElement * TauMothers
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsHpm_eX
MonitorElement * TauSpinEffectsHpm_eX
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsZ_X100to120
MonitorElement * TauSpinEffectsZ_X100to120
Definition: TauValidation.h:92
TauValidation::FindPhotosFSR
void FindPhotosFSR(const reco::GenParticle *p, std::vector< const reco::GenParticle * > &ListofFSR, double &BosonScale)
Definition: TauValidation.cc:313
TauValidation::TauSpinEffectsW_muX
MonitorElement * TauSpinEffectsW_muX
Definition: TauValidation.h:92
TauValidation::TauBremPhotonsN
MonitorElement * TauBremPhotonsN
Definition: TauValidation.h:92
TauValidation::spinEffectsWHpm
void spinEffectsWHpm(const reco::GenParticle *, int, int, std::vector< const reco::GenParticle * > &part, double weight)
Definition: TauValidation.cc:408
MonitorElement::setBinLabel
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)
Definition: MonitorElement.cc:942
TauValidation::TauSpinEffectsZ_Xf
MonitorElement * TauSpinEffectsZ_Xf
Definition: TauValidation.h:92
TauValidation::isLastTauinChain
bool isLastTauinChain(const reco::GenParticle *tau)
Definition: TauValidation.cc:278
TauValidation::TauBremPhotonsPtSum
MonitorElement * TauBremPhotonsPtSum
Definition: TauValidation.h:92
genParticleCandidates2GenParticles_cfi.genParticles
tuple genParticles
Definition: genParticleCandidates2GenParticles_cfi.py:3
TauValidation::TauSpinEffectsH_pipiAcoplanarity
MonitorElement * TauSpinEffectsH_pipiAcoplanarity
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsHpm_UpsilonRho
MonitorElement * TauSpinEffectsHpm_UpsilonRho
Definition: TauValidation.h:92
DQMHelper::book1dHisto
MonitorElement * book1dHisto(std::string name, std::string title, int n, double xmin, double xmax, std::string xaxis, std::string yaxis)
Definition: DQMHelper.cc:12
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
TauValidation::leadingPionMomentum
double leadingPionMomentum(const reco::GenParticle *, double weight)
Definition: TauValidation.cc:678
AlCaHLTBitMon_QueryRunRegistry.string
string string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:255
edm::EventSetup::getData
void getData(T &iHolder) const
Definition: EventSetup.h:67
TauValidation::TauSpinEffectsZ_eX
MonitorElement * TauSpinEffectsZ_eX
Definition: TauValidation.h:92
MonitorElement::Fill
void Fill(long long x)
Definition: MonitorElement.h:166
TauValidation::TauProngs
MonitorElement * TauProngs
Definition: TauValidation.h:92
TauValidation::motherP4
TLorentzVector motherP4(const reco::GenParticle *)
Definition: TauValidation.cc:694
double_binary
Definition: DDStreamer.cc:227
TauDecay::isTauFinalStateParticle
bool isTauFinalStateParticle(int pdgid)
Definition: TauDecay.cc:36
TauValidation::dqmBeginRun
virtual void dqmBeginRun(const edm::Run &r, const edm::EventSetup &c) override
Definition: TauValidation.cc:31
TauValidation::TauSpinEffectsH_pipiAcollinearityzoom
MonitorElement * TauSpinEffectsH_pipiAcollinearityzoom
Definition: TauValidation.h:92
cond::to_string
std::string to_string(const T &t)
Definition: Logger.cc:26
LV
math::XYZTLorentzVectorD LV
Definition: HLTTauDQMPlotter.h:16
TauValidation::photons
void photons(const reco::GenParticle *, double weight)
Definition: TauValidation.cc:712
reco::LeafCandidate::energy
virtual double energy() const
energy
Definition: LeafCandidate.h:159
iEvent
int iEvent
Definition: GenABIO.cc:230
TauValidation::TauSpinEffectsW_UpsilonRho
MonitorElement * TauSpinEffectsW_UpsilonRho
Definition: TauValidation.h:92
TauValidation::GetMother
const reco::GenParticle * GetMother(const reco::GenParticle *tau)
Definition: TauValidation.cc:255
TauValidation::TauSpinEffectsZ_X88to100
MonitorElement * TauSpinEffectsZ_X88to100
Definition: TauValidation.h:92
TauDecay_GenParticle.h
TauValidation::TauSpinEffectsH_rhorhoAcoplanarityminus
MonitorElement * TauSpinEffectsH_rhorhoAcoplanarityminus
Definition: TauValidation.h:92
TauValidation::analyze
virtual void analyze(edm::Event const &, edm::EventSetup const &) override
Definition: TauValidation.cc:156
reco::CompositeRefCandidateT::numberOfMothers
virtual size_t numberOfMothers() const
number of mothers
TauValidation::TauSpinEffectsZ_X50to75
MonitorElement * TauSpinEffectsZ_X50to75
Definition: TauValidation.h:92
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:48
p4
double p4[4]
Definition: TauolaWrapper.h:92
reco::CompositeRefCandidateT::numberOfDaughters
virtual size_t numberOfDaughters() const
number of daughters
reco::CompositeRefCandidateT::daughter
virtual const Candidate * daughter(size_type) const
return daughter at a given position, i = 0, ... numberOfDaughters() - 1 (read only mode) ...
TauValidation::TauSpinEffectsH_Xf
MonitorElement * TauSpinEffectsH_Xf
Definition: TauValidation.h:92
TauValidation::TauFSRPhotonsPt
MonitorElement * TauFSRPhotonsPt
Definition: TauValidation.h:92
funct::abs
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
j
int j
Definition: DBlmapReader.cc:9
TauValidation::TauSpinEffectsZ_muX
MonitorElement * TauSpinEffectsZ_muX
Definition: TauValidation.h:92
TauValidation::TauDecayChannels
MonitorElement * TauDecayChannels
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsH_Zs
MonitorElement * TauSpinEffectsH_Zs
Definition: TauValidation.h:92
TauValidation::NMODEID
unsigned int NMODEID
Definition: TauValidation.h:104
TauValidation::TauFSRPhotonsPtSum
MonitorElement * TauFSRPhotonsPtSum
Definition: TauValidation.h:92
ParticleData
HepPDT::ParticleData ParticleData
Definition: ParticleDataTable.h:9
TauValidation::tauMother
int tauMother(const reco::GenParticle *, double weight)
Definition: TauValidation.cc:334
TauValidation::bookHistograms
virtual void bookHistograms(DQMStore::IBooker &i, edm::Run const &, edm::EventSetup const &) override
Definition: TauValidation.cc:35
TauDecay_GenParticle::Get_TauDecayProducts
std::vector< const reco::GenParticle * > Get_TauDecayProducts()
Definition: TauDecay_GenParticle.h:31
TauValidation::GetMothers
const std::vector< const reco::GenParticle * > GetMothers(const reco::GenParticle *boson)
Definition: TauValidation.cc:264
TauValidation::TauSpinEffectsZ_MVis
MonitorElement * TauSpinEffectsZ_MVis
Definition: TauValidation.h:92
reco::Candidate::pdgId
virtual int pdgId() const =0
PDG identifier.
TauValidation::TauSpinEffectsZ_X75to88
MonitorElement * TauSpinEffectsZ_X75to88
Definition: TauValidation.h:92
TauValidation::visibleTauEnergy
double visibleTauEnergy(const reco::GenParticle *)
Definition: TauValidation.cc:699
TauDecay::DecayMode
static std::string DecayMode(unsigned int &MODE_ID)
Definition: TauDecay.cc:258
TauValidation::tauDecayChannel
int tauDecayChannel(const reco::GenParticle *tau, int jak_id, unsigned int TauBitMask, double weight)
Definition: TauValidation.cc:354
DQMStore::IBooker::setCurrentFolder
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:274
TauValidation::TauSpinEffectsH_pipiAcollinearity
MonitorElement * TauSpinEffectsH_pipiAcollinearity
Definition: TauValidation.h:92
reco::LeafCandidate::px
virtual double px() const
x coordinate of momentum vector
Definition: LeafCandidate.h:172
part
part
Definition: HCALResponse.h:20
sysUtil.pid
tuple pid
Definition: sysUtil.py:22
TauValidation::TauSpinEffectsZ_X
MonitorElement * TauSpinEffectsZ_X
Definition: TauValidation.h:92
TauValidation::fPDGTable
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
PDT table.
Definition: TauValidation.h:89
reco::LeafCandidate::pz
virtual double pz() const
z coordinate of momentum vector
Definition: LeafCandidate.h:176
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
TauValidation::TauSpinEffectsH_eX
MonitorElement * TauSpinEffectsH_eX
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsW_eX
MonitorElement * TauSpinEffectsW_eX
Definition: TauValidation.h:92
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
TauValidation::genparticleCollectionToken_
edm::EDGetTokenT< reco::GenParticleCollection > genparticleCollectionToken_
Definition: TauValidation.h:111
TauValidation::Zstoa
double Zstoa(double zs)
Definition: TauValidation.cc:669
TauValidation::findTauList
void findTauList(const reco::GenParticle *tau, std::vector< const reco::GenParticle * > &TauList)
Definition: TauValidation.cc:285
TauValidation::findFSRandBrem
void findFSRandBrem(const reco::GenParticle *p, bool doBrem, std::vector< const reco::GenParticle * > &ListofFSR, std::vector< const reco::GenParticle * > &ListofBrem)
Definition: TauValidation.cc:295
TauValidation::LifeTime
MonitorElement * LifeTime
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsH_muX
MonitorElement * TauSpinEffectsH_muX
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsW_UpsilonA1
MonitorElement * TauSpinEffectsW_UpsilonA1
Definition: TauValidation.h:92
a
double a
Definition: hdecay.h:121
TauValidation::countParticles
void countParticles(const reco::GenParticle *p, int &allCount, int &eCount, int &muCount, int &pi0Count, int &piCount, int &rhoCount, int &a1Count, int &KCount, int &KstarCount)
Definition: TauValidation.cc:388
TauValidation::TauEta
MonitorElement * TauEta
Definition: TauValidation.h:92
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
TauValidation::genparticleCollection_
edm::InputTag genparticleCollection_
Definition: TauValidation.h:86
TauValidation::TauSpinEffectsZ_Zs
MonitorElement * TauSpinEffectsZ_Zs
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsHpm_X
MonitorElement * TauSpinEffectsHpm_X
Definition: TauValidation.h:92
benchmark_cfg.pdgId
tuple pdgId
Definition: benchmark_cfg.py:18
histoStyle.weight
int weight
Definition: histoStyle.py:50
TauValidation::TauSpinEffectsHpm_muX
MonitorElement * TauSpinEffectsHpm_muX
Definition: TauValidation.h:92
TauDecay_GenParticle::AnalyzeTau
bool AnalyzeTau(const reco::GenParticle *Tau, unsigned int &MODE_ID, unsigned int &TauBitMask, bool dores, bool dopi0)
Definition: TauDecay_GenParticle.cc:17
TauValidation::TauSpinEffectsZ_X120UP
MonitorElement * TauSpinEffectsZ_X120UP
Definition: TauValidation.h:92
TauValidation::TauSpinEffectsH_Xb
MonitorElement * TauSpinEffectsH_Xb
Definition: TauValidation.h:92
reco::CompositeRefCandidateT::mother
virtual const Candidate * mother(size_type=0) const
return mother at a given position, i = 0, ... numberOfMothers() - 1 (read only mode) ...
TauValidation::nPrimeTaus
MonitorElement * nPrimeTaus
Definition: TauValidation.h:91
TauValidation::TauBremPhotonsPt
MonitorElement * TauBremPhotonsPt
Definition: TauValidation.h:92
reco::LeafCandidate::py
virtual double py() const
y coordinate of momentum vector
Definition: LeafCandidate.h:174
TauValidation::spinEffectsZH
void spinEffectsZH(const reco::GenParticle *boson, double weight)
Definition: TauValidation.cc:456
edm::Run
Definition: Run.h:41
TauDecay::nProng
unsigned int nProng(unsigned int &TauBitMask)
Definition: TauDecay.cc:218
TauValidation::TauSpinEffectsH_X
MonitorElement * TauSpinEffectsH_X
Definition: TauValidation.h:92
TauValidation::leadingPionP4
TLorentzVector leadingPionP4(const reco::GenParticle *)
Definition: TauValidation.cc:682
Generated for CMSSW Reference Manual by doxygen 1.8.5