FEDRA emulsion software from the OPERA Collaboration
EdbViewRec Class Reference

#include <EdbViewRec.h>

Inheritance diagram for EdbViewRec:
Collaboration diagram for EdbViewRec:

Public Member Functions

float CalculateSegmentChi2 (EdbSegment &seg, float sx, float sy, float sz)
 
int CheckFramesShift ()
 
int CheckSeedThres ()
 
float Chi2Seg (EdbSegment &s1, EdbSegment &s2)
 
 EdbViewRec ()
 
 EdbViewRec (EdbViewDef &vd)
 
int FillGrainsTree ()
 
int FindGrains (int option=0)
 
int FindSeeds ()
 
bool GoodSegment (EdbSegment &s, int wkey=0)
 
bool Init ()
 
void InitGrainsTree (const char *file="grains.root")
 
void InitR ()
 
int MergeSegments ()
 
int ReconstructGrains ()
 
int ReconstructSegments ()
 
int RefillSegment (EdbSegment &s)
 
int RefitSegments (int wkey=0)
 
void ResetClustersSeg ()
 
bool SaveToOutputView (EdbView &vout, int do_h=1, int do_c=2, int do_s=2, int do_tr=0, int do_f=2)
 
int SelectSegments ()
 
void SetAddGrainsToView (bool yesno)
 
void SetClThres (int mina, int maxa)
 
void SetNclGrLim (int mincl, int maxcl)
 
void SetNgrMax0 (Int_t ngr)
 
void SetNSeedMax (int nt, int th[])
 
void SetNSeedMax0 (int n)
 
void SetPrimary ()
 
void SetPulsThres (int minp, int maxp=500)
 
void SetRmax (float rmax)
 
void SetSeedsLim (int nseedslim=100000, int seedlim=48)
 
void SetSeedThres (int nt, int th[])
 
void SetSeedThres0 (Short_t mins)
 
void SetSigmaThres (float smin, float smax)
 
int SetStep (int sfrom, int sto, int step=1)
 
void SetThetaLim (float t)
 
bool SetView (EdbView *v)
 
float SPhiGr (float theta, float phi, float dz, float sx, float sy, float sz)
 
float SThetaGr (float theta, float phi, float dz, float sx, float sy, float sz)
 
virtual ~EdbViewRec ()
 
- Public Member Functions inherited from EdbViewDef
 EdbViewDef ()
 
void Print ()
 
void SetDef ()
 
 ~EdbViewDef ()
 

Static Public Member Functions

static int FitSegment (EdbSegment &s, int wkey=0)
 
static int FitSegmentToCl (EdbSegment &s, EdbCluster &c, int wkey=0)
 

Public Attributes

bool eAddGrainsToView
 
bool eCheckSeedThres
 
bool eDoGrainsProcessing
 
Float_t eGrainNbin
 
bool ePropagateToBase
 
- Public Attributes inherited from EdbViewDef
Float_t eClaSX
 
Float_t eClaSY
 
Float_t eClaSZ
 
Float_t eClaSZvar
 
Float_t eDZdead
 
Float_t eFogDens
 
Float_t eFogGrainArea
 
Float_t eGrainArea
 
Float_t eGrainSX
 
Float_t eGrainSY
 
Float_t eGrainSZ
 
Int_t eNframes
 
Float_t eTXopt
 
Float_t eTYopt
 
Float_t eX0
 
Float_t eX0opt
 
Float_t eXmax
 
Float_t eXmin
 
Float_t eY0
 
Float_t eY0opt
 
Float_t eYmax
 
Float_t eYmin
 
Float_t eZdead
 
Float_t eZmax
 
Float_t eZmin
 
Float_t eZxy
 

Private Attributes

TClonesArray * eCL
 [eNgrMax] array of grains represented as clusters More...
 
Int_t eClMaxA
 
Int_t eClMinA
 pointer to array of segments (output of tracking) More...
 
Float_t eDmax
 pointers to clusters More...
 
Float_t eDZmin
 
Float_t eFact
 
TClonesArray * eG
 cells with grains More...
 
TClonesArray * eGCla
 pointer to eView->GetSegments() or to eGSeg as the output for grain search More...
 
TTree * eGrainsTree
 
Short_t * ehX
 [enYtot] - pointers to the first x[iy] More...
 
int eNclGrMax
 
int eNclGrMin
 debug tree More...
 
Int_t eNgr
 
Int_t eNgrMax
 pointer to the input view currently in processing More...
 
TArrayI enP
 
Int_t enPtot
 
TArrayI eNseedMax
 
Int_t eNseedMax0
 
Int_t enSeedsLim
 
Int_t eNsegMax
 
Int_t enT
 
TArrayI enX
 
Int_t enXtot
 
TArrayI enY
 
Int_t enYtot
 
EdbCluster ** epC
 pointers to seeds list More...
 
Short_t *** epP
 [enT] - pointers to the first phi[it] More...
 
EdbCluster *** epS
 [enXtot] - phase histogram More...
 
Short_t **** epT
 
Int_t ePulsMax
 
Int_t ePulsMin
 
Short_t ** epY
 [enPtot] - pointers to the first y[ip] More...
 
TArrayF eR
 
Float_t eRmax
 
TClonesArray * eSA
 pointer to eView->GetClusters() or to eGCla as the input for tracking More...
 
Short_t eSeedLim
 
TArrayI eSeedThres
 
Int_t eSeedThres0
 
Float_t eSigmaMax
 
Float_t eSigmaMin
 
TArrayF esP
 
Int_t eStep
 
Int_t eStepFrom
 
Int_t eStepTo
 
TArrayF esX
 
TArrayF esY
 
TArrayF eTheta
 
Float_t eThetaLim
 
EdbViewCelleVC
 
EdbViewCell eVCC
 pointer to eVCC or eVCG More...
 
EdbViewCell eVCG
 cells with raw clusters More...
 
EdbVieweView
 
Float_t eZcenter
 

Constructor & Destructor Documentation

◆ EdbViewRec() [1/2]

EdbViewRec::EdbViewRec ( )
167{
168 SetDef();
169 SetPrimary();
170}
void SetDef()
Definition: EdbViewDef.cxx:29
void SetPrimary()
Definition: EdbViewRec.cxx:198

◆ EdbViewRec() [2/2]

EdbViewRec::EdbViewRec ( EdbViewDef vd)
inline
void vd(int trmin=2, float amin=.0)
Definition: check_vertex.C:217

◆ ~EdbViewRec()

EdbViewRec::~EdbViewRec ( )
virtual
174{
175 printf("EdbViewRec::~EdbViewRec()\n");
176 if(eGCla) { delete eGCla; eGCla=0; }
177 printf("1\n");
178 if(eGrainsTree) { delete eGrainsTree; eGrainsTree=0; }
179 printf("2\n");
180
181 if(epT) { delete[] epT; epT=0; }
182 if(epP) { delete[] epP; epP=0; }
183 if(epY) { delete[] epY; epY=0; }
184 if(ehX) { delete[] ehX; ehX=0; }
185 if(epS) { delete[] epS; epS=0; }
186 if(epC) { delete[] epC; epC=0; }
187 printf("3\n");
188
189 if(eNsegMax>0) {
190 if(eSA!=eG) if(eSA) {eSA->Clear(); delete eSA; eSA=0; }
191 printf("4\n");
192 if(eG) { eG->Clear(); delete eG; eG=0; }
193 }
194 printf("5\n");
195}
Short_t ** epY
[enPtot] - pointers to the first y[ip]
Definition: EdbViewRec.h:177
TClonesArray * eGCla
pointer to eView->GetSegments() or to eGSeg as the output for grain search
Definition: EdbViewRec.h:124
TClonesArray * eSA
pointer to eView->GetClusters() or to eGCla as the input for tracking
Definition: EdbViewRec.h:127
TTree * eGrainsTree
Definition: EdbViewRec.h:134
Int_t eNsegMax
Definition: EdbViewRec.h:115
Short_t *** epP
[enT] - pointers to the first phi[it]
Definition: EdbViewRec.h:176
TClonesArray * eG
cells with grains
Definition: EdbViewRec.h:122
EdbCluster ** epC
pointers to seeds list
Definition: EdbViewRec.h:181
Short_t * ehX
[enYtot] - pointers to the first x[iy]
Definition: EdbViewRec.h:178
EdbCluster *** epS
[enXtot] - phase histogram
Definition: EdbViewRec.h:180
Short_t **** epT
Definition: EdbViewRec.h:175

Member Function Documentation

◆ CalculateSegmentChi2()

float EdbViewRec::CalculateSegmentChi2 ( EdbSegment s,
float  sx,
float  sy,
float  sz 
)

816{
817 //assumed that clusters are attached to segments
818 TObjArray *clusters = s.GetElements();
819 if(!clusters) return 0;
820 int ncl = clusters->GetLast()+1;
821 if(ncl<=0) return 0;
822
823 // segment line parametrized as 2 points
824 float xyz1[3] = { s.GetX0() /sx,
825 s.GetY0() /sy,
826 s.GetZ0() /sz };
827 float xyz2[3] = { (s.GetX0() + s.GetDz()*s.GetTx()) /sx,
828 (s.GetY0() + s.GetDz()*s.GetTy()) /sy,
829 (s.GetZ0() + s.GetDz()) /sz };
830
831 bool inside=true;
832 double d, chi2=0;
833 EdbCluster *cl=0;
834 float xyz[3];
835 for(int i=0; i<ncl; i++ ) {
836 cl = (EdbCluster*)clusters->At(i);
837 xyz[0] = cl->GetX()/sx;
838 xyz[1] = cl->GetY()/sy;
839 xyz[2] = cl->GetZ()/sz;
840 d = EdbMath::DistancePointLine3(xyz,xyz1,xyz2, inside);
841 chi2 += d*d;
842 }
843 return TMath::Sqrt(chi2/ncl);
844}
void d()
Definition: RecDispEX.C:381
Definition: EdbCluster.h:19
Float_t GetX() const
Definition: EdbCluster.h:51
Float_t GetY() const
Definition: EdbCluster.h:52
Float_t GetZ() const
Definition: EdbCluster.h:53
static double DistancePointLine3(float Point[3], float LineStart[3], float LineEnd[3], bool inside)
Definition: EdbMath.cxx:61
s
Definition: check_shower.C:55
Float_t chi2
Definition: testBGReduction_By_ANN.C:14

◆ CheckFramesShift()

int EdbViewRec::CheckFramesShift ( )
679{
680 if( eNgr<1 ) return 0;
681 int nfr = eView->GetNframes();
682 if( nfr<2 ) return 0;
683
684 TArrayD Xshift(nfr),Yshift(nfr),Zshift(nfr);
685 TArrayI Nshift(nfr);
686 EdbSegment *s=0;
687 EdbCluster *cl0=0,*cl1=0;
688 TObjArray *arr=0;
689 int ncl=0;
690 int ifr=0;
691
692 //TNtuple *fshift = new TNtuple("fshift","frames shift","x0:y0:z0:x1:y1:z1:ifr:ncl:gr");
693 for(int i=0; i<eNgr; i++) {
694 s = (EdbSegment*)(eG->UncheckedAt(i));
695 arr = s->GetElements();
696 ncl = arr->GetEntriesFast();
697 if(ncl<2) continue;
698 for(int j=1; j<ncl; j++) {
699 cl0 = (EdbCluster*)s->GetElements()->UncheckedAt(j-1);
700 cl1 = (EdbCluster*)s->GetElements()->UncheckedAt(j);
701 ifr = cl1->eFrame;
702 Nshift[ifr]++;
703 Xshift[ifr] += (cl1->eX - cl0->eX);
704 Yshift[ifr] += (cl1->eY - cl0->eY);
705 Zshift[ifr] += (cl1->eZ - cl0->eZ);
706 // fshift->Fill( cl0->eX,cl0->eY,cl0->eZ, cl1->eX,cl1->eY,cl1->eZ,1.*cl1->eFrame,1.*ncl, 1.*cl0->eSegment );
707 }
708 }
709
710 for(int i=1; i<nfr; i++) {
711 Xshift[i] /= Nshift[i];
712 Yshift[i] /= Nshift[i];
713 Zshift[i] /= Nshift[i];
714 printf("Frame shift(%2d): %+6.4f %+6.4f %f\t by %d shadows\n",i, Xshift[i],Yshift[i],Zshift[i],Nshift[i] );
715 }
716
717 printf("------------------------------------------------------\n");
718 printf("Total shift : %+4.2f %+4.2f \t by %d shadows; dz=%6.2f\n",
719 Xshift.GetSum(),Yshift.GetSum(),(int)Nshift.GetSum(),Zshift.GetSum() );
720 printf("Mean angle : %+4.3f %+4.3f \n\n",
721 Xshift.GetSum()/Zshift.GetSum(),
722 Yshift.GetSum()/Zshift.GetSum() );
723
724 return (int)Nshift.GetSum();
725}
Float_t eZ
cluster coordinates in pixels(?)
Definition: EdbCluster.h:25
Float_t eY
cluster coordinates in pixels(?)
Definition: EdbCluster.h:24
Float_t eX
cluster coordinates in pixels(?)
Definition: EdbCluster.h:23
segment of the track
Definition: EdbSegment.h:63
EdbView * eView
Definition: EdbViewRec.h:112
Int_t eNgr
Definition: EdbViewRec.h:133
Int_t GetNframes() const
Definition: EdbView.h:206

◆ CheckSeedThres()

int EdbViewRec::CheckSeedThres ( )
1133{
1134 Short_t *p0,*pnext, *plast = ehX+enXtot;
1135 Int_t nfill;
1136 Int_t ntot = 0;
1137
1138 for(int i=0; i<enT; i++) {
1140 }
1141
1142 for(int it=0; it<enT; it++) {
1143 TArrayL spt(256); //phase occupancy histogram
1144 p0 = epT[it][0][0];
1145 if(it==enT-1) pnext = plast;
1146 else pnext = epT[it+1][0][0];
1147 ntot = (int)((pnext-p0)/eFact);
1148 nfill=0;
1149 while(p0<pnext) {
1150 spt[*p0]++;
1151 if(*p0>0) nfill++;
1152 p0++;
1153 }
1154 printf("\nit = %3d occup: %d / %d = %f %% \n",it, nfill, ntot, 100.*nfill/ntot );
1155
1156 int thres=0, sum=0, sumthres=0;
1157 for (int i=255; i>0; i--) {
1158 if(spt[i]<1) continue;
1159 printf("%3d %ld\n",i,spt[i]);
1160 if(!thres) sumthres=sum;
1161 sum += spt[i];
1162 if(!thres&&sum>eNseedMax[it]) thres = i+1;
1163 }
1164 if(thres>eSeedThres[it]) eSeedThres[it]=thres;
1165 printf("sumthres(%d) = %d eSeedThres[%d]= %d\n\n",thres,sumthres, it, eSeedThres[it]);
1166 }
1167
1168
1169 //TArrayL sp(256); //phase occupancy histogram
1170 //for (int i=0; i<enXtot; i++)
1171 //if(ehX[i]>0) sp[ehX[i]]++;
1172 //printf("\n total occup: %d / %d = %f %% \n", (int)(sp.GetSum()), enXtot, 100.*eFact*sp.GetSum()/enXtot );
1173
1174// TNtuple *ntpix = new TNtuple("ntpix","pixels","i:pix:ix:iy:t:p");
1175// float p,t;
1176// int iphi,ix,iy;
1177// int jxy=0;
1178// for (int i=0; i<eNtot; i++) {
1179// if(ePhist[i]<eSeedThres) continue;
1180// iphi = i%eNphiTot;
1181// t = eTmask[iphi];
1182// p = ePmask[iphi];
1183// jxy = i/eNphiTot;
1184// iy = jxy/eVC->eNx;
1185// ix = jxy%eVC->eNx;
1186// ntpix->Fill(i,ePhist[i],ix,iy,t,p);
1187// }
1188// return thres;
1189
1190 return 0;
1191}
TArrayI eNseedMax
Definition: EdbViewRec.h:153
Int_t enT
Definition: EdbViewRec.h:161
Float_t eFact
Definition: EdbViewRec.h:183
TArrayI eSeedThres
Definition: EdbViewRec.h:148
Int_t eSeedThres0
Definition: EdbViewRec.h:147
Int_t enXtot
Definition: EdbViewRec.h:167

◆ Chi2Seg()

float EdbViewRec::Chi2Seg ( EdbSegment s1,
EdbSegment s2 
)
1127{
1128 return 1.;
1129}

◆ FillGrainsTree()

int EdbViewRec::FillGrainsTree ( )

877{
878 using namespace GRAINS_TREE;
879 if( eNgr<1 ) return 0;
880
883 eXview = eView->GetXview();
884 eYview = eView->GetYview();
885 EdbCluster *cl=0;
886 int ncl=0;
887 for(int i=0; i<eNgr; i++) {
888 eSeg = (EdbSegment*)(eG->UncheckedAt(i));
889 //printf("id=%d\n",s->GetID());
890 ncl = eSeg->GetElements()->GetEntriesFast();
891 if(ncl>1) {
892 for(int j=0; j<ncl; j++) {
893 cl = (EdbCluster*)eSeg->GetElements()->UncheckedAt(j);
894 new((*eClust)[j]) EdbCluster( *cl );
895 }
896 eGrainsTree->Fill();
897 }
898 eClust->Clear();
899 }
900 eGrainsTree->AutoSave();
901 return eNgr;
902}
TObjArray * GetElements() const
Definition: EdbSegment.h:117
Float_t GetXview() const
Definition: EdbView.h:193
Int_t GetViewID() const
Definition: EdbView.h:190
Int_t GetAreaID() const
Definition: EdbView.h:191
Float_t GetYview() const
Definition: EdbView.h:194
Definition: EdbViewRec.cxx:848
Float_t eYview
Definition: EdbViewRec.cxx:854
Int_t eViewID
Definition: EdbViewRec.cxx:851
Float_t eXview
Definition: EdbViewRec.cxx:853
TClonesArray * eClust
Definition: EdbViewRec.cxx:849
Int_t eAreaID
Definition: EdbViewRec.cxx:852
EdbSegment * eSeg
Definition: EdbViewRec.cxx:850

◆ FindGrains()

int EdbViewRec::FindGrains ( int  option = 0)
568{
569 // if option = 0 (default) grains reconstructed as segments and added to eGSeg array
570 // if option = 1 grains reconstructed as clusters and added to eGCla array
571
572 if(eVC->eNcl<1) return 0;
573 if(eVC->eNfr<2) return 0;
574 printf("FindGrains with %d clusters, %d frames:\n\n",eVC->eNcl,eVC->eNfr);
575 if(eNgr) {eG->Delete(); eNgr=0;}
576
577 EdbCluster **arr0=0, **arr1=0;
578 EdbCluster *cl0=0, *cl1=0;
579 float dx,dy;
580 int iz1=0, jn=0;
581 int n0=0,n1=0;
582 int jcell=0, jcell1=0;
583
584 // set to clusters eSegment attribute:
585
586 for(int iz=0; iz<eVC->eNfr-1; iz++) {
587 iz1=iz+1;
588 for(int j=0; j<eVC->eNcellXY; j++) {
589 jcell = eVC->Jcell(j,iz);
590 n0 = eVC->eNC[jcell];
591 if(n0<1) continue;
592 arr0 = eVC->GetCell(jcell);
593
594 for(int i0=0; i0<n0; i0++) {
595 cl0 = arr0[i0];
596
597 for(int in=0; in<9; in++) { //neighbours
598 jn = j + eVC->Jneib(in);
599 if(jn<0) continue;
600 if(jn>=eVC->eNcellXY) continue;
601
602 jcell1 = eVC->Jcell(jn,iz1);
603 n1 = eVC->eNC[jcell1];
604 if(n1<1) continue;
605 arr1 = eVC->GetCell(jcell1);
606
607 for(int i1=0; i1<n1; i1++) {
608 cl1 = arr1[i1];
609
610 dx = cl1->eX - (cl0->eX + eTXopt*(cl0->eX-eX0opt)*(cl1->eZ - cl0->eZ)); // with off-axis correction
611 if( Abs(dx) > eGrainNbin*eGrainSX ) continue;
612 dy = cl1->eY - (cl0->eY + eTYopt*(cl0->eY-eY0opt)*(cl1->eZ - cl0->eZ));
613 if( Abs(dy) > eGrainNbin*eGrainSY ) continue;
614
615 if( cl0->eSegment < 0 ) cl0->eSegment = eNgr++;
616 cl1->eSegment = cl0->eSegment;
617 }
618 }
619 }
620 }
621 }
622
623 eVC->CalcStat();
624
625 // assemble grains lists (as a segments)
626
627 for(int i=0; i<eVC->eNcl; i++) {
628 cl0 = (EdbCluster*)eCL->UncheckedAt(i);
629 if(!eG) printf("\n***** eG is 0!!!!!!\n\n");
630 if( cl0->eSegment<0 ) cl0->eSegment = eNgr++;
631
632 //printf("eNgr = %d eNgrMax = %d \n",eNgr, eNgrMax);
633 //cl0->Print();
634
635 if(eNgr>eNgrMax) {printf("ERROR: too many grains: %d - stop tracking!!\n",eNgr); eNgr--; return -1; }
636 if(!eG->UncheckedAt(cl0->eSegment)) new((*eG)[cl0->eSegment]) EdbSegment();
637 ((EdbSegment*)eG->UncheckedAt(cl0->eSegment))->AddElement(cl0);
638 }
639 printf("ncl: %d ngr: %d ( %3.1f %% )\n", eVC->eNcl,eNgr, 100.*eNgr/eVC->eNcl );
640
641 // fit grains to segments or to clusters:
642
643 int iseg;
644 EdbSegment *s=0;
645 EdbCluster c;
646 for(int i=0; i<eNgr; i++) {
647 s = (EdbSegment*)eG->UncheckedAt(i);
648 s->SetID(i);
649 if( s->GetNelements() < eNclGrMin) s->UnSetIDE(); // release clusters
650 else {
651 if(option==1) { // grains preprocessing before tracking
652 if( s->GetNelements() > eNclGrMax ) { // long grain is a segment
653 GoodSegment(*s,1);
654 iseg = eSA->GetLast()+1;
655 s->SetIDE(iseg);
656 new((*eSA)[iseg]) EdbSegment( *s );
657 } else { // short grain is a cluster
658 FitSegmentToCl(*s,c, 1);
659 new((*eGCla)[eGCla->GetLast()+1]) EdbCluster(c);
660 s->UnSetIDE();
661 }
662 }
663 else FitSegment(*s,1);
664 }
665 s->SetPuls(s->GetNelements());
666 }
667
668 if(option==1) {
669 printf("%d segments has been converted from grains\n",eSA->GetLast()+1);
671 printf("WARNING: FindGrains: eNclGrMax(%d) < ePulsMin(%d) : can loose segments!\n",eNclGrMax,ePulsMin);
672 }
673
674 return eNgr;
675}
Int_t eSegment
segment id to be attached (-1 if no segment)
Definition: EdbCluster.h:30
Int_t * eNC
Definition: EdbViewRec.h:40
Int_t eNfr
Definition: EdbViewRec.h:29
EdbCluster ** GetCell(int j) const
Definition: EdbViewRec.h:95
int Jcell(float x, float y, int ifr) const
Definition: EdbViewRec.h:89
int Jneib(int i) const
Definition: EdbViewRec.h:93
void CalcStat()
Definition: EdbViewRec.cxx:145
Int_t eNcl
Definition: EdbViewRec.h:35
Int_t eNcellXY
Definition: EdbViewRec.h:33
Float_t eTYopt
Definition: EdbViewDef.h:25
Float_t eX0opt
Definition: EdbViewDef.h:24
Float_t eGrainSX
Definition: EdbViewDef.h:30
Float_t eGrainSY
Definition: EdbViewDef.h:30
Float_t eY0opt
Definition: EdbViewDef.h:24
Float_t eTXopt
Definition: EdbViewDef.h:25
int eNclGrMax
Definition: EdbViewRec.h:136
Float_t eGrainNbin
Definition: EdbViewRec.h:109
bool GoodSegment(EdbSegment &s, int wkey=0)
Definition: EdbViewRec.cxx:950
Int_t ePulsMin
Definition: EdbViewRec.h:154
static int FitSegment(EdbSegment &s, int wkey=0)
Definition: EdbViewRec.cxx:728
EdbViewCell * eVC
Definition: EdbViewRec.h:117
static int FitSegmentToCl(EdbSegment &s, EdbCluster &c, int wkey=0)
Definition: EdbViewRec.cxx:776
Int_t eNgrMax
pointer to the input view currently in processing
Definition: EdbViewRec.h:114
int eNclGrMin
debug tree
Definition: EdbViewRec.h:135
TClonesArray * eCL
[eNgrMax] array of grains represented as clusters
Definition: EdbViewRec.h:125

◆ FindSeeds()

int EdbViewRec::FindSeeds ( )
1195{
1196
1197 if(eVC->eNcl<1) return 0;
1198 if(eVC->eNfr<2) return 0;
1199 printf("FindSeed with %d clusters, %d frames \t eZcenter[%6.1f:%6.1f] = %6.1f eZxy = %6.1f\n",
1201
1202#ifdef WIN32
1203 memset(ehX,'\0',enXtot*sizeof(Short_t));
1204#else
1205 bzero(ehX,enXtot*sizeof(Short_t));
1206#endif
1207
1208 EdbCluster **arr0=0, **arr1=0;
1209 EdbCluster *cl0=0, *cl1=0;
1210 int jcell0=0,jcell1=0;
1211 int n0=0,n1=0;
1212 int iz1=0, jn=0;
1213 long ncycle=0, ict=0, nseeds=0;
1214 float r2max=eVC->eBinX*eVC->eBinX;
1215 float r2d3max = eRmax*eRmax;
1216 float dx,dy,dz,r2, r2d3;
1217 float r,phi,t;
1218 float zratio,xZxy,yZxy;
1219 int it,ip,ix,iy;
1220 float x0,x1,y0,y1, phirot;
1221 EdbCluster **pps;
1222
1223 for(int istep=eStepFrom; istep<eStepTo; istep+=eStep) {
1224 for(int iz=0; iz<eVC->eNfr-istep; iz++) {
1225 iz1=iz+istep;
1226
1227 for(int j=0; j<eVC->eNcellXY; j++) {
1228 jcell0 = eVC->Jcell(j,iz);
1229 n0 = eVC->eNC[jcell0];
1230 if(n0<1) continue;
1231 arr0 = eVC->GetCell(jcell0);
1232
1233 for(int i0=0; i0<n0; i0++) {
1234 cl0 = arr0[i0];
1235
1236 for(int in=0; in<9; in++) { //neighbours
1237 jn = j + eVC->Jneib(in);
1238 if(jn<0) continue;
1239 if(jn>=eVC->eNcellXY) continue;
1240
1241 jcell1 = eVC->Jcell(jn,iz1);
1242 n1 = eVC->eNC[jcell1];
1243 if(n1<1) continue;
1244 arr1 = eVC->GetCell(jcell1);
1245
1246 for(int i1=0; i1<n1; i1++) {
1247 cl1 = arr1[i1];
1248 ncycle++;
1249
1250
1251 x0 = cl0->eX;
1252 x1 = cl1->eX;
1253 y0 = cl0->eY;
1254 y1 = cl1->eY;
1255 dx = x1 - x0;
1256 dy = y1 - y0;
1257 dz = cl1->eZ - cl0->eZ;
1258 r2 = dx*dx+dy*dy;
1259 r2d3 = r2+dz*dz;
1260 if(r2>r2max) continue;
1261 if(r2d3>r2d3max) continue;
1262 r = Sqrt(r2);
1263 t = ATan2(r,Abs(dz));
1264 if(t>eTheta[enT-1]) continue;
1265 phi = ATan2(dy,dx); // defined in [ -pi : pi ]
1266
1267 for(it=0; it<enT-1; it++) if(t<eTheta[it]) break;
1268 ip = (int)((Pi()+phi-0.00001)/esP[it]+0.00001);
1269
1270 if(it>0) { //rotation:
1271 phirot = -((ip+0.5)*esP[it] - Pi()); // rotation angle is the same for all entries of the same phi-bin
1272 x0 = cl0->eX*Cos(phirot) - cl0->eY*Sin(phirot);
1273 y0 = cl0->eX*Sin(phirot) + cl0->eY*Cos(phirot);
1274 x1 = cl1->eX*Cos(phirot) - cl1->eY*Sin(phirot);
1275 y1 = cl1->eX*Sin(phirot) + cl1->eY*Cos(phirot);
1276 dx = x1-x0;
1277 dy = y1-y0;
1278 }
1279
1280 //printf("dx,dy,dz: %f %f %f\n",dx,dy,dz);
1281
1282 zratio = (eZcenter - cl0->eZ)/dz;
1283 xZxy = eDmax/2. + x0 + dx*zratio;
1284 if(xZxy<0) continue;
1285 if(xZxy>eDmax) continue;
1286
1287 yZxy = eDmax/2. + y0 + dy*zratio;
1288 if(yZxy<0) continue;
1289 if(yZxy>eDmax) continue;
1290
1291 iy = (int)(yZxy/esY[it]);
1292 ix = (int)(xZxy/esX[it]);
1293
1294 //printf("ict[ %d %d %d %d ]\n",it,ip,iy,ix);
1295 ict = epT[it][ip][iy] + ix -ehX;
1296 //printf("%d\n",ict);
1297
1298 if(ict<0||ict>=enXtot) {
1299 printf("ict[ %d %d %d ] = %ld\n",it,ip,iy,ict);
1300 printf("esX, esY: %f %f\n",esX[it], esY[it]);
1301 printf("eYmin, eYmax: %f %f\n",eYmin,eYmax );
1302 printf("eXmin, eXmax: %f %f\n",eXmin,eXmax );
1303 printf("yZxy, xZxy: %f %f\n",yZxy,xZxy);
1304 return 0;
1305 }
1306
1307 if(ehX[ict]<1) {
1308 epS[ict] = epC + nseeds*eSeedLim;
1309 nseeds++;
1310 if( nseeds >= enSeedsLim) {
1311 printf("ERROR: too many seeds: %ld, abandon tracking! \n",nseeds);
1312 return 0;
1313 }
1314 }
1315 pps = epS[ict] + ehX[ict];
1316
1317 if(ehX[ict]<eSeedLim) *pps = cl0;
1318 pps++;
1319 //printf("ict = %d bin = %d\n",ict, ehX[ict]);
1320 ehX[ict]++;
1321 if(ehX[ict]<eSeedLim) *pps=cl1;
1322 ehX[ict]++;
1323
1324 }
1325 }
1326 }
1327 }
1328 }
1329 }
1330 printf("occupancy: %ld / %d * %4.3f = %5.2f %% \n", nseeds,enXtot, eFact, 100.*eFact*nseeds/enXtot);
1331 return nseeds;
1332}
brick dz
Definition: RecDispMC.C:107
Float_t eBinX
Definition: EdbViewRec.h:27
Float_t eYmin
Definition: EdbViewDef.h:21
Float_t eXmax
Definition: EdbViewDef.h:20
Float_t eZmax
Definition: EdbViewDef.h:17
Float_t eZxy
Definition: EdbViewDef.h:18
Float_t eZmin
Definition: EdbViewDef.h:17
Float_t eYmax
Definition: EdbViewDef.h:21
Float_t eXmin
Definition: EdbViewDef.h:20
Short_t eSeedLim
Definition: EdbViewRec.h:151
Float_t eRmax
Definition: EdbViewRec.h:184
TArrayF esP
Definition: EdbViewRec.h:171
Int_t eStep
Definition: EdbViewRec.h:143
Float_t eZcenter
Definition: EdbViewRec.h:159
TArrayF eTheta
Definition: EdbViewRec.h:169
Float_t eDmax
pointers to clusters
Definition: EdbViewRec.h:182
Int_t enSeedsLim
Definition: EdbViewRec.h:150
TArrayF esY
Definition: EdbViewRec.h:172
TArrayF esX
Definition: EdbViewRec.h:173
Int_t eStepFrom
Definition: EdbViewRec.h:144
Int_t eStepTo
Definition: EdbViewRec.h:145
TTree * t
Definition: check_shower.C:4
void r(int rid=2)
Definition: test.C:201

◆ FitSegment()

int EdbViewRec::FitSegment ( EdbSegment s,
int  wkey = 0 
)
static
729{
730 // if wkey = 0 - equal weights (default)
731 // = 1 - weighted with cluster area
732
733 EdbCluster *c1=0;
734 TObjArray *arr = s.GetElements();
735 if(!arr) return -1;
736 int ncl = arr->GetEntriesFast();
737 if(!ncl) return 0;
738
739 if(ncl==1) {
740 c1 = (EdbCluster*)arr->UncheckedAt(0);
741 s.Set( c1->X(), c1->Y(), c1->Z(), 0.,0., 0., s.GetSide(), ncl, s.GetID() );
742 return ncl;
743 }
744 else {
745 float *X = new float[ncl];
746 float *Y = new float[ncl];
747 float *Z = new float[ncl];
748 float *W = new float[ncl];
749 for(int i=0; i<ncl; i++) {
750 c1 = (EdbCluster*)arr->UncheckedAt(i);
751 if(wkey==1) W[i] = c1->eArea;
752 else W[i] = 1.;
753 X[i] = c1->eX;
754 Y[i] = c1->eY;
755 Z[i] = c1->eZ;
756 }
757 float zmin=Z[0],zmax=Z[0];
758 for(int i=1; i<ncl; i++) {
759 if(zmin>Z[i]) zmin=Z[i];
760 if(zmax<Z[i]) zmax=Z[i];
761 }
762 float x0,y0,z0,tx,ty,ex,ey;
763 EdbMath::LFIT3( X, Y, Z, W, ncl,
764 x0, y0, z0, tx, ty, ex, ey );
765 s.Set( x0,y0,z0, tx,ty, zmax-zmin, s.GetSide(), ncl, s.GetID());
766 s.SetSigma(ex,ey);
767 delete[] X;
768 delete[] Y;
769 delete[] Z;
770 delete[] W;
771 }
772 return ncl;
773}
brick z0
Definition: RecDispMC.C:106
static int LFIT3(float *X, float *Y, float *Z, float *W, int L, float &X0, float &Y0, float &Z0, float &TX, float &TY, float &EX, float &EY)
Definition: EdbMath.cxx:313
TCanvas * c1
Definition: energy.C:13
Double_t X
Definition: tlg2couples.C:76
Double_t Y
Definition: tlg2couples.C:76
Double_t Z
Definition: tlg2couples.C:104
Int_t W
Definition: testBGReduction_By_ANN.C:15

◆ FitSegmentToCl()

int EdbViewRec::FitSegmentToCl ( EdbSegment s,
EdbCluster c,
int  wkey = 0 
)
static
777{
778 // Special function to calculate the grain parameters and fill the output cluster c
779 // if wkey = 0 - equal weights (default)
780 // = 1 - weighted with cluster area
781
782 EdbCluster *c1=0;
783 TObjArray *arr = s.GetElements();
784 if(!arr) return -1;
785 int ncl = arr->GetEntriesFast();
786 if(!ncl) return 0;
787
788 if(ncl==1) {
789 c1 = (EdbCluster*)arr->UncheckedAt(0);
790 c.Set( c1->eX, c1->eY, c1->eZ, c1->eArea, c1->eVolume, c1->eFrame, c1->eSide, c1->eSegment );
791 return ncl;
792 }
793 else {
794 double SX=0,SY=0,SZ=0;
795 float SA=0,SV=0,SF=0,SW=0;
796 float w;
797 for(int i=0; i<ncl; i++) {
798 c1 = (EdbCluster*)arr->UncheckedAt(i);
799 if(wkey==1) w = c1->eArea;
800 else w = 1.;
801 SX += w*c1->eX;
802 SY += w*c1->eY;
803 SZ += w*c1->eZ;
804 SF += w*c1->eFrame;
805 SW += w;
806 SA += c1->eArea;
807 SV += c1->eVolume;
808 }
809 c.Set( (float)(SX/SW), (float)(SY/SW), (float)(SZ/SW), SA, SV, Nint(SF/SW), c1->eSide, c1->eSegment );
810 }
811 return ncl;
812}
void Set(float x, float y, float z, float a, float v, int f, int s, int seg=-1)
Definition: EdbCluster.h:43
Double_t SY
Definition: tlg2couples.C:106
Double_t SX
Definition: tlg2couples.C:105
Double_t SZ
Definition: tlg2couples.C:107
void w(int rid=2, int nviews=2)
Definition: test.C:27

◆ GoodSegment()

bool EdbViewRec::GoodSegment ( EdbSegment s,
int  wkey = 0 
)
951{
952 // TODO:
953 // - probability to be a straight line (chi2)
954 // - probability to be a grain (big&long)
955 // - probability to be a casual coincidence of more grains
956
957 if(!s.GetElements()) return false;
958 s.SetPuls(s.GetNelements());
959 if( s.GetPuls() < ePulsMin ) return false;
960 s.SetSide(((EdbCluster*)(s.GetElements()->At(0)))->eSide);
961 FitSegment( s, wkey );
962 if( Sqrt(s.GetTx()*s.GetTx()+s.GetTy()*s.GetTy()) > eThetaLim ) return false;
963 s.SetSigma( Sqrt(s.GetSigmaX()*s.GetSigmaX() + s.GetSigmaY()*s.GetSigmaY()),
965
966 if( s.GetSigmaY() > eSigmaMax ) return false;
967
968 return true;
969}
Float_t eGrainSZ
Definition: EdbViewDef.h:30
Float_t eSigmaMax
Definition: EdbViewRec.h:157
Float_t eThetaLim
Definition: EdbViewRec.h:141
float CalculateSegmentChi2(EdbSegment &seg, float sx, float sy, float sz)
Definition: EdbViewRec.cxx:815

◆ Init()

bool EdbViewRec::Init ( void  )
268{
269 // allocate memory here, return true if ok
270 // all settings must be done before!
271
272 printf("\nInit reconstructor:\n\n");
273
274 if( eSeedThres.GetSize() != enT ) { // set default threshold for the seed value
275 eSeedThres.Set(enT);
276 for(int i=0; i<enT; i++) eSeedThres[i]=eSeedThres0;
277 }
278 if( eNseedMax.GetSize() != enT ) { // set default threshold for the seed value
279 eNseedMax.Set(enT);
280 for(int i=0; i<enT; i++) eNseedMax[i]=eNseedMax0;
281 }
282
283 enP.Set(enT); // number of Phi divisions [enT]
284 eTheta.Set(enT); // limits of theta divisions [enT]
285 eR.Set(enT); // limits of theta divisions [enT]
286 esP.Set(enT); // number of Phi divisions [enT]
287 esY.Set(enT); // step of Y divisions [enT]
288 esX.Set(enT); // step of X divisions [enT]
289
290 printf("enT \t= %d\n",enT);
291 float st=0, sp=0, t0=0, dt = 0.;
292 for(int it=0; it<enT; it++) {
293 st = 1.8*SThetaGr( t0+dt , 0., eDZmin , eGrainSX, eGrainSY, eGrainSZ);
294 dt = st/2.;
295 eTheta[it] = t0 + st;
296 t0 = eTheta[it];
297 eR[it] = eDZmin*Tan(eTheta[it]);
298 sp = Min( TwoPi(), 2.*SPhiGr( eTheta[it]-st , 0., eDZmin , eGrainSX, eGrainSY, eGrainSZ) );
299 enP[it] = (int)(TwoPi()/sp + 0.001);
300 esP[it] = TwoPi()/enP[it]; // reset Sphi for having the exactly equal sectors
301
302 esX[it] = st/Cos(eTheta[it])*20. * 2.;
303 esY[it] = eTheta[0]/Cos(eTheta[0])*20. * 2.;
304 }
305
306 enPtot = (int)(enP.GetSum()); // total number of phi divisions
307 enY.Set(enPtot); // number of Y divisions [enPtot]
308 printf("enPtot \t= %d\n",enPtot);
309
310 eDmax = Sqrt( (eYmax-eYmin)*(eYmax-eYmin) + (eXmax-eXmin)*(eXmax-eXmin) );
311 eFact = (eDmax*eDmax) / ((eYmax-eYmin)*(eXmax-eXmin));
312
313 int kp=0;
314 for(int it=0; it<enT; it++) {
315 for(int ip=0; ip<enP[it]; ip++) {
316 enY[kp] = (int)(eDmax/esY[it] + 1);
317 kp++;
318 }
319 }
320 enYtot = (int)(enY.GetSum()); // total number of Y divisions
321 enX.Set(enYtot); // number of X divisions [enYtot]
322 printf("enYtot \t= %d\n",enYtot);
323
324 kp=0;
325 int ky=0;
326 for(int it=0; it<enT; it++) {
327 for(int ip=0; ip<enP[it]; ip++) {
328 for(int iy=0; iy<enY[kp]; iy++) {
329 enX[ky] = (int)(eDmax/esX[it] + 1);
330 ky++;
331 }
332 kp++;
333 }
334 }
335 enXtot = (int)(enX.GetSum()); // total number of X divisions
336 printf("enXtot \t= %d\n",enXtot);
337 printf("eDmax \t= %f eFact = %f \n",eDmax, eFact);
338
339 ehX = new Short_t[enXtot]; //- phase histogram
340 epY = new Short_t*[enYtot]; //- pointers to the first x[iy]
341 epP = new Short_t**[enPtot]; //- pointers to the first y[ip]
342 epT = new Short_t***[enT]; //- pointers to the first phi[it]
343
344 epS = new EdbCluster**[enXtot]; //- pointers to seeds lists (in epC)
345 epC = new EdbCluster*[enSeedsLim*eSeedLim]; //- pointers to clusters
346
347 Short_t *px = ehX;
348 epY[0] = px;
349 for(int iy=1; iy<enYtot; iy++) { px += enX[iy-1]; epY[iy]=px; }
350
351 Short_t **py = epY;
352 epP[0] = py;
353 for(int ip=1; ip<enPtot; ip++) { py += enY[ip-1]; epP[ip]=py; }
354
355 Short_t ***pp = epP;
356 epT[0] = &(epP[0]);
357 for(int it=1; it<enT; it++) { pp += enP[it-1]; epT[it]=pp; }
358
359 for(int it=0; it<enT; it++)
360 printf("%d theta = %5.4f (r= %4.2f)\t esP = %f x %d \t sx = %f sy = %f \tseedthres = %d \tNseedsmax = %d\n",
361 it, eTheta[it], eR[it],
362 esP[it], enP[it],
363 esX[it], esY[it],
364 eSeedThres[it],
365 eNseedMax[it]);
366
367 eVCC.Init();
368 if(eDoGrainsProcessing) { // segments reconstruction with grains preprocessing
369 eVCG.Init();
370 eGCla = new TClonesArray("EdbCluster",eNgrMax);
371 }
372
373 // if(!eAddGrainsToView)
374 if(eNgrMax>0) eG = new TClonesArray("EdbSegment",eNgrMax);
375
376 if(eNsegMax>0) if(!eSA) eSA = new TClonesArray("EdbSegment",eNsegMax);
377
378 return true;
379}
void Init()
Definition: EdbViewRec.cxx:105
float SThetaGr(float theta, float phi, float dz, float sx, float sy, float sz)
Definition: EdbViewRec.cxx:394
TArrayI enP
Definition: EdbViewRec.h:162
EdbViewCell eVCC
pointer to eVCC or eVCG
Definition: EdbViewRec.h:118
Int_t enPtot
Definition: EdbViewRec.h:163
EdbViewCell eVCG
cells with raw clusters
Definition: EdbViewRec.h:119
Int_t enYtot
Definition: EdbViewRec.h:165
Int_t eNseedMax0
Definition: EdbViewRec.h:152
TArrayI enX
Definition: EdbViewRec.h:166
Float_t eDZmin
Definition: EdbViewRec.h:140
bool eDoGrainsProcessing
Definition: EdbViewRec.h:106
TArrayF eR
Definition: EdbViewRec.h:149
float SPhiGr(float theta, float phi, float dz, float sx, float sy, float sz)
Definition: EdbViewRec.cxx:403
TArrayI enY
Definition: EdbViewRec.h:164

◆ InitGrainsTree()

void EdbViewRec::InitGrainsTree ( const char *  file = "grains.root")

859{
860 using namespace GRAINS_TREE;
861
862 TFile *f;
863 f = new TFile(file,"RECREATE");
864 eGrainsTree = new TTree("grains","grains");
865 eClust = new TClonesArray("EdbCluster");
866 eSeg = 0;
867 eGrainsTree->Branch("view",&eViewID,"view/I");
868 eGrainsTree->Branch("area",&eAreaID,"area/I");
869 eGrainsTree->Branch("xv",&eXview,"xv");
870 eGrainsTree->Branch("yv",&eYview,"xv");
871 eGrainsTree->Branch("c",&eClust);
872 eGrainsTree->BranchOld("s","EdbSegment",&eSeg);
873}
FILE * f
Definition: RecDispMC.C:150
TFile * file
Definition: write_pvr.C:3

◆ InitR()

void EdbViewRec::InitR ( )
383{
384 eR[0] = .5; //microns
385 eR[1] = 1.;
386 eR[2] = 2.;
387 eR[3] = 4.;
388 eR[4] = 6.;
389 eR[5] = 8.;
390 eR[6] = 10.;
391}

◆ MergeSegments()

int EdbViewRec::MergeSegments ( )
1069{
1070 int nseg = eSA->GetLast()+1;
1071 if(nseg<2) return nseg;
1072
1073 int nmerge=0;
1074 TObject *c;
1075 EdbSegment *s;
1076 EdbSegment *s1,*s2;
1077 for(int i=0; i<nseg-1; i++) {
1078 s1 = (EdbSegment*)(eSA->UncheckedAt(i));;
1079 for(int j=i+1; j<nseg; j++) {
1080 s2 = (EdbSegment*)(eSA->UncheckedAt(j));
1081 if( s1->GetPuls() < 1 ) continue;
1082 if( s2->GetPuls() < 1 ) continue;
1083 if( Abs(s1->GetX0() - s2->GetX0()) > 40.) continue;
1084 if( Abs(s1->GetY0() - s2->GetY0()) > 40.) continue;
1085 if( Abs(s1->GetTx() - s2->GetTx()) > 0.12) continue;
1086 if( Abs(s1->GetTy() - s2->GetTy()) > 0.12) continue;
1087 if( Chi2Seg( *s1,*s2 ) > 3. ) continue;
1088
1089 s = new EdbSegment();
1090 s->Copy(*s1);
1091 int nsame=0;
1092 for(int ii=0; ii<s2->GetNelements(); ii++) {
1093 c = s2->GetElements()->UncheckedAt(ii);
1094 if( !(s->GetElements()->FindObject(c)) ) s->AddElement( c );
1095 else nsame++;
1096 }
1097 if(nsame > s->GetNelements()/2) { // same clusters - discard the second segment
1098 s2->UnSetIDE();
1099 s2->SetPuls(0);
1100 nmerge++;
1101 } else if( GoodSegment(*s) ) { // s1 and s2 are parts of the same segment
1102 s->SetID(s1->GetID());
1103 s1->Copy(*s);
1104 s2->UnSetIDE();
1105 s2->SetPuls(0);
1106 nmerge++;
1107 }
1108 delete s;
1109 }
1110 }
1111
1112 int good=0;
1113 nseg = eSA->GetLast()+1;
1114 for(int i=0; i<nseg; i++) {
1115 s = (EdbSegment*)(eSA->UncheckedAt(i));
1116 if(s->GetPuls()<1) continue;
1117 s->SetIDE(i);
1118 good++;
1119 }
1120
1121 printf("%d segments after %d merging\n", good, nmerge);
1122 return good;
1123}
void Copy(const EdbSegP &s)
Definition: EdbSegP.cxx:105
float Chi2Seg(EdbSegment &s1, EdbSegment &s2)
Definition: EdbViewRec.cxx:1126
EdbSegP * s1
Definition: tlg2couples.C:29
EdbSegP * s2
Definition: tlg2couples.C:30

◆ ReconstructGrains()

int EdbViewRec::ReconstructGrains ( )
456{
457 // stearing function for grains reconstruction
458 // options: - save grains as the segments into eView
459 // - keep them as an independent array and then save to separate grains tree
460
461 int ngr=0;
462
463 ngr = FindGrains();
464
465 if(eSA!=eG) { delete eSA; eSA=eG; } // no tracking in this routine!
466 //FillGrainsTree();
467
468 return ngr;
469}
int FindGrains(int option=0)
Definition: EdbViewRec.cxx:567

◆ ReconstructSegments()

int EdbViewRec::ReconstructSegments ( )
473{
474 // stearing function for segments reconstruction
475 // options: - with grains preprocessing
476 // - without grains preprocessing
477
478 int nseg=0;
479 int ngr=0;
480
481 eSA->Clear(); // clear segments array
482
484 ngr = FindGrains(1); // 1- fit grains to clusters
485 eVC = &eVCG;
486 eCL = eGCla; // switch processing to grains
487 eVC->FillCell( *(eCL) );
488 }
489
490 FindSeeds();
491 printf("FindSeeds ok\n");
495
496 if(eDoGrainsProcessing) { // switch processing back to clusters
497 eCL->Clear();
498 eVC = &eVCC;
499 eCL = eView->GetClusters();
500 }
501
502 RefitSegments(1);
503
504 return nseg;
505}
int FillCell(TClonesArray &v)
Definition: EdbViewRec.cxx:121
int MergeSegments()
Definition: EdbViewRec.cxx:1068
int SelectSegments()
Definition: EdbViewRec.cxx:905
int FindSeeds()
Definition: EdbViewRec.cxx:1194
bool eCheckSeedThres
Definition: EdbViewRec.h:107
int RefitSegments(int wkey=0)
Definition: EdbViewRec.cxx:1038
int CheckSeedThres()
Definition: EdbViewRec.cxx:1132
TClonesArray * GetClusters() const
Definition: EdbView.h:164

◆ RefillSegment()

int EdbViewRec::RefillSegment ( EdbSegment s)
973{
974 // attach to the segment all clusters(grains) close to him
975 // assumed the segment is fitted
976
977 float SX = Sqrt(eGrainSX*eGrainSX + eClaSX*eClaSX);
978 float SY = Sqrt(eGrainSY*eGrainSY + eClaSY*eClaSY);
979 float SZ = Sqrt(eGrainSZ*eGrainSZ + eClaSZ*eClaSZ);
980
981 float dz0 = s.GetZ0() - eZmin;
982 float dz1 = eZmax - s.GetZ0();
983 float v[3];
984 // segment line in a normalized space:
985 float start[3]= { (s.GetX0()-dz0*s.GetTx())/SX,
986 (s.GetY0()-dz0*s.GetTy())/SY,
987 eZmin/SZ };
988 float end[3] = { (s.GetX0()+dz1*s.GetTx())/SX,
989 (s.GetY0()+dz1*s.GetTy())/SY,
990 eZmax/SZ };
991
992 float sx = Sqrt( eGrainSZ*eGrainSZ*s.GetTx()*s.GetTx() + eGrainSX*eGrainSX );
993 float sy = Sqrt( eGrainSZ*eGrainSZ*s.GetTy()*s.GetTy() + eGrainSY*eGrainSY );
994
995 EdbCluster **cla=0;
996 TObjArray *ela=s.GetElements();
997 ela->Clear();
998 EdbCluster *c=0;
999 float x0,y0,x,y,dz;
1000 int n1, lx, ly, lx0,ly0;
1001 int jcell=0;
1002 for(int i=0; i<eVC->eNfr; i++) {
1003 dz = eView->GetFrame(i)->GetZ() - s.GetZ0();
1004 x0 = s.GetX0() + s.GetTx()*dz;
1005 y0 = s.GetY0() + s.GetTy()*dz;
1006 lx0 = eVC->IXcell(x0);
1007 ly0 = eVC->IYcell(y0);
1008 for(int ix=-1; ix<2; ix++)
1009 for(int iy=-1; iy<2; iy++) {
1010 x = x0 + ix*sx;
1011 if(x<eVC->eXmin||x>eVC->eXmax) continue;
1012 y = y0 + iy*sy;
1013 if(y<eVC->eYmin||y>eVC->eYmax) continue;
1014 lx = eVC->IXcell(x);
1015 ly = eVC->IYcell(y);
1016 if(ix!=0&&iy!=0)
1017 if( lx == lx0 && ly == ly0 ) continue;
1018 jcell = eVC->Jcell(lx,ly,i);
1019 n1 = eVC->eNC[jcell];
1020 if(!n1) continue;
1021 cla = eVC->GetCell(jcell);
1022 for(int i1=0; i1<n1; i1++) {
1023 c = cla[i1];
1024 v[0] = c->eX/SX;
1025 v[1] = c->eY/SY;
1026 v[2] = c->eZ/SZ;
1027 if( EdbMath::DistancePointLine3(v,start,end, true) < 2. )
1028 if( !(ela->FindObject(c)) ) ela->Add( c );
1029 }
1030 }
1031 }
1032 int nn = ela->GetEntriesFast();
1033 s.SetPuls(nn);
1034 return nn;
1035}
float GetZ() const
Definition: EdbFrame.h:42
Float_t eYmax
Definition: EdbViewRec.h:24
Float_t eXmax
Definition: EdbViewRec.h:23
int IYcell(float y) const
Definition: EdbViewRec.h:87
int IXcell(float x) const
Definition: EdbViewRec.h:86
Float_t eClaSX
Definition: EdbViewDef.h:28
Float_t eClaSZ
Definition: EdbViewDef.h:28
Float_t eClaSY
Definition: EdbViewDef.h:28
EdbFrame * GetFrame(int i) const
Definition: EdbView.h:221

◆ RefitSegments()

int EdbViewRec::RefitSegments ( int  wkey = 0)
1039{
1040 int good=0;
1041 int nseg = eSA->GetLast()+1;
1042 EdbSegment *s;
1043 for(int i=0; i<nseg; i++) {
1044 s = (EdbSegment*)(eSA->UncheckedAt(i));
1045 if(s->GetPuls()<1) continue;
1046 s->UnSetIDE();
1047 RefillSegment(*s);
1048 if( !GoodSegment(*s,wkey) ) {
1049 s->UnSetIDE();
1050 s->SetPuls(0);
1051 continue;
1052 }
1053 if( s->GetSigmaY() < eSigmaMin || s->GetSigmaY() > eSigmaMax ) {
1054 s->UnSetIDE();
1055 s->SetPuls(0);
1056 continue;
1057 }
1058
1059 s->SetIDE(i);
1060 good++;
1061 }
1062
1063 printf("%d segments after refitting\n", good);
1064 return good;
1065}
int RefillSegment(EdbSegment &s)
Definition: EdbViewRec.cxx:972
Float_t eSigmaMin
Definition: EdbViewRec.h:156

◆ ResetClustersSeg()

void EdbViewRec::ResetClustersSeg ( )
558{
559 // set sl.eSegment to the default value (-1) - important only for simulation
560
561 int ncl=eCL->GetLast()+1;
562 if(ncl>0)
563 for(int i=0; i<ncl; i++) ((EdbCluster*)eCL->UncheckedAt(i))->eSegment=-1;
564}

◆ SaveToOutputView()

bool EdbViewRec::SaveToOutputView ( EdbView vout,
int  do_h = 1,
int  do_c = 2,
int  do_s = 2,
int  do_tr = 0,
int  do_f = 2 
)
509{
510 // Put the result into the output view; return true if successful
511 //
512 // do_h : 0/1 (1) - header info: 1-save; 0 - do not save;
513 // do_c : 0/1/2 (2) - cluster info: 1-save all; 2-save attached to segments; 0- do not save;
514 // do_s : 0/1/2 (2) - segments info: 1-save all; 2-save selected; 0- do not save;
515 // do_tr: (0) - track info:
516 // do_f : 0/1/2 (2) - frames info: 1-save all; 2-do not save images(if any); 0-do not save;
517
518 if(do_h>0) vout.GetHeader()->Copy(eView->GetHeader());
519
520 if(do_s>0) {
521 int nseg = eSA->GetLast()+1;
522 EdbSegment *s;
523 for(int i=0; i<nseg; i++) {
524 s = (EdbSegment*)(eSA->UncheckedAt(i));
525 if(do_s==2) if(s->GetPuls()<1) continue;
526 if(ePropagateToBase) {
527 s->SetX0( s->GetX0() + s->GetTx()*(eZxy-s->GetZ0()) );
528 s->SetY0( s->GetY0() + s->GetTy()*(eZxy-s->GetZ0()) );
529 s->SetZ0(eZxy);
530 }
531 vout.AddSegment(s);
532 }
533 printf( "%d segments saved (ot of %d)\n", vout.Nsegments(),nseg );
534 }
535
536 if(do_c>0) {
537 int ncl = eCL->GetLast()+1;
538 EdbCluster *c=0;
539 for(int i=0; i<ncl; i++) {
540 c = (EdbCluster*)(eCL->UncheckedAt(i));
541 if(do_c==2) if(c->eSegment<0) continue;
542 vout.AddCluster(c);
543 }
544 }
545
546 if(do_f>0) {
547 int nfr = eView->GetNframes();
548 for(int i=0; i<nfr; i++) {
549 vout.AddFrame(eView->GetFrame(i));
550 }
551 }
552
553 return true;
554}
void Copy(EdbViewHeader *h)
Definition: EdbView.h:65
bool ePropagateToBase
Definition: EdbViewRec.h:108
EdbViewHeader * GetHeader() const
Definition: EdbView.h:163
EdbSegment * AddSegment(float x, float y, float z, float tx, float ty, float dz=0, int side=0, int puls=0, int id=-1)
Definition: EdbView.h:231
Int_t Nsegments() const
Definition: EdbView.h:216
EdbCluster * AddCluster(EdbCluster *c)
Definition: EdbView.h:225
void AddFrame(int id, float z, int ncl=0, int npix=0)
Definition: EdbView.cxx:268

◆ SelectSegments()

int EdbViewRec::SelectSegments ( )
906{
907 int good = 0;
908 int seeds = 0;
909
910 EdbSegment *s=0;
911 EdbCluster *c=0;
912 EdbCluster **pps;
913 int jmax, iseg;
914 Short_t *p0,*pnext, *plast = ehX+enXtot;
915 for(int it=0; it<enT; it++) {
916 p0 = epT[it][0][0];
917 if(it==enT-1) pnext = plast;
918 else pnext = epT[it+1][0][0];
919 while(p0<pnext) {
920 jmax = *p0;
921 if( jmax>=eSeedThres[it] ) {
922 int ict = p0-ehX;
923 pps = epS[ict];
924 jmax = (jmax<eSeedLim)?jmax:eSeedLim;
925 seeds++;
926 s = new EdbSegment();
927 for(int j=0; j<jmax; j++) {
928 c = *pps;
929 if( !(s->GetElements()) ) s->AddElement(c);
930 else if( !(s->GetElements()->FindObject(c)) ) s->AddElement(c);
931 pps++;
932 }
933 if( GoodSegment(*s) ) {
934 good++;
935 iseg = eSA->GetLast()+1;
936 s->SetIDE(iseg);
937 new((*eSA)[iseg]) EdbSegment( *s );
938 }
939 else delete s;
940 }
941 p0++;
942 }
943 }
944
945 printf("%d ( %d ) segments (seeds) are selected\n", good,seeds);
946 return good;
947}

◆ SetAddGrainsToView()

void EdbViewRec::SetAddGrainsToView ( bool  yesno)
inline
198{eAddGrainsToView=yesno;}
bool eAddGrainsToView
Definition: EdbViewRec.h:105

◆ SetClThres()

void EdbViewRec::SetClThres ( int  mina,
int  maxa 
)
inline
201{ eClMinA=mina; eClMaxA=maxa; }
Int_t eClMaxA
Definition: EdbViewRec.h:130
Int_t eClMinA
pointer to array of segments (output of tracking)
Definition: EdbViewRec.h:129

◆ SetNclGrLim()

void EdbViewRec::SetNclGrLim ( int  mincl,
int  maxcl 
)
inline
194{ eNclGrMin=mincl; eNclGrMax=maxcl; }

◆ SetNgrMax0()

void EdbViewRec::SetNgrMax0 ( Int_t  ngr)
inline
197{eNgrMax=ngr;}

◆ SetNSeedMax()

void EdbViewRec::SetNSeedMax ( int  nt,
int  th[] 
)
inline
207{ if(nt>0){ enT=nt; eNseedMax.Set(nt,th);} };

◆ SetNSeedMax0()

void EdbViewRec::SetNSeedMax0 ( int  n)
inline
206{ eNseedMax0=n; }

◆ SetPrimary()

void EdbViewRec::SetPrimary ( )
199{
200 // define default primary values used later for the calculations in Init() function
201 // before Init them can be modified by correspondent setters
202
203 eClMinA = 1; // by default accept all clusters
204 eClMaxA = 10000; // by default accept almoust all clusters
205
207 eVCC.SetBin(10,10); // cell size
208 eVCC.SetCellLimits(100000,20);
209
211 eVCG.SetBin(10,10); // cell size (TODO!)
212 eVCG.SetCellLimits(100000,20);
213
214 eDZmin = 5.87; // TODO!!!
215 enT = 7; // number of Theta divisions!
216
217 eThetaLim = 1.5; // max theta angle
218 enSeedsLim = 100000; // limit for the total number of seeds
219 eSeedLim = 48; // limit for the seed value
220 eNseedMax0 = 10; // starting limit for the good seeds (segments) to be analysed (used in CheckSeedThres)/theta
221 eStep = 1;
222 eStepFrom = 2;
223 eStepTo = 3;
224 eRmax = 11.; // limit on 3-dim distance between points [microns]
225 eSeedThres0 = 5; // limit for the seed value
226 eFact=1.; // occupancy correction factor
227
228 ePulsMin = 4; //default threshold for the segment puls value
229 ePulsMax = 500; //default threshold for the segment puls value
230
231 eNsegMax=10000; // limit for the segments number
232
233 eNgr=0;
234 eNgrMax=200000; // limit for the grains number
235 // eAddGrainsToView = false; // use a dedicated array for eG
236 eNclGrMin = 1;
237 eNclGrMax = 6;
238
239 eGrainNbin = 2.; // acceptance for grain preprocessing = eGrainNbin*eGrainSX
240
241 eDoGrainsProcessing = true; // when reconstruct segments do grains preprocessing
242 eCheckSeedThres = false; // if true: use adaptive seeds threshold (based on eNseedMax)
243 ePropagateToBase = false; // if true: segments are propagated to the base position
244
245// set to 0 all pointers
246 eView = 0;
247 eGrainsTree = 0;
248
249 eG = 0; // grains
250 eVC = 0;
251 // eGSeg = 0;
252 eG = 0;
253 eGCla = 0;
254 eCL = 0;
255 eSA = 0;
256
257 epT = 0;
258 epP = 0;
259 epY = 0;
260 ehX = 0;
261 epS = 0;
262 epC = 0;
263
264}
void SetLimits(float xmin, float xmax, float ymin, float ymax, float zmin, float zmax)
Definition: EdbViewRec.h:56
void SetCellLimits(int ncell, int ncl)
Definition: EdbViewRec.h:60
void SetBin(float binx, float biny, float binz=-1)
Definition: EdbViewRec.h:58
Int_t ePulsMax
Definition: EdbViewRec.h:155

◆ SetPulsThres()

void EdbViewRec::SetPulsThres ( int  minp,
int  maxp = 500 
)
inline
195{ ePulsMin=minp; ePulsMax=maxp; }

◆ SetRmax()

void EdbViewRec::SetRmax ( float  rmax)
inline
204{ eRmax=rmax; }

◆ SetSeedsLim()

void EdbViewRec::SetSeedsLim ( int  nseedslim = 100000,
int  seedlim = 48 
)
inline
202{ enSeedsLim=nseedslim; eSeedLim=seedlim; }

◆ SetSeedThres()

void EdbViewRec::SetSeedThres ( int  nt,
int  th[] 
)
inline
209{ if(nt>0){ enT=nt; eSeedThres.Set(nt,th);} };

◆ SetSeedThres0()

void EdbViewRec::SetSeedThres0 ( Short_t  mins)
inline
208{ eSeedThres0=mins;}

◆ SetSigmaThres()

void EdbViewRec::SetSigmaThres ( float  smin,
float  smax 
)
inline
196{ eSigmaMin=smin; eSigmaMax=smax; }

◆ SetStep()

int EdbViewRec::SetStep ( int  sfrom,
int  sto,
int  step = 1 
)
inline
203{ eStepFrom=sfrom; eStepTo=sto; eStep=step; return (sto-sfrom)/step; }

◆ SetThetaLim()

void EdbViewRec::SetThetaLim ( float  t)
inline
205{ eThetaLim=t; }

◆ SetView()

bool EdbViewRec::SetView ( EdbView v)
414{
415 if(!v) return false;
416 eView = v;
417 eCL = v->GetClusters();
418 if(!eCL) return false;
419
420 if ( eView->GetNframesTop()==0 && eView->GetNframesBot()>0 ) { //bottom side
421 //eZmin = eView->GetZ1();
422 //eZmax = eView->GetZ2();
423 eZmin = eView->ZFrameMin();
424 eZmax = eView->ZFrameMax();
425 eZxy = eZmax; // base point
426 }
427 else if( eView->GetNframesBot()==0 && eView->GetNframesTop()>0 ) { //top side
428 //eZmin = eView->GetZ3();
429 //eZmax = eView->GetZ4();
430 eZmin = eView->ZFrameMin();
431 eZmax = eView->ZFrameMax();
432 eZxy = eZmin; // base point
433 } else return false;
434 eZcenter = 0.5*(eZmin+eZmax);
435// if(eAddGrainsToView) {
436// eG = eView->GetSegments();
437// if(!eG) printf("\n***** eG is 0!!!!!!\n\n");
438// printf("***** eG size is is %d \n", eG->GetSize() );
439// if(eG->GetSize()<eNgrMax) eG->Expand(eNgrMax);
440// }
441
442
443 //eVCC.SetNfr( eView->GetNframes(), eZmin, eZmax );
445 eVCG.SetNfr( 16, eZmin, eZmax, 1 );
446
447 eVC = &eVCC;
448 eVC->Print();
449 eVC->FillCell( *(eCL) );
450
451 return true;
452}
void Print()
Definition: EdbViewRec.cxx:153
void SetNfr(int nfr, float zmin, float zmax, int ifz=0)
Definition: EdbViewRec.cxx:73
Int_t GetNframesTop() const
Definition: EdbView.h:207
float ZFrameMin() const
Definition: EdbView.cxx:222
Int_t GetNframesBot() const
Definition: EdbView.h:208
float ZFrameMax() const
Definition: EdbView.cxx:234

◆ SPhiGr()

float EdbViewRec::SPhiGr ( float  theta,
float  phi,
float  dz,
float  sx,
float  sy,
float  sz 
)
404{
405 // cone aperture phi angle for grain of size sx,sy,sz visible at theta,phi,dz
406 float rz = Abs(dz*Tan(theta));
407 float rg = Sqrt( sx*Sin(phi)*sx*Sin(phi) + sy*Cos(phi)*sy*Cos(phi) );
408 if(rz<rg) return TwoPi();
409 return ATan(rg/rz);
410}

◆ SThetaGr()

float EdbViewRec::SThetaGr ( float  theta,
float  phi,
float  dz,
float  sx,
float  sy,
float  sz 
)
395{
396 // cone aperture theta angle for grain of size sx,sy,sz visible at theta,phi,dz
397 return ATan( Sqrt( sz*Sin(theta)*sz*Sin(theta) +
398 sx*Sin(phi)*sx*Sin(phi) + sy*Cos(phi)*sy*Cos(phi) ) /
399 (dz/Cos(theta)) );
400}

Member Data Documentation

◆ eAddGrainsToView

bool EdbViewRec::eAddGrainsToView

◆ eCheckSeedThres

bool EdbViewRec::eCheckSeedThres

◆ eCL

TClonesArray* EdbViewRec::eCL
private

[eNgrMax] array of grains represented as clusters

◆ eClMaxA

Int_t EdbViewRec::eClMaxA
private

◆ eClMinA

Int_t EdbViewRec::eClMinA
private

pointer to array of segments (output of tracking)

◆ eDmax

Float_t EdbViewRec::eDmax
private

pointers to clusters

◆ eDoGrainsProcessing

bool EdbViewRec::eDoGrainsProcessing

◆ eDZmin

Float_t EdbViewRec::eDZmin
private

◆ eFact

Float_t EdbViewRec::eFact
private

◆ eG

TClonesArray* EdbViewRec::eG
private

cells with grains

◆ eGCla

TClonesArray* EdbViewRec::eGCla
private

pointer to eView->GetSegments() or to eGSeg as the output for grain search

◆ eGrainNbin

Float_t EdbViewRec::eGrainNbin

◆ eGrainsTree

TTree* EdbViewRec::eGrainsTree
private

◆ ehX

Short_t* EdbViewRec::ehX
private

[enYtot] - pointers to the first x[iy]

◆ eNclGrMax

int EdbViewRec::eNclGrMax
private

◆ eNclGrMin

int EdbViewRec::eNclGrMin
private

debug tree

◆ eNgr

Int_t EdbViewRec::eNgr
private

◆ eNgrMax

Int_t EdbViewRec::eNgrMax
private

pointer to the input view currently in processing

◆ enP

TArrayI EdbViewRec::enP
private

◆ enPtot

Int_t EdbViewRec::enPtot
private

◆ eNseedMax

TArrayI EdbViewRec::eNseedMax
private

◆ eNseedMax0

Int_t EdbViewRec::eNseedMax0
private

◆ enSeedsLim

Int_t EdbViewRec::enSeedsLim
private

◆ eNsegMax

Int_t EdbViewRec::eNsegMax
private

◆ enT

Int_t EdbViewRec::enT
private

◆ enX

TArrayI EdbViewRec::enX
private

◆ enXtot

Int_t EdbViewRec::enXtot
private

◆ enY

TArrayI EdbViewRec::enY
private

◆ enYtot

Int_t EdbViewRec::enYtot
private

◆ epC

EdbCluster** EdbViewRec::epC
private

pointers to seeds list

◆ epP

Short_t*** EdbViewRec::epP
private

[enT] - pointers to the first phi[it]

◆ ePropagateToBase

bool EdbViewRec::ePropagateToBase

◆ epS

EdbCluster*** EdbViewRec::epS
private

[enXtot] - phase histogram

◆ epT

Short_t**** EdbViewRec::epT
private

◆ ePulsMax

Int_t EdbViewRec::ePulsMax
private

◆ ePulsMin

Int_t EdbViewRec::ePulsMin
private

◆ epY

Short_t** EdbViewRec::epY
private

[enPtot] - pointers to the first y[ip]

◆ eR

TArrayF EdbViewRec::eR
private

◆ eRmax

Float_t EdbViewRec::eRmax
private

◆ eSA

TClonesArray* EdbViewRec::eSA
private

pointer to eView->GetClusters() or to eGCla as the input for tracking

◆ eSeedLim

Short_t EdbViewRec::eSeedLim
private

◆ eSeedThres

TArrayI EdbViewRec::eSeedThres
private

◆ eSeedThres0

Int_t EdbViewRec::eSeedThres0
private

◆ eSigmaMax

Float_t EdbViewRec::eSigmaMax
private

◆ eSigmaMin

Float_t EdbViewRec::eSigmaMin
private

◆ esP

TArrayF EdbViewRec::esP
private

◆ eStep

Int_t EdbViewRec::eStep
private

◆ eStepFrom

Int_t EdbViewRec::eStepFrom
private

◆ eStepTo

Int_t EdbViewRec::eStepTo
private

◆ esX

TArrayF EdbViewRec::esX
private

◆ esY

TArrayF EdbViewRec::esY
private

◆ eTheta

TArrayF EdbViewRec::eTheta
private

◆ eThetaLim

Float_t EdbViewRec::eThetaLim
private

◆ eVC

EdbViewCell* EdbViewRec::eVC
private

◆ eVCC

EdbViewCell EdbViewRec::eVCC
private

pointer to eVCC or eVCG

◆ eVCG

EdbViewCell EdbViewRec::eVCG
private

cells with raw clusters

◆ eView

EdbView* EdbViewRec::eView
private

◆ eZcenter

Float_t EdbViewRec::eZcenter
private

The documentation for this class was generated from the following files: