00001 #include "TrackSegmentAtNu.h"
00002
00003 #include <cmath>
00004
00005 ClassImp(TrackSegmentAtNu)
00006
00007
00008
00009 TrackSegmentAtNu::TrackSegmentAtNu()
00010 {
00011
00012 }
00013
00014 TrackSegmentAtNu::TrackSegmentAtNu(ClusterAtNu* clrm, ClusterAtNu* clr0, ClusterAtNu* clrp) :
00015 fClrList(0),
00016 fBegTriList(0), fEndTriList(0),
00017 fBegAssocTriList(0), fEndAssocTriList(0),
00018 fBegSegList(0), fEndSegList(0),
00019 fSegment(0),
00020 fTrkFlag(0),fTmpTrkFlag(0),fNPlanes(0)
00021 {
00022
00023 this->AddCluster(clrm); this->AddCluster(clr0); this->AddCluster(clrp);
00024 }
00025
00026 TrackSegmentAtNu::TrackSegmentAtNu(const TrackSegmentAtNu& rhs) :
00027 ObjTrackAtNu(rhs),
00028 fClrList(rhs.fClrList),
00029 fBegTriList(rhs.fBegTriList),
00030 fEndTriList(rhs.fEndTriList),
00031 fBegAssocTriList(rhs.fBegAssocTriList),
00032 fEndAssocTriList(rhs.fEndAssocTriList),
00033 fBegSegList(rhs.fBegSegList),
00034 fEndSegList(rhs.fEndSegList),
00035 fSegment(rhs.fSegment),
00036 fTrkFlag(rhs.fTrkFlag),
00037 fTmpTrkFlag(rhs.fTmpTrkFlag),
00038 fNPlanes(rhs.fNPlanes)
00039 {
00040
00041 }
00042
00043 TrackSegmentAtNu::~TrackSegmentAtNu()
00044 {
00045 if(fClrList) delete fClrList;
00046 if(fBegTriList) delete fBegTriList;
00047 if(fEndTriList) delete fEndTriList;
00048 if(fBegAssocTriList) delete fBegAssocTriList;
00049 if(fEndAssocTriList) delete fEndAssocTriList;
00050 if(fBegSegList) delete fBegSegList;
00051 if(fEndSegList) delete fEndSegList;
00052 }
00053
00054 void TrackSegmentAtNu::AddSegment(TrackSegmentAtNu* segment)
00055 {
00056 for(Int_t i=0;i<1+segment->GetClrLast();i++){
00057 ClusterAtNu* clr = (ClusterAtNu*)(segment->GetClrAt(i));
00058 this->AddCluster(clr);
00059 }
00060 }
00061
00062 void TrackSegmentAtNu::AddCluster(ClusterAtNu* cluster)
00063 {
00064 if(!fClrList) fClrList = new TObjArray();
00065 if(!this->ContainsClr(cluster)){
00066 fClrList->Add(cluster);
00067 for(Int_t i=0;i<1+cluster->GetHitLast();i++){
00068 HitAtNu* hit = (HitAtNu*)(cluster->GetHitAt(i));
00069 this->AddHitToTrack(hit);
00070 }
00071 }
00072 return;
00073 }
00074
00075 Bool_t TrackSegmentAtNu::ContainsClr(ClusterAtNu* cluster)
00076 {
00077 for(Int_t i=0;i<1+fClrList->GetLast();i++){
00078 if(cluster==(ClusterAtNu*)(fClrList->At(i))) return 1;
00079 }
00080 return 0;
00081 }
00082
00083 ClusterAtNu* TrackSegmentAtNu::GetClrAt(Int_t i)
00084 {
00085 if(fClrList) return (ClusterAtNu*)(fClrList->At(i)); else return 0;
00086 }
00087
00088 Int_t TrackSegmentAtNu::GetClrLast()
00089 {
00090 if(fClrList) return fClrList->GetLast(); else return -1;
00091 }
00092
00093 void TrackSegmentAtNu::AddTriToBeg(TrackSegmentAtNu* segment)
00094 {
00095 if(!fBegTriList) fBegTriList = new TObjArray();
00096 fBegTriList->Add(segment);
00097 return;
00098 }
00099
00100 TrackSegmentAtNu* TrackSegmentAtNu::GetTriBegAt(Int_t i)
00101 {
00102 if(fBegTriList) return (TrackSegmentAtNu*)(fBegTriList->At(i)); else return 0;
00103 }
00104
00105 Int_t TrackSegmentAtNu::GetTriBegLast()
00106 {
00107 if(fBegTriList) return fBegTriList->GetLast(); else return -1;
00108 }
00109
00110 void TrackSegmentAtNu::AddTriToEnd(TrackSegmentAtNu* segment)
00111 {
00112 if(!fEndTriList) fEndTriList = new TObjArray();
00113 fEndTriList->Add(segment);
00114 return;
00115 }
00116
00117 TrackSegmentAtNu* TrackSegmentAtNu::GetTriEndAt(Int_t i)
00118 {
00119 if(fEndTriList) return (TrackSegmentAtNu*)(fEndTriList->At(i)); else return 0;
00120 }
00121
00122 Int_t TrackSegmentAtNu::GetTriEndLast()
00123 {
00124 if(fEndTriList) return fEndTriList->GetLast(); else return -1;
00125 }
00126
00127 void TrackSegmentAtNu::AddAssocTriToBeg(TrackSegmentAtNu* segment)
00128 {
00129 if(!fBegAssocTriList) fBegAssocTriList = new TObjArray();
00130 fBegAssocTriList->Add(segment);
00131 return;
00132 }
00133
00134 TrackSegmentAtNu* TrackSegmentAtNu::GetAssocTriBegAt(Int_t i)
00135 {
00136 if(fBegAssocTriList) return (TrackSegmentAtNu*)(fBegAssocTriList->At(i)); else return 0;
00137 }
00138
00139 Int_t TrackSegmentAtNu::GetAssocTriBegLast()
00140 {
00141 if(fBegAssocTriList) return fBegAssocTriList->GetLast(); else return -1;
00142 }
00143
00144 void TrackSegmentAtNu::AddAssocTriToEnd(TrackSegmentAtNu* segment)
00145 {
00146 if(!fEndAssocTriList) fEndAssocTriList = new TObjArray();
00147 fEndAssocTriList->Add(segment);
00148 return;
00149 }
00150
00151 TrackSegmentAtNu* TrackSegmentAtNu::GetAssocTriEndAt(Int_t i)
00152 {
00153 if(fEndAssocTriList) return (TrackSegmentAtNu*)(fEndAssocTriList->At(i)); else return 0;
00154 }
00155
00156 Int_t TrackSegmentAtNu::GetAssocTriEndLast()
00157 {
00158 if(fEndAssocTriList) return fEndAssocTriList->GetLast(); else return -1;
00159 }
00160
00161 Int_t TrackSegmentAtNu::IsBegAssoc(TrackSegmentAtNu* segment)
00162 {
00163 Int_t i,assoc=0;
00164 Int_t begflag=0,endflag=0;
00165 if(segment->GetEndPlane()>=this->GetBegPlane()){
00166 begflag=1; endflag=1;
00167 for(i=0;i<1+segment->GetClrLast();i++){
00168 ClusterAtNu* clr = (ClusterAtNu*)(segment->GetClrAt(i));
00169 if(clr->GetPlane()>=this->GetBegPlane()){ if(!(this->ContainsClr(clr))) endflag=0; }
00170 }
00171 for(i=0;i<1+this->GetClrLast();i++){
00172 ClusterAtNu* clr = (ClusterAtNu*)(this->GetClrAt(i));
00173 if(clr->GetPlane()<=segment->GetEndPlane()){ if(!(segment->ContainsClr(clr))) begflag=0; }
00174 }
00175 if(segment->GetEndPlane()>this->GetBegPlane()){ if(begflag && endflag) assoc=1; }
00176 }
00177
00178 if(segment->GetEndPlane()<=this->GetBegPlane()){
00179 Int_t idb=0,idb0=0;
00180 Int_t bpln=this->GetEndPlane()+1;
00181 for(i=0;i<1+this->GetClrLast();i++){
00182 ClusterAtNu* clr = (ClusterAtNu*)(this->GetClrAt(i));
00183 if(clr->GetPlane()<bpln && clr->GetPlane()>this->GetBegPlane()+2){
00184 bpln=clr->GetPlane(); idb=i;
00185 }
00186 else{ if(clr->GetPlane()==this->GetBegPlane()) idb0=i; }
00187 }
00188 ClusterAtNu* clrb = (ClusterAtNu*)(this->GetClrAt(idb));
00189 ClusterAtNu* clrb0 = (ClusterAtNu*)(this->GetClrAt(idb0));
00190
00191 Int_t ide=0,ide0=0;
00192 Int_t epln=segment->GetBegPlane()-1;
00193 for(i=0;i<1+segment->GetClrLast();i++){
00194 ClusterAtNu* clr = (ClusterAtNu*)(segment->GetClrAt(i));
00195 if(clr->GetPlane()>epln && clr->GetPlane()<segment->GetEndPlane()-2){
00196 epln=clr->GetPlane(); ide=i;
00197 }
00198 else{ if(clr->GetPlane()==segment->GetEndPlane()) ide0=i; }
00199 }
00200 ClusterAtNu* clre = (ClusterAtNu*)(segment->GetClrAt(ide));
00201 ClusterAtNu* clre0 = (ClusterAtNu*)(segment->GetClrAt(ide0));
00202
00203 if( ( clrb->GetEndStrip()-clrb0->GetBegStrip()>-2 && clrb0->GetEndStrip()-clrb->GetBegStrip()>-2
00204 && clre->GetEndStrip()-clre0->GetBegStrip()>-2 && clre0->GetEndStrip()-clre->GetBegStrip()>-2 )
00205 || ( ( clrb->GetBegStrip()-clrb0->GetBegStrip()>-2 && clre0->GetBegStrip()-clre->GetBegStrip()>-2 )
00206 || ( clrb->GetEndStrip()-clrb0->GetEndStrip()>-2 && clre0->GetEndStrip()-clre->GetEndStrip()>-2 ) )
00207 || ( ( clrb->GetBegStrip()-clrb0->GetBegStrip()<2 && clre0->GetBegStrip()-clre->GetBegStrip()<2 )
00208 || ( clrb->GetEndStrip()-clrb0->GetEndStrip()<2 && clre0->GetEndStrip()-clre->GetEndStrip()<2 ) ) ){
00209 if(this->GetBegPlane()==segment->GetEndPlane()){ if(begflag && endflag) assoc=1; }
00210 else{
00211 Double_t z1,z2,t1,t2,dt1,dt2,dir1,dir2,dirtmp,win;
00212 z1=segment->GetEndVtxZ(); z2=this->GetBegVtxZ();
00213 t1=segment->GetEndVtxT(); t2=this->GetBegVtxT();
00214 dir1=segment->GetEndDir(); dir2=this->GetBegDir();
00215 dt1=t2-t1-dir1*(z2-z1); dt2=t2-t1-dir2*(z2-z1); win=0.1+0.35*(z2-z1);
00216 dirtmp=(1+dir1*dir2)/(sqrt(1+dir1*dir1)*sqrt(1+dir2*dir2));
00217 if( ( (dt1>-win&&dt1<win)||(dt2>-win&&dt2<win) ) && ( dirtmp>0.65 ) ) assoc=1;
00218 }
00219 }
00220 }
00221 return assoc;
00222 }
00223
00224 Int_t TrackSegmentAtNu::IsEndAssoc(TrackSegmentAtNu* segment)
00225 {
00226 Int_t i,assoc=0;
00227 Int_t begflag=0,endflag=0;
00228 if(segment->GetBegPlane()<=this->GetEndPlane()){
00229 begflag=1; endflag=1;
00230 for(i=0;i<1+segment->GetClrLast();i++){
00231 ClusterAtNu* clr = (ClusterAtNu*)(segment->GetClrAt(i));
00232 if(clr->GetPlane()<=this->GetEndPlane()){ if(!(this->ContainsClr(clr))) begflag=0; }
00233 }
00234 for(i=0;i<1+this->GetClrLast();i++){
00235 ClusterAtNu* clr = (ClusterAtNu*)(this->GetClrAt(i));
00236 if(clr->GetPlane()>=segment->GetBegPlane()){ if(!(segment->ContainsClr(clr))) endflag=0; }
00237 }
00238 if(segment->GetBegPlane()<this->GetEndPlane()){ if(begflag && endflag) assoc=1; }
00239 }
00240
00241 if(segment->GetBegPlane()>=this->GetEndPlane()){
00242 Int_t idb=0,idb0=0;
00243 Int_t bpln=segment->GetEndPlane()+1;
00244 for(i=0;i<1+segment->GetClrLast();i++){
00245 ClusterAtNu* clr = (ClusterAtNu*)(segment->GetClrAt(i));
00246 if(clr->GetPlane()<bpln && clr->GetPlane()>segment->GetBegPlane()+2){
00247 bpln=clr->GetPlane(); idb=i;
00248 }
00249 else{ if(clr->GetPlane()==segment->GetBegPlane()) idb0=i; }
00250 }
00251 ClusterAtNu* clrb = (ClusterAtNu*)(segment->GetClrAt(idb));
00252 ClusterAtNu* clrb0 = (ClusterAtNu*)(segment->GetClrAt(idb0));
00253
00254 Int_t ide=0,ide0=0;
00255 Int_t epln=this->GetBegPlane()-1;
00256 for(i=0;i<1+this->GetClrLast();i++){
00257 ClusterAtNu* clr = (ClusterAtNu*)(this->GetClrAt(i));
00258 if(clr->GetPlane()>epln && clr->GetPlane()<this->GetEndPlane()-2){
00259 epln=clr->GetPlane(); ide=i;
00260 }
00261 else{ if(clr->GetPlane()==this->GetEndPlane()) ide0=i; }
00262 }
00263 ClusterAtNu* clre = (ClusterAtNu*)(this->GetClrAt(ide));
00264 ClusterAtNu* clre0 = (ClusterAtNu*)(this->GetClrAt(ide0));
00265
00266 if( ( clrb->GetEndStrip()-clrb0->GetBegStrip()>-2 && clrb0->GetEndStrip()-clrb->GetBegStrip()>-2
00267 && clre->GetEndStrip()-clre0->GetBegStrip()>-2 && clre0->GetEndStrip()-clre->GetBegStrip()>-2 )
00268 || ( ( clrb->GetBegStrip()-clrb0->GetBegStrip()>-2 && clre0->GetBegStrip()-clre->GetBegStrip()>-2 )
00269 || ( clrb->GetEndStrip()-clrb0->GetEndStrip()>-2 && clre0->GetEndStrip()-clre->GetEndStrip()>-2 ) )
00270 || ( ( clrb->GetBegStrip()-clrb0->GetBegStrip()<2 && clre0->GetBegStrip()-clre->GetBegStrip()<2 )
00271 || ( clrb->GetEndStrip()-clrb0->GetEndStrip()<2 && clre0->GetEndStrip()-clre->GetEndStrip()<2 ) ) ){
00272 if(segment->GetBegPlane()==this->GetEndPlane()){ if(begflag && endflag) assoc=1; }
00273 else{
00274 Double_t z1,z2,t1,t2,dt1,dt2,dir1,dir2,dirtmp,win;
00275 z1=this->GetEndVtxZ(); z2=segment->GetBegVtxZ();
00276 t1=this->GetEndVtxT(); t2=segment->GetBegVtxT();
00277 dir1=this->GetEndDir(); dir2=segment->GetBegDir();
00278 dt1=t2-t1-dir1*(z2-z1); dt2=t2-t1-dir2*(z2-z1); win=0.1+0.35*(z2-z1);
00279 dirtmp=(1+dir1*dir2)/(sqrt(1+dir1*dir1)*sqrt(1+dir2*dir2));
00280 if( ( (dt1>-win&&dt1<win)||(dt2>-win&&dt2<win) ) && ( dirtmp>0.65 ) ) assoc=1;
00281 }
00282 }
00283 }
00284 return assoc;
00285 }
00286
00287 void TrackSegmentAtNu::AddSegToBeg(TrackSegmentAtNu* segment)
00288 {
00289 if(!fBegSegList) fBegSegList = new TObjArray();
00290 fBegSegList->Add(segment);
00291 return;
00292 }
00293
00294 TrackSegmentAtNu* TrackSegmentAtNu::GetSegBegAt(Int_t i)
00295 {
00296 if(fBegSegList) return (TrackSegmentAtNu*)(fBegSegList->At(i)); else return 0;
00297 }
00298
00299 Int_t TrackSegmentAtNu::GetSegBegLast()
00300 {
00301 if(fBegSegList) return fBegSegList->GetLast(); else return -1;
00302 }
00303
00304 void TrackSegmentAtNu::AddSegToEnd(TrackSegmentAtNu* segment)
00305 {
00306 if(!fEndSegList) fEndSegList = new TObjArray();
00307 fEndSegList->Add(segment);
00308 return;
00309 }
00310
00311 TrackSegmentAtNu* TrackSegmentAtNu::GetSegEndAt(Int_t i)
00312 {
00313 if(fEndSegList) return (TrackSegmentAtNu*)(fEndSegList->At(i)); else return 0;
00314 }
00315
00316 Int_t TrackSegmentAtNu::GetSegEndLast()
00317 {
00318 if(fEndSegList) return fEndSegList->GetLast(); else return -1;
00319 }
00320
00321 void TrackSegmentAtNu::SetSegment(TrackSegmentAtNu* segment)
00322 {
00323 fSegment = segment;
00324 }
00325
00326 TrackSegmentAtNu* TrackSegmentAtNu::GetSegment()
00327 {
00328 return (TrackSegmentAtNu*)(fSegment);
00329 }
00330
00331 void TrackSegmentAtNu::SetTrkFlag(Int_t flag)
00332 {
00333 fTrkFlag=flag;
00334 }
00335
00336 Int_t TrackSegmentAtNu::GetTrkFlag()
00337 {
00338 return fTrkFlag;
00339 }
00340
00341 void TrackSegmentAtNu::SetTmpTrkFlag(Int_t flag)
00342 {
00343 fTmpTrkFlag=flag;
00344 }
00345
00346 Int_t TrackSegmentAtNu::GetTmpTrkFlag()
00347 {
00348 return fTmpTrkFlag;
00349 }
00350
00351 void TrackSegmentAtNu::SetNPlanes(Int_t plns)
00352 {
00353 fNPlanes=plns;
00354 }
00355
00356 Int_t TrackSegmentAtNu::GetNPlanes()
00357 {
00358 return fNPlanes;
00359 }
00360
00361
00362 Double_t TrackSegmentAtNu::GetScore(TObjArray* objbeg, TObjArray* objend)
00363 {
00364
00365 Int_t k,k1;
00366 Int_t bpln2,bpln1,epln1,epln2;
00367 Int_t pln,npln;
00368 TObjArray tmparr;
00369
00370 bpln1=this->GetBegPlane();
00371 if(objbeg){
00372 TrackSegmentAtNu* segb = (TrackSegmentAtNu*)(objbeg->Last());
00373 bpln2=segb->GetBegPlane();
00374 }
00375 else{
00376 bpln2=this->GetBegPlane();
00377 }
00378
00379 epln1=this->GetEndPlane();
00380 if(objend){
00381 TrackSegmentAtNu* sege = (TrackSegmentAtNu*)(objend->Last());
00382 epln2=sege->GetEndPlane();
00383 }
00384 else{
00385 epln2=this->GetEndPlane();
00386 }
00387
00388
00389 for(k=0;k<1+fClrList->GetLast();k++){
00390 ClusterAtNu* clr = (ClusterAtNu*)(fClrList->At(k));
00391 tmparr.Add(clr);
00392 }
00393
00394 if(objbeg){
00395 for(k=0;k<1+objbeg->GetLast();k++){
00396 TrackSegmentAtNu* segb = (TrackSegmentAtNu*)(objbeg->At(k));
00397 for(k1=0;k1<1+segb->GetClrLast();k1++){
00398 ClusterAtNu* clr = (ClusterAtNu*)(segb->GetClrAt(k1));
00399 if(!tmparr.Contains(clr)) tmparr.Add(clr);
00400 }
00401 }
00402 }
00403
00404 if(objend){
00405 for(k=0;k<1+objend->GetLast();k++){
00406 TrackSegmentAtNu* sege = (TrackSegmentAtNu*)(objend->At(k));
00407 for(k1=0;k1<1+sege->GetClrLast();k1++){
00408 ClusterAtNu* clr = (ClusterAtNu*)(sege->GetClrAt(k1));
00409 if(!tmparr.Contains(clr)) tmparr.Add(clr);
00410 }
00411 }
00412 }
00413
00414 npln=1+(epln2-bpln2)/2;
00415
00416 Double_t* T = new Double_t[npln];
00417 Double_t* Z = new Double_t[npln];
00418 Double_t* W = new Double_t[npln];
00419
00420 for(k=0;k<npln;k++){
00421 T[k]=0.0; Z[k]=0.0; W[k]=0.0;
00422 }
00423
00424 Int_t km,kp;
00425 Double_t m,c;
00426 Double_t dt2,sn;
00427 Double_t dscr,scr,dtot,tot,dplnscr,plnscr;
00428 Double_t sw,swz,swt,swzt,swzz;
00429 for(k=0;k<1+tmparr.GetLast();k++){
00430 ClusterAtNu* clr = (ClusterAtNu*)(tmparr.At(k));
00431 pln=(clr->GetPlane()-bpln2)/2;
00432 sw=0.0; swz=0.0; swt=0.0;
00433 for(k1=0;k1<1+clr->GetHitLast();k1++){
00434 HitAtNu* hit = (HitAtNu*)(clr->GetHitAt(k1));
00435 swz+=hit->GetCharge()*hit->GetZPos();
00436 swt+=hit->GetCharge()*hit->GetTPos();
00437 sw+=hit->GetCharge();
00438 }
00439 if(sw>0.0){
00440 Z[pln]=swz/sw; T[pln]=swt/sw; W[pln]=1.0;
00441 if( pln+1>(bpln1-bpln2)/2
00442 && pln-1<(epln1-bpln2)/2 ) W[pln]=5.0*W[pln]; else W[pln]=0.5*W[pln];
00443 }
00444 }
00445
00446 scr=0.0; tot=0.0; plnscr=0.0;
00447 for(k=0;k<npln;k++){
00448 dscr=0.01; dtot=0.01; dplnscr=0.01;
00449 if(W[k]>0.0){
00450 swz=0.0; swt=0.0; swzz=0.0; swzt=0.0; sw=0.0; sn=0.0;
00451 km=k-4; kp=k+4;
00452 if(W[k]>1.0){
00453 if( km<(bpln1-bpln2)/2 ) km=(bpln1-bpln2)/2;
00454 if( kp>(epln1-bpln2)/2 ) kp=(epln1-bpln2)/2;
00455 }
00456 if(W[k]<1.0){
00457 if( km>(epln1-bpln2)/2 ) km=-4+(epln1-bpln2)/2;
00458 if( kp<(bpln1-bpln2)/2 ) kp=+4+(bpln1-bpln2)/2;
00459 }
00460 if(km<0) km=0; if(kp>npln-1) kp=npln-1;
00461 for(k1=km;k1<kp+1;k1++){
00462 if(W[k1]>0.0){
00463 swz+=W[k1]*Z[k1]; swt+=W[k1]*T[k1];
00464 swzz+=W[k1]*Z[k1]*Z[k1]; swzt+=W[k1]*Z[k1]*T[k1];
00465 sw+=W[k1]; sn+=1.0;
00466 }
00467 }
00468 if(sn>2.0){
00469 m=(sw*swzt-swz*swt)/(sw*swzz-swz*swz);
00470 c=(swt*swzz-swz*swzt)/(sw*swzz-swz*swz);
00471 dt2=(T[k]-(m*Z[k]+c))*(T[k]-(m*Z[k]+c));
00472 if(dt2<dscr) dscr=dt2; else dscr+=0.01*(1.0-dscr/dt2);
00473 }
00474 }
00475 scr+=dscr; tot+=dtot; plnscr+=dplnscr;
00476 }
00477
00478 scr = 1.0 + 4.0*plnscr - sqrt(scr/tot);
00479 if(scr<0.0) scr=0.0;
00480
00481 delete [] T;
00482 delete [] Z;
00483 delete [] W;
00484
00485 return scr;
00486
00487 }
00488
00489
00490
