11 #include "StEmcPosition.h"
13 #include "SystemOfUnits.h"
14 #include "PhysicalConstants.h"
15 #include "StPhysicalHelixD.hh"
17 #include "StMcEvent.hh"
18 #include "StMcEventTypes.hh"
19 #include "StMcEventMaker/StMcEventMaker.h"
21 #include "StEventTypes.h"
22 #include "StEmcUtil/geometry/StEmcGeom.h"
25 #include "StMuDSTMaker/COMMON/StMuTrack.h"
32 mGeom[0] = StEmcGeom::getEmcGeom(
"bemc");
33 mGeom[1] = StEmcGeom::getEmcGeom(
"bprs");
34 mGeom[2] = StEmcGeom::getEmcGeom(
"bsmde");
35 mGeom[3] = StEmcGeom::getEmcGeom(
"bsmdp");
38 StEmcPosition::~StEmcPosition()
43 const StMuTrack*
const track,
double magField,
double radius,
int option)
const
47 *momentumAtFinal=Zero;
58 double charge = track->
charge();
62 s1 = pathLength.first;
63 s2 = pathLength.second;
68 if (finite(s1) == 0 && finite(s2) == 0) {
return kFALSE;}
73 if (s1 >= 0 && s2 >= 0) {s = s1; goProj = kTRUE; }
74 if (s1 >= 0 && s2 < 0) { s = s1; goProj = kTRUE; }
75 if (s1 < 0 && s2 >= 0) { s = s2; goProj = kTRUE; }
81 if (s1 <= 0 && s2 <= 0) { s = s2; goProj = kTRUE; }
82 if (s1 <= 0 && s2 > 0) { s = s1; goProj = kTRUE; }
83 if (s1 > 0 && s2 <= 0) { s = s2; goProj = kTRUE; }
88 *atFinal = helix.at( s );
89 *momentumAtFinal = helix.momentumAt( s, magField*tesla );
90 if (charge == 0) *momentumAtFinal = momentum;
96 const StPhysicalHelixD*
const helix, Double_t magField, Double_t radius, Int_t option)
const
100 *momentumAtFinal=Zero;
106 s1 = pathLength.first;
107 s2 = pathLength.second;
112 if (finite(s1) == 0 && finite(s2) == 0) {
return kFALSE;}
117 if (s1 >= 0 && s2 >= 0) {s = s1; goProj = kTRUE; }
118 if (s1 >= 0 && s2 < 0) { s = s1; goProj = kTRUE; }
119 if (s1 < 0 && s2 >= 0) { s = s2; goProj = kTRUE; }
125 if (s1 <= 0 && s2 <= 0) { s = s2; goProj = kTRUE; }
126 if (s1 <= 0 && s2 > 0) { s = s1; goProj = kTRUE; }
127 if (s1 > 0 && s2 <= 0) { s = s2; goProj = kTRUE; }
132 *atFinal = helix->at( s );
133 *momentumAtFinal = helix->momentumAt( s, magField*tesla );
140 const StTrack*
const track, Double_t magField, Double_t radius, Int_t option)
const
144 *momentumAtFinal=Zero;
147 const StThreeVectorF& momentum = track->outerGeometry()->momentum();
148 Double_t charge = track->outerGeometry()->charge();
154 s1 = pathLength.first;
155 s2 = pathLength.second;
160 if (finite(s1) == 0 && finite(s2) == 0) {
return kFALSE;}
165 if (s1 >= 0 && s2 >= 0) {s = s1; goProj = kTRUE; }
166 if (s1 >= 0 && s2 < 0) { s = s1; goProj = kTRUE; }
167 if (s1 < 0 && s2 >= 0) { s = s2; goProj = kTRUE; }
173 if (s1 <= 0 && s2 <= 0) { s = s2; goProj = kTRUE; }
174 if (s1 <= 0 && s2 > 0) { s = s1; goProj = kTRUE; }
175 if (s1 > 0 && s2 <= 0) { s = s2; goProj = kTRUE; }
180 *atFinal = helix.at( s );
181 *momentumAtFinal = helix.momentumAt( s, magField*tesla );
182 if (charge == 0) *momentumAtFinal = momentum;
188 const StMcTrack*
const mcTrack, Double_t magField, Double_t radius, Int_t option)
const
192 *momentumAtFinal=Zero;
194 const StThreeVectorF& origin = mcTrack->startVertex()->position();
196 Double_t charge = mcTrack->particleDefinition()->charge();
202 s1 = pathLength.first;
203 s2 = pathLength.second;
208 if (finite(s1) == 0 && finite(s2) == 0) {
return kFALSE;}
213 if (s1 >= 0 && s2 >= 0) { s = s1; goProj = kTRUE; }
214 if (s1 >= 0 && s2 < 0) { s = s1; goProj = kTRUE; }
215 if (s1 < 0 && s2 >= 0) { s = s2; goProj = kTRUE; }
221 if (s1 <= 0 && s2 <= 0) { s = s2; goProj = kTRUE; }
222 if (s1 <= 0 && s2 > 0) { s = s1; goProj = kTRUE; }
223 if (s1 > 0 && s2 <= 0) { s = s2; goProj = kTRUE; }
228 *atFinal = helix.at( s );
229 *momentumAtFinal = helix.momentumAt( s, magField*tesla );
230 if (charge == 0) *momentumAtFinal = momentum;
247 float xO = origin.x();
248 float yO = origin.y();
249 float distToOrigin = ::sqrt( ::pow(xO, 2) + ::pow(yO, 2) );
250 if ( distToOrigin < emcRadius )
253 Bool_t projTrackOk =
projTrack( position, momentum, track, magField, emcRadius );
257 int m = 0, e = 0, s = 0;
258 float phi = position->phi();
259 float eta = position->pseudoRapidity();
263 if ( mGeom[0]->getBin(phi, eta, m, e, s) == 0 && s != -1 )
return kTRUE;
271 const StTrack*
const track, Double_t magField, Double_t emcRadius )
const
275 if (!track->outerGeometry())
return kFALSE;
278 Float_t xO = origin.x();
279 Float_t yO = origin.y();
280 Float_t distToOrigin = ::sqrt( ::pow(xO, 2) + ::pow(yO, 2) );
281 if ( distToOrigin < emcRadius )
283 Bool_t projTrackOk =
projTrack( position, momentum, track, magField, emcRadius );
286 Int_t m = 0, e = 0, s = 0;
287 Float_t phi = position->phi();
288 Float_t eta = position->pseudoRapidity();
289 if ( mGeom[0]->getBin(phi, eta, m, e, s) == 0 && s != -1 )
return kTRUE;
297 const StMcTrack*
const mcTrack, Double_t magField, Double_t emcRadius )
const
299 Float_t startVertexX = mcTrack->startVertex()->position().x();
300 Float_t startVertexY = mcTrack->startVertex()->position().y();
301 Float_t startVtxToOrigin = ::sqrt( ::pow( startVertexX, 2 ) + ::pow( startVertexY, 2 ) );
303 if ( !mcTrack->stopVertex() && startVtxToOrigin < emcRadius )
305 Bool_t projTrackOk =
projTrack( position, momentum, mcTrack, magField, emcRadius );
308 Int_t m = 0, e = 0, s = 0;
309 Float_t phi = position->phi();
310 Float_t eta = position->pseudoRapidity();
311 if ( mGeom[0]->getBin(phi, eta, m, e, s) == 0 && s != -1 )
return kTRUE;
316 Float_t stopVtxToOrigin = -1;
317 if ( mcTrack->stopVertex() )
319 Float_t stopVertexX = mcTrack->stopVertex()->position().x();
320 Float_t stopVertexY = mcTrack->stopVertex()->position().y();
321 stopVtxToOrigin = ::sqrt( ::pow( stopVertexX,2 ) + ::pow(stopVertexY,2) );
324 if (stopVtxToOrigin >= emcRadius)
326 Bool_t projTrackOk =
projTrack( position, momentum, mcTrack, magField, emcRadius );
329 Int_t m = 0, e = 0, s = 0;
330 Float_t phi = position->phi();
331 Float_t eta = position->pseudoRapidity();
332 if ( mGeom[0]->getBin(phi, eta, m, e, s) == 0 && s != -1 )
return kTRUE;
340 Float_t* towerEta, Float_t* towerPhi )
const
342 *towerEta = 0; *towerPhi = 0;
343 Float_t tempTowerEta = 0, tempTowerPhi = 0;
344 Int_t m = 0, e = 0, s = 0, towerId = -1;
346 mGeom[0]->
getBin(phi, eta, m, e, s);
349 mGeom[0]->getId(m, e, s, towerId);
350 mGeom[0]->getEtaPhi(towerId, tempTowerEta, tempTowerPhi);
351 *towerEta = tempTowerEta;
352 *towerPhi = tempTowerPhi;
359 mGeom[0]->
getBin( phi, eta, m, e, s );
370 if(softId<1 || softId>4800)
return 0;
372 mGeom[0]->
getBin(softId,m,e,s);
378 if(m<1 || m>120)
return 0;
379 if(e<1 || e>20)
return 0;
380 if(s<1 || s>2)
return 0;
381 return getNextId(1,m,e,s,nTowersdEta,nTowersdPhi);
384 Int_t
StEmcPosition::getNextId(
const Int_t det,
const Int_t m,
const Int_t e,
const Int_t s,
const Int_t nEta,
const Int_t nPhi)
const
386 if(det<1 || det>4)
return 0;
387 if(m<1 || m>120)
return 0;
388 if(s<1 || s>mGeom[det-1]->NSub())
return 0;
389 if(e<1 || e>mGeom[det-1]->NEta())
return 0;
395 Int_t NE=mGeom[det-1]->NEta();
396 Int_t NS=mGeom[det-1]->NSub();
398 if(abs(ef)>NE)
return 0;
416 }
while(sf<=0 || sf>NS);
425 mGeom[det-1]->getId(mf, ef, sf, rid);
426 mGeom[det-1]->getEtaPhi(rid, eta, phi);
427 mGeom[det-1]->
getBin(phi,-eta,mf,etmp,stmp);
431 if(mf<1 || mf>120)
return 0;
432 if(ef<1 || ef>NE)
return 0;
433 if(sf<1 || sf>NS)
return 0;
434 mGeom[det-1]->getId(mf, ef, sf, rid);
441 if(det<1 || det>4)
return 0;
443 if(softId<1)
return 0;
444 if((det == 1 || det == 2) && softId > 4800)
return 0;
445 if((det == 3 || det == 4) && softId > 18000)
return 0;
446 mGeom[det-1]->
getBin(softId,m,e,s);
456 Int_t nTowersdEta, Int_t nTowersdPhi )
const
460 Float_t towerEta = 0, towerToTrackdEta = 0;
461 Float_t towerPhi = 0, towerToTrackdPhi = 0;
462 Float_t mdistTowerToTrack = 0;
464 towerId =
getNextTowerId( trackEta, trackPhi, nTowersdEta, nTowersdPhi );
468 mGeom[0]->getEtaPhi(towerId, towerEta, towerPhi);
469 towerToTrackdEta = towerEta-trackEta;
470 towerToTrackdPhi = towerPhi-trackPhi;
472 mdistTowerToTrack = ::sqrt( ::pow(towerToTrackdEta, 2) + ::pow(towerToTrackdPhi, 2) );
474 return mdistTowerToTrack;
483 if(TowerId<1 || TowerId>4800)
return Zero;
485 Float_t xTower,yTower,zTower;
487 mGeom[0]->getXYZ(TowerId, xTower, yTower, zTower);
491 return PositionFromVertex;
497 if(TowerId<1 || TowerId>4800)
return Zero;
499 float xTower,yTower,zTower;
501 mGeom[0]->getXYZ(TowerId, xTower, yTower, zTower);
505 return PositionFromVertex;
511 if(TowerId<1 || TowerId>4800)
return Zero;
513 Float_t xTower,yTower,zTower;
515 mGeom[0]->getXYZ(TowerId, xTower, yTower, zTower);
519 return PositionFromVertex;
543 return p.pseudoRapidity();
549 return p.pseudoRapidity();
555 return p.pseudoRapidity();
Float_t getThetaFromVertex(const StVertex *const vertex, Int_t TowerId) const
Return theta of the tower considering the collision vertex.
Bool_t trackOnEmc(StThreeVectorD *position, StThreeVectorD *momentum, const StTrack *const track, Double_t magField, Double_t emcRadius=225.405) const
Track projection utility.
Float_t getDistTowerToTrack(Double_t trackEta, Double_t trackPhi, Int_t nTowersdEta, Int_t nTowersdPhi) const
Return distance from track to center of one tower.
Int_t getNextId(Int_t det, Int_t m, Int_t e, Int_t s, Int_t nEta, Int_t nPhi) const
Return neighbor id (works for all detectors 1=bemc, 2=bprs, 3=bsmde, 4=bsmdp)
Monte Carlo Track class All information on a simulated track is stored in this class: kinematics...
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Short_t charge() const
Returns charge.
Int_t getTowerEtaPhi(Double_t eta, Double_t phi, Float_t *towerEta, Float_t *towerPhi) const
Return tower eta/phi.
Int_t getNextTowerId(Float_t eta, Float_t phi, Int_t nTowersdEta, Int_t nTowersdPhi) const
Return neighbor tower id's.
const StThreeVector< double > & origin() const
-sign(q*B);
StThreeVectorF getPosFromVertex(const StVertex *const vertex, Int_t TowerId) const
Return Position from collision vertex.
const StThreeVectorF & momentum() const
Returns 3-momentum at dca to primary vertex.
Float_t getEtaFromVertex(const StVertex *const vertex, Int_t TowerId) const
Return eta of the tower considering the collision vertex.
Float_t getPhiFromVertex(const StVertex *const vertex, Int_t TowerId) const
Return phi of the tower considering the collision vertex.
Int_t getBin(const Float_t phi, const Float_t eta, Int_t &m, Int_t &e, Int_t &s) const
Bool_t projTrack(StThreeVectorD *atFinal, StThreeVectorD *momentumAtFinal, const StTrack *const track, Double_t magField, Double_t radius=225.405, Int_t option=1) const
Track projection utility.
StPhysicalHelixD outerHelix() const
Returns outer helix (last measured point)