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 
8 
9 #include "CLHEP/Units/defs.h"
10 #include "CLHEP/Units/PhysicalConstants.h"
11 
13 
17 
18 using namespace edm;
19 
21  _wmanager(iPSet)
22  ,hepmcCollection_(iPSet.getParameter<edm::InputTag>("hepmcCollection"))
23  ,tauEtCut(iPSet.getParameter<double>("tauEtCutForRtau"))
24  ,NJAKID(22)
25  ,zsbins(20)
26  ,zsmin(-0.5)
27  ,zsmax(0.5)
28 {
29  dbe = 0;
31 }
32 
34 
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);
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}");
63  TauDecayChannels->setBinLabel(1+stable,"Stable");
64 
65  TauMothers = dbe->book1D("TauMothers","Tau mother particles", 10 ,0,10);
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  TauRtauW = dbe->book1D("TauRtauW","W->Tau p(leading track)/E(visible tau)", 50 ,0,1); TauRtauW->setAxisTitle("rtau");
78  TauRtauHpm = dbe->book1D("TauRtauHpm","Hpm->Tau p(leading track)/E(visible tau)", 50 ,0,1); TauRtauHpm->setAxisTitle("rtau");
79 
80  TauSpinEffectsW_X = dbe->book1D("TauSpinEffectsWX","Pion energy in W rest frame", 50 ,0,1); TauSpinEffectsW_X->setAxisTitle("X");
81  TauSpinEffectsHpm_X = dbe->book1D("TauSpinEffectsHpmX","Pion energy in Hpm rest frame", 50 ,0,1); TauSpinEffectsHpm_X->setAxisTitle("X");
82 
83  TauSpinEffectsW_eX = dbe->book1D("TauSpinEffectsWeX","e energy in W rest frame", 50 ,0,1); TauSpinEffectsW_eX->setAxisTitle("X");
84  TauSpinEffectsHpm_eX = dbe->book1D("TauSpinEffectsHpmeX","e energy in Hpm rest frame", 50 ,0,1); TauSpinEffectsHpm_eX->setAxisTitle("X");
85 
86  TauSpinEffectsW_muX = dbe->book1D("TauSpinEffectsWmuX","mu energy in W rest frame", 50 ,0,1); TauSpinEffectsW_muX->setAxisTitle("X");
87  TauSpinEffectsHpm_muX = dbe->book1D("TauSpinEffectsHpmmuX","mu energy in Hpm rest frame", 50 ,0,1); TauSpinEffectsHpm_muX->setAxisTitle("X");
88 
89  TauSpinEffectsW_UpsilonRho = dbe->book1D("TauSpinEffectsWUpsilonRho","#Upsilon for #rho", 50 ,-1,1); TauSpinEffectsW_UpsilonRho->setAxisTitle("#Upsilon");
90  TauSpinEffectsHpm_UpsilonRho = dbe->book1D("TauSpinEffectsHpmUpsilonRho","#Upsilon for #rho", 50 ,-1,1); TauSpinEffectsHpm_UpsilonRho->setAxisTitle("#Upsilon");
91 
92  TauSpinEffectsW_UpsilonA1 = dbe->book1D("TauSpinEffectsWUpsilonA1","#Upsilon for a1", 50 ,-1,1); TauSpinEffectsW_UpsilonA1->setAxisTitle("#Upsilon");
93  TauSpinEffectsHpm_UpsilonA1 = dbe->book1D("TauSpinEffectsHpmUpsilonA1","#Upsilon for a1", 50 ,-1,1); TauSpinEffectsHpm_UpsilonA1->setAxisTitle("#Upsilon");
94 
95  TauSpinEffectsZ_MVis = dbe->book1D("TauSpinEffectsZMVis","Mass of pi+ pi-", 25 ,0,1.1); TauSpinEffectsZ_MVis->setAxisTitle("M_{#pi^{+}#pi^{-}}");
96  TauSpinEffectsH_MVis = dbe->book1D("TauSpinEffectsHMVis","Mass of pi+ pi-", 25 ,0,1.1); TauSpinEffectsZ_MVis->setAxisTitle("M_{#pi^{+}#pi^{-}}");
97 
98  TauSpinEffectsZ_Zs = dbe->book1D("TauSpinEffectsZZs","Z_{s}", zsbins ,zsmin,zsmax); TauSpinEffectsZ_Zs->setAxisTitle("Z_{s}");
99  TauSpinEffectsH_Zs = dbe->book1D("TauSpinEffectsHZs","Z_{s}", zsbins ,zsmin,zsmax); TauSpinEffectsZ_Zs->setAxisTitle("Z_{s}");
100 
101  TauSpinEffectsZ_Xf = dbe->book1D("TauSpinEffectsZXf","X of forward emitted #tau^{-}", 25 ,0,1.0); TauSpinEffectsZ_Xf->setAxisTitle("X_{f}");
102  TauSpinEffectsH_Xf = dbe->book1D("TauSpinEffectsHXf","X of forward emitted #tau^{-}", 25 ,0,1.0); TauSpinEffectsZ_Xf->setAxisTitle("X_{f}");
103 
104  TauSpinEffectsZ_Xb = dbe->book1D("TauSpinEffectsZXb","X of backward emitted #tau^{-}", 25 ,0,1.0); TauSpinEffectsZ_Xb->setAxisTitle("X_{b}");
105  TauSpinEffectsH_Xb = dbe->book1D("TauSpinEffectsHXb","X of backward emitted #tau^{-}", 25 ,0,1.0); TauSpinEffectsZ_Xb->setAxisTitle("X_{b}");
106 
107  TauSpinEffectsZ_eX = dbe->book1D("TauSpinEffectsZeX","e energy in Z rest frame", 50 ,0,1); TauSpinEffectsZ_eX->setAxisTitle("X");
108  TauSpinEffectsH_eX = dbe->book1D("TauSpinEffectsHeX","e energy in H rest frame", 50 ,0,1); TauSpinEffectsH_eX->setAxisTitle("X");
109 
110  TauSpinEffectsZ_muX = dbe->book1D("TauSpinEffectsZmuX","mu energy in z rest frame", 50 ,0,1); TauSpinEffectsZ_muX->setAxisTitle("X");
111  TauSpinEffectsH_muX = dbe->book1D("TauSpinEffectsHmuX","mu energy in H rest frame", 50 ,0,1); TauSpinEffectsH_muX->setAxisTitle("X");
112 
113 
114  TauFSRPhotonsN=dbe->book1D("TauFSRPhotonsN","FSR Photons radiating from/with tau (Gauge Boson)", 5 ,-0.5,4.5);
115  TauFSRPhotonsN->setAxisTitle("N FSR Photons radiating from/with tau");
116  TauFSRPhotonsPt=dbe->book1D("TauFSRPhotonsPt","Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
117  TauFSRPhotonsPt->setAxisTitle("P_{t} of FSR Photons radiating from/with tau [per tau]");
118  TauFSRPhotonsPtSum=dbe->book1D("TauFSRPhotonsPtSum","Pt of FSR Photons radiating from/with tau (Gauge Boson)", 100 ,0,100);
119  TauFSRPhotonsPtSum->setAxisTitle("P_{t} of FSR Photons radiating from/with tau [per tau]");
120 
121  TauBremPhotonsN=dbe->book1D("TauBremPhotonsN","Brem. Photons radiating in tau decay", 5 ,-0.5,4.5);
122  TauBremPhotonsN->setAxisTitle("N FSR Photons radiating from/with tau");
123  TauBremPhotonsPt=dbe->book1D("TauBremPhotonsPt","Sum Brem Pt ", 100 ,0,100);
124  TauFSRPhotonsPt->setAxisTitle("P_{t} of Brem. Photons radiating in tau decay");
125  TauBremPhotonsPtSum =dbe->book1D("TauBremPhotonsPtSum","Sum of Brem Pt ", 100 ,0,100);
126  TauFSRPhotonsPtSum->setAxisTitle("Sum P_{t} of Brem. Photons radiating in tau decay");
127 
128  JAKID =dbe->book1D("JAKID","JAK ID",NJAKID+1,-0.5,NJAKID+0.5);
129  for(unsigned int i=0; i<NJAKID+1;i++){
130  JAKInvMass.push_back(std::vector<MonitorElement *>());
131  TString tmp="JAKID";
132  tmp+=i;
133  JAKInvMass.at(i).push_back(dbe->book1D("M"+tmp,"M_{"+tmp+"} (GeV)", 80 ,0,2.0));
134  if(i==TauDecay::JAK_A1_3PI ||
137  JAKInvMass.at(i).push_back(dbe->book1D("M13"+tmp,"M_{13,"+tmp+"} (GeV)", 80 ,0,2.0));
138  JAKInvMass.at(i).push_back(dbe->book1D("M23"+tmp,"M_{23,"+tmp+"} (GeV)", 80 ,0,2.0));
139  JAKInvMass.at(i).push_back(dbe->book1D("M12"+tmp,"M_{12,"+tmp+"} (GeV)", 80 ,0,2.0));
140  }
141  }
142  }
143 
144  return;
145 }
146 
148  return;
149 }
150 
151 void TauValidation::beginRun(const edm::Run& iRun,const edm::EventSetup& iSetup)
152 {
154  iSetup.getData( fPDGTable );
155  return;
156 }
157 void TauValidation::endRun(const edm::Run& iRun,const edm::EventSetup& iSetup){return;}
159 {
162  iEvent.getByLabel(hepmcCollection_, evt);
163 
164  //Get EVENT
165  HepMC::GenEvent *myGenEvent = new HepMC::GenEvent(*(evt->GetEvent()));
166 
167  double weight = _wmanager.weight(iEvent);
168 
170  /*
171  edm::Handle<double> WT;
172  iEvent.getByLabel(edm::InputTag("TauSpinnerGen","TauSpinnerWT"),WT);
173  weight = 1.0;
174  if(*(WT.product())>1e-3 && *(WT.product())<=10.0) weight=(*(WT.product()));
175  else {weight=1.0;}
176  */
178 
179  // find taus
180  for(HepMC::GenEvent::particle_const_iterator iter = myGenEvent->particles_begin(); iter != myGenEvent->particles_end(); iter++) {
181  if(abs((*iter)->pdg_id())==23){
182  spinEffectsZ(*iter,weight);
183  }
184  if(abs((*iter)->pdg_id())==15){
185  if(isLastTauinChain(*iter)){
186  nTaus->Fill(0.5,weight);
187  int mother = tauMother(*iter,weight);
188  int decaychannel = tauDecayChannel(*iter,weight);
189  tauProngs(*iter, weight);
190  if(mother>-1){ // exclude B, D and other non-signal decay modes
191  nPrimeTaus->Fill(0.5,weight);
192  TauPt->Fill((*iter)->momentum().perp(),weight);
193  TauEta->Fill((*iter)->momentum().eta(),weight);
194  TauPhi->Fill((*iter)->momentum().phi(),weight);
195  rtau(*iter,mother,decaychannel,weight);
196  photons(*iter,weight);
197  }
199  //Adding JAKID and Mass information
200  //
201  TauDecay_CMSSW TD;
202  unsigned int jak_id, TauBitMask;
203  if(TD.AnalyzeTau((*iter),jak_id,TauBitMask,false,false)){
204  JAKID->Fill(jak_id,weight);
205  if(jak_id<=NJAKID){
206  int tcharge=(*iter)->pdg_id()/abs((*iter)->pdg_id());
207  std::vector<HepMC::GenParticle*> part=TD.Get_TauDecayProducts();
208  spinEffects(*iter,mother,jak_id,part,weight);
209  TLorentzVector LVQ(0,0,0,0);
210  TLorentzVector LVS12(0,0,0,0);
211  TLorentzVector LVS13(0,0,0,0);
212  TLorentzVector LVS23(0,0,0,0);
213  bool haspart1=false;
214  for(unsigned int i=0;i<part.size();i++){
215  if(TD.isTauFinalStateParticle(part.at(i)->pdg_id()) &&
216  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_e &&
217  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_mu &&
218  abs(part.at(i)->pdg_id())!=PdtPdgMini::nu_tau ){
219  TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
220  LVQ+=LV;
221  if(jak_id==TauDecay::JAK_A1_3PI ||
222  jak_id==TauDecay::JAK_KPIK ||
223  jak_id==TauDecay::JAK_KPIPI
224  ){
225  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) ){
226  LVS13+=LV;
227  LVS23+=LV;
228  }
229  else{
230  LVS12+=LV;
231  if(!haspart1 && ((jak_id==TauDecay::JAK_A1_3PI) || (jak_id!=TauDecay::JAK_A1_3PI && abs(part.at(i)->pdg_id())==PdtPdgMini::K_plus) )){
232  LVS13+=LV;
233  haspart1=true;
234  }
235  else{
236  LVS23+=LV;
237  }
238  }
239  }
240  }
241  }
242  part.clear();
243  JAKInvMass.at(jak_id).at(0)->Fill(LVQ.M(),weight);
244  if(jak_id==TauDecay::JAK_A1_3PI ||
245  jak_id==TauDecay::JAK_KPIK ||
246  jak_id==TauDecay::JAK_KPIPI
247  ){
248  JAKInvMass.at(jak_id).at(1)->Fill(LVS13.M(),weight);
249  JAKInvMass.at(jak_id).at(2)->Fill(LVS23.M(),weight);
250  JAKInvMass.at(jak_id).at(3)->Fill(LVS12.M(),weight);
251  }
252  }
253  }
254  else{
255  JAKID->Fill(jak_id,weight);
256  }
257  }
258  }
259  }
260  delete myGenEvent;
261 }//analyze
262 
264  if ( tau->production_vertex() ) {
265  HepMC::GenVertex::particle_iterator mother;
266  for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
267  if((*mother)->pdg_id() == tau->pdg_id()) return GetMother(*mother);
268  return (*mother);
269  }
270  }
271  return tau;
272 }
273 
274 
275 const std::vector<HepMC::GenParticle*> TauValidation::GetMothers(const HepMC::GenParticle* boson){
276  std::vector<HepMC::GenParticle*> mothers;
277  if ( boson->production_vertex() ) {
278  HepMC::GenVertex::particle_iterator mother;
279  for (mother = boson->production_vertex()->particles_begin(HepMC::parents); mother!= boson->production_vertex()->particles_end(HepMC::parents); mother++ ) {
280  if((*mother)->pdg_id() == boson->pdg_id()) return GetMothers(*mother);
281  mothers.push_back(*mother);
282  }
283  }
284  return mothers;
285 }
286 
287 
289  int mother_pid = 0;
290  if ( tau->production_vertex() ) {
291  HepMC::GenVertex::particle_iterator mother;
292  for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents); mother++ ) {
293  mother_pid = (*mother)->pdg_id();
294  if(mother_pid == tau->pdg_id()) return findMother(*mother); //mother_pid = -1; Make recursive to look for last tau in chain
295  }
296  }
297  return mother_pid;
298 }
299 
300 
302  if ( tau->end_vertex() ) {
303  HepMC::GenVertex::particle_iterator dau;
304  for (dau = tau->end_vertex()->particles_begin(HepMC::children); dau!= tau->end_vertex()->particles_end(HepMC::children); dau++ ) {
305  int dau_pid = (*dau)->pdg_id();
306  if(dau_pid == tau->pdg_id()) return false;
307  }
308  }
309  return true;
310 }
311 
312 
313 void TauValidation::findTauList(const HepMC::GenParticle* tau,std::vector<const HepMC::GenParticle*> &TauList){
314  TauList.insert(TauList.begin(),tau);
315  if ( tau->production_vertex() ) {
316  HepMC::GenVertex::particle_iterator mother;
317  for (mother = tau->production_vertex()->particles_begin(HepMC::parents); mother!= tau->production_vertex()->particles_end(HepMC::parents);mother++) {
318  if((*mother)->pdg_id() == tau->pdg_id()){
319  findTauList(*mother,TauList);
320  }
321  }
322  }
323 }
324 
325 void TauValidation::findFSRandBrem(const HepMC::GenParticle* p, bool doBrem, std::vector<const HepMC::GenParticle*> &ListofFSR,
326  std::vector<const HepMC::GenParticle*> &ListofBrem){
327  // 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.
328  if(abs(p->pdg_id())==15){
329  if(isLastTauinChain(p)){ doBrem=true;}
330  else{ doBrem=false;}
331  }
332  int photo_ID=22;
333  if ( p->end_vertex() ) {
334  HepMC::GenVertex::particle_iterator dau;
335  for (dau = p->end_vertex()->particles_begin(HepMC::children); dau!= p->end_vertex()->particles_end(HepMC::children); dau++ ) {
336  //if(doBrem) std::cout << "true " << (*dau)->pdg_id() << " " << std::endl;
337  //else std::cout << "false " << (*dau)->pdg_id() << " " << std::endl;
338  if(abs((*dau)->pdg_id()) == abs(photo_ID) && !doBrem){ListofFSR.push_back(*dau);}
339  if(abs((*dau)->pdg_id()) == abs(photo_ID) && doBrem){ListofBrem.push_back(*dau);}
340  if((*dau)->end_vertex() && (*dau)->end_vertex()->particles_out_size()>0 && abs((*dau)->pdg_id()) != 111 && abs((*dau)->pdg_id()) != 221/* remove pi0 and eta decays*/){
341  findFSRandBrem(*dau,doBrem,ListofFSR,ListofBrem);
342  }
343  }
344  }
345 }
346 
347 
348 
349 void TauValidation::FindPhotosFSR(const HepMC::GenParticle* p,std::vector<const HepMC::GenParticle*> &ListofFSR,double &BosonScale){
350  BosonScale=0.0;
351  const HepMC::GenParticle* m=GetMother(p);
352  double mother_pid=m->pdg_id();
353  if(m->end_vertex() && mother_pid!=p->pdg_id()){
354  HepMC::GenVertex::particle_iterator dau;
355  for (dau = m->end_vertex()->particles_begin(HepMC::children); dau!= m->end_vertex()->particles_end(HepMC::children); dau++ ) {
356  int dau_pid = (*dau)->pdg_id();
357  if(fabs(dau_pid) == 22) {
358  ListofFSR.push_back(*dau);
359  }
360  }
361  }
362  if(abs(mother_pid) == 24) BosonScale=1.0; // W
363  if(abs(mother_pid) == 23) BosonScale=2.0; // Z;
364  if(abs(mother_pid) == 22) BosonScale=2.0; // gamma;
365  if(abs(mother_pid) == 25) BosonScale=2.0; // HSM;
366  if(abs(mother_pid) == 35) BosonScale=2.0; // H0;
367  if(abs(mother_pid) == 36) BosonScale=2.0; // A0;
368  if(abs(mother_pid) == 37) BosonScale=1.0; //Hpm;
369 }
370 
371 
373 
374  if(abs(tau->pdg_id()) != 15 ) return -3;
375 
376  int mother_pid = findMother(tau);
377  if(mother_pid == -2) return -2;
378 
379  int label = other;
380  if(abs(mother_pid) == 24) label = W;
381  if(abs(mother_pid) == 23) label = Z;
382  if(abs(mother_pid) == 22) label = gamma;
383  if(abs(mother_pid) == 25) label = HSM;
384  if(abs(mother_pid) == 35) label = H0;
385  if(abs(mother_pid) == 36) label = A0;
386  if(abs(mother_pid) == 37) label = Hpm;
387 
388  int mother_shortpid=(abs(mother_pid)%10000);
389  if(mother_shortpid>500 && mother_shortpid<600 )label = B;
390  if(mother_shortpid>400 && mother_shortpid<500)label = D;
391  TauMothers->Fill(label,weight);
392  if(label==B || label == D || label == other) return -1;
393 
394  return mother_pid;
395 }
396 
398  int nProngs = 0;
399  if ( tau->end_vertex() ) {
400  HepMC::GenVertex::particle_iterator des;
401  for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
402  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
403  int pid = (*des)->pdg_id();
404  if(abs(pid) == 15) return tauProngs(*des, weight);
405  if((*des)->status() != 1) continue; // dont count unstable particles
406 
407  const HepPDT::ParticleData* pd = fPDGTable->particle((*des)->pdg_id ());
408  int charge = (int) pd->charge();
409  if(charge == 0) continue;
410  nProngs++;
411  }
412  }
413  TauProngs->Fill(nProngs,weight);
414  return nProngs;
415 }
416 
418 
419  int channel = undetermined;
420  if(tau->status() == 1) channel = stable;
421  int allCount = 0,
422  eCount = 0,
423  muCount = 0,
424  pi0Count = 0,
425  piCount = 0,
426  rhoCount = 0,
427  a1Count = 0,
428  KCount = 0,
429  KstarCount = 0;
430 
431  if ( tau->end_vertex() ) {
432  HepMC::GenVertex::particle_iterator des;
433  for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
434  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
435  int pid = (*des)->pdg_id();
436  if(abs(pid) == 15) return tauDecayChannel(*des,weight);
437 
438  allCount++;
439  if(abs(pid) == 11) eCount++;
440  if(abs(pid) == 13) muCount++;
441  if(abs(pid) == 111) pi0Count++;
442  if(abs(pid) == 211) piCount++;
443  if(abs(pid) == 213) rhoCount++;
444  if(abs(pid) == 20213) a1Count++;
445  if(abs(pid) == 321) KCount++;
446  if(abs(pid) == 323) KstarCount++;
447 
448  }
449  }
450  // resonances
451  if(KCount >= 1) channel = K;
452  if(KstarCount >= 1) channel = Kstar;
453  if(a1Count >= 1) channel = a1;
454  if(rhoCount >= 1) channel = rho;
455  if(channel!=undetermined && weight!=0.0) TauDecayChannels->Fill(channel,weight);
456 
457  // final state products
458  if(piCount == 1 && pi0Count == 0) channel = pi;
459  if(piCount == 1 && pi0Count == 1) channel = pi1pi0;
460  if(piCount == 1 && pi0Count > 1) channel = pinpi0;
461  if(piCount == 3 && pi0Count == 0) channel = tripi;
462  if(piCount == 3 && pi0Count > 0) channel = tripinpi0;
463  if(eCount == 1) channel = electron;
464  if(muCount == 1) channel = muon;
465  if(weight!=0.0) TauDecayChannels->Fill(channel,weight);
466  return channel;
467 }
468 
469 
470 void TauValidation::rtau(const HepMC::GenParticle* tau,int mother, int decay, double weight){
471 
472  if(decay != pi1pi0) return; // polarization only for 1-prong hadronic taus with one neutral pion to make a clean case
473 
474  if(tau->momentum().perp() < tauEtCut) return; // rtau visible only for boosted taus
475 
476  double rTau = 0;
477  double ltrack = leadingPionMomentum(tau, weight);
478  double visibleTauE = visibleTauEnergy(tau);
479 
480  if(visibleTauE != 0) rTau = ltrack/visibleTauE;
481 
482  if(abs(mother) == 24) TauRtauW->Fill(rTau,weight);
483  if(abs(mother) == 37) TauRtauHpm->Fill(rTau,weight);
484 }
485 
486 void TauValidation::spinEffects(const HepMC::GenParticle* tau,int mother, int decay, std::vector<HepMC::GenParticle*> &part,double weight){
487  if(decay == TauDecay::JAK_PION || decay == TauDecay::JAK_MUON || decay == TauDecay::JAK_ELECTRON){ // polarization only for 1-prong hadronic taus with no neutral pions
488  TLorentzVector momP4 = motherP4(tau);
489  TLorentzVector pionP4 = leadingPionP4(tau);
490  pionP4.Boost(-1*momP4.BoostVector());
491  double energy = pionP4.E()/(momP4.M()/2);
492  if(decay == TauDecay::JAK_PION){
493  if(abs(mother) == 24) TauSpinEffectsW_X->Fill(energy,weight);
494  if(abs(mother) == 37) TauSpinEffectsHpm_X->Fill(energy,weight);
495  }
496  if(decay == TauDecay::JAK_MUON){
497  if(abs(mother) == 24) TauSpinEffectsW_muX->Fill(energy,weight);
498  if(abs(mother) == 37) TauSpinEffectsHpm_muX->Fill(energy,weight);
499  }
500  if(decay == TauDecay::JAK_ELECTRON){
501  if(abs(mother) == 24) TauSpinEffectsW_eX->Fill(energy,weight);
502  if(abs(mother) == 37) TauSpinEffectsHpm_eX->Fill(energy,weight);
503  }
504 
505  }
506  else if(decay==TauDecay::JAK_RHO_PIPI0){
507  TLorentzVector rho(0,0,0,0),pi(0,0,0,0);
508  for(unsigned int i=0;i<part.size();i++){
509  TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
510  if(abs(part.at(i)->pdg_id())==PdtPdgMini::pi_plus){pi+=LV; rho+=LV;}
511  if(abs(part.at(i)->pdg_id())==PdtPdgMini::pi0){rho+=LV;}
512  }
513  if(abs(mother) == 24) TauSpinEffectsW_UpsilonRho->Fill(2*pi.Pt()/rho.Pt()-1,weight);
514  if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonRho->Fill(2*pi.Pt()/rho.Pt()-1,weight);
515  }
516  else if(decay==TauDecay::JAK_A1_3PI){ // only for pi2pi0 for now
517  TLorentzVector a1(0,0,0,0),pi_p(0,0,0,0),pi_m(0,0,0,0);
518  int nplus(0),nminus(0);
519  for(unsigned int i=0;i<part.size();i++){
520  TLorentzVector LV(part.at(i)->momentum().px(),part.at(i)->momentum().py(),part.at(i)->momentum().pz(),part.at(i)->momentum().e());
521  if(part.at(i)->pdg_id()==PdtPdgMini::pi_plus){ pi_p+=LV; a1+=LV; nplus++;}
522  if(part.at(i)->pdg_id()==PdtPdgMini::pi_minus){pi_m+=LV; a1+=LV; nminus++;}
523  }
524  double gamma=0;
525  if(nplus+nminus==3 && nplus==1) gamma=2*pi_p.Pt()/a1.Pt()-1;
526  if(nplus+nminus==3 && nminus==1) gamma=2*pi_m.Pt()/a1.Pt()-1;
527  else{
528  pi_p+=pi_m; gamma=2*pi_p.Pt()/a1.Pt()-1;
529  }
530  if(abs(mother) == 24) TauSpinEffectsW_UpsilonA1->Fill(gamma,weight);
531  if(abs(mother) == 37) TauSpinEffectsHpm_UpsilonA1->Fill(gamma,weight);
532  }
533 }
534 
536 
537  TLorentzVector tautau(0,0,0,0);
538  TLorentzVector pipi(0,0,0,0);
539  TLorentzVector taum(0,0,0,0);
540  int nSinglePionDecays(0),nSingleMuonDecays(0),nSingleElectronDecays(0);
541  double x1(0),x2(0);
542  TLorentzVector Zboson(boson->momentum().px(),boson->momentum().py(),boson->momentum().pz(),boson->momentum().e());
543  if ( boson->end_vertex() ) {
544  HepMC::GenVertex::particle_iterator des;
545  for(des = boson->end_vertex()->particles_begin(HepMC::children); des!= boson->end_vertex()->particles_end(HepMC::children);++des ) {
546  int pid = (*des)->pdg_id();
547  if(abs(findMother(*des)) != 15 && abs(pid) == 15 && (tauDecayChannel(*des) == pi || tauDecayChannel(*des) == muon || tauDecayChannel(*des) == electron )){
548  if(tauDecayChannel(*des) == pi)nSinglePionDecays++;
549  if(tauDecayChannel(*des) == muon)nSingleMuonDecays++;
550  if(tauDecayChannel(*des) == electron)nSingleElectronDecays++;
551  TLorentzVector LVtau((*des)->momentum().px(),(*des)->momentum().py(),(*des)->momentum().pz(),(*des)->momentum().e());
552  tautau += LVtau;
553  TLorentzVector LVpi=leadingPionP4(*des);
554  pipi+=LVpi;
555  const HepPDT::ParticleData* pd = fPDGTable->particle((*des)->pdg_id ());
556  int charge = (int) pd->charge();
557  LVtau.Boost(-1*Zboson.BoostVector());
558  LVpi.Boost(-1*Zboson.BoostVector());
559  if(charge<0){x1=LVpi.P()/LVtau.E(); taum=LVtau;}
560  else{ x2=LVpi.P()/LVtau.E();}
561  }
562  }
563  }
564  if(nSingleMuonDecays==2){
565  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_muX->Fill(x1,weight);
566  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_muX->Fill(x1,weight);
567  }
568  if(nSingleElectronDecays==2){
569  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_eX->Fill(x1,weight);
570  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_eX->Fill(x1,weight);
571  }
572  if(nSinglePionDecays == 2 && tautau.M()!= 0) {
573  for(int i=0;i<zsbins;i++){
574  double zslow=((double)i)*(zsmax-zsmin)/((double)zsbins)+zsmin;
575  double zsup=((double)i+1)*(zsmax-zsmin)/((double)zsbins)+zsmin;
576  double aup=Zstoa(zsup), alow=Zstoa(zslow);
577  if(x2-x1>alow && x2-x1<aup){
578  double zs=(zsup+zslow)/2;
579  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_Zs->Fill(zs,weight);
580  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Zs->Fill(zs,weight);
581  break;
582  }
583  }
584  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_MVis->Fill(pipi.M()/tautau.M(),weight);
585  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_MVis->Fill(pipi.M()/tautau.M(),weight);
586 
587  if(x1!=0){
588  const std::vector<HepMC::GenParticle*> m=GetMothers(boson);
589  int q(0),qbar(0);
590  TLorentzVector Z(0,0,0,0);
591  for(unsigned int i=0;i<m.size();i++){
592  if(m.at(i)->pdg_id()==PdtPdgMini::d || m.at(i)->pdg_id()==PdtPdgMini::u ){q++;}
593  if(m.at(i)->pdg_id()==PdtPdgMini::anti_d || m.at(i)->pdg_id()==PdtPdgMini::anti_u ){qbar++;}
594  }
595  if(q==1 && qbar==1){// assume q has largest E (valence vs see quarks)
596  if(taum.Vect().Dot(Zboson.Vect())/(Zboson.P()*taum.P())>0){
597  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_Xf->Fill(x1,weight);
598  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xf->Fill(x1,weight);
599  }
600  else{
601  if(abs(boson->pdg_id())==PdtPdgMini::Z0) TauSpinEffectsZ_Xb->Fill(x1,weight);
602  if(abs(boson->pdg_id())==PdtPdgMini::Higgs0) TauSpinEffectsH_Xb->Fill(x1,weight);
603  }
604  }
605  }
606  }
607 }
608 
609 double TauValidation::Zstoa(double zs){
610  double a=1-sqrt(fabs(1.0-2*fabs(zs)));
611  if(zs<0){
612  a*=-1.0;
613  }
614  return a;
615 }
616 
617 
619  return leadingPionP4(tau).P();
620 }
621 
623 
624  TLorentzVector p4(0,0,0,0);
625 
626  if ( tau->end_vertex() ) {
627  HepMC::GenVertex::particle_iterator des;
628  for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
629  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
630  int pid = (*des)->pdg_id();
631 
632  if(abs(pid) == 15) return leadingPionP4(*des);
633 
634  if(!(abs(pid)==211 || abs(pid)==13 || abs(pid)==11)) continue;
635 
636  if((*des)->momentum().rho() > p4.P()) {
637  p4 = TLorentzVector((*des)->momentum().px(),
638  (*des)->momentum().py(),
639  (*des)->momentum().pz(),
640  (*des)->momentum().e());
641  }
642  }
643  }
644 
645  return p4;
646 }
647 
649  const HepMC::GenParticle* m=GetMother(tau);
650  return TLorentzVector(m->momentum().px(),m->momentum().py(),m->momentum().pz(),m->momentum().e());
651 }
652 
654  TLorentzVector p4(tau->momentum().px(),
655  tau->momentum().py(),
656  tau->momentum().pz(),
657  tau->momentum().e());
658 
659  if ( tau->end_vertex() ) {
660  HepMC::GenVertex::particle_iterator des;
661  for(des = tau->end_vertex()->particles_begin(HepMC::descendants);
662  des!= tau->end_vertex()->particles_end(HepMC::descendants);++des ) {
663  int pid = (*des)->pdg_id();
664 
665  if(abs(pid) == 15) return visibleTauEnergy(*des);
666 
667  if(abs(pid) == 12 || abs(pid) == 14 || abs(pid) == 16) {
668  p4 -= TLorentzVector((*des)->momentum().px(),
669  (*des)->momentum().py(),
670  (*des)->momentum().pz(),
671  (*des)->momentum().e());
672  }
673  }
674  }
675 
676  return p4.E();
677 }
678 
680  // Find First tau in chain
681  std::vector<const HepMC::GenParticle*> TauList;
682  findTauList(tau,TauList);
683 
684  // Get List of Gauge Boson to tau(s) FSR and Brem
685  bool passedW=false;
686  std::vector<const HepMC::GenParticle*> ListofFSR; ListofFSR.clear();
687  std::vector<const HepMC::GenParticle*> ListofBrem; ListofBrem.clear();
688  std::vector<const HepMC::GenParticle*> FSR_photos; FSR_photos.clear();
689  double BosonScale(1);
690  if(TauList.size()>0){
691  TauValidation::findFSRandBrem(TauList.at(0),passedW,ListofFSR,ListofBrem);
692  TauValidation::FindPhotosFSR(TauList.at(0),FSR_photos,BosonScale);
693 
694  // Add the Tau Brem. information
695  TauBremPhotonsN->Fill(ListofBrem.size(),weight);
696  double photonPtSum=0;
697  for(unsigned int i=0;i<ListofBrem.size();i++){
698  photonPtSum+=ListofBrem.at(i)->momentum().perp();
699  TauBremPhotonsPt->Fill(ListofBrem.at(i)->momentum().perp(),weight);
700  }
701  TauBremPhotonsPtSum->Fill(photonPtSum,weight);
702 
703  // Now add the Gauge Boson FSR information
704  if(BosonScale!=0){
705  TauFSRPhotonsN->Fill(ListofFSR.size(),weight);
706  photonPtSum=0;
707  for(unsigned int i=0;i<ListofFSR.size();i++){
708  photonPtSum+=ListofFSR.at(i)->momentum().perp();
709  TauFSRPhotonsPt->Fill(ListofFSR.at(i)->momentum().perp(),weight);
710  }
711  double FSR_photosSum(0);
712  for(unsigned int i=0;i<FSR_photos.size();i++){
713  FSR_photosSum+=FSR_photos.at(i)->momentum().perp();
714  TauFSRPhotonsPt->Fill(FSR_photos.at(i)->momentum().perp()/BosonScale,weight*BosonScale);
715  }
716  TauFSRPhotonsPtSum->Fill(photonPtSum+FSR_photosSum/BosonScale,weight);
717  }
718  }
719 }
720 
void findTauList(const HepMC::GenParticle *tau, std::vector< const HepMC::GenParticle * > &TauList)
MonitorElement * TauFSRPhotonsN
int i
Definition: DBlmapReader.cc:9
virtual ~TauValidation()
MonitorElement * TauPhi
MonitorElement * TauSpinEffectsW_X
TPRegexp parents
Definition: eve_filter.cc:24
MonitorElement * TauSpinEffectsHpm_UpsilonA1
std::vector< std::vector< MonitorElement * > > JAKInvMass
MonitorElement * TauSpinEffectsZ_Xb
MonitorElement * nTaus
MonitorElement * book1D(const char *name, const char *title, int nchX, double lowX, double highX)
Book 1D histogram.
Definition: DQMStore.cc:722
MonitorElement * TauRtauW
TauValidation(const edm::ParameterSet &)
MonitorElement * TauPt
void spinEffects(const HepMC::GenParticle *, int, int, std::vector< HepMC::GenParticle * > &part, double weight)
MonitorElement * TauSpinEffectsH_MVis
MonitorElement * TauMothers
MonitorElement * TauSpinEffectsHpm_eX
MonitorElement * TauSpinEffectsW_muX
MonitorElement * TauBremPhotonsN
void setBinLabel(int bin, const std::string &label, int axis=1)
set bin label for x, y or z axis (axis=1, 2, 3 respectively)
#define abs(x)
Definition: mlp_lapack.h:159
MonitorElement * JAKID
MonitorElement * TauSpinEffectsZ_Xf
int tauProngs(const HepMC::GenParticle *, double weight)
MonitorElement * TauBremPhotonsPtSum
virtual void endJob()
void spinEffectsZ(const HepMC::GenParticle *boson, double weight)
MonitorElement * TauSpinEffectsHpm_UpsilonRho
double charge(const std::vector< uint8_t > &Ampls)
virtual void analyze(const edm::Event &, const edm::EventSetup &)
void getData(T &iHolder) const
Definition: EventSetup.h:67
MonitorElement * TauSpinEffectsZ_eX
int findMother(const HepMC::GenParticle *)
void Fill(long long x)
MonitorElement * TauProngs
int tauMother(const HepMC::GenParticle *, double weight)
int tauDecayChannel(const HepMC::GenParticle *, double weight=0.0)
bool isTauFinalStateParticle(int pdgid)
Definition: TauDecay.cc:38
const std::vector< HepMC::GenParticle * > GetMothers(const HepMC::GenParticle *boson)
math::XYZTLorentzVectorD LV
int iEvent
Definition: GenABIO.cc:243
MonitorElement * TauSpinEffectsW_UpsilonRho
void photons(const HepMC::GenParticle *, double weight)
T sqrt(T t)
Definition: SSEVec.h:48
double p4[4]
Definition: TauolaWrapper.h:92
double leadingPionMomentum(const HepMC::GenParticle *, double weight)
virtual void beginJob()
MonitorElement * TauSpinEffectsH_Xf
MonitorElement * TauFSRPhotonsPt
MonitorElement * TauSpinEffectsZ_muX
MonitorElement * TauDecayChannels
edm::InputTag hepmcCollection_
Definition: TauValidation.h:98
MonitorElement * TauSpinEffectsH_Zs
MonitorElement * TauFSRPhotonsPtSum
HepPDT::ParticleData ParticleData
double visibleTauEnergy(const HepMC::GenParticle *)
bool getByLabel(InputTag const &tag, Handle< PROD > &result) const
Definition: Event.h:361
TLorentzVector leadingPionP4(const HepMC::GenParticle *)
bool isLastTauinChain(const HepMC::GenParticle *tau)
MonitorElement * TauSpinEffectsZ_MVis
TLorentzVector motherP4(const HepMC::GenParticle *)
void rtau(const HepMC::GenParticle *, int, int, double weight)
WeightManager _wmanager
Definition: TauValidation.h:96
part
Definition: HCALResponse.h:20
void FindPhotosFSR(const HepMC::GenParticle *p, std::vector< const HepMC::GenParticle * > &ListofFSR, double &BosonScale)
virtual void endRun(const edm::Run &, const edm::EventSetup &)
edm::ESHandle< HepPDT::ParticleDataTable > fPDGTable
PDT table.
MonitorElement * TauSpinEffectsH_eX
MonitorElement * TauSpinEffectsW_eX
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
double Zstoa(double zs)
const HepMC::GenParticle * GetMother(const HepMC::GenParticle *tau)
MonitorElement * TauSpinEffectsH_muX
MonitorElement * TauSpinEffectsW_UpsilonA1
double a
Definition: hdecay.h:121
MonitorElement * TauEta
DQMStore * dbe
ME&#39;s &quot;container&quot;.
void findFSRandBrem(const HepMC::GenParticle *p, bool doBrem, std::vector< const HepMC::GenParticle * > &ListofFSR, std::vector< const HepMC::GenParticle * > &ListofBrem)
unsigned int NJAKID
MonitorElement * TauSpinEffectsZ_Zs
MonitorElement * TauSpinEffectsHpm_X
virtual void beginRun(const edm::Run &, const edm::EventSetup &)
int weight
Definition: histoStyle.py:50
MonitorElement * TauSpinEffectsHpm_muX
double weight(const edm::Event &)
void setAxisTitle(const std::string &title, int axis=1)
set x-, y- or z-axis title (axis=1, 2, 3 respectively)
MonitorElement * TauSpinEffectsH_Xb
std::vector< HepMC::GenParticle * > Get_TauDecayProducts()
bool AnalyzeTau(HepMC::GenParticle *Tau, unsigned int &JAK_ID, unsigned int &TauBitMask, bool dores=true, bool dopi0=true)
MonitorElement * nPrimeTaus
MonitorElement * TauBremPhotonsPt
MonitorElement * TauRtauHpm
void setCurrentFolder(const std::string &fullpath)
Definition: DQMStore.cc:434
Definition: Run.h:36
unsigned int nProng(unsigned int &TauBitMask)
Definition: TauDecay.h:104