2 #include "StPrimaryVertex.h"
3 #include "StEventInfo.h"
4 #include "StEventSummary.h"
6 #include "StTrackNode.h"
7 #include "StPrimaryTrack.h"
8 #include "StGlobalTrack.h"
9 #include "StTrackDetectorInfo.h"
10 #include "StTrackGeometry.h"
13 #include "TGeoMatrix.h"
14 #include "StarRoot/THelixTrack.h"
19 #include "TDirectory.h"
23 #include "TRSymMatrix.h"
24 #include "StSvtBarrelHitCollection.h"
25 #include "StSvtHitCollection.h"
26 #include "StSvtLadderHitCollection.h"
27 #include "StSsdLadderHitCollection.h"
28 #include "StSsdWaferHitCollection.h"
29 #include "StSsdHitCollection.h"
30 #include "StDbUtilities/St_svtRDOstrippedC.h"
31 #include "StDedxPidTraits.h"
35 #include "StDbUtilities/St_svtHybridDriftVelocityC.h"
41 TClonesArray *EventT::fgTracks = 0;
42 TClonesArray *EventT::fgHits = 0;
43 THashList *EventT::fRotList = 0;
46 static Int_t _debug = 0;
48 EventT::EventT() : fIsValid(kFALSE)
54 if (!fgTracks) fgTracks =
new TClonesArray(
"TrackT", 1000);
57 if (!fgHits) fgHits =
new TClonesArray(
"HitT", 1000);
70 Int_t EventT::Build(
StEvent *pEventT, UInt_t MinNoHits, Double_t pCut) {
71 static const Int_t NoFitPointCutForGoodTrackT = 15;
74 if (! pEventT)
return iok;
75 UInt_t NprimVtx = pEventT->numberOfPrimaryVertices();
76 if (! NprimVtx)
return iok;
79 Int_t nBestTracks = -1;
81 for (UInt_t ipr=0; ipr < NprimVtx ; ipr++) {
82 pVertex = pEventT->primaryVertex(ipr);
83 if (! pVertex)
continue;
84 UInt_t nDaughters = pVertex->numberOfDaughters();
86 for (UInt_t i=0; i < nDaughters; i++) {
87 StTrack* pTrackT = pVertex->daughter(i);
88 if ( pTrackT->fitTraits().numberOfFitPoints(kTpcId) >= NoFitPointCutForGoodTrackT) nGoodTpcTracks++;
90 if (nBestTracks < nGoodTpcTracks) {nBestTracks = nGoodTpcTracks; ibest = ipr;}
92 if (ibest < 0)
return iok;
93 pVertex = pEventT->primaryVertex(ibest);
97 if (! SvtHitCollection && ! SsdHitCollection) { cout <<
"No SVT & SSD HitT Collections" << endl;
return iok;}
99 fVertex[0] = xyzP.x();
100 fVertex[1] = xyzP.y();
101 fVertex[2] = xyzP.z();
102 StMatrixF vCM = pVertex->covariantMatrix();
103 fCovariantMatrix[0] = vCM(1,1);
104 fCovariantMatrix[1] = vCM(1,2);
105 fCovariantMatrix[2] = vCM(2,2);
106 fCovariantMatrix[3] = vCM(1,3);
107 fCovariantMatrix[4] = vCM(2,3);
108 fCovariantMatrix[5] = vCM(3,3);
109 fNPTracks = pVertex->numberOfDaughters();
114 Int_t ev = 0, run = 0, time = 0;
121 Double32_t field = 0;
122 if (summary) field = summary->magneticField();
123 SetHeader(ev,run,time,field);
126 for (UInt_t t = 0; t < fNPTracks; t++) {
127 StTrack *pTrackT = pVertex->daughter(t);
128 if (! pTrackT)
continue;
130 if (! node)
continue;
131 #ifdef __USE_GLOBAL__
133 if (! gTrackT)
continue;
136 if (! dInfo)
continue;
137 static StDetectorId ids[2] = {kSvtId, kSsdId};
138 UInt_t Nsp = dInfo->numberOfPoints(ids[0]) + dInfo->numberOfPoints(ids[1]);
139 if (MinNoHits > 0 && Nsp < MinNoHits)
continue;
140 UInt_t npoints = dInfo->numberOfPoints() + 100*(dInfo->numberOfPoints(ids[0]) + 10*dInfo->numberOfPoints(ids[1]));
142 pTrackT->numberOfPossiblePoints() +
143 100*(pTrackT->numberOfPossiblePoints(ids[0]) + 10*pTrackT->numberOfPossiblePoints(ids[1]));
145 #ifdef __USE_GLOBAL__
148 Double_t pT = g3.perp();
149 Double_t pMom = g3.mag();
150 if (pMom < pCut)
continue;
153 Double_t TanL = 999999;
154 if (TMath::Abs(pT) > 1.e-7) {
155 InvpT = pTrackT->geometry()->charge()/pT;
158 track->SetInvpT(InvpT);
159 track->SetPhi(TMath::ATan2(g3.y(),g3.x()));
160 track->SetTanL(TanL);
161 static const Double_t EC = 2.9979251E-4;
162 Double_t Rho = - EC*InvpT*field;
165 Double_t TrackLength70 = 0;
166 StSPtrVecTrackPidTraits &traits = pTrackT->pidTraits();
167 UInt_t size = traits.size();
169 for (UInt_t i = 0; i < size; i++) {
170 if (! traits[i])
continue;
171 if ( traits[i]->IsZombie())
continue;
173 if (! pid || pid->method() != kTruncatedMeanId)
continue;
175 TrackLength70 = pid->length();
177 track->SetdEdx(I70,TrackLength70);
178 #ifdef __USE_GLOBAL__
179 Double_t pTGl = g3Gl.perp();
181 Double_t InvpTGl = 0;
182 Double_t TanLGl = 999999;
183 if (TMath::Abs(pTGl) > 1.e-7) {
184 InvpTGl = gTrackT->geometry()->charge()/pTGl;
185 TanLGl = g3Gl.z()/pTGl;
187 track->SetInvpTGl(InvpTGl);
188 track->SetPhiGl(TMath::ATan2(g3Gl.y(),g3Gl.x()));
189 track->SetTanLGl(TanLGl);
190 Double_t RhoGl = - EC*InvpT*field;
191 track->SetRhoGl(RhoGl);
192 track->SetNPpoint(nPpoints);
196 track->SetNpoint(npoints);
197 track->SetNPpoint(nPpoints);
199 const Double_t XyzDirRho[6] = {fVertex[0], fVertex[1], fVertex[2], track->GetTanL(), track->GetPhi(), track->GetRho()};
200 const Double_t XyzDirRhoGl[6] = {fVertex[0], fVertex[1], fVertex[2], track->GetTanL(), track->GetPhi(), track->GetRho()};
202 THashList *fRotList = RotMatrices();
203 if (! fRotList)
continue;
204 TIter next(fRotList);
205 TGeoHMatrix *comb = 0;
206 while ((comb = (TGeoHMatrix *) next())) {
207 TString combName(comb->GetName());
208 if (! combName.BeginsWith(
"R"))
continue;
210 sscanf(comb->GetName()+1,
"%04i",&Id);
211 UInt_t Ladder = Id%100;
212 UInt_t Layer = Id/1000;
if (Layer > 7) Layer = 7;
213 UInt_t Wafer = (Id - 1000*Layer)/100;
214 UInt_t Barrel = (Layer - 1)/2 + 1;
216 cout << comb->GetName() <<
"\tLayer/Ladder/Wafer = " << Layer <<
"/" << Ladder <<
"/" << Wafer << endl;
219 static Double_t dz[2] = {3.00, 2.10};
220 static Double_t dx[2] = {3.00, 3.65};
222 if (Layer > 6) k = 1;
225 const Double_t *X = &XyzDirRho[0];
227 comb->MasterToLocal(X,x0);
228 Double_t CosL = 1./TMath::Sqrt(1. + XyzDirRho[3]*XyzDirRho[3]);
229 Double_t SinL = XyzDirRho[3]*CosL;
230 Double_t CosP = TMath::Cos(XyzDirRho[4]);
231 Double_t SinP = TMath::Sin(XyzDirRho[4]);
232 Double_t dir[3], d[3];
236 comb->MasterToLocalVect(dir,d);
237 if (TMath::Abs(d[2]) < 1.e-7)
continue;
238 Double_t s = - x0[2]/d[2];
if (_debug) cout <<
"straight line s " << s << endl;;
241 TRVector tuvP(2,d[0]/d[2],d[1]/d[2]);
244 uvP += DD;
if (_debug) cout <<
"straight line uvP\t" << uvP <<
"\ttuvP \t" << tuvP << endl;
246 if (TMath::Abs(uvP[0]) > dx[k] + 1.0)
continue;
247 if (TMath::Abs(uvP[1]) > dz[k] + 1.0)
continue;
250 Double_t *rot = comb->GetRotationMatrix();
251 Double_t *tra = comb->GetTranslation();
255 Double_t sh = helixI.
pathLength(middle, normal);
if (_debug) cout <<
"StHelix sh " << sh << endl;
256 if (sh <= 0 || sh > 1e3)
continue;
257 StThreeVectorD xyzG = helixI.at(sh);
if (_debug) cout <<
"StHelix xyzG\t" << xyzG << endl;
258 Double_t xyzGPred[3] = {xyzG.x(),xyzG.y(),xyzG.z()};
260 comb->MasterToLocal(xyzGPred,uvPred);
261 TRVector xyzL(3,uvPred);
if (_debug) cout <<
"StHelix xyzL\t" << xyzL << endl;
262 Double_t dirGPred[3] = {helixI.
cx(sh),helixI.cy(sh),helixI.cz(sh)};
264 comb->MasterToLocalVect(dirGPred,dxyzL);
265 Double_t tuvPred[2] = {dxyzL[0]/dxyzL[2], dxyzL[1]/dxyzL[2]};
266 if (_debug) cout <<
"StHelix tU/tV = " << tuvPred[0] <<
"\t" << tuvPred[1] << endl;
267 #ifdef __USE_GLOBAL__
270 sh = helixG.
pathLength(middle, normal);
if (_debug) cout <<
"StHelix sh " << sh << endl;
271 if (sh < 0 || sh > 1e3)
continue;
272 xyzG = helixG.at(sh);
if (_debug) cout <<
"StHelix xyzG\t" << xyzG << endl;
273 Double_t xyzGPredGl[3] = {xyzG.x(),xyzG.y(),xyzG.z()};
274 Double_t uvPredGl[3];
275 comb->MasterToLocal(xyzGPredGl,uvPredGl);
276 TRVector xyzLGl(3,uvPredGl);
if (_debug) cout <<
"StHelix xyzL\t" << xyzL << endl;
277 Double_t dirGPredGl[3] = {helixG.
cx(sh),helixG.cy(sh),helixG.cz(sh)};
278 comb->MasterToLocalVect(dirGPredGl,dxyzL);
279 Double_t tuvPredGl[2] = {dxyzL[0]/dxyzL[2], dxyzL[1]/dxyzL[2]};
280 if (_debug) cout <<
"StHelix tU/tV = " << tuvPredGl[0] <<
"\t" << tuvPredGl[1] << endl;
283 if (TMath::Abs(uvPred[0]) > dx[k] + 1.0)
continue;
284 if (TMath::Abs(uvPred[1]) > dz[k] + 1.0)
continue;
285 StPtrVecHit &hvec = pTrackT->detectorInfo()->hits();
286 UInt_t NoHits = hvec.size();
288 for (UInt_t j = 0; j < NoHits; j++) {
289 StHit *hitc = hvec[j];
290 if (! hitc)
continue;
291 if (hitc->detector() == kSvtId) {
293 if (htSvt->barrel() == Barrel &&
294 htSvt->ladder() == Ladder &&
295 htSvt->wafer() == Wafer) {hit = hitc;
break;}
297 if (hitc->detector() == kSsdId) {
299 if (htSsd->ladder() == Ladder &&
300 htSsd->wafer() == Wafer) {hit = hitc;
break;}
303 HitT *ht = AddHitT();
305 if (uvPred[0] >= 0) Hybrid = 2;
306 ht->SetId(Barrel,Layer,Ladder,Wafer,Hybrid);
308 ht->SetUVPred (uvPred[0],uvPred[1]);
309 ht->SettUVPred(tuvPred[0],tuvPred[1]);
310 ht->SetXyzG(xyzGPred);
311 ht->SetDirG(dirGPred);
312 #ifdef __USE_GLOBAL__
313 ht->SetUVPredGl (uvPredGl[0],uvPredGl[1]);
314 ht->SettUVPredGl(tuvPredGl[0],tuvPredGl[1]);
315 ht->SetXyzGl(xyzGPredGl);
316 ht->SetDirGl(dirGPredGl);
319 SetHitT(ht, hit, comb, track);
320 ht->SetisFitted(t+1);
322 StSPtrVecSvtHit *hitsvt = 0;
323 StSPtrVecSsdHit *hitssd = 0;
324 Int_t NoHitPerTrack = 0;
326 if (! SvtHitCollection)
continue;
328 if (! barrelCollection)
continue;
330 if (! ladderCollection)
continue;
331 if (svtRDOs && svtRDOs->svtRDOstrippedStatus(Barrel,Ladder,Wafer))
continue;
333 if (! waferCollection)
continue;
334 hitsvt = &waferCollection->hits();
336 if (! SsdHitCollection)
continue;
338 if (! ladderCollection)
continue;
341 if (! waferCollection)
continue;
342 hitssd = &waferCollection->hits();
344 if (! hitsvt && ! hitssd)
continue;
345 if (hitsvt) NoHits = hitsvt->size();
346 if (hitssd) NoHits = hitssd->size();
347 ht->SetNofHits(NoHits);
348 if (! NoHits)
continue;
350 for (UInt_t l = 0; l < NoHits; l++) {
352 if (hitsvt) hit = (*hitsvt)[l];
353 if (hitssd) hit = (*hitssd)[l];
359 h->SetHitFlag(UInt_t(hit->flag()));
360 h->SetUVPred (uvPred[0],uvPred[1]);
361 h->SettUVPred(tuvPred[0],tuvPred[1]);
362 h->SetXyzG(xyzGPred);
363 h->SetDirG(dirGPred);
364 #ifdef __USE_GLOBAL__
366 h->SetUVPredGl (uvPredGl[0],uvPredGl[1]);
367 h->SettUVPredGl(tuvPredGl[0],tuvPredGl[1]);
368 h->SetXyzGl(xyzGPredGl);
369 h->SetDirGl(dirGPredGl);
371 SetHitT(h, hit, comb, track);
373 h->SetHitPerTrack(NoHitPerTrack);
384 TrackT *EventT::AddTrackT()
392 TClonesArray &tracks = *fTracks;
398 HitT *EventT::AddHitT()
406 TClonesArray &hits = *fHits;
407 HitT *hit =
new(hits[fNhit++])
HitT();
413 void EventT::Clear(Option_t * )
420 void EventT::Reset(Option_t * )
425 delete fgTracks; fgTracks = 0;
426 delete fgHits; fgHits = 0;
430 void EventT::SetHeader(Int_t i, Int_t run, Int_t date, Double32_t field)
434 fEvtHdr.Set(i, run, date, field);
437 void EventT::Print(Option_t *opt)
const {
438 cout <<
"Run/EventT\t" << fEvtHdr.GetRun() <<
"/" << fEvtHdr.GetEvtNum() <<
"\tDate " << fEvtHdr.GetDate()
439 <<
"\tField " << fEvtHdr.GetField() << endl;
440 cout <<
"Total no. tracks " << GetTotalNoTracks() <<
"\tRecorded tracks " << GetNtrack()
441 <<
"\tRecorded hits " << GetNhit() << endl;
444 cout <<
"Primary vertex " <<
vertex << endl;
445 cout <<
"Its cov. matrix " << cov << endl;
446 for (UInt_t i = 0; i < GetNtrack(); i++) {cout << i <<
"\t"; GetTrackT(i)->Print();}
447 for (UInt_t i = 0; i < GetNhit(); i++) {cout << i <<
"\t"; GetHitT(i)->Print();}
454 Int_t ladder, barrel;
457 const Int_t NRDOS = 36;
458 const RDO_t RDOS[NRDOS] = {
459 {
"L01B1", 1,1, 1},{
"L02B1", 2,1, 2},{
"L03B1", 3,1, 4},{
"L04B1", 4,1, 5},{
"L05B1", 5,1, 7},
460 {
"L06B1", 6,1, 8},{
"L07B1", 7,1,10},{
"L08B1", 8,1,11},{
"L01B2", 1,2, 3},{
"L02B2", 2,2, 3},
461 {
"L03B2", 3,2, 3},{
"L04B2", 4,2, 6},{
"L05B2", 5,2, 6},{
"L06B2", 6,2, 6},{
"L07B2", 7,2, 9},
462 {
"L08B2", 8,2, 9},{
"L09B2", 9,2, 9},{
"L10B2",10,2,12},{
"L11B2",11,2,12},{
"L12B2",12,2,12},
463 {
"L01B3", 1,3, 1},{
"L02B3", 2,3, 1},{
"L03B3", 3,3, 2},{
"L04B3", 4,3, 2},{
"L05B3", 5,3, 4},
464 {
"L06B3", 6,3, 4},{
"L07B3", 7,3, 5},{
"L08B3", 8,3, 5},{
"L09B3", 9,3, 7},{
"L10B3",10,3, 7},
465 {
"L11B3",11,3, 8},{
"L12B3",12,3, 8},{
"L13B3",13,3,10},{
"L14B3",14,3,10},{
"L15B3",15,3,11},
468 UInt_t B = 0, L = 0, l = 0, W = 0, H = 0;
471 if (hit->detector() == kSvtId) {
479 h->SetuvD(ht->localPosition(0), ht->localPosition(1));
480 h->SetAnode(ht->anode());
481 h->SetTimeB(ht->timebucket());
482 for (Int_t r = 0; r < NRDOS; r++) {
483 if (h->Ladder() == RDOS[r].ladder && h->Barrel() == RDOS[r].barrel) {
488 if (rdo) h->SetRDO(rdo);
490 if (d && d->p(B,l,W,H)) {
491 h->SetLM(d->CalcU(B,l,W,H,ht->timebucket(),ht->anode()),
492 d->CalcV(H,ht->anode()));
493 h->SetuHat(d->uHat(B,l,W,H,ht->timebucket()));
496 h->SetUsedInFit(hit->usedInFit());
497 if (hit->detector() == kSsdId) {
504 h->SetLM(-ht->localPosition(0), ht->localPosition(1));
505 h->SetuHat(-ht->localPosition(0));
510 Double_t xyzG[3] = {position.x(),position.y(),position.z()};
511 h->SetGC(xyzG[0],xyzG[1],xyzG[2]);
512 Double_t xyzL[3] = {0,0,0};
513 comb->MasterToLocal(xyzG,xyzL);
515 Double_t uvw[3] = {h->GetU(),h->GetV(),0};
516 comb->LocalToMaster(uvw,xyzG);
518 Double_t *rot = comb->GetRotationMatrix();
519 h->SetWG(rot[2],rot[5],rot[8]);
521 Int_t IdH = fNhit - 1;
522 track->SetHitTId(IdH);
523 Double_t invpT = track->GetInvpT();
524 if (TMath::Abs(invpT) < 1e-7) invpT = 1e-7;
526 h->SetMom(track->GetMomentum());
527 h->SetWG(rot[2],rot[5],rot[8]);
528 TGeoHMatrix *rotL = (TGeoHMatrix *) RotMatrices()->FindObject(Form(
"WL%s",comb->GetName()+1));
529 Double_t xyzLadder[3] = {0,0,0};
532 rotL->LocalToMaster(uvw,xyzLadder);
533 h->SetL(xyzLadder[0],xyzLadder[1],xyzLadder[2]);
534 Double_t uvwP[3] = {h->GetPredU(),h->GetPredV(),0};
535 rotL->LocalToMaster(uvwP,xyzLadder);
536 h->SetXyzL(xyzLadder);
537 #ifdef __USE_GLOBAL__
539 Double_t uvwPGl[3] = {h->GetPredGlU(),h->GetPredGlV(),0};
540 rotL->LocalToMaster(uvwPGl,xyzLadder);
541 h->SetXyzGlL(xyzLadder);
545 cout << Form(
"WL%s",comb->GetName()+1) <<
" has not been found" << endl;
546 h->SetL(xyzLadder[0],xyzLadder[1],xyzLadder[2]);
547 h->SetXyzL(xyzLadder);
548 #ifdef __USE_GLOBAL__
549 h->SetXyzGlL(xyzLadder);
555 void TrackT::Print(Option_t *opt)
const {
556 cout <<
"TrackT: InvpT " << fInvpT <<
"\tTanL " << fTanL
557 <<
"\tPhi " << fPhi <<
"\tRho " << fRho
558 <<
"\tNpoint " << fNpoint <<
"\tNsp " << fNsp << endl;
559 for (UInt_t i = 0; i < fNsp; i++) cout <<
"\t" << fIdHitT[i];
563 void HitT::SetId(Int_t B, Int_t L, Int_t l, Int_t W, Int_t H) {
570 const Int_t NoLayers = 7;
572 static const Geom_t SvtSsdConfig[NoLayers] =
581 static const Int_t ssdSector[20] = {
583 203, 204, 205, 206, 207, 208, 209,
585 413, 414, 415, 416, 417, 418, 419,
588 barrel = B; layer = L; ladder = l; wafer = W; hybrid = H;
589 if (barrel == 0) Id = 7000 + 100*wafer + ladder;
590 else Id = 1000*layer + 100*wafer + ladder;
594 if (ladder > SvtSsdConfig[layer-1].NoLadders/2) sector = 1;
595 }
else {sector = ssdSector[ladder-1]/100 + 1; barrel = 4;}
599 void HitT::SetId(
StHit *shit) {
601 if (shit->detector() == kSvtId) {
609 if (shit->detector() == kSsdId) {
621 void HitT::Print(Option_t *opt)
const {
622 cout <<
"HitT: Id " << Id <<
"\tpT = " << pT <<
"\tmomentum " << pMom << endl;
623 TRVector glob(3,&xG); cout <<
"Global :" << glob << endl;
624 cout <<
"Local u/v/w " << u <<
"/ " << v <<
"/ " << w << endl;
625 cout <<
"Prediction uP/vP " << uP <<
"/ " << vP <<
"\ttuP/tvP " << tuP <<
"/ " << tvP << endl;
628 void EventT::RestoreListOfRotations() {
629 if (fRotList)
return;
630 if (! gDirectory)
return;
631 fRotList =
new THashList(100,0);
632 fRotList->SetOwner();
633 TIter nextkey(gDirectory->GetListOfKeys() );
635 while ((key = (TKey*) nextkey())) {
636 TObject *obj = key->ReadObj();
637 if ( obj->IsA()->InheritsFrom(
"TGeoHMatrix" ) ) {
double cx(double s) const
pointing vector of helix at point s
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)