StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StDijetFilter.cxx
1 #include "StDijetFilter.h"
2 #include "StGenParticle.h"
3 #include <string>
4 #include <iostream>
5 #include <fstream>
6 #include <cmath>
7 #include <vector>
8 #include <algorithm>
9 
10 static StDijetFilter dijetFilter;
11 
12 using namespace std;
13 //_______________________________________________________________
14 StDijetFilter::StDijetFilter():StMCFilter("dijet")
15 {
16 
17  mRBTOW = 225.405;
18  mRcone = 0.7;
19  mSeed = 0.5;
20  mAssoc = 0.0;
21  mSplitfraction = 0.5;
22  mAddmidpoints = true;
23  mStablemidpoints = true;
24  mSplitmerge = true;
25  nEvents = new int;
26  (*nEvents) = 0;
27 
28  mParticleEtaRange = 3.5;
29  mJetEtaHigh = 1.3;
30  mJetEtaLow = -1.3;
31  mPtLow = 7.0;
32  mPtHigh = 10.0;
33  mDPhi = 2.0;
34  mMinJetPt = 4.0;
35  mDEta = mJetEtaHigh - mJetEtaLow;
36 
37  mRecohadron = 1.1;
38  mRecolepton = 1.5;
39 
40  mVertex[0] = 0.;
41  mVertex[1] = 0.;
42  mVertex[2] = 0.;
43 
44  readConfig();
45 
46 }
47 //_______________________________________________________________
48 StDijetFilter::~StDijetFilter()
49 {
50  delete nEvents;
51 }
52 //_______________________________________________________________
53 int StDijetFilter::RejectEG(const StGenParticleMaster &ptl) const
54 {
55  (*nEvents)++;
56  return 0;
57 }//_______________________________________________________________
58 int StDijetFilter::RejectGT(const StGenParticleMaster &ptl) const
59 {
60  if(ptl.Size() <= 0)return 1;
61 
62  vector<JetFourVec*> finalparticles;
63  vector<JetFourVec*> seeds;
64  double v[3];
65 
66  //apply cuts on particles
67  //get seeds
68  for(int i = 0; i < ptl.Size(); i++){
69  if(ptl(i)->GetStatusCode() != 1)continue;
70  if(fabs(ptl(i)->Eta()) > mParticleEtaRange)continue;
71  JetFourVec* p = new JetFourVec((StGenParticle*)ptl(i));
72  finalparticles.push_back(p);
73  ptl(i)->Vertex(v);
74 
75  if(p->Pt() < mSeed)continue;
76  seeds.push_back(p);
77  }
78 
79 
80  if(mAddmidpoints){
81  seeds = addMidPoints(seeds);
82  }
83 
84  int nJets = 0;
85  vector< vector<JetFourVec*> > jetFour;
86  vector< vector<JetFourVec*> > recoJetFour;
87  jetFour.clear();
88  recoJetFour.clear();
89  for(unsigned int k = 0; k < seeds.size(); k++){
90  JetFourVec* j = new JetFourVec(seeds[k]);
91  vector<JetFourVec*> jf;
92  int nChange = 1;
93  int nIter = 0;
94  while(nChange != 0 && nIter < 10){
95  nChange = 0;
96  vector<JetFourVec*>::iterator jfiter;
97  for(unsigned int l = 0; l < finalparticles.size(); l++){
98  if(finalparticles[l]->getEn()*sin(finalparticles[l]->Theta()) < mAssoc)continue;
99 
100 
101  jfiter = find(jf.begin(),jf.end(),finalparticles[l]);
102  if(dR(finalparticles[l],j) > mRcone && jfiter == jf.end())continue;
103  if(dR(finalparticles[l],j) > mRcone && jfiter != jf.end()){
104  jf.erase(jfiter);
105  nChange++;
106  continue;
107  }
108  if(jfiter != jf.end())continue;
109  jf.push_back(finalparticles[l]);
110  nChange++;
111  }
112  if (j) delete j; // plug leak
113  j = combineTracks(jf);
114  nIter++;
115  }
116  if(jf.size() == 0){
117  if(j) delete j;
118  continue;
119  }
120  jetFour.push_back(jf);
121 
122  if (j) delete j; // plug leak
123  }
124 
125  jetFour = EtOrderedList(jetFour);
126  jetFour = RemoveDuplicates(jetFour);
127 
128  if(mSplitmerge){
129  jetFour = doSplitMerge(jetFour);
130  jetFour = RemoveDuplicates(jetFour);
131  }
132 
133 
134  nJets = 0;
135  JetFourVec* j = 0;
136  vector<JetFourVec*> dijetFour;
137  vector<JetFourVec*> dijetFourReco;
138  JetFourVec* dj = 0;
139  for(unsigned int k = 0; k < jetFour.size(); k++){
140  if(j)delete j;
141  j = combineTracks(jetFour[k]);
142  //cout<<"jet "<<k<<" "<<j->Pt()<<" "<<j->Eta()<<" "<<j->Phi()<<endl;
143  if(j->Pt() < mMinJetPt)continue;
144  if(fabs(j->Eta()) > mJetEtaHigh)continue;
145  nJets++;
146  if(j->Eta() < mJetEtaLow)continue;
147  dj = combineTracks(jetFour[k]);
148  dijetFour.push_back(dj);
149  }
150  if(j)delete j;
151  int nTrJets = 0;
152  dj = 0;
153  for(unsigned int k = 0; k < jetFour.size(); k++){
154  JetFourVec* j = recoJet(jetFour[k],v);
155  if(j->getEn()*sin(j->Theta()) > 4 && fabs(j->Eta()) < 2.5) nTrJets++;
156  if(j->Eta() < mJetEtaHigh && j->Eta() > mJetEtaLow){
157  dj = new JetFourVec(j);
158  dijetFourReco.push_back(dj);
159  }
160  delete j;
161  }
162 
163  dijetFourReco = EtOrderedList(dijetFourReco);
164  int yesdijet = 0;
165  if(dijetFour.size() > 1){
166  dj = new JetFourVec;
167  (*dj) = *(dijetFour[0]) + *(dijetFour[1]);
168  float dphi = fabs(dijetFour[0]->Phi() - dijetFour[1]->Phi());
169  float deta = fabs(dijetFour[0]->Eta() - dijetFour[1]->Eta());
170  if(dphi > acos(-1))dphi = fabs(2*acos(-1) - dphi);
171  if(dphi > mDPhi && deta < mDEta && dj->M() > 10. && dijetFour[0]->Pt() > mPtHigh && dijetFour[1]->Pt() > mPtLow)yesdijet += 1;
172  delete dj;
173  dj = 0;
174  }
175 
176  if(dijetFourReco.size() > 1){
177  dj = new JetFourVec;
178  (*dj) = *(dijetFourReco[0]) + *(dijetFourReco[1]);
179  float dphi = fabs(dijetFourReco[0]->Phi() - dijetFourReco[1]->Phi());
180  float deta = fabs(dijetFourReco[0]->Eta() - dijetFourReco[1]->Eta());
181  if(dphi > acos(-1))dphi = fabs(2*acos(-1) - dphi);
182  if(dphi > mDPhi && deta < mDEta && dj->M() > 10. && dijetFourReco[0]->Pt() > mPtHigh && dijetFourReco[1]->Pt() > mPtLow)yesdijet += 2;
183  delete dj;
184  dj = 0;
185  }
186  cout << "\n" <<"agrdl: event "<<(*nEvents)<<" "<<nJets<<" "<<nTrJets<<" "<<yesdijet<<endl;
187  for(unsigned int k = 0; k < finalparticles.size(); k++)
188  delete finalparticles[k];
189  finalparticles.clear();
190 
191  if(yesdijet == 0) return 1;
192 
193  return 0;
194 }
195 //_______________________________________________________________
196 int StDijetFilter::RejectGE(const StGenParticleMaster &ptl) const
197 {
198  return 0;
199 }
200 //_______________________________________________________________
201 void StDijetFilter::readConfig()
202 {
203  ifstream cfile("dijet.cnf");
204  while(1){
205  string attr;
206  float val;
207  cfile >> attr >> val;
208  if(!cfile.good())break;
209  parseConfig(attr,val);
210  }
211 }
212 //_______________________________________________________________
213 void StDijetFilter::parseConfig(string attr, float val)
214 {
215  if(attr == "RBTOW"){
216  mRBTOW = val;
217  cout<<"RBTOW changed to "<<val<<endl;
218  return;
219  }
220  if(attr == "RCONE"){
221  mRcone = val;
222  cout<<"RCONE changed to "<<val<<endl;
223  return;
224  }
225  if(attr == "SEED"){
226  mSeed = val;
227  cout<<"SEED changed to "<<val<<endl;
228  return;
229  }
230  if(attr == "ASSOC"){
231  mAssoc = val;
232  cout<<"ASSOC changed to "<<val<<endl;
233  return;
234  }
235  if(attr == "SPLIT"){
236  mSplitfraction = val;
237  cout<<"SPLIT changed to "<<val<<endl;
238  return;
239  }
240  if(attr == "MIDPOINT"){
241  if(val == 1.0){
242  mAddmidpoints = true;
243  cout<<"MIDPOINT changed to TRUE"<<endl;
244  return;
245  }else if(val == 0.0){
246  mAddmidpoints = false;
247  cout<<"MIDPOINT changed to FALSE"<<endl;
248  return;
249  }else{
250  cout<<"IMPRORER INPUT TO MIDPOINTS"<<endl;
251  }
252  }
253  if(attr == "STABLEMIDPOINTS"){
254  if(val == 1.0){
255  mStablemidpoints = true;
256  cout<<"STABLEMIDPOINTS changed to TRUE"<<endl;
257  return;
258  }else if(val == 0.0){
259  mStablemidpoints = false;
260  cout<<"STABLEMIDPOINTS changed to FALSE"<<endl;
261  return;
262  }else{
263  cout<<"IMPRORER INPUT TO STABLEMIDPOINTS"<<endl;
264  }
265  }
266  if(attr == "SPLITMERGE"){
267  if(val == 1.0){
268  mSplitmerge = true;
269  cout<<"SPLITMERGE changed to TRUE"<<endl;
270  return;
271  }else if(val == 0.0){
272  mSplitmerge = false;
273  cout<<"SPLITMERGE changed to FALSE"<<endl;
274  return;
275  }else{
276  cout<<"IMPRORER INPUT TO SPLITMERGE"<<endl;
277  }
278  }
279  if(attr == "PARTICLEETA"){
280  mParticleEtaRange = val;
281  cout<<"PARTICLEETA changed to "<<val<<endl;
282  return;
283  }
284  if(attr == "JETETAHIGH"){
285  mJetEtaHigh = val;
286  cout<<"JETETAHIGH changed to "<<val<<endl;
287  return;
288  }
289  if(attr == "JETETALOW"){
290  mJetEtaLow = val;
291  cout<<"JETETALOW changed to "<<val<<endl;
292  return;
293  }
294  if(attr == "DIJETPTLOW"){
295  mPtLow = val;
296  cout<<"DIJETPTLOW changed to "<<val<<endl;
297  return;
298  }
299  if(attr == "DIJETPTHIGH"){
300  mPtHigh = val;
301  cout<<"DIJETPTHIGH changed to "<<val<<endl;
302  return;
303  }
304  if(attr == "DPHI"){
305  mDPhi = val;
306  cout<<"DPHI changed to "<<val<<endl;
307  return;
308  }
309  if(attr == "MINJETPT"){
310  mMinJetPt = val;
311  cout<<"MINJETPT changed to "<<val<<endl;
312  return;
313  }
314  if(attr == "DETA"){
315  mDEta = val;
316  cout<<"DETA changed to "<<val<<endl;
317  return;
318  }
319  if(attr == "RECOHADRON"){
320  mRecohadron = val;
321  cout<<"RECOHADRON changed to "<<val<<endl;
322  return;
323  }
324  if(attr == "RECOLEPTON"){
325  mRecolepton = val;
326  cout<<"RECOLEPTON changed to "<<val<<endl;
327  return;
328  }
329 }
330 //_______________________________________________________________
331 JetFourVec* StDijetFilter::recoJet(vector<JetFourVec*> v, double* vert) const
332 {
333  JetFourVec* j = new JetFourVec;
334  for(unsigned int i = 0; i < v.size(); i++){
335  int id = v[i]->getCode();
336  if(abs(id) != 11 && id != -2212 && id != -2122){
337  v[i]->setPtEtaPhiM(mRecohadron*v[i]->Pt(),v[i]->Eta(),v[i]->Phi(),v[i]->M());
338  }else{
339  v[i]->setPtEtaPhiM(mRecolepton*v[i]->Pt(),v[i]->Eta(),v[i]->Phi(),v[i]->M());
340  }
341  (*j) = (*j) + (*v[i]);
342  }
343 
344  return j;
345 }
346 //_______________________________________________________________
347 vector<JetFourVec*> StDijetFilter::addMidPoints(vector<JetFourVec*> seeds) const
348 {
349  unsigned int n = seeds.size();
350  float pi = acos(-1.0);
351  for(unsigned int k = 0; k < n; k++){
352  for(unsigned int l = 0; l < k; l++){
353  if(dR(seeds[k],seeds[l]) > 2*mRcone)continue;
354  JetFourVec* s = new JetFourVec;
355  float maxphi = max(seeds[k]->Phi(),seeds[l]->Phi());
356  float minphi = min(seeds[k]->Phi(),seeds[l]->Phi());
357  float dphi = fabs(maxphi - minphi);
358  if(dphi > pi)minphi += 2*pi;
359  float neta = 0.5*(seeds[k]->Eta() + seeds[l]->Eta());
360  float nphi = 0.5*(maxphi + minphi);
361  if(nphi > pi)nphi -= 2*pi;
362  if(nphi < -pi) nphi += 2*pi;
363  s->setPtEtaPhiM(0.001,neta,nphi,0);
364  seeds.push_back(s);
365  }
366  }
367  return seeds;
368 }
369 //_______________________________________________________________
370 void StDijetFilter::split(vector<JetFourVec*> &v1, vector<JetFourVec*> &v2) const
371 {
372  JetFourVec* j1 = combineTracks(v1);
373  JetFourVec* j2 = combineTracks(v2);
374 
375  vector<JetFourVec*>::iterator dit;
376  for(vector<JetFourVec*>::iterator it = v1.begin(); it != v1.end(); it++){
377  dit = find(v2.begin(),v2.end(),*it);
378  if(dit == v2.end())continue;
379  float d1 = dR(j1,*it);
380  float d2 = dR(j2,*it);
381  if(d2 > d1){
382  v2.erase(dit);
383  }
384  }
385  for(vector<JetFourVec*>::iterator it = v2.begin(); it != v2.end(); it++){
386  dit = find(v1.begin(),v1.end(),*it);
387  if(dit == v1.end())continue;
388  float d1 = dR(j1,*it);
389  float d2 = dR(j2,*it);
390  if(d1 > d2){
391  v1.erase(dit);
392  }
393  }
394 
395  delete j1;
396  delete j2;
397 }
398 //_______________________________________________________________
399 vector<JetFourVec*> StDijetFilter::merge(vector<JetFourVec*> v1, vector<JetFourVec*> v2) const
400 {
401  vector<JetFourVec*> nj = v1;
402  vector<JetFourVec*>::iterator iter;
403  for(iter = v2.begin(); iter != v2.end(); iter++){
404  vector<JetFourVec*>::iterator it2 = find(v1.begin(),v1.end(),*iter);
405  if(it2 == v1.end()) nj.push_back(*iter);
406  }
407 
408  return nj;
409 }
410 //_______________________________________________________________
411 vector< vector<JetFourVec*> > StDijetFilter::doSplitMerge(vector< vector<JetFourVec*> >orig) const
412 {
413  vector< vector<JetFourVec*> > njf;
414  vector< vector<JetFourVec* > > jetFour = orig;
415  int nChange = 0;
416  int nIter = 0;
417  while(1){
418  if(nIter > 200)break;
419  nChange = 0;
420  njf.clear();
421  njf = jetFour;
422  vector< vector<JetFourVec*> >::iterator iter1;
423  vector< vector<JetFourVec*> >::iterator iter2;
424  for(iter1 = njf.begin(); iter1 != njf.end(); iter1++){
425  if(nChange)continue;
426  if((*iter1).size() == 0)continue;
427  JetFourVec* j = combineTracks(*iter1);
428  for(iter2 = iter1+1; iter2 != njf.end(); iter2++){
429  if(nChange)continue;
430  if((*iter2).size() == 0)continue;
431  JetFourVec* nj = combineTracks(*iter2);
432  if(dR(nj,j) > 2*mRcone){
433  delete nj;
434  continue;
435  }
436  float oe = overlapEnergy(*iter1,*iter2);
437  if(oe == 0){
438  delete nj;
439  continue;
440  }
441  nChange = 1;
442  vector< vector<JetFourVec*> >::iterator njit1;
443  vector< vector<JetFourVec*> >::iterator njit2;
444  njit1 = find(jetFour.begin(),jetFour.end(),*iter1);
445  njit2 = find(jetFour.begin(),jetFour.end(),*iter2);
446  if ( njit1 != jetFour.end() && njit2 != jetFour.end() ) {
447  if(oe > mSplitfraction){
448  vector<JetFourVec*> mj = merge(*iter1,*iter2);
449  (*njit1).clear();
450  (*njit2).clear();
451  jetFour.erase(njit2);
452  jetFour.erase(njit1);
453  jetFour.insert(jetFour.begin(),mj);
454  if(nj) delete nj;
455  continue;
456  }else{
457  split(*njit1,*njit2);
458  if(nj) delete nj;
459  continue;
460  }}
461  delete nj;
462  }
463  delete j;
464  }
465  jetFour = EtOrderedList(jetFour);
466  if(!nChange)break;
467  nIter++;
468  }
469  return njf;
470 
471 }
472 //_______________________________________________________________
473 float StDijetFilter::overlapEnergy(vector<JetFourVec*> v1,vector<JetFourVec*> v2) const
474 {
475  float e = 0;
476  JetFourVec* j1 = combineTracks(v1);
477  JetFourVec* j2 = combineTracks(v2);
478  float high = min(j1->getEn()*sin(j1->Theta()),j2->getEn()*sin(j2->Theta()));
479 
480  for(vector<JetFourVec*>::iterator iter1 = v1.begin(); iter1 != v1.end(); iter1++){
481  vector<JetFourVec*>::iterator iter2 = find(v2.begin(),v2.end(),*iter1);
482  if(iter2 == v2.end())continue;
483  e += (*iter2)->getEn()*sin((*iter2)->Theta());
484  }
485 
486  delete j1;
487  delete j2;
488 
489  return e/high;
490 
491 }
492 //_______________________________________________________________
493 vector<JetFourVec*> StDijetFilter::EtOrderedList(vector<JetFourVec*> jetFour) const
494 {
495  vector<JetFourVec*> njf;
496  vector<JetFourVec*>::iterator jfvIter;
497  vector<JetFourVec*>::iterator njfvIter;
498  for(jfvIter = jetFour.begin(); jfvIter != jetFour.end(); jfvIter++){
499  if((*jfvIter)->Pt() == 0)continue;
500  int added = 0;
501  for(njfvIter = njf.begin(); njfvIter != njf.end(); njfvIter++){
502  if((*jfvIter)->Pt() > (*njfvIter)->Pt()){
503  njf.insert(njfvIter,*jfvIter);
504  added = 1;
505  break;
506  }
507  }
508  if(!added){
509  njf.push_back(*jfvIter);
510  }
511  }
512 
513  return njf;
514 }
515 //_______________________________________________________________
516 vector< vector<JetFourVec*> > StDijetFilter::RemoveDuplicates(vector< vector<JetFourVec*> >jetFour) const
517 {
518  vector< vector<JetFourVec*> > njf;
519  vector< vector<JetFourVec*> >::iterator jfvIter;
520  vector< vector<JetFourVec*> >::iterator njfvIter;
521  for(jfvIter = jetFour.begin(); jfvIter != jetFour.end(); jfvIter++){
522  if((*jfvIter).size() == 0)continue;
523  JetFourVec* j = combineTracks(*jfvIter);
524  int dupe = 0;
525  for(njfvIter = njf.begin(); njfvIter != njf.end(); njfvIter++){
526  JetFourVec* nj = combineTracks(*njfvIter);
527  if(*nj == *j)dupe = 1;
528  delete nj;
529  }
530  if(!dupe)njf.push_back(*jfvIter);
531  delete j;
532  }
533  return njf;
534 }
535 //_______________________________________________________________
536 vector< vector<JetFourVec*> > StDijetFilter::EtOrderedList(vector< vector<JetFourVec*> >jetFour) const
537 {
538  vector< vector<JetFourVec*> > njf;
539  vector< vector<JetFourVec*> >::iterator jfvIter;
540  vector< vector<JetFourVec*> >::iterator njfvIter;
541  for(jfvIter = jetFour.begin(); jfvIter != jetFour.end(); jfvIter++){
542  if((*jfvIter).size() == 0)continue;
543  JetFourVec* j = combineTracks(*jfvIter);
544  if(j->Pt() == 0){ delete j;continue;}
545  int added = 0;
546  for(njfvIter = njf.begin(); njfvIter != njf.end(); njfvIter++){
547  JetFourVec* nj = combineTracks(*njfvIter);
548  if(j->Pt() > nj->Pt()){
549  njf.insert(njfvIter,*jfvIter);
550  delete nj;
551  added = 1;
552  break;
553  }
554  delete nj;
555  }
556  if(!added){
557  njf.push_back(*jfvIter);
558  }
559  delete j;
560  }
561 
562  return njf;
563 
564 }
565 //_______________________________________________________________
566 float StDijetFilter::dR(StGenParticle* p1,StGenParticle* p2) const
567 {
568  float dphi = fabs(p1->Phi() - p2->Phi());
569  float deta = p1->Eta() - p2->Eta();
570  float pi = acos(-1.0);
571  if(dphi > pi)dphi = 2*pi - dphi;
572  return sqrt(pow(deta,2) + pow(dphi,2));
573 }
574 //_______________________________________________________________
575 float StDijetFilter::dR(StGenParticle* p,JetFourVec* v) const
576 {
577  float dphi = fabs(p->Phi() - v->Phi());
578  float deta = p->Eta() - v->Eta();
579  float pi = acos(-1.0);
580  if(dphi > pi)dphi = 2*pi - dphi;
581  return sqrt(pow(deta,2) + pow(dphi,2));
582 }
583 //_______________________________________________________________
584 float StDijetFilter::dR(JetFourVec* v,StGenParticle* p) const
585 {
586  return dR(p,v);
587 }
588 //_______________________________________________________________
589 float StDijetFilter::dR(JetFourVec* v1, JetFourVec* v2) const
590 {
591  float dphi = fabs(v1->Phi() - v2->Phi());
592  float deta = v1->Eta() - v2->Eta();
593  float pi = acos(-1.0);
594  if(dphi > pi)dphi = 2*pi - dphi;
595  return sqrt(pow(deta,2) + pow(dphi,2));
596 }
597 //_______________________________________________________________
598 JetFourVec* StDijetFilter::combineTracks(vector<JetFourVec*> jf) const
599 {
600  JetFourVec* j = new JetFourVec;
601  if(jf.size() == 0)return j;
602  for(unsigned int i = 0; i < jf.size(); i++){
603  (*j) = (*j) + (*jf[i]);
604  }
605 
606  return j;
607 
608 }
609 //_______________________________________________________________
610 JetFourVec::JetFourVec()
611 {
612  px = 0,py = 0,pz = 0,en = 0;code = 0;
613 }
614 //_______________________________________________________________
615 JetFourVec::JetFourVec(JetFourVec* v)
616 {
617  code = v->getCode();
618  px = v->getPx();
619  py = v->getPy();
620  pz = v->getPz();
621  en = v->getEn();
622 }
623 //_______________________________________________________________
624 JetFourVec::JetFourVec(StGenParticle* g)
625 {
626  double p[4];
627  g->Momentum(p);
628  px = p[0];
629  py = p[1];
630  pz = p[2];
631  en = p[3];
632  code = g->GetPdgCode();
633 }
634 //_______________________________________________________________
635 void JetFourVec::setPxPyPzEn(float x, float y, float z, float e)
636 {
637  setPx(x);
638  setPy(y);
639  setPz(z);
640  setEn(e);
641 }
642 //_______________________________________________________________
643 float JetFourVec::P()
644 {
645  return sqrt(px*px + py*py + pz*pz);
646 }
647 //_______________________________________________________________
648 float JetFourVec::Pt()
649 {
650  return sqrt(px*px + py*py);
651 }
652 //_______________________________________________________________
653 float JetFourVec::Eta()
654 {
655  float pmom = P();
656  if(pmom > fabs(pz))return -log(tan(Theta()/2.));
657  else return 1.e30;
658 }
659 //_______________________________________________________________
660 float JetFourVec::Phi()
661 {
662  return atan2(py,px);
663 }
664 //_______________________________________________________________
665 float JetFourVec::Theta()
666 {
667  return acos(pz/P());
668 }
669 //_______________________________________________________________
670 float JetFourVec::M()
671 {
672  if(en*en - P()*P() < 0)return 0;
673  return sqrt(en*en-P()*P());
674 }
675 //_______________________________________________________________
676 void JetFourVec::setPtEtaPhiM(float pt, float eta, float phi, float m)
677 {
678  float x = pt * cos(phi);
679  float y = pt * sin(phi);
680  float p = pt*cosh(eta);
681  float z = sqrt(p*p-pt*pt)*eta/fabs(eta);
682  float e = sqrt(m*m + p*p);
683  setPxPyPzEn(x,y,z,e);
684 }
685 //_______________________________________________________________
686 void JetFourVec::setPtEtaPhiE(float pt, float eta, float phi, float e)
687 {
688  float x = pt * cos(phi);
689  float y = pt * sin(phi);
690  float p = pt*cosh(eta);
691  float z = sqrt(p*p-pt*pt)*eta/fabs(eta);
692  setPxPyPzEn(x,y,z,e);
693 }
694 //_______________________________________________________________
695 JetFourVec JetFourVec::operator +(JetFourVec x)
696 {
697  JetFourVec t;
698  t.setCode(0);
699  t.setPx(x.getPx() + px);
700  t.setPy(x.getPy() + py);
701  t.setPz(x.getPz() + pz);
702  t.setEn(x.getEn() + en);
703 
704  return t;
705 }
706 //_______________________________________________________________
707 bool JetFourVec::operator ==(JetFourVec x)
708 {
709  if(code != x.getCode())return false;
710  if(fabs(px - x.getPx()) > 1e-4)return false;
711  if(fabs(py - x.getPy()) > 1e-4)return false;
712  if(fabs(pz - x.getPz()) > 1e-4)return false;
713  if(fabs(en - x.getEn()) > 1e-4)return false;
714 
715  return true;
716 }
JetFourVec::getPz
float getPz()
getter for py
Definition: StDijetFilter.h:55
StDijetFilter::~StDijetFilter
virtual ~StDijetFilter()
constructor
Definition: StDijetFilter.cxx:46
JetFourVec::getEn
float getEn()
getter for pz
Definition: StDijetFilter.h:56
JetFourVec::Eta
float Eta()
calculate vector pt
Definition: StDijetFilter.cxx:647
StGenParticle
Abstract base class for particles related to common /HEPEVT/.
Definition: StGenParticle.h:133
JetFourVec::getPx
float getPx()
getter for pdg code
Definition: StDijetFilter.h:53
JetFourVec::JetFourVec
JetFourVec()
pdg code of four vector
Definition: StDijetFilter.cxx:606
JetFourVec::M
float M()
calculate vector theta
Definition: StDijetFilter.cxx:664
JetFourVec::setPx
void setPx(float x)
setter for pdg code
Definition: StDijetFilter.h:44
JetFourVec::setPxPyPzEn
void setPxPyPzEn(float, float, float, float)
setter for en
Definition: StDijetFilter.cxx:629
JetFourVec::setCode
void setCode(int x)
equality operator
Definition: StDijetFilter.h:43
JetFourVec::operator+
JetFourVec operator+(JetFourVec)
destructor
Definition: StDijetFilter.cxx:689
JetFourVec::Theta
float Theta()
calculate vector phi
Definition: StDijetFilter.cxx:659
JetFourVec::P
float P()
calculate vector mass
Definition: StDijetFilter.cxx:637
JetFourVec::setPz
void setPz(float x)
setter for py
Definition: StDijetFilter.h:46
JetFourVec::setPy
void setPy(float x)
setter for px
Definition: StDijetFilter.h:45
JetFourVec::setPtEtaPhiM
void setPtEtaPhiM(float, float, float, float)
four element setter
Definition: StDijetFilter.cxx:670
JetFourVec::operator==
bool operator==(JetFourVec)
addition operator
Definition: StDijetFilter.cxx:701
JetFourVec::getCode
int getCode()
alternative four element setter
Definition: StDijetFilter.h:52
StDijetFilter::RejectEG
int RejectEG(const StGenParticleMaster &ptl) const
destructor
Definition: StDijetFilter.cxx:53
StDijetFilter::parseConfig
void parseConfig(std::string, float)
Definition: StDijetFilter.cxx:207
JetFourVec::Phi
float Phi()
calculate vector pseudorapidity
Definition: StDijetFilter.cxx:654
JetFourVec::setPtEtaPhiE
void setPtEtaPhiE(float, float, float, float)
alternative four element setter
Definition: StDijetFilter.cxx:680
StDijetFilter::RejectGT
int RejectGT(const StGenParticleMaster &ptl) const
Rejection of GEANT Tracking.
Definition: StDijetFilter.cxx:58
JetFourVec::setEn
void setEn(float x)
setter for pz
Definition: StDijetFilter.h:47
JetFourVec::Pt
float Pt()
getter for en
Definition: StDijetFilter.cxx:642
JetFourVec::getPy
float getPy()
getter for px
Definition: StDijetFilter.h:54
StDijetFilter::StDijetFilter
StDijetFilter()
read a config file to adjust parameters
Definition: StDijetFilter.cxx:12
StDijetFilter::RejectGE
int RejectGE(const StGenParticleMaster &ptl) const
Rejection at GEANT End, No GEANT output.
Definition: StDijetFilter.cxx:196
Generated by   doxygen 1.8.5