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