00001 #include <stdio.h>
00002 #include <stdlib.h>
00003 #include <assert.h>
00004 #include "TROOT.h"
00005 #include "TClass.h"
00006 #if ROOT_VERSION_CODE < 331013
00007 #include "TCL.h"
00008 #else
00009 #include "TCernLib.h"
00010 #endif
00011 #include "TMath.h"
00012 #include "TBaseClass.h"
00013 #include "TDataMember.h"
00014 #include "TMethod.h"
00015 #include "TMethodArg.h"
00016 #include "TDataType.h"
00017 #include "Api.h"
00018 #include "TMemberInspector.h"
00019 #include "TExMap.h"
00020 #include "TCollection.h"
00021 #include "TRegexp.h"
00022 #include "TRandom.h"
00023 #include "TError.h"
00024 #include "TPoints3DABC.h"
00025 #include "TSystem.h"
00026 #include "TPad.h"
00027 #include "TView.h"
00028
00029 #include "StEvent.h"
00030 #include "StObject.h"
00031 #include "StHit.h"
00032 #include "StTrack.h"
00033 #include "StTrackNode.h"
00034 #include "StVertex.h"
00035 #include "StTrackGeometry.h"
00036 #include "StTrackDetectorInfo.h"
00037
00038 #include "StarClassLibrary/BetheBloch.h"
00039 #include "StEvent/StDedxPidTraits.h"
00040
00041 #include "StEventTypes.h"
00042 #include "StProbPidTraits.h"
00043 #include "StTpcDedxPidAlgorithm.h"
00044 #include "THelixTrack.h"
00045
00046 #define __EVENTHELPER_ONLY__
00047 #include "StEventHelper.h"
00048 #undef __EVENTHELPER_ONLY__
00049 #include <map>
00050
00051 #include "StRnDHitCollection.h"
00052 #include "StEtrHitCollection.h"
00053
00054 void Break(){printf("InBreak\n");}
00055
00056 std::map<long,long> myMap;
00057 typedef std::pair <long,long> MyPair;
00058 std::map <long,long> :: const_iterator myFinder;
00059
00060
00061
00062
00063 class StEventInspector : public TMemberInspector {
00064 public:
00065 StEventInspector(TExMap *map,Int_t &count,const char *opt="");
00066 virtual ~StEventInspector(){delete fSkip;};
00067 virtual void Inspect(TClass* cl, const char* parent, const char* name, const void* addr);
00068 void CheckIn(TObject *obj,const char *bwname="");
00069
00070 Int_t &fCount;
00071 TExMap *fMap;
00072 TRegexp *fSkip;
00073 TString fOpt;
00074
00075 };
00076
00077
00078 StEventInspector::StEventInspector(TExMap *map,Int_t &count,const char *opt):fCount(count)
00079 {
00080 fMap = map;
00081 fSkip = 0;
00082 fOpt = opt;
00083 if (fOpt.Length()) fSkip = new TRegexp(fOpt.Data());
00084 }
00085
00086 void StEventInspector::Inspect(TClass* kl, const char* tit , const char* name, const void* addr)
00087 {
00088 if(tit && strchr(tit,'.')) return ;
00089
00090 TString ts;
00091
00092 if (!kl) return;
00093 if (name[0] == '*') name++;
00094 int ln = strcspn(name,"[ ");
00095 TString iname(name,ln);
00096 const char *iName=iname.Data();
00097 if (iName[1]=='P' && strcmp(iName,"fParent" )==0) return;
00098 if (iName[0]=='G' && strcmp(iName,"G__virtualinfo")==0) return;
00099
00100 G__ClassInfo *classInfo = (G__ClassInfo *)kl->GetClassInfo();
00101 if (!classInfo) return;
00102 G__ClassInfo &cl = *classInfo;
00103
00104
00105
00106 G__DataMemberInfo m(cl);
00107 int found=0;
00108 const char *mName=0;
00109 while (m.Next()) {
00110 mName = m.Name();
00111 if (mName[1] != iName[1]) continue;
00112 if (strcmp(mName,iName) ) continue;
00113 found = 1; break;
00114 }
00115 assert(found);
00116
00117
00118
00119
00120 long prop = m.Property() | m.Type()->Property();
00121 if (prop & G__BIT_ISFUNDAMENTAL) return;
00122 if (prop & G__BIT_ISSTATIC) return;
00123 if (prop & G__BIT_ISENUM) return;
00124 if (strcmp(m.Type()->Fullname(),"TObject") && !m.Type()->IsBase("TObject"))
00125 return;
00126
00127 int size = sizeof(void*);
00128 if (!(prop&G__BIT_ISPOINTER)) size = m.Type()->Size();
00129
00130 int nmax = 1;
00131 if (prop & G__BIT_ISARRAY) {
00132 for (int dim = 0; dim < m.ArrayDim(); dim++) nmax *= m.MaxIndex(dim);
00133 }
00134
00135 for(int i=0; i<nmax; i++) {
00136 char *ptr = (char*)addr + i*size;
00137 TObject *obj = (prop&G__BIT_ISPOINTER) ? *((TObject**)ptr) : (TObject*)ptr;
00138 if (!obj) continue;
00139 const char *bwname = obj->ClassName();
00140 if (!bwname[0] || strcmp(bwname,obj->ClassName())==0) {
00141 bwname = name;
00142 int l = strcspn(bwname,"[ ");
00143 if (bwname[l]=='[') {
00144 char cbuf[12]; sprintf(cbuf,"[%02d]",i);
00145 ts.Replace(0,999,bwname,l);
00146 ts += cbuf;
00147 bwname = (const char*)ts;
00148 }
00149 }
00150
00151 CheckIn(obj,bwname);
00152
00153 }
00154
00155 }
00156
00157 void StEventInspector::CheckIn(TObject *obj,const char *bwname)
00158 {
00159 if (!obj) return;
00160 TObject *inobj=0;
00161 if (obj->InheritsFrom(StRefArray::Class())) return;
00162 if (obj->InheritsFrom( StObjLink::Class())) return;
00163 int n;
00164 if (fSkip && (fSkip->Index(obj->ClassName(),&n)>=0)) return;
00165
00166 if (obj->InheritsFrom(TCollection::Class())){
00167 TCollection *tcol = (TCollection*)obj;
00168 TIter next(tcol);
00169 while ((inobj=next())) {CheckIn(inobj);}
00170 return;
00171 }
00172
00173 if (obj->InheritsFrom(StXRef::Class())){
00174 LongKey_t &inmap = (*fMap)(TMath::Hash(&obj,sizeof(void*)),(Long_t)obj);
00175 myFinder = myMap.find((long)obj);
00176 assert((inmap==0) == (myFinder==myMap.end()));
00177
00178 if (inmap) return;
00179 myMap.insert(MyPair((long)obj,1));
00180 inmap = 1;fCount++;
00181 }
00182
00183 if (obj->InheritsFrom(StStrArray::Class())){
00184
00185 if (((StStrArray*)obj)->size()) {
00186
00187
00188 LongKey_t &inmap = (*fMap)(TMath::Hash(&obj,sizeof(void*)),(Long_t)obj);
00189 myFinder = myMap.find((long)obj);
00190 assert((inmap==0) == (myFinder==myMap.end()));
00191 if (inmap) return;
00192 myMap.insert(MyPair((long)obj,2));
00193
00194 inmap = 2; fCount++;
00195 int vecobj = ( obj->IsA() == StSPtrVecObject::Class());
00196
00197 StStrArray *arr = (StStrArray*)obj;
00198 int sz = arr->size();
00199 for (int idx=0;idx<sz; idx++) {
00200 inobj = arr->at(idx);
00201 Int_t count = fCount;
00202 CheckIn(inobj);
00203 if (count==fCount && !vecobj) break;
00204 }
00205 return;
00206 } }
00207
00208 char cbuf[1000];*cbuf=0;
00209 StEventInspector insp(fMap,fCount);
00210 obj->ShowMembers(insp,cbuf);
00211 }
00212
00213 ClassImp(StEventHelper)
00214
00215 StEventHelper::StEventHelper(const TObject *evt,const char *opt)
00216 {
00217 fMap = new TExMap(10000);
00218 myMap.clear();
00219 fObject = 0;
00220 Reset(evt,opt);
00221 }
00222
00223 StEventHelper::~StEventHelper()
00224 {
00225 myMap.clear();
00226 delete fMap; fMap=0;
00227 Clear();
00228 }
00229
00230 void StEventHelper::Clear(Option_t *opt)
00231 {
00232 }
00233
00234 void StEventHelper::Reset(const TObject *evt,const char *opt)
00235 {
00236 fObject = (TObject *)evt;
00237 Clear();
00238 myMap.clear();
00239 fMap->Delete();
00240 if (!fObject) return;
00241 int kount=0;
00242 StEventInspector insp(fMap,kount,opt);
00243 char cbuf[1024];
00244
00245 fObject->ShowMembers(insp,cbuf);
00246 }
00247
00248 int StEventHelper::Kind(const TObject *to)
00249 {
00250 static TClass *klass=0;
00251 static int who=0;
00252 int kind = 0;
00253 TClass *myClass=to->IsA();
00254 if (myClass!=klass) {
00255 klass = myClass; who = 0;
00256 if (klass->InheritsFrom( StHit::Class())) { who=kHIT;}
00257 else if (klass->InheritsFrom( StTrack::Class())) { who=kTRK;}
00258 else if (klass->InheritsFrom(StPtrVecHit::Class())) { who=kHRR;}
00259 else if (klass->InheritsFrom( StVertex::Class())) { who=kVTX;}
00260 else if (klass->InheritsFrom( TObjArray::Class())) { who=kTRR;}
00261 }
00262 kind = who;
00263 if (kind==kHIT) {
00264 StHitHelper hh((StHit*)to);
00265 if (hh.IsUsed()) {kind|=kUSE;} else {kind|=kUNU;}
00266 if (hh.IsFit ()) kind|=kFIT;
00267 return kind;
00268 }
00269 return kind;
00270 }
00271
00272
00273 void StEventHelper::ls(Option_t* option) const
00274 {
00275 typedef struct { int nb; int sz; const char *tenant; } QWE;
00276 QWE *qwe=0;
00277
00278 TExMap map;
00279 TExMapIter it(fMap);
00280 LongKey_t key,val;
00281 while( it.Next(key,val) ) {
00282 if (val != 2) continue;
00283 StStrArray *a = (StStrArray *)(key);
00284 LongKey_t &cnt = map((Long_t)a->IsA());
00285
00286 if (!cnt) {
00287 qwe = new QWE;
00288 cnt = (Long_t)qwe;
00289 qwe->nb=0; qwe->sz=0;qwe->tenant=0;
00290 }
00291 qwe = (QWE*)cnt;
00292 qwe->nb++; qwe->sz += a->size();
00293 if (qwe->tenant==0 && a->size()) {
00294 TObject *to = a->front();
00295 if (to) qwe->tenant = to->ClassName();
00296 }
00297
00298 }
00299 TExMapIter itt(&map);
00300 printf("\n StEvent(%p)\n",(void*)fObject);
00301
00302 while( itt.Next(key,val) ) {
00303 TObject *kl = (TObject *)key;
00304 qwe = (QWE*)val;
00305 printf ("%8d(%8d) - %s (%s)\n",qwe->nb,qwe->sz,kl->GetName(),qwe->tenant);
00306 delete qwe;
00307 }
00308 printf("\n");
00309
00310 }
00311
00312 TObjArray *StEventHelper::SelConts(const char *sel)
00313 {
00314 TObjArray *tarr = new TObjArray;
00315 TRegexp reg(sel);
00316
00317 TExMapIter it(fMap);
00318 LongKey_t key,val;
00319 while( it.Next(key,val) ) {
00320 if (val == 1) continue;
00321 StStrArray *a = (StStrArray *)(key);
00322 if(a->size()==0) continue;
00323 int n =0;
00324 if (reg.Index(a->ClassName(),&n)<0) continue;
00325 tarr->Add(a);
00326 }
00327 return tarr;
00328 }
00329
00330 TObjArray *StEventHelper::SelTracks(const char*,int flag)
00331 {
00332 int trackTypes[]= {global, primary, tpt, secondary, estGlobal, estPrimary,-1};
00333
00334 TObjArray *conts = SelConts("^StSPtrVecTrackNode$");
00335 TObjArray *traks = new TObjArray();
00336 Int_t ilast = conts->GetLast();
00337 for (int idx=0;idx<=ilast;idx++) {
00338 StObjArray *arr = (StObjArray *)conts->At(idx);
00339 if (!arr) continue;
00340 int ntrk = arr->size();
00341 if (!ntrk) continue;
00342 for (int itrk=0;itrk<ntrk;itrk++) {
00343 StTrackNode *tn = (StTrackNode*)arr->at(itrk);
00344 if (!tn) continue;
00345 StTrack *trk = 0;int ity; int bty=kTGB;
00346 for (int jty=0;(ity=trackTypes[jty])>=0;jty++,bty<<=1){
00347 if (!(flag&bty)) continue;
00348 trk=tn->track(ity);
00349 if (!trk) continue;
00350 if (trk->IsZombie()) continue;
00351
00352 if ((flag&kMark2Draw) && !trk->TestBit(kMark2Draw)) continue;
00353 traks->Add(trk);
00354 }
00355 }
00356 }
00357 delete conts;
00358 return traks;
00359 }
00360
00361 TObjArray *StEventHelper::SelHits(const char *RegEx, Int_t un, Int_t flag)
00362 {
00363
00364
00365 TObjArray *conts = SelConts(RegEx);
00366 TObjArray *hits = new TObjArray();
00367 Int_t ilast = conts->GetLast();
00368
00369 for (int idx=0;idx<=ilast;idx++) {
00370 StObjArray *arr = (StObjArray *)conts->At(idx);
00371 if (!arr) continue;
00372 int sz = arr->size();
00373 if (!sz) continue;
00374 if (!arr->at(0)->InheritsFrom(StHit::Class())) continue;
00375 for(int ih=0;ih<sz; ih++) {
00376 StHit *hit = (StHit*)arr->at(ih);
00377 if (!hit) continue;
00378 if (hit->IsZombie()) continue;
00379
00380 if ((flag&kMark2Draw) && !hit->TestBit(kMark2Draw)) continue;
00381 int used = (hit->trackReferenceCount()!=0);
00382 int take = 0;
00383 if ( used && (un&kUSE)) take++;
00384 if (!used && (un&kUNU)) take++;
00385 if (take) hits->Add(hit);
00386 }
00387 }
00388 delete conts;
00389 return hits;
00390 }
00391
00392 TObjArray *StEventHelper::SelVertex(const char *sel,Int_t flag)
00393 {
00394
00395 TObjArray *conts = SelConts(sel);
00396 TObjArray *verts = new TObjArray();
00397 Int_t ilast = conts->GetLast();
00398 int nvtx =0;
00399 for (int idx=0;idx<=ilast;idx++) {
00400 StObjArray *arr = (StObjArray *)conts->At(idx);
00401 if (!arr) continue;
00402 int sz = arr->size();
00403 if (!sz) continue;
00404 for (int ivx=0; ivx<sz; ivx++) {
00405 StVertex *vx = (StVertex*)arr->at(ivx);
00406 if (!vx) continue;
00407 if (vx->IsZombie()) continue;
00408
00409 if ((flag&kMark2Draw) && !vx->TestBit(kMark2Draw)) continue;
00410 verts->Add(vx);nvtx++;
00411 }
00412 }
00413 delete conts;
00414 return verts;
00415 }
00416
00417 TObjArray *StEventHelper::ExpandAndFilter(const TObject *eObj, int flag, TObjArray *out)
00418 {
00419
00420 if (!out) {out = new TObjArray;}
00421 TObject *eobj = (TObject*)eObj;
00422
00423 int kind = Kind(eobj);
00424 if (kind&kHIT) {
00425 if (!(flag&kHIT)) return out;
00426 int take=kind & (kUSE|kUNU|kFIT) &flag;
00427 if (take) out->Add(eobj);
00428 return out;
00429 }
00430
00431
00432 if (kind&kTRK) {
00433 if (flag&kTRK) out->Add(eobj);
00434 if (!(flag&kHRR)) return out;
00435 StTrack *trk = (StTrack*)eobj;
00436 StTrackHelper trkh(trk);
00437 out->Add((TObject*)trkh.GetHits());
00438 return out;
00439 }
00440
00441
00442 if (kind&kVTX) {
00443 if (flag&kVTX) out->Add(eobj);
00444 StVertex *vtx = (StVertex*)eobj;
00445 if (!(flag&(kTRK|kHRR))) return out;
00446 StVertexHelper vtxh(vtx);
00447 int n = vtxh.GetNTracks();
00448 for (int i=-1;i<n;i++) {
00449 const TObject *to = vtxh.GetTrack(i);
00450 if (!to) continue;
00451 ExpandAndFilter(to,flag,out);
00452 }
00453 return out;
00454 }
00455
00456
00457 if (kind&kTRR) {
00458 TObjArray *inp = (TObjArray *)eobj;
00459 inp->Compress();
00460 int nbjs = inp->GetLast()+1;
00461 if (!nbjs) return 0;
00462 for (int i=0;i<nbjs;i++) {
00463 ExpandAndFilter(inp->At(i),flag,out);
00464 }
00465 return out;
00466 }
00467
00468
00469 if (kind&kHRR) {
00470 if (!(flag&kHRR)) return out;
00471 out->Add(eobj);
00472 return out;
00473 }
00474 return 0;
00475 }
00476
00477 TObjArray *StEventHelper::MakePoints(TObjArray *inp, int flag)
00478 {
00479 static const Color_t plitra[]={kRed,kGreen,kBlue,kMagenta, kCyan};
00480 static const int nlitra = sizeof(plitra)/sizeof(Color_t);
00481 int ilitra=0;
00482 inp->Compress();
00483 int nbjs = inp->GetLast()+1;
00484 if (!nbjs) return 0;
00485 TObjArray *out = new TObjArray;out->SetOwner();
00486 StPoints3DABC *p[3] ; int np=0;
00487 for (int i=0;i<nbjs;i++) {
00488 TObject *to = inp->At(i);
00489 int kind = Kind(to);
00490 if (!(kind&kHRR)) { ilitra++; ilitra = ilitra%nlitra; }
00491 int take = (kind&flag);
00492 if (!take) continue;
00493
00494 if (!take) continue;
00495
00496 np = 0;
00497 if (kind&kHIT) {np=1;p[0] = new StHitPoints ((StHit *)to );}
00498 else if (kind&kHRR && ((StPtrVecHit*)to)->size())
00499 {np=1;p[0] = new StHitPoints ((StPtrVecHit*)to );}
00500
00501 else if (kind&kTRK) {np=3;p[0] = new StTrackPoints ((StTrack *)to );
00502 p[1] = new StInnOutPoints((StTrack *)to,0);
00503 p[2] = new StInnOutPoints((StTrack *)to,1);}
00504
00505 else if (kind&kVTX) {np=1;p[0] = new StVertexPoints((StVertex *)to );}
00506
00507 for (int j=0;j<np;j++){p[j]->SetUniqueID(plitra[ilitra]); out->Add(p[j]);}
00508 }
00509
00510 return out;
00511 }
00512
00513 void StEventHelper::Break(int kase)
00514 {
00515 fprintf(stderr,"Break(%d)\n",kase);
00516 }
00517
00518 void StEventHelper::Remove(StEvent *ev,const char *className)
00519 {
00520 StSPtrVecObject& V = ev->content();
00521 int n = V.size();
00522 for (int i=0; i<n; i++) {
00523 StObject *to = V[i];
00524 if (!to) continue;
00525 if (!strstr(to->ClassName(),className)) continue;
00526 V[i] = 0; delete to;
00527 break;
00528 }
00529 }
00530
00531
00532 ClassImp(StPoints3DABC)
00533 void StPoints3DABC::Add(StPoints3DABC *add)
00534 {
00535 int n = add->fSize + fSize;
00536 if (n > fN) {
00537 if (n < fN*2) n = fN*2;
00538 Float_t *arr = new Float_t[n*3];
00539 memcpy(arr, fXYZ, fSize*3*sizeof(Float_t));
00540 delete [] fXYZ; fXYZ = arr; fN = n;
00541 }
00542 memcpy(fXYZ+fSize*3,add->fXYZ,add->fSize*3*sizeof(Float_t));
00543 fSize+=add->fSize;
00544 }
00545
00546 ClassImp(StTrackPoints)
00547
00548 StTrackPoints::StTrackPoints(const StTrack *st,const char *name,const char *title)
00549 :StPoints3DABC(name,title,st)
00550 {
00551 Init();
00552 }
00553
00554 void StTrackPoints::Init()
00555 {
00556 if (fXYZ) return;
00557 StTrack *trk = ((StTrack*)fObj);
00558 StTrackHelper th(trk);
00559 fXYZ = th.GetPoints(fSize);
00560 fN = fSize;
00561 if (!fSize) { MakeZombie(); return;}
00562 }
00563
00564 Int_t StTrackPoints::DistancetoPrimitive(Int_t px, Int_t py)
00565 {
00566
00567
00568
00569
00570
00571
00572
00573
00574
00575
00576 enum {inaxis = 7,mindist=10};
00577 Float_t dist = 999999;
00578
00579 Int_t puxmin = gPad->XtoAbsPixel(gPad->GetUxmin());
00580 Int_t puymin = gPad->YtoAbsPixel(gPad->GetUymin());
00581 Int_t puxmax = gPad->XtoAbsPixel(gPad->GetUxmax());
00582 Int_t puymax = gPad->YtoAbsPixel(gPad->GetUymax());
00583
00584 TView *view = 0;
00585
00586 if (px < puxmin - inaxis) goto END;
00587 if (py > puymin + inaxis) goto END;
00588 if (px > puxmax + inaxis) goto END;
00589 if (py < puymax - inaxis) goto END;
00590
00591 view = gPad->GetView();
00592 if (!view) goto END;
00593
00594 {Int_t i;
00595 Float_t alfa;
00596 Float_t xndc[3];
00597 Int_t x1,y1,x0,y0;
00598 Int_t pointSize = fN*3;
00599 view->WCtoNDC(fXYZ, xndc);
00600 x0 = gPad->XtoAbsPixel(xndc[0]);
00601 y0 = gPad->YtoAbsPixel(xndc[1]);
00602
00603 float dif[2],difdif,cur[2],curcur,difcur;
00604 for (i=3;i<pointSize;i+=3) {
00605 view->WCtoNDC(fXYZ+i, xndc);
00606 x1 = gPad->XtoAbsPixel(xndc[0]);
00607 y1 = gPad->YtoAbsPixel(xndc[1]);
00608 dif[0] = x1-x0; dif[1]=y1-y0;
00609 cur[0] = x0-px; cur[1]=y0-py;
00610 difdif = (dif[0]*dif[0]+dif[1]*dif[1]);
00611 difcur = (dif[0]*cur[0]+dif[1]*cur[1]);
00612 curcur = cur[0]*cur[0]+cur[1]*cur[1];
00613 if (difdif<mindist*mindist) {
00614 if ((i+3)<pointSize) continue;
00615 dist = curcur; break;
00616 }
00617 alfa = -difcur/difdif;
00618
00619 if (alfa<0.) {dist = curcur; break;}
00620
00621 x0=x1; y0=y1;
00622 if (alfa > 1.) {
00623 if (i+3 < pointSize) continue;
00624 dist = (px-x1)*(px-x1) + (py-y1)*(py-y1); break;
00625 }
00626 dist = curcur+alfa*(2*difcur+difdif*alfa);
00627 break;
00628 }}
00629 END:
00630 dist = TMath::Sqrt(dist);
00631
00632 if (dist <= mindist) { dist = 0; gPad->SetSelected(this);}
00633
00634 return Int_t(dist);
00635 }
00636
00637
00638
00639 ClassImp(StVertexPoints)
00640
00641 StVertexPoints::StVertexPoints(const StVertex *sv,const char *name,const char *title)
00642 :StPoints3DABC(name,title,sv)
00643 {
00644 SetBit(1);
00645 fSize = 1; fN =1;
00646 fXYZ = new Float_t[3];
00647 fXYZ[0] = ((StVertex*)fObj)->position().x();
00648 fXYZ[1] = ((StVertex*)fObj)->position().y();
00649 fXYZ[2] = ((StVertex*)fObj)->position().z();
00650 }
00651
00652 ClassImp(StVertexPoints)
00653
00654 StInnOutPoints::StInnOutPoints(const StTrack *st,int innout,const char *name,const char *title)
00655 :StPoints3DABC(name,title,st)
00656 {
00657 fSize = 1; fN =1; fInnOut=innout;
00658 const StTrackGeometry *geo = (fInnOut==0) ? st->geometry():st->outerGeometry();
00659 fXYZ = new Float_t[3];
00660 fXYZ[0] = geo->origin().x();
00661 fXYZ[1] = geo->origin().y();
00662 fXYZ[2] = geo->origin().z();
00663 }
00664
00665
00666 ClassImp(StHitPoints)
00667
00668 StHitPoints::StHitPoints(const StHit *sh,const char *name,const char *title)
00669 :StPoints3DABC(name,title,sh)
00670 {
00671 fSize = 1; fN =1;
00672 Init();
00673 }
00674
00675 StHitPoints::StHitPoints(const StRefArray *ar,const char *name,const char *title)
00676 :StPoints3DABC(name,title,ar)
00677 {
00678 fSize = ar->size();
00679 fN = fSize;
00680 if (!fSize) return;
00681 fObj = (fSize==1) ? ar->front() : ar;
00682 Init();
00683 }
00684
00685 void StHitPoints::Init()
00686 {
00687 if (fXYZ) return;
00688 fXYZ = new Float_t[fN*3];
00689
00690 int n=0;
00691 for (int i =0;i<fSize;i++)
00692 {
00693 StHit *hit= (fSize==1) ? (StHit*)fObj: (StHit*)((StRefArray*)fObj)->at(i);
00694 if (fSize>1 && !hit->trackReferenceCount()) continue;
00695 if (fSize>1 && !hit->usedInFit()) continue;
00696 StThreeVectorF v3 = hit->position();
00697 fXYZ[n*3+0] = v3.x();
00698 fXYZ[n*3+1] = v3.y();
00699 fXYZ[n*3+2] = v3.z();
00700 n++;
00701 }
00702 fN=n; fSize=n;
00703 }
00704
00705
00706 ClassImp(StFilterABC)
00707
00708 int StFilterABC::fgDial=0;
00709
00710 StFilterABC::StFilterABC(const char *name,bool active):TNamed(name,""),fActive(active)
00711 {
00712 }
00713
00714 void StFilterABC::SetDefs()
00715 {
00716 for (int i=0;GetNams() && GetNams()[i]; i++) {GetPars()[i]=GetDefs()[i];}
00717 }
00718
00719
00720 ClassImp(StFilterDef)
00721 StFilterDef::StFilterDef(const char *name,bool active):StFilterABC(name,active)
00722 {
00723 SetDefs();
00724
00725 }
00726
00727 const char **StFilterDef::GetNams() const
00728 {
00729 static const char *nams[] = {
00730 " RandomSelect ",
00731 " RxyMin ",
00732 " RxyMax ",
00733 " ZMin ",
00734 " ZMax ",
00735 " PhiMin ",
00736 " PhiMax ",
00737 " LenMin ",
00738 " LenMax ",
00739 " PtMin ",
00740 " PtMax ",
00741 " PseudoMin ",
00742 " PseudoMax ",
00743 " QMin ",
00744 " QMax ",
00745 " EncodedMethod",
00746 0};
00747 return nams;
00748 }
00749
00750 const float *StFilterDef::GetDefs() const
00751 {
00752 static const float defs[] = {
00753 1.00,
00754 0.00,
00755 900.00,
00756 -900.00,
00757 +900.00,
00758 -180.01,
00759 +181.01,
00760 +0.00,
00761 +999.00,
00762 0.00,
00763 999.00,
00764 -999.00,
00765 999.00,
00766 -1 ,
00767 +1 ,
00768 -1 ,
00769
00770 0};
00771 return defs;
00772 }
00773
00774
00775 Int_t StFilterDef::Accept(StPoints3DABC *pnt,Color_t &color, Size_t&, Style_t&)
00776 {
00777 static TRandom rrr;
00778 float x,y,z,r2xy,phid,len,pt,ps,q;
00779 const TObject *to;
00780 const StTrack *trk;
00781
00782 color = (((color-kRed)+1)%6)+kRed;
00783
00784 int cut = 1;
00785 if (fRandomSelect < 1. && fRandomSelect < rrr.Rndm())return 0;
00786
00787 z = pnt->GetZ(0);
00788 cut++;
00789 if (fZMin >z || z > fZMax) goto SKIP;
00790
00791 x = pnt->GetX(0);
00792 y = pnt->GetY(0);
00793 r2xy = x*x+y*y;
00794 cut++;
00795 if (fRxyMin*fRxyMin > r2xy || r2xy > fRxyMax*fRxyMax)goto SKIP;
00796 phid = atan2(y,x)*(180./M_PI);
00797 cut++;
00798 if (fPhiMin > phid || phid > fPhiMax) goto SKIP;
00799 to = pnt->GetObject();
00800 if (!to) return 1;
00801 if (!to->InheritsFrom(StTrack::Class())) return 1;
00802
00803
00804
00805
00806 trk = (StTrack*)to;
00807 len = trk->length();
00808 cut++;
00809 if (fLenMin >len || len > fLenMax) goto SKIP;
00810 pt = trk->geometry()->momentum().perp();
00811 cut++;
00812 if (fPtMin >pt || pt > fPtMax) goto SKIP;
00813 ps = trk->geometry()->momentum().pseudoRapidity();
00814 cut++;
00815 if (fPsMin >ps || ps > fPsMax) goto SKIP;
00816 q = trk->geometry()->charge();
00817 cut++;
00818 if (fQMin >q || q > fQMax) goto SKIP;
00819 cut++;
00820 if ( (int(fEncodedMethod) != -1) && (trk->encodedMethod() != int(fEncodedMethod)) )
00821 goto SKIP;
00822 return 1;
00823
00824 SKIP: return 0;
00825
00826 }
00827
00828
00829
00830 ClassImp(StMuDstFilterHelper)
00831 StMuDstFilterHelper::StMuDstFilterHelper(const char *name,bool active):StFilterABC(name,active)
00832 {
00833 mBB = new BetheBloch();
00834 SetDefs();
00835
00836 }
00837
00838 StMuDstFilterHelper::~StMuDstFilterHelper()
00839 { delete mBB;}
00840
00841 const char **StMuDstFilterHelper::GetNams() const
00842 {
00843 static const char *nams[] = {
00844 " pCutHigh ",
00845 " nHitsCutHighP ",
00846 " pCutLow ",
00847 " nHitsCutLowP ",
00848 " chargeForLowP ",
00849 " dEdxMassCutHigh ",
00850 " dEdxFractionCutHigh ",
00851 " dEdxMassCutLow ",
00852 " dEdxFractionCutLow ",
00853 0
00854 };
00855 return nams;
00856 }
00857
00858 const float *StMuDstFilterHelper::GetDefs() const
00859 {
00860 static const float defs[] = {
00861 2.0,
00862 10,
00863 0.2,
00864 15,
00865 -1,
00866 0.939,
00867 0.6,
00868 0.494,
00869 1.1,
00870 0
00871 };
00872 return defs;
00873 }
00874
00875 Int_t StMuDstFilterHelper::Accept(const StTrack* track) {
00876
00877 float pCutHigh = fpCutHigh;
00878 int nHitsCutHighP = int(fnHitsCutHighP);
00879
00880
00881 float pCutLow = fpCutLow;
00882 int nHitsCutLowP = int(fnHitsCutLowP);
00883 int chargeForLowP = int(fchargeForLowP);
00884 float dEdxMassCutHigh = fdEdxMassCutHigh;
00885 float dEdxFractionCutHigh = fdEdxFractionCutHigh;
00886 float dEdxMassCutLow = fdEdxMassCutLow;
00887 float dEdxFractionCutLow = fdEdxFractionCutLow;
00888
00889 int iret = 0;
00890 int chargeOK = 0;
00891 int dedxOK = 0;
00892
00893 float magnitude = track->geometry()->momentum().magnitude();
00894 int nPoints = track->detectorInfo()->numberOfPoints();
00895
00896 if ( magnitude > pCutHigh && nPoints >= nHitsCutHighP) iret = 1;
00897 else {
00898 if ( magnitude > pCutLow && nPoints >= nHitsCutLowP )
00899 {
00900
00901 if (chargeForLowP==0)
00902 chargeOK = 1;
00903 else if (track->geometry()->charge() == chargeForLowP)
00904 chargeOK = 1;
00905
00906
00907
00908 float dedxHigh = dEdxFractionCutHigh * mBB->Sirrf(magnitude/dEdxMassCutHigh);
00909 float dedxLow = dEdxFractionCutLow * mBB->Sirrf(magnitude/dEdxMassCutLow);
00910 float dedx = 0;
00911
00912
00913 const StSPtrVecTrackPidTraits& traits = track->pidTraits();
00914 StDedxPidTraits* dedxPidTr;
00915 for (unsigned int itrait = 0; itrait < traits.size(); itrait++){
00916 dedxPidTr = 0;
00917 if (traits[itrait]->detector() == kTpcId) {
00918 StTrackPidTraits* thisTrait = traits[itrait];
00919 dedxPidTr = dynamic_cast<StDedxPidTraits*>(thisTrait);
00920 if (dedxPidTr && dedxPidTr->method() == kTruncatedMeanId) {
00921
00922 dedx = 2 * dedxPidTr->mean();
00923 }
00924 }
00925 }
00926 if (dedx > dedxHigh && dedx > dedxLow)
00927 dedxOK = 1;
00928
00929 iret = chargeOK * dedxOK;
00930 }
00931 }
00932 return iret;
00933 }
00934
00935
00936 Int_t StMuDstFilterHelper::Accept(StPoints3DABC *pnt)
00937 {
00938 const TObject *to;
00939 const StTrack *trk;
00940 to = pnt->GetObject();
00941 if (!to) return 1;
00942 if (!to->InheritsFrom(StTrack::Class())) return 1;
00943 trk = (StTrack*)to;
00944 return Accept(trk);
00945 }
00946
00947
00948 ClassImp(StColorFilterHelper)
00949 StColorFilterHelper::StColorFilterHelper(const char *name,bool active):StFilterABC(name,active)
00950 {
00951 fPidAlgorithm = new StTpcDedxPidAlgorithm();;
00952 fElectron = StElectron::instance();
00953 fPion = StPionPlus::instance();
00954 fKaon = StKaonPlus::instance();
00955 fProton = StProton::instance();
00956
00957 SetDefs();
00958
00959 }
00960
00961 StColorFilterHelper::~StColorFilterHelper()
00962 { delete fPidAlgorithm;}
00963
00964 const char **StColorFilterHelper::GetNams() const
00965 {
00966 static const char *nams[] = {
00967 " Electron sigma ",
00968 " Electron color ",
00969 " Pion sigma ",
00970 " Pion color ",
00971 " Kaon sigma ",
00972 " Kaon color ",
00973 " Proton sigma ",
00974 " Proton color ",
00975 " others sigma ",
00976 " others color ",
00977 0
00978 };
00979 return nams;
00980 }
00981
00982 const float *StColorFilterHelper::GetDefs() const
00983 {
00984 static const float defs[] = {
00985 1 ,
00986 2 ,
00987 1 ,
00988 3 ,
00989 1 ,
00990 4 ,
00991 1 ,
00992 6 ,
00993
00994 -1,
00995 0,
00996 0
00997 };
00998 return defs;
00999 }
01000
01001 Int_t StColorFilterHelper::Accept(const StTrack* track, Color_t &color, Size_t&size, Style_t&) {
01002
01003 float sigmaElectron = fNSigmaElectron ;
01004 Color_t colorElectron = (Color_t)fNColorElectron ;
01005
01006 float sigmaPion = fNSigmaPion ;
01007 Color_t colorPion = (Color_t)fNColorPion ;
01008
01009 float sigmaKaon = fNSigmaKaon ;
01010 Color_t colorKaon = (Color_t)fNColorKaon ;
01011
01012 float sigmaProton = fNSigmaProton ;
01013 Color_t colorProton = (Color_t)fNColorProton ;
01014
01015
01016 Color_t colorOther = (Color_t)fNColorOther ;
01017
01018
01019
01020 track->pidTraits(*fPidAlgorithm);
01021
01022 color = colorOther;
01023 size = 1;
01024
01025 if (TMath::Abs(fPidAlgorithm->numberOfSigma(fElectron)) < sigmaElectron)
01026 { color = colorElectron; size = 2; }
01027
01028 if (TMath::Abs(fPidAlgorithm->numberOfSigma(fKaon)) < sigmaKaon)
01029 { color = colorKaon; size = 4; }
01030
01031 if (TMath::Abs(fPidAlgorithm->numberOfSigma(fPion)) < sigmaPion)
01032 { color = colorPion; size = 5; }
01033
01034 if (TMath::Abs(fPidAlgorithm->numberOfSigma(fProton)) < sigmaProton)
01035 { color = colorProton; size = 3; }
01036
01037 return 1;
01038 }
01039
01040
01041 Int_t StColorFilterHelper::Accept(StPoints3DABC *pnt, Color_t&color, Size_t&size, Style_t&style)
01042 {
01043 const TObject *to;
01044 const StTrack *trk;
01045 to = pnt->GetObject();
01046 if (!to) return 1;
01047 if (!to->InheritsFrom(StTrack::Class())) return 1;
01048 trk = (StTrack*)to;
01049 return Accept(trk,color,size,style);
01050 }
01051
01052
01053
01054 ClassImp(StVertexHelper)
01055 StVertexHelper::StVertexHelper(const StVertex *vtx)
01056 { SetVertex(vtx);}
01057 StVertexHelper::StVertexHelper(const StEvent *evt)
01058 { SetVertex(evt->primaryVertex(0));}
01059
01060 void StVertexHelper::SetVertex(const StVertex *vtx){fVtx = vtx;}
01061 int StVertexHelper::GetType() {return (int)fVtx->type();}
01062 int StVertexHelper::GetFlag() {return fVtx->flag();};
01063 int StVertexHelper::GetNTracks() {return fVtx->numberOfDaughters();}
01064
01065 const StThreeVectorF &StVertexHelper::GetPoint()
01066 {
01067 return fVtx->position();
01068 }
01069
01070 const StTrack *StVertexHelper::GetTrack(int idx)
01071 {
01072 if (idx==-1) return fVtx->parent();
01073 if (idx>= GetNTracks()) return 0;
01074 return fVtx->daughter((UInt_t)idx);
01075 }
01076
01077 const float *StVertexHelper::GetErrMtx()
01078 {
01080
01081 StMatrixF mxF = fVtx->covariantMatrix();
01082 int jj=0;
01083 for (int i=0;i< 3;i++) {
01084 for (int j=0;j<=i;j++) {
01085 fErrMtx[jj++] = mxF(i+1,j+1);}}
01086 return fErrMtx;
01087 }
01088
01089
01090
01091 ClassImp(StTrackHelper)
01092 StTrackHelper::StTrackHelper(const StTrack *trk) :
01093 StHelixHelper(trk->geometry()->helix(),
01094 trk->outerGeometry()->helix(),trk->length())
01095 , fTrk(trk), fHits(0)
01096 { GetNHits(); }
01097
01098 StTrackHelper::~StTrackHelper()
01099 { }
01100 int StTrackHelper::GetType() const {return fTrk->type();}
01101 int StTrackHelper::GetFlag() const {return fTrk->flag();}
01102 int StTrackHelper::GetCharge() const {return fTrk->geometry()->charge();}
01103 const StVertex *StTrackHelper::GetParent() const {return fTrk->vertex();}
01104 float StTrackHelper::GetImpact() const {return fTrk->impactParameter();}
01105 float StTrackHelper::GetCurv() const {return GetTHelix(0)->GetRho() ;}
01106 const StThreeVectorF &StTrackHelper::GetFirstPoint() const {return fTrk->geometry()->origin();}
01107 const StThreeVectorF &StTrackHelper::GetLastPoint() const {return fTrk->outerGeometry()->origin();}
01108 const StThreeVectorF &StTrackHelper::GetMom() const {return fTrk->geometry()->momentum();}
01109
01110
01111 const StPtrVecHit *StTrackHelper::GetHits() const
01112 {
01113 if (fHits) return fHits;
01114 const StTrackDetectorInfo *tdi = fTrk->detectorInfo();
01115 if (!tdi) return 0;
01116 fHits = &tdi->hits();
01117 return fHits;
01118 }
01119
01120 int StTrackHelper::GetNHits() const
01121 {
01122 if (fHits) return fHits->size();
01123 GetHits();
01124 return (fHits)? fHits->size():0;
01125 }
01126
01127 const StHit *StTrackHelper::GetHit(int idx) const
01128 {
01129 if (idx<0) return 0;
01130 if (idx>=GetNHits()) return 0;
01131 if (!fHits) return 0;
01132 return fHits->at(idx);
01133 }
01134
01135 int StTrackHelper::numberOfFitPoints(int det) const
01136 {
01137 const StTrackFitTraits& trait = fTrk->fitTraits();
01138 return (det)? trait.numberOfFitPoints((StDetectorId)det): trait.numberOfFitPoints();
01139 }
01140
01141 StMCTruth StTrackHelper::GetTruth(int byNumb,double rXYMin,double rXYMax) const
01142 {
01143 StMCPivotTruth pivo(1);
01144 int nHits = GetNHits();
01145 int nUsed=0;
01146 for (int jh=0;jh<nHits;jh++) {
01147 const StHit *hit = GetHit(jh);
01148 double r = sqrt(pow(hit->position().x(),2)+pow(hit->position().y(),2));
01149 if (r<rXYMin) continue;
01150 if (r>rXYMax) continue;
01151 int idTruth=hit->idTruth();
01152 int wtTruth=hit->qaTruth();
01153 if (!wtTruth) wtTruth=1;
01154
01155
01156
01157 nUsed++; pivo.Add(idTruth,wtTruth);
01158 }
01159 if (!nUsed) return 0;
01160 return pivo.Get(byNumb);
01161 }
01162
01163
01164
01165
01166 ClassImp(StHitHelper)
01167 StHitHelper::StHitHelper(const StHit *hit){fHit = hit;}
01168 void StHitHelper::SetHit(const StHit *hit) {fHit = hit;}
01169 int StHitHelper::GetDetId() {return fHit->detector();}
01170 int StHitHelper::GetFlag() {return fHit->flag();}
01171 float StHitHelper::GetCharge() {return fHit->charge();}
01172 int StHitHelper::IsUsed() {return fHit->trackReferenceCount();}
01173 int StHitHelper::IsFit() {return fHit->usedInFit();}
01174 const StThreeVectorF &StHitHelper::GetPoint() {return fHit->position();}
01175
01176
01177
01178
01179 ClassImp(StErrorHelper)
01180
01181 StErrorHelper::StErrorHelper()
01182 {
01183 fNErr=0; fNTot=0; fKErr=0;
01184 fMap = new TExMap;
01185 fArr = new TArrayI;
01186 }
01187
01188 StErrorHelper::~StErrorHelper()
01189 {
01190 delete fMap;
01191 delete fArr;
01192 }
01193
01194 void StErrorHelper::Add(int errn)
01195 {
01196 fNTot++;
01197 if(!errn) return;
01198 fNErr++;
01199 (*fMap)(errn)++;
01200 }
01201
01202 void StErrorHelper::MakeArray()
01203 {
01204 if (!fNErr) return;
01205 fKErr = fMap->GetSize();
01206 fArr->Set(fKErr*3);
01207 TExMapIter it(fMap);
01208 LongKey_t lerr,lnum;
01209
01210 int idx=0;
01211 while(it.Next(lerr,lnum)) {
01212 (*fArr)[idx+ 0] = lnum;
01213 (*fArr)[idx+fKErr] = lerr;
01214 idx++;
01215 }
01216
01217 TMath::Sort(fKErr, fArr->GetArray(), fArr->GetArray()+2*fKErr);
01218 }
01219
01220 void StErrorHelper::Print(const char* txt) const
01221 {
01222 StErrorHelper *This = (StErrorHelper *)this;
01223 This->MakeArray();
01224 if (!txt) txt="";
01225 printf("StEvent Error Summary:%s\n",txt);
01226
01227 printf("%4d -%8d(%4d)\n",0,0,fNTot-fNErr);
01228 int *nrr=fArr->GetArray();
01229 int *krr=nrr+fKErr;
01230 int *idx=krr+fKErr;
01231 for (int i=0;i<fKErr;i++) {
01232 int j = idx[i];
01233 printf("%4d -%8d(%4d) //%s\n",i+1,krr[j],nrr[j],Say(krr[j]).Data());
01234 }
01235 }
01236
01237
01238 TString StErrorHelper::Say(int ierr,const char *klass)
01239 {
01240 static const char *TabErr[] =
01241 {
01242 "StTrack" ,"mFlag" ,"1","2","is Negative",
01243 "StTrack" ,"mFlag" ,"1","3","is Zero",
01244
01245 "StTrack" ,"mImpactParameter" ,"2","1","is NaN",
01246 "StTrack" ,"mImpactParameter" ,"2","2","is huge",
01247
01248 "StTrack" ,"mLength" ,"3","1","is NaN",
01249 "StTrack" ,"mLength" ,"3","2","is huge",
01250 "StTrack" ,"mLength" ,"3","3","is too small",
01251 "StTrack" ,"mLength" ,"3","4","contradicts to In/Out distance",
01252 "StTrack" ,"mLength" ,"3","5","helix out of Zmax",
01253 "StTrack" ,"mLength" ,"3","6","helix out of Rmax",
01254
01255 "StTrack" ,"mGeometry" ,"4","2","iz zero",
01256 "StTrack" ,"mGeometry" ,"4","0","StTrackGeometry",
01257
01258 "StTrack" ,"mOuterGeometry" ,"5","2","iz zero",
01259 "StTrack" ,"mOuterGeometry" ,"5","0","StTrackGeometry",
01260
01261 "StTrack" ,"mDetectorInfo" ,"6","2","iz zero",
01262 "StTrack" ,"mDetectorInfo" ,"6","0","StTrackDetectorInfo",
01263
01264 "StTrackGeometry" ,"Helix" ,"1","0","StPhysicalHelixD",
01265 "StTrackGeometry" ,"Helix" ,"1","2","out of zMax",
01266 "StTrackGeometry" ,"Helix" ,"1","3","out of rMax",
01267
01268 "StTrackDetectorInfo" ,"mFirstPoint" ,"1","0","StThreeVectorF",
01269 "StTrackDetectorInfo" ,"mFirstPoint" ,"1","2","out of zMax",
01270 "StTrackDetectorInfo" ,"mFirstPoint" ,"1","3","out of rMax",
01271
01272 "StTrackDetectorInfo" ,"mLastPoint" ,"2","0","StThreeVectorF",
01273 "StTrackDetectorInfo" ,"mLastPoint" ,"2","2","out of zMax",
01274 "StTrackDetectorInfo" ,"mFLastPoint" ,"2","3","out of rMax",
01275
01276
01277 "StPhysicalHelixD" ,"mDipAngle" ,"1","1","is NaN",
01278 "StPhysicalHelixD" ,"mDipAngle" ,"1","2","> Py/2",
01279 "StPhysicalHelixD" ,"mDipAngle" ,"1","3","== Py/2",
01280
01281 "StPhysicalHelixD" ,"mCurvature" ,"2","1","is NaN",
01282 "StPhysicalHelixD" ,"mCurvature" ,"2","2","too big",
01283 "StPhysicalHelixD" ,"mCurvature" ,"2","3","is Negaive",
01284
01285 "StPhysicalHelixD" ,"mOrigin" ,"3","0","StThreeVectorD",
01286 "StPhysicalHelixD" ,"mH" ,"4","2","!= 1 or -1",
01287
01288 "StThreeVectorD" ,"mX1" ,"1","1","is NaN",
01289 "StThreeVectorD" ,"mX1" ,"1","2","too big",
01290 "StThreeVectorD" ,"mX2" ,"2","1","is NaN",
01291 "StThreeVectorD" ,"mX2" ,"2","2","too big",
01292 "StThreeVectorD" ,"mX3" ,"3","1","is NaN",
01293 "StThreeVectorD" ,"mX3" ,"3","2","too big",
01294
01295
01296 "StThreeVectorF" ,"mX1" ,"1","1","is NaN",
01297 "StThreeVectorF" ,"mX1" ,"1","2","too big",
01298 "StThreeVectorF" ,"mX2" ,"2","1","is NaN",
01299 "StThreeVectorF" ,"mX2" ,"2","2","too big",
01300 "StThreeVectorF" ,"mX3" ,"3","1","is NaN",
01301 "StThreeVectorF" ,"mX3" ,"3","2","too big",
01302 0};
01303 TString ts;
01304
01305 int jmm = ierr%10;
01306 int jrr = (ierr/10)%10;
01307 for (const char **jt=TabErr;*jt;jt+=5) {
01308 if (strcmp(klass,*jt) ) continue;
01309 if (atoi(jt[2]) != jmm) continue;
01310 if (atoi(jt[3]) != jrr) continue;
01311 ts+=jt[1];
01312 if (jrr) { ts+=": "; ts+=jt[4];}
01313 else { ts+="."; ts+=Say(ierr/100,jt[4]);}
01314 return ts;
01315 }
01316 ts="***Unknown***";
01317 return ts;
01318 }
