00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00035
00036
00037
00039 #include "TTree.h"
00040 #include "StEvent/StEvent.h"
00041 #include "StMcEventMaker/StMcEventMaker.h"
00042 #include "StAssociationMaker/StAssociationMaker.h"
00043 #include "StAssociationMaker/StTrackPairInfo.hh"
00044 #include "StTrack.h"
00045 #include "StGlobalTrack.h"
00046 #include "StXiVertex.h"
00047 #include "StXiMuDst.hh"
00048 #include "StXiMc.hh"
00049 #include "StMcEventTypes.hh"
00050 #include "StParticleDefinition.hh"
00051 #include "StTrackDetectorInfo.h"
00052
00053 #include "StStrangeControllerInclude.h"
00054
00055 class StStrangeEvMuDst;
00056
00057
00058 StXiController::StXiController() : StStrangeControllerBase(xiT) {
00059 }
00060
00061 StXiController::~StXiController() {
00062 }
00063
00064 Int_t StXiController::MakeReadDst() {
00065
00066 Int_t j;
00067 StStrangeEvMuDst* ev = masterMaker->GetEvent();
00068 entries = GetN();
00069 for (j = 0; j<entries; j++) {
00070 StXiMuDst* xi = (StXiMuDst*) (*dataArray)[j];
00071 xi->SetEvent(ev);
00072 }
00073 PrintNumCand("read",entries);
00074 nEntries += entries;
00075
00076 if (doMc) {
00077 ev = masterMaker->GetMcEvent();
00078 Int_t mc_entries = GetNMc();
00079 for (j = 0; j<mc_entries; j++) {
00080 StXiMc* mc_xi = (StXiMc*) (*mcArray)[j];
00081 mc_xi->SetEvent(ev);
00082 }
00083 }
00084
00085 return kStOK;
00086 }
00087
00088 Int_t StXiController::MakeCreateDst(StEvent& event) {
00089
00090
00091 Int_t nBad = 0;
00092 StSPtrVecXiVertex& xiVertices = event.xiVertices();
00093 entries = xiVertices.size();
00094 Int_t asize = dataArray->GetSize();
00095 if (entries > asize) dataArray->Expand(entries+increment);
00096 StStrangeEvMuDst* ev = masterMaker->GetEvent();
00097 Int_t j=0;
00098 for (Int_t i=0; i<entries; i++) {
00099 StXiVertex* xiVertex = xiVertices[i];
00100 if (xiVertex) {
00101 StV0Vertex* v0Vertex = xiVertex->v0Vertex();
00102 if (v0Vertex) {
00103 new((*dataArray)[j++]) StXiMuDst(xiVertex,v0Vertex,ev);
00104 } else nBad++;
00105 }
00106 }
00107 entries = j;
00108 PrintNumCand("found",entries);
00109 if (nBad) gMessMgr->Warning() << "StXiController: " << nBad
00110 << "with missing V0 vertices" << endm;
00111 nEntries += entries;
00112
00113 return kStOK;
00114 }
00115
00116 Int_t StXiController::MakeCreateMcDst(StMcVertex* mcVert) {
00117
00118 mcXiMapType* theMcXiMap = 0;
00119 mcTrackMapType* theMcTrackMap = 0;
00120 if (assocMaker) {
00121 theMcXiMap = assocMaker->mcXiMap();
00122 theMcTrackMap = assocMaker->mcTrackMap();
00123 }
00124 if (!((assocMaker)&&(theMcXiMap)&&(theMcTrackMap))) return kStOk;
00125 StStrangeEvMuDst* ev = masterMaker->GetMcEvent();
00126 const StXiVertex* rcXiPartner = 0;
00127 StMcTrack* Bach = 0;
00128 StMcTrack* V0daughter = 0;
00129 Int_t indexRecoArray = -1;
00130 Int_t count = theMcXiMap->count(mcVert);
00131 StSPtrVecMcTrack& Daughters = mcVert->daughters();
00132
00133 for (StMcTrackIterator DTrackIt = Daughters.begin();
00134 DTrackIt != Daughters.end(); DTrackIt++) {
00135 if ((Int_t)(*DTrackIt)->particleDefinition()->charge())
00136 Bach = (*DTrackIt);
00137 else
00138 V0daughter = (*DTrackIt);
00139 }
00140
00141 if (Bach) {
00142 StXiMc* xiMc = new((*mcArray)[mcEntries++]) StXiMc(mcVert,Bach,ev);
00143 if (V0daughter) {
00144 StStrangeControllerBase* v0Cont = masterMaker->Get(v0T);
00145 if (v0Cont) {
00146 StMcVertex* v0Vertex = V0daughter->stopVertex();
00147 if (v0Vertex) {
00148 Int_t before = v0Cont->GetNMc();
00149 v0Cont->MakeCreateMcDst(v0Vertex);
00150 Int_t after = v0Cont->GetNMc();
00151 if (!(before==after)) xiMc->SetV0Index(before);
00152 }
00153 }
00154 }
00155 if (count>0) {
00156 pair<mcXiMapIter,mcXiMapIter> mcXiBounds =
00157 theMcXiMap->equal_range(mcVert);
00158 indexRecoArray = -1;
00159 for(mcXiMapIter mcXiMapIt = mcXiBounds.first;
00160 mcXiMapIt != mcXiBounds.second; ++mcXiMapIt) {
00161 rcXiPartner = (*mcXiMapIt).second;
00162 float x, y, z, delta, xd, yd, zd;
00163 x = mcVert->position().x();
00164 y = mcVert->position().y();
00165 z = mcVert->position().z();
00166 xd = x - rcXiPartner->position().x();
00167 yd = y - rcXiPartner->position().y();
00168 zd = z - rcXiPartner->position().z();
00169 delta = xd*xd + yd*yd + zd*zd;
00170
00171
00172 for(mcXiMapIter mcXiMapIt = mcXiBounds.first;
00173 mcXiMapIt != mcXiBounds.second; ++mcXiMapIt) {
00174 const StXiVertex *temp = (*mcXiMapIt).second;
00175 if (temp != rcXiPartner) {
00176 xd = x - temp->position().x();
00177 yd = y - temp->position().y();
00178 zd = z - temp->position().z();
00179 float delta2 = xd*xd + yd*yd + zd*zd;
00180 if (delta2 < delta) { rcXiPartner = temp; delta = delta2; }
00181 }
00182 }
00183 x = rcXiPartner->position().x();
00184 y = rcXiPartner->position().y();
00185 z = rcXiPartner->position().z();
00186
00187 for(Int_t i = 0; i <= GetN(); i++) {
00188 StXiMuDst* tmpXi = (StXiMuDst*) dataArray->At(i);
00189 if( fabs(x - tmpXi->decayVertexXiX()) < 0.00001 &&
00190 fabs(y - tmpXi->decayVertexXiY()) < 0.00001 &&
00191 fabs(z - tmpXi->decayVertexXiZ()) < 0.00001 )
00192 { indexRecoArray = i; break; }
00193 }
00194 }
00195 new((*assocArray)[assocEntries++])
00196 StStrangeAssoc(indexRecoArray,(mcEntries-1));
00197 if(indexRecoArray!=-1) {
00198 pair<mcTrackMapIter,mcTrackMapIter> mcTrackBounds =
00199 theMcTrackMap->equal_range(Bach);
00200 StTrackPairInfo* bestPairInfo = (*mcTrackBounds.first).second;
00201 for(mcTrackMapIter mcMapIt = mcTrackBounds.first;
00202 mcMapIt != mcTrackBounds.second; ++mcMapIt) {
00203 if ((*mcMapIt).second->commonTpcHits() > bestPairInfo->commonTpcHits())
00204 bestPairInfo = (*mcMapIt).second;
00205 }
00206 if (mcTrackBounds.first != mcTrackBounds.second) {
00207 xiMc->SetHitInfo(bestPairInfo->commonTpcHits());
00208 }
00209 }
00210 }
00211 }
00212
00213 return kStOK;
00214 }
00215
00216 ClassImp(StXiController)
