#include "halignmentmeta.h"
#include "TMinuit.h"
#include "TFile.h"
#include "TString.h"
#include "TNtuple.h"
#include "hmdcsizescells.h"
#include "hades.h"
#include "hruntimedb.h"
#include "hspecgeompar.h"
#include "hgeomvolume.h"
#include <iostream>
using namespace std;
ClassImp(HAlignmentMeta)
HAlignmentMeta::HAlignmentMeta() {
trackSelecCutX = 1.5;
trackSelecCutY = 1.0;
filterFlag = kTRUE;
nMetaModules = 0;
metaDetector = 3;
fitTofModYPos = kFALSE;
xShitfRpc = 0.;
calcCellXOffset = kFALSE;
nCells = 0;
nCellsTot = 0;
}
void HAlignmentMeta::setTofDetector(Double_t cutX,Double_t cutY) {
metaDetector = 0;
nMetaModules = 8;
setCuts(cutX,cutY);
nCells = 8;
nCellsTot = nMetaModules*nCells;
}
void HAlignmentMeta::setShowerDetector(Double_t cutX,Double_t cutY) {
metaDetector = 1;
nMetaModules = 1;
setCuts(cutX,cutY);
nCells = 32;
nCellsTot = nCells;
}
void HAlignmentMeta::setRpcDetector(Double_t cutX,Double_t cutY) {
metaDetector = 2;
nMetaModules = 1;
setCuts(cutX,cutY);
nCells = 32;
nCellsTot = 6*nCells;
}
void HAlignmentMeta::setCuts(Double_t cutX,Double_t cutY) {
trackSelecCutX = cutX;
trackSelecCutY = cutY;
if(trackSelecCutX > 0. && trackSelecCutY > 0.) filterFlag = kTRUE;
else filterFlag = kFALSE;
}
void HAlignmentMeta::setNtuple(TNtuple *ntp) {
nt = ntp;
nt -> SetBranchAddress("s",&sec);
nt -> SetBranchAddress("x1s",&x1);
nt -> SetBranchAddress("y1s",&y1);
nt -> SetBranchAddress("z1s",&z1);
nt -> SetBranchAddress("x2s",&x2);
nt -> SetBranchAddress("y2s",&y2);
nt -> SetBranchAddress("z2s",&z2);
nt -> SetBranchAddress("metaMod",&metaModule);
nt -> SetBranchAddress("col",&metaColumn);
nt -> SetBranchAddress("cell",&metaCell);
nt -> SetBranchAddress("xMeta",&xMetaLocal);
nt -> SetBranchAddress("yMeta",&yMetaLocal);
nt -> SetBranchAddress("zMeta",&zMetaLocal);
nt -> SetBranchAddress("sigX",&xRMS);
nt -> SetBranchAddress("sigY",&yRMS);
nt -> SetBranchAddress("sigZ",&zRMS);
}
HAlignmentMeta::~HAlignmentMeta() {
if(gMinuit != NULL) delete gMinuit;
}
void HAlignmentMeta::fcnMeta(Int_t &npar, Double_t *gin, Double_t &fn, Double_t *par, Int_t iflag) {
static int count = 1;
count++;
HAlignmentMeta *myObject = (HAlignmentMeta*)(gMinuit -> GetObjectFit());
fn = myObject->getMinFunction(par);
if(count%100 == 0) cout<<"iter: "<<count<<" function "<<fn<<" "<<par[0]<<" "<<par[1]<<" "
<<par[2]<<" "<<par[3]<<" "<<par[4]<<" "<<par[5]<< endl;
}
void HAlignmentMeta::alignMeta(Int_t sec,TNtuple *ntp) {
if(metaDetector<0 || metaDetector>2) {
Error("alignMeta","Meta detector type is not setted! Stop!");
exit(1);
}
alignSec = sec;
setNtuple(ntp);
HSpecGeomPar* fSpecGeomPar = (HSpecGeomPar*)gHades->getRuntimeDb()->getContainer("SpecGeomPar");
const HGeomTransform& secLabTrans = fSpecGeomPar->getSector(alignSec)->getTransform();
xShitfRpc = 0.;
for(Int_t im=0;im<nMetaModules;im++) {
transMetaModSecOld[im] = transMetaModLabOld[im];
transMetaModSecOld[im].transTo(secLabTrans);
}
Int_t ierflg = 0;
if(gMinuit == NULL) {
if(metaDetector == 0) new TMinuit(6+8);
else new TMinuit(6);
}
gMinuit->mncler();
gMinuit->SetFCN(fcnMeta);
gMinuit->SetObjectFit(this);
Double_t arglist[6];
fillArray();
gMinuit->mnparm(0,"x", 0., 1.,0,0,ierflg);
gMinuit->mnparm(1,"y", 0., 1.,0,0,ierflg);
gMinuit->mnparm(2,"z", 0., 1.,0,0,ierflg);
gMinuit->mnparm(3,"alpha", 0., 0.1,0,0,ierflg);
gMinuit->mnparm(4,"beta", 0., 0.1,0,0,ierflg);
gMinuit->mnparm(5,"gamma", 0., 0.1,0,0,ierflg);
if(metaDetector != 1) {
gMinuit->FixParameter(0);
gMinuit->FixParameter(4);
if(metaDetector == 0) {
if(fitTofModYPos) {
for(Int_t ip=0;ip<8;ip++) {
TString parName("Ymod");
parName += ip;
gMinuit->mnparm(6+ip,parName.Data(),0.,1.,0,0,ierflg);
}
for(Int_t ip=0;ip<8;ip++) gMinuit->FixParameter(6+ip);
for(Int_t ip=0;ip<8;ip++) tofModYSh[ip] = 0.;
}
}
}
gMinuit->SetPrintLevel(1);
arglist[0] = 2;
gMinuit->mnexcm("SET STR",arglist,1,ierflg);
arglist[1] = 0.001;
arglist[0] = 20000;
for(Int_t i=0;i<5;i++) {
calcMinDist();
if(filterFlag) {
if(metaDetector == 1) selectTracks(i==0 ? 2. : 1.);
else calcXOffset(i==0 ? 2. : 1.);
}
setWeights();
if(i==0) for(Int_t tr = 0; tr < nTracks; tr++) {
tracks[tr].xMinDistInit = tracks[tr].xMinDist;
tracks[tr].yMinDistInit = tracks[tr].yMinDist;
}
gMinuit->mnexcm("MINI",arglist,2,ierflg);
if(i<2) continue;
if(ierflg == 0) break;
}
calcMinDist();
if(metaDetector == 0 && fitTofModYPos) {
printf(" ------ Fit Yshifts --------------\n");
gMinuit->mnfree(0);
for(Int_t ip=0;ip<6;ip++) gMinuit->FixParameter(ip);
gMinuit->mnexcm("MINI",arglist,2,ierflg);
selectTracks(1.);
setWeights();
gMinuit->mnfree(0);
for(Int_t ip=0;ip<8;ip++) gMinuit->FixParameter(6+ip);
gMinuit->FixParameter(0);
gMinuit->FixParameter(4);
gMinuit->mnexcm("MINI",arglist,2,ierflg);
calcMinDist();
}
if(metaDetector != 1) {
calcXOffset(1.0);
calcMinDist();
}
for(Int_t im=0;im<nMetaModules;im++) {
transMetaModLabNew[im] = transMetaModSecNew[im];
transMetaModLabNew[im].transFrom(secLabTrans);
}
}
void HAlignmentMeta::calcMinDist(void) {
Double_t out[6];
Double_t err;
for(Int_t ip=0;ip<6;ip++) gMinuit->GetParameter(ip,out[ip],err);
if(metaDetector==0 && fitTofModYPos)
for(Int_t ip=0;ip<8;ip++) gMinuit->GetParameter(6+ip,tofModYSh[ip],err);
HGeomTransform trans;
HMdcSizesCells::setTransform(out,trans);
calcMinDist(trans);
}
void HAlignmentMeta::calcMinDist(Double_t *par) {
HGeomTransform trans;
HMdcSizesCells::setTransform(par,trans);
if(metaDetector == 0 && fitTofModYPos) for(Int_t ip=0;ip<8;ip++) tofModYSh[ip] = par[6+ip];
calcMinDist(trans);
}
void HAlignmentMeta::fillArray(void) {
nTracks = 0;
for(Int_t c=0;c<nCellsTot;c++) {
cellStat[c] = 0;
cellXCorr[c] = 0.;
}
Int_t nentries = nt->GetEntries();
HGeomVector trPoint1pr,trPoint2pr;
Int_t strInd = 0;
for(Int_t i = 0; i < nentries; i++) {
nt->GetEntry(i);
if(sec != alignSec) continue;
if(TMath::IsNaN(xMetaLocal+yMetaLocal+zMetaLocal)) {
printf("NaN!!! Entry %i: xMetaLocal=%f yMetaLocal=%f zMetaLocal=%f\n",
i,xMetaLocal,yMetaLocal,zMetaLocal);
continue;
}
TrackMdcMeta &t = tracks[nTracks];
t.mdcPnt1Sec.setXYZ(x1, y1, z1);
t.mdcPnt2Sec.setXYZ(x2, y2, z2);
t.metaMod = Short_t(metaModule + 0.1);
t.column = Short_t(metaColumn+0.1);
t.cell = Short_t(metaCell+0.1);
if (metaDetector == 0) t.cellInd = t.metaMod*nCells + t.cell;
else if(metaDetector == 1) t.cellInd = t.cell;
else if(metaDetector == 2) t.cellInd = t.column*nCells + t.cell;
t.xMeta = xMetaLocal;
t.yMeta = yMetaLocal;
t.zMeta = zMetaLocal;
t.sigmaX = xRMS;
t.sigmaY = yRMS;
t.sigmaZ = zRMS;
t.useIt = kTRUE;
t.wt = 1.;
HGeomVector dv1(trPoint1pr-t.mdcPnt1Sec);
HGeomVector dv2(trPoint2pr-t.mdcPnt2Sec);
if(dv1.length()+dv2.length()<0.01) t.startTrInd = strInd;
else {
t.startTrInd = nTracks;
strInd = nTracks;
trPoint1pr = t.mdcPnt1Sec;
trPoint2pr = t.mdcPnt2Sec;
}
if(nTracks==0 || yMetaLocal<yMinMetaLocal) yMinMetaLocal = yMetaLocal;
if(nTracks==0 || yMetaLocal>yMaxMetaLocal) yMaxMetaLocal = yMetaLocal;
cellStat[t.cellInd]++;
nTracks++;
if(nTracks==1000000) return;
}
if(metaDetector > 0) {
Double_t binSize = (yMaxMetaLocal-yMinMetaLocal)/8.;
for(Int_t tr = 0; tr < nTracks; tr++) {
TrackMdcMeta &t = tracks[tr];
t.binScWt = Short_t((t.yMeta-yMinMetaLocal)/binSize);
if(t.binScWt<0) t.binScWt = 0;
else if(t.binScWt>7) t.binScWt = 7;
}
} else {
for(Int_t tr = 0; tr < nTracks; tr++) {
TrackMdcMeta &t = tracks[tr];
t.binScWt = t.metaMod;
}
}
}
Double_t HAlignmentMeta::getMinFunction(Double_t *par) {
calcMinDist(par);
Double_t dist = 0;
for(Int_t tr = 0; tr < nTracks; tr++) if(tracks[tr].useIt) {
dist += tracks[tr].minDist2*tracks[tr].wt;
}
return dist;
}
void HAlignmentMeta::calcMinDist(HGeomTransform& trans) {
for(Int_t im=0;im<nMetaModules;im++) {
transMetaModSecNew[im] = transMetaModSecOld[im];
transMetaModSecNew[im].transTo(trans);
if(metaDetector == 2 && xShitfRpc != 0.) {
HGeomTransform shifterTrans;
HGeomVector tr(xShitfRpc,0.,0.);
shifterTrans.setTransVector(tr);
shifterTrans.transFrom(transMetaModSecNew[im]);
transMetaModSecNew[im] = shifterTrans;
}
if(metaDetector == 0 && fitTofModYPos && tofModYSh[im] != 0.) {
HGeomTransform trans;
HGeomVector tr(0.,tofModYSh[im],0.);
trans.setTransVector(tr) ;
trans.transFrom(transMetaModSecNew[im]);
transMetaModSecNew[im] = trans;
}
}
for(Int_t tr = 0; tr < nTracks; tr++) {
TrackMdcMeta &t = tracks[tr];
HGeomVector point1(transMetaModSecNew[t.metaMod].transTo(t.mdcPnt1Sec));
HGeomVector point2(transMetaModSecNew[t.metaMod].transTo(t.mdcPnt2Sec));
point1.setZ(point1.getZ() - t.zMeta);
point2.setZ(point2.getZ() - t.zMeta);
Double_t diffZ = point2.getZ() - point1.getZ();
Double_t diffX = point2.getX() - point1.getX();
Double_t diffY = point2.getY() - point1.getY();
t.xMinDist = (point1.getX()*diffZ - point1.getZ()*diffX)/diffZ - t.xMeta;
t.yMinDist = (point1.getY()*diffZ - point1.getZ()*diffY)/diffZ - t.yMeta;
t.xMinDist -= cellXCorr[t.cellInd];
Double_t dYNr = t.dYNorm();
t.minDist2 = dYNr*dYNr;
if(metaDetector == 1) {
Double_t dXNr = t.dXNorm();
t.minDist2 += dXNr*dXNr;
}
}
}
void HAlignmentMeta::setWeights(void) {
Double_t stat[8];
for(Int_t i=0;i<8;i++) stat[i] = 0;
for(Int_t tr = 0; tr < nTracks; tr++) {
TrackMdcMeta &t = tracks[tr];
t.wt = t.useIt ? 1.:0.;
if(t.useIt && t.startTrInd != tr) {
Int_t nTr = 1;
Int_t str = t.startTrInd;
for(Int_t p=str;p<tr;p++) if(tracks[p].useIt && t.oneLay(tracks[p])) nTr++;
if(nTr == 1) continue;
Double_t w = 1./Double_t(nTr);
t.wt = w;
for(Int_t p=str;p<tr;p++) if(tracks[p].useIt && t.oneLay(tracks[p])) tracks[p].wt = w;
}
stat[t.binScWt] += t.wt;
}
Double_t statMax = 0.;
for(Int_t i=0;i<8;i++) if(statMax < stat[i]) statMax = stat[i];
for(Int_t tr = 0; tr < nTracks; tr++) {
TrackMdcMeta &t = tracks[tr];
if(t.wt>0. && stat[t.binScWt]>100.) t.wt *= statMax/stat[t.binScWt];
}
}
void HAlignmentMeta::selectTracks(Double_t nSigmasCut) {
meanX = 0.;
meanY = 0.;
meanZ = 0.;
sigmX = 0.;
sigmY = 0.;
sigmZ = 0.;
isFirstSIter = kTRUE;
while(selectTracksIter(nSigmasCut));
}
Bool_t HAlignmentMeta::selectTracksIter(Double_t nSigmasCut) {
Double_t dXM = 0.;
Double_t dX2M = 0.;
Double_t dYM = 0.;
Double_t dY2M = 0.;
Int_t nTr = 0;
Bool_t doNextIter = kFALSE;
for(Int_t tr = 0; tr < nTracks; tr++) {
TrackMdcMeta &t = tracks[tr];
Double_t dXm = t.dXNorm();
Double_t dYm = t.dYNorm();
if(isFirstSIter) {
t.useIt = kTRUE;
doNextIter = kTRUE;
} else if(TMath::Abs(dXm-meanX) <= sigmX*nSigmasCut &&
TMath::Abs(dYm-meanY) <= sigmY*nSigmasCut) {
if(!t.useIt) doNextIter = kTRUE;
t.useIt = kTRUE;
} else {
if(t.useIt) doNextIter = kTRUE;
t.useIt = kFALSE;
continue;
}
dXM += dXm;
dX2M += dXm*dXm;
dYM += dYm;
dY2M += dYm*dYm;
nTr++;
}
meanX = dXM/nTr;
meanY = dYM/nTr;
sigmX = TMath::Max(TMath::Sqrt(dX2M/nTr - meanX*meanX),trackSelecCutX);
sigmY = TMath::Max(TMath::Sqrt(dY2M/nTr - meanY*meanY),trackSelecCutY);
if(doNextIter) printf(" * From %i tracks %i selected. <X>=%f sX=%f <Y>=%f sY=%f\n",
nTracks,nTr,meanX,sigmX,meanY,sigmY);
isFirstSIter = kFALSE;
return doNextIter;
}
void HAlignmentMeta::checkAlignment(void) {
TString fileName("align");
if(metaDetector == 0) fileName += "Tof_S";
else if(metaDetector == 1) fileName += "Shower_S";
else if(metaDetector == 2) fileName += "Rpc_S";
fileName += alignSec;
fileName += ".root";
TFile *f = new TFile(fileName.Data(),"recreate");
TNtuple *chnt = new TNtuple("chnt","chnt",
"alignSec:wt:xMeta:yMeta:zMeta:dXold:dYold:dX:dY:sigX:sigY:chi2:mod:col:cell");
Float_t data[15];
for(Int_t tr = 0; tr < nTracks; tr++) {
TrackMdcMeta &t = tracks[tr];
data[ 0] = sec;
data[ 1] = t.wt;
data[ 2] = t.xMeta + cellXCorr[t.cellInd];
data[ 3] = t.yMeta;
data[ 4] = t.zMeta;
data[ 5] = t.xMinDistInit;
data[ 6] = t.yMinDistInit;
data[ 7] = t.xMinDist;
data[ 8] = t.yMinDist;
data[ 9] = t.sigmaX;
data[10] = t.sigmaY;
data[11] = t.minDist2;
data[12] = t.metaMod;
data[13] = t.column;
data[14] = t.cell;
chnt->Fill(data);
}
f->cd();
chnt->Write();
f->Close();
delete f;
}
void HAlignmentMeta::calcXOffset(Double_t nSigmasCut) {
Int_t nCellsC = 0;
xShitfRpc = 0.;
for(Short_t c=0;c<nCellsTot;c++) {
meanX = 0.;
meanY = 0.;
meanZ = 0.;
sigmX = 0.;
sigmY = 0.;
sigmZ = 0.;
isFirstSIter = kTRUE;
if(cellStat[c] > 100) {
while(calcXOffset(nSigmasCut,c));
nCellsC++;
if(metaDetector==2 && calcCellXOffset) {
xShitfRpc += cellXCorr[c];
nCellsC++;
}
}
}
if(metaDetector==2 && calcCellXOffset ) {
if(nCellsC>0) xShitfRpc /= nCellsC;
printf("xShitf =%f nCells=%i !!!!!!!!!!!\n",xShitfRpc,nCellsC);
for(Short_t c=0;c<nCellsTot;c++) if(cellStat[c] > 100) cellXCorr[c] -= xShitfRpc;
}
}
Bool_t HAlignmentMeta::calcXOffset(Double_t nSigmasCut,Short_t cellInd) {
Double_t dXM = 0.;
Double_t dX2M = 0.;
Double_t dYM = 0.;
Double_t dY2M = 0.;
Double_t xShift = 0.;
Int_t nTr = 0;
Bool_t doNextIter = kFALSE;
Int_t nTot = 0;
for(Int_t tr = 0; tr < nTracks; tr++) {
TrackMdcMeta &t = tracks[tr];
if(cellInd != t.cellInd) continue;
nTot++;
Double_t dXm = t.dXNorm();
Double_t dYm = t.dYNorm();
if(isFirstSIter) {
t.useIt = kTRUE;
doNextIter = kTRUE;
} else if(TMath::Abs(dXm-meanX) <= sigmX*nSigmasCut &&
TMath::Abs(dYm-meanY) <= sigmY*nSigmasCut) {
if(!t.useIt) doNextIter = kTRUE;
t.useIt = kTRUE;
} else {
if(t.useIt) doNextIter = kTRUE;
t.useIt = kFALSE;
continue;
}
dXM += dXm;
dX2M += dXm*dXm;
dYM += dYm;
dY2M += dYm*dYm;
xShift += t.xMinDist;
nTr++;
}
meanX = dXM/nTr;
meanY = dYM/nTr;
sigmX = TMath::Max(TMath::Sqrt(dX2M/nTr - meanX*meanX),trackSelecCutX);
sigmY = TMath::Max(TMath::Sqrt(dY2M/nTr - meanY*meanY),trackSelecCutY);
isFirstSIter = kFALSE;
if(!doNextIter) {
cellXCorr[cellInd] += xShift/nTr;
printf(" * %icolumn %2icell: From %i tracks %i selected. xOffset=%f\n",
cellInd/nCells,cellInd%nCells,nTot,nTr,cellXCorr[cellInd]);
}
return doNextIter;
}