using namespace std;
#include "hmdctrackfitter.h"
#include "hdebug.h"
#include "hades.h"
#include "hevent.h"
#include "hmdctrackfitpar.h"
#include "hruntimedb.h"
#include "hmdcsizescells.h"
#include "hgeomvector.h"
#include "hmdcgetcontainers.h"
#include "hmdclistcells.h"
#include "hmdcwirefitsim.h"
#include "hmdcclusfitsim.h"
#include "hmdcclussim.h"
#include "hmdctrackdset.h"
#include "hmdcdigitpar.h"
ClassImp(HMdcTrackFitInOut)
ClassImp(HMdcTrackFitter)
HMdcTrackFitInOut::HMdcTrackFitInOut(void) {
fSizesCells=0;
setDefaultFitParam();
useRndbPar = kTRUE;
}
Bool_t HMdcTrackFitInOut::init(void) {
version=HMdcTrackDSet::getFitVersion();
if(version==0) {
Warning("init","track fit version 0 is not supported more, version 1 will used");
version=1;
}
if(!useRndbPar) fitPar = 0;
else fitPar = (HMdcTrackFitPar*)gHades->getRuntimeDb()->getContainer("MdcTrackFitPar");
fprint = HMdcTrackDSet::fPrint();
HMdcGetContainers* fGetCont = HMdcGetContainers::getObject();
fDriftTimePar = HMdcDriftTimePar::getObject();
HMdcWireData::setDriftTimePar(fDriftTimePar);
wireOffsetFlag=HMdcTrackDSet::getUseWireOffset();
if(wireOffsetFlag) {
fDigitPar=(HMdcDigitPar*)gHades->getRuntimeDb()->getContainer("MdcDigitPar");
if(!fDigitPar) {
Error("init:","Zero pointer for HMdcDigitPar recieved!");
return kFALSE;
}
} else fDigitPar=0;
fSizesCells=HMdcSizesCells::getObject();
if (!fSizesCells) {
Error("init","HMdcSizesCells is absent");
return kFALSE;
}
geantFlag=HMdcGetContainers::isGeant();
if(geantFlag && HMdcTrackDSet::fNTuple()) {
fGeantKineCat = fGetCont->getCatGeantKine();
fGeantMdcCat = fGetCont->getCatGeantMdc();
} else {
fGeantKineCat = 0;
fGeantMdcCat = 0;
}
fCalCat = fGetCont->getCatMdcCal1();
if (!fCalCat) return kFALSE;
if(HMdcTrackDSet::fNTuple()) {
fClusFitCat = fGetCont->getCatMdcClusFit(kTRUE);
fWireFitCat = fGetCont->getCatMdcWireFit(kTRUE);
if(!fClusFitCat || !fWireFitCat) return kFALSE;
} else {
fClusFitCat = 0;
fWireFitCat = 0;
}
locClusFit.set(1,0);
locWireFit.set(1,0);
loc.set(4,0,0,0,0);
return kTRUE;
}
Bool_t HMdcTrackFitInOut::reinit(void) {
if(!fSizesCells->initContainer()) return kFALSE;
if(!fDriftTimePar->initContainer()) return kFALSE;
for(Int_t sec=0; sec<6; sec++) {
HMdcSizesCellsSec& fSCSec=(*fSizesCells)[sec];
for(Int_t mod=0; mod<4; mod++)
fSCModAr[sec][mod]=(fSizesCells->modStatus(sec,mod)) ? &fSCSec[mod] : 0;
}
if( fitPar ) {
if( !HMdcGetContainers::isInited(fitPar) ) return kFALSE;
else {
if(useRtdbTofFlag) tofFlag = fitPar->getTofFlag();
cutWeight = fitPar->getCutWeight ();
if(useRtdbTScFlag) doTargScan = fitPar->getDoTargScan();
minTimeOffset = fitPar->getMinTimeOffset ();
maxTimeOffset = fitPar->getMaxTimeOffset ();
minCellsNum = fitPar->getMinCellsNum ();
chi2CutFlag = fitPar->getChi2CutFlag ();
totalChi2Cut = fitPar->getTotalChi2Cut ();
chi2PerNdfCut = fitPar->getChi2PerNdfCut ();
if(useRtdbTFlag) useTukeyFlag = fitPar->getUseTukeyFlag();
cnWs = fitPar->getCnWs ();
cn2s = fitPar->getCn2s ();
cn4s = fitPar->getCn4s ();
cn2s = cn4s*cn4s/cnWs;
minSig2 = fitPar->getMinSig2 ();
maxNFilterIter= fitPar->getMaxNFilterIter();
minWeight = fitPar->getMinWeight ();
maxChi2 = fitPar->getMaxChi2 ();
minTOffsetIter= fitPar->getMinTOffsetIter();
funCt1 = fitPar->getFunCt1 ();
stepD1 = fitPar->getStepD1 ();
funCt2 = fitPar->getFunCt2 ();
stepD2 = fitPar->getStepD2 ();
stepD3 = fitPar->getStepD3 ();
}
}
if(wireOffsetFlag) {
if( !HMdcGetContainers::isInited(fDigitPar) ) return kFALSE;
signalSpeed=fDigitPar->getSignalSpeed();
} else signalSpeed=0.;
return kTRUE;
}
void HMdcTrackFitInOut::setDefaultFitParam(void) {
useRndbPar = kFALSE;
fprint = kFALSE;
printStartEv = 0;
nEventPrint = 1000000000;
version = 1;
geantFlag = kFALSE;
wireOffsetFlag = 0;
signalSpeed = 0.005;
cutWeight = 0.01;
useRtdbTofFlag = kTRUE;
tofFlag = 3;
doTargScan = kTRUE;
calcInitValue = HMdcTrackDSet::getCalcInitValueFlag();
minTimeOffset = -30.;
maxTimeOffset = 60.;
minCellsNum = 5;
chi2CutFlag = kTRUE;
totalChi2Cut = 300.;
chi2PerNdfCut = 50.;
useTukeyFlag = kTRUE;
useRtdbTFlag = kTRUE;
useRtdbTScFlag = kTRUE;
cnWs = 6.4516;
cn4s = 10.6276;
cn2s = cn4s*cn4s/cnWs;
minSig2 = 2.5*2.5;
tukeyScale = 1.;
maxNFilterIter = 4;
minWeight = 0.5;
maxChi2 = 25.; ;
minTOffsetIter = -50.;
funCt1 = 500.;
stepD1 = 0.0001;
funCt2 = 10000.0;
stepD2 = 0.001;
stepD3 = 0.01;
for(Int_t s=0;s<6;s++) for(Int_t m=0;m<4;m++) for(Int_t l=0;l<6;l++) exclLay[s][m][l] = kFALSE;
}
void HMdcTrackFitInOut::setNEventsPrint(Int_t start,Int_t nev) {
printStartEv = start;
nEventPrint = nev;
}
void HMdcTrackFitInOut::setTukeyConstants(Double_t cw,Double_t c2,Double_t c4) {
cnWs = cw*cw;
cn4s = c4*c4;
cn2s = cn4s*cn4s/cnWs;
}
void HMdcTrackFitInOut::getTukeyConstants(Double_t& cw,Double_t& c2,
Double_t& c4) const {
cw = sqrt(cnWs);
c2 = sqrt(cn2s);
c4 = sqrt(cn4s);
}
HMdcWireFit* HMdcTrackFitInOut::getNewWireFitSlot(Int_t* indWireFit) {
TObject* fWireFit=fWireFitCat->getNewSlot(locWireFit,indWireFit);
if(!fWireFit) {
Warning("getNewWireFitSlot","No slot HMdcWireFit available");
return 0;
}
if(geantFlag) return new(fWireFit) HMdcWireFitSim;
else return new(fWireFit) HMdcWireFit;
}
HMdcClusFit* HMdcTrackFitInOut::getNewClusFitSlot(Int_t* indClusFit) {
TObject* fClusFit=fClusFitCat->getNewSlot(locClusFit,indClusFit);
if(!fClusFit) {
Warning("getNewClusFitSlot","No slot HMdcClusFit available");
return 0;
}
if(geantFlag) return new(fClusFit) HMdcClusFitSim;
else return new(fClusFit) HMdcClusFit;
}
Double_t HMdcTrackFitInOut::getStepDer(Double_t funct) const {
return (funct<funCt1) ? stepD1:((funct<funCt2) ? stepD2:stepD3);
}
void HMdcTrackFitInOut::setPrintFlag(Int_t currentEvent) {
fprint = currentEvent>=printStartEv && currentEvent<printStartEv+nEventPrint;
}
void HMdcTrackFitInOut::setChi2PerNdfCut(Double_t cut) {
chi2CutFlag = kFALSE;
chi2PerNdfCut = cut;
}
void HMdcTrackFitInOut::setTotalChi2Cut(Double_t cut) {
chi2CutFlag = kTRUE;
totalChi2Cut = cut;
}
void HMdcTrackFitInOut::excludeLayer(UInt_t s,UInt_t m,UInt_t l) {
if(s<6 && m<4 && l<6) exclLay[s][m][l] = kTRUE;
}
void HMdcTrackFitInOut::excludeModule(UInt_t s,UInt_t m) {
if(s<6 && m<4) for(UInt_t l=0;l<6;l++) exclLay[s][m][l] = kTRUE;
}
Bool_t HMdcTrackFitInOut::isLayerExcluded(Int_t s,Int_t m,Int_t l) const {
return s<6 && m<4 && l<6 ? exclLay[s][m][l] : kFALSE;
}
HMdcTrackFitter::HMdcTrackFitter(HMdcTrackFitInOut* fIO) {
fitInOut=fIO;
init();
}
void HMdcTrackFitter::init(void) {
fprint = fitInOut->getPrintFlag();
tofFlag = fitInOut->getTofFlag();
wires.setPrintFlag(fprint);
wires.setTrackFitInOut(fitInOut);
}
void HMdcTrackFitter::setPrintFlag(Bool_t prnt) {
fprint=prnt;
wires.setPrintFlag(fprint);
}
Bool_t HMdcTrackFitter::fillListHits(HMdcClus* cl1,HMdcClus* cl2) {
segIndex=-1;
indClusFit=-1;
if(!wires.fillListHits(cl1,cl2)) return kFALSE;
setPlanes();
return kTRUE;
}
Bool_t HMdcTrackFitter::fillListHits(HMdcClus* cl1,HMdcClus* cl2,HMdcClus* cl3,HMdcClus* cl4,
HMdcClus* cl5,HMdcClus* cl6,HMdcClus* cl7,HMdcClus* cl8) {
segIndex=-1;
indClusFit=-1;
if( !wires.fillListHits(cl1,cl2,cl3,cl4,cl5,cl6,cl7,cl8) ) return kFALSE;
setPlanes();
return kTRUE;
}
Bool_t HMdcTrackFitter::fillListHits(HMdcEvntListCells* store,HMdcClus* clus1,HMdcClus* clus2) {
segIndex=-1;
indClusFit=-1;
if(!wires.fillListHits(store,clus1,clus2)) return kFALSE;
setPlanes();
return kTRUE;
}
Bool_t HMdcTrackFitter::fillListHits(HMdcEvntListCells* store) {
segIndex=-1;
indClusFit=-1;
if(!wires.fillListHits(store)) return kFALSE;
setPlanes();
return kTRUE;
}
void HMdcTrackFitter::setPlanes(void) {
Int_t sec = wires.getSector();
HMdcSizesCellsMod** fSCModAr = fitInOut->getSCellsModArr(sec);
Int_t nClTimes = wires.getNDrTmFirstSec();
initParam.setFirstPlane(&((*(fSCModAr[wires[0].getModule()]))[0]));
initParam.setSecondPlane(&((*(fSCModAr[wires[nClTimes-1].getModule()]))[5]));
initParam.setCoorSys(sec);
}
Bool_t HMdcTrackFitter::fillClusFitCont(void) {
HMdcClusFit* fClusFit=fitInOut->getNewClusFitSlot(&indClusFit);
if(!fClusFit) return kFALSE;
fClusFit->setFitStatus(fitStatus);
finalParam.fillClusFit(fClusFit);
fClusFit->setExitFlag(exitFlag);
wires.calcDistanceSign(finalParam);
wires.fillClusFitSim(fClusFit,finalParam);
wires.fillClusFitAndWireFit(fClusFit);
return kTRUE;
}
Bool_t HMdcTrackFitter::fitCluster(Int_t fittingMod) {
if(wires.getNCellsInInput(fittingMod) < 5) return kFALSE;
wires.setHitStatM1toP1();
wires.subtractWireOffset(initParam);
if(HMdcTrackDSet::getCalcInitValueFlag()==1) wires.calcInitialValue(initParam);
wires.fillLookupTableForDer(initParam);
fitStatus = fit(fittingMod);
if(HMdcTrackDSet::getCalcInitValueFlag()==2 && wires.getNWiresInFit()<30) {
if(!fitStatus) refitCluster(fittingMod);
else {
Int_t nInputLayers = wires.getInputListCells().getNLayers();
Int_t nFittedLayers = wires.getOutputListCells().getNLayers();
if(nInputLayers-nFittedLayers >= 3) {
refitCluster(fittingMod);
nFittedLayers = wires.getOutputListCells().getNLayers();
}
}
}
if(fitInOut->getClusFitCat() && HMdcTrackDSet::fNTuple()) fillClusFitCont();
return fitStatus;
}
void HMdcTrackFitter::refitCluster(Int_t fittingMod) {
if(wires.calcInitialValue(initParam)) {
wires.reinitWtSt();
fitStatus = fit(fittingMod);
}
}
Bool_t HMdcTrackFitter::fit(Int_t fittingMod) {
Int_t iter = 0;
while(kTRUE) {
Int_t exit = minimize(iter++);
Double_t delW = wires.testFitResult();
Int_t nCells = wires.getNCellsInOutput(fittingMod);
if( delW<0.5 || nCells<6 ) {
if(exit == 0) return kFALSE;
if(delW > 0.)
if(finalParam.calcChi2PerDF(wires.getNumOfGoodWires())<0. || !testChi2Cut()) return kFALSE;
if(finalParam.testParameters(fitInOut->getMinTimeOffset(),fitInOut->getMaxTimeOffset()) &&
nCells >= fitInOut->getMinCellsNum()) return kTRUE;
break;
}
if(fprint) printf("TestFit: num.of deleted cells=%.1f, refit this!\n",delW);
}
return kFALSE;
}
Bool_t HMdcTrackFitter::testChi2Cut(void) {
if(fitInOut->getChi2CutFlag()) {
if(finalParam.functional()<fitInOut->getTotalChi2Cut() ) return kTRUE;
} else if(finalParam.getChi2() <fitInOut->getChi2PerNdfCut()) return kTRUE;
return kFALSE;
}
void HMdcTrackFitter::setRegionOfWires(Int_t mod) {
wires.setRegionOfWires(mod);
}
Bool_t HMdcTrackFitter::setClustAndFill(HMdcClus* cl1, HMdcClus* cl2) {
if(fprint != fitInOut->getPrintFlag()) setPrintFlag(fitInOut->getPrintFlag());
if(fprint) {
cl1->print();
if(cl2) cl2->print();
}
if(!fillListHits(cl1,cl2)) return kFALSE;
initParam.setParam(cl1->getXTarg(),cl1->getYTarg(),cl1->getZTarg(),
cl1->getX(), cl1->getY(), cl1->getZ());
return kTRUE;
}
Bool_t HMdcTrackFitter::setClustAndFill(HMdcClus* cl1,HMdcClus* cl2,HMdcClus* cl3,HMdcClus* cl4,
HMdcClus* cl5,HMdcClus* cl6,HMdcClus* cl7,HMdcClus* cl8) {
if(fprint != fitInOut->getPrintFlag()) setPrintFlag(fitInOut->getPrintFlag());
if(fprint) {
cl1->print();
if(cl2 != NULL) cl2->print();
if(cl3 != NULL) cl3->print();
if(cl4 != NULL) cl4->print();
if(cl5 != NULL) cl5->print();
if(cl6 != NULL) cl6->print();
if(cl7 != NULL) cl7->print();
if(cl8 != NULL) cl8->print();
}
Int_t nSecs = 2;
if(cl5 != NULL) nSecs = 3;
if(cl7 != NULL) nSecs = 4;
initParam.setNMods(nSecs);
finalParam.setNMods(nSecs);
if(!fillListHits(cl1,cl2,cl3,cl4,cl5,cl6,cl7,cl8)) return kFALSE;
initParam.setParam(cl1->getXTarg(),cl1->getYTarg(),cl1->getZTarg(),
cl1->getX(), cl1->getY(), cl1->getZ());
return kTRUE;
}
Bool_t HMdcTrackFitter::setClustAndFill(HMdcEvntListCells* store,
Double_t x1, Double_t y1, Double_t z1,
Double_t x2, Double_t y2, Double_t z2) {
if(fprint != fitInOut->getPrintFlag()) setPrintFlag(fitInOut->getPrintFlag());
if(fprint) printf("x1=%f y1=%f z1=%f x2=%f y2=%f z2=%f \n",x1,y1,z1,x2,y2,z2);
if(!fillListHits(store)) return kFALSE;
initParam.setParam(x1, y1, z1, x2, y2, z2);
wires.setXYZClust(x2, y2, z2);
return kTRUE;
}
void HMdcTrackFitter::fillOutput() {
finalParam.copyLine(initParam);
wires.calcNGoodWiresAndChi2(finalParam);
wires.valueOfFunctional(finalParam,2);
wires.calculateErrors(finalParam);
}