// ----------------------------------------------------------------------------
// $Id: LCDVTopl.cxx,v 1.3 2001/06/19 03:46:40 toshi Exp $
// ----------------------------------------------------------------------------
//
// $Log: LCDVTopl.cxx,v $
// Revision 1.3 2001/06/19 03:46:40 toshi
// Small bug fixes.
//
// Revision 1.2 2001/05/11 21:34:51 toshi
// Minar changes. MS VC++ 6 picked up some multi-difinition error of
// Int_t i .
//
// Revision 1.1 2001/05/10 19:32:09 toshi
// C++ file name convention has been changed from .C to .cxx .
//
//
////////////////////////////////////////////////////////////
//
// LCDVTopl
//
// $Log :$
//
#include "LCDVTopl.h"
#include <stdio.h>
#include <iostream.h>
ClassImp(LCDVTopl)
//__________________________________________________________________________
//Constructor and Deconstructor
LCDVTopl::LCDVTopl() {
Init();
CalcIPErrorMatrix();
}
//__________________________________________________________________________
//__________________________________________________________________________
LCDVTopl::LCDVTopl(TObjArray* track_list,
Int_t ntrk,
Int_t* track_indecies,
const TVector3& ip,
const TVector3& bfld) {
Init();
SetUpLCDVTopl(track_list,ntrk,track_indecies,ip,bfld);
}
//__________________________________________________________________________
//__________________________________________________________________________
LCDVTopl::LCDVTopl(TObjArray* track_list,
Int_t ntrk,
Int_t* track_indecies,
const TVector3& ip,
const TVector3& bfld,
const TVector3& pjet) {
Init();
SetUpLCDVTopl(track_list,ntrk,track_indecies,ip,bfld,pjet);
}
//__________________________________________________________________________
//__________________________________________________________________________
LCDVTopl::LCDVTopl(TObjArray* track_list,
const TVector3& ip,
const TVector3& bfld,
const TVector3& pjet) {
Init();
SetUpLCDVTopl(track_list,ip,bfld,pjet);
}
//__________________________________________________________________________
//__________________________________________________________________________
LCDVTopl::~LCDVTopl(){
// fclose(fp);
if (vtrack_list) {
vtrack_list->Delete();
delete vtrack_list;
vtrack_list=0;
}
if (vertex_list) {
vertex_list->Delete();
delete vertex_list;
vertex_list=0;
}
if (vrttrack_list) {
delete vrttrack_list;
vrttrack_list=0;
}
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::SetUpLCDVTopl(TObjArray* track_list,
Int_t ntrk,
Int_t* track_indecies,
const TVector3& ip,
const TVector3& bfld) {
Clean();
ipv =ip;
bfield=bfld;
SetLCDVToplTRKs(track_list,ntrk,track_indecies);
//This loop to
// make sure position and momentum in trk at 2d POCA
// and to make LCDVToplTRK.
ptks.SetXYZ(0.0,0.0,0.0);
Double_t l =0.0;
LCDVToplTRK* trk;
for (Int_t i=0; i < ntrk ; i++) {
trk= GetTrackAt(i);
//Accumurate track momentum.
ptks += trk->GetMomentumVector(l);
} //loop over zxtrks
pja=ptks;
CalcIPErrorMatrix();
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::SetUpLCDVTopl(TObjArray* track_list,
Int_t ntrk,
Int_t* track_indecies,
const TVector3& ip,
const TVector3& bfld,
const TVector3& pjet) {
Clean();
ipv =ip;
bfield=bfld;
SetLCDVToplTRKs(track_list,ntrk,track_indecies);
//This loop to
// make sure position and momentum in trk at 2d POCA
// and to make LCDVToplTRK.
ptks.SetXYZ(0.0,0.0,0.0);
Double_t l =0.0;
LCDVToplTRK* trk;
for( Int_t i=0; i < ntrk ; i++ ) {
trk= GetTrackAt(i);
//Accumurate track momentum.
ptks += trk->GetMomentumVector(l);
} //loop over zxtrks
pja=pjet;
CalcIPErrorMatrix();
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::SetUpLCDVTopl(TObjArray* track_list,
const TVector3& ip,
const TVector3& bfld,
const TVector3& pjet) {
Clean();
ipv =ip;
bfield=bfld;
LCDVToplTRK* trk;
Int_t i;
for(i=0 ; i < track_list->GetEntries() ; i++) {
trk= new((*vtrack_list)[i]) LCDVToplTRK(track_list,i,ipv);
} //loop over zxtrks
//This loop to
// make sure position and momentum in trk at 2d POCA
// and to make LCDVToplTRK.
ptks.SetXYZ(0.0,0.0,0.0);
Double_t l =0.0;
for(i=0; i < track_list->GetEntries() ; i++ ) {
trk= GetTrackAt(i);
//Accumurate track momentum.
ptks += trk->GetMomentumVector(l);
} //loop over zxtrks
pja=pjet;
CalcIPErrorMatrix();
}
//__________________________________________________________________________
//__________________________________________________________________________
//Init...
void LCDVTopl::Init() {
iper.ResizeTo(3,3);
iper_inv.ResizeTo(3,3);
vertex_list = new TClonesArray("LCDVToplVRT",1000);
vtrack_list = new TClonesArray("LCDVToplTRK",1000);
vrttrack_list = new TObjArray(100);
InitControlParameter();
dblv=0;
// Initialize jet total momentum
pja.SetXYZ(0.0,0.0,0.0);
ptks.SetXYZ(0.0,0.0,0.0);
//Set Beam Position (Temporary)
ipv.SetXYZ(0.0,0.0,0.0);
}
//__________________________________________________________________________
//__________________________________________________________________________
//Delete LCDVToplVRT, LCDVToplTRK...
void LCDVTopl::Clean() {
if (vtrack_list) vtrack_list->Delete();
if (vertex_list) vertex_list->Delete();
if (vrttrack_list) vrttrack_list->Clear();
}
//__________________________________________________________________________
//__________________________________________________________________________
//Get IP probabirity
Double_t LCDVTopl::GetIPProbabirity(TVector3 xvtx) {
TMatrixD ip2xvtx(3,1);
for (Int_t i=0 ; i < 3 ; i++) { ip2xvtx(i,0)=xvtx[i] - ipv[i]; }
TMatrixD tmp0(iper_inv,TMatrixD::kMult,ip2xvtx);
TMatrixD tmp1(ip2xvtx,TMatrixD::kTransposeMult,tmp0);
Double_t prfu = 0.0;
// protect against floating overflow
if (tmp1(0,0) < 100.0) {
prfu = TMath::Exp(-0.5*tmp1(0,0));
}
return prfu;
}
//__________________________________________________________________________
//__________________________________________________________________________
//Get Initial vertexposition rrack and track
TVector3 LCDVTopl::GetInitialVertexT2T(Int_t itrk1, Int_t itrk2) {
//initial vertex point for given tracks
LCDVToplTRK* trk1= GetTrackAt(itrk1);
TVector3 posv1 = trk1->GetPositionVector(0.0);
// Double_t phi0pr1 = trk1->GetTrackParameter(1)-TMath::Pi();
Double_t phi0pr1 = trk1->GetTrackParameter(1)-TMath::Pi()/2.0;
Double_t cspr1 = TMath::Cos(phi0pr1);
Double_t snpr1 = TMath::Sin(phi0pr1);
Double_t z0pr1 = trk1->GetTrackParameter(3);
Double_t x0pr1 = cspr1*(posv1[0] - ipv[0]) + snpr1*(posv1[1] - ipv[1]);
Double_t tanlpr1 = trk1->GetTrackParameter(4);
Double_t s1spr1 = trk1->GetErrorTsiApprox(0.0);
Double_t s2spr1 = trk1->GetErrorEtaApprox(0.0)*(1.0 + tanlpr1*tanlpr1);
LCDVToplTRK* trk2= GetTrackAt(itrk2);
TVector3 posv2 = trk2->GetPositionVector(0.0);
// Double_t phi0pr2 = trk2->GetTrackParameter(1)-TMath::Pi();
Double_t phi0pr2 = trk2->GetTrackParameter(1)-TMath::Pi()/2.0;
Double_t cspr2 = TMath::Cos(phi0pr2);
Double_t snpr2 = TMath::Sin(phi0pr2);
Double_t z0pr2 = trk2->GetTrackParameter(3);
Double_t x0pr2 = cspr2*(posv2[0] - ipv[0]) + snpr2*(posv2[1] - ipv[1]);
Double_t tanlpr2 = trk2->GetTrackParameter(4);
Double_t s1spr2 = trk2->GetErrorTsiApprox(0.0);
Double_t s2spr2 = trk2->GetErrorEtaApprox(0.0)*(1.0 + tanlpr2*tanlpr2);
TMatrixD m_rot(3,3);
m_rot(0,0) = cspr1*cspr1/s1spr1 + tanlpr1*tanlpr1*snpr1*snpr1/s2spr1 +
cspr2*cspr2/s1spr2 + tanlpr2*tanlpr2*snpr2*snpr2/s2spr2;
m_rot(1,0) = cspr1*snpr1/s1spr1 - tanlpr1*tanlpr1*snpr1*cspr1/s2spr1 +
cspr2*snpr2/s1spr2 - tanlpr2*tanlpr2*snpr2*cspr2/s2spr2;
m_rot(2,0) = tanlpr1*snpr1/s2spr1 + tanlpr2*snpr2/s2spr2;
m_rot(0,1) = m_rot(1,0);
m_rot(1,1) = snpr1*snpr1/s1spr1 + tanlpr1*tanlpr1*cspr1*cspr1/s2spr1 +
snpr2*snpr2/s1spr2 + tanlpr2*tanlpr2*cspr2*cspr2/s2spr2;
m_rot(2,1) = -(tanlpr1*cspr1/s2spr1 + tanlpr2*cspr2/s2spr2);
m_rot(0,2) = m_rot(2,0);
m_rot(1,2) = m_rot(2,1);
m_rot(2,2) = -(z0pr1/s2spr1 + z0pr2/s2spr2);
TMatrixD m_tmp(3,1);
m_tmp(0,0) = -(cspr1*x0pr1/s1spr1 + z0pr1*tanlpr1*snpr1/s2spr1 +
cspr2*x0pr2/s1spr2 + z0pr2*tanlpr2*snpr2/s2spr2);
m_tmp(1,0) = -snpr1*x0pr1/s1spr1 + z0pr1*tanlpr1*cspr1/s2spr1 -
snpr2*x0pr2/s1spr2 + z0pr2*tanlpr2*cspr2/s2spr2;
m_tmp(2,0) = -(z0pr1/s2spr1 + z0pr1/s2spr2);
TMatrixD m_rot_inv(TMatrixD::kInverted,m_rot);
TVector3 xdv=ipv;
TMatrixD m_tmp1(3,1);
m_tmp1.Mult(m_rot_inv,m_tmp);
// xdv.SetXYZ(m_tmp1(0,0)+ipv[0],m_tmp1(1,0)+ipv[1],m_tmp1(2,0));
xdv.SetXYZ(-m_tmp1(0,0)+ipv[0],-m_tmp1(1,0)+ipv[1],-m_tmp1(2,0));
// avoid being too far off
TVector3 xdv2(xdv[0],xdv[1],xdv[2]-ipv[2]);
if ( xdv2.Mag2() > 16.0 ) {
xdv=ipv;
}
return xdv;
}
//__________________________________________________________________________
//__________________________________________________________________________
//Get Initial vertexposition ip and track
TVector3 LCDVTopl::GetInitialVertexIP2T(Int_t itrk) {
//initial vertex point for given tracks
LCDVToplTRK* trk=GetTrackAt(itrk);
TVector3 xpt=trk->GetPositionVector(0.0);
Double_t stk2=trk->GetErrorDistance(0.0,ipv[0],ipv[1],ipv[2]);
TVector3 dx=xpt-ipv;
Double_t adx=dx.Mag();
TMatrixD dddxT(3,1);
for (Int_t i=0 ; i < 3 ; i++) { dddxT(i,0)=dx[i]/adx; }
TMatrixD tmp0(iper,TMatrixD::kMult,dddxT);
TMatrixD tmp1(dddxT,TMatrixD::kTransposeMult,tmp0);
Double_t sip2=tmp1(0,0);
Double_t wtot=1.0/sip2 + 1.0/stk2;
return ipv + (1.0/wtot/stk2)*dx;
}
//__________________________________________________________________________
//__________________________________________________________________________
//Get vertex significance
Double_t LCDVTopl::GetVertexSignificance(TVector3 xvrt) {
//return vertex significance V(r) for NTRK tracks at x,y,z
TVector3 ip2xvrt=xvrt-ipv;
Double_t dlong=pja.Mag()> 0 ? (ip2xvrt*pja)/pja.Mag() : ip2xvrt.Mag();
// cut out if behind ip or > 10.0 cm in front
if (dlong > 10.0 || dlong < -0.01) {
return -1.0;
}
Double_t dmag = ip2xvrt.Mag();
Double_t dtran= 0.0;
if (TMath::Abs(dlong) < TMath::Abs(dmag)) {
dtran = TMath::Sqrt(dmag*dmag - dlong*dlong);
}
Double_t prsumu = 0.0;
Double_t prsum2u = 0.0;
Double_t prfu=0.0;
Int_t ntrk=GetTrackEntries();
LCDVToplTRK* trk=0;
Int_t i=0;
for (i=0 ; i < ntrk ; i++) {
trk =GetTrackAt(i);
if (trk->GetTrackID() < 0) { continue; } //skip IP track...
//prfu=trk->GetTrackProbabirity(xvrt);
prfu=trk->GetTrackProbabirity(xvrt,ipv);
prsumu += prfu;
prsum2u += prfu*prfu;
} //End of the track loop.
// add in IP information if required
prfu = GetIPProbabirity(xvrt);
prsumu += kipw*prfu;
prsum2u += kipw*prfu*prfu;
Double_t vsg = 0.0;
if (prsumu > 0.0) {
vsg = prsumu - prsum2u/prsumu;
};
//include extra factors such as jet core weighting
// What are 0.01 and 0.005?
if (dtran > 0.005) {
Double_t alpha = TMath::ACos((dlong + 0.01)/
TMath::Sqrt(TMath::Power(dlong + 0.01,2)
+ TMath::Power(dtran - 0.005,2)));
//What is 0.01 and 0.005????
vsg *= TMath::Exp(-kang*alpha*alpha);
}
return TMath::Max(vsg,0.0);
}
//__________________________________________________________________________
//__________________________________________________________________________
//Get min(Vr/min(V))
Double_t LCDVTopl::GetMinVrOverMinVs(LCDVToplVRT* vrt1, LCDVToplVRT* vrt2, Int_t nstp) {
Double_t v1 =vrt1->GetVsig();
Double_t v2 =vrt2->GetVsig();
Double_t vmin =TMath::Min(v1,v2);
TVector3 xvrt1=vrt1->GetVertexVector();
TVector3 xvrt2=vrt2->GetVertexVector();
TVector3 dxvrt=xvrt2-xvrt1;
// prevent having to find vmin if vertices are within 10 um
Double_t dstep =dxvrt.Mag();
Double_t vrat =0;
Double_t vmaxi =0;
Double_t vmini =0;
Int_t mstep =nstp;
TVector3 xpt(0.0,0.0,0.0);
if (dstep < dcut) {
vrat = 1.0;
} else {
dstep /= nstp + 1;
if (dstep < dcut/2.0) {
mstep=(Int_t)((nstp+1)*dstep/dcut*2.0) - 1;
dstep = (nstp + 1)*dstep/(mstep + 1);
}
dxvrt.SetMag(dstep);
vmini=vmin;
mstep=TMath::Max(mstep,1);
for (Int_t step=1 ; step <= mstep ; step++) {
xpt =xvrt1 + step*dxvrt;
vmaxi=GetVertexSignificance(xpt);
if (vmaxi < vmini) { vmini = vmaxi; }
} //endloop over steps
vrat = vmini/vmin;
}
return vrat;
}
//__________________________________________________________________________
//__________________________________________________________________________
Double_t LCDVTopl::GetMinVrOverMinVs(LCDVToplVRT* vrt1, LCDVToplVRT* vrt2) {
return GetMinVrOverMinVs(vrt1,vrt2,nstep);
}
//__________________________________________________________________________
//__________________________________________________________________________
//Get Vertex by simple iteration
TVector3 LCDVTopl::GetVertexBySimpleIteration(LCDVToplVRT* vrt) {
// ************ 3D iterative bit ****************
TVector3 spx = vrt->GetVertexVector();
Double_t vm3d = vrt->GetVsig();
//start with 4 um step size
Double_t stepsize = 0.0005;
Bool_t f_iteration = kFALSE;
Bool_t f_again = kFALSE;
Double_t vsg = 0.0;
while(kTRUE) {
f_iteration = kFALSE;
while(1) {
spx += TVector3(stepsize,0.0,0.0);
vsg=GetVertexSignificance(spx);
if (vsg > vm3d) {
vm3d = vsg;
f_iteration = kTRUE;
} else {
spx -= TVector3(stepsize,0.0,0.0);
break;
}
}
if (!f_iteration) {
while(1) {
spx -= TVector3(stepsize,0.0,0.0);
vsg=GetVertexSignificance(spx);
if (vsg > vm3d) {
vm3d = vsg;
f_iteration = kTRUE;
} else {
spx += TVector3(stepsize,0.0,0.0);
break;
}
}
}
f_again=kFALSE;
//now y co-ordinate
f_iteration = kFALSE;
while(1) {
spx += TVector3(0.0,stepsize,0.0);
vsg=GetVertexSignificance(spx);
if (vsg > vm3d) {
f_again = kTRUE;
vm3d = vsg;
f_iteration = kTRUE;
} else {
spx -= TVector3(0.0,stepsize,0.0);
break;
}
}
if (!f_iteration) {
while(1) {
spx -= TVector3(0.0,stepsize,0.0);
vsg=GetVertexSignificance(spx);
if (vsg > vm3d) {
f_again = kTRUE;
vm3d = vsg;
} else {
spx += TVector3(0.0,stepsize,0.0);
break;
}
}
}
//now z co-ordinate
f_iteration = kFALSE;
while(1) {
spx += TVector3(0.0,0.0,stepsize);
vsg=GetVertexSignificance(spx);
if (vsg > vm3d) {
f_again = kTRUE;
vm3d = vsg;
f_iteration = kTRUE;
} else {
spx -= TVector3(0.0,0.0,stepsize);
break;
}
}
if (!f_iteration) {
while(1) {
spx -= TVector3(0.0,0.0,stepsize);
vsg=GetVertexSignificance(spx);
if (vsg > vm3d) {
f_again = kTRUE;
vm3d = vsg;
} else {
spx += TVector3(0.0,0.0,stepsize);
break;
}
}
}
if (!f_again) {
if (stepsize > dcut) {
stepsize = stepsize/2.0;
} else {
break;
}
}
}
vrt->SetVsig(vm3d);
vrt->SetVertex(spx);
return spx;
}
//__________________________________________________________________________
//__________________________________________________________________________
//Get Vertex Charge
Double_t LCDVTopl::GetVertexCharge() {
Int_t nvtrks=GetVrtTrackEntries();
LCDVToplTRK* trk=0;
Double_t q_vertex=0;
for (Int_t itrk=0 ; itrk < nvtrks ; itrk++) {
trk=GetVrtTrackAt(itrk);
q_vertex += trk->GetCharge();
}
return q_vertex;
}
//__________________________________________________________________________
//__________________________________________________________________________
//Get # of tracks
Int_t LCDVTopl::GetTrackEntries(){
return vtrack_list->GetEntries();
}
//__________________________________________________________________________
//__________________________________________________________________________
LCDVToplTRK* LCDVTopl::GetTrackAt(Int_t itrk) {
return (LCDVToplTRK*)(vtrack_list->At(itrk));
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::GetVertexEntries(){
return vertex_list->GetEntries();
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::GetVertexEntriesExceptV0(){
Int_t nvrt=GetVertexEntries();
Int_t nvrtnov0=0;
LCDVToplVRT* vrt;
for (Int_t ivrt=0 ; ivrt < nvrt ; ivrt++) {
vrt=GetVertexAt(ivrt);
if (CheckKs0Vertex(vrt)) continue;
if (CheckLambda0Vertex(vrt)) continue;
nvrtnov0++;
}
return nvrtnov0;
}
//__________________________________________________________________________
//__________________________________________________________________________
LCDVToplVRT* LCDVTopl::GetVertexAt(Int_t ivrt) {
return (LCDVToplVRT*)(vertex_list->At(ivrt));
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::GetVrtTrackEntries(){
return vrttrack_list->GetEntries();
}
//__________________________________________________________________________
//__________________________________________________________________________
LCDVToplTRK* LCDVTopl::GetVrtTrackAt(Int_t itrk) {
return (LCDVToplTRK*)(vrttrack_list->At(itrk));
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::InitControlParameter(){
// ripe = 3.0;
// zipe = 10.0;
ripe = 3.0e-4;
zipe = 6.0e-4;
rcut = 0.6;
xcut = 10.0;
kang = 5.0;
kipw = 1.0;
nstep = 7;
dcut = 0.0002;
m_SeedVrtRMax =2.3 ;
m_SeedVrtK0Mass=0.015;
m_SeedVrtLambda0Mass=0.015;
m_SeedVrtSigMin=2.0 ;
m_TrdiMax=0.1;
m_LodiMin=0.05;
m_LodrMin=0.25;
//m_LodrMax=2.0;
m_LodrMax=100.0;
//m_AntaMax=0.5;
m_AntaMax=TMath::Pi()/2.0;
//m_PTIP[0]=0.0010;
//m_PTIP[1]=0.0040;
m_PTIP[0]=1.5*ripe;
m_PTIP[1]=1.5*zipe;
m_NVsigma=1.0;
m_Rmass =0.0;
m_PtMass =0.0;
m_PtMassEr=0.0;
m_Masstag =0.0;
m_MissPt =0.0;
m_UseNN =kFALSE;
m_SeedVrtNN =0.7;
m_TrackAttNN=0.4;
m_SelectNN =0.5;
m_p4Vertex.SetXYZT(0.,0.,0.,0.);
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::SetLCDVToplTRKs(TObjArray* track_list,
Int_t ntrk,
Int_t* track_indecies) {
//This loop to
// make sure position and momentum in trk at 2d POCA
// and to make LCDVToplTRK.
// pja.SetXYZ(0.0,0.0,0.0);
for( Int_t i=0; i < ntrk ; i++ ) {
new((*vtrack_list)[i]) LCDVToplTRK(track_list,track_indecies[i],ipv);
} //loop over zxtrks
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::CalcIPErrorMatrix(){
//Double_t ripe2=ripe/1e8*ripe;
//Double_t zipe2=zipe/1e8*zipe;
Double_t ripe2=ripe*ripe;
Double_t zipe2=zipe*zipe;
iper(0,0)=iper_inv(0,0)=ripe2;
iper(1,0)=iper_inv(1,0)=0.0 ;
iper(2,0)=iper_inv(2,0)=0.0 ;
iper(0,1)=iper_inv(0,1)=0.0 ;
iper(1,1)=iper_inv(1,1)=ripe2;
iper(2,1)=iper_inv(2,1)=0.0 ;
iper(0,2)=iper_inv(0,2)=0.0 ;
iper(1,2)=iper_inv(1,2)=0.0 ;
iper(2,2)=iper_inv(2,2)=zipe2;
Double_t det_iper;
iper_inv.Invert(&det_iper);
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::DeleteTracks(){
vtrack_list->Delete();
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::DeleteVertices(){
vertex_list->Delete();
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::AddVertex(LCDVToplVRT* vrt) {
Int_t nvrt=GetVertexEntries();
Int_t ivrt=0;
Bool_t f_no_same_vertex=kTRUE;
for (ivrt=0 ; ivrt < nvrt-1 ; ivrt++) {
if (vrt == GetVertexAt(ivrt)) {
f_no_same_vertex=kFALSE;
break;
}
}
if (f_no_same_vertex) {
SortVertexByVsig();
} else {
vertex_list->Remove(vrt);
}
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::AddVertexIntoTrack(LCDVToplTRK* trk, LCDVToplVRT* vrt) {
Int_t nvrt=trk->GetVertexEntries();
Int_t ivrt=0;
Bool_t f_no_same_vertex=kTRUE;
LCDVToplVRT* tvrt=0;
for (ivrt=0; ivrt < nvrt ; ivrt++) {
tvrt=(LCDVToplVRT*)trk->GetVertexAt(ivrt);
if (tvrt == vrt) {
f_no_same_vertex=kFALSE;
break;
}
}
if (f_no_same_vertex) {
trk->AddVertex(vrt);
SortVertexByVsigInTrack(trk);
}
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::AddVrtTrack(LCDVToplTRK* trk) {
vrttrack_list->Add(trk);
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::CompressVerticesList() {
vertex_list->Compress();
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::RemoveVertexAt(Int_t i) {
RemoveVertex(GetVertexAt(i));
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::RemoveVertex(LCDVToplVRT* vrt) {
vertex_list->Remove(vrt);
CompressVerticesList();
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::MergeVertex(LCDVToplVRT* vrti,LCDVToplVRT* vrtj,Int_t f_fit) {
if (vrti == vrtj) { return; }
LCDVToplTRK* trk=0;
Int_t ntrk=vrtj->GetTrackEntries();
Int_t itrk=0;
for (itrk=0; itrk < ntrk ; itrk++) {
trk =vrtj->GetTrackAt(itrk);
vrti->AddTrack(trk,f_fit);
trk->RemoveVertex(vrtj);
AddVertexIntoTrack(trk,vrti);
trk->CompressVerticesList();
}
if (f_fit > 1) {
vrti->SetVsig(GetVertexSignificance(vrti->GetVertexVector()));
}
RemoveVertex(vrtj);
CompressVerticesList();
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::MergeVertex(Int_t i, Int_t j, Int_t f_fit) {
if (i == j) { return; }
LCDVToplVRT* vrti=GetVertexAt(i);
LCDVToplVRT* vrtj=GetVertexAt(j);
MergeVertex(vrti,vrtj,f_fit);
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::SortVertexByVsig() {
vertex_list->Sort();
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::SortVertexByDecayLength() {
Int_t nvrt=GetVertexEntries();
Int_t ivrti,ivrtj;
Double_t disti,distj;
LCDVToplVRT* vrti;
LCDVToplVRT* vrtj;
LCDVToplVRT vrtk;
for (ivrti=0 ; ivrti < nvrt-1 ; ivrti++) {
for (ivrtj=ivrti+1 ; ivrtj < nvrt ; ivrtj++) {
vrti=GetVertexAt(ivrti);
vrtj=GetVertexAt(ivrtj);
disti=vrti->GetDecayLength(ipv);
distj=vrtj->GetDecayLength(ipv);
if (disti > distj) {
vrtk = *vrtj;
*vrtj = *vrti;
*vrti = vrtk;
}
}
}
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::SortVertexByVsigInTrack(LCDVToplTRK* trk) {
Int_t nvrt=trk->GetVertexEntries();
Int_t ivrti=0,ivrtj=0;
LCDVToplVRT* vrti=0;
LCDVToplVRT* vrtj=0;
for (ivrti=0 ; ivrti < nvrt-1 ; ivrti++) {
for (ivrtj=ivrti+1 ; ivrtj < nvrt ; ivrtj++) {
vrti=(LCDVToplVRT*)(trk->GetVertexAt(ivrti));
vrtj=(LCDVToplVRT*)(trk->GetVertexAt(ivrtj));
if (vrti->GetVsig() < vrtj->GetVsig()) {
// trk->AddVertexAt((TObject*)vrtj,ivrti);
// trk->AddVertexAt((TObject*)vrti,ivrtj);
trk->AddVertexAt(vrtj,ivrti);
trk->AddVertexAt(vrti,ivrtj);
}
}
}
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::SortVertexByVsigInAllTrack() {
Int_t ntrk=GetTrackEntries();
for (Int_t itrk=0 ; itrk < ntrk ; itrk++) {
SortVertexByVsigInTrack(GetTrackAt(itrk));
}
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::FindVertices() {
return FindVertices(0);
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::FindVertices(Int_t f_v0check) {
//-------------------------------------------------------
//
// FindVertices --
// determine the TOPOlogical Vertex structure of a jet
//------------------------------------------------------
//
// resolve_algo -- resolve track V(r) into inclusive vertices
//-----------------------------------------------------
if (dblv > 0)
cout << "=========<FindVertices started!>============"<<endl;
Int_t ntrk=GetTrackEntries();
if (dblv > 0)
cout << " ---> Check 1: N of LCDVToplTRKs= " << ntrk << endl;
// first loop over track pairs to find Vs and interactions
if (ntrk < 2) {
//not enough tracks!
return 0;
}
if (dblv > 0)
cout <<" ---> Check 2: CalcIPErrorMatrix() ... ";
CalcIPErrorMatrix();
if (dblv > 0)
cout <<"done"<<endl;
// Addd IP-Track...
if (dblv > 0)
cout << " ---> Check 3: AddIPTrack ... ";
LCDVToplTRK* iptrk=new((*vtrack_list)[ntrk]) LCDVToplTRK();
if (dblv > 0)
cout <<" done"<<endl;
//----------------------------------------------------
// now loop over track pairs and IP
//----------------------------------------------------
if (dblv > 0)
cout << " ---> Check 4: Making Vertices with IP-Track ... ";
TVector3 xvinit(0.0,0.0,0.0);
TVector3 xvrt(0.0,0.0,0.0);
LCDVToplTRK* trkj = 0;
LCDVToplTRK* trk = 0;
LCDVToplVRT* vrt = 0;
Double_t vmaxi= 0.0;
Double_t prfi = 0;
Double_t prfj = 0;
Double_t lprfi= 0;
Double_t lprfj= 0;
Int_t itrk = ntrk;
Int_t jtrk = 0;
Int_t nvrt = 0;
for (jtrk=0 ; jtrk < ntrk ;jtrk++) {
trkj=GetTrackAt(jtrk);
//ip-track combination
xvinit=GetInitialVertexIP2T(jtrk);
vrt=new((*vertex_list)[nvrt]) LCDVToplVRT(xvinit);
vrt->AddTrack(iptrk);
vrt->AddTrack(trkj);
vrt->SetVertex(xvinit);
xvrt=vrt->GetVertexVector();
vmaxi=GetVertexSignificance(xvrt);
//Vertex significance check.
if (vmaxi <= 0.001) {
vertex_list->Remove(vrt);
continue;
}
vrt->SetVsig(vmaxi);
prfi=GetIPProbabirity(xvrt);
prfj=trkj->GetTrackProbabirity(xvrt,ipv);
//Track Probabirity Check.
if (prfi*prfj <= 0.0) {
vertex_list->Remove(vrt);
continue;
}
lprfi = -2.0*TMath::Log(prfi);
lprfj = -2.0*TMath::Log(prfj);
//Chi**2 contribution of tracks check.
if (lprfi > xcut || lprfj > xcut) {
vertex_list->Remove(vrt);
continue;
}
nvrt++;
AddVertexIntoTrack(iptrk,vrt);
AddVertexIntoTrack(trkj ,vrt);
} // end do loop over track pairs
//Int_t nvrt=GetVertexEntries();
if (dblv > 0)
cout <<nvrt<<" done"<<endl;
//----------------------------------------------------
// now loop over track pairs
//----------------------------------------------------
if (dblv > 0)
cout << " ---> Check 5: Making Vertices by Track-Track ... ";
LCDVToplTRK* trki;
for (itrk=1 ; itrk < ntrk ; itrk++) {
trki=GetTrackAt(itrk);
for (jtrk=0 ; jtrk < itrk ;jtrk++) {
trkj=GetTrackAt(jtrk);
//track-track combination
xvinit=GetInitialVertexT2T(itrk,jtrk);
vrt=new((*vertex_list)[nvrt]) LCDVToplVRT(xvinit);
vrt->AddTrack(trki);
vrt->AddTrack(trkj);
xvrt=vrt->GetVertexVector();
prfi=trki->GetTrackProbabirity(xvrt,ipv);
prfj=trkj->GetTrackProbabirity(xvrt,ipv);
//Track Probabirity Check.
if (prfi*prfj <= 0.0) {
vertex_list->Remove(vrt);
continue;
}
lprfi = -2.0*TMath::Log(prfi);
lprfj = -2.0*TMath::Log(prfj);
//Chi**2 contribution of tracks check.
if (lprfi > xcut || lprfj > xcut ) {
vertex_list->Remove(vrt);
continue;
}
vmaxi=GetVertexSignificance(xvrt);
vrt->SetVsig(vmaxi);
//Vertex significance check.
if (vmaxi <= 0.001) {
vertex_list->Remove(vrt);
continue;
}
//V0 check...
if (f_v0check != 0 && vrt->GetChi2() < 16.0) {
vrt->CalcFourVector();
if (CheckKs0Vertex(vrt)) {
trki->SetV0flag(1);
trkj->SetV0flag(1);
}
if (CheckLambda0Vertex(vrt)) {
trki->SetV0flag(2);
trkj->SetV0flag(2);
}
}
// xvrt=GetVertexBySimpleIteration(vrt);
nvrt++;
AddVertexIntoTrack(trki,vrt);
AddVertexIntoTrack(trkj,vrt);
} // end do loop over track pairs
} //end do loop over track pairs
SortVertexByVsig();
SortVertexByVsigInAllTrack();
//nvrt=GetVertexEntries();
if (dblv > 0)
cout <<nvrt<<" done"<<endl;
Int_t ivrt;
Int_t nvtrk = 0;
PrintOutForDebug();
// -----------> now resolve vmax for each track
if (dblv > 0)
cout << " ---> Check 6: Resolve Vertices per track ... ";
Double_t gmax= 0.0;
Double_t vmax=-1.0;
Int_t itsp=0,ntsp=0;
Int_t irsp=0,nrsp=0;
TVector3 xvrt_prev(0.0,0.0,0.0);
TVector3 dvrt(0.0,0.0,0.0);
Bool_t f_resolved=kTRUE;
LCDVToplVRT* vrt_prev=0;
Double_t vrat=0.0;
for (itrk=0 ; itrk < ntrk+1 ; itrk++ ) {
//loop over tracks + ip
trk=GetTrackAt(itrk);
gmax = 0.0; // global max on track
vmax =-1.0;
nrsp = 0;
ntsp = trk->GetVertexEntries();
itsp=0;
while (itsp < ntsp) { //loop over track's spatial points
vrt=(LCDVToplVRT*)(trk->GetVertexAt(itsp));
vmax = vrt->GetVsig();
if (nrsp == 0) { gmax = vmax; }
if (vmax > 0.001 && vmax > 0.1*gmax) { // must be >frac of big max
xvrt=vrt->GetVertexVector();
f_resolved=kTRUE;
// no. of res. spatial points on track
for (irsp=0 ; irsp < nrsp ; irsp++) {
vrt_prev=(LCDVToplVRT*)(trk->GetVertexAt(irsp));
xvrt_prev=vrt_prev->GetVertexVector();
dvrt=xvrt - xvrt_prev;
if (dvrt.Mag() < dcut) {
f_resolved=kFALSE;
break;
}
vrat =GetMinVrOverMinVs(vrt,vrt_prev,7);
// vrat =GetMinVrOverMinVs(vrt,vrt_prev,15);
if (vrat > rcut) {
f_resolved=kFALSE;
break;
}
} //endloop over resolved max.
if (f_resolved) {
nrsp++;
itsp++;
} else {
MergeVertex(vrt_prev,vrt,0);
ntsp=trk->GetVertexEntries();
} //endif f_resolved
} else {
trk->RemoveVertexAt(itsp);
vrt->RemoveTrack(trk,0);
ntsp=trk->GetVertexEntries();
} // end if vmax > 0.001 && vmax > 0.1*gmax
} //endloop over spatial points on track.
} //endloop over tracks
nvrt=GetVertexEntries();
if (dblv > 0)
cout <<nvrt<<" done"<<endl;
PrintOutForDebug();
//fill resolved spatial max. arrays
if (dblv > 0)
cout << " ---> Check 7: Resolve Vertices by Track info. ";
Int_t nsp = GetVertexEntries();
Int_t isp = 0;
Int_t nresp = 0;
Int_t iresp = 0;
LCDVToplVRT* tvrt = 0;
Bool_t f_not_found_yet=kTRUE;
// Int_t k=0;
isp=0;
while(isp < nsp) {
// for (isp=0 ; isp < nsp ; isp++) {
vrt=GetVertexAt(isp);
nvtrk=vrt->GetTrackEntries();
if (nvtrk > 1) {
//use logical f_resolved for convenience
f_resolved=kFALSE;
for (itrk=0 ; itrk < ntrk+1 ; itrk++) {
trk =GetTrackAt(itrk);
nrsp=trk->GetVertexEntries();
for (irsp=0 ; irsp < nrsp ; irsp++) {
tvrt=(LCDVToplVRT*)(trk->GetVertexAt(irsp));
if (vrt == tvrt) {
f_resolved=kTRUE;
break;
}
}
if (f_resolved) { break; }
}
if (f_resolved && vrt->GetVsig() > 0.001) {
xvrt=GetVertexBySimpleIteration(vrt);
//has this spatial max already been found
f_not_found_yet=kTRUE;
for (iresp=0 ; iresp < nresp ; iresp++) {
vrt_prev =GetVertexAt(iresp);
xvrt_prev=vrt_prev->GetVertexVector();
dvrt =xvrt - xvrt_prev;
if (dvrt.Mag() < dcut) {
//if old spatial point, assign tracks to it
nvtrk=vrt->GetTrackEntries();
for (itrk=0 ; itrk < nvtrk ; itrk++) {
trk =vrt->GetTrackAt(itrk);
nrsp=trk->GetVertexEntries();
for (irsp=0; irsp < nrsp ; irsp++) {
tvrt = (LCDVToplVRT*)(trk->GetVertexAt(irsp));
if (vrt == tvrt) {
trk->RemoveVertexAt(irsp);
AddVertexIntoTrack(trk,vrt_prev);
//T.Abe Add 3/16/2001
vrt_prev->AddTrack(trk);
}
}
}
RemoveVertexAt(isp);
nsp=GetVertexEntries();
f_not_found_yet=kFALSE;
break;
}
}
if (f_not_found_yet) {
nresp++;
isp++;
}
} else {
RemoveVertexAt(isp);
nsp=GetVertexEntries();
}// endif ....
} else { // # of tracks in the vertex < 2
if (nvtrk > 0) {
trk=vrt->GetTrackAt(0);
// trk->RemoveVertex((TObject*)vrt);
trk->RemoveVertex(vrt);
}
RemoveVertexAt(isp);
nsp=GetVertexEntries();
} // endif nvtrk > 1
} // endloop of isp
nvrt=GetVertexEntries();
SortVertexByVsig();
SortVertexByVsigInAllTrack();
if (dblv > 0)
cout <<nvrt<<" done"<<endl;
PrintOutForDebug();
//
// loop over spatial points, fill resolution matrix Vminsp
//
// now find resolved vertices within rcut
if (dblv > 0)
cout << " ---> Check 8: Resolve Vertices by Vsig... ";
nresp=GetVertexEntries();
LCDVToplVRT* vrti=0;
LCDVToplVRT* vrtj=0;
iresp=0;
Int_t jresp=0;
while(iresp < nresp-1) {
vrti=GetVertexAt(iresp);
jresp = iresp + 1;
while(jresp < nresp) {
vrtj=GetVertexAt(jresp);
vrat=GetMinVrOverMinVs(vrti,vrtj,7);
// vrat=GetMinVrOverMinVs(vrti,vrtj,15);
if (vrat > rcut) {
MergeVertex(iresp,jresp,0);
nresp=GetVertexEntries();
} else {
jresp++;
}
}
iresp++;
}
nvrt=GetVertexEntries();
SortVertexByVsig();
SortVertexByVsigInAllTrack();
if (dblv > 0)
cout <<nvrt<<" done"<<endl;
PrintOutForDebug();
//
//now fit vertices and do chi**2 trimming
//
if (dblv > 0)
cout << " ---> Check 9: Chi**2 triming ... ";
nvrt=GetVertexEntries();
Int_t ntracks_in_vertex =0;
//Int_t nvertices_in_track=0;
Int_t itrkmx =0;
Double_t maxchi2 =0.0;
Double_t pchi2 =0.0;
Int_t itvrt =0;
ivrt = 0;
Bool_t f_inc_iptrk;
while (ivrt < nvrt) {
vrt=GetVertexAt(ivrt);
xvrt=vrt->GetVertexVector();
vrt->VertexFit(xvrt);
while(1) {
ntracks_in_vertex=vrt->GetTrackEntries();
f_inc_iptrk = (vrt->IndexOfTrack(iptrk) >= 0) ? kTRUE : kFALSE;
if (ntracks_in_vertex < 2) {
for (itrk=0 ; itrk < ntracks_in_vertex ; itrk++) {
trk=vrt->GetTrackAt(itrk);
// trk->RemoveVertex((TObject*)vrt);
trk->RemoveVertex(vrt);
}
RemoveVertex(vrt);
nvrt=GetVertexEntries();
break;
} else {
// xvinit = ipv;
itrkmx =0;
maxchi2=0.0;
if (ntracks_in_vertex ==2 && f_inc_iptrk) {
trk=0;
for (itrk=0; itrk < ntracks_in_vertex ; itrk++) {
trk=vrt->GetTrackAt(itrk);
if (trk != iptrk) {
break;
}
}
xvrt=GetInitialVertexIP2T(vtrack_list->IndexOf(trk));
vmaxi=GetVertexSignificance(xvrt);
vrt->SetVertex(xvrt);
vrt->SetVsig(vmaxi);
} else {
xvrt=vrt->GetVertexVector();
}
for (itrk=0; itrk < ntracks_in_vertex ; itrk++) {
trk=vrt->GetTrackAt(itrk);
if (trk != iptrk) {
if ( (ntracks_in_vertex >= 3 && !f_inc_iptrk) ||
(ntracks_in_vertex > 3 && f_inc_iptrk)) {
pchi2=vrt->GetPchi2(itrk);
} else {
pchi2=-2.0*
TMath::Log(TMath::Max(1e-40,
trk->GetTrackProbabirity(xvrt,ipv)));
}
} else {
if ( (ntracks_in_vertex >= 3 && !f_inc_iptrk) ||
(ntracks_in_vertex > 3 && f_inc_iptrk)) {
pchi2=TMath::Power(vrt->GetDecayLength(ipv),2)
/vrt->GetErrorDecayLength(ipv,iper);
} else {
pchi2=-2.0*TMath::Log(TMath::Max(1e-40,GetIPProbabirity(xvrt)));
}
}
if (pchi2 > maxchi2) {
maxchi2=pchi2;
itrkmx =itrk;
}
}
if (maxchi2 > xcut) {
//Remove this vertex from track info.
trk=vrt->GetTrackAt(itrkmx);
// trk->RemoveVertex((TObject*)vrt);
trk->RemoveVertex(vrt);
//Drop this track
vrt->RemoveTrackAt(itrkmx);
} else {
ivrt++;
break;
} // ifmaxchi2....
} //more than 1 track in vertex
} // loop over for chi**2 triming
} //loop over vertices
nvrt=GetVertexEntries();
for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
vrt=GetVertexAt(ivrt);
vmaxi=GetVertexSignificance(vrt->GetVertexVector());
vrt->SetVsig(vmaxi);
}
SortVertexByVsig();
SortVertexByVsigInAllTrack();
if (dblv > 0)
cout <<nvrt<<" done"<<endl;
PrintOutForDebug();
//
// now decide on track ambiguities by starting with highest spat. max in vrt
//
if (dblv > 0)
cout << " ---> Check 10 : Solve track ambiguities ... ";
for (itrk=0 ; itrk < ntrk+1 ; itrk++ ) {
//loop over tracks + ip
trk = GetTrackAt(itrk);
vrti = 0;
nvrt = trk->GetVertexEntries();
for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
vrt=(LCDVToplVRT*)(trk->GetVertexAt(ivrt));
ntracks_in_vertex=vrt->GetTrackEntries();
f_inc_iptrk = (vrt->IndexOfTrack(iptrk) >= 0) ? kTRUE : kFALSE;
//require at least 2 trk in the vertex or in ip...
if ( (ntracks_in_vertex > 1 && !f_inc_iptrk) || f_inc_iptrk) {
vrti=vrt;
break;
}
}
itvrt=nvrt;
while(--itvrt > 0) {
vrtj=(LCDVToplVRT*)(trk->GetVertexAt(itvrt));
if (vrtj == vrti) {
vrtj = (LCDVToplVRT*)(trk->GetVertexAt(0));
trk->AddVertexAt(vrtj,itvrt);
trk->AddVertexAt(vrti,0 );
}
vrtj->RemoveTrack(trk,0);
trk->RemoveVertexAt(itvrt);
}
}
nvrt=GetVertexEntries();
if (dblv > 0)
cout<<nvrt<<" done"<<endl;
PrintOutForDebug();
//
//re-sort vertices if at least two tracks or ip
//
if (dblv > 0)
cout << " ---> Check 11: Resort Vertices ...";
nresp = 0;
iresp = 0;
ivrt=0;
while(ivrt < nvrt) {
vrt=GetVertexAt(ivrt);
nvtrk=vrt->GetTrackEntries();
f_inc_iptrk = (vrt->IndexOfTrack(iptrk) >= 0) ? kTRUE : kFALSE;
//require at least 2 trk in the vertex or in ip...
if ( (nvtrk > 1 && !f_inc_iptrk) || f_inc_iptrk) {
vrt->VertexFit(ipv);
ivrt++;
} else {
for (itrk=0 ; itrk < nvtrk ; itrk++) {
trk=vrt->GetTrackAt(itrk);
// trk->RemoveVertex((TObject*)vrt);
trk->RemoveVertex(vrt);
}
RemoveVertex(vrt);
nvrt=GetVertexEntries();
}
} // endloop of ivrt
nvrt=GetVertexEntries();
if (dblv > 0)
cout <<nvrt<<" done"<<endl;
PrintOutForDebug();
// now find resolved vertices within rcut
if (dblv > 0)
cout << " ---> Check 12: Resolve Vertices by Vsig (again)... ";
nresp=GetVertexEntries();
iresp=0;
jresp=0;
while(iresp < nresp-1) {
vrti=GetVertexAt(iresp);
jresp = iresp + 1;
while(jresp < nresp) {
vrtj=GetVertexAt(jresp);
vrat=GetMinVrOverMinVs(vrti,vrtj,7);
// vrat=GetMinVrOverMinVs(vrti,vrtj,15);
if (vrat > rcut) {
MergeVertex(iresp,jresp,0);
nresp=GetVertexEntries();
} else {
jresp++;
}
}
iresp++;
}
Int_t nvrt_prev=nvrt;
nvrt=GetVertexEntries();
if (nvrt < nvrt_prev) {
for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
vrt=GetVertexAt(ivrt);
xvrt=vrt->GetVertexVector();
vrt->VertexFit(xvrt);
}
SortVertexByVsig();
//SortVertexByVsigInAllTrack();
}
if (dblv > 0)
cout <<nvrt<<" done"<<endl;
PrintOutForDebug();
//
//remove IP track...
//
if (dblv > 0)
cout << " ---> Check 13: Check IP ... ";
Bool_t f_no_ipvrt=kTRUE;
for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
vrt=GetVertexAt(ivrt);
ntracks_in_vertex=vrt->GetTrackEntries();
//remove IP-track from vertices.
itrk=0;
while(itrk < ntracks_in_vertex) {
trk=vrt->GetTrackAt(itrk);
if (trk == iptrk) {
f_no_ipvrt=kFALSE;
vrt->RemoveTrackAt(itrk);
ntracks_in_vertex=vrt->GetTrackEntries();
} else {
itrk++;
}
}
}
vtrack_list->Remove(iptrk);
vtrack_list->Compress();
//if (nvrt <= 0) {
if (f_no_ipvrt || nvrt <= 0) {
vrt=new((*vertex_list)[nvrt]) LCDVToplVRT(ipv);
vmaxi=GetVertexSignificance(ipv);
vrt->SetVsig(vmaxi);
for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
vrti=GetVertexAt(ivrt);
vrat =GetMinVrOverMinVs(vrt,vrti,7);
if (vrat > rcut) {
vrt=0;
break;
}
}
if (vrt == 0)
RemoveVertex(vrt);
}
SortVertexByDecayLength();
nvrt=GetVertexEntries();
if (dblv > 0)
cout <<nvrt<<" done"<<endl;
PrintOutForDebug();
//
//save info.
//
if (dblv > 0)
cout << " ---> Check 14: Save info ... ";
for (ivrt=0 ; ivrt < nvrt ; ivrt++) {
vrt=GetVertexAt(ivrt);
vrt->CalcFourVector();
}
nvrt=GetVertexEntries();
if (dblv > 0)
cout <<nvrt<<" done"<<endl;
if (nvrt == 0) { cout << " nvrt == 0 @ Check 12 " << ntrk << endl; }
return 1;
}
//__________________________________________________________________________
//-------------------------------------------
// PT corrected mass relates...
//-------------------------------------------
//__________________________________________________________________________
Int_t LCDVTopl::GetSeedVertex() {
Int_t seedvertex_id=0;
if (m_UseNN) {
seedvertex_id=GetSeedVertexByNNMethod();
} else {
seedvertex_id=GetSeedVertexByClassicMethod();
}
return seedvertex_id;
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::GetSeedVertexByClassicMethod() {
Int_t nseedvrt = -1;
Int_t nvrt = GetVertexEntries();
if (nvrt < 2) return nseedvrt;
Double_t maxvsig = -1.;
LCDVToplVRT* vrt=GetVertexAt(0);
for(Int_t ivrt =1; ivrt < nvrt ; ivrt++){
vrt = GetVertexAt(ivrt);
// Apply cuts to select seed vertex candidates
// Cut 1
Double_t decayL = vrt->GetDecayLength(ipv);
if (decayL > m_SeedVrtRMax) continue;
// Cut 2
Int_t ntrks = vrt->GetTrackEntries();
Double_t netCharge = vrt->GetVertexCharge();
if (ntrks == 2 && netCharge == 0.) {
Double_t reconMass = vrt->GetMass();
if ( TMath::Abs(reconMass - 0.49767) < m_SeedVrtK0Mass ) continue;
}
// Cut 3
TVector3 vrtpos(vrt->GetVertexVector() - ipv);
if (vrtpos*pja < 0.) continue;
// Cut 4
TVector3 vrtmon(vrt->GetFourVector().Vect());
if (vrtpos*vrtmon < 0.) continue;
Double_t decayLerr = TMath::Sqrt(vrt->GetErrorDecayLength(ipv,iper));
Double_t vrtsig = decayL/decayLerr;
if (vrtsig < m_SeedVrtSigMin) continue;
// Select the closest vertex
Double_t vsig=vrt->GetVsig();
if (nseedvrt < 0 || vsig > maxvsig) {
maxvsig = vsig; nseedvrt = ivrt;
}
}
return nseedvrt;
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::GetSeedVertexByNNMethod() {
Int_t nseedvrt = -1;
Int_t nvrt = GetVertexEntries();
if (nvrt < 2) return nseedvrt;
Double_t maxnn = -1.;
LCDVToplVRT* vrt=GetVertexAt(0);
Float_t nninp[3];
Float_t nnout;
for(Int_t ivrt =1; ivrt < nvrt ; ivrt++){
vrt = GetVertexAt(ivrt);
// Apply cuts to select seed vertex candidates
// Cut 1
Double_t decayL = vrt->GetDecayLength(ipv);
TVector3 vrtpos(vrt->GetVertexVector() - ipv);
TVector3 vrtmon(vrt->GetFourVector().Vect());
Double_t decayLerr = TMath::Sqrt(vrt->GetErrorDecayLength(ipv,iper));
Double_t vrtsig = decayL/decayLerr;
Double_t pdang = vrtpos.Angle(vrtmon);
decayL=TMath::Log10(decayL);
vrtsig=TMath::Log10(vrtsig);
pdang =TMath::Log10(TMath::Max(pdang,1e-8));
nninp[0]=TMath::Min(TMath::Max((decayL+3.0)/3.5,0.0),1.0);
nninp[1]=TMath::Min(TMath::Max((vrtsig+0.5)/3.5,0.0),1.0);
nninp[2]=TMath::Min(TMath::Max((pdang +3.5)/4.0,0.0),1.0);
LCDNNSeedVertex(nninp,&nnout,0);
if (nnout < m_SeedVrtNN) continue;
// Cut 2
Int_t ntrks = vrt->GetTrackEntries();
Double_t netCharge = vrt->GetVertexCharge();
if (ntrks == 2 && netCharge == 0.) {
Double_t reconMass = vrt->GetMass();
if ( TMath::Abs(reconMass - 0.49767) < m_SeedVrtK0Mass ) continue;
}
// Select the closest vertex
if (nseedvrt < 0 || nnout > maxnn) {
maxnn = nnout; nseedvrt = ivrt;
}
}
return nseedvrt;
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::CalcPTMass(Int_t iseedvrt) {
m_Rmass = 0.;
m_PtMass = 0.;
m_PtMassEr = 0.;
m_Masstag = 0.;
m_MissPt = 0.;
if (iseedvrt < 0 ) { return; }
TVector3 p3trk(0.0,0.0,0.0);
const Double_t mpi = 0.139567;
LCDVToplTRK* trk = 0;
LCDVToplVRT* seedvrt = GetVertexAt(iseedvrt);
Int_t nvtrk=seedvrt->GetTrackEntries();
Int_t ivtrk=0;
for (ivtrk=0 ; ivtrk < nvtrk ; ivtrk++) {
trk = seedvrt->GetTrackAt(ivtrk);
trk->SetTrkVal(1);
AddVrtTrack(trk);
}
m_p4Vertex=seedvrt->GetFourVector();
TVector3 xseedvrt(seedvrt->GetVertexVector());
TVector3 ipseed(xseedvrt - ipv);
Double_t dist = ipseed.Mag();
// Attatch tracks to the Seed Vertex
Float_t nninp[5];
Float_t nnout=0;
Int_t ntrks = GetTrackEntries();
Double_t l ;
Double_t lodi;
Double_t trdi;
Double_t anta;
Double_t ldist;
for (Int_t itrk = 0; itrk < ntrks ; itrk++){
trk = GetTrackAt(itrk);
if (trk->GetTrkVal() != 0) continue;
//Double_t l=trk->GetLminTrdiApprox(ipv,xseedvrt);
//Double_t l=trk->GetLminTrdi(ipv,xseedvrt);
//Double_t trdi = trk->GetTrdi(l,ipv,xseedvrt);
//Double_t lodi = trk->GetLodi(l,ipv,xseedvrt);
l = 0;
lodi = 0;
trdi = 0;
anta = 0;
trk->GetTrdiLodiAnta(l,ipv,xseedvrt,&trdi,&lodi,&anta);
ldist= lodi/dist;
if (m_UseNN) {
Double_t sbds =
TMath::Log10(trk->GetDistanceNormByError(l,ipv,iper));
nninp[0]=
TMath::Max(0.0,TMath::Min(1.0,trdi/0.3));
nninp[1]=
TMath::Max(0.0,TMath::Min(1.0,(lodi-0.3)/3.0));
nninp[2]=
TMath::Max(0.0,TMath::Min(1.0,(ldist-0.5)/3.0));
nninp[3]=
TMath::Max(0.0,TMath::Min(1.0,anta/TMath::Pi()));
nninp[4]=
TMath::Max(0.0,TMath::Min(1.0,sbds/3.0));
LCDNNTrackAttach(nninp,&nnout,0);
if (nnout < m_TrackAttNN) continue;
} else {
if (trdi > m_TrdiMax ) continue;
if (lodi < m_LodiMin ) continue;
if (ldist < m_LodrMin ) continue;
if (ldist > m_LodrMax ) continue;
if (anta > m_AntaMax ) continue;
if (trdi > lodi ) continue;
}
AddVrtTrack(trk);
trk->SetTrkVal(2);
p3trk = trk->GetMomentumVector(l); // Get momentum at Seed Vertex.
m_p4Vertex += TLorentzVector(p3trk,TMath::Sqrt(p3trk.Mag2() + mpi*mpi));
if (dblv > 9) {
cout<<" "<<" "
<<"itrk="<<itrk<<" "
<<"trdi="<<trdi<<" "
<<"lodi="<<lodi<<" "
<<"lovd="<<ldist<<" "
<<endl;
}
}
// Raw Mass
m_Rmass = m_p4Vertex.Mag();
// Momentum vector of reconstructed Vertex
TVector3 p3Vertex=m_p4Vertex.Vect();
Double_t pt = p3Vertex.Perp(ipseed);
// Missing Pt
m_MissPt = pt;
// Pt corrected mass
m_PtMass = TMath::Sqrt(m_Rmass*m_Rmass + pt*pt) + TMath::Abs(pt);
//Find delta phi and delta lambda.
Double_t sedxy = ipseed.Pt();
Double_t pxy = p3Vertex.Pt();
Double_t dPhi = p3Vertex.DeltaPhi(ipseed);
Double_t dLam = TMath::ATan((p3Vertex.Z()/pxy - ipseed.Z()/sedxy)/
(1.+p3Vertex.Z()*ipseed.Z()/(pxy*sedxy)));
//Rotate vertex error matrix
TMatrixD evrt(3,3); seedvrt->GetErrorMatrixVertex(evrt);
//TVector3 pvn = p3Vertex;
TVector3 pvn = ipseed;
Double_t pvnphi= pvn.Phi();
Double_t csphi = TMath::Cos(pvnphi);
Double_t snphi = TMath::Sin(pvnphi);
Double_t csthe = pvn.CosTheta();
Double_t snthe = TMath::Sqrt(1.0 - csthe*csthe);
TMatrixD dxpdxT(3,3);
dxpdxT(0,0)=+snphi; dxpdxT(0,1)=-csphi*csthe; dxpdxT(0,2)= csphi*snthe;
dxpdxT(1,0)=-csphi; dxpdxT(1,1)=-snphi*csthe; dxpdxT(1,2)= snphi*snthe;
dxpdxT(2,0)=+0.0 ; dxpdxT(2,1)= snthe; dxpdxT(2,2)= csthe;
//dxpdxT(0,0)=+snphi; dxpdxT(1,0)=-csphi*csthe; dxpdxT(2,0)= csphi*snthe;
//dxpdxT(0,1)=-csphi; dxpdxT(1,1)=-snphi*csthe; dxpdxT(2,1)= snphi*snthe;
//dxpdxT(0,2)=+0.0 ; dxpdxT(1,2)= snthe; dxpdxT(2,2)= csthe;
TMatrixD tmp0(evrt,TMatrixD::kMult,dxpdxT);
TMatrixD verrot(tmp0,TMatrixD::kTransposeMult,dxpdxT);
//Take errors into accound
//Allow sigma IP error to make conservative Min missing P
//Also allow seed vertex error n sigma Min missing pt
Double_t dPhiM = TMath::Abs(dPhi) -
TMath::Sqrt(m_PTIP[0]*m_PTIP[0] + m_NVsigma*m_NVsigma*verrot(0,0)
+ m_NVsigma*m_NVsigma*verrot(2,2)
*TMath::Sin(dPhi)*TMath::Sin(dPhi) )/sedxy;
if (dPhiM < 0.) dPhiM = 0.;
Double_t dLamM = TMath::Abs(dLam) -
TMath::Sqrt(m_PTIP[1]*m_PTIP[1]*sedxy*sedxy/(dist*dist)
+ m_NVsigma*m_NVsigma*verrot(1,1)
+ m_NVsigma*m_NVsigma
*verrot(2,2)*TMath::Sin(dLam)*TMath::Sin(dLam))/dist;
if (dLamM < 0.) dLamM = 0.;
pt = TMath::Sqrt(TMath::Power(p3Vertex.Mag()*TMath::Sin(dLamM),2)
+TMath::Power(pxy*TMath::Sin(dPhiM),2));
// Ptcorrected mass with IP errors
m_PtMassEr = TMath::Sqrt(m_Rmass*m_Rmass + pt*pt) + TMath::Abs(pt);
// Pt corrected IP err with 2*raw mass protection
if (m_PtMassEr > 2.*m_Rmass) {
m_Masstag = 2.*m_Rmass;
} else {
m_Masstag = m_PtMassEr;
}
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::CalcPTMass() {
// Get Seed Vertex
Int_t iseedvrt = GetSeedVertex();
CalcPTMass(iseedvrt);
}
//__________________________________________________________________________
//__________________________________________________________________________
Double_t LCDVTopl::GetHFTagNN() {
// Get Seed Vertex
Int_t iseedvrt = GetSeedVertex();
return GetHFTagNN(iseedvrt);
}
//__________________________________________________________________________
//__________________________________________________________________________
Double_t LCDVTopl::GetHFTagNN(Int_t iseed) {
if (iseed < 0 ) { return -1.0; } //no seed vertex...
if (m_Masstag == 0.0) {
CalcPTMass(iseed);
}
LCDVToplVRT* seedvrt = GetVertexAt(iseed);
TVector3 ipseed(seedvrt->GetVertexVector()-ipv);
Double_t dist=ipseed.Mag();
Double_t mvrt=GetMassTag();
Double_t pvrt=15*mvrt-GetP4Vertex().Rho();
Double_t dvrt=TMath::Log10(dist);
Double_t tvrt=GetVrtTrackEntries();
Double_t nvrt=GetVertexEntriesExceptV0();
Float_t nninp[5];
nninp[0]=TMath::Max(0.0,TMath::Min(1.0,mvrt/ 5.0));
nninp[1]=TMath::Max(0.0,TMath::Min(1.0,(pvrt+20.0)/100.0));
nninp[2]=TMath::Max(0.0,TMath::Min(1.0,(dvrt+2.0)/2.5));
nninp[3]=TMath::Max(0.0,TMath::Min(1.0,tvrt/10.0));
nninp[4]=TMath::Max(0.0,TMath::Min(1.0,(nvrt-1.0)/3.0));
Float_t nnout;
LCDNNHFSelect(nninp,&nnout,0);
return nnout;
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::GetNsig(Double_t nsigma, TVector3 axis, TVector3 xfrom,
TMatrixD& exfrom) {
LCDVToplTRK* trk=0;
Int_t itrk=0;
Int_t ntrk=GetTrackEntries();
Int_t nsig=0;
Double_t l=0;
for (itrk=0 ; itrk < ntrk ; itrk++) {
trk=GetTrackAt(itrk);
if (trk->GetV0flag()) continue;
l=trk->CalcPOCA3(xfrom);
if (trk->GetSignedDistanceNormByError(l,xfrom,axis,exfrom) > nsigma) {
nsig++;
}
}
return nsig;
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::GetNsig(Double_t nsigma) {
return GetNsig(nsigma,pja,ipv,iper);
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::GetNsig(Double_t nsigma, TVector3 axis, TVector3 xfrom,
TMatrixD& exfrom,TObjArray* nsig_list) {
LCDVToplTRK* trk=0;
Int_t itrk=0;
Int_t ntrk=GetTrackEntries();
Int_t nsig=0;
Double_t l=0;
for (itrk=0 ; itrk < ntrk ; itrk++) {
trk=GetTrackAt(itrk);
if (trk->GetV0flag()) continue;
l=trk->CalcPOCA3(xfrom);
if (trk->GetSignedDistanceNormByError(l,xfrom,axis,exfrom) > nsigma) {
nsig_list->Add((TObject*)trk);
nsig++;
}
}
return nsig;
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::GetNsig(Double_t nsigma,TObjArray* nsig_list) {
return GetNsig(nsigma,pja,ipv,iper,nsig_list);
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::GetNsig(Double_t nsigma, TVector3 axis, TVector3 xfrom,
TMatrixD& exfrom,TObjArray* nsig_list,
Double_t rnsigcut, Double_t znsigcut) {
LCDVToplTRK* trk=0;
Int_t itrk=0;
Int_t ntrk=GetTrackEntries();
Int_t nsig=0;
Double_t l=0;
TVector3 ip2ca;
Double_t rnsig;
Double_t znsig;
for (itrk=0 ; itrk < ntrk ; itrk++) {
trk=GetTrackAt(itrk);
if (trk->GetV0flag()) continue;
//r & z check...
ip2ca=trk->GetPositionVector(0.0)-xfrom;
rnsig=ip2ca.Perp();
znsig=TMath::Abs(ip2ca.Z());
if (rnsig > rnsigcut) continue;
if (znsig > znsigcut) continue;
l=trk->CalcPOCA3(xfrom);
if (trk->GetSignedDistanceNormByError(l,xfrom,axis,exfrom) > nsigma) {
nsig_list->Add((TObject*)trk);
nsig++;
}
}
return nsig;
}
//__________________________________________________________________________
//__________________________________________________________________________
Int_t LCDVTopl::GetNsig(Double_t nsigma,TObjArray* nsig_list,
Double_t rnsigcut, Double_t znsigcut) {
return GetNsig(nsigma,pja,ipv,iper,nsig_list,rnsigcut,znsigcut);
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::CheckV0Vertex() {
Int_t ntrk=GetTrackEntries();
LCDVToplTRK* trki=0;
LCDVToplTRK* trkj=0;
Int_t itrk=0;
Int_t jtrk=0;
TVector3 xvinit;
LCDVToplVRT vrt;
for (itrk=1 ; itrk < ntrk ; itrk++) {
trki=GetTrackAt(itrk);
for (jtrk=0 ; jtrk < itrk ;jtrk++) {
trkj=GetTrackAt(jtrk);
//track-track combination
xvinit=GetInitialVertexT2T(itrk,jtrk);
vrt.Clean();
vrt.SetVertex(xvinit);
vrt.AddTrack(trki);
vrt.AddTrack(trkj);
//V0 check...
if (vrt.GetChi2() < 16.0) {
vrt.CalcFourVector();
if (CheckKs0Vertex(&vrt)) {
trki->SetV0flag(1);
trkj->SetV0flag(1);
}
if (CheckLambda0Vertex(&vrt)) {
trki->SetV0flag(2);
trkj->SetV0flag(2);
}
}
}
}
}
//__________________________________________________________________________
//__________________________________________________________________________
Bool_t LCDVTopl::CheckKs0Vertex(LCDVToplVRT* vrt) {
Bool_t fk0s=kFALSE;
Int_t ntrks = vrt->GetTrackEntries();
Double_t netCharge = vrt->GetVertexCharge();
if (ntrks == 2 && netCharge == 0.) {
Double_t reconMass = vrt->GetMass();
if ( TMath::Abs(reconMass - 0.49767) < m_SeedVrtK0Mass ) fk0s=kTRUE;
}
return fk0s;
}
//__________________________________________________________________________
//__________________________________________________________________________
Bool_t LCDVTopl::CheckLambda0Vertex(LCDVToplVRT* vrt) {
Bool_t flambda=kFALSE;
Double_t mpi = 0.13957;
Double_t mp = 0.93827;
Int_t ntrks = vrt->GetTrackEntries();
Double_t netCharge = vrt->GetVertexCharge();
if (ntrks == 2 && netCharge == 0.) {
TVector3 p3trk1=vrt->GetTrackMomentumVector(0);
TVector3 p3trk2=vrt->GetTrackMomentumVector(1);
TLorentzVector p4trk1a(p3trk1,TMath::Sqrt(p3trk1.Mag2()+mpi*mpi));
TLorentzVector p4trk2a(p3trk2,TMath::Sqrt(p3trk1.Mag2()+mp *mp ));
TLorentzVector p4vrta=p4trk1a+p4trk2a;
if ( TMath::Abs(p4vrta.M() - 1.1157) < m_SeedVrtLambda0Mass ) {
flambda=kTRUE;
} else {
TLorentzVector p4trk1b(p3trk1,TMath::Sqrt(p3trk1.Mag2()+mp *mp ));
TLorentzVector p4trk2b(p3trk2,TMath::Sqrt(p3trk1.Mag2()+mpi*mpi));
TLorentzVector p4vrtb=p4trk1b+p4trk2b;
if ( TMath::Abs(p4vrtb.M() - 1.1157) < m_SeedVrtLambda0Mass )
flambda=kTRUE;
}
}
return flambda;
}
//__________________________________________________________________________
//__________________________________________________________________________
void LCDVTopl::PrintOutForDebug() {
if (dblv > 19) {
Int_t nvrt=GetVertexEntries();
LCDVToplVRT* vrt;
Int_t nvtrk;
for (Int_t ivrt=0 ; ivrt < nvrt ; ivrt++) {
vrt=GetVertexAt(ivrt);
nvtrk=vrt->GetTrackEntries();
cout<<" "
<<"ivrt="<<ivrt<<" "
<<"vrt="<<vrt<<" "
<<"ntrk="<<nvtrk<<" "
<<"decayL="<<vrt->GetDecayLength(ipv)<<" "
<<"chi2="<<vrt->GetChi2()<<" "
<<"Vsig="<<vrt->GetVsig()<<" "
<<endl;
if (dblv > 24) {
cout<<" "
<<"vpos="
<<vrt->GetVertexX()<<" "
<<vrt->GetVertexY()<<" "
<<vrt->GetVertexZ()<<" "
<<endl;
}
if (dblv > 29) {
for (Int_t itrk=0 ; itrk < nvtrk ; itrk++)
cout<<" "
<<"itrk="<<itrk<<" "
<<"trk="<<vrt->GetTrackAt(itrk)<<" "
<<endl;
}
}
if (dblv > 34) {
LCDVToplTRK* trk;
for (Int_t itrk=0 ; itrk < GetTrackEntries() ; itrk++) {
trk=(LCDVToplTRK*)GetTrackAt(itrk);
cout<<" "
<<"itrk="<<itrk<<" "
<<"trk="<<trk<<" "
<<endl;
for (Int_t itvrt=0 ; itvrt < trk->GetVertexEntries() ; itvrt++) {
vrt=(LCDVToplVRT*)(trk->GetVertexAt(itvrt));
cout<<" "
<<" "<<" "
<<"itvrt="<<itvrt<<" "
<<"vrt="<<vrt<<" "
<<endl;
}
}
}
}
}
//__________________________________________________________________________
ROOT 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.