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 
8 #include "CLHEP/Units/defs.h"
9 #include "CLHEP/Units/PhysicalConstants.h"
10 
11 #include "DataFormats/Math/interface/LorentzVector.h"
12 
13 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
14 #include "Validation/EventGenerator/interface/TauDecay_CMSSW.h"
15 #include "Validation/EventGenerator/interface/PdtPdgMini.h"
16 
17 using namespace edm;
18 
19 TauValidation::TauValidation(const edm::ParameterSet& iPSet):
20  _wmanager(iPSet)
21  ,hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection"))
22  ,tauEtCut(iPSet.getParameter<double>("tauEtCutForRtau"))
23  ,NJAKID(22)
24  ,zsbins(20)
25  ,zsmin(-0.5)
26  ,zsmax(0.5)
27  ,n_taupinu(200,0)
28 {
29  dbe = 0;
30  dbe = edm::Service<DQMStore>().operator->();
31 }
32 
33 TauValidation::~TauValidation() {}
34 
35 void TauValidation::beginJob()
36 {
37  if(dbe){
39  dbe->setCurrentFolder("Generator/Tau");
40 
41  // Number of analyzed events
42  nTaus = dbe->book1D("nTaus", "n analyzed Taus", 1, 0., 1.);
43  nPrimeTaus = dbe->book1D("nPrimeTaus", "n analyzed prime Taus", 1, 0., 1.);
44 
45  //Kinematics
46  TauPt = dbe->book1D("TauPt","Tau pT", 100 ,0,100);
47  TauEta = dbe->book1D("TauEta","Tau eta", 100 ,-2.5,2.5);
48  TauPhi = dbe->book1D("TauPhi","Tau phi", 100 ,-3.14,3.14);
49  TauProngs = dbe->book1D("TauProngs","Tau n prongs", 7 ,0,7);
50  TauDecayChannels = dbe->book1D("TauDecayChannels","Tau decay channels", 13 ,0,13);
51  TauDecayChannels->setBinLabel(1+undetermined,"?");
52  TauDecayChannels->setBinLabel(1+electron,"e");
53  TauDecayChannels->setBinLabel(1+muon,"mu");
54  TauDecayChannels->setBinLabel(1+pi,"#pi^{#pm}");
55  TauDecayChannels->setBinLabel(1+rho,"#rho^{#pm}");
56  TauDecayChannels->setBinLabel(1+a1,"a_{1}^{#pm}");
57  TauDecayChannels->setBinLabel(1+pi1pi0,"#pi^{#pm}#pi^{0}");
58  TauDecayChannels->setBinLabel(1+pinpi0,"#pi^{#pm}n#pi^{0}");
59  TauDecayChannels->setBinLabel(1+tripi,"3#pi^{#pm}");
60  TauDecayChannels->setBinLabel(1+tripinpi0,"3#pi^{#pm}n#pi^{0}");
61  TauDecayChannels->setBinLabel(1+K,"K");
62  TauDecayChannels->setBinLabel(1+Kstar,"K^{*}");
63  TauDecayChannels->setBinLabel(1+stable,"Stable");
64 
65  TauMothers = dbe->book1D("TauMothers","Tau mother particles", 10 ,0,10);
66  TauMothers->setBinLabel(1+other,"?");
67  TauMothers->setBinLabel(1+B,"B Decays");
68  TauMothers->setBinLabel(1+D,"D Decays");
69  TauMothers->setBinLabel(1+gamma,"#gamma");
70  TauMothers->setBinLabel(1+Z,"Z");
71  TauMothers->setBinLabel(1+W,"W");
72  TauMothers->setBinLabel(1+HSM,"H_{SM}/h^{0}");
73  TauMothers->setBinLabel(1+H0,"H^{0}");
74  TauMothers->setBinLabel(1+A0,"A^{0}");
75  TauMothers->setBinLabel(1+Hpm,"H^{#pm}");
76 
77  DecayLength = dbe->book1D("DecayLength","#tau Decay Length", 100 ,0,20); DecayLength->setAxisTitle("L_{#tau} (mm)");
78  LifeTime = dbe->book1D("LifeTime","#tau LifeTime ", 100 ,0,1000E-15); LifeTime->setAxisTitle("#tau_{#tau}s)");
79 
80  TauRtauW = dbe->book1D("TauRtauW","W->Tau p(leading track)/E(visible tau)", 50 ,0,1); TauRtauW->setAxisTitle("rtau");
81  TauRtauHpm = dbe->book1D("TauRtauHpm","Hpm->Tau p(leading track)/E(visible tau)", 50 ,0,1); TauRtauHpm->setAxisTitle("rtau");
82 
83  TauSpinEffectsW_X = dbe->book1D("TauSpinEffectsWX","Pion energy in W rest frame", 50 ,0,1); TauSpinEffectsW_X->setAxisTitle("X");
84  TauSpinEffectsHpm_X = dbe->book1D("TauSpinEffectsHpmX","Pion energy in Hpm rest frame", 50 ,0,1); TauSpinEffectsHpm_X->setAxisTitle("X");
85 
86  TauSpinEffectsW_eX = dbe->book1D("TauSpinEffectsWeX","e energy in W rest frame", 50 ,0,1); TauSpinEffectsW_eX->setAxisTitle("X");
87  TauSpinEffectsHpm_eX = dbe->book1D("TauSpinEffectsHpmeX","e energy in Hpm rest frame", 50 ,0,1); TauSpinEffectsHpm_eX->setAxisTitle("X");
88 
89  TauSpinEffectsW_muX = dbe->book1D("TauSpinEffectsWmuX","mu energy in W rest frame", 50 ,0,1); TauSpinEffectsW_muX->setAxisTitle("X");
90  TauSpinEffectsHpm_muX = dbe->book1D("TauSpinEffectsHpmmuX","mu energy in Hpm rest frame", 50 ,0,1); TauSpinEffectsHpm_muX->setAxisTitle("X");
91 
92  TauSpinEffectsW_UpsilonRho = dbe->book1D("TauSpinEffectsWUpsilonRho","#Upsilon for #rho", 50 ,-1,1); TauSpinEffectsW_UpsilonRho->setAxisTitle("#Upsilon");
93  TauSpinEffectsHpm_UpsilonRho = dbe->book1D("TauSpinEffectsHpmUpsilonRho","#Upsilon for #rho", 50 ,-1,1); TauSpinEffectsHpm_UpsilonRho->setAxisTitle("#Upsilon");
94 
95  TauSpinEffectsW_UpsilonA1 = dbe->book1D("TauSpinEffectsWUpsilonA1","#Upsilon for a1", 50 ,-1,1); TauSpinEffectsW_UpsilonA1->setAxisTitle("#Upsilon");
96  TauSpinEffectsHpm_UpsilonA1 = dbe->book1D("TauSpinEffectsHpmUpsilonA1","#Upsilon for a1", 50 ,-1,1); TauSpinEffectsHpm_UpsilonA1->setAxisTitle("#Upsilon");
97 
98  TauSpinEffectsH_pipiAcoplanarity = dbe->book1D("TauSpinEffectsH_pipiAcoplanarity","H Acoplanarity for #pi^{-}#pi^{+}", 50 ,0,2*TMath::Pi()); TauSpinEffectsH_pipiAcoplanarity->setAxisTitle("Acoplanarity");
99 
100  TauSpinEffectsH_pipiAcollinearity = dbe->book1D("TauSpinEffectsH_pipiAcollinearity","H Acollinearity for #pi^{-}#pi^{+}", 50 ,0,TMath::Pi()); TauSpinEffectsH_pipiAcollinearity->setAxisTitle("Acollinearity");
101  TauSpinEffectsH_pipiAcollinearityzoom = dbe->book1D("TauSpinEffectsH_pipiAcollinearityzoom","H Acollinearity for #pi^{-}#pi^{+}", 50 ,3,TMath::Pi()); TauSpinEffectsH_pipiAcollinearityzoom->setAxisTitle("Acollinearity");
102 
103  TauSpinEffectsZ_MVis = dbe->book1D("TauSpinEffectsZMVis","Mass of pi+ pi-", 25 ,0,1.1); TauSpinEffectsZ_MVis->setAxisTitle("M_{#pi^{+}#pi^{-}}");
104  TauSpinEffectsH_MVis = dbe->book1D("TauSpinEffectsHMVis","Mass of pi+ pi-", 25 ,0,1.1); TauSpinEffectsZ_MVis->setAxisTitle("M_{#pi^{+}#pi^{-}}");
105 
106  TauSpinEffectsZ_Zs = dbe->book1D("TauSpinEffectsZZs","Z_{s}", zsbins ,zsmin,zsmax); TauSpinEffectsZ_Zs->setAxisTitle("Z_{s}");
107  TauSpinEffectsH_Zs = dbe->book1D("TauSpinEffectsHZs","Z_{s}", zsbins ,zsmin,zsmax); TauSpinEffectsZ_Zs->setAxisTitle("Z_{s}");
108 
109  TauSpinEffectsZ_X= dbe->book1D("TauSpinEffectsZX","X of #tau^{-}", 25 ,0,1.0); TauSpinEffectsZ_X->setAxisTitle("X");
110  TauSpinEffectsH_X= dbe->book1D("TauSpinEffectsH_X","X of #tau^{-}", 25 ,0,1.0); TauSpinEffectsH_X->setAxisTitle("X");
111 
112  TauSpinEffectsZ_Xf = dbe->book1D("TauSpinEffectsZXf","X of forward emitted #tau^{-}", 25 ,0,1.0); TauSpinEffectsZ_Xf->setAxisTitle("X_{f}");
113  TauSpinEffectsH_Xf = dbe->book1D("TauSpinEffectsHXf","X of forward emitted #tau^{-}", 25 ,0,1.0); TauSpinEffectsZ_Xf->setAxisTitle("X_{f}");
114 
115  TauSpinEffectsZ_Xb = dbe->book1D("TauSpinEffectsZXb","X of backward emitted #tau^{-}", 25 ,0,1.0); TauSpinEffectsZ_Xb->setAxisTitle("X_{b}");
116  TauSpinEffectsH_Xb = dbe->book1D("TauSpinEffectsHXb","X of backward emitted #tau^{-}", 25 ,0,1.0); TauSpinEffectsZ_Xb->setAxisTitle("X_{b}");
117 
118  TauSpinEffectsZ_eX = dbe->book1D("TauSpinEffectsZeX","e energy in Z rest frame", 50 ,0,1); TauSpinEffectsZ_eX->setAxisTitle("X");
119  TauSpinEffectsH_eX = dbe->book1D("TauSpinEffectsHeX","e energy in H rest frame", 50 ,0,1); TauSpinEffectsH_eX->setAxisTitle("X");
120 
121  TauSpinEffectsZ_muX = dbe->book1D("TauSpinEffectsZmuX","mu energy in z rest frame", 50 ,0,1); TauSpinEffectsZ_muX->setAxisTitle("X");
122  TauSpinEffectsH_muX = dbe->book1D("TauSpinEffectsHmuX","mu energy in H rest frame", 50 ,0,1); TauSpinEffectsH_muX->setAxisTitle("X");
123 
124  TauSpinEffectsH_rhorhoAcoplanarityminus = dbe->book1D("TauSpinEffectsH_rhorhoAcoplanarityminus","#phi^{*-} (acoplanarity) for Higgs #rightarrow #rho-#rho (y_{1}*y_{2}<0)", 32 ,0,2*TMath::Pi()); TauSpinEffectsH_rhorhoAcoplanarityminus->setAxisTitle("#phi^{*-} (Acoplanarity)");
125  TauSpinEffectsH_rhorhoAcoplanarityplus = dbe->book1D("TauSpinEffectsH_rhorhoAcoplanarityplus","#phi^{*+} (acoplanarity) for Higgs #rightarrow #rho-#rho (y_{1}*y_{2}>0)", 32 ,0,2*TMath::Pi()); TauSpinEffectsH_rhorhoAcoplanarityplus->setAxisTitle("#phi^{*+} (Acoplanarity)");
126 
127 
128 
129  TauSpinEffectstautau_polvxM = dbe->book1D("TauSpinEffectsZ_polvxM","#tau Polarization ", 40,0,200); TauSpinEffectstautau_polvxM->setAxisTitle("Mass (GeV)");
130 
131 
132  TauFSRPhotonsN=dbe->book1D("TauFSRPhotonsN","FSR Photons radiating from/with tau (Gauge Boson)", 5 ,-0.5,4.5);
133  TauFSRPhotonsN->setAxisTitle("N FSR Photons radiating from/with tau");
134  TauFSRPhotonsPt=dbe->book1D("TauFSRPhotonsPt","Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
135  TauFSRPhotonsPt->setAxisTitle("P_{t} of FSR Photons radiating from/with tau [per tau]");
136  TauFSRPhotonsPtSum=dbe->book1D("TauFSRPhotonsPtSum","Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
137  TauFSRPhotonsPtSum->setAxisTitle("P_{t} of FSR Photons radiating from/with tau [per tau]");
138 
139  TauBremPhotonsN=dbe->book1D("TauBremPhotonsN","Brem. Photons radiating in tau decay", 5 ,-0.5,4.5);
140  TauBremPhotonsN->setAxisTitle("N FSR Photons radiating from/with tau");
141  TauBremPhotonsPt=dbe->book1D("TauBremPhotonsPt","Sum Brem Pt ", 100 ,0,100);
142  TauFSRPhotonsPt->setAxisTitle("P_{t} of Brem. Photons radiating in tau decay");
143  TauBremPhotonsPtSum =dbe->book1D("TauBremPhotonsPtSum","Sum of Brem Pt ", 100 ,0,100);
144  TauFSRPhotonsPtSum->setAxisTitle("Sum P_{t} of Brem. Photons radiating in tau decay");
145 
146  JAKID =dbe->book1D("JAKID","JAK ID",NJAKID+1,-0.5,NJAKID+0.5);
147  for(unsigned int i=0; i<NJAKID+1;i++){
148  JAKInvMass.push_back(std::vector<MonitorElement *>());
149  TString tmp="JAKID";
150  tmp+=i;
151  JAKInvMass.at(i).push_back(dbe->book1D("M"+tmp,"M_{"+tmp+"} (GeV)", 80 ,0,2.0));
152  if(i==TauDecay::JAK_A1_3PI ||
153  i==TauDecay::JAK_KPIK ||
154  i==TauDecay::JAK_KPIPI ){
155  JAKInvMass.at(i).push_back(dbe->book1D("M13"+tmp,"M_{13,"+tmp+"} (GeV)", 80 ,0,2.0));
156  JAKInvMass.at(i).push_back(dbe->book1D("M23"+tmp,"M_{23,"+tmp+"} (GeV)", 80 ,0,2.0));
157  JAKInvMass.at(i).push_back(dbe->book1D("M12"+tmp,"M_{12,"+tmp+"} (GeV)", 80 ,0,2.0));
158  }
159  }
160  }
161 
162  return;
163 }
164 
165 void TauValidation::endJob(){
166  return;
167 }
168 
169 void TauValidation::beginRun(const edm::Run& iRun,const edm::EventSetup& iSetup)
170 {
172  iSetup.getData( fPDGTable );
173  return;
174 }
175 void TauValidation::endRun(const edm::Run& iRun,const edm::EventSetup& iSetup){return;}
176 void TauValidation::analyze(const edm::Event& iEvent,const edm::EventSetup& iSetup)
177 {
179  edm::Handle<HepMCProduct> evt;
180  iEvent.getByLabel(hepmcCollection_, evt);
181 
182  //Get EVENT
183  HepMC::GenEvent *myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
184 
185  double weight = _wmanager.weight(iEvent);
186 
188  /*
189  edm::Handle<double> WT;
190  iEvent.getByLabel(edm::InputTag("TauSpinnerGen","TauSpinnerWT"),WT);
191  weight = 1.0;
192  if(*(WT.product())>1e-3 && *(WT.product())<=10.0) weight=(*(WT.product()));
193  else {weight=1.0;}
194  */
196 
197  // find taus
198  for(HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin(); iter != myGenEvent->particles_end(); iter++) {
199  if(abs((*iter)->pdg_id())==PdtPdgMini::Z0 || abs((*iter)->pdg_id())==PdtPdgMini::Higgs0){
200  spinEffectsZH(*iter,weight);
201  }
202  if(abs((*iter)->pdg_id())==15){
203  if(isLastTauinChain(*iter)){
204  nTaus->Fill(0.5,weight);
205  int mother = tauMother(*iter,weight);
206  int decaychannel = tauDecayChannel(*iter,weight);
207  tauProngs(*iter, weight);
208  if(mother>-1){ // exclude B, D and other non-signal decay modes
209  nPrimeTaus->Fill(0.5,weight);
210  TauPt->Fill((*iter)->momentum().perp(),weight);
211  TauEta->Fill((*iter)->momentum().eta(),weight);
212  TauPhi->Fill((*iter)->momentum().phi(),weight);
213  rtau(*iter,mother,decaychannel,weight);
214  photons(*iter,weight);
215  }
217  //Adding JAKID and Mass information
218  //
219  TauDecay_CMSSW TD;
220  unsigned int jak_id, TauBitMask;
221  if(TD.AnalyzeTau((*iter),jak_id,TauBitMask,false,false)){
222  JAKID->Fill(jak_id,weight);
223  if(jak_id<=NJAKID){
224  int tcharge=(*iter)->pdg_id()/abs((*iter)->pdg_id());
225  std::vector<HepMC::GenParticle*> part=TD.Get_TauDecayProducts();
226  spinEffectsWHpm(*iter,mother,jak_id,part,weight);
227  TLorentzVector LVQ(0,0,0,0);
228  TLorentzVector LVS12(0,0,0,0);
229  TLorentzVector LVS13(0,0,0,0);
230  TLorentzVector LVS23(0,0,0,0);
231  bool haspart1=false;
232  TVector3 PV,SV;
233  bool hasDL(false);
234  for(unsigned int i=0;i<part.size();i++){
235  if(abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_tau && TD.isTauFinalStateParticle(part.at(i)->pdg_id()) && !hasDL){
236  PV=TVector3((*iter)->production_vertex()->point3d().x(),(*iter)->production_vertex()->point3d().y(),(*iter)->production_vertex()->point3d().z());
237  SV=TVector3(part.at(i)->production_vertex()->point3d().x(),part.at(i)->production_vertex()->point3d().y(),part.at(i)->production_vertex()->point3d().z());
238  TVector3 DL=SV-PV;
239  DecayLength->Fill(DL.Mag(),weight);
240  double c(2.99792458E8),Ltau(DL.Mag()/1000),beta((*iter)->momentum().rho()/(*iter)->momentum().m());
241  LifeTime->Fill( Ltau/(c*beta),weight);
242  }
243 
244  if(TD.isTauFinalStateParticle(part.at(i)->pdg_id()) &&
245  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_e &&
246  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_mu &&
247  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_tau ){
248  TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
249  LVQ+=LV;
250  if(jak_id==TauDecay::JAK_A1_3PI ||
251  jak_id==TauDecay::JAK_KPIK ||
252  jak_id==TauDecay::JAK_KPIPI
253  ){
254  if((tcharge==part.at(i)->pdg_id()/abs(part.at(i)->pdg_id()) && TD.nProng(TauBitMask)==3) || (jak_id==TauDecay::JAK_A1_3PI && TD.nProng(TauBitMask)==1 && abs(part.at(i)->pdg_id())==PdtPdgMini::pi_plus) ){
255  LVS13+=LV;
256  LVS23+=LV;
257  }
258  else{
259  LVS12+=LV;
260  if(!haspart1 && ((jak_id==TauDecay::JAK_A1_3PI) || (jak_id!=TauDecay::JAK_A1_3PI && abs(part.at(i)->pdg_id())==PdtPdgMini::K_plus) )){
261  LVS13+=LV;
262  haspart1=true;
263  }
264  else{
265  LVS23+=LV;
266  }
267  }
268  }
269  }
270  }
271  part.clear();
272  JAKInvMass.at(jak_id).at(0)->Fill(LVQ.M(),weight);
273  if(jak_id==TauDecay::JAK_A1_3PI ||
274  jak_id==TauDecay::JAK_KPIK ||
275  jak_id==TauDecay::JAK_KPIPI
276  ){
277  JAKInvMass.at(jak_id).at(1)->Fill(LVS13.M(),weight);
278  JAKInvMass.at(jak_id).at(2)->Fill(LVS23.M(),weight);
279  JAKInvMass.at(jak_id).at(3)->Fill(LVS12.M(),weight);
280  }
281  }
282  }
283  else{
284  JAKID->Fill(jak_id,weight);
285  }
286  }
287  }
288  }
289  delete myGenEvent;
290 }//analyze
291 
292 const HepMC::GenParticle* TauValidation::GetMother(const HepMC::GenParticle* tau){
293  if ( tau->production_vertex() ) {
294  HepMC::GenVertex::particle_iterator mother;
295  for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
296  if((*mother)->pdg_id() == tau->pdg_id()) return GetMother(*mother);
297  return (*mother);
298  }
299  }
300  return tau;
301 }
302 
303 
304 const std::vector<HepMC::GenParticle*> TauValidation::GetMothers(const HepMC::GenParticle* boson){
305  std::vector<HepMC::GenParticle*> mothers;
306  if ( boson->production_vertex() ) {
307  HepMC::GenVertex::particle_iterator mother;
308  for (mother = boson->production_vertex()->particles_begin(HepMC::parents); mother!= boson->production_vertex()->particles_end(HepMC::parents); mother++ ) {
309  if((*mother)->pdg_id() == boson->pdg_id()) return GetMothers(*mother);
310  mothers.push_back(*mother);
311  }
312  }
313  return mothers;
314 }
315 
316 
317 int TauValidation::findMother(const HepMC::GenParticle* tau){
318  int mother_pid = 0;
319  if ( tau->production_vertex() ) {
320  HepMC::GenVertex::particle_iterator mother;
321  for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
322  mother_pid = (*mother)->pdg_id();
323  if(mother_pid == tau->pdg_id()) return findMother(*mother); //mother_pid = -1; Make recursive to look for last tau in chain
324  }
325  }
326  return mother_pid;
327 }
328 
329 
330 bool TauValidation::isLastTauinChain(const HepMC::GenParticle* tau){
331  if ( tau->end_vertex() ) {
332  HepMC::GenVertex::particle_iterator dau;
333  for (dau = tau->end_vertex()->particles_begin(HepMC::children); dau!= tau->end_vertex()->particles_end(HepMC::children); dau++ ) {
334  int dau_pid = (*dau)->pdg_id();
335  if(dau_pid == tau->pdg_id()) return false;
336  }
337  }
338  return true;
339 }
340 
341 
342 void TauValidation::findTauList(const HepMC::GenParticle* tau,std::vector<const HepMC::GenParticle*> &TauList){
343  TauList.insert(TauList.begin(),tau);
344  if ( tau->production_vertex() ) {
345  HepMC::GenVertex::particle_iterator mother;
346  for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents);mother++) {
347  if((*mother)->pdg_id() == tau->pdg_id()){
348  findTauList(*mother,TauList);
349  }
350  }
351  }
352 }
353 
354 void TauValidation::findFSRandBrem(const HepMC::GenParticle* p, bool doBrem, std::vector<const HepMC::GenParticle*> &ListofFSR,
355  std::vector<const HepMC::GenParticle*> &ListofBrem){
356  // 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.
357  if(abs(p->pdg_id())==15){
358  if(isLastTauinChain(p)){ doBrem=true;}
359  else{ doBrem=false;}
360  }
361  int photo_ID=22;
362  if ( p->end_vertex() ) {
363  HepMC::GenVertex::particle_iterator dau;
364  for (dau = p->end_vertex()->particles_begin(HepMC::children); dau!= p->end_vertex()->particles_end(HepMC::children); dau++ ) {
365  //if(doBrem) std::cout << "true " << (*dau)->pdg_id() << " " << std::endl;
366  //else std::cout << "false " << (*dau)->pdg_id() << " " << std::endl;
367  if(abs((*dau)->pdg_id()) == abs(photo_ID) && !doBrem){ListofFSR.push_back(*dau);}
368  if(abs((*dau)->pdg_id()) == abs(photo_ID) && doBrem){ListofBrem.push_back(*dau);}
369  if((*dau)->end_vertex() && (*dau)->end_vertex()->particles_out_size()>0 && abs((*dau)->pdg_id()) != 111 && abs((*dau)->pdg_id()) != 221/* remove pi0 and eta decays*/){
370  findFSRandBrem(*dau,doBrem,ListofFSR,ListofBrem);
371  }
372  }
373  }
374 }
375 
376 
377 
378 void TauValidation::FindPhotosFSR(const HepMC::GenParticle* p,std::vector<const HepMC::GenParticle*> &ListofFSR,double &BosonScale){
379  BosonScale=0.0;
380  const HepMC::GenParticle* m=GetMother(p);
381  double mother_pid=m->pdg_id();
382  if(m->end_vertex() && mother_pid!=p->pdg_id()){
383  HepMC::GenVertex::particle_iterator dau;
384  for (dau = m->end_vertex()->particles_begin(HepMC::children); dau!= m->end_vertex()->particles_end(HepMC::children); dau++ ) {
385  int dau_pid = (*dau)->pdg_id();
386  if(fabs(dau_pid) == 22) {
387  ListofFSR.push_back(*dau);
388  }
389  }
390  }
391  if(abs(mother_pid) == 24) BosonScale=1.0; // W
392  if(abs(mother_pid) == 23) BosonScale=2.0; // Z;
393  if(abs(mother_pid) == 22) BosonScale=2.0; // gamma;
394  if(abs(mother_pid) == 25) BosonScale=2.0; // HSM;
395  if(abs(mother_pid) == 35) BosonScale=2.0; // H0;
396  if(abs(mother_pid) == 36) BosonScale=2.0; // A0;
397  if(abs(mother_pid) == 37) BosonScale=1.0; //Hpm;
398 }
399 
400 
401 int TauValidation::tauMother(const HepMC::GenParticle* tau, double weight){
402 
403  if(abs(tau->pdg_id()) != 15 ) return -3;
404 
405  int mother_pid = findMother(tau);
406  if(mother_pid == -2) return -2;
407 
408  int label = other;
409  if(abs(mother_pid) == 24) label = W;
410  if(abs(mother_pid) == 23) label = Z;
411  if(abs(mother_pid) == 22) label = gamma;
412  if(abs(mother_pid) == 25) label = HSM;
413  if(abs(mother_pid) == 35) label = H0;
414  if(abs(mother_pid) == 36) label = A0;
415  if(abs(mother_pid) == 37) label = Hpm;
416 
417  int mother_shortpid=(abs(mother_pid)%10000);
418  if(mother_shortpid>500 && mother_shortpid<600 )label = B;
419  if(mother_shortpid>400 && mother_shortpid<500)label = D;
420  TauMothers->Fill(label,weight);
421  if(label==B || label == D || label == other) return -1;
422 
423  return mother_pid;
424 }
425 
426 int TauValidation::tauProngs(const HepMC::GenParticle* tau, double weight){
427  int nProngs = 0;
428  if ( tau->end_vertex() ) {
429  HepMC::GenVertex::particle_iterator des;
430  for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
431  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
432  int pid = (*des)->pdg_id();
433  if(abs(pid) == 15) return tauProngs(*des, weight);
434  if((*des)->status() != 1) continue; // dont count unstable particles
435 
436  const HepPDT::ParticleData* pd = fPDGTable->particle((*des)->pdg_id ());
437  int charge = (int) pd->charge();
438  if(charge == 0) continue;
439  nProngs++;
440  }
441  }
442  TauProngs->Fill(nProngs,weight);
443  return nProngs;
444 }
445 
446 int TauValidation::tauDecayChannel(const HepMC::GenParticle* tau, double weight){
447 
448  int channel = undetermined;
449  if(tau->status() == 1) channel = stable;
450  int allCount = 0,
451  eCount = 0,
452  muCount = 0,
453  pi0Count = 0,
454  piCount = 0,
455  rhoCount = 0,
456  a1Count = 0,
457  KCount = 0,
458  KstarCount = 0;
459 
460  if ( tau->end_vertex() ) {
461  HepMC::GenVertex::particle_iterator des;
462  for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
463  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
464  int pid = (*des)->pdg_id();
465  if(abs(pid) == 15) return tauDecayChannel(*des,weight);
466 
467  allCount++;
468  if(abs(pid) == 11) eCount++;
469  if(abs(pid) == 13) muCount++;
470  if(abs(pid) == 111) pi0Count++;
471  if(abs(pid) == 211) piCount++;
472  if(abs(pid) == 213) rhoCount++;
473  if(abs(pid) == 20213) a1Count++;
474  if(abs(pid) == 321) KCount++;
475  if(abs(pid) == 323) KstarCount++;
476 
477  }
478  }
479  // resonances
480  if(KCount >= 1) channel = K;
481  if(KstarCount >= 1) channel = Kstar;
482  if(a1Count >= 1) channel = a1;
483  if(rhoCount >= 1) channel = rho;
484  if(channel!=undetermined && weight!=0.0) TauDecayChannels->Fill(channel,weight);
485 
486  // final state products
487  if(piCount == 1 && pi0Count == 0) channel = pi;
488  if(piCount == 1 && pi0Count == 1) channel = pi1pi0;
489  if(piCount == 1 && pi0Count > 1) channel = pinpi0;
490  if(piCount == 3 && pi0Count == 0) channel = tripi;
491  if(piCount == 3 && pi0Count > 0) channel = tripinpi0;
492  if(eCount == 1) channel = electron;
493  if(muCount == 1) channel = muon;
494  if(weight!=0.0) TauDecayChannels->Fill(channel,weight);
495  return channel;
496 }
497 
498 
499 void TauValidation::rtau(const HepMC::GenParticle* tau,int mother, int decay, double weight){
500 
501  if(decay != pi1pi0) return; // polarization only for 1-prong hadronic taus with one neutral pion to make a clean case
502 
503  if(tau->momentum().perp() < tauEtCut) return; // rtau visible only for boosted taus
504 
505  double rTau = 0;
506  double ltrack = leadingPionMomentum(tau, weight);
507  double visibleTauE = visibleTauEnergy(tau);
508 
509  if(visibleTauE != 0) rTau = ltrack/visibleTauE;
510 
511  if(abs(mother) == 24) TauRtauW->Fill(rTau,weight);
512  if(abs(mother) == 37) TauRtauHpm->Fill(rTau,weight);
513 }
514 
515  void TauValidation::spinEffectsWHpm(const HepMC::GenParticle* tau,int mother, int decay, std::vector<HepMC::GenParticle*> &part,double weight){
516  if(decay == TauDecay::JAK_PION || decay == TauDecay::JAK_MUON || decay == TauDecay::JAK_ELECTRON){ // polarization only for 1-prong hadronic taus with no neutral pions
517  TLorentzVector momP4 = motherP4(tau);
518  TLorentzVector pionP4 = leadingPionP4(tau);
519  pionP4.Boost(-1*momP4.BoostVector());
520  double energy = pionP4.E()/(momP4.M()/2);
521  if(decay == TauDecay::JAK_PION){
522  if(abs(mother) == 24) TauSpinEffectsW_X->Fill(energy,weight);
523  if(abs(mother) == 37) TauSpinEffectsHpm_X->Fill(energy,weight);
524  }
525  if(decay == TauDecay::JAK_MUON){
526  if(abs(mother) == 24) TauSpinEffectsW_muX->Fill(energy,weight);
527  if(abs(mother) == 37) TauSpinEffectsHpm_muX->Fill(energy,weight);
528  }
529  if(decay == TauDecay::JAK_ELECTRON){
530  if(abs(mother) == 24) TauSpinEffectsW_eX->Fill(energy,weight);
531  if(abs(mother) == 37) TauSpinEffectsHpm_eX->Fill(energy,weight);
532  }
533 
534  }
535  else if(decay==TauDecay::JAK_RHO_PIPI0){
536  TLorentzVector rho(0,0,0,0),pi(0,0,0,0);
537  for(unsigned int i=0;i<part.size();i++){
538  TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
539  if(abs(part.at(i)->pdg_id())==PdtPdgMini::pi_plus){pi+=LV; rho+=LV;}
540  if(abs(part.at(i)->pdg_id())==PdtPdgMini::pi0){rho+=LV;}
541  }
542  if(abs(mother) == 24) TauSpinEffectsW_UpsilonRho->Fill(2*pi.Pt()/rho.Pt()-1,weight);
543  if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonRho->Fill(2*pi.Pt()/rho.Pt()-1,weight);
544  }
545  else if(decay==TauDecay::JAK_A1_3PI){ // only for pi2pi0 for now
546  TLorentzVector a1(0,0,0,0),pi_p(0,0,0,0),pi_m(0,0,0,0);
547  int nplus(0),nminus(0);
548  for(unsigned int i=0;i<part.size();i++){
549  TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
550  if(part.at(i)->pdg_id()==PdtPdgMini::pi_plus){ pi_p+=LV; a1+=LV; nplus++;}
551  if(part.at(i)->pdg_id()==PdtPdgMini::pi_minus){pi_m+=LV; a1+=LV; nminus++;}
552  }
553  double gamma=0;
554  if(nplus+nminus==3 && nplus==1) gamma=2*pi_p.Pt()/a1.Pt()-1;
555  if(nplus+nminus==3 && nminus==1) gamma=2*pi_m.Pt()/a1.Pt()-1;
556  else{
557  pi_p+=pi_m; gamma=2*pi_p.Pt()/a1.Pt()-1;
558  }
559  if(abs(mother) == 24) TauSpinEffectsW_UpsilonA1->Fill(gamma,weight);
560  if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonA1->Fill(gamma,weight);
561  }
562 }
563 
564  void TauValidation::spinEffectsZH(const HepMC::GenParticle* boson, double weight){
565  TLorentzVector tautau(0,0,0,0);
566  TLorentzVector pipi(0,0,0,0);
567  TLorentzVector taum(0,0,0,0);
568  TLorentzVector taup(0,0,0,0);
569  TLorentzVector rho_plus,rho_minus,pi_rhominus,pi0_rhominus,pi_rhoplus,pi0_rhoplus,pi_plus,pi_minus;
570  bool hasrho_minus(false),hasrho_plus(false),haspi_minus(false),haspi_plus(false);
571  int nSinglePionDecays(0),nSingleMuonDecays(0),nSingleElectronDecays(0);
572  double x1(0),x2(0);
573  TLorentzVector Zboson(boson->momentum().px(),boson->momentum().py(),boson->momentum().pz(),boson->momentum().e());
574  if ( boson->end_vertex() ) {
575  HepMC::GenVertex::particle_iterator des;
576  for(des = boson->end_vertex()->particles_begin(HepMC::children); des!= boson->end_vertex()->particles_end(HepMC::children);++des ) {
577  int pid = (*des)->pdg_id();
578  if(abs(findMother(*des)) != 15 && abs(pid) == 15 && (tauDecayChannel(*des) == pi || tauDecayChannel(*des) == muon || tauDecayChannel(*des) == electron )){
579  if(tauDecayChannel(*des) == pi)nSinglePionDecays++;
580  if(tauDecayChannel(*des) == muon)nSingleMuonDecays++;
581  if(tauDecayChannel(*des) == electron)nSingleElectronDecays++;
582  TLorentzVector LVtau((*des)->momentum().px(),(*des)->momentum().py(),(*des)->momentum().pz(),(*des)->momentum().e());
583  tautau += LVtau;
584  TLorentzVector LVpi=leadingPionP4(*des);
585  pipi+=LVpi;
586  const HepPDT::ParticleData* pd = fPDGTable->particle((*des)->pdg_id ());
587  int charge = (int) pd->charge();
588  LVtau.Boost(-1*Zboson.BoostVector());
589  LVpi.Boost(-1*Zboson.BoostVector());
590  if(tauDecayChannel(*des) == pi){
591  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_X->Fill(LVpi.P()/LVtau.E(),weight);
592  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_X->Fill(LVpi.P()/LVtau.E(),weight);
593  if(tautau.M()>2.0){
594  // for both z and h compute tau polarization
595  TH1* histo = TauSpinEffectstautau_polvxM->getTH1();
596  int n=histo->FindBin(tautau.M());
597  double mean=histo->GetBinContent(n);
598  double err=histo->GetBinError(n);
599  double x=(6.0*(LVpi.P()/LVtau.E()-0.5));
600  double M2=err*err*(n_taupinu[n]-1)*n_taupinu[n];
601  //use Knuth algorithm: Donald E. Knuth (1998). The Art of Computer Programming,
602  //volume 2: Seminumerical Algorithms, 3rd edn., p. 232. Boston: Addison-Wesley.
603  n_taupinu[n]+=1.0;
604  double d=x-mean;
605  mean+=d/n_taupinu[n];
606  //std::cout << n << M2 << " " << d << " " << (x-mean) << std::endl;
607  M2+=d*(x-mean);
608  if(n_taupinu[n]>1){err=sqrt((M2/(n_taupinu[n]-1))/n_taupinu[n]);}
609  else{ err=sqrt(M2);}
610  //save values
611  //std::cout << n << " " << tautau.M() << " " << mean << "+/-" << err << std::endl;
612  histo->SetBinContent(n,mean);
613  histo->SetBinError(n,err);
614  }
615  }
616  if(charge<0){x1=LVpi.P()/LVtau.E(); taum=LVtau;}
617  else{ x2=LVpi.P()/LVtau.E();}
618  }
619  if(abs(findMother(*des)) != 15 && abs(pid) == 15 && (tauDecayChannel(*des) == rho || tauDecayChannel(*des) == pi1pi0) ){
620  if ( (*des)->end_vertex() ) {
621  HepMC::GenVertex::particle_iterator tauprod;
622  TLorentzVector LVtau((*des)->momentum().px(),(*des)->momentum().py(),(*des)->momentum().pz(),(*des)->momentum().e());
623  if(pid == 15)taum=LVtau;
624  if(pid ==-15)taup=LVtau;
625  for(tauprod = (*des)->end_vertex()->particles_begin(HepMC::descendants); tauprod!= (*des)->end_vertex()->particles_end(HepMC::descendants);++tauprod ) {
626  int pid_d = (*tauprod)->pdg_id();
627  if(abs(pid_d)==211 || abs(pid_d)==111){
628  TLorentzVector LV((*tauprod)->momentum().px(),(*tauprod)->momentum().py(),(*tauprod)->momentum().pz(),(*tauprod)->momentum().e());
629  if(pid==15){
630  hasrho_minus=true;
631  if(pid_d==-211 ){ pi_rhominus=LV;}
632  if(abs(pid_d)==111 ){ pi0_rhominus=LV;}
633  }
634  if(pid==-15){
635  hasrho_plus=true;
636  if(pid_d==211 ){pi_rhoplus=LV;}
637  if(abs(pid_d)==111 ){pi0_rhoplus=LV;}
638  }
639  }
640  }
641  }
642  }
643  if(abs(findMother(*des)) != 15 && abs(pid) == 15 && (tauDecayChannel(*des) == pi) ){
644  if ( (*des)->end_vertex() ) {
645  HepMC::GenVertex::particle_iterator tauprod;
646  TLorentzVector LVtau((*des)->momentum().px(),(*des)->momentum().py(),(*des)->momentum().pz(),(*des)->momentum().e());
647  if(pid == 15)taum=LVtau;
648  if(pid ==-15)taup=LVtau;
649  for(tauprod = (*des)->end_vertex()->particles_begin(HepMC::descendants); tauprod!= (*des)->end_vertex()->particles_end(HepMC::descendants);++tauprod ) {
650  int pid_d = (*tauprod)->pdg_id();
651  if(abs(pid_d)==211 || abs(pid_d)==111){
652  TLorentzVector LV((*tauprod)->momentum().px(),(*tauprod)->momentum().py(),(*tauprod)->momentum().pz(),(*tauprod)->momentum().e());
653  if(pid==15){
654  haspi_minus=true;
655  if(pid_d==-211 ){ pi_minus=LV;}
656  }
657  if(pid==-15){
658  haspi_plus=true;
659  if(pid_d==211 ){pi_plus=LV;}
660  }
661  }
662  }
663  }
664  }
665  }
666  }
667  if(hasrho_minus && hasrho_plus){
668  //compute rhorho
669  rho_minus=pi_rhominus;
670  rho_minus+=pi0_rhominus;
671  rho_plus=pi_rhoplus;
672  rho_plus+=pi0_rhoplus;
673  TLorentzVector rhorho=rho_minus;rhorho+=rho_plus;
674 
675  // boost to rhorho cm
676  TLorentzVector pi_rhoplusb=pi_rhoplus; pi_rhoplusb.Boost(-1*rhorho.BoostVector());
677  TLorentzVector pi0_rhoplusb=pi0_rhoplus; pi0_rhoplusb.Boost(-1*rhorho.BoostVector());
678  TLorentzVector pi_rhominusb=pi_rhominus; pi_rhominusb.Boost(-1*rhorho.BoostVector());
679  TLorentzVector pi0_rhominusb=pi0_rhominus; pi0_rhominusb.Boost(-1*rhorho.BoostVector());
680 
681  // compute n+/-
682  TVector3 n_plus=pi_rhoplusb.Vect().Cross(pi0_rhoplusb.Vect());
683  TVector3 n_minus=pi_rhominusb.Vect().Cross(pi0_rhominusb.Vect());
684 
685  // compute the acoplanarity
686  double Acoplanarity=acos(n_plus.Dot(n_minus)/(n_plus.Mag()*n_minus.Mag()));
687  if(pi_rhominus.Vect().Dot(n_plus)>0){Acoplanarity*=-1;Acoplanarity+=2*TMath::Pi();}
688 
689  // now boost to tau frame
690  pi_rhoplus.Boost(-1*taup.BoostVector());
691  pi0_rhoplus.Boost(-1*taup.BoostVector());
692  pi_rhominus.Boost(-1*taum.BoostVector());
693  pi0_rhominus.Boost(-1*taum.BoostVector());
694 
695  // compute y1 and y2
696  double y1=(pi_rhoplus.E()-pi0_rhoplus.E())/(pi_rhoplus.E()+pi0_rhoplus.E());
697  double y2=(pi_rhominus.E()-pi0_rhominus.E())/(pi_rhominus.E()+pi0_rhominus.E());
698 
699  // fill histograms
700  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0 && y1*y2<0) TauSpinEffectsH_rhorhoAcoplanarityminus->Fill(Acoplanarity,weight);
701  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0 && y1*y2>0) TauSpinEffectsH_rhorhoAcoplanarityplus->Fill(Acoplanarity,weight);
702  }
703  if(haspi_minus && haspi_plus){
704  TLorentzVector tauporig=taup;
705  TLorentzVector taumorig=taum;
706 
707  // now boost to Higgs frame
708  pi_plus.Boost(-1*Zboson.BoostVector());
709  pi_minus.Boost(-1*Zboson.BoostVector());
710 
711  taup.Boost(-1*Zboson.BoostVector());
712  taum.Boost(-1*Zboson.BoostVector());
713 
714  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0){
715  TauSpinEffectsH_pipiAcollinearity->Fill(acos(pi_plus.Vect().Dot(pi_minus.Vect())/(pi_plus.P()*pi_minus.P())));
716  TauSpinEffectsH_pipiAcollinearityzoom->Fill(acos(pi_plus.Vect().Dot(pi_minus.Vect())/(pi_plus.P()*pi_minus.P())));
717  }
718 
719  double proj_m=taum.Vect().Dot(pi_minus.Vect())/(taum.P()*taup.P());
720  double proj_p=taup.Vect().Dot(pi_plus.Vect())/(taup.P()*taup.P());
721  TVector3 Tau_m=taum.Vect();
722  TVector3 Tau_p=taup.Vect();
723  Tau_m*=proj_m;
724  Tau_p*=proj_p;
725  TVector3 Pit_m=pi_minus.Vect()-Tau_m;
726  TVector3 Pit_p=pi_plus.Vect()-Tau_p;
727 
728  double Acoplanarity=acos(Pit_m.Dot(Pit_p)/(Pit_p.Mag()*Pit_m.Mag()));
729  TVector3 n=Pit_p.Cross(Pit_m);
730  if(n.Dot(Tau_m)/Tau_m.Mag()>0){Acoplanarity*=-1; Acoplanarity+=2*TMath::Pi();}
731  // fill histograms
732  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_pipiAcoplanarity->Fill(Acoplanarity,weight);
733  taup=tauporig;
734  taum=taumorig;
735 
736  }
737 
738 
739  if(nSingleMuonDecays==2){
740  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_muX->Fill(x1,weight);
741  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_muX->Fill(x1,weight);
742  }
743  if(nSingleElectronDecays==2){
744  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_eX->Fill(x1,weight);
745  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_eX->Fill(x1,weight);
746  }
747  if(nSinglePionDecays == 2 && tautau.M()!= 0) {
748  for(int i=0;i<zsbins;i++){
749  double zslow=((double)i)*(zsmax-zsmin)/((double)zsbins)+zsmin;
750  double zsup=((double)i+1)*(zsmax-zsmin)/((double)zsbins)+zsmin;
751  double aup=Zstoa(zsup), alow=Zstoa(zslow);
752  if(x2-x1>alow && x2-x1<aup){
753  double zs=(zsup+zslow)/2;
754  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_Zs->Fill(zs,weight);
755  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Zs->Fill(zs,weight);
756  break;
757  }
758  }
759  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_MVis->Fill(pipi.M()/tautau.M(),weight);
760  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_MVis->Fill(pipi.M()/tautau.M(),weight);
761 
762  if(x1!=0){
763  const std::vector<HepMC::GenParticle*> m=GetMothers(boson);
764  int q(0),qbar(0);
765  TLorentzVector Z(0,0,0,0);
766  for(unsigned int i=0;i<m.size();i++){
767  if(m.at(i)->pdg_id()==PdtPdgMini::d || m.at(i)->pdg_id()==PdtPdgMini::u ){q++;}
768  if(m.at(i)->pdg_id()==PdtPdgMini::anti_d || m.at(i)->pdg_id()==PdtPdgMini::anti_u ){qbar++;}
769  }
770  if(q==1 && qbar==1){// assume q has largest E (valence vs see quarks)
771  if(taum.Vect().Dot(Zboson.Vect())/(Zboson.P()*taum.P())>0){
772  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_Xf->Fill(x1,weight);
773  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xf->Fill(x1,weight);
774  }
775  else{
776  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_Xb->Fill(x1,weight);
777  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xb->Fill(x1,weight);
778  }
779  }
780  }
781  }
782 }
783 
784 double TauValidation::Zstoa(double zs){
785  double a=1-sqrt(fabs(1.0-2*fabs(zs)));
786  if(zs<0){
787  a*=-1.0;
788  }
789  return a;
790 }
791 
792 
793 double TauValidation::leadingPionMomentum(const HepMC::GenParticle* tau, double weight){
794  return leadingPionP4(tau).P();
795 }
796 
797 TLorentzVector TauValidation::leadingPionP4(const HepMC::GenParticle* tau){
798 
799  TLorentzVector p4(0,0,0,0);
800 
801  if ( tau->end_vertex() ) {
802  HepMC::GenVertex::particle_iterator des;
803  for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
804  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
805  int pid = (*des)->pdg_id();
806 
807  if(abs(pid) == 15) return leadingPionP4(*des);
808 
809  if(!(abs(pid)==211 || abs(pid)==13 || abs(pid)==11)) continue;
810 
811  if((*des)->momentum().rho() > p4.P()) {
812  p4 = TLorentzVector((*des)->momentum().px(),
813  (*des)->momentum().py(),
814  (*des)->momentum().pz(),
815  (*des)->momentum().e());
816  }
817  }
818  }
819 
820  return p4;
821 }
822 
823 TLorentzVector TauValidation::motherP4(const HepMC::GenParticle* tau){
824  const HepMC::GenParticle* m=GetMother(tau);
825  return TLorentzVector(m->momentum().px(),m->momentum().py(),m->momentum().pz(),m->momentum().e());
826 }
827 
828 double TauValidation::visibleTauEnergy(const HepMC::GenParticle* tau){
829  TLorentzVector p4(tau->momentum().px(),
830  tau->momentum().py(),
831  tau->momentum().pz(),
832  tau->momentum().e());
833 
834  if ( tau->end_vertex() ) {
835  HepMC::GenVertex::particle_iterator des;
836  for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
837  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
838  int pid = (*des)->pdg_id();
839 
840  if(abs(pid) == 15) return visibleTauEnergy(*des);
841 
842  if(abs(pid) == 12 || abs(pid) == 14 || abs(pid) == 16) {
843  p4 -= TLorentzVector((*des)->momentum().px(),
844  (*des)->momentum().py(),
845  (*des)->momentum().pz(),
846  (*des)->momentum().e());
847  }
848  }
849  }
850 
851  return p4.E();
852 }
853 
854 void TauValidation::photons(const HepMC::GenParticle* tau, double weight){
855  // Find First tau in chain
856  std::vector<const HepMC::GenParticle*> TauList;
857  findTauList(tau,TauList);
858 
859  // Get List of Gauge Boson to tau(s) FSR and Brem
860  bool passedW=false;
861  std::vector<const HepMC::GenParticle*> ListofFSR; ListofFSR.clear();
862  std::vector<const HepMC::GenParticle*> ListofBrem; ListofBrem.clear();
863  std::vector<const HepMC::GenParticle*> FSR_photos; FSR_photos.clear();
864  double BosonScale(1);
865  if(TauList.size()>0){
866  TauValidation::findFSRandBrem(TauList.at(0),passedW,ListofFSR,ListofBrem);
867  TauValidation::FindPhotosFSR(TauList.at(0),FSR_photos,BosonScale);
868 
869  // Add the Tau Brem. information
870  TauBremPhotonsN->Fill(ListofBrem.size(),weight);
871  double photonPtSum=0;
872  for(unsigned int i=0;i<ListofBrem.size();i++){
873  photonPtSum+=ListofBrem.at(i)->momentum().perp();
874  TauBremPhotonsPt->Fill(ListofBrem.at(i)->momentum().perp(),weight);
875  }
876  TauBremPhotonsPtSum->Fill(photonPtSum,weight);
877 
878  // Now add the Gauge Boson FSR information
879  if(BosonScale!=0){
880  TauFSRPhotonsN->Fill(ListofFSR.size(),weight);
881  photonPtSum=0;
882  for(unsigned int i=0;i<ListofFSR.size();i++){
883  photonPtSum+=ListofFSR.at(i)->momentum().perp();
884  TauFSRPhotonsPt->Fill(ListofFSR.at(i)->momentum().perp(),weight);
885  }
886  double FSR_photosSum(0);
887  for(unsigned int i=0;i<FSR_photos.size();i++){
888  FSR_photosSum+=FSR_photos.at(i)->momentum().perp();
889  TauFSRPhotonsPt->Fill(FSR_photos.at(i)->momentum().perp()/BosonScale,weight*BosonScale);
890  }
891  TauFSRPhotonsPtSum->Fill(photonPtSum+FSR_photosSum/BosonScale,weight);
892  }
893  }
894 }
895 
TauValidation::findTauList
void findTauList(const HepMC::GenParticle *tau, std::vector< const HepMC::GenParticle * > &TauList)
Definition: TauValidation.cc:342
TauValidation::TauFSRPhotonsN
MonitorElement * TauFSRPhotonsN
Definition: TauValidation.h:109
beta
const double beta
Definition: CosmicMuonParameters.h:69
Pi
const double Pi
Definition: CosmicMuonParameters.h:18
i
int i
Definition: DBlmapReader.cc:9
TauValidation::~TauValidation
virtual ~TauValidation()
Definition: TauValidation.cc:33
TauValidation::n_taupinu
std::vector< int > n_taupinu
Definition: TauValidation.h:127
TauValidation::TauSpinEffectsH_rhorhoAcoplanarityplus
MonitorElement * TauSpinEffectsH_rhorhoAcoplanarityplus
Definition: TauValidation.h:109
TauValidation::TauPhi
MonitorElement * TauPhi
Definition: TauValidation.h:109
TauValidation::TauSpinEffectsW_X
MonitorElement * TauSpinEffectsW_X
Definition: TauValidation.h:109
parents
TPRegexp parents
Definition: eve_filter.cc:24
TauValidation::TauSpinEffectsHpm_UpsilonA1
MonitorElement * TauSpinEffectsHpm_UpsilonA1
Definition: TauValidation.h:109
TauValidation::JAKInvMass
std::vector< std::vector< MonitorElement * > > JAKInvMass
Definition: TauValidation.h:123
TauValidation::TauSpinEffectsZ_Xb
MonitorElement * TauSpinEffectsZ_Xb
Definition: TauValidation.h:109
TauValidation::nTaus
MonitorElement * nTaus
Definition: TauValidation.h:108
DQMStore::book1D
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:717
TauValidation::DecayLength
MonitorElement * DecayLength
Definition: TauValidation.h:109
TauValidation::TauRtauW
MonitorElement * TauRtauW
Definition: TauValidation.h:109
TauValidation::TauValidation
TauValidation(const edm::ParameterSet &)
Definition: TauValidation.cc:19
TauValidation::TauPt
MonitorElement * TauPt
Definition: TauValidation.h:109
TauValidation::TauSpinEffectsH_MVis
MonitorElement * TauSpinEffectsH_MVis
Definition: TauValidation.h:109
TauValidation::spinEffectsZH
void spinEffectsZH(const HepMC::GenParticle *boson, double weight)
Definition: TauValidation.cc:564
TauValidation::TauMothers
MonitorElement * TauMothers
Definition: TauValidation.h:109
TauValidation::TauSpinEffectsHpm_eX
MonitorElement * TauSpinEffectsHpm_eX
Definition: TauValidation.h:109
TauValidation::TauSpinEffectsW_muX
MonitorElement * TauSpinEffectsW_muX
Definition: TauValidation.h:109
TauValidation::TauBremPhotonsN
MonitorElement * TauBremPhotonsN
Definition: TauValidation.h:109
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:889
abs
#define abs(x)
Definition: mlp_lapack.h:159
TauValidation::JAKID
MonitorElement * JAKID
Definition: TauValidation.h:122
TauValidation::TauSpinEffectsZ_Xf
MonitorElement * TauSpinEffectsZ_Xf
Definition: TauValidation.h:109
TauValidation::tauProngs
int tauProngs(const HepMC::GenParticle *, double weight)
Definition: TauValidation.cc:426
TauDecay_CMSSW
Definition: TauDecay_CMSSW.h:31
TauValidation::TauBremPhotonsPtSum
MonitorElement * TauBremPhotonsPtSum
Definition: TauValidation.h:109
TauValidation::endJob
virtual void endJob()
Definition: TauValidation.cc:165
TauValidation::TauSpinEffectsH_pipiAcoplanarity
MonitorElement * TauSpinEffectsH_pipiAcoplanarity
Definition: TauValidation.h:109
TauValidation::TauSpinEffectstautau_polvxM
MonitorElement * TauSpinEffectstautau_polvxM
Definition: TauValidation.h:109
TauDecay_CMSSW.h
TauValidation::TauSpinEffectsHpm_UpsilonRho
MonitorElement * TauSpinEffectsHpm_UpsilonRho
Definition: TauValidation.h:109
DeDxDiscriminatorTools::charge
double charge(const std::vector< uint8_t > &Ampls)
Definition: DeDxDiscriminatorTools.cc:43
TauValidation::analyze
virtual void analyze(const edm::Event &, const edm::EventSetup &)
Definition: TauValidation.cc:176
edm::EventSetup::getData
void getData(T &iHolder) const
Definition: EventSetup.h:67
TauValidation::TauSpinEffectsZ_eX
MonitorElement * TauSpinEffectsZ_eX
Definition: TauValidation.h:109
TauValidation::findMother
int findMother(const HepMC::GenParticle *)
Definition: TauValidation.cc:317
MonitorElement::Fill
void Fill(long long x)
Definition: MonitorElement.h:152
TauValidation::TauProngs
MonitorElement * TauProngs
Definition: TauValidation.h:109
TauValidation::tauMother
int tauMother(const HepMC::GenParticle *, double weight)
Definition: TauValidation.cc:401
double_binary
Definition: DDStreamer.cc:227
TauValidation::tauDecayChannel
int tauDecayChannel(const HepMC::GenParticle *, double weight=0.0)
Definition: TauValidation.cc:446
TauDecay::isTauFinalStateParticle
bool isTauFinalStateParticle(int pdgid)
Definition: TauDecay.cc:38
TauValidation::TauSpinEffectsH_pipiAcollinearityzoom
MonitorElement * TauSpinEffectsH_pipiAcollinearityzoom
Definition: TauValidation.h:109
TauValidation::GetMothers
const std::vector< HepMC::GenParticle * > GetMothers(const HepMC::GenParticle *boson)
Definition: TauValidation.cc:304
LV
math::XYZTLorentzVectorD LV
Definition: HLTTauDQMPlotter.h:22
iEvent
int iEvent
Definition: GenABIO.cc:243
TauValidation::TauSpinEffectsW_UpsilonRho
MonitorElement * TauSpinEffectsW_UpsilonRho
Definition: TauValidation.h:109
TauValidation::photons
void photons(const HepMC::GenParticle *, double weight)
Definition: TauValidation.cc:854
TauValidation::TauSpinEffectsH_rhorhoAcoplanarityminus
MonitorElement * TauSpinEffectsH_rhorhoAcoplanarityminus
Definition: TauValidation.h:109
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:46
p4
double p4[4]
Definition: TauolaWrapper.h:92
TauValidation::leadingPionMomentum
double leadingPionMomentum(const HepMC::GenParticle *, double weight)
Definition: TauValidation.cc:793
TauValidation::beginJob
virtual void beginJob()
Definition: TauValidation.cc:35
TauValidation::TauSpinEffectsH_Xf
MonitorElement * TauSpinEffectsH_Xf
Definition: TauValidation.h:109
TauValidation::TauFSRPhotonsPt
MonitorElement * TauFSRPhotonsPt
Definition: TauValidation.h:109
MonitorElement::getTH1
TH1 * getTH1(void) const
Definition: MonitorElement.cc:1479
TauValidation::TauSpinEffectsZ_muX
MonitorElement * TauSpinEffectsZ_muX
Definition: TauValidation.h:109
TauValidation::TauDecayChannels
MonitorElement * TauDecayChannels
Definition: TauValidation.h:109
TauValidation::hepmcCollection_
edm::InputTag hepmcCollection_
Definition: TauValidation.h:98
TauValidation::TauSpinEffectsH_Zs
MonitorElement * TauSpinEffectsH_Zs
Definition: TauValidation.h:109
TauValidation::TauFSRPhotonsPtSum
MonitorElement * TauFSRPhotonsPtSum
Definition: TauValidation.h:109
ParticleData
HepPDT::ParticleData ParticleData
Definition: ParticleDataTable.h:10
TauValidation::visibleTauEnergy
double visibleTauEnergy(const HepMC::GenParticle *)
Definition: TauValidation.cc:828
edm::Event::getByLabel
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:356
TauValidation::leadingPionP4
TLorentzVector leadingPionP4(const HepMC::GenParticle *)
Definition: TauValidation.cc:797
TauValidation::isLastTauinChain
bool isLastTauinChain(const HepMC::GenParticle *tau)
Definition: TauValidation.cc:330
TauValidation::TauSpinEffectsZ_MVis
MonitorElement * TauSpinEffectsZ_MVis
Definition: TauValidation.h:109
TauValidation::TauSpinEffectsH_pipiAcollinearity
MonitorElement * TauSpinEffectsH_pipiAcollinearity
Definition: TauValidation.h:109
TauValidation::motherP4
TLorentzVector motherP4(const HepMC::GenParticle *)
Definition: TauValidation.cc:823
TauValidation::rtau
void rtau(const HepMC::GenParticle *, int, int, double weight)
Definition: TauValidation.cc:499
TauValidation::_wmanager
WeightManager _wmanager
Definition: TauValidation.h:96
part
part
Definition: HCALResponse.h:21
TauValidation::FindPhotosFSR
void FindPhotosFSR(const HepMC::GenParticle *p, std::vector< const HepMC::GenParticle * > &ListofFSR, double &BosonScale)
Definition: TauValidation.cc:378
TauValidation::TauSpinEffectsZ_X
MonitorElement * TauSpinEffectsZ_X
Definition: TauValidation.h:109
TauValidation::endRun
virtual void endRun(const edm::Run &, const edm::EventSetup &)
Definition: TauValidation.cc:175
TauValidation::fPDGTable
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
PDT table.
Definition: TauValidation.h:103
AlCaHLTBitMon_ParallelJobs.p
tuple p
Definition: AlCaHLTBitMon_ParallelJobs.py:152
TauValidation::TauSpinEffectsH_eX
MonitorElement * TauSpinEffectsH_eX
Definition: TauValidation.h:109
TauValidation::TauSpinEffectsW_eX
MonitorElement * TauSpinEffectsW_eX
Definition: TauValidation.h:109
tmp
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
TauValidation::Zstoa
double Zstoa(double zs)
Definition: TauValidation.cc:784
TauValidation::GetMother
const HepMC::GenParticle * GetMother(const HepMC::GenParticle *tau)
Definition: TauValidation.cc:292
TauValidation::spinEffectsWHpm
void spinEffectsWHpm(const HepMC::GenParticle *, int, int, std::vector< HepMC::GenParticle * > &part, double weight)
Definition: TauValidation.cc:515
TauValidation::LifeTime
MonitorElement * LifeTime
Definition: TauValidation.h:109
TauValidation::TauSpinEffectsH_muX
MonitorElement * TauSpinEffectsH_muX
Definition: TauValidation.h:109
TauValidation::TauSpinEffectsW_UpsilonA1
MonitorElement * TauSpinEffectsW_UpsilonA1
Definition: TauValidation.h:109
a
double a
Definition: hdecay.h:121
TauValidation::TauEta
MonitorElement * TauEta
Definition: TauValidation.h:109
TauValidation::dbe
DQMStore * dbe
ME&#39;s &quot;container&quot;.
Definition: TauValidation.h:106
TauValidation::findFSRandBrem
void findFSRandBrem(const HepMC::GenParticle *p, bool doBrem, std::vector< const HepMC::GenParticle * > &ListofFSR, std::vector< const HepMC::GenParticle * > &ListofBrem)
Definition: TauValidation.cc:354
relval_parameters_module.energy
string energy
Definition: relval_parameters_module.py:29
TauValidation::NJAKID
unsigned int NJAKID
Definition: TauValidation.h:121
TauValidation::TauSpinEffectsZ_Zs
MonitorElement * TauSpinEffectsZ_Zs
Definition: TauValidation.h:109
TauValidation::TauSpinEffectsHpm_X
MonitorElement * TauSpinEffectsHpm_X
Definition: TauValidation.h:109
TauValidation::beginRun
virtual void beginRun(const edm::Run &, const edm::EventSetup &)
Definition: TauValidation.cc:169
vdt::x
x
Definition: VDTMath.h:216
TauValidation::TauSpinEffectsHpm_muX
MonitorElement * TauSpinEffectsHpm_muX
Definition: TauValidation.h:109
WeightManager::weight
double weight(const edm::Event &)
Definition: WeightManager.cc:21
MonitorElement::setAxisTitle
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
Definition: MonitorElement.cc:916
TauValidation::TauSpinEffectsH_Xb
MonitorElement * TauSpinEffectsH_Xb
Definition: TauValidation.h:109
TauDecay_CMSSW::Get_TauDecayProducts
std::vector< HepMC::GenParticle * > Get_TauDecayProducts()
Definition: TauDecay_CMSSW.h:39
TauDecay_CMSSW::AnalyzeTau
bool AnalyzeTau(HepMC::GenParticle *Tau, unsigned int &JAK_ID, unsigned int &TauBitMask, bool dores=true, bool dopi0=true)
Definition: TauDecay_CMSSW.cc:17
TauValidation::nPrimeTaus
MonitorElement * nPrimeTaus
Definition: TauValidation.h:108
TauValidation::TauBremPhotonsPt
MonitorElement * TauBremPhotonsPt
Definition: TauValidation.h:109
TauValidation::TauRtauHpm
MonitorElement * TauRtauHpm
Definition: TauValidation.h:109
DQMStore::setCurrentFolder
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:429
edm::Run
Definition: Run.h:33
TauDecay::nProng
unsigned int nProng(unsigned int &TauBitMask)
Definition: TauDecay.h:104
TauValidation::TauSpinEffectsH_X
MonitorElement * TauSpinEffectsH_X
Definition: TauValidation.h:109
Generated for CMSSW Reference Manual by doxygen 1.8.5