00001
00002
00003
00004
00006
00007 #include "StRedoTracks.h"
00008 #include "StDbUtilities/StMagUtilities.h"
00009 #include "StEventTypes.h"
00010 #include "StMessMgr.h"
00011 #include "StTpcDb/StTpcDb.h"
00012 #include "StTpcDb/StTpcDbMaker.h"
00013 #include "St_db_Maker/St_db_Maker.h"
00014
00015
00016 ClassImp(StRedoTracks)
00017
00018
00019 StRedoTracks::StRedoTracks(const char *name, StTpcDbMaker* mkr):StMaker(name),
00020 m_ExB(0), tpcDbMaker(mkr), redo(kTRUE) {
00021 }
00022
00023 StRedoTracks::~StRedoTracks() {}
00024
00025 Int_t StRedoTracks::Init(){
00026
00027 if (!tpcDbMaker) {
00028 StMakerIter iter(GetParentChain());
00029 StMaker* mkr;
00030 while ((mkr = iter.NextMaker())) {
00031 if (mkr->IsA() == StTpcDbMaker::Class()) {
00032 tpcDbMaker = (StTpcDbMaker*) mkr;
00033 break;
00034 }
00035 }
00036 if (!tpcDbMaker) gMessMgr->Warning("StRedoTracks: No StTpcDbMaker found.");
00037 }
00038 return StMaker::Init();
00039 }
00040
00041 Int_t StRedoTracks::Make(){
00042
00043 int option = 0;
00044 if (!m_ExB) {
00045 TDataSet *RunLog = GetDataBase("RunLog");
00046 if (!RunLog) gMessMgr->Warning("StRedoTracks: No RunLog found.");
00047 m_ExB = new StMagUtilities(tpcDbMaker->tpcDbInterface(),RunLog,option);
00048 }
00049
00050 StEvent* event = (StEvent*) GetInputDS("StEvent");
00051 if (!event) {
00052 gMessMgr->Warning("StRedoTracks: no StEvent; skipping event.");
00053 return kStWarn;
00054 }
00055
00056 unsigned int i,j,k;
00057 float x[3],p[3],x_new[3],p_new[3];
00058 StTrackGeometry* triGeom = 0;
00059 StThreeVectorD ooo;
00060 StTrackType typ;
00061 Prime typ2;
00062
00063 StPrimaryVertex* pvtx = event->primaryVertex();
00064 if (!pvtx) return kStOk;
00065 ooo = pvtx->position();
00066 StSPtrVecTrackNode& theNodes = event->trackNodes();
00067
00068
00069
00070 UInt_t nPrims = pvtx->numberOfDaughters();
00071 Float_t pv_err = TMath::Sqrt(0.0004 + (0.11/nPrims));
00072
00073 for (i=0; i<theNodes.size(); i++) {
00074 typ = global; typ2 = IsGlobal;
00075 Bool_t iterate = kTRUE;
00076 while (iterate) {
00077 for (j=0; j<theNodes[i]->entries(typ); j++) {
00078 StTrack* tri = theNodes[i]->track(typ,j);
00079 const StTrackTopologyMap& map = tri->topologyMap();
00080 for (k=0; k<2; k++) {
00081 if (k) triGeom = tri->outerGeometry();
00082 else triGeom = tri->geometry();
00083
00084 StThreeVectorF xvec = triGeom->origin();
00085 if (!(xvec.x() || xvec.y() || xvec.z())) continue;
00086 StThreeVectorF pvec = triGeom->momentum();
00087 if (!(pvec.x() || pvec.y())) continue;
00088
00089 float oldPt = pvec.perp();
00090 if (oldPt < 0.0001) continue;
00091 p[0] = pvec.x();
00092 p[1] = pvec.y();
00093 p[2] = pvec.z();
00094 x[0] = xvec.x();
00095 x[1] = xvec.y();
00096 x[2] = xvec.z();
00097
00098 if (!redo) continue;
00099 m_ExB->FixSpaceChargeDistortion(triGeom->charge(),x,p,
00100 typ2,x_new,p_new,map.data(0),map.data(1),pv_err);
00101
00102 StThreeVectorF npvec(p_new);
00103 StThreeVectorF nxvec(x_new);
00104 float newPt = npvec.perp();
00105 float inv_newPt = 1.0/newPt;
00106 float psi = TMath::ACos(npvec.x()*inv_newPt);
00107 if (npvec.y() < 0) psi = TMath::TwoPi() - psi;
00108
00109 triGeom->setMomentum(npvec);
00110 triGeom->setOrigin(nxvec);
00111 triGeom->setCurvature(triGeom->curvature()*pvec.perp()*inv_newPt);
00112 triGeom->setDipAngle(TMath::ATan(npvec.z()*inv_newPt));
00113 triGeom->setPsi(psi);
00114
00115 }
00116 }
00117 if (typ == global) { typ = primary; typ2 = IsPrimary; }
00118 else iterate = kFALSE;
00119 }
00120 }
00121
00122
00123 return kStOK;
00124 }
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134
00135
00136
00137