StRoot: StMCFilter/StDijetFilter.cxx Source File

00001 #include "StDijetFilter.h"
00002 #include "StGenParticle.h"
00003 #include <string>
00004 #include <iostream>
00005 #include <fstream>
00006 #include <cmath>
00007 #include <vector>
00008 #include <algorithm>
00009 
00010 static StDijetFilter dijetFilter;
00011 
00012 using namespace std;
00013 //_______________________________________________________________
00014 StDijetFilter::StDijetFilter():StMCFilter("dijet")
00015 {
00016 
00017   mRBTOW = 225.405;
00018   mRcone = 0.7;
00019   mSeed = 0.5;
00020   mAssoc = 0.0;
00021   mSplitfraction = 0.5;
00022   mAddmidpoints = true;
00023   mStablemidpoints = true;
00024   mSplitmerge = true;
00025   nEvents = new int;
00026   (*nEvents) = 0;
00027 
00028   mParticleEtaRange = 3.5;
00029   mJetEtaHigh = 1.3;
00030   mJetEtaLow = -1.3;
00031   mPtLow = 7.0;
00032   mPtHigh = 10.0;
00033   mDPhi = 2.0;
00034   mMinJetPt = 4.0;
00035   mDEta = mJetEtaHigh - mJetEtaLow;
00036 
00037   mRecohadron = 1.1;
00038   mRecolepton = 1.5;
00039 
00040   mVertex[0] = 0.;
00041   mVertex[1] = 0.;
00042   mVertex[2] = 0.;
00043 
00044   readConfig();
00045 
00046 }
00047 //_______________________________________________________________
00048 StDijetFilter::~StDijetFilter()
00049 {
00050   delete nEvents;
00051 }
00052 //_______________________________________________________________
00053 int StDijetFilter::RejectEG(const StGenParticleMaster &ptl) const
00054 {
00055   (*nEvents)++;
00056   return 0;
00057 }//_______________________________________________________________
00058 int StDijetFilter::RejectGT(const StGenParticleMaster &ptl) const
00059 {
00060   if(ptl.Size() <= 0)return 1;
00061 
00062   vector<JetFourVec*> finalparticles;
00063   vector<JetFourVec*> seeds;
00064   double v[3];
00065 
00066   //apply cuts on particles
00067   //get seeds
00068   for(int i = 0; i < ptl.Size(); i++){
00069     if(ptl(i)->GetStatusCode() != 1)continue;
00070     if(fabs(ptl(i)->Eta()) > mParticleEtaRange)continue;
00071     JetFourVec* p = new JetFourVec((StGenParticle*)ptl(i));
00072     finalparticles.push_back(p);
00073     ptl(i)->Vertex(v);
00074 
00075     if(p->Pt() < mSeed)continue;
00076     seeds.push_back(p);
00077   }
00078 
00079 
00080   if(mAddmidpoints){
00081     seeds = addMidPoints(seeds);
00082   }
00083 
00084   int nJets = 0;
00085   vector< vector<JetFourVec*> > jetFour;
00086   vector< vector<JetFourVec*> > recoJetFour;
00087   jetFour.clear();
00088   recoJetFour.clear();
00089   for(unsigned int k = 0; k < seeds.size(); k++){
00090     JetFourVec* j = new JetFourVec(seeds[k]);
00091     vector<JetFourVec*> jf;
00092     int nChange = 1;
00093     int nIter = 0;
00094     while(nChange != 0 && nIter < 10){
00095       nChange = 0;
00096       vector<JetFourVec*>::iterator jfiter;
00097       for(unsigned int l = 0; l < finalparticles.size(); l++){
00098         if(finalparticles[l]->getEn()*sin(finalparticles[l]->Theta()) < mAssoc)continue;
00099 
00100         
00101         jfiter = find(jf.begin(),jf.end(),finalparticles[l]);
00102         if(dR(finalparticles[l],j) > mRcone && jfiter == jf.end())continue;
00103         if(dR(finalparticles[l],j) > mRcone && jfiter != jf.end()){
00104           jf.erase(jfiter);
00105           nChange++;
00106           continue;
00107         }
00108         if(jfiter != jf.end())continue;
00109         jf.push_back(finalparticles[l]);
00110         nChange++;
00111       }
00112       j = combineTracks(jf);
00113       nIter++;
00114     }
00115     if(jf.size() == 0)continue;
00116     jetFour.push_back(jf);
00117   }
00118 
00119   jetFour = EtOrderedList(jetFour);
00120   jetFour = RemoveDuplicates(jetFour);
00121 
00122   if(mSplitmerge){
00123     jetFour = doSplitMerge(jetFour);
00124     jetFour = RemoveDuplicates(jetFour);
00125   }
00126 
00127 
00128   nJets = 0;
00129   JetFourVec* j = 0;
00130   vector<JetFourVec*> dijetFour;
00131   vector<JetFourVec*> dijetFourReco;
00132   JetFourVec* dj = 0;
00133   for(unsigned int k = 0; k < jetFour.size(); k++){
00134     if(j)delete j;
00135     j = combineTracks(jetFour[k]);
00136     //cout<<"jet "<<k<<" "<<j->Pt()<<" "<<j->Eta()<<" "<<j->Phi()<<endl;
00137     if(j->Pt() < mMinJetPt)continue;
00138     if(fabs(j->Eta()) > mJetEtaHigh)continue;
00139     nJets++;
00140     if(j->Eta() < mJetEtaLow)continue;
00141     dj = combineTracks(jetFour[k]);
00142     dijetFour.push_back(dj);
00143   }
00144   if(j)delete j;
00145   int nTrJets = 0;
00146   dj = 0;
00147   for(unsigned int k = 0; k < jetFour.size(); k++){
00148     JetFourVec* j = recoJet(jetFour[k],v);
00149     if(j->getEn()*sin(j->Theta()) > 4 && fabs(j->Eta()) < 2.5) nTrJets++;
00150     if(j->Eta() < mJetEtaHigh && j->Eta() > mJetEtaLow){
00151       dj = new JetFourVec(j);
00152       dijetFourReco.push_back(dj);
00153     }
00154     delete j;
00155   }
00156 
00157   dijetFourReco = EtOrderedList(dijetFourReco);
00158   int yesdijet = 0;
00159   if(dijetFour.size() > 1){
00160     dj = new JetFourVec;
00161     (*dj) = *(dijetFour[0]) + *(dijetFour[1]);
00162     float dphi = fabs(dijetFour[0]->Phi() - dijetFour[1]->Phi());
00163     float deta = fabs(dijetFour[0]->Eta() - dijetFour[1]->Eta());
00164     if(dphi > acos(-1))dphi = fabs(2*acos(-1) - dphi);
00165     if(dphi > mDPhi && deta < mDEta && dj->M() > 10. && dijetFour[0]->Pt() > mPtHigh && dijetFour[1]->Pt() > mPtLow)yesdijet += 1;
00166     delete dj;
00167     dj = 0;
00168   }
00169 
00170   if(dijetFourReco.size() > 1){
00171     dj = new JetFourVec;
00172     (*dj) = *(dijetFourReco[0]) + *(dijetFourReco[1]);
00173     float dphi = fabs(dijetFourReco[0]->Phi() - dijetFourReco[1]->Phi());
00174     float deta = fabs(dijetFourReco[0]->Eta() - dijetFourReco[1]->Eta());
00175     if(dphi > acos(-1))dphi = fabs(2*acos(-1) - dphi);
00176     if(dphi > mDPhi && deta < mDEta && dj->M() > 10. && dijetFourReco[0]->Pt() > mPtHigh && dijetFourReco[1]->Pt() > mPtLow)yesdijet += 2;
00177     delete dj;
00178     dj = 0;
00179   }
00180   cout << "\n" <<"agrdl: event "<<(*nEvents)<<" "<<nJets<<" "<<nTrJets<<" "<<yesdijet<<endl;
00181   for(unsigned int k = 0; k < finalparticles.size(); k++)
00182     delete finalparticles[k];
00183   finalparticles.clear();
00184  
00185   if(yesdijet == 0) return 1;
00186 
00187   return 0;
00188 }
00189 //_______________________________________________________________
00190 int StDijetFilter::RejectGE(const StGenParticleMaster &ptl) const
00191 {
00192   return 0;
00193 }
00194 //_______________________________________________________________
00195 void StDijetFilter::readConfig()
00196 {
00197   ifstream cfile("dijet.cnf");
00198   while(1){
00199     string attr;
00200     float val;
00201     cfile >> attr >> val;
00202     if(!cfile.good())break;
00203     parseConfig(attr,val);
00204   }
00205 }
00206 //_______________________________________________________________
00207 void StDijetFilter::parseConfig(string attr, float val)
00208 {
00209   if(attr == "RBTOW"){
00210     mRBTOW = val;
00211     cout<<"RBTOW changed to "<<val<<endl;
00212     return;
00213   }
00214   if(attr == "RCONE"){
00215     mRcone = val;
00216     cout<<"RCONE changed to "<<val<<endl;
00217     return;
00218   }
00219   if(attr == "SEED"){
00220     mSeed = val;
00221     cout<<"SEED changed to "<<val<<endl;
00222     return;
00223   }
00224   if(attr == "ASSOC"){
00225     mAssoc = val;
00226     cout<<"ASSOC changed to "<<val<<endl;
00227     return;
00228   }
00229   if(attr == "SPLIT"){
00230     mSplitfraction = val;
00231     cout<<"SPLIT changed to "<<val<<endl;
00232     return;
00233   }
00234   if(attr == "MIDPOINT"){
00235     if(val == 1.0){
00236       mAddmidpoints = true;
00237       cout<<"MIDPOINT changed to TRUE"<<endl;
00238       return;
00239     }else if(val == 0.0){
00240       mAddmidpoints = false;
00241       cout<<"MIDPOINT changed to FALSE"<<endl;
00242       return;
00243     }else{
00244       cout<<"IMPRORER INPUT TO MIDPOINTS"<<endl;
00245     }
00246   }
00247   if(attr == "STABLEMIDPOINTS"){
00248       if(val == 1.0){
00249         mStablemidpoints = true;
00250         cout<<"STABLEMIDPOINTS changed to TRUE"<<endl;
00251         return;
00252       }else if(val == 0.0){
00253         mStablemidpoints = false;
00254         cout<<"STABLEMIDPOINTS changed to FALSE"<<endl;
00255         return;
00256       }else{
00257         cout<<"IMPRORER INPUT TO STABLEMIDPOINTS"<<endl;
00258       }
00259   }
00260   if(attr == "SPLITMERGE"){
00261     if(val == 1.0){
00262       mSplitmerge = true;
00263       cout<<"SPLITMERGE changed to TRUE"<<endl;
00264       return;
00265     }else if(val == 0.0){
00266       mSplitmerge = false;
00267       cout<<"SPLITMERGE changed to FALSE"<<endl;
00268       return;
00269     }else{
00270       cout<<"IMPRORER INPUT TO SPLITMERGE"<<endl;
00271     }
00272   }
00273   if(attr == "PARTICLEETA"){
00274     mParticleEtaRange = val;
00275     cout<<"PARTICLEETA changed to "<<val<<endl;
00276     return;
00277   }
00278   if(attr == "JETETAHIGH"){
00279     mJetEtaHigh = val;
00280     cout<<"JETETAHIGH changed to "<<val<<endl;
00281     return;
00282   }
00283   if(attr == "JETETALOW"){
00284     mJetEtaLow = val;
00285     cout<<"JETETALOW changed to "<<val<<endl;
00286     return;
00287   }
00288   if(attr == "DIJETPTLOW"){
00289     mPtLow = val;
00290     cout<<"DIJETPTLOW changed to "<<val<<endl;
00291     return;
00292   }
00293   if(attr == "DIJETPTHIGH"){
00294     mPtHigh = val;
00295     cout<<"DIJETPTHIGH changed to "<<val<<endl;
00296     return;
00297   }
00298   if(attr == "DPHI"){
00299     mDPhi = val;
00300     cout<<"DPHI changed to "<<val<<endl;
00301     return;
00302   }
00303   if(attr == "MINJETPT"){
00304     mMinJetPt = val;
00305     cout<<"MINJETPT changed to "<<val<<endl;
00306     return;
00307   }
00308   if(attr == "DETA"){
00309     mDEta = val;
00310     cout<<"DETA changed to "<<val<<endl;
00311     return;
00312   }
00313   if(attr == "RECOHADRON"){
00314     mRecohadron = val;
00315     cout<<"RECOHADRON changed to "<<val<<endl;
00316     return;
00317   }
00318   if(attr == "RECOLEPTON"){
00319     mRecolepton = val;
00320     cout<<"RECOLEPTON changed to "<<val<<endl;
00321     return;
00322   }
00323 }
00324 //_______________________________________________________________
00325 JetFourVec* StDijetFilter::recoJet(vector<JetFourVec*> v, double* vert) const
00326 {
00327   JetFourVec* j = new JetFourVec;
00328   for(unsigned int i = 0; i < v.size(); i++){
00329     int id = v[i]->getCode();
00330     if(abs(id) != 11 && id != -2212 && id != -2122){
00331       v[i]->setPtEtaPhiM(mRecohadron*v[i]->Pt(),v[i]->Eta(),v[i]->Phi(),v[i]->M());
00332     }else{
00333       v[i]->setPtEtaPhiM(mRecolepton*v[i]->Pt(),v[i]->Eta(),v[i]->Phi(),v[i]->M());
00334     }
00335     (*j) = (*j) + (*v[i]);
00336   }
00337 
00338   return j;
00339 }
00340 //_______________________________________________________________
00341 vector<JetFourVec*> StDijetFilter::addMidPoints(vector<JetFourVec*> seeds) const
00342 {
00343   unsigned int n = seeds.size();
00344   float pi = acos(-1.0);
00345   for(unsigned int k = 0; k < n; k++){
00346     for(unsigned int l = 0; l < k; l++){
00347       if(dR(seeds[k],seeds[l]) > 2*mRcone)continue;
00348       JetFourVec* s = new JetFourVec;
00349       float maxphi = max(seeds[k]->Phi(),seeds[l]->Phi());
00350       float minphi = min(seeds[k]->Phi(),seeds[l]->Phi());
00351       float dphi = fabs(maxphi - minphi);
00352       if(dphi > pi)minphi += 2*pi;
00353       float neta = 0.5*(seeds[k]->Eta() + seeds[l]->Eta());
00354       float nphi = 0.5*(maxphi + minphi);
00355       if(nphi > pi)nphi -= 2*pi;
00356       if(nphi < -pi) nphi += 2*pi;
00357       s->setPtEtaPhiM(0.001,neta,nphi,0);
00358       seeds.push_back(s);
00359     }
00360   }
00361   return seeds;
00362 }
00363 //_______________________________________________________________
00364 void StDijetFilter::split(vector<JetFourVec*> &v1, vector<JetFourVec*> &v2) const
00365 {
00366   JetFourVec* j1 = combineTracks(v1);
00367   JetFourVec* j2 = combineTracks(v2);
00368 
00369   vector<JetFourVec*>::iterator dit;
00370   for(vector<JetFourVec*>::iterator it = v1.begin(); it != v1.end(); it++){
00371     dit = find(v2.begin(),v2.end(),*it);
00372     if(dit == v2.end())continue;
00373     float d1 = dR(j1,*it);
00374     float d2 = dR(j2,*it);
00375     if(d2 > d1){
00376       v2.erase(dit);
00377     }
00378   }
00379   for(vector<JetFourVec*>::iterator it = v2.begin(); it != v2.end(); it++){
00380     dit = find(v1.begin(),v1.end(),*it);
00381     if(dit == v1.end())continue;
00382     float d1 = dR(j1,*it);
00383     float d2 = dR(j2,*it);
00384     if(d1 > d2){
00385       v1.erase(dit);
00386     }
00387   }
00388 
00389   delete j1;
00390   delete j2;
00391 }
00392 //_______________________________________________________________
00393 vector<JetFourVec*> StDijetFilter::merge(vector<JetFourVec*> v1, vector<JetFourVec*> v2) const
00394 {
00395   vector<JetFourVec*> nj = v1;
00396   vector<JetFourVec*>::iterator iter;
00397   for(iter = v2.begin(); iter != v2.end(); iter++){
00398     vector<JetFourVec*>::iterator it2 = find(v1.begin(),v1.end(),*iter);
00399     if(it2 == v1.end()) nj.push_back(*iter);
00400   }
00401 
00402   return nj;
00403 }
00404 //_______________________________________________________________
00405 vector< vector<JetFourVec*> > StDijetFilter::doSplitMerge(vector< vector<JetFourVec*> >orig) const
00406 {
00407   vector< vector<JetFourVec*> > njf;
00408   vector< vector<JetFourVec* > > jetFour = orig;
00409   int nChange = 0;
00410   int nIter = 0;
00411   while(1){
00412     if(nIter > 200)break;
00413     nChange = 0;
00414     njf.clear();
00415     njf = jetFour;
00416     vector< vector<JetFourVec*> >::iterator iter1;
00417     vector< vector<JetFourVec*> >::iterator iter2;
00418     for(iter1 = njf.begin(); iter1 != njf.end(); iter1++){
00419       if(nChange)continue;
00420       if((*iter1).size() == 0)continue;
00421       JetFourVec* j = combineTracks(*iter1);
00422       for(iter2 = iter1+1; iter2 != njf.end(); iter2++){
00423         if(nChange)continue;
00424         if((*iter2).size() == 0)continue;
00425         JetFourVec* nj = combineTracks(*iter2);
00426         if(dR(nj,j) > 2*mRcone){
00427           delete nj;
00428           continue;
00429         }
00430         float oe = overlapEnergy(*iter1,*iter2);
00431         if(oe == 0){
00432           delete nj;
00433           continue;
00434         }
00435         nChange = 1;
00436         vector< vector<JetFourVec*> >::iterator njit1;
00437         vector< vector<JetFourVec*> >::iterator njit2;
00438         njit1 = find(jetFour.begin(),jetFour.end(),*iter1);
00439         njit2 = find(jetFour.begin(),jetFour.end(),*iter2);
00440         if(oe > mSplitfraction){
00441           vector<JetFourVec*> mj = merge(*iter1,*iter2);
00442           (*njit1).clear();
00443           (*njit2).clear();
00444           jetFour.erase(njit2);
00445           jetFour.erase(njit1);
00446           jetFour.insert(jetFour.begin(),mj);
00447           delete nj;
00448           continue;
00449         }else{
00450           split(*njit1,*njit2);
00451           delete nj;
00452           continue;
00453         }
00454       }
00455       delete j;
00456     }
00457     jetFour = EtOrderedList(jetFour);
00458     if(!nChange)break;
00459     nIter++;
00460   }
00461   return njf;
00462 
00463 }
00464 //_______________________________________________________________
00465 float StDijetFilter::overlapEnergy(vector<JetFourVec*> v1,vector<JetFourVec*> v2) const
00466 {
00467   float e = 0;
00468   JetFourVec* j1 = combineTracks(v1);
00469   JetFourVec* j2 = combineTracks(v2);
00470   float high = min(j1->getEn()*sin(j1->Theta()),j2->getEn()*sin(j2->Theta()));
00471 
00472   for(vector<JetFourVec*>::iterator iter1 = v1.begin(); iter1 != v1.end(); iter1++){
00473     vector<JetFourVec*>::iterator iter2 = find(v2.begin(),v2.end(),*iter1);
00474     if(iter2 == v2.end())continue;
00475     e += (*iter2)->getEn()*sin((*iter2)->Theta());
00476   }
00477 
00478   delete j1;
00479   delete j2;
00480 
00481   return e/high;
00482 
00483 }
00484 //_______________________________________________________________
00485 vector<JetFourVec*> StDijetFilter::EtOrderedList(vector<JetFourVec*> jetFour) const
00486 {
00487   vector<JetFourVec*> njf;
00488   vector<JetFourVec*>::iterator jfvIter;
00489   vector<JetFourVec*>::iterator njfvIter;
00490   for(jfvIter = jetFour.begin(); jfvIter != jetFour.end(); jfvIter++){
00491     if((*jfvIter)->Pt() == 0)continue;
00492     int added = 0;
00493     for(njfvIter = njf.begin(); njfvIter != njf.end(); njfvIter++){
00494       if((*jfvIter)->Pt() > (*njfvIter)->Pt()){
00495         njf.insert(njfvIter,*jfvIter);
00496         added = 1;
00497         break;
00498       }
00499     }
00500     if(!added){
00501       njf.push_back(*jfvIter);
00502     }
00503   }
00504 
00505   return njf;
00506 }
00507 //_______________________________________________________________
00508 vector< vector<JetFourVec*> > StDijetFilter::RemoveDuplicates(vector< vector<JetFourVec*> >jetFour) const
00509 {
00510   vector< vector<JetFourVec*> > njf;
00511   vector< vector<JetFourVec*> >::iterator jfvIter;
00512   vector< vector<JetFourVec*> >::iterator njfvIter;
00513   for(jfvIter = jetFour.begin(); jfvIter != jetFour.end(); jfvIter++){
00514     if((*jfvIter).size() == 0)continue;
00515     JetFourVec* j = combineTracks(*jfvIter);
00516     int dupe = 0;
00517     for(njfvIter = njf.begin(); njfvIter != njf.end(); njfvIter++){
00518       JetFourVec* nj = combineTracks(*njfvIter);
00519       if(*nj == *j)dupe = 1;
00520       delete nj;
00521     }
00522     if(!dupe)njf.push_back(*jfvIter);
00523     delete j;
00524   }
00525   return njf;
00526 }
00527 //_______________________________________________________________
00528 vector< vector<JetFourVec*> > StDijetFilter::EtOrderedList(vector< vector<JetFourVec*> >jetFour) const
00529 {
00530   vector< vector<JetFourVec*> > njf;
00531   vector< vector<JetFourVec*> >::iterator jfvIter;
00532   vector< vector<JetFourVec*> >::iterator njfvIter;
00533   for(jfvIter = jetFour.begin(); jfvIter != jetFour.end(); jfvIter++){
00534     if((*jfvIter).size() == 0)continue;
00535     JetFourVec* j = combineTracks(*jfvIter);
00536     if(j->Pt() == 0){ delete j;continue;}
00537     int added = 0;
00538     for(njfvIter = njf.begin(); njfvIter != njf.end(); njfvIter++){
00539       JetFourVec* nj = combineTracks(*njfvIter);
00540       if(j->Pt() > nj->Pt()){
00541         njf.insert(njfvIter,*jfvIter);
00542         delete nj;
00543         added = 1;
00544         break;
00545       }
00546       delete nj;
00547     }
00548     if(!added){
00549       njf.push_back(*jfvIter);
00550     }
00551     delete j;
00552   }
00553 
00554   return njf;
00555 
00556 }
00557 //_______________________________________________________________
00558 float StDijetFilter::dR(StGenParticle* p1,StGenParticle* p2) const
00559 {
00560   float dphi = fabs(p1->Phi() - p2->Phi());
00561   float deta = p1->Eta() - p2->Eta();
00562   float pi = acos(-1.0);
00563   if(dphi > pi)dphi = 2*pi - dphi;
00564   return sqrt(pow(deta,2) + pow(dphi,2));
00565 }
00566 //_______________________________________________________________
00567 float StDijetFilter::dR(StGenParticle* p,JetFourVec* v) const
00568 {
00569   float dphi = fabs(p->Phi() - v->Phi());
00570   float deta = p->Eta() - v->Eta();
00571   float pi = acos(-1.0);
00572   if(dphi > pi)dphi = 2*pi - dphi;
00573   return sqrt(pow(deta,2) + pow(dphi,2));
00574 }
00575 //_______________________________________________________________
00576 float StDijetFilter::dR(JetFourVec* v,StGenParticle* p) const
00577 {
00578   return dR(p,v);
00579 }
00580 //_______________________________________________________________
00581 float StDijetFilter::dR(JetFourVec* v1, JetFourVec* v2) const
00582 {
00583   float dphi = fabs(v1->Phi() - v2->Phi());
00584   float  deta = v1->Eta() - v2->Eta();
00585   float pi = acos(-1.0);
00586   if(dphi > pi)dphi = 2*pi - dphi;
00587   return sqrt(pow(deta,2) + pow(dphi,2));
00588 }
00589 //_______________________________________________________________
00590 JetFourVec* StDijetFilter::combineTracks(vector<JetFourVec*> jf) const
00591 {
00592   JetFourVec* j = new JetFourVec;
00593   if(jf.size() == 0)return j;
00594   for(unsigned int i = 0; i < jf.size(); i++){
00595     (*j) = (*j) + (*jf[i]);
00596   }
00597 
00598   return j;
00599 
00600 }
00601 //_______________________________________________________________
00602 JetFourVec::JetFourVec()
00603 {
00604   px = 0,py = 0,pz = 0,en = 0;code = 0;
00605 }
00606 //_______________________________________________________________
00607 JetFourVec::JetFourVec(JetFourVec* v)
00608 {
00609   code = v->getCode();
00610   px = v->getPx();
00611   py = v->getPy();
00612   pz = v->getPz();
00613   en = v->getEn();
00614 }
00615 //_______________________________________________________________
00616 JetFourVec::JetFourVec(StGenParticle* g)
00617 {
00618   double p[4];
00619   g->Momentum(p);
00620   px = p[0];
00621   py = p[1];
00622   pz = p[2];
00623   en = p[3];
00624   code = g->GetPdgCode();
00625 }
00626 //_______________________________________________________________
00627 void JetFourVec::setPxPyPzEn(float x, float y, float z, float e)
00628 {
00629   setPx(x);
00630   setPy(y);
00631   setPz(z);
00632   setEn(e);
00633 }
00634 //_______________________________________________________________
00635 float JetFourVec::P()
00636 {
00637   return sqrt(px*px + py*py + pz*pz);
00638 }
00639 //_______________________________________________________________
00640 float JetFourVec::Pt()
00641 {
00642   return sqrt(px*px + py*py);
00643 }
00644 //_______________________________________________________________
00645 float JetFourVec::Eta()
00646 {
00647   float pmom = P();
00648   if(pmom > fabs(pz))return -log(tan(Theta()/2.));
00649   else return 1.e30;
00650 }
00651 //_______________________________________________________________
00652 float JetFourVec::Phi()
00653 {
00654   return atan2(py,px);
00655 }
00656 //_______________________________________________________________
00657 float JetFourVec::Theta()
00658 {
00659   return acos(pz/P());
00660 }
00661 //_______________________________________________________________
00662 float JetFourVec::M()
00663 {
00664   if(en*en - P()*P() < 0)return 0;
00665   return sqrt(en*en-P()*P());
00666 }
00667 //_______________________________________________________________
00668 void JetFourVec::setPtEtaPhiM(float pt, float eta, float phi, float m)
00669 {
00670   float x = pt * cos(phi);
00671   float y = pt * sin(phi);
00672   float p = pt*cosh(eta);
00673   float z = sqrt(p*p-pt*pt)*eta/fabs(eta);
00674   float e = sqrt(m*m + p*p);
00675   setPxPyPzEn(x,y,z,e);
00676 }
00677 //_______________________________________________________________
00678 void JetFourVec::setPtEtaPhiE(float pt, float eta, float phi, float e)
00679 {
00680   float x = pt * cos(phi);
00681   float y = pt * sin(phi);
00682   float p = pt*cosh(eta);
00683   float z = sqrt(p*p-pt*pt)*eta/fabs(eta);
00684   setPxPyPzEn(x,y,z,e);
00685 }
00686 //_______________________________________________________________
00687 JetFourVec JetFourVec::operator +(JetFourVec x)
00688 {
00689   JetFourVec t;
00690   t.setCode(0);
00691   t.setPx(x.getPx() + px);
00692   t.setPy(x.getPy() + py);
00693   t.setPz(x.getPz() + pz);
00694   t.setEn(x.getEn() + en);
00695 
00696   return t;
00697 }
00698 //_______________________________________________________________
00699 bool JetFourVec::operator ==(JetFourVec x)
00700 {
00701   if(code != x.getCode())return false;
00702   if(fabs(px - x.getPx()) > 1e-4)return false;
00703   if(fabs(py - x.getPy()) > 1e-4)return false;
00704   if(fabs(pz - x.getPz()) > 1e-4)return false;
00705   if(fabs(en - x.getEn()) > 1e-4)return false;
00706 
00707   return true;
00708 }