//*LastUpdate : 9-June-2000 Akiya Miyamoto
//*-- Author : Akiya Miyamoto 9-June-2000
///////////////////////////////////////////////////////////////////
//
// Do topological vertex finding.
//
//$Id: JSFZVTOP3.cxx,v 1.4 2003/05/25 10:44:51 miyamoto Exp $
//
//////////////////////////////////////////////////////////////////
// #define debug
#include "JSFZVTOP3.h"
#include "JSFSteer.h"
#include "JSFSIMDST.h"
#include "JSFVertexing.h"
#include "JSFGeoCFit.h"
#include "JSFGeneratorParticle.h"
#include <TNtuple.h>
#include "JSFLCFULL.h"
ClassImp(JSFZVTOP3Track)
ClassImp(JSFZVTOP3Vertex)
ClassImp(JSFZVTOP3)
const Int_t kMaxTrks=100;
JSFCDCTrack **JSFZVTOP3::fTrackPointers=(JSFCDCTrack**)(new UInt_t[kMaxTrks]);
extern "C" {
void zxswim_(Int_t *mode, Float_t *ok,Float_t vpos[], Int_t *id){
JSFZVTOP3::ZXSWIM(mode, ok, vpos, id);
};
Int_t zxfit_(Int_t *mode, Int_t *jj, Int_t vid[], Float_t vpos[],
Float_t vpossg[],Float_t *chisq, Float_t xvtx[], Float_t xvtxsg[],
Float_t chisqtk[],Float_t pxyz[][3], Int_t *ier){
return JSFZVTOP3::ZXFIT(mode,jj,vid,vpos,vpossg,chisq,xvtx,xvtxsg,chisqtk,pxyz,ier);
};
extern Int_t zvres3_(Int_t *ntrk, Int_t id[],Float_t *pxj, Float_t *pyj,
Float_t *pzj,Float_t *ptotj, Int_t *jetno);
extern Int_t zvtop3_(Int_t *ntrk, Int_t id[], Int_t *jetno, Int_t *errflg);
Int_t jsfzvtop3_gufld_(Float_t pos[], Float_t b0[]){
return JSFZVTOP3::GUFLD(pos, b0);
}
extern Int_t zvptm_(Int_t *mode, Float_t axi[3], Float_t eaxi[6],
Float_t axf[3], Float_t eaxf[6], Float_t p[3],
Float_t *sigmax, Float_t *p3norm, Float_t *pt,
Float_t *ptmin, Float_t *ptmax);
extern void uconfl_(Float_t *prob, Int_t *ndf, Float_t *chisq);
Float_t prob_(Float_t *chisq, Int_t *ndf){
Float_t pr;
uconfl_(&pr, ndf, chisq);
return pr;
}
};
typedef struct {
Int_t ibnk[9];
Float_t rbnk[16];
Int_t itrk[kMaxTrks][12];
Float_t rtrk[kMaxTrks][52];
} COMMON_ZXTRKS;
extern COMMON_ZXTRKS zxtrks_bank_;
typedef struct {
Int_t tidy,mc,algo, mode,nvreq;
Float_t prbv, b3min, disne, mkmin,mkmax,mlmin,mlmax,
mgmax,dkmin,dlmin,rimin,i3min,icutp,icutn,ripe,
zipe,pxja,pyja,pzja,rcut,xcut,kang,kipw,pwid,swim,
gini,gmin,cmax,pcut,momf;
} COMMON_ZVKON3;
extern COMMON_ZVKON3 zvkon3_bank_;
const Int_t kMaxTopl3Jets=20;
const Int_t kMaxTopl3Trks=40;
typedef struct {
Int_t ibnk[kMaxTopl3Jets][15];
Int_t ivrt[kMaxTopl3Jets][7][2];
Int_t itrk[kMaxTopl3Jets][kMaxTopl3Trks][4];
Float_t rbnk[kMaxTopl3Jets][22];
Float_t rvrt[kMaxTopl3Jets][7][14];
Float_t rtrk[kMaxTopl3Jets][kMaxTopl3Trks][17];
} COMMON_ZVTOPL3;
extern COMMON_ZVTOPL3 zvtopl3_;
typedef struct {
Float_t ecm, decm,beampos[3],dbeamps[6];
Float_t pol,dpol,pdtime;
Float_t xyzpri[3],dxyzpri[6],rmspri;
Int_t rtrnpri;
} COMMON_PHBM ;
extern COMMON_PHBM phbm_bank_;
TNtuple *ntrlodi=0;
//_____________________________________________________________________________
JSFZVTOP3::JSFZVTOP3()
{
fNTRK = 0 ; // number of input tracks.
fNVRT = 0 ; // number of found vertices.
fNISO = 0 ; // number of tracks left isolated.
for(Int_t i=0;i<4;i++) { fPTKS[i]=0.0; }
for(Int_t i=0;i<3;i++) { fIPX[i]=0.0; }
for(Int_t i=0;i<3;i++) { fAXI[i]=0.0; }
fGWID = 0.0 ; // Fitted ghost track gaussian width /um
fMSPTM = 0.0 ; //
fMSPT = 0.0 ;
fPTM2 = 0.0 ;
fChisq = 0.0 ; //
fProb = 0.0 ; //
fProb2 = 0.0 ; // Vfit probability of first fit.
fProbo = 0.0 ;// Probability of original vertex
fNseco = 0 ;
fNsec = 0 ; //
fDecayLength = 0.0 ; //
fDecayLength0 = 0.0 ; //
fSecondaryMass = 0.0 ; //
fSecondaryMomentum = 0.0 ; //
// if( ntrlodi == 0 ) {
// ntrlodi=new TNtuple("ntrlodi","Ntuple","trdi:lodi:disv:lodidisv");
// }
zvkon3_bank_.tidy=0 ; // set TIDY=0 to use all tracks
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.TIDY","0"),"%d",&zvkon3_bank_.tidy);
zvkon3_bank_.prbv=0.01 ; // consider V if fit prob. > PRBV
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.PRBV","0.01"),"%g",&zvkon3_bank_.prbv);
zvkon3_bank_.b3min=3.0 ; // both V tracks must be B3MIN sigma from IP
// (for Kshort and lambda only)
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.B3MIN","3.0"),"%g",&zvkon3_bank_.b3min);
zvkon3_bank_.disne=250.0; // min DIStance to NEarest other track from V
// in microns (for Kshort, lambda and gamma)
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.DISNE","250.0"),"%g",&zvkon3_bank_.disne);
zvkon3_bank_.mkmin=0.485; // lower Mass range for Kshort identification
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.MKMIN","0.485"),"%g",&zvkon3_bank_.mkmin);
zvkon3_bank_.mkmax=0.515; // upper Mass...
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.MKMAX","0.515"),"%g",&zvkon3_bank_.mkmax);
zvkon3_bank_.mlmin=1.112; // lower Mass range for Lambda identification
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.MLMIN","1.112"),"%g",&zvkon3_bank_.mlmin);
zvkon3_bank_.mlmax=1.120; // upper Mass (..both pion-proton combinations)
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.MLMAX","1.120"),"%g",&zvkon3_bank_.mlmax);
zvkon3_bank_.mgmax=0.025; // MAX Mass for e+e- conversion
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.MGMAX","0.025"),"%g",&zvkon3_bank_.mgmax);
zvkon3_bank_.dkmin=0.5; // MIN distance/cm from IP to V for Kshort
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.DKMIN","0.5"),"%g",&zvkon3_bank_.dkmin);
zvkon3_bank_.dlmin=1.0; // MIN distance/cm from IP to V for Lambda
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.DLMIN","1.0"),"%g",&zvkon3_bank_.dlmin);
zvkon3_bank_.rimin=2.2; // MIN Radius/cm for Interaction vertex
// (for gamma conversion, detector interaction)
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.RIMIN","2.2"),"%g",&zvkon3_bank_.rimin);
zvkon3_bank_.i3min=0.15; // remove if both tracks 3D Imp. para. > I3MIN/cm
// (detector interaction selection only)
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.I3MIN","0.15"),"%g",&zvkon3_bank_.i3min);
zvkon3_bank_.icutp=0.3; // remove ANY track 3D Imp. para > ICUTP/cm
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.ICUTP","0.3"),"%g",&zvkon3_bank_.icutp);
zvkon3_bank_.icutn=-0.1; // remove ANY track 3D Imp. para < ICUTN/cm
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.ICUTN","-0.1"),"%g",&zvkon3_bank_.icutn);
//
//! general ZVTOP input
//
zvkon3_bank_.ripe=7.0; // xy error on IP ellipsoid in ZVTOP, microns.
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.RIPE","7.0"),"%g",&zvkon3_bank_.ripe);
zvkon3_bank_.zipe=30.0; // z error on IP ellipsoid in ZVTOP, microns.
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.ZIPE","30.0"),"%g",&zvkon3_bank_.zipe);
zvkon3_bank_.pxja=0.0; // Jet Axis direction used to calculate angle
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.PXJA","0.0"),"%g",&zvkon3_bank_.pxja);
zvkon3_bank_.pyja=0.0; // A for KANG weight for ALGO=1 and for initial
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.PYJA","0.0"),"%g",&zvkon3_bank_.pyja);
zvkon3_bank_.pzja=0.0; // ghost axis in ALGO=2. By default (all P*JA = 0)
// the sum of the input track momenta is used.
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.PZJA","0.0"),"%g",&zvkon3_bank_.pzja);
zvkon3_bank_.mc=0; // MC=1 : include Monte Carlo information in
// the output.
zvkon3_bank_.mc=gJSF->Env()->GetValue("JSFZVTOP3.MC",0);
//! select ZVTOP algorithm
zvkon3_bank_.algo=1; // 1 = ZVRES original pure vertex RESolubility
// 2 = ZVKIN KINematic info and ghost track
zvkon3_bank_.algo=gJSF->Env()->GetValue("JSFZVTOP3.ALGO",1);
//! tunable inputs for ZVRES
zvkon3_bank_.rcut=0.6; // resolubility cut, can vary as 0.0 < VCUT < 1.0,
// larger values produce more vertices.
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.RCUT","0.6"),"%g",&zvkon3_bank_.rcut);
zvkon3_bank_.xcut=10.0; // cut on the chi**2 contribution of a track to
// a vertex.
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.XCUT","10.0"),"%g",&zvkon3_bank_.xcut);
zvkon3_bank_.kang=5.0; // Weights vertex finding probability by the
// factor exp(-KANG*A**2) where A is the angle
// in radians between the line joining the IP
// to the 3D spatial co-ordinate and the jet
// axis. (+ve values ~10 increase vertex
// finding efficiency in core of jet).
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.KANG","5.0"),"%g",&zvkon3_bank_.kang);
zvkon3_bank_.kipw=1.0; // i.p track weight, weights the significance
// of the IP, large values will tend to absorb
// nearby vertices/fakes into the primary vertex.
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.KIPW","1.0"),"%g",&zvkon3_bank_.kipw);
zvkon3_bank_.mode=0; // 0 -usual ZVTOP/ZVRES
// 1 -lepton or high impact track mode - Note 2
// 2 -like 1, but remove track for vertex finding
// 3 -...and remove space on IP side of lepton
zvkon3_bank_.mode=gJSF->Env()->GetValue("JSFZVTOP3.MODE",0);
zvkon3_bank_.pwid=100.0; // Gaussian WIDth of Plane for MODE>=1 in microns
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.PWID","100.0"),"%g",&zvkon3_bank_.pwid);
zvkon3_bank_.swim=0.0; // 1.0 - use Aaron's track swim errors in ZVRES
// SWIM>0.0 multiply swam track by SWIM
// SWIM<0.0 multiply non-swam track by |SWIM|
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.SWIM","0.0"),"%g",&zvkon3_bank_.swim);
//
//! tunable inputs for ZVKIN
//
zvkon3_bank_.gini=25.0; // Initial ghost track width/um
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.GINI","25.0"),"%g",&zvkon3_bank_.gini);
zvkon3_bank_.gmin=25.0; // Minimum ghost track width/um
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.GMIN","25.0"),"%g",&zvkon3_bank_.gmin);
zvkon3_bank_.cmax=1.0; // Max chi**2 of +ve LODI track to ghost
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.CMAX","1.0"),"%g",&zvkon3_bank_.cmax);
zvkon3_bank_.pcut=0.01; // Min prob. to cluster track vertices
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.PCUT","0.01"),"%g",&zvkon3_bank_.pcut);
zvkon3_bank_.nvreq=0; // Number of reconstructed Vertices REQuested
// NVREQ = 0, no request, PCUT is used
// NVREQ > 0, vertex clustering stops when
// vertex candidates reduced to NVREQ.
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.NVREQ","0"),"%d",&zvkon3_bank_.nvreq);
zvkon3_bank_.momf=0.0; // MOMentum Factor IP/BD discrimination
// MOMF=1.0 is fairly large momentum weight
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.MOMF","0.0"),"%g",&zvkon3_bank_.momf);
// Note 1
// TIDY > 0 initially loops over all track pairs and find the
// furthest V vertex from IP with fit probability > PRBV.
// If the 2 satisfies the above cuts for Kshort, lambda
// or gamma (and is 0 charge vertex), or interaction
// with detector material, the two tracks are flagged
// and the next furthest 2-prong is then considered
// iteratively until no more tracks are flagged. Flagged
// tracks are not used further by ZVTOP, but they are
// output with their flag in ZVTOPL.
//
// TIDY = 2 if N >= 2, (where N is the number of secondary tracks
// (ALGO=2) found by ZVKIN) they are fit to a common vertex distance
// D along vertex axis. Tracks that initially failed ICUTP
// but pass close to the vertex (in TRDI and LODI - D) are
// 'unflagged' and the whole algorithm is called again.
// Similarly, if N < 2, all tracks initially failing ICUTP
// are put back into the jet and the algorithm rerun.
// (ZVTOPL3.GWID is set -ve to flag such a second iteration).
//
// Note 2
// MODE > 0 the first track passed to ZVTOP, id(1), together
// (ALGO=1) with the IP location, is used to define a plane.
// Locations not in the plane are suppressed as
// a function of distance from the plane by weighing
// V(r) with Gaussian factor tunable via PWID, smaller
// values to PWID confine vertex search closer to plane.
// - For track id(1) a reasonable 3D impact distance
// to the IP, > 2 or 3 sigma, is recommended.
phbm_bank_.beampos[0]=0;
phbm_bank_.beampos[1]=0;
phbm_bank_.beampos[2]=0;
for(Int_t i=0;i<6;i++){ phbm_bank_.dbeamps[i]=0.0; }
phbm_bank_.dbeamps[0]=1.0e-4;
phbm_bank_.dbeamps[2]=1.0e-4;
phbm_bank_.dbeamps[3]=1.0e-4;
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.BEAMSIGX","0.0001"),"%g",&phbm_bank_.dbeamps[0]);
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.BEAMSIGY","0.0001"),"%g",&phbm_bank_.dbeamps[1]);
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.BEAMSIGZ","0.0001"),"%g",&phbm_bank_.dbeamps[2]);
}
//_____________________________________________________________________________
JSFZVTOP3::~JSFZVTOP3()
{
Delete();
}
// ---------------------------------------------------------------
void JSFZVTOP3::ZXSWIM(Int_t *mode, Float_t *ok, Float_t vpos[], Int_t *id)
{
//
printf(" ZXSWIM is called. .. mode=%d",*mode);
printf(" id=%d",*id);
printf(" vpos=%g %g %g",vpos[0],vpos[1],vpos[2]);
*ok=0.5;
return ;
}
// ---------------------------------------------------------------
Int_t JSFZVTOP3::ZXFIT(Int_t *mode, Int_t *jj, Int_t vid[], Float_t vpos[],
Float_t vpossg[], Float_t *chisq, Float_t xvtx[], Float_t xvtxsg[],
Float_t chisqtk[],Float_t pxyz[][3], Int_t *ier)
{
// Does vertex fitting...
//(Inputs)
// mode ;
// jj ; Number of tracks
// vid[]; An array of track id whose vertex is fitted.
// vpos[3] : Initial vertex position
// vpossg[3] : error of vertex position.
//(Output)
// chisq : Chisq of vertex fit
// xvtx(3) : Fitted vertex (x,y,z)
// xvtxsq(6) : Error matrix of vertex position
// chisqtk() : Chi-square of each track.
// pxyz[][3] : three momentum at the vertex
Int_t ntrk=*jj;
if( *jj < 2 ) {
printf(" JSFZVTOP3::ZXFIT was called, but the number of tracks is less than 2, *jj=%dn",*jj);
printf(" mode=%dn",*mode);
printf(" vpos=%g %g %gn",vpos[0],vpos[1],vpos[2]);
return 4;
}
Float_t b0[3];
JSFZVTOP3::GUFLD(vpos,b0);
Float_t bf=b0[2];
TClonesArray fitbuf("JSFHelicalTrack",ntrk);
JSFVertexing vtx(2);
for(Int_t i=0;i<ntrk;i++){
if( vid[i] < 1 ) {
printf(" JSFZVTOP3::ZXFIT .. Error: invalid index to track pointersn");
printf(" Program may be broken somewhere.n");
return 4;
}
JSFCDCTrack *t=JSFZVTOP3::fTrackPointers[vid[i]-1];
if( !t ) {
printf(" JSFZVTOP3::ZXFIT .. Error: Pointer to CDC track is not available.n");
return 4;
}
Int_t ndf=t->GetNDF();
Double_t errmtr[15];
Float_t helix[5], pivot[3];
t->GetHelix(helix);
t->GetPivot(pivot);
t->GetError(errmtr);
new( fitbuf[i] ) JSFHelicalTrack(helix, pivot, ndf, 1.0, 0.0, errmtr);
JSFHelicalTrack *ht=(JSFHelicalTrack*)fitbuf.At(i);
ht->SetBfield(bf);
if( ntrk == 2 ) { vtx.SetTrack(i,*ht); }
}
TVector3 v(vpos[0], vpos[1], vpos[2]);
JSFGeoCFit fit(ntrk, &fitbuf, v);
Bool_t err=fit.Fit();
// printf(" Fit err =%d ntry=%dn",err,fit.GetNtry());
#if __ROOT_FULLVERSION__ < 30000
JSFDMatrix ans(fit.GetParameter());
JSFDMatrix errmatrix(fit.GetErrorMatrix());
#else
TMatrixD ans(fit.GetParameter());
TMatrixD errmatrix(fit.GetErrorMatrix());
#endif
*chisq=fit.GetChisq();
xvtx[0]=ans(0,0);
xvtx[1]=ans(1,0);
xvtx[2]=ans(2,0);
fit.GetTrackChisq(chisqtk);
Int_t ip=0;
for(Int_t i=0;i<3;i++){ for(Int_t j=0;j<i+1;j++){
xvtxsg[ip++]=errmatrix(i,j);
}}
for(Int_t i=0;i<ntrk;i++){
Int_t ip=3*(i+1);
Float_t pt=1.0/TMath::Abs(ans(ip+1,0));
Float_t phi0=ans(ip,0);
Float_t tanl=ans(ip+2,0);
pxyz[i][0]=-pt*TMath::Cos(phi0);
pxyz[i][1]= pt*TMath::Sin(phi0);
pxyz[i][2]= pt*tanl;
}
fitbuf.Delete();
if( err ) return 0;
return 4;
}
// ---------------------------------------------------------------
Int_t JSFZVTOP3::GUFLD(Float_t pos[3], Float_t b0[3])
{
// Return magnetic field.
static JSFQuickSimParam *fpar=0;
if( fpar == 0 ) {
JSFSIMDST *sds=(JSFSIMDST*)gJSF->FindModule("JSFSIMDST");
fpar=sds->Param();
}
b0[0]=0;
b0[1]=0;
b0[2]=fpar->GetBField();
return 0;
}
// ---------------------------------------------------------------
void JSFZVTOP3::MakeVertices(TObjArray *tracks)
{
// Do vertexing using ZVTOP Algorithm
//(Input)
// tracks : Pointers to JSFCDCTrack class which belongs to this track
//(Output)
// Results are saved as data member of the class.
// Loop over jets: select charged track parameter and
// move its pivot to IP.
Double_t err[15];
Float_t trk[30];
// Clear array for vertex finding
Int_t nt=0;
for(Int_t i=0;i<9;i++){ zxtrks_bank_.ibnk[i]=0; }
for(Int_t i=0;i<16;i++){ zxtrks_bank_.rbnk[i]=0; }
for(Int_t ic=0;ic<kMaxTrks;ic++){
for(Int_t i=0;i<12;i++){ zxtrks_bank_.itrk[ic][i]=0; }
for(Int_t i=0;i<52;i++){ zxtrks_bank_.rtrk[ic][i]=0; }
}
JSFSIMDST *sds=(JSFSIMDST*)gJSF->FindModule("JSFSIMDST");
JSFQuickSimParam *fpar=sds->Param();
// Loop over particles in jet to prepare zvtrks common.
Int_t trkid[kMaxTrks];
Int_t numtrk=0;
Float_t pntip[3]={0.0, 0.0, 0.0};
Float_t bf=fpar->GetBField();
TIter next(tracks);
JSFCDCTrack *ct;
while( (ct=(JSFCDCTrack*)next())) {
ct->GetTrackParameter(trk);
Float_t rclosest=TMath::Sqrt(trk[4]*trk[4]+trk[5]*trk[5]);
// if( rclosest > 1.0 ) continue;
if( rclosest > 2.0 ) {
ct->MovePivot(pntip,bf);
}
else {
Float_t pv[3];
ct->GetPivot(pv);
// If current pivot position is more than 1 cm away from IP,
// move it to IP.
if( pv[0]*pv[0]+pv[1]*pv[1]+pv[2]*pv[2] > 1.0 ) {
ct->MovePivotToIP(fpar);
}
}
ct->GetErrorMatrix(err);
ct->GetTrackParameter(trk);
zxtrks_bank_.rtrk[nt][0]=ct->GetCharge();
zxtrks_bank_.rtrk[nt][1]=1;
trkid[numtrk]=nt+1;
if( numtrk > kMaxTrks ) {
printf(" JSFZVTOP3::MakeVertices .. Number of CDC tracks(%d)",numtrk);
printf(" exceeds buffer size. Track#%d is neglected.n",nt);
continue;
}
numtrk++;
JSFZVTOP3::fTrackPointers[nt]=ct;
Float_t hlx[5];
Float_t pivot[3];
ct->GetHelix(hlx);
ct->GetPivot(pivot);
zxtrks_bank_.rtrk[nt][23+1]=hlx[1]+TMath::Pi()/2;
zxtrks_bank_.rtrk[nt][23+2]=TMath::Abs(hlx[2]);
zxtrks_bank_.rtrk[nt][23+3]=hlx[4];
#ifdef debug
Float_t cosct=ct->GetPz()/TMath::Sqrt(ct->GetPx()*ct->GetPx()
+ ct->GetPy()*ct->GetPy()+ct->GetPz()*ct->GetPz());
printf(" GenID=%d",ct->GetGenID());
printf(" nt=%d px,py,pz=%g,%g,%g costh=%g\n",
numtrk,ct->GetPx(),ct->GetPy(),ct->GetPz(),cosct);
#endif
Double_t xip=pivot[0]+hlx[0]*TMath::Cos(hlx[1]);
Double_t yip=pivot[1]+hlx[0]*TMath::Sin(hlx[1]);
Double_t zip=pivot[2]+hlx[3];
zxtrks_bank_.rtrk[nt][23+4]=xip;
zxtrks_bank_.rtrk[nt][23+5]=yip;
zxtrks_bank_.rtrk[nt][23+6]=zip;
#ifdef debug
printf(" track %d ip %g %g %g ",nt,xip,yip,zip);
#endif
zxtrks_bank_.rtrk[nt][23]=TMath::Sqrt(xip*xip+yip*yip);
zxtrks_bank_.rtrk[nt][22]=zxtrks_bank_.rtrk[nt][23]/
TMath::Sqrt(err[0]);
zxtrks_bank_.rtrk[nt][48]=TMath::Abs(zip)*TMath::Sqrt(xip*xip+yip*yip)/
TMath::Sqrt(xip*xip+yip*yip+zip*zip);
zxtrks_bank_.rtrk[nt][49]=
TMath::Sqrt( (xip*xip+yip*yip)/err[0]+zip*zip/err[9]);
#ifdef debug
printf(" 2d drp=%g drp/s=%g ",zxtrks_bank_.rtrk[nt][23],zxtrks_bank_.rtrk[nt][22]);
printf(" 3d drp=%g drp/s=%g ",zxtrks_bank_.rtrk[nt][48],zxtrks_bank_.rtrk[nt][49]);
printf("\n");
#endif
zxtrks_bank_.rtrk[nt][29+10]=err[2];
zxtrks_bank_.rtrk[nt][29+15]=err[9];
nt++;
zxtrks_bank_.ibnk[8]=nt;
}
// Call ZVRES3
Int_t errflg;
Int_t jetno=1;
#ifdef debug
printf(" \n*********************************\n");
printf(" start ZVTOP3 .. jetno=%d no of trackis %d\n",jetno,numtrk);
#endif
Int_t iret=zvtop3_(&numtrk, &trkid[0], &jetno, &errflg);
if( iret != 1 ) {
printf(" error return from ZVTOP3 .. iret=%d\n",iret);
}
// Save results in Class members
fNTRK=zvtopl3_.ibnk[0][0];
fNVRT=zvtopl3_.ibnk[0][1];
fNISO=zvtopl3_.ibnk[0][2];
#ifdef debug
printf(" fNTRK=%d fNVRT=%d fNISO=%d\n",fNTRK, fNVRT, fNISO);
#endif
for(Int_t i=0;i<4;i++){ fPTKS[i]=zvtopl3_.ibnk[0][i+3]; }
for(Int_t i=0;i<3;i++){ fIPX[i] =zvtopl3_.ibnk[0][i+7]; }
for(Int_t i=0;i<3;i++){ fAXI[i] =zvtopl3_.ibnk[0][i+10]; }
fGWID=zvtopl3_.ibnk[0][13];
// Save vertex information
for(Int_t iv=0;iv<fNVRT;iv++){
Add(new JSFZVTOP3Vertex(iv));
}
CalculateMSPTM();
return ;
}
//_________________________________________________________
void JSFZVTOP3::DebugPrint()
{
// Print a vertexing result.
Int_t ij=0;
printf(" jet %dn",ij);
printf(" #Tracks %d # Vertex %d #iso=%dn",
zvtopl3_.ibnk[ij][0],zvtopl3_.ibnk[ij][1],zvtopl3_.ibnk[ij][2]);
printf(" RBNK_ZVTOPL3 informationn");
printf(" PTKS(XYZ magnitude of summed track momentum):%g %g %g %gn",
zvtopl3_.rbnk[ij][0], zvtopl3_.rbnk[ij][1],
zvtopl3_.rbnk[ij][2], zvtopl3_.rbnk[ij][3]);
printf(" IPX(XYZ location of IP from PHBM:%g %g %gn",
zvtopl3_.rbnk[ij][4], zvtopl3_.rbnk[ij][5],zvtopl3_.rbnk[ij][6]);
printf(" MSPTM=%g",GetMSPTM());
printf(" MSEC=%g ",GetSecondaryMass());
printf(" DecayLength=%g",GetDecayLength());
printf("n");
Int_t nvrt=zvtopl3_.ibnk[ij][1];
printf(" nvrt=%dn",nvrt);
for(Int_t iv=0;iv<nvrt;iv++){
printf(" Vertex Number %dn",iv);
printf(" Number of tracks in this vertex : %dn",
zvtopl3_.ivrt[ij][iv][0]);
printf(" vertex charge : %dn" ,zvtopl3_.ivrt[ij][iv][1]);
printf(" Vertex significance func. value:%gn",zvtopl3_.rvrt[ij][iv][0]);
printf(" Vertex position:%g %g %gn",zvtopl3_.rvrt[ij][iv][1],
zvtopl3_.rvrt[ij][iv][2],zvtopl3_.rvrt[ij][iv][3]);
printf(" Chisq of vertex fit:%gn",zvtopl3_.rvrt[ij][iv][10]);
printf(" Probability of vertex fit:%gn",zvtopl3_.rvrt[ij][iv][11]);
printf(" vertex mass:%gn",zvtopl3_.rvrt[ij][iv][12]);
}
printf(" Track id and Vertex number in ths jetn");
for(Int_t k=0;k< zvtopl3_.ibnk[ij][0];k++){
printf(" (Id,Vno)=(%d,%d)",zvtopl3_.itrk[ij][k][0],zvtopl3_.itrk[ij][k][1]);
printf(" chisq=%g",zvtopl3_.rtrk[ij][k][0]);
printf(" trdi,lodi=%g,%g",zvtopl3_.rtrk[ij][k][15],zvtopl3_.rtrk[ij][k][16]);
printf("n");
}
}
//_________________________________________________________
void JSFZVTOP3::CalculateMSPTM()
{
Int_t jetno=0;
Int_t idsec[50];
fMSPTM=0.0;
fProbo=0.0;
fNseco=0;
fProb2=0.0;
fDecayLength=0.0;
fDecayLength0=0.0;
fNsec =0;
fChisq=0.0;
fSecondaryMomentum=0.0;
fSecondaryMass=0.0;
fMSPT = 0.0;
fPTM2 = 0.0;
Int_t nvrt=fNVRT;
#ifdef debug
printf("************************************\n");
printf(" Start MSPTM... nvrt=%d\n",nvrt);
#endif
if ( nvrt < 2 ) return ;
// Vertex axis
Int_t ivrt=1;
Float_t vax = zvtopl3_.rvrt[jetno][ivrt][1];
Float_t vay = zvtopl3_.rvrt[jetno][ivrt][2];
Float_t vaz = zvtopl3_.rvrt[jetno][ivrt][3];
Float_t chisq=zvtopl3_.rvrt[jetno][ivrt][10];
Int_t nseco=zvtopl3_.ivrt[jetno][ivrt][0];
Float_t probo;
Int_t ndfo=5*nseco -3*(nseco+1);
uconfl_(&probo, &ndfo, &chisq);
// take furthest rec. vertex as seed to add isolated tracks.
Float_t disv = TMath::Sqrt(vax*vax + vay*vay + vaz*vaz);
#ifdef debug
printf(" Vertex position is %g %g %g disv=%g \n",vax,vay,vaz,disv);
#endif
#if debug
Double_t verr[3][3];
verr[0][0]=zvtopl3_.rvrt[jetno][ivrt][4];
verr[1][0]=zvtopl3_.rvrt[jetno][ivrt][5];
verr[1][1]=zvtopl3_.rvrt[jetno][ivrt][6];
verr[2][0]=zvtopl3_.rvrt[jetno][ivrt][7];
verr[2][1]=zvtopl3_.rvrt[jetno][ivrt][8];
verr[2][2]=zvtopl3_.rvrt[jetno][ivrt][9];
verr[0][1]=verr[1][0];
verr[0][2]=verr[2][0];
verr[1][2]=verr[2][1];
Double_t xyz[3]={vax,vay,vaz};
Double_t errw=0.0;
Int_t i,j;
for(i=0;i<3;i++){ for(j=0;j<3;j++) {
errw+= xyz[i]*verr[i][j]*xyz[j];
}}
Float_t disverr=disv*disv/TMath::Sqrt(errw);
printf(" disv=%g sqrt(errw)=%g ",disv,TMath::Sqrt(errw)/disv);
printf(" disverr = %g \n",disverr);
#endif
Int_t nsec=0;
Float_t ptot=0.0;
Float_t pxt =0.0;
Float_t pyt =0.0;
Float_t pzt =0.0;
Float_t etot=0.0;
#if debug
Int_t ievt=gJSF->GetEventNumber();
#endif
TVector3 vtx1(vax, vay, vaz);
TVector3 vtx2;
Float_t vdis=100.0;
if( nvrt > 2 ) {
vtx2=TVector3( zvtopl3_.rvrt[jetno][2][1],
zvtopl3_.rvrt[jetno][2][2], zvtopl3_.rvrt[jetno][2][3]);
vdis=(vtx1-vtx2).Mag();
}
#ifdef debug
printf("\n nvrt=%d disv=%g\n",nvrt,disv);
printf(" dl1=%g ",vtx1.Mag());
if( nvrt > 2 ) { printf(" dl2=%g ",vtx2.Mag()) ; }
printf("\n");
#endif
Float_t lodimin;
Float_t trdimax;
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.MSPTM.LODIMIN","0.3"),"%g",&lodimin);
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.MSPTM.TRDIMAX","0.1"),"%g",&trdimax);
Int_t mtvers=gJSF->Env()->GetValue("JSFZVTOP3.MSPTM.MTVERS",0);
for(Int_t it=0;it<GetNTRK();it++){
Float_t trdi = zvtopl3_.rtrk[jetno][it][15];
Float_t lodi = zvtopl3_.rtrk[jetno][it][16];
// trdi : distance in cm from the track to the IP ( no secondary)
// or closest approach to the track from the IP pasing through ...
// lodi : distrance from IP along 2.5cm line to TRDI
#ifdef debug
printf(" it =%d trdi=%g lodi=%g lodi/disv=%g \n",it,trdi,lodi,lodi/disv);
#endif
/*
Float_t array[3];
array[0]=trdi;
array[1]=lodi;
array[2]=disv;
array[3]=lodi/disv;
ntrlodi->Fill(array);
*/
if( (mtvers == 0 && (zvtopl3_.itrk[jetno][it][1] == 2 ||
( zvtopl3_.itrk[jetno][it][1] > 2 &&
lodi/disv > lodimin && trdi < trdimax ))) ||
(mtvers == 3 && ( zvtopl3_.itrk[jetno][it][1] == 2 ||
( zvtopl3_.itrk[jetno][it][1] < 2 &&
lodi/disv > lodimin && trdi < trdimax ))) ||
(mtvers == 4 && ( zvtopl3_.itrk[jetno][it][1] == nvrt ||
( zvtopl3_.itrk[jetno][it][1] < nvrt &&
lodi/disv > lodimin && trdi < trdimax ))) ) {
Int_t ids=zvtopl3_.itrk[jetno][it][0]-1;
Float_t cosl=TMath::Cos(TMath::ATan(zxtrks_bank_.rtrk[ids][23+3]));
Float_t sinl=TMath::Sin(TMath::ATan(zxtrks_bank_.rtrk[ids][23+3]));
// ptot=1.0/(zxtrks_bank_.rtrk[ids][23+2]*cosl);
Float_t tanl=zxtrks_bank_.rtrk[ids][23+3];
ptot=TMath::Sqrt(1.0+tanl*tanl)/zxtrks_bank_.rtrk[ids][23+2];
Float_t px=ptot*cosl*TMath::Cos(zxtrks_bank_.rtrk[ids][23+1]);
Float_t py=ptot*cosl*TMath::Sin(zxtrks_bank_.rtrk[ids][23+1]);
Float_t pz=ptot*sinl;
#ifdef debug
printf(" px,py,pz=%g %g %g\n",px,py,pz);
#endif
pxt += px ;
pyt += py ;
pzt += pz ;
etot += TMath::Sqrt(ptot*ptot + 0.01948);
idsec[nsec++]=zvtopl3_.itrk[jetno][it][0];
} // Track is in secondary
} // Loop over quality tracks in hemi
if( nsec < 2 ) {
printf(" in JSFZVTOP3:: nsec=%d\n",nsec);
printf(" nvrt=%d\n",nvrt);
return ;
}
fProbo=probo;
fNseco=nseco;
fProb2 =zvtopl3_.rvrt[jetno][ivrt][11];
Float_t vpos[3]={0.0, 0.0, 0.0};
Float_t vpossg[3]={10.0, 10.0, 10.0};
Int_t mode=1;
Float_t xvtx[3], xvtxsg[6];
Float_t chisqtk[50], pxyz[50][3];
Int_t ierr;
ZXFIT(&mode, &nsec, idsec, vpos , vpossg,
&chisq, xvtx, xvtxsg, chisqtk, pxyz, &ierr);
Float_t dsec=TMath::Sqrt(
TMath::Power((xvtx[0]-zvtopl3_.rbnk[jetno][4]),2) +
TMath::Power((xvtx[1]-zvtopl3_.rbnk[jetno][5]),2) +
TMath::Power((xvtx[2]-zvtopl3_.rbnk[jetno][6]),2));
Int_t ndf=5*nsec -3*(nsec+1);
uconfl_(&fProb, &ndf, &chisq);
fDecayLength=dsec;
fDecayLength0=disv;
fNsec = nsec;
fChisq= chisq;
Float_t msec = TMath::Sqrt(etot*etot-pxt*pxt - pyt*pyt - pzt*pzt);
Float_t moms = TMath::Sqrt(pxt*pxt+pyt*pyt+pzt*pzt);
Float_t plse = (pxt*vax + pyt*vay + pzt*vaz)/disv;
Float_t pt3d = TMath::Sqrt( TMath::Abs((moms-plse)*(moms+plse)));
Float_t mspt = TMath::Sqrt(msec*msec + pt3d*pt3d) + pt3d;
fSecondaryMomentum=moms;
fSecondaryMass=msec;
Float_t axi[3];
axi[0]=zvtopl3_.rbnk[jetno][4];
axi[1]=zvtopl3_.rbnk[jetno][5];
axi[2]=zvtopl3_.rbnk[jetno][6];
Float_t eaxi[6];
eaxi[0]=zvkon3_bank_.ripe*zvkon3_bank_.ripe*1.0E-8;
eaxi[2]=zvkon3_bank_.ripe*zvkon3_bank_.ripe*1.0E-8;
eaxi[5]=zvkon3_bank_.zipe*zvkon3_bank_.zipe*1.0E-8;
// eaxi[0]=0.0007*0.0007 ; eaxi[2]=0.0007*0.0007; eaxi[5]=0.007*0.007;
// eaxi[0]=0.0007*0.0007 ; eaxi[2]=0.0007*0.0007; eaxi[5]=0.003*0.003;
eaxi[0]=0.0001*0.0001 ; eaxi[2]=0.0001*0.0001; eaxi[5]=0.0004*0.0004;
eaxi[1]=0.0; eaxi[3]=0; eaxi[4]=0.0;
Float_t p[3];
p[0]=pxt ; p[1]=pyt ; p[2]=pzt;
mode=3;
Float_t sigmax=1.0;
sscanf(gJSF->Env()->GetValue("JSFZVTOP3.MSPTM.SIGMAX","1.0"),"%g",&sigmax);
if( gJSF->Env()->GetValue("JSFZVTOP3.MSPT.VTXSRC",0) != 0 ) {
for(Int_t i=0;i<3;i++){ xvtx[i]=zvtopl3_.rvrt[jetno][1][i+1]; }
for(Int_t i=0;i<6;i++){ xvtxsg[i]=zvtopl3_.rvrt[jetno][1][i+4]; }
}
Float_t p3norm, ptm0, ptm2, ptmax;
zvptm_(&mode, axi, eaxi, xvtx, xvtxsg, p, &sigmax, &p3norm,
&ptm0, &ptm2, &ptmax);
Float_t msptm = TMath::Sqrt(msec*msec + ptm2*ptm2) + TMath::Abs(ptm2);
#ifdef debug
printf(" refitted vertex position is %g %g %g\n",xvtx[0],xvtx[1],xvtx[2]);
#endif
#ifdef debug
if( msptm > 2.0 && dsec > 0.03 && disv > 0.03 ) {
// if( msptm > 2.0 && disverr > 3.0 ) {
printf(" prob=%g ndf=%d chisq=%gn",fProb,ndf,chisq);
printf(" msec=%g dsec=%g ",msec,dsec);
// printf(" disverr=%g",disverr);
printf(" msptm=%gn",msptm);
}
#endif
fMSPTM=msptm;
fMSPT = mspt;
fPTM2 = ptm2;
return ;
}
//_____________________________________________________________________
JSFZVTOP3Track::JSFZVTOP3Track(Int_t id, Int_t jetno)
{
//(Function)
// Copy contents of ZVTOP3 common to data member
//(Input)
// id : Track number ( 0,1,2,....)
// jetno: jet number (0,1,2,....)
fZID=zvtopl3_.itrk[jetno][id][0];
fVNO=zvtopl3_.itrk[jetno][id][1];
fXIST=zvtopl3_.rtrk[jetno][id][0];
for(Int_t i=0;i<7;i++){
fXISG[i]=zvtopl3_.rtrk[jetno][id][i+1];
fPRBG[i]=zvtopl3_.rtrk[jetno][id][i+8];
}
fTRDI = zvtopl3_.rtrk[jetno][id][15];
fLODI = zvtopl3_.rtrk[jetno][id][16];
}
//_____________________________________________________________________
JSFZVTOP3Vertex::~JSFZVTOP3Vertex()
{
fTracks->Delete();
}
//_____________________________________________________________________
JSFZVTOP3Vertex::JSFZVTOP3Vertex(Int_t vno, Int_t jetno)
{
//(Function)
// Copy contents of ZVTOP3 common to data member
//(Input)
// vno : vertex number (0,1,2,....)
// jetno: jet number (0,1,2,....)
fNcharged=zvtopl3_.ivrt[jetno][vno][0];
fCharge=zvtopl3_.ivrt[jetno][vno][1];
fVSIG =zvtopl3_.rvrt[jetno][vno][0];
fPos[0]=zvtopl3_.rvrt[jetno][vno][1];
fPos[1]=zvtopl3_.rvrt[jetno][vno][2];
fPos[2]=zvtopl3_.rvrt[jetno][vno][3];
for(Int_t i=0;i<6;i++){
fDPOS[i]=zvtopl3_.rvrt[jetno][vno][4+i];
}
fXISQ =zvtopl3_.rvrt[jetno][vno][10];
fPRBV =zvtopl3_.rvrt[jetno][vno][11];
fMass =zvtopl3_.rvrt[jetno][vno][12];
fALPHA =zvtopl3_.rvrt[jetno][vno][13];
if( zvkon3_bank_.algo == 1 ) {
Int_t ndf=(2*fNcharged-3);
fPRBV=prob_(&fXISQ,&ndf);
zvtopl3_.rvrt[jetno][vno][11]=fPRBV;
}
fTracks = new TObjArray();
for(Int_t i=0;i<zvtopl3_.ibnk[jetno][0];i++){
if( zvtopl3_.itrk[jetno][i][1] == vno+1 ) {
fTracks->Add(new JSFZVTOP3Track(i,jetno));
}
}
}
#undef debug
ROOT page - Home page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.