00001 #include <cmath>
00002 #include <TArrayD.h>
00003 #include "StEmcMath.h"
00004 #include "StEvent/StMeasuredPoint.h"
00005 #include "StarClassLibrary/StThreeVectorF.hh"
00006 #include "StEmcUtil/others/emcInternalDef.h"
00007 #include "StEmcUtil/geometry/StEmcGeom.h"
00008
00009 ClassImp(StEmcMath)
00010
00011 Bool_t
00012 StEmcMath::etaPhi(StMeasuredPoint* point,StMeasuredPoint* vertex, Double_t &eta, Double_t &phi)
00013 {
00014
00015
00016
00017 const StThreeVectorF *p1, *p2;
00018 p1 = NULL;
00019 StThreeVectorF vt;
00020
00021 if(point) {
00022 if(vertex) p1 = &vertex->position();
00023 p2 = &point->position();
00024
00025 if(p1) vt = (*p2) - (*p1);
00026 else vt = (*p2);
00027
00028 eta = vt.pseudoRapidity();
00029 phi = vt.phi();
00030 return kTRUE;
00031 }
00032 return kFALSE;
00033 }
00034
00035 Double_t
00036 StEmcMath::pseudoRapidity(StMeasuredPoint* point,StMeasuredPoint* vertex)
00037 {
00038 Double_t eta, phi;
00039 if(etaPhi(point,vertex, eta,phi)) return eta;
00040 else return -9999.;
00041 }
00042
00043 Double_t
00044 StEmcMath::phi(StMeasuredPoint* point,StMeasuredPoint* vertex)
00045 {
00046 Double_t eta, phi;
00047 if(etaPhi(point,vertex, eta,phi)) return phi;
00048 else return -9999.;
00049 }
00050
00051 UInt_t
00052 StEmcMath::detectorId(const StDetectorId stId)
00053 {
00054
00055 Int_t id = stId - kBarrelEmcTowerIdentifier + 1;
00056 if(id<BEMC || id> ESMDP) return 0;
00057 else return UInt_t(id);
00058 }
00059
00060 StDetectorId
00061 StEmcMath::detectorId(const UInt_t id)
00062 {
00063
00064 StDetectorId stId = StDetectorId(id + kBarrelEmcTowerIdentifier - 1);
00065 if(stId<kBarrelEmcTowerIdentifier || stId >kEndcapSmdVStripIdentifier)
00066 return kUnknownId;
00067 else return stId;
00068 }
00069
00070 Double_t
00071 StEmcMath::getPhiPlusMinusPi(const Double_t phi)
00072 {
00073
00074 Double_t phiw = phi;
00075 while(phiw >= M_PI) phiw -= 2.*M_PI;
00076 while(phiw < -M_PI) phiw += 2.*M_PI;
00077 return phiw;
00078 }
00079
00080 TArrayD *StEmcMath::binForSmde(Bool_t kprint)
00081 {
00082
00083 StEmcGeom* geom = StEmcGeom::getEmcGeom(3);
00084 const Int_t neta = geom->NEta();
00085 Int_t iw1, iw2;
00086 const Float_t *eb = geom->Eta();
00087 TArrayD *xb = new TArrayD(2*neta+1);
00088 (*xb)[neta] = 0.0;
00089 for(Int_t ik=0; ik<neta; ik++){
00090 iw1 = neta + 1 + ik;
00091 iw2 = neta-ik-1;
00092 Float_t x1 = eb[ik], x2, xw;
00093 if(ik<neta-1) {
00094 x2 = eb[ik+1];
00095 xw = (x1+x2)*0.5;
00096 }
00097 else xw = 0.99;
00098 (*xb)[iw1] = +xw;
00099 (*xb)[iw2] = -xw;
00100 if(kprint)
00101 printf(" iw1 %i %f => iw2 %i %f => eta %f\n", iw1,(*xb)[iw1], iw2,(*xb)[iw2], eb[ik]);
00102 }
00103 return xb;
00104 }