EdbEDADecaySearch Class Reference

#include <EdbEDADecaySearch.h>

Inheritance diagram for EdbEDADecaySearch:

Collaboration diagram for EdbEDADecaySearch:

Public Types
enum	{ kParentSearch =0x01 , kBaseTrackSearch =0x02 , kDaughterSearch =0x04 }

Public Member Functions
TObjArray *	CheckInTrackKink (EdbTrackP *trk)
TObjArray *	CheckInTrackKinks (TObjArray *tracks=NULL)
TObjArray *	DoSearch ()
	EdbEDADecaySearch (int DSversion=2)
EdbVertex *	FindPrimaryVertex ()
EdbEDATrackP *	FindTrack (EdbTrackP *t)
TObjArray *	FindUpstreamVertices ()
EdbEDADecayVertex *	GetDecayVertex (int i)
float	GetIPCut (EdbVertex *v, EdbSegP *s)
TObjArray *	GetKinks (void)
EdbVertex *	GetPrimaryVertex ()
EdbPVRec *	GetPVR ()
EdbEDATrackP *	GetTrack (int i)
EdbEDATrackSet *	GetTrackSet ()
double	GetTSDauIP1 (double dz)
double	GetTSDauIP2 (double dz)
void	KinkSearch ()
EdbVertex *	MakeFakeVertex (EdbTrackP *t, double dz=1150)
void	MTSearch2seg (TObjArray *tracks)
void	MTSearchAll (TObjArray *tracks)
int	NDecayVertices ()
int	Ntracks ()
void	PrintRunTracking ()
void	PrintTracks ()
void	SetBTSearch (int do_search=1, int npl_up=0, int npl_down=2, float ipcut=20, float phcut=17)
void	SetIPHist1 (TH1F *h)
void	SetIPHist2 (TH1F *h)
void	SetParentSearch (int do_search=1, int npl=2, float ipcut=20, float phcut=17)
void	SetPVR (EdbPVRec *pvr)
void	SetSmallKinkSearch (int do_search=1, int npl_down=5, float Rmin=5.)
void	SetTracks (TObjArray *array)
void	SetTrackSet (EdbEDATrackSet *set)
void	SetVertex (EdbVertex *v)
void	ShortDecaySearch ()
void	SmallKinkSearch ()
TObjArray *	TSBaseTracks (int ipl)
TObjArray *	TSDaughterTracks (TObjArray *base)
void	TSDaughterTracks2 ()
TObjArray *	TSParentTracks (TObjArray *base)
virtual	~EdbEDADecaySearch ()

Public Attributes
int	eBT
double	eBTIP
double	eBTPH
int	eBTPlateDown
int	eBTPlateUp
TObjArray *	eDecayVertices
int	eDSVer
TObjArray *	eKinks
double	eKinkSigma
int	eMT2seg
int	eMTAll
EdbPVRec *	ePVR
TObjArray *	eSegments
EdbEDATrackSet *	eSet
int	eSmallKink
int	eSmallKinkNpl
float	eSmallKinkRmin
TObjArray *	eTracks
int	eTSDau
int	eTSDau2
TH1F *	eTSDauIPHist1
TH1F *	eTSDauIPHist2
int	eTSDauNseg
int	eTSDauNseg2
double	eTSDauPH
int	eTSPar
double	eTSParIP
int	eTSParNseg
double	eTSParPH
int	eTSParPlate
EdbVertex *	eVertex
int	eVertexPlatePID
TObjArray *	eVertices
int	eVtxUp
int	eVtxUpDZ
int	eVtxUpIP

Member Enumeration Documentation

anonymous enum

anonymous enum

Enumerator
kParentSearch
kBaseTrackSearch
kDaughterSearch

        {
        kParentSearch       =0x01,
        kBaseTrackSearch    =0x02,
        kDaughterSearch     =0x04
    };

Constructor & Destructor Documentation

EdbEDADecaySearch()

EdbEDADecaySearch::EdbEDADecaySearch ( int  DSversion = 2 )

inline

                                      :
        eDSVer     (DSversion),
        eVertex    (NULL),
        eTracks    (new TObjArray),
        eSegments  (new TObjArray),
        eVertices  (new TObjArray),
        eSet       (NULL),
        
        eTSDau      (1),
        eTSDauNseg  (2),
        eTSDauPH    (17),
        
        eMT2seg     (0),
        eMTAll      (1),
        
        eTSDau2     (1),
        eTSDauNseg2 (3),
        
        eSmallKink    (0),
        eSmallKinkNpl (5),
        eSmallKinkRmin(5.),
        
        eTSPar     (0),
        eTSParPlate(2),
        eTSParNseg (1),
        eTSParIP   (20),
        eTSParPH   (17),
        
        eBT            (0),
        eBTPlateUp     (0),
        eBTPlateDown   (2),
        eBTIP          (20),
        eBTPH          (17),
        
        eVtxUp   (0),
        eVtxUpDZ (3000),
        eVtxUpIP (500),
        
        eKinkSigma(3.),
        eKinks (new TObjArray),
        eDecayVertices(new TObjArray)
        
        {
            
            float ipcutdau_lt_1 = eDSVer==1 ? 500 : 300;
            float ipcutdau_gt_1 = eDSVer==1 ? 800 : 500;
            
            if(gDirectory->Get("eTSDauIPHist1")) gDirectory->Delete("eTSDauIPHist1");
            eTSDauIPHist1 = new TH1F ("eTSDauIPHist1","IP cut function along dZ, 1st.", 121, -100, 12000);
            // Fill function
            TH1F *h = eTSDauIPHist1;
            h->SetFillColor(kBlue);
            h->SetXTitle("dZ (#mum)");
            h->SetYTitle("IP cut (#mum)");
            for(int i=1;i<=100;i++){
                if(0<=h->GetBinLowEdge(i)&&h->GetBinLowEdge(i)<1000) h->SetBinContent(i, ipcutdau_lt_1);
                else if(1000<=h->GetBinLowEdge(i)&&h->GetBinLowEdge(i)<11700) h->SetBinContent(i, ipcutdau_gt_1);
                else if(11700<=h->GetBinLowEdge(i)) h->SetBinContent(i,0);
            }
            if(gDirectory->Get("eTSDauIPHist1")) gDirectory->Delete("eTSDauIPHist2");
            eTSDauIPHist2 = new TH1F ("eTSDauIPHist2","IP cut function along dZ, 2nd", 121, -100, 12000);
            // Fill function
            h = eTSDauIPHist2;
            h->SetFillColor(kBlue);
            h->SetXTitle("dZ (#mum)");
            h->SetYTitle("IP cut (#mum)");
            for(int i=1;i<=100;i++){
                if(0<=h->GetBinLowEdge(i)&&h->GetBinLowEdge(i)<1000) h->SetBinContent(i, ipcutdau_lt_1);
                else if(1000<=h->GetBinLowEdge(i)&&h->GetBinLowEdge(i)<6000) h->SetBinContent(i, ipcutdau_gt_1);
                else if(6000<=h->GetBinLowEdge(i)) h->SetBinContent(i,0);
            }
        }

~EdbEDADecaySearch()

virtual EdbEDADecaySearch::~EdbEDADecaySearch ( )

inlinevirtual

{}

Member Function Documentation

CheckInTrackKink()

TObjArray * EdbEDADecaySearch::CheckInTrackKink

(

EdbTrackP *

trk

)

                                                            {
    // Search kink. return TObjArray of EdbEDASmallKink
    
    if(trk->N()<2) {
        printf("Track %d nseg=%d, too short for kink search\n", trk->ID(), trk->N());
        return NULL;
    }
    EdbTrackP *t = CleanTrack(trk); // Remove fake segment/microtracks. 
 
    double rms,rmst,rmsl;
    DTRMSTL1Kink(t, &rms, &rmst, &rmsl); // DTRMS for each projection.
 
    if(t->N()<=2) {
        printf("Track %d nseg=%d. too short to calculate delta-theta rms. use 2.0mrad.\n", t->ID(), t->N());
        rms=rmst=rmsl=2.0;
    }
    
    TObjArray *kinks = new TObjArray;
    
    int vtxpl;
    if(gEDA) vtxpl=gEDA->GetIPLZ(eVertex->Z());
    else {
        vtxpl=0; // TODO
        printf("VERTEX plate is not set!!!!! \n");
    }
    
    for(int i=0; i<ePVR->Npatterns(); i++){
        EdbPattern *pat = ePVR->GetPattern(i);
        if( eVertex->Z()<pat->Z()) {
            vtxpl = pat->Plate();
            break;
        }
    }
    
    int n=t->N();
    for(int j=0;j<n-1;j++){
        EdbSegP *s1 = t->GetSegment(j);
        EdbSegP *s2 = t->GetSegment(j+1);
        
        // if s1 is more downstream than 5 plates from vertex, skip. search only kink between 1-x,2-x,3-x,4-x.
        if(s1->Plate()>vtxpl+eSmallKinkNpl) continue;
        
        // calculate kink angle in transverse and longitudinal projection.
        double dtt, dtl;
        CalcDTTransLongi(s1,s2, &dtt, &dtl);
        double dt = sqrt(dtt*dtt+dtl*dtl);
        
        // calculate position difference from 1st segment.
        double dxt, dxl;
        CalcDXTransLongi(t->GetSegmentFirst(),s2, &dxt, &dxl);
        
        if( fabs(dtt)>rmst*eSmallKinkRmin || fabs(dtl)>rmsl*eSmallKinkRmin ) {
            // if there is kink.
            // if kink angle is bigger than 5 times delta-theta rms in one of the 2 projection.
 
            // momentum calculation using downstream from s2.
            double p,pmin,pmax;
            CalcPPartial(t,s2,t->GetSegmentLast(), p, pmin, pmax);
            
            // calculate Pt at the kink, using downstream momentum.
            double Pt = p>0 ? p*dt : -1.;
 
            // calculate vertex point with 2 segments.
            TObjArray segs;
            segs.Add(s1);
            segs.Add(s2);
            EdbVertex *v = CalcVertex(&segs);
            
            int ndau = n-j-1;
            // put the data in a struct.
            EdbEDASmallKink *kink = new EdbEDASmallKink(v, t, s1, s2, dtt, dtl, dxt, dxl, ndau, p, pmin, pmax, Pt, rmst, rmsl);
            
            if(gEDA) gEDA->AddVertex(v);
            kinks->Add(kink); // for this track
            eKinks->Add(kink); // for all tracks.
            
            {
                // Make Decay vertex
                // this is a temporary solution.
                EdbEDADecayVertex *vdecay = new EdbEDADecayVertex;
                vdecay->SetXYZ( v->X(), v->Y(), v->Z());
                vdecay->SetPrimaryVertex(eVertex);
                
                EdbTrackP *tp = new EdbTrackP(*t);
                tp->Clear();
                for(int i=0; i<=j; i++) tp->AddSegment( t->GetSegment(i));
                tp->SetCounters();
                vdecay->SetParent( new EdbEDATrackP(tp, eVertex));
                
                EdbTrackP *td = new EdbTrackP(*t);
                td->Clear();
                for(int i=j+1; i<t->N(); i++) td->AddSegment( t->GetSegment(i));
                td->SetCounters();
                vdecay->SetDaughter( new EdbEDATrackP(td, eVertex));
                vdecay->SetType( EdbEDADecayVertex::kSmallKink | EdbEDADecayVertex::kLong);
                eDecayVertices->Add(vdecay);
            }
            
            
            printf("Kink candidate. itrk %d plate %d - %d kink angle %.4lf P %.3lf Pt %.3lf ", 
                t->ID(), s1->PID(), s2->PID(), dt, p, Pt);
            printf("(Trans, Longi) = ( %.4lf, %.4lf ) thr ( %.4lf, %.4lf )\n",
                dtt, dtl, rmst*3, rmsl*3);
            //gEDA->AddDrawObject(kink);
        }
    }
    
    printf("%d kink candidates are found.\n", kinks->GetEntriesFast());
    
    return kinks;
}

CheckInTrackKinks()

TObjArray* EdbEDADecaySearch::CheckInTrackKinks ( TObjArray *  tracks = NULL )

inline

                                                           { 
        if(NULL==tracks) tracks=eTracks; 
        for(int i=0;i<tracks->GetEntriesFast();i++) {
            CheckInTrackKink((EdbTrackP *) tracks->At(i)); 
        }
        return eKinks;
    }

DoSearch()

TObjArray * EdbEDADecaySearch::DoSearch

(

)

                                       {
    // Do track search
    
    
    // clear previous results.
    eTracks->Clear();
    eSegments->Clear();
    
    int nAll=-1;
    int nTSDau1=-1;
    int nTSDau2=-1;
    if(gEve) gEve->SetStatusLine("Start Track search");
    printf("######### Track Search #########\n");
    if(eVertex==NULL) {
        printf("Vertex is not selected!! End.\n");
        return eTracks;
    }
    if(ePVR==NULL){
        printf("No PVRec set\n");
        return eTracks;
    }
    
    int hall[11],hdau1[11];
    int hdau2[11],hmt[11];
    for(int i=0;i<11;i++) hall[i]=hdau1[i]=hdau2[i]=hmt[i]=0;
    
    printf(" target vertex = ( %8.1f, %8.1f, %8.1f) on pid%d.\n", eVertex->X(), eVertex->Y(), eVertex->Z(), eVertexPlatePID);
    printf("  with %d track\n", ePVR->Ntracks());
    
    TObjArray *base = ePVR->eTracks;
    
    for(int i=0;i<base->GetEntriesFast();i++) hall[((EdbTrackP *) base->At(i))->N()>=10?10:((EdbTrackP *) base->At(i))->N()]++;
    nAll=base->GetEntriesFast();
    
    // Daughter search 1st
    if(eTSDau) TSDaughterTracks(base);
    
    for(int i=0;i<eTracks->GetEntriesFast();i++) {
        EdbTrackP *t = ((EdbEDATrackP *) eTracks->At(i))->GetOriginal();
        int nseg= t->N();
        if( nseg>10) nseg=10;
        hdau1[nseg]++;
    }
    nTSDau1 = eTracks->GetEntriesFast();
    printf("%d\n", nTSDau1);
    
    // Microtrack search
    if(eMT2seg&&!eMTAll) MTSearch2seg(eTracks);
    if(eMTAll) {
        MTSearchAll(eTracks);
    }
    for(int i=0;i<eTracks->GetEntriesFast();i++) hmt[((EdbEDATrackP *) eTracks->At(i))->GetOriginal()->N()>=10?10:((EdbEDATrackP *) eTracks->At(i))->GetOriginal()->N()]++;
    
    // Daughter cut 2nd.
    if(eTSDau2){
        TObjArray tmp = *eTracks;
        TSDaughterTracks2();
    }
    for(int i=0;i<eTracks->GetEntriesFast();i++) hdau2[((EdbEDATrackP *) eTracks->At(i))->N()>=10?10:((EdbEDATrackP *) eTracks->At(i))->GetOriginal()->N()]++;
 
    nTSDau2 = eTracks->GetEntriesFast();
    
    // small kin search
    if(eSmallKink) SmallKinkSearch();
    
    if(eTSPar) TSParentTracks(base);
    
    if(eBT){
        printf(" Base track search. ");
        printf("Upstream %d plates, Downstream %d plates ip<%.1lf w>%.1lf\n", 
                eBTPlateUp, eBTPlateDown, eBTIP, eBTPH);
        for(int dPlate=-eBTPlateUp+1; dPlate<=eBTPlateDown; dPlate++){
            int pid = eVertexPlatePID + dPlate -1;
            if(pid<0||ePVR->Npatterns()-1<pid) continue;
            TSBaseTracks(pid);
        }
    }
    
    PrintTracks();
    
    printf("------  SUMMARY  --------------\n");
    printf("                         nseg ");
    for(int i=1;i<11;i++) printf("%3d ", i);
    printf("\n");
    printf(" Total            %4d tracks ", nAll);
    for(int i=1;i<11;i++) printf("%3d ", hall[i]);
    printf("\n");
    printf(" Daughter search  %4d tracks ", nTSDau1);
    for(int i=1;i<11;i++) printf("%3d ", hdau1[i]);
    printf("\n");
    printf(" Microtrack search            ");
    for(int i=1;i<11;i++) printf("%3d ", hmt[i]);
    printf("\n");
    printf(" Nseg>=3 (inc.MT) %4d tracks ", nTSDau2);
    for(int i=1;i<11;i++) printf("%3d ", hdau2[i]);
    printf("\n");
    printf("-------------------------------\n");
    
    return eTracks;
}

FindPrimaryVertex()

EdbVertex * EdbEDADecaySearch::FindPrimaryVertex

(

)

                                                {
    // Find Best Vertex.
    // select best one from pvr->eVertex
    
    printf("FindBestVertex(): \n");
    if( ePVR->Nvtx()==0 ){
        printf("No vertex is set. make a fake vertex.");
        
        EdbTrackP *tup=NULL;
        for(int i=0; i<ePVR->Ntracks(); i++){
            EdbTrackP *t = ePVR->GetTrack(i);
            if(t->N()<3) continue;
            if(EdbEDAUtil::DTRMS(t)>0.02) continue;
            if(tup==NULL) tup=t;
            if(t->Z()<tup->Z()) tup=t;
        }
        if(tup==NULL) return NULL;
        
        EdbVertex *v = MakeFakeVertex( tup, -650);
        
        return v;
    }
    
    // select most upstream vertex
    int zero=0;
    EdbVertex *v_most_upstream = ePVR->GetVertex(zero);
    for(int i=1; i<ePVR->Nvtx(); i++){
        EdbVertex *v = ePVR->GetVertex(i);
        if(v->Z()<v_most_upstream->Z()) v_most_upstream=v;
    }
    
    // gather all vertex near the most upstream vertex
    // vertices within 200 micron in Z.
    
    printf("TODO: FindPrimaryVertex... find real vertex, especially QE like kink event.\n");
    printf("memo: currently mostly 2 track vertex is selected.");
    printf("TODO: use SB info.\n");
    printf("TODO: constraint for vertex Z < SB stop Z\n");
    TObjArray *vertices = new TObjArray;
    for(int i=0; i<ePVR->Nvtx(); i++){
        EdbVertex *v = ePVR->GetVertex(i);
        if(v->Z()<v_most_upstream->Z()+200) vertices->Add(v);
    }
    
    // select most high probability vertex
    EdbVertex *v_highest_prob= (EdbVertex *) vertices->At(0);
    for(int i=1; i<vertices->GetEntriesFast(); i++){
        EdbVertex *v= (EdbVertex *) vertices->At(i);
        if(v->V()->prob()>v_highest_prob->V()->prob()) v_highest_prob = v;
    }
    
    
    printf("select upstream and high probable vertex: %d vertex\n", ePVR->Nvtx());
    for(int i=0; i<ePVR->Nvtx();i++){
        EdbVertex *v = ePVR->GetVertex(i);
        int flag = 0;
        for(int j=0; j<vertices->GetEntriesFast(); j++) { if(v==vertices->At(j)) flag=1;}
        printf("vertex %2d, %2d tracks z=%8.1f prob=%7.5f %s %s, %s\n", i, v->N(), v->Z(), v->V()->prob(), 
                v==v_most_upstream?"most":"    ", 
                flag?"upstream": "         ", 
                v==v_highest_prob?"highest prob":"");
    }
    
    delete vertices;
    return v_highest_prob;
    //return v_most_upstream;
    
}

FindTrack()

EdbEDATrackP* EdbEDADecaySearch::FindTrack ( EdbTrackP *  t )

inline

                                         {
        for(int i=0;i<Ntracks();i++){
            if(GetTrack(i)->ID()==t->ID()) return GetTrack(i);
        }
        return NULL;
    }

FindUpstreamVertices()

TObjArray * EdbEDADecaySearch::FindUpstreamVertices

(

)

                                                  {
    
    printf("Find upstream vertices. Nvtx_base = %d\n", ePVR->Nvtx());
    TObjArray *vertices= new TObjArray;
    
    for(int i=0;i<ePVR->Nvtx();i++){
        EdbVertex *v = ePVR->GetVertex(i);
        double dz = v->Z() - eVertex->Z();
        if( dz < -eVtxUpDZ || 0 < dz ) continue;
        
        int flag=0;
        for(int j=0;j<eVertex->N();j++){
            EdbSegP *s = eVertex->GetTrack(j)->GetSegmentFirst();
            if(CalcIP(s,v)<eVtxUpIP) flag++;
        }
        if(flag ==0) continue;
        vertices->Add(v);
    }
    
    return vertices;
}

GetDecayVertex()

EdbEDADecayVertex* EdbEDADecaySearch::GetDecayVertex ( int  i )

inline

{ return (EdbEDADecayVertex *)eDecayVertices->At(i);}

GetIPCut()

float EdbEDADecaySearch::GetIPCut	(	EdbVertex *	v,
		EdbSegP *	s
	)

                                                         {
    // calculate IP cut for given 1ry vertex and segment.
    // if segment is BT, Z of segment is assumed Z=s->Z()-150;
    // otherwise, use s->Z();
    
    float Z = s->Side()==0? s->Z()-150 : s->Z();
    float dZ = fabs( Z-v->Z());
    return 5.0+0.01*dZ;
}

GetKinks()

TObjArray* EdbEDADecaySearch::GetKinks ( void  )

inline

{ return eKinks;}

GetPrimaryVertex()

EdbVertex* EdbEDADecaySearch::GetPrimaryVertex ( )

inline

{ return eVertex;}

GetPVR()

EdbPVRec* EdbEDADecaySearch::GetPVR ( )

inline

{ return ePVR;}

GetTrack()

EdbEDATrackP* EdbEDADecaySearch::GetTrack ( int  i )

inline

{ return (EdbEDATrackP *) eTracks->At(i);}

GetTrackSet()

EdbEDATrackSet* EdbEDADecaySearch::GetTrackSet ( )

inline

{ return eSet;}

GetTSDauIP1()

double EdbEDADecaySearch::GetTSDauIP1

(

double

dz

)

                                              {
    TH1F *h = eTSDauIPHist1;
    
    if(dz<h->GetXaxis()->GetXmin()) return h->GetBinContent(1);
    if(dz>h->GetXaxis()->GetXmax()) return h->GetBinContent(h->GetNbinsX());
    for(int i=1;i<h->GetNbinsX();i++){
        if(dz>=h->GetBinLowEdge(i)&&dz<h->GetBinLowEdge(i)+h->GetBinWidth(i)) return h->GetBinContent(i);
    }
    return 0.0;
}

GetTSDauIP2()

double EdbEDADecaySearch::GetTSDauIP2

(

double

dz

)

                                              {
    TH1F *h = eTSDauIPHist2;
    
    if(dz<h->GetXaxis()->GetXmin()) return h->GetBinContent(1);
    if(dz>h->GetXaxis()->GetXmax()) return h->GetBinContent(h->GetNbinsX());
    for(int i=1;i<h->GetNbinsX();i++){
        if(dz>=h->GetBinLowEdge(i)&&dz<h->GetBinLowEdge(i)+h->GetBinWidth(i)) return h->GetBinContent(i);
    }
    return 0.0;
}

KinkSearch()

void EdbEDADecaySearch::KinkSearch

(

)

                                  {
    printf("------------------------------\n");
    printf("Kink search in selected tracks\n");
    for(int i=0; i<Ntracks(); i++){
        for(int j=0; j<Ntracks(); j++){
            
            if(i==j) continue;
            
            EdbEDATrackP *t1 = GetTrack(i);
            EdbEDATrackP *t2 = GetTrack(j);
            
            EdbSegP *s1 = t1->GetSegmentLast();
            EdbSegP *s2 = t2->GetSegmentFirst();
            
            if(s1->Plate() >= s2->Plate()) continue;
            // if last segment of t1 is upstream of first segment of t2.
            
            // ignore kink angle < 20 mrad.
            float kinkangle = sqrt((s1->TX()-s2->TX())*(s1->TX()-s2->TX())+(s1->TY()-s2->TY())*(s1->TY()-s2->TY()));
            if( kinkangle<0.02) continue;
            
            // remove very low momentum
            double rmss, rmst, rmsl;
            DTRMSTL(t2, &rmss, &rmst, &rmsl);
            if( rmst>0.015 && rmsl>0.015) continue;
            
            TObjArray segs;
            segs.Add(s1);
            segs.Add(s2);
            EdbVertex *v = EdbEDAUtil::CalcVertex(&segs);
            float ip1 = EdbEDAUtil::CalcIP(s1, v);
            float ip2 = EdbEDAUtil::CalcIP(s2, v);
            
            // cut with IP and Z. (Z should between 2 segments.)
            if(ip1<20 && ip2<20 && s1->Z()-500<v->Z() && v->Z()<s2->Z()+500){
                printf("Decay vertex t1 %d t2 %d ip1 %f ip2 %f\n", t1->ID(), t2->ID(), ip1, ip2);
                EdbEDADecayVertex *dv = new EdbEDADecayVertex;
                dv->SetXYZ( v->X(), v->Y(), v->Z());
                dv->SetParent(t1);
                dv->SetDaughter(t2);
                dv->SetType( EdbEDADecayVertex::kLong);
                printf("kink found, %d %f %f %f\n", dv->N(), dv->X(), dv->Y(), dv->Z());
                eDecayVertices->Add(dv);
            }
            
            delete v;
        }
    }
}

MakeFakeVertex()

EdbVertex * EdbEDADecaySearch::MakeFakeVertex	(	EdbTrackP *	t,
		double	dz = 1150
	)

                                                                   {
    EdbSegP *s0 = t->GetSegmentFirst();
    
    EdbSegP s(*s0);
    s.PropagateTo(s0->Z()+dz);
 
    EdbVertex *v = new EdbVertex;
    v->SetID(t->ID());
    v->SetXYZ(s.X(),s.Y(),s.Z());
    
    EdbVTA *vta = new EdbVTA(t, v);
    vta->SetZpos(1);
    vta->SetFlag(2);
    vta->SetImp(EdbEDAUtil::CalcIP(&s,v));
    v->AddVTA(vta);
    return v;
}

MTSearch2seg()

void EdbEDADecaySearch::MTSearch2seg

(

TObjArray *

tracks

)

                                                     {
    // Microtrack search for 2 segment track.
    // if eSet is not set, do nothing.
    
    // if eSet==NULL or Preparation is failed, do nothing.
    int ret = eSet ? eSet->PrepareScanSetForMT() : -1; 
    if(ret==-1) {
        printf("ScanSet error. stop.\n");
        return;
    }
    for(int i=0;i<tracks->GetEntriesFast();i++){
        EdbTrackP *t = ((EdbEDATrackP *) tracks->At(i))->GetOriginal();
        if(t->N()==2) eSet->MicroTrackSearch(t);
    }
}

MTSearchAll()

void EdbEDADecaySearch::MTSearchAll

(

TObjArray *

tracks

)

                                                    {
    // Microtrack search for all tracks.
    // if eSet is not set, do nothing.
    
    // if eSet==NULL or Preparation is failed, do nothing.
    int ret = eSet ? eSet->PrepareScanSetForMT() : -1; 
    if(ret==-1) {
        printf("ScanSet error. stop.\n");
        return;
    }
    for(int i=0;i<tracks->GetEntriesFast();i++){
        gEve->SetStatusLine(Form("MicroTrackSearch %d/%d", i+1, tracks->GetEntriesFast()));
        EdbTrackP *t = ((EdbEDATrackP *) tracks->At(i))->GetOriginal();
        eSet->MicroTrackSearch(t);
    }
}

NDecayVertices()

int EdbEDADecaySearch::NDecayVertices ( )

inline

{ return eDecayVertices->GetEntriesFast();}

Ntracks()

int EdbEDADecaySearch::Ntracks ( )

inline

{ return eTracks->GetEntriesFast();}

PrintRunTracking()

void EdbEDADecaySearch::PrintRunTracking ( )

inline

                           {
        if(eSet) eSet->eRunTracking.Print();
        else printf("No EdbEDATrackSet is set\n");
    }

PrintTracks()

void EdbEDADecaySearch::PrintTracks

(

)

                                   {
    printf("ibrick trkid  flag pl ns x        y        z        tx      ty      pl pa cs dk ip1ry class PBDLU comment\n");
    for(int i=0;i<Ntracks();i++){
        EdbEDATrackP *t = GetTrack(i);
        t->Print();
    }
}

SetBTSearch()

void EdbEDADecaySearch::SetBTSearch ( int  do_search = 1, int  npl_up = 0, int  npl_down = 2, float  ipcut = 20, float  phcut = 17  )

inline

                                                                                                    {
        eBT = do_search; eBTPlateUp = npl_up; eBTPlateDown = npl_down; eBTIP = ipcut; eBTPH = phcut;}

SetIPHist1()

void EdbEDADecaySearch::SetIPHist1 ( TH1F *  h )

inline

{ eTSDauIPHist1=h;}

SetIPHist2()

void EdbEDADecaySearch::SetIPHist2 ( TH1F *  h )

inline

{ eTSDauIPHist2=h;}

SetParentSearch()

void EdbEDADecaySearch::SetParentSearch ( int  do_search = 1, int  npl = 2, float  ipcut = 20, float  phcut = 17  )

inline

                                                                                     {
        eTSPar = do_search; eTSParPlate = npl; eTSParIP = ipcut; eTSParPH = phcut;}

SetPVR()

void EdbEDADecaySearch::SetPVR ( EdbPVRec *  pvr )

inline

{ ePVR = pvr;}

SetSmallKinkSearch()

void EdbEDADecaySearch::SetSmallKinkSearch ( int  do_search = 1, int  npl_down = 5, float  Rmin = 5.  )

inline

                                                                            {
        eSmallKink = do_search; eSmallKinkNpl = npl_down; eSmallKinkRmin = Rmin;}

SetTracks()

void EdbEDADecaySearch::SetTracks ( TObjArray *  array )

inline

{eTracks=array;}

SetTrackSet()

void EdbEDADecaySearch::SetTrackSet ( EdbEDATrackSet *  set )

inline

                                         { 
        eSet = set;
        SetPVR(set->GetPVRec());
    }

SetVertex()

void EdbEDADecaySearch::SetVertex

(

EdbVertex *

v

)

                                               { 
    eVertex = v;
    if(NULL==v) return;
    double vertexplz = 1e9;
    
    printf("vertex z = %f\n", v->Z());
    for(int i=0;i<ePVR->Npatterns();i++){
        EdbPattern *pat = ePVR->GetPattern(i);
        if(pat==NULL) continue;
        if(pat->Z()<v->Z()) continue;
        if(pat->Z()<vertexplz) {
            vertexplz    = pat->Z();
            eVertexPlatePID = pat->ID();   // Set eVertexPlatePID
        }
    }
    printf("VertexPlatePID=%d\n", eVertexPlatePID);
}

ShortDecaySearch()

void EdbEDADecaySearch::ShortDecaySearch

(

)

                                        {
    printf(" --------------------------\n");
    printf("// Short decay search    //\n");
    printf("-------------------------- \n");
    printf("Using tracks which start from vertex plate and Nseg>=3\n");
    printf("TODO: Microtrack search!!\n");
    
    // select tracks in vertex plate
    TObjArray tracks;
    for(int i=0; i<Ntracks(); i++){
        EdbEDATrackP *t = GetTrack(i);
        EdbSegP *s1 = t->GetSegmentFirst();
        
        // Nseg cut
        if(t->N()<=2) continue;
        
        // use only tracks which start from vertex plate.
        if(s1->PID()!=eVertexPlatePID) continue;
        
        tracks.Add(t);
    }
    printf("Tracks selection for short decay : %d -> %d\n", Ntracks(), tracks.GetEntriesFast());
    
    if(tracks.GetEntriesFast()==0) return;
    
    // Calculate vertex with all tracks.  First time.
    TObjArray segs;
    for(int i=0; i<tracks.GetEntriesFast(); i++) segs.Add( ((EdbEDATrackP *)tracks.At(i))->GetSegmentFirst());
    EdbVertex *v = CalcVertex(&segs);
    float ipcut = GetIPCut(v, (EdbSegP *)segs.At(0));
    printf("Vertex with all tracks: %.1f %.1f %.1f\n", v->X(), v->Y(), v->Z());
    printf("dZ=%.1f microns -> IP cut %.1f\n", ((EdbSegP *)segs.At(0))->Z()-v->Z(), ipcut);
    
    int flag=0;
    // check if short decay candidate exists or not.
    
    if(tracks.GetEntriesFast()==2){ // 2 track case.
        EdbEDATrackP *t1 = (EdbEDATrackP *) tracks.At(0);
        EdbSegP *s1 = t1->GetSegmentFirst();
        EdbEDATrackP *t2 = (EdbEDATrackP *) tracks.At(1);
        EdbSegP *s2 = t2->GetSegmentFirst();
        float ip = CalcDmin(s1,s2);
        printf("track %5d x track %5d, dmin %4.1f %s\n", t1->ID(), t2->ID(), ip, ip>ipcut?">ipcut":"");
        if(ip>ipcut) flag=2;
    }
    else { // more tracks case
        float ipmean=0.0;
        for(int i=0; i<tracks.GetEntriesFast(); i++){
            EdbEDATrackP *t = (EdbEDATrackP *) tracks.At(i);
            EdbSegP *s = t->GetSegmentFirst();
            float ip = CalcIP( s, v);
            
            printf("track %5d (%7.4f, %7.4f) ip %4.1f %s\n", t->ID(), t->TX(), t->TY(), ip, ip>ipcut?">ipcut":"");
            ipmean +=ip;
            if(ip>ipcut) flag++;
        }
        ipmean /= tracks.GetEntriesFast();
        printf("IP mean %.1lf microns with %d tracks\n", ipmean, tracks.GetEntriesFast());
    }
    if(flag==0) {
        printf("No short decay. stop.\n");
        return;
    }
    printf("%d short decay candidates\n", flag);
    
    printf("Vertex refinement by removing one or more tracks.\n");
    
    while(flag&&tracks.GetEntriesFast()!=2){
        // loop to make vertex removing one track.
        int imax=-1;
        float ipmax = 0;
        for(int i=0; i<tracks.GetEntriesFast(); i++){
            EdbEDATrackP *t = (EdbEDATrackP *) tracks.At(i);
            EdbSegP *s = t->GetSegmentFirst();
            
            TObjArray segs2 = segs;
            segs2.Remove( s);
            segs2.Sort();
            
            EdbVertex *v2 = CalcVertex(&segs2);
            float ip = CalcIP(s, v2);
            
            if(ip>ipmax) {
                imax=i;
                ipmax=ip;
            }
        }
        
        // if ipmax>ipcut.
        if(ipmax>ipcut) {
            EdbEDATrackP *tsd = (EdbEDATrackP *) tracks.At(imax);
            
            // Calculate new vertex position.
            tracks.Remove(tsd);
            tracks.Sort();
            segs.Clear();
            for(int i=0; i<tracks.GetEntriesFast(); i++) segs.Add( ((EdbEDATrackP *)tracks.At(i))->GetSegmentFirst());
            EdbVertex *v2 = CalcVertex(&segs);
            
            // Set as Primary vertex.
            eVertex = v2;
            
            // remove very low momentum
            double rmss, rmst, rmsl;
            DTRMSTL(tsd, &rmss, &rmst, &rmsl);
            if( rmst>0.015 && rmsl>0.015) continue;
            
            // Make Short Decay vertex
            EdbVertex *vsd = MakeFakeVertex( tsd, -150);
            EdbEDADecayVertex *vdecay = new EdbEDADecayVertex;
            vdecay->SetXYZ( vsd->X(), vsd->Y(), vsd->Z());
            vdecay->SetPrimaryVertex(v2);
            vdecay->SetDaughter( tsd);
            vdecay->SetType( EdbEDADecayVertex::kShort);
 
            eDecayVertices->Add(vdecay);
            
            
            printf("Short decay of track %d is found. ip = %5.1f.  \n", tsd->ID(), ipmax);
            
            // Check if no other tracks has ip>ipcut
            flag=0;
            if(tracks.GetEntriesFast()==2){
                EdbEDATrackP *t1 = (EdbEDATrackP *) tracks.At(0);
                EdbSegP *s1 = t1->GetSegmentFirst();
                EdbEDATrackP *t2 = (EdbEDATrackP *) tracks.At(1);
                EdbSegP *s2 = t2->GetSegmentFirst();
                float ip = CalcDmin(s1,s2);
                printf("track %5d x track %5d, dmin %4.1f %s\n", t1->ID(), t2->ID(), ip, ip>ipcut?">ipcut":"");
                if(ip>ipcut) flag=2;
            }
            else {
                float ipmean=0.0;
                for(int i=0; i<segs.GetEntriesFast(); i++){
                    EdbSegP *s = (EdbSegP *) segs.At(i);
                    float ip = CalcIP(s, v2);
                    if(ip>ipcut) flag++;
                    ipmean += ip;
                }
                ipmean /= segs.GetEntriesFast();
                printf("IP mean %.1lf microns with %d tracks, removing track %d\n",  ipmean, tracks.GetEntriesFast(), tsd->ID());
            }
        } else flag=0;
        
    }
    
    
    if( flag && tracks.GetEntriesFast()==2){
        // Short decay with only 1 partner track. (without vertex)
        for(int i=0; i<2; i++){
            EdbEDATrackP *t1 = (EdbEDATrackP *) tracks.At(i);
            EdbSegP *s1 = t1->GetSegmentFirst();
            EdbEDATrackP *t2 = (EdbEDATrackP *) tracks.At(!i);
            EdbSegP *s2 = t2->GetSegmentFirst();
            
            // remove very low momentum
            double rmss, rmst, rmsl;
            DTRMSTL(t1, &rmss, &rmst, &rmsl);
            if( rmst>0.015 && rmsl>0.015) continue;
            
            EdbVertex *vsd = MakeFakeVertex( t1, -150);
            EdbEDADecayVertex *vdecay = new EdbEDADecayVertex;
            vdecay->SetXYZ( vsd->X(), vsd->Y(), vsd->Z());
            vdecay->SetDaughter( t1);
            vdecay->SetPartner( t2);
            vdecay->SetType(EdbEDADecayVertex::kShort);
            eDecayVertices->Add(vdecay);
            
        }
    }
    
    for(int i=0; i<NDecayVertices(); i++){
        EdbEDADecayVertex *v = GetDecayVertex(i);
        if(v->GetPrimaryVertex()) v->SetPrimaryVertex(eVertex);
    }
}

SmallKinkSearch()

void EdbEDADecaySearch::SmallKinkSearch

(

)

                                       {
    for(int i=0;i<Ntracks();i++){
        EdbEDATrackP *t = GetTrack(i);
        TObjArray *kinks=CheckInTrackKink(t);
        if(kinks->GetEntriesFast()!=0) {
            if(gEDA==NULL) continue;
            EdbEDAPlotTab *p = gEDA->GetPlotTab();
            if(p!=NULL) p->CheckKink(t);
        }
    }
}

TSBaseTracks()

TObjArray * EdbEDADecaySearch::TSBaseTracks

(

int

ipl

)

                                                  {
    // Search Tracks and set eTracks. return TObjArray filled with selected segments.
    TObjArray *filtered = new TObjArray();
    printf("  pid=%2d ", pid);
    EdbPattern *pat = ePVR->GetPattern(pid);
    
    if(pat==NULL) {
        printf("No pattern is available for pid=%d\n", pid);
        return filtered;
    }
 
    int found=0;
    int found_as_track=0;
    for(int i=0;i<pat->N();i++){
        EdbSegP *s = pat->GetSegment(i);
        if(s==NULL) continue;
        
        // IP cut
        double ip=CalcIP(s, eVertex);
        if( ip > eBTIP )  continue;
        
        // PH cut
        if( s->W() < eBTPH ) continue;
        
        found++;
        // already found as tracks to be drawn?
        int flag=0;
        for(int j=0;j<eTracks->GetEntriesFast();j++){
            EdbTrackP *t = (EdbTrackP *)eTracks->At(j);
            for(int k=0;k<t->N();k++){
                EdbSegP *ss = t->GetSegment(k);
                if(ss->PID()==s->PID()&&ss->ID()==s->ID()){
                    flag++;
                    found_as_track++;
                }
            }
        }
        if(flag) continue;
        
        // add to array
        filtered->Add(s);
        eSegments->Add(s);
        // add to array
        EdbTrackP *t = new EdbTrackP(s);
        t->SetCounters();
 
        EdbEDATrackP *t2 = new EdbEDATrackP(t, eVertex);
        t2->SetTrackSearch( t2->TrackSearch()|kBaseTrackSearch);
        t2->SetIPUp(ip);
        eTracks->Add(t2);
 
        printf("   Track id=%6d ip=%5.1lf pid%02d->%02d %7.4lf %7.4lf\n", t->ID(), ip, 
            t->GetSegmentFirst()->PID(),
            t->GetSegmentLast()->PID(),
            t->TX(),t->TY());
    }
    printf("%5d segments from pid %2d. (%d found, %d already selected as track)\n", filtered->GetEntriesFast(), pid, found, found_as_track);
    return filtered;
}

TSDaughterTracks()

TObjArray * EdbEDADecaySearch::TSDaughterTracks

(

TObjArray *

base

)

                                                              {
    // Search Tracks and set eTracks. return TObjArray filled with selected tracks.
    printf(" Daughter track search. ");
    printf("nseg>=%d, ip=func of dz, w>%.1lf\n", eTSDauNseg,eTSDauPH);
//  printf("nseg>=%d, start within %d plates, ip<%.1lf w>%.1lf\n", 
//          eTSDauNseg, eTSDauPlate, 0.0 /*eTSDauIP*/, eTSDauPH);
    
    TObjArray *filtered = new TObjArray;
    
    int inew=0;
    
    for(int i=0;i<base->GetEntriesFast();i++){
        EdbTrackP *t = (EdbTrackP *) base->At(i);
        if(t==NULL) continue;
        EdbSegP *s = t->GetSegmentFirst();
        
        // Nseg cut
        if( t->N() < eTSDauNseg ) continue;
 
        // dPlate cut
        //int dPlate = s->PID() - eVertexPlatePID + 1; // 
        //if( dPlate < 1 || eTSDauPlate < dPlate ) continue; // 0 <= dPlate <= eTSDauPlate
        
        // IP cut
        double ip = CalcIP(s, eVertex);
        
        double dz = s->Z() - eVertex->Z();
        
        if( ip > GetTSDauIP1(dz) )  continue;
        
        // PH cut
        double ph = t->Wgrains()/t->N();
        if( ph < eTSDauPH ) continue;
        
 
        // add to array
        filtered->Add(t);
        EdbEDATrackP *t2 = FindTrack(t);
        if(t2) t2->SetTrackSearch( t2->TrackSearch()|kDaughterSearch);
        else {
            inew++;
            printf("   Track id=%6d ip=%5.1lf pid%02d->%02d %7.4lf %7.4lf\n", t->ID(), ip, 
            t->GetSegmentFirst()->PID(),
            t->GetSegmentLast()->PID(),
            t->TX(),t->TY());
            
            t2 = new EdbEDATrackP(t,eVertex);
            t2->SetTrackSearch( t2->TrackSearch()|kDaughterSearch);
            t2->SetIPUp(ip);
            eTracks->Add(t2);
        }
    }
    printf("%5d tracks are selected (new %2d tracks)\n", filtered->GetEntriesFast(), inew);
    return filtered;
}

TSDaughterTracks2()

void EdbEDADecaySearch::TSDaughterTracks2

(

)

                                         {
    // Filter 2segment tracks after microtrack search.
    // sevior dz cut will be applied.
    
    printf(" Daughter track search2. ");
    printf("nseg>=%d, ip=func of dz\n", eTSDauNseg2);
    
    int ntrk=eTracks->GetEntriesFast();
    for(int i=0;i<ntrk;i++){
        EdbEDATrackP *tt = (EdbEDATrackP *) eTracks->At(i);
        EdbTrackP *t = tt->GetOriginal();
        if(t==NULL) continue;
        
        EdbSegP *s = t->GetSegmentFirst();
        
        // Nseg cut
        // IP cut
        double ip = CalcIP(s, eVertex);
        double dz = s->Z() - eVertex->Z();
        
        if( ip < GetTSDauIP2(dz) && t->N() >= eTSDauNseg2)  continue;
        
        // remove this track
        eTracks->Remove(tt);
    }
    eTracks->Sort();
    printf("%5d tracks are remained. original %d tracks\n", eTracks->GetEntriesFast(), ntrk);
}

TSParentTracks()

TObjArray * EdbEDADecaySearch::TSParentTracks

(

TObjArray *

base

)

                                                            {
    // Search Tracks and set eTracks. return TObjArray filled with selected tracks.
    
    printf(" Parent track search. ");
    printf("nseg>=%d, start within %d plates, ip<%.1lf w>%.1lf\n", 
            eTSParNseg, eTSParPlate, eTSParIP, eTSParPH);
    
    printf("eVertexPlatePID=%d\n", eVertexPlatePID);
    
    TObjArray *filtered = new TObjArray;
    int inew=0;
    for(int i=0;i<base->GetEntriesFast();i++){
        EdbTrackP *t = (EdbTrackP *) base->At(i);
        if(t==NULL) continue;
        EdbSegP *s = t->GetSegmentFirst();
        
        // Nseg cut
        if( t->N() < eTSParNseg ) continue;
 
        // dPlate cut
        int dPlate = s->PID() - eVertexPlatePID + 1;
        if( dPlate <= -eTSParPlate || eTSParPlate < dPlate ) continue; // 0 <= dPlate <= eTSParPlate
        
        // IP cut
        double ip = CalcIP(s, eVertex);
        if( ip > eTSParIP )  continue;
        
        // PH cut
        double ph = t->Wgrains()/t->N();
        if( ph < eTSParPH ) continue;
        
 
        // add to array
        filtered->Add(t);
        EdbEDATrackP *t2 = FindTrack(t);
        if(t2) t2->SetTrackSearch( t2->TrackSearch()|kParentSearch);
        else {
            inew++;
            printf("   Track id=%6d ip=%5.1lf pid%02d->%02d %7.4lf %7.4lf\n", t->ID(), ip, 
                t->GetSegmentFirst()->PID(),
                t->GetSegmentLast()->PID(),
                t->TX(),t->TY());
            EdbEDATrackP *t2 = new EdbEDATrackP( t, eVertex);
            t2->SetTrackSearch( t2->TrackSearch()|kParentSearch);
            eTracks->Add(t2);
        }
    }
    
    printf("%5d tracks are selected (new %2d tracks)\n", filtered->GetEntriesFast(), inew);
    return filtered;
}

Member Data Documentation

eBT

int EdbEDADecaySearch::eBT

eBTIP

double EdbEDADecaySearch::eBTIP

eBTPH

double EdbEDADecaySearch::eBTPH

eBTPlateDown

int EdbEDADecaySearch::eBTPlateDown

eBTPlateUp

int EdbEDADecaySearch::eBTPlateUp

eDecayVertices

TObjArray* EdbEDADecaySearch::eDecayVertices

eDSVer

int EdbEDADecaySearch::eDSVer

eKinks

TObjArray* EdbEDADecaySearch::eKinks

eKinkSigma

double EdbEDADecaySearch::eKinkSigma

eMT2seg

int EdbEDADecaySearch::eMT2seg

eMTAll

int EdbEDADecaySearch::eMTAll

ePVR

EdbPVRec* EdbEDADecaySearch::ePVR

eSegments

TObjArray* EdbEDADecaySearch::eSegments

eSet

EdbEDATrackSet* EdbEDADecaySearch::eSet

eSmallKink

int EdbEDADecaySearch::eSmallKink

eSmallKinkNpl

int EdbEDADecaySearch::eSmallKinkNpl

eSmallKinkRmin

float EdbEDADecaySearch::eSmallKinkRmin

eTracks

TObjArray* EdbEDADecaySearch::eTracks

eTSDau

int EdbEDADecaySearch::eTSDau

eTSDau2

int EdbEDADecaySearch::eTSDau2

eTSDauIPHist1

TH1F* EdbEDADecaySearch::eTSDauIPHist1

eTSDauIPHist2

TH1F* EdbEDADecaySearch::eTSDauIPHist2

eTSDauNseg

int EdbEDADecaySearch::eTSDauNseg

eTSDauNseg2

int EdbEDADecaySearch::eTSDauNseg2

eTSDauPH

double EdbEDADecaySearch::eTSDauPH

eTSPar

int EdbEDADecaySearch::eTSPar

eTSParIP

double EdbEDADecaySearch::eTSParIP

eTSParNseg

int EdbEDADecaySearch::eTSParNseg

eTSParPH

double EdbEDADecaySearch::eTSParPH

eTSParPlate

int EdbEDADecaySearch::eTSParPlate

eVertex

EdbVertex* EdbEDADecaySearch::eVertex

eVertexPlatePID

int EdbEDADecaySearch::eVertexPlatePID

eVertices

TObjArray* EdbEDADecaySearch::eVertices

eVtxUp

int EdbEDADecaySearch::eVtxUp

eVtxUpDZ

int EdbEDADecaySearch::eVtxUpDZ

eVtxUpIP

int EdbEDADecaySearch::eVtxUpIP

---

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

• src/libEDA/EdbEDADecaySearch.h
• src/libEDA/EdbEDADecaySearch.C