#include <iostream>
#include "hrktrackBF2.h"
#include "hrktrackB.h"
#include "hades.h"
#include "hruntimedb.h"
#include "hevent.h"
#include "hdetector.h"
#include "hiterator.h"
#include "hcategory.h"
#include "hmdctrkcand.h"
#include "hmdcgetcontainers.h"
#include "hrungekutta.h"
#include "hmetamatch2.h"
#include "htofcluster.h"
#include "htofclustersim.h"
#include "hgeomvolume.h"
#include "hspectrometer.h"
#include "hmdcdetector.h"
#include "hmdcgeompar.h"
#include "htofgeompar.h"
#include "hshowergeometry.h"
#include "hemcgeompar.h"
#include "hspecgeompar.h"
#include "tofdef.h"
#include "hmdcseg.h"
#include "hmdchit.h"
#include "showerdef.h"
#include "emcdef.h"
#include "hmdctrackgfieldpar.h"
#include "hmatrixcategory.h"
#include "hgeomtransform.h"
#include "hmdctrackddef.h"
#include "hmdctrackgdef.h"
#include "hmdctrackgspline.h"
#include "hgeomvector.h"
#include "hsplinetrack.h"
#include "hmagnetpar.h"
#include "hshowerhitsim.h"
#include "hemcclustersim.h"
#include "hmdcsizescells.h"
#include "hrpccluster.h"
#include "rpcdef.h"
#include "hrpcgeompar.h"
#include "hmetamatchpar.h"
using namespace std;
ClassImp(HRKTrackBF2)
HRKTrackBF2::HRKTrackBF2() {
mode=0;
clear();
}
HRKTrackBF2::HRKTrackBF2(const Text_t *name, Short_t m) : HReconstructor(name,name) {
mode=m;
if(mode==1) {
Error("HRKTrackBF2","mode=%i is not supported any more!",mode);
mode = 2;
}
clear();
}
void HRKTrackBF2::clear()
{
fMatchPar = NULL;
pMagnet = NULL;
field = NULL;
fieldFactor = -1000.;
fSpecGeomPar = NULL;
fGetCont = NULL;
fTofGeometry = NULL;
fRpcGeometry = NULL;
fShowerGeometry = NULL;
fEmcGeometry = NULL;
fCatMetaMatch = NULL;
fMetaMatchIter = NULL;
fCatMdcTrkCand = NULL;
fCatRpcCluster = NULL;
fCatMdcSeg = NULL;
fCatMdcHit = NULL;
fSplineTrack = NULL;
fCatKine = NULL;
fCatShower = NULL;
fCatEmc = NULL;
fCatTof = NULL;
fCatTofCluster = NULL;
fCatRKTrack = NULL;
pMetaMatch = NULL;
pMdcTrkCand = NULL;
pSplineTrack = NULL;
pTofCluster = NULL;
pShowerHit = NULL;
pEmcCluster = NULL;
pRungeKutta = NULL;
pMdcSeg1 = NULL;
pMdcSeg2 = NULL;
qualityRpc = -1.;
qualityShower = -1.;
qualityTof = -1.;
for(Int_t i = 0; i < 3; ++i)
{
pTofHit[i] = NULL;
}
for(Int_t m = 0;m < 4; m++)
{
for(Int_t s = 0;s < 6; s++)
{
mdcInstalled[m][s]=kFALSE;
}
}
}
HRKTrackBF2::~HRKTrackBF2() {
if (pRungeKutta) {
delete pRungeKutta;
pRungeKutta=0;
}
HMdcSizesCells::deleteCont();
}
Bool_t HRKTrackBF2::init(){
if (gHades) {
HRuntimeDb *rtdb=gHades->getRuntimeDb();
HSpectrometer *spec=gHades->getSetup();
HEvent *event=gHades->getCurrentEvent();
if (!event) return kFALSE;
field=(HMdcTrackGFieldPar*)(rtdb->getContainer("MdcTrackGFieldPar"));
fSpecGeomPar=(HSpecGeomPar*)(rtdb->getContainer("SpecGeomPar"));
pMagnet=(HMagnetPar*)( rtdb->getContainer("MagnetPar"));
fGetCont=HMdcGetContainers::getObject();
pMSizesCells = HMdcSizesCells::getObject();
if (spec->getDetector("Tof")) {
fTofGeometry = (HTofGeomPar *)rtdb->getContainer("TofGeomPar");
}
if (spec->getDetector("Shower")) {
fShowerGeometry = (HShowerGeometry *)rtdb->getContainer("ShowerGeometry");
}
if (spec->getDetector("Emc")) {
fEmcGeometry = (HEmcGeomPar *)rtdb->getContainer("EmcGeomPar");
}
if (spec->getDetector("Rpc"))
{
fRpcGeometry = (HRpcGeomPar*)rtdb->getContainer("RpcGeomPar");
}
fMatchPar=(HMetaMatchPar*)rtdb->getContainer("MetaMatchPar");
fCatMetaMatch=event->getCategory(catMetaMatch);
if (!fCatMetaMatch) return kFALSE;
fMetaMatchIter=(HIterator*)fCatMetaMatch->MakeIterator();
if (!fMetaMatchIter) return kFALSE;
fCatMdcTrkCand=event->getCategory(catMdcTrkCand);
if (!fCatMdcTrkCand) return kFALSE;
fCatMdcSeg=event->getCategory(catMdcSeg);
if (!fCatMdcSeg) return kFALSE;
fCatMdcHit=event->getCategory(catMdcHit);
fCatRpcCluster = event->getCategory(catRpcCluster);
if(!fCatRpcCluster) Warning("init","NO catRpcCluster in input! \n");
if (mode==2) {
fSplineTrack=event->getCategory(catSplineTrack);
if (!fSplineTrack) {
Error("init",
"Spline is used as initial momentum guess, but was not initialized before Runge-Kutta");
return kFALSE;
}
}
fCatTof=event->getCategory(catTofHit);
if (!fCatTof) Warning("init","No catTofHit in input!");;
fCatTofCluster=event->getCategory(catTofCluster);
if (!fCatTofCluster) Warning("init","No catTofCluster in input!");
fCatShower = event->getCategory(catShowerHit);
fCatEmc = event->getCategory(catEmcCluster);
if(fCatShower == NULL && fCatEmc == NULL) {
Warning("init","NO catShowerHit and catEmcCluster in input!");
}
fCatRKTrack=event->getCategory(catRKTrackB);
if (!fCatRKTrack) {
Int_t size[2]={6,8000};
fCatRKTrack=new HMatrixCategory("HRKTrackB",2,size,0.5);
if (fCatRKTrack) event->addCategory(catRKTrackB,fCatRKTrack,"RKTrackB");
}
return kTRUE;
}
return kFALSE;
}
Bool_t HRKTrackBF2::reinit()
{
#warning will be changed: Anar
if (!pRungeKutta) {
pRungeKutta = new HRungeKutta();
}
if (pMagnet->hasChanged() || fieldFactor ==-1000) {
fieldFactor=pMagnet->getScalingFactor();
pRungeKutta->setField(field->getPointer(),fieldFactor);
}
if (!pMSizesCells->initContainer()) {
Error("reinit()","HMdcSizesCells is not initialzed!");
return kFALSE;
}
HGeomTransform gtrMod2Sec;
for(Int_t is=0; is<6; is++) {
HMdcSizesCellsSec& fSCSec = (*pMSizesCells)[is];
if(&fSCSec == 0) continue;
secTrans[is]=*(fSCSec.getLabTrans());
for(Int_t im=0; im<4; im++) {
HMdcSizesCellsMod& fSCMod=fSCSec[im];
if (&fSCMod) {
gtrMod2Sec = *(fSCMod.getSecTrans());
pRungeKutta->setMdcPosition(is,im,gtrMod2Sec);
mdcInstalled[im][is]=kTRUE;
}
}
}
if ( !fMatchPar ) return kFALSE;
for(Int_t is=0; is<6; is++) {
if (fRpcGeometry) {
HModGeomPar *module = fRpcGeometry -> getModule(is,0);
HGeomTransform modTrans( module -> getLabTransform());
modTrans.transTo(secTrans[is]);
rpcSM[is] = modTrans;
HGeomVector r0_mod(0.0, 0.0, 0.0);
centerRpc[is] = modTrans.transFrom(r0_mod);
HGeomVector rz_mod(0.0, 0.0, 1.0);
normVecRpc[is] = modTrans.transFrom(rz_mod) - centerRpc[is];
}
if (fShowerGeometry) {
HGeomTransform modTrans(fShowerGeometry->getTransform(is));
modTrans.transTo(secTrans[is]);
showerSM[is] = modTrans;
HGeomVector r0_mod(0.0, 0.0, 0.);
HGeomVector rz_mod(0.0, 0.0, 1.);
normVecShower[is] = modTrans.transFrom(rz_mod) - modTrans.transFrom(r0_mod);
}
if (fEmcGeometry) {
HModGeomPar *pmodgeom = fEmcGeometry->getModule(is);
HGeomTransform modTrans(pmodgeom->getLabTransform());
modTrans.transTo(secTrans[is]);
emcSM[is] = modTrans;
HGeomVector r0_mod(0.0, 0.0, 0.);
HGeomVector rz_mod(0.0, 0.0, 1.);
normVecEmc[is] = modTrans.transFrom(rz_mod) - modTrans.transFrom(r0_mod);
}
if (fTofGeometry) {
for(Int_t im=0; im<8; im++) {
HModGeomPar *module = fTofGeometry->getModule(is,im);
HGeomTransform modTrans(module->getLabTransform());
modTrans.transTo(secTrans[is]);
tofSM[is][im] = modTrans;
HGeomVector r0_mod(0.0, 0.0, 0.0);
HGeomVector rz_mod(0.0, 0.0, 1.0);
normVecTof[is][im] = modTrans.transFrom(rz_mod) - modTrans.transFrom(r0_mod);
}
}
setMatchingParams(is);
}
return kTRUE;
}
void HRKTrackBF2::setMatchingParams(Int_t s)
{
if(s>5||s<0) {
Error("setMatchingParams()","s = %i out of range!",s);
return;
}
sigma2MdcInRpcX[s] = fMatchPar -> getRpcSigmaXMdc(s)*fMatchPar -> getRpcSigmaXMdc(s);
sigma2MdcInRpcY[s] = fMatchPar -> getRpcSigmaYMdc(s)*fMatchPar -> getRpcSigmaYMdc(s);
sRpcX[s] = fMatchPar -> getRpcSigmaXOffset(s);
sRpcY[s] = fMatchPar -> getRpcSigmaYOffset(s);
quality2RPCCut[s] = fMatchPar -> getRpcQualityCut(s)*fMatchPar -> getRpcQualityCut(s);
sigma2MdcInShrX[s] = fMatchPar -> getShowerSigmaXMdc(s)*fMatchPar -> getShowerSigmaXMdc(s);
sigma2MdcInShrY[s] = fMatchPar -> getShowerSigmaYMdc(s)*fMatchPar -> getShowerSigmaYMdc(s);
sShowerX[s] = fMatchPar -> getShowerSigmaXOffset(s);
sShowerY[s] = fMatchPar -> getShowerSigmaYOffset(s);
quality2SHOWERCut[s] = fMatchPar -> getShowerQualityCut(s)*fMatchPar -> getShowerQualityCut(s);
sigma2MdcInEmcX[s] = fMatchPar -> getEmcSigmaXMdc(s)*fMatchPar -> getEmcSigmaXMdc(s);
sigma2MdcInEmcY[s] = fMatchPar -> getEmcSigmaYMdc(s)*fMatchPar -> getEmcSigmaYMdc(s);
sEmcX[s] = fMatchPar -> getEmcSigmaXOffset(s);
sEmcY[s] = fMatchPar -> getEmcSigmaYOffset(s);
quality2EMCCut[s] = fMatchPar -> getEmcQualityCut(s)*fMatchPar -> getEmcQualityCut(s);
sigma2TofX[s] = fMatchPar -> getTofSigmaX(s)*fMatchPar -> getTofSigmaX(s);
sigma2TofY[s] = fMatchPar -> getTofSigmaY(s)*fMatchPar -> getTofSigmaY(s);
sTofX[s] = fMatchPar -> getTofSigmaXOffset(s);
sTofY[s] = fMatchPar -> getTofSigmaYOffset(s);
quality2TOFCut[s] = fMatchPar -> getTofQualityCut(s)*fMatchPar -> getTofQualityCut(s);
}
Int_t HRKTrackBF2::execute()
{
success = kTRUE;
fMetaMatchIter->Reset();
while ((pMetaMatch=(HMetaMatch2*)(fMetaMatchIter->Next()))!=0)
{
sector = pMetaMatch->getSector();
pMdcTrkCand = (HMdcTrkCand*)fCatMdcTrkCand -> getObject(pMetaMatch->getTrkCandInd());
if(!pMdcTrkCand) { continue; }
pMdcSeg1 = (HMdcSeg*)fCatMdcSeg -> getObject(pMdcTrkCand->getSeg1Ind());
pMdcSeg2 = (HMdcSeg*)fCatMdcSeg -> getObject(pMdcTrkCand->getSeg2Ind());
if (!pMdcSeg1 || !pMdcSeg2) { continue; }
sectorloc.set(1,sector);
Double_t u1pos[3], u1dir[3], u2pos[3], u2dir[3];
Bool_t doFakeHit=kFALSE;
Int_t nHitsCheck[2]={0};
if(pMdcSeg1->getHitInd(0)!=-1) nHitsCheck[0]++;
if(pMdcSeg1->getHitInd(1)!=-1) nHitsCheck[0]++;
if(pMdcSeg2->getHitInd(0)!=-1) nHitsCheck[1]++;
if(pMdcSeg2->getHitInd(1)!=-1) nHitsCheck[1]++;
if(mdcInstalled[0][pMdcSeg1->getSec()] &&
mdcInstalled[1][pMdcSeg1->getSec()] &&
nHitsCheck[0]==1 && nHitsCheck[1]==1)
{
doFakeHit=kTRUE;
}
Int_t nHits=calcPosDirFromSegment(pMdcSeg1,0,&u1pos[0],&u1dir[0],doFakeHit);
if(nHits==1 &&
mdcInstalled[2][pMdcSeg2->getSec()] &&
mdcInstalled[3][pMdcSeg2->getSec()] &&
nHitsCheck[0]==1 && nHitsCheck[1]==1)
{
doFakeHit=kTRUE;
}
else {
doFakeHit=kFALSE;
}
nHits +=calcPosDirFromSegment(pMdcSeg2,1,&u2pos[0],&u2dir[0],doFakeHit);
if (nHits<3) {
Warning("execute()","Less than 3 points for RK detected. This should never happen!");
continue;
}
Int_t splineTrackIndex=-1;
momentumGuess=-1000000.0;
pRK=-1000000.0;
qRK=0;
if (mode==2) {
splineTrackIndex = pMetaMatch->getSplineInd();
if (splineTrackIndex>=0) {
pSplineTrack = (HSplineTrack*)fSplineTrack->getObject(splineTrackIndex);
pRK = pSplineTrack -> getP();
qRK = pSplineTrack -> getPolarity();
} else {
Warning("execute()","SplineTrack did not provide information for RK");
}
}
if (pRK>50.) {
momentumGuess=qRK*pRK;
} else if (pRK!=-1000000.0) {
momentumGuess=qRK*50.;
}
if (nHits == 4) {
success = pRungeKutta->fit4Hits(u1pos,u1dir,u2pos,u2dir,multSig,sector,momentumGuess);
} else {
if (momentumGuess != -1000000.0) {
success = pRungeKutta->fit3Hits(u1pos,u1dir,u2pos,u2dir,multSig,sector,momentumGuess);
} else {
Warning("execute()",
"HRKTrackF: MDC123 mode, but user did not provide momentum guess - RK failed!");
success = kFALSE;
}
}
if(fieldFactor==0.){
success = kTRUE;
doMassStuff();
fillData(pMdcSeg1,pMdcSeg2,pSplineTrack,indexRK);
pMetaMatch -> setRungeKuttaInd((Short_t)indexRK);
success = kFALSE;
continue;
}
if (!success) {
fillData(pMdcSeg1,pMdcSeg2,pSplineTrack,indexRK);
pMetaMatch -> setRungeKuttaInd((Short_t)indexRK);
continue;
}
chiqRK=pRungeKutta->getChi2();
pRK=pRungeKutta->getPfit();
if (pRK>0) {
qRK=+1;
} else {
pRK*=-1.;
qRK=-1;
}
if (pRK>10000.) pRK=10000.0;
pRungeKutta->traceToVertex(pMSizesCells);
if ( !doMassStuff() )
{
tof = -1.0;
trackLength = -1.0;
beta = -1.0;
mass2 = -1.0;
metaeloss = -1.0;
RKxyzMETA[0] = -1.0;
RKxyzMETA[1] = -1.0;
RKxyzMETA[2] = -1.0;
pRungeKutta -> traceToMETA(centerRpc[sector],normVecRpc[sector]);
fillData(pMdcSeg1,pMdcSeg2,pSplineTrack,indexRK);
pMetaMatch -> setRungeKuttaInd((Short_t)indexRK);
}
}
return 0;
}
Bool_t HRKTrackBF2::doMassStuff()
{
Bool_t isMatched = kFALSE;
if( pMetaMatch -> getSystem() == -1) return kFALSE;
if( pMetaMatch -> getNRpcClusters() )
{
isMatched = kTRUE;
matchWithRpc();
}
if( fShowerGeometry!=NULL && pMetaMatch -> getNShrHits() )
{
isMatched = kTRUE;
matchWithShower();
}
else if(fEmcGeometry!=NULL && pMetaMatch -> getNEmcClusters() )
{
isMatched = kTRUE;
matchWithEmc();
}
if ( pMetaMatch -> getNTofHits() )
{
isMatched = kTRUE;
matchWithTof();
}
if(!isMatched) return kFALSE;
return kTRUE;
}
void HRKTrackBF2::matchWithRpc()
{
metaNormVec = normVecRpc[sector];
transMetaSM = rpcSM[sector];
for(UChar_t n = 0; n < pMetaMatch -> getNRpcClusters(); ++n )
{
indRpc = pMetaMatch -> getRpcClstInd(n);
if( indRpc < 0 )
{
Warning("matchWithRpc","Rpc index is < 0, DISASTER!");
continue;
}
pRpc = (HRpcCluster*)fCatRpcCluster -> getObject(indRpc);
if(!pRpc)
{
Warning("matchWithRpc","Pointer to Rpc is NULL, DISASTER!");
continue;
}
tof = pRpc -> getTof();
metaeloss = pRpc -> getCharge();
xTof = pRpc -> getXSec();
yTof = pRpc -> getYSec();
zTof = pRpc -> getZSec();
zMod = 0;
pointMeta.setXYZ(xTof,yTof,zTof);
qualityRpc = pMetaMatch -> getRpcClstQuality(n);
RKxyzMETA[0] = pMetaMatch -> getRpcClstDX(n);
RKxyzMETA[1] = pMetaMatch -> getRpcClstDY(n);
RKxyzMETA[2] = 0.;
calcBeta(zMod,0);
HRKTrackB* pRK = fillData(pMdcSeg1,pMdcSeg2,pSplineTrack,indexRK);
pRK->setRpcClustInd(indRpc);
pMetaMatch -> setRungeKuttaIndRpcClst(n, (Short_t)indexRK);
}
}
void HRKTrackBF2::matchWithShower()
{
if(!fCatShower) return;
metaNormVec = normVecShower[sector];
transMetaSM = showerSM[sector];
for(UChar_t n = 0; n < pMetaMatch -> getNShrHits(); ++n )
{
indShower = pMetaMatch -> getShowerHitInd(n);
if( indShower < 0 )
{
Warning("matchWithShower","Index of shower is < 0, DISASTER!");
continue;
}
pShowerHit = (HShowerHit*)fCatShower -> getObject(indShower);
if(!pShowerHit)
{
Warning("matchWithShower","Pointer to Shower is NULL, DISASTER!");
continue;
}
tof = -1.0;
beta = -1.0;
mass2 = -1.0;
metaeloss = -1.0;
pShowerHit -> getLabXYZ(&xTof,&yTof,&zTof);
pointMeta.setXYZ(xTof,yTof,zTof);
pointMeta = secTrans[sector].transTo(pointMeta);
zMod = pShowerHit-> getZ();
qualityShower = pMetaMatch -> getShowerHitQuality(n);
RKxyzMETA[0] = pMetaMatch -> getShowerHitDX(n);
RKxyzMETA[1] = pMetaMatch -> getShowerHitDY(n);
RKxyzMETA[2] = 0.;
calcBeta(zMod,1, kFALSE);
HRKTrackB* pRK = fillData(pMdcSeg1,pMdcSeg2,pSplineTrack,indexRK);
pRK->setShowerHitInd(indShower);
pMetaMatch -> setRungeKuttaIndShowerHit(n,(Short_t)indexRK);
}
}
void HRKTrackBF2::matchWithEmc()
{
if(!fCatEmc) return;
metaNormVec = normVecEmc[sector];
transMetaSM = emcSM[sector];
for(UChar_t n = 0; n < pMetaMatch -> getNEmcClusters(); ++n )
{
indEmc = pMetaMatch -> getEmcClusterInd(n);
if( indEmc < 0 )
{
Warning("matchWithEmc","Index of Emc is < 0, DISASTER!");
continue;
}
pEmcCluster = (HEmcCluster*)fCatEmc -> getObject(indEmc);
if(!pEmcCluster)
{
Warning("matchWithEmc","Pointer to Emc is NULL, DISASTER!");
continue;
}
tof = pEmcCluster->getTime();
beta = -1.0;
mass2 = -1.0;
metaeloss = -1.0;
pEmcCluster -> getXYZLab(xTof,yTof,zTof);
pointMeta.setXYZ(xTof,yTof,zTof);
pointMeta = secTrans[sector].transTo(pointMeta);
zMod = 0.;
qualityEmc = pMetaMatch -> getEmcClusterQuality(n);
RKxyzMETA[0] = pMetaMatch -> getEmcClusterDX(n);
RKxyzMETA[1] = pMetaMatch -> getEmcClusterDY(n);
RKxyzMETA[2] = 0.;
calcBeta(zMod,3, kTRUE);
HRKTrackB* pRK = fillData(pMdcSeg1,pMdcSeg2,pSplineTrack,indexRK);
pRK->setShowerHitInd(indEmc);
pMetaMatch -> setRungeKuttaIndEmcCluster(n,(Short_t)indexRK);
}
}
void HRKTrackBF2::calcBeta(Float_t zMod, Int_t mode, Bool_t option)
{
if(fieldFactor==0.) {
HMdcSizesCellsSec& fSCSec=(*pMSizesCells)[sector];
const HGeomVector& target=fSCSec.getTargetMiddlePoint();
HGeomVector distance = pointMeta - target;
trackLength = distance.length();
beta = 1.0e6*trackLength/tof/TMath::C();
mass2 = 1.0e6;
chiqRK = 1.;
return;
}
pRungeKutta -> traceToMETA(pointMeta, metaNormVec);
HGeomVector localMeta1(pRungeKutta -> getXtrackOnMETA(),
pRungeKutta -> getYtrackOnMETA(),
pRungeKutta -> getZtrackOnMETA());
HGeomVector localMeta2(pointMeta);
localMeta1 = transMetaSM.transTo(localMeta1);
localMeta2 = transMetaSM.transTo(localMeta2);
HGeomVector shiftz(0.,0.,zMod);
localMeta2 -= shiftz;
HGeomVector localMeta = localMeta2 - localMeta1;
RKxyzMETA[0] = localMeta.getX();
RKxyzMETA[1] = localMeta.getY();
RKxyzMETA[2] = localMeta.getZ();
qualityRpc = qualityShower = qualityTof = -1.;
if(mode == 0)
{
dXrms2 = pRpc -> getXRMS();
dYrms2 = pRpc -> getYRMS();
dX = localMeta.getX() - sRpcX[sector];
dY = localMeta.getY() - sRpcY[sector];
qualityRpc = getQuality(dX,dY,dXrms2*dXrms2 + sigma2MdcInRpcX[sector],dYrms2*dYrms2 + sigma2MdcInRpcY[sector]);
}
else if (mode == 1)
{
dXrms2 = pShowerHit -> getSigmaX();
dYrms2 = pShowerHit -> getSigmaY();
dX = localMeta.getX() - sShowerX[sector];
dY = localMeta.getY() - sShowerY[sector];
qualityShower = getQuality(dX, dY,dXrms2*dXrms2 + sigma2MdcInShrX[sector],dYrms2*dYrms2 + sigma2MdcInShrY[sector]);
}
else if(mode == 2)
{
dX = localMeta.getX() - sTofX[sector];
dY = localMeta.getY() - sTofY[sector];
qualityTof = getQuality(dX, dY, sigma2TofX[sector], sigma2TofY[sector]);
}
else if (mode == 3)
{
dXrms2 = pEmcCluster -> getSigmaXMod();
dYrms2 = pEmcCluster -> getSigmaYMod();
dX = localMeta.getX() - sEmcX[sector];
dY = localMeta.getY() - sEmcY[sector];
qualityEmc = getQuality(dX, dY,dXrms2*dXrms2 + sigma2MdcInEmcX[sector],dYrms2*dYrms2 + sigma2MdcInEmcY[sector]);
}
else
{
Warning("calcBeta", "This option does not exist");
}
trackLength = pRungeKutta->getTrackLength();
if(option)
{
beta = 1.0e6*trackLength/tof/TMath::C();
mass2 = pRK*pRK*(1-beta*beta)/beta/beta;
}
}
Float_t HRKTrackBF2::getQuality(Float_t dx, Float_t dy, Float_t s2x, Float_t s2y)
{
return sqrt(dx*dx/s2x + dy*dy/s2y);
}
void HRKTrackBF2::matchWithTof()
{
for(UChar_t n = 0; n < pMetaMatch -> getNTofHits(); ++n )
{
indTof[0] = pMetaMatch -> getTofHit1Ind(n);
indTof[1] = pMetaMatch -> getTofHit2Ind(n);
indTof[2] = pMetaMatch -> getTofClstInd(n);
for(Int_t i = 0; i < 3; ++i )
{
if( indTof[i] < 0 ) continue;
if (i == 2)
pTofHit[i] = (HTofHit*)fCatTofCluster->getObject(indTof[i]);
else
pTofHit[i] = (HTofHit*)fCatTof->getObject(indTof[i]);
if( !pTofHit[i] )
{
Warning("matchWithTof","Pointer to Tof is NULL, DISASTER!");
continue;
}
metaNormVec = normVecTof[sector][(Int_t)pTofHit[i] -> getModule()];
transMetaSM = tofSM[sector][(Int_t)pTofHit[i] -> getModule()];
tof = pTofHit[i] -> getTof();
metaeloss = pTofHit[i] -> getEdep();
pTofHit[i] -> getXYZLab(xTof, yTof, zTof);
pointMeta.setXYZ(xTof,yTof,zTof);
pointMeta = secTrans[sector].transTo(pointMeta);
zMod = 0;
if( i == 2){
qualityTof = pMetaMatch -> getTofClstQuality(n);
RKxyzMETA[0] = pMetaMatch -> getTofClstDX(n);
RKxyzMETA[1] = pMetaMatch -> getTofClstDY(n);
} else if (i == 1){
qualityTof = pMetaMatch -> getTofHit2Quality(n);
RKxyzMETA[0] = pMetaMatch -> getTofHit2DX(n);
RKxyzMETA[1] = pMetaMatch -> getTofHit2DY(n);
} else {
qualityTof = pMetaMatch -> getTofHit1Quality(n);
RKxyzMETA[0] = pMetaMatch -> getTofHit1DX(n);
RKxyzMETA[1] = pMetaMatch -> getTofHit1DY(n);
}
RKxyzMETA[2] = 0.;
calcBeta(zMod,2);
HRKTrackB* pRK = fillData(pMdcSeg1,pMdcSeg2,pSplineTrack,indexRK);
if( i == 0) { pRK->setTofHitInd (indTof[0]); pMetaMatch -> setRungeKuttaIndTofHit1(n, (Short_t)indexRK); }
else if( i == 1) { pRK->setTofHitInd (indTof[1]); pMetaMatch -> setRungeKuttaIndTofHit2(n, (Short_t)indexRK); }
else { pRK->setTofClustInd(indTof[2]); pMetaMatch -> setRungeKuttaIndTofClst(n, (Short_t)indexRK); }
}
}
}
Int_t HRKTrackBF2::calcPosDirFromSegment(HMdcSeg* pSeg,Int_t ioseg, Double_t* pos, Double_t* dir, Bool_t flag) {
Int_t index[2]={pSeg->getHitInd(0),pSeg->getHitInd(1)};
Int_t nHits=2;
if (ioseg==0) {
for (Int_t i=0;i<2;i++) {
if (index[i]==-1 && flag==kFALSE) {
multSig[2*i] = multSig[2*i+1] = -1.;
nHits--;
} else multSig[2*i] = multSig[2*i+1] = 1.;
}
} else {
for (Int_t i=0;i<2;i++) {
if (index[i]==-1) {
if (flag==kFALSE) {
multSig[4+2*i] = multSig[4+2*i+1] = -1.;
nHits--;
} else {
if(pSeg->getChi2()>=0) multSig[4+2*i] = multSig[4+2*i+1] = 1.;
else multSig[4+2*i] = multSig[4+2*i+1] = 5.;
}
} else if (pSeg->getChi2()>=0) {
multSig[4+2*i] = multSig[4+2*i+1] = 1.;
} else {
multSig[4+2*i] = multSig[4+2*i+1] = 5.;
}
}
}
if (nHits>0) {
Double_t hpi=TMath::Pi()/2.;
pos[0] = pSeg->getR()*cos(pSeg->getPhi()+hpi);
pos[1] = pSeg->getR()*sin(pSeg->getPhi()+hpi);
pos[2] = pSeg->getZ();
dir[0] = cos(pSeg->getPhi())*sin(pSeg->getTheta());
dir[1] = sin(pSeg->getPhi())*sin(pSeg->getTheta());
dir[2] = cos(pSeg->getTheta());
}
return nHits;
}
HRKTrackB* HRKTrackBF2::fillData(HMdcSeg* segment1,HMdcSeg* segment2,HSplineTrack* spline, Int_t &indexRK)
{
TObject* slot=0;
HRKTrackB* rkt=0;
slot=fCatRKTrack->getNewSlot(sectorloc,&indexRK);
if (!slot) {
Error("fillData","No slot available");
} else {
rkt=new(slot) HRKTrackB;
rkt->setSector(sector);
rkt->setZ( segment1->getZ() ,segment1->getErrZ() );
rkt->setR( segment1->getR() ,segment1->getErrR() );
rkt->setTheta(segment1->getTheta() ,segment1->getErrTheta() );
rkt->setPhi( segment1->getPhi() ,segment1->getErrPhi() );
rkt->setP(pRK,0.);
rkt->setPolarity(qRK);
rkt->setTof(tof);
rkt->setMetaEloss(metaeloss);
rkt->setTarDist(spline->getTarDist());
rkt->setIOMatching(spline->getIOMatching());
if (success) {
rkt->setChiq(chiqRK);
if(fieldFactor==0.) {
rkt->setZSeg1RK(segment1->getZ());
rkt->setRSeg1RK(segment1->getR());
rkt->setThetaSeg1RK(segment1->getTheta());
rkt->setPhiSeg1RK(segment1->getPhi());
rkt->setZSeg2RK(segment2->getZ());
rkt->setRSeg2RK(segment2->getR());
rkt->setThetaSeg2RK(segment2->getTheta());
rkt->setPhiSeg2RK(segment2->getPhi());
} else {
rkt->setZSeg1RK(pRungeKutta->getZSeg1());
rkt->setRSeg1RK(pRungeKutta->getRSeg1());
rkt->setThetaSeg1RK(pRungeKutta->getThetaSeg1());
rkt->setPhiSeg1RK(pRungeKutta->getPhiSeg1());
rkt->setZSeg2RK(pRungeKutta->getZSeg2());
rkt->setRSeg2RK(pRungeKutta->getRSeg2());
rkt->setThetaSeg2RK(pRungeKutta->getThetaSeg2());
rkt->setPhiSeg2RK(pRungeKutta->getPhiSeg2());
}
rkt->setTofDist(trackLength);
rkt->setBeta(beta);
rkt->setMass2(mass2, 0.);
rkt->setMETAdx(RKxyzMETA[0]);
rkt->setMETAdy(RKxyzMETA[1]);
rkt->setMETAdz(RKxyzMETA[2]);
rkt->setQualityRpc(qualityRpc);
if(fCatShower != NULL) rkt->setQualityShower(qualityShower);
else rkt->setQualityShower(qualityEmc);
rkt->setQualityTof(qualityTof);
} else {
rkt->setZSeg1RK(segment1->getZ());
rkt->setRSeg1RK(segment1->getR());
rkt->setThetaSeg1RK(segment1->getTheta());
rkt->setPhiSeg1RK(segment1->getPhi());
rkt->setZSeg2RK(segment2->getZ());
rkt->setRSeg2RK(segment2->getR());
rkt->setThetaSeg2RK(segment2->getTheta());
rkt->setPhiSeg2RK(segment2->getPhi());
rkt->setTofDist(spline->getTofDist());
rkt->setBeta( spline->getBeta());
rkt->setMass2( spline->getMass2(), 0.);
}
}
return rkt;
}