StRoot: StSpinPool/StJetEvent/StJetCandidate.cxx Source File

00001 //
00002 // Pibero Djawotho <pibero@tamu.edu>
00003 // Texas A&M University
00004 // 3 September 2009
00005 //
00006 
00007 #include "StJetTrack.h"
00008 #include "StJetTower.h"
00009 #include "StJetCandidate.h"
00010 
00011 ClassImp(StJetCandidate);
00012 
00013 StJetCandidate::StJetCandidate(const TVector3& vertex, const TLorentzVector& fourMomentum)
00014   : mPt(fourMomentum.Pt())
00015   , mEta(fourMomentum.Eta())
00016   , mPhi(fourMomentum.Phi())
00017   , mE(fourMomentum.E())
00018 {
00019   mDetEta = detEta(vertex);
00020 }
00021 
00022 float StJetCandidate::detEta(const TVector3& vertex) const
00023 {
00024   float DetEta;
00025   if (getBarrelDetectorEta(vertex,DetEta)) return DetEta;
00026   if (getEndcapDetectorEta(vertex,DetEta)) return DetEta;
00027   return -999;
00028 }
00029 
00030 bool StJetCandidate::getBarrelDetectorEta(const TVector3& vertex, float& detEta) const
00031 {
00032   // Front plate of BEMC or BPRS layer
00033   // See StEmcGeom/geometry/StEmcGeom.cxx
00034   static const double BEMC_RADIUS = 225.405; // cm
00035   TVector3 pos = momentum();
00036   pos.SetMag(BEMC_RADIUS/pos.Unit().Perp());
00037   pos += vertex;
00038   detEta = pos.Eta();
00039   return fabs(detEta) < 1;
00040 }
00041 
00042 bool StJetCandidate::getEndcapDetectorEta(const TVector3& vertex, float& detEta) const
00043 {
00044   // EEMC preshower1 layer
00045   // See StEEmcUtil/EEmcGeom/EEmcGeomDefs.h
00046   static const double EEMC_Z = 270.190; // cm
00047   TVector3 pos = momentum();
00048   if (pos.z() > 0) {
00049     // Real endcap
00050     pos.SetMag((EEMC_Z-vertex.z())/pos.Unit().z());
00051     pos += vertex;
00052     detEta = pos.Eta();
00053   }
00054   else {
00055     // Mirror endcap (flip about xy-plane, solve and flip sign)
00056     TVector3 pos2(pos.x(),pos.y(),-pos.z());
00057     TVector3 vertex2(vertex.x(),vertex.y(),-vertex.z());
00058     pos2.SetMag((EEMC_Z-vertex2.z())/pos2.Unit().z());
00059     pos2 += vertex2;
00060     detEta = -pos2.Eta();
00061   }
00062   return fabs(detEta) > 1 && fabs(detEta) < 2;
00063 }
00064 
00065 StJetTrack* StJetCandidate::leadingChargedParticle() const
00066 {
00067   StJetTrack* lcp = 0;
00068   for (int i = 0; i < numberOfTracks(); ++i) {
00069     StJetTrack* t = track(i);
00070     if (!lcp || t->pt() > lcp->pt()) lcp = t;
00071   }
00072   return lcp;
00073 }
00074 
00075 float StJetCandidate::deltaR(const StJetElement* element) const
00076 {
00077   return momentum().DeltaR(element->momentum());
00078 }
00079 
00080 float StJetCandidate::sumTrackPt() const
00081 {
00082   float s = 0;
00083   for (int i = 0; i < numberOfTracks(); ++i) s += track(i)->pt();
00084   return s;
00085 }
00086 
00087 float StJetCandidate::sumTrackPt(float radius) const
00088 {
00089   float s = 0;
00090   for (int i = 0; i < numberOfTracks(); ++i) {
00091     const StJetTrack* t = track(i);
00092     if (deltaR(t) < radius) s += t->pt();
00093   }
00094   return s;
00095 }
00096 
00097 float StJetCandidate::sumTowerPt() const
00098 {
00099   float s = 0;
00100   for (int i = 0; i < numberOfTowers(); ++i) s += tower(i)->pt();
00101   return s;
00102 }
00103 
00104 float StJetCandidate::sumTowerPt(float radius) const
00105 {
00106   float s = 0;
00107   for (int i = 0; i < numberOfTowers(); ++i) {
00108     const StJetTower* t = tower(i);
00109     if (deltaR(t) < radius) s += t->pt();
00110   }
00111   return s;
00112 }
00113 
00114 StJetTrack* StJetCandidate::getTrackById(int id) const
00115 {
00116   for (int i = 0; i < numberOfTracks(); ++i) {
00117     StJetTrack* t = track(i);
00118     if (t->id() == id) return t;
00119   }
00120   return 0;
00121 }
00122 
00123 StJetTower* StJetCandidate::getTowerById(int id) const
00124 {
00125   for (int i = 0; i < numberOfTowers(); ++i) {
00126     StJetTower* t = tower(i);
00127     if (t->id() == id) return t;
00128   }
00129   return 0;
00130 }
00131 
00132 float StJetCandidate::getJetPatchPhi(int jetPatch)
00133 {
00134   return TVector2::Phi_mpi_pi((150 - (jetPatch % 6) * 60) * TMath::DegToRad());
00135 }
00136 
00137 bool StJetCandidate::getBarrelJetPatchEtaPhi(int jetPatch, float& eta, float& phi)
00138 {
00139   //
00140   // Pibero Djawotho <pibero@tamu.edu>
00141   // Texas A&M University
00142   // 5 September 2009
00143   // Revised 9 Feruary 2010 to include middle jet patches
00144   //
00145 
00146   // Sanity check
00147   if (jetPatch < 0 || jetPatch >= 18) return false;
00148 
00149   //
00150   // The jet patches are numbered starting with JP0 centered at 150 degrees
00151   // looking from the West into the IR (intersection region) and increasing
00152   // clockwise, i.e. JP1 at 90 degrees, JP2 at 30 degrees, etc. On the East
00153   // side the numbering picks up at JP6 centered again at 150 degrees and
00154   // increasing clockwise (again as seen from the *West* into the IR). Thus
00155   // JP0 and JP6 are in the same phi location in the STAR coordinate system.
00156   // So are JP1 and JP7, etc.
00157   //
00158   // Jet Patch#  Eta  Phi  Quadrant
00159   //
00160   //          0  0.5  150       10'
00161   //          1  0.5   90       12'
00162   //          2  0.5   30        2'
00163   //          3  0.5  -30        4'
00164   //          4  0.5  -90        6'
00165   //          5  0.5 -150        8'
00166   //          6 -0.5  150       10'
00167   //          7 -0.5   90       12'
00168   //          8 -0.5   30        2'
00169   //          9 -0.5  -30        4'
00170   //         10 -0.5  -90        6'
00171   //         11 -0.5 -150        8'
00172   //
00173   // http://www.nikhef.nl/~ogrebeny/emc/files/Towers%20Layout.pdf
00174   // http://www.nikhef.nl/~ogrebeny/emc/files/BEMC.pdf
00175   // http://drupal.star.bnl.gov/STAR/system/files/BEMC_y2004.pdf
00176   //
00177 
00178   if (jetPatch >=  0 && jetPatch <  6) eta =  0.5;
00179   if (jetPatch >=  6 && jetPatch < 12) eta = -0.5;
00180   if (jetPatch >= 12 && jetPatch < 18) eta = -0.1;
00181 
00182   phi = getJetPatchPhi(jetPatch);
00183 
00184   return true;
00185 }
00186 
00187 bool StJetCandidate::getEndcapJetPatchEtaPhi(int jetPatch, float& eta, float& phi)
00188 {
00189   if (jetPatch >= 0 && jetPatch < 6)
00190     eta = 1.5;
00191   else
00192     return false;
00193 
00194   phi = getJetPatchPhi(jetPatch);
00195 
00196   return true;
00197 }
00198 
00199 bool StJetCandidate::getOverlapJetPatchEtaPhi(int jetPatch, float& eta, float& phi)
00200 {
00201   if (jetPatch >= 0 && jetPatch < 6)
00202     eta = 0.9;
00203   else
00204     return false;
00205 
00206   phi = getJetPatchPhi(jetPatch);
00207 
00208   return true;
00209 }
00210 
00211 bool StJetCandidate::getBarrelJetPatchId(float eta, float phi, int& id)
00212 {
00213   //
00214   // Pibero Djawotho <pibero@tamu.edu>
00215   // Texas A&M University
00216   // 5 September 2009
00217   // To do: Need to include possibility of returning middle jet patch
00218   //
00219 
00220   // Jet patch id is left at -1 on failure
00221   id = -1;
00222 
00223   // Check range of eta
00224   if (eta < -1 || eta > 1) return false;
00225 
00226   // Map phi into [-pi,pi] if necessary
00227   if (phi < -M_PI || phi > M_PI) phi = TVector2::Phi_mpi_pi(phi);
00228 
00229   // Get jet patch id
00230   static const double PI_OVER_3 = M_PI/3;
00231 
00232   if (0 <= eta && eta <= 1) {
00233     if ( 2*PI_OVER_3 <= phi && phi <         M_PI) id = 0;
00234     if (   PI_OVER_3 <= phi && phi <  2*PI_OVER_3) id = 1;
00235     if (           0 <= phi && phi <    PI_OVER_3) id = 2;
00236     if (  -PI_OVER_3 <= phi && phi <            0) id = 3;
00237     if (-2*PI_OVER_3 <= phi && phi <   -PI_OVER_3) id = 4;
00238     if (       -M_PI <= phi && phi < -2*PI_OVER_3) id = 5;
00239   }
00240 
00241   if (-1 <= eta && eta < 0) {
00242     if ( 2*PI_OVER_3 <= phi && phi <         M_PI) id = 6;
00243     if (   PI_OVER_3 <= phi && phi <  2*PI_OVER_3) id = 7;
00244     if (           0 <= phi && phi <    PI_OVER_3) id = 8;
00245     if (  -PI_OVER_3 <= phi && phi <            0) id = 9;
00246     if (-2*PI_OVER_3 <= phi && phi <   -PI_OVER_3) id = 10;
00247     if (       -M_PI <= phi && phi < -2*PI_OVER_3) id = 11;
00248   }
00249 
00250   return (0 <= id && id < 12);
00251 }
00252 
00253 bool StJetCandidate::getEndcapJetPatchId(float eta, float phi, int& id)
00254 {
00255   // Jet patch id is left at -1 on failure
00256   id = -1;
00257 
00258   // Check range of eta
00259   if (eta < 1 || eta > 2) return false;
00260 
00261   // Map phi into [-pi,pi] if necessary
00262   if (phi < -M_PI || phi > M_PI) phi = TVector2::Phi_mpi_pi(phi);
00263 
00264   // Get jet patch id
00265   static const double PI_OVER_3 = M_PI/3;
00266 
00267   if (1 <= eta && eta <= 2) {
00268     if ( 2*PI_OVER_3 <= phi && phi <         M_PI) id = 0;
00269     if (   PI_OVER_3 <= phi && phi <  2*PI_OVER_3) id = 1;
00270     if (           0 <= phi && phi <    PI_OVER_3) id = 2;
00271     if (  -PI_OVER_3 <= phi && phi <            0) id = 3;
00272     if (-2*PI_OVER_3 <= phi && phi <   -PI_OVER_3) id = 4;
00273     if (       -M_PI <= phi && phi < -2*PI_OVER_3) id = 5;
00274   }
00275 
00276   return (0 <= id && id < 6);
00277 }
00278 
00279 bool StJetCandidate::getOverlapJetPatchId(float eta, float phi, int& id)
00280 {
00281   // Jet patch id is left at -1 on failure
00282   id = -1;
00283 
00284   // Check range of eta
00285   if (eta < 0.4 || eta > 1.4) return false;
00286 
00287   // Map phi into [-pi,pi] if necessary
00288   if (phi < -M_PI || phi > M_PI) phi = TVector2::Phi_mpi_pi(phi);
00289 
00290   // Get jet patch id
00291   static const double PI_OVER_3 = M_PI/3;
00292 
00293   if (0.4 <= eta && eta <= 1.4) {
00294     if ( 2*PI_OVER_3 <= phi && phi <         M_PI) id = 0;
00295     if (   PI_OVER_3 <= phi && phi <  2*PI_OVER_3) id = 1;
00296     if (           0 <= phi && phi <    PI_OVER_3) id = 2;
00297     if (  -PI_OVER_3 <= phi && phi <            0) id = 3;
00298     if (-2*PI_OVER_3 <= phi && phi <   -PI_OVER_3) id = 4;
00299     if (       -M_PI <= phi && phi < -2*PI_OVER_3) id = 5;
00300   }
00301 
00302   return (0 <= id && id < 6);
00303 }