Docker-in-Docker (DinD) capabilities of public runners deactivated. More info

Commit 61bd0452 authored by dino's avatar dino
Browse files

Some more work on inverting xtalk in PreprocessingFilterPSA

git-svn-id: svn://gal-serv.lnl.infn.it/agata/trunk/narval_emulator@1127 170316e4-aea8-4b27-aad4-0380ec0519c9
parent 58a703ca
......@@ -179,6 +179,11 @@ int main(int argc, char **argv)
cout << "# Using Narval actors without Narval #" << endl;
cout << "######################################" << endl;
getparams(argc, argv);
if(listKeys)
listKeysAndExit();
#ifdef WCT_THREADED
cout << "\nBoost C++ Library Version " << BOOST_VERSION/100000 << "_" << (BOOST_VERSION/100)%1000 << "_" << BOOST_VERSION%100 << endl;
#endif
......@@ -189,10 +194,6 @@ int main(int argc, char **argv)
if(!ok)
exit(EXIT_FAILURE);
getparams(argc, argv);
if(listKeys)
listKeysAndExit();
ok = topologyRead(topologyFile);
if(!ok)
......@@ -318,6 +319,10 @@ void getparams(int argc, char *argv[])
else if(!strcmp(argv[1], "-h")|| !strcmp(argv[1], "--help")) {
ok = -1;
}
else if(!strcmp(argv[1], "-k") || !strcmp(argv[1], "-keys")) {
listKeys = true;
return;
}
else {
topologyFile = argv[1];
ok = (argc > 1) ? 1 : 0;
......@@ -348,6 +353,7 @@ void getparams(int argc, char *argv[])
}
else if(!strcmp(cmd, "-k") || !strcmp(cmd, "-keys")) {
listKeys = true;
return;
}
else {
ok = 0;
......
......@@ -1873,11 +1873,10 @@ void PreprocessingFilterPSA::CoreEnergyFromXtalk(int segMult)
cxtalk += tt*tt;
}
}
extalk /= (float)cxtalk;
extalk /= cxtalk;
}
#if 1
// not really worthwhile
#if 0
else if(segMult == 2) {
double d00 = 0, d01 = 0;
double d10 = 0, d11 = 0;
......@@ -1906,7 +1905,37 @@ void PreprocessingFilterPSA::CoreEnergyFromXtalk(int segMult)
#endif
#if 1
//not really worthwhile
else if(segMult == 2) {
double dd[4]; memset(dd, 0, sizeof(dd)); // row-wise
double xd[2]; memset(xd, 0, sizeof(xd));
// 0 1
// 2 3
for(int ns = 0; ns < CC.nSG; ns++) {
float vv = SegE[ns];
if(vv < CC.pSG[ns].emink) {
float *px = xTalkDir[ns];
double tt0 = px[hitSeg[0]];
double tt1 = px[hitSeg[1]];
xd[0] += vv*tt0;
xd[1] += vv*tt1;
dd[0] += tt0*tt0;
dd[1] += tt0*tt1;
dd[3] += tt1*tt1;
}
}
dd[ 2] = dd[ 1];
double DD[4];
if(!InvertMatrix2(dd, DD))
return; // ??
double eSeg0 = DD[0]*xd[0] + DD[1]*xd[1];
double eSeg1 = DD[2]*xd[0] + DD[3]*xd[1];
extalk = float(eSeg0 + eSeg1);
}
#endif
#if 0
else if(segMult == 3) {
double d00 = 0, d01 = 0, d02 = 0;
double d10 = 0, d11 = 0, d12 = 0;
......@@ -1944,6 +1973,93 @@ void PreprocessingFilterPSA::CoreEnergyFromXtalk(int segMult)
}
#endif
#if 1
else if(segMult == 3) {
double dd[9]; memset(dd, 0, sizeof(dd)); // row-wise
double xd[3]; memset(xd, 0, sizeof(xd));
// 0 1 2
// 3 4 5
// 6 7 8
for(int ns = 0; ns < CC.nSG; ns++) {
float vv = SegE[ns];
if(vv < CC.pSG[ns].emink) {
float *px = xTalkDir[ns];
double tt0 = px[hitSeg[0]];
double tt1 = px[hitSeg[1]];
double tt2 = px[hitSeg[2]];
xd[0] += vv*tt0;
xd[1] += vv*tt1;
xd[2] += vv*tt2;
dd[0] += tt0*tt0;
dd[1] += tt0*tt1;
dd[2] += tt0*tt2;
dd[4] += tt1*tt1;
dd[5] += tt1*tt2;
dd[8] += tt2*tt2;
}
}
dd[ 3] = dd[ 1];
dd[ 6] = dd[ 2]; dd[ 7] = dd[ 5];
double DD[9];
if(!InvertMatrix3(dd, DD))
return; // ??
double eSeg0 = DD[0]*xd[0] + DD[1]*xd[1] + DD[2]*xd[2];
double eSeg1 = DD[3]*xd[0] + DD[4]*xd[1] + DD[5]*xd[2];
double eSeg2 = DD[6]*xd[0] + DD[7]*xd[1] + DD[8]*xd[2];
extalk = float(eSeg0 + eSeg1 + eSeg2);
}
#endif
#if 1
else if(segMult == 4) {
double dd[16]; memset(dd, 0, sizeof(dd)); // row-wise
double xd[ 4]; memset(xd, 0, sizeof(xd));
// 0 1 2 3
// 4 5 6 7
// 8 9 10 11
//12 13 14 15
for(int ns = 0; ns < CC.nSG; ns++) {
float vv = SegE[ns];
if(vv < CC.pSG[ns].emink) {
float *px = xTalkDir[ns];
double tt0 = px[hitSeg[0]];
double tt1 = px[hitSeg[1]];
double tt2 = px[hitSeg[2]];
double tt3 = px[hitSeg[3]];
xd[ 0] += vv*tt0;
xd[ 1] += vv*tt1;
xd[ 2] += vv*tt2;
xd[ 3] += vv*tt3;
dd[ 0] += tt0*tt0;
dd[ 1] += tt0*tt1;
dd[ 2] += tt0*tt2;
dd[ 3] += tt0*tt3;
dd[ 5] += tt1*tt1;
dd[ 6] += tt1*tt2;
dd[ 7] += tt1*tt3;
dd[10] += tt2*tt2;
dd[11] += tt2*tt3;
dd[15] += tt3*tt3;
}
}
dd[ 4] = dd[ 1];
dd[ 8] = dd[ 2]; dd[ 9] = dd[ 6];
dd[12] = dd[ 3]; dd[13] = dd[ 7]; dd[14] = dd[11];
double DD[16];
if(!InvertMatrix4(dd, DD))
return; // ??
double eSeg0 = DD[ 0]*xd[0] + DD[ 1]*xd[1] + DD[ 2]*xd[2] + DD[ 3]*xd[3];
double eSeg1 = DD[ 4]*xd[0] + DD[ 5]*xd[1] + DD[ 6]*xd[2] + DD[ 7]*xd[3];
double eSeg2 = DD[ 8]*xd[0] + DD[ 9]*xd[1] + DD[10]*xd[2] + DD[11]*xd[3];
double eSeg3 = DD[12]*xd[0] + DD[13]*xd[1] + DD[14]*xd[2] + DD[15]*xd[3];
extalk = float(eSeg0 + eSeg1 + eSeg2 + eSeg3);
}
#endif
else {
int nxtalk = 0;
for(int ns = 0; ns < CC.nSG; ns++) {
......@@ -1957,7 +2073,7 @@ void PreprocessingFilterPSA::CoreEnergyFromXtalk(int segMult)
}
}
if(nxtalk < CC.nSG/4) // ??
return; // ??
return;
extalk /= xTalkAveraged*nxtalk;
}
......@@ -1978,4 +2094,240 @@ void PreprocessingFilterPSA::CoreEnergyFromXtalk(int segMult)
}
bool PreprocessingFilterPSA::InvertMatrix2(const double m[4], double invOut[4])
{
double d00 = m[0];
double d01 = m[1];
double d10 = m[2];
double d11 = m[3];
double Det = d00*d11 - d01*d10;
if (Det == 0)
return false;
double D00 = d11; double D01 = -d01;
double D10 = -d10; double D11 = d00;
Det = 1/Det;
invOut[0] = D00 * Det;
invOut[1] = D01 * Det;
invOut[2] = D10 * Det;
invOut[3] = D11 * Det;
#if 0
#define indx(n1, n2) (n1*2+n2)
double unit[4]; memset(unit, 0, sizeof(unit));
for(int i1 = 0; i1 < 2; i1++) {
for(int i2 = 0; i2 < 2; i2++) {
double accu = 0;
for(int ii = 0; ii < 2; ii++) {
accu += m[indx(i1,ii)]*invOut[indx(ii,i2)];
}
unit[indx(i1,i2)] = accu;
}
}
#undef indx
#endif
return true;
}
bool PreprocessingFilterPSA::InvertMatrix3(const double m[9], double invOut[9])
{
double d00 = m[0];
double d01 = m[1];
double d02 = m[2];
double d10 = m[3];
double d11 = m[4];
double d12 = m[5];
double d20 = m[6];
double d21 = m[7];
double d22 = m[8];
double Det = d00*(d11*d22 - d12*d21) + d01*(d12*d20 - d10*d22) + d02*(d10*d21 - d11*d20);
if (Det == 0)
return false;
double D00 = d11*d22 - d12*d21; double D01 = d12*d20 - d10*d22; double D02 = d10*d21 - d11*d20;
double D10 = d21*d02 - d22*d01; double D11 = d22*d00 - d20*d02; double D12 = d20*d01 - d21*d00;
double D20 = d01*d12 - d02*d11; double D21 = d02*d10 - d00*d12; double D22 = d00*d11 - d01*d10;
Det = 1/Det;
invOut[0] = D00 * Det;
invOut[1] = D01 * Det;
invOut[2] = D02 * Det;
invOut[3] = D10 * Det;
invOut[4] = D11 * Det;
invOut[5] = D12 * Det;
invOut[6] = D20 * Det;
invOut[7] = D21 * Det;
invOut[8] = D22 * Det;
#if 0
#define indx(n1, n2) (n1*3+n2)
double unit[9]; memset(unit, 0, sizeof(unit));
for(int i1 = 0; i1 < 3; i1++) {
for(int i2 = 0; i2 < 3; i2++) {
double accu = 0;
for(int ii = 0; ii < 3; ii++) {
accu += m[indx(i1,ii)]*invOut[indx(ii,i2)];
}
unit[indx(i1,i2)] = accu;
}
}
#undef indx
#endif
return true;
}
// this is gluInvertMatrix from Mesa
bool PreprocessingFilterPSA::InvertMatrix4(const double m[16], double invOut[16])
{
double inv[16], det;
int i;
inv[0] = m[5] * m[10] * m[15] -
m[5] * m[11] * m[14] -
m[9] * m[6] * m[15] +
m[9] * m[7] * m[14] +
m[13] * m[6] * m[11] -
m[13] * m[7] * m[10];
inv[4] = -m[4] * m[10] * m[15] +
m[4] * m[11] * m[14] +
m[8] * m[6] * m[15] -
m[8] * m[7] * m[14] -
m[12] * m[6] * m[11] +
m[12] * m[7] * m[10];
inv[8] = m[4] * m[9] * m[15] -
m[4] * m[11] * m[13] -
m[8] * m[5] * m[15] +
m[8] * m[7] * m[13] +
m[12] * m[5] * m[11] -
m[12] * m[7] * m[9];
inv[12] = -m[4] * m[9] * m[14] +
m[4] * m[10] * m[13] +
m[8] * m[5] * m[14] -
m[8] * m[6] * m[13] -
m[12] * m[5] * m[10] +
m[12] * m[6] * m[9];
inv[1] = -m[1] * m[10] * m[15] +
m[1] * m[11] * m[14] +
m[9] * m[2] * m[15] -
m[9] * m[3] * m[14] -
m[13] * m[2] * m[11] +
m[13] * m[3] * m[10];
inv[5] = m[0] * m[10] * m[15] -
m[0] * m[11] * m[14] -
m[8] * m[2] * m[15] +
m[8] * m[3] * m[14] +
m[12] * m[2] * m[11] -
m[12] * m[3] * m[10];
inv[9] = -m[0] * m[9] * m[15] +
m[0] * m[11] * m[13] +
m[8] * m[1] * m[15] -
m[8] * m[3] * m[13] -
m[12] * m[1] * m[11] +
m[12] * m[3] * m[9];
inv[13] = m[0] * m[9] * m[14] -
m[0] * m[10] * m[13] -
m[8] * m[1] * m[14] +
m[8] * m[2] * m[13] +
m[12] * m[1] * m[10] -
m[12] * m[2] * m[9];
inv[2] = m[1] * m[6] * m[15] -
m[1] * m[7] * m[14] -
m[5] * m[2] * m[15] +
m[5] * m[3] * m[14] +
m[13] * m[2] * m[7] -
m[13] * m[3] * m[6];
inv[6] = -m[0] * m[6] * m[15] +
m[0] * m[7] * m[14] +
m[4] * m[2] * m[15] -
m[4] * m[3] * m[14] -
m[12] * m[2] * m[7] +
m[12] * m[3] * m[6];
inv[10] = m[0] * m[5] * m[15] -
m[0] * m[7] * m[13] -
m[4] * m[1] * m[15] +
m[4] * m[3] * m[13] +
m[12] * m[1] * m[7] -
m[12] * m[3] * m[5];
inv[14] = -m[0] * m[5] * m[14] +
m[0] * m[6] * m[13] +
m[4] * m[1] * m[14] -
m[4] * m[2] * m[13] -
m[12] * m[1] * m[6] +
m[12] * m[2] * m[5];
inv[3] = -m[1] * m[6] * m[11] +
m[1] * m[7] * m[10] +
m[5] * m[2] * m[11] -
m[5] * m[3] * m[10] -
m[9] * m[2] * m[7] +
m[9] * m[3] * m[6];
inv[7] = m[0] * m[6] * m[11] -
m[0] * m[7] * m[10] -
m[4] * m[2] * m[11] +
m[4] * m[3] * m[10] +
m[8] * m[2] * m[7] -
m[8] * m[3] * m[6];
inv[11] = -m[0] * m[5] * m[11] +
m[0] * m[7] * m[9] +
m[4] * m[1] * m[11] -
m[4] * m[3] * m[9] -
m[8] * m[1] * m[7] +
m[8] * m[3] * m[5];
inv[15] = m[0] * m[5] * m[10] -
m[0] * m[6] * m[9] -
m[4] * m[1] * m[10] +
m[4] * m[2] * m[9] +
m[8] * m[1] * m[6] -
m[8] * m[2] * m[5];
det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
if (det == 0)
return false;
det = 1.0 / det;
for (i = 0; i < 16; i++)
invOut[i] = inv[i] * det;
#if 0
#define indx(n1, n2) (n1*4+n2)
double unit[16]; memset(unit, 0, sizeof(unit));
for(int i1 = 0; i1 < 4; i1++) {
for(int i2 = 0; i2 < 4; i2++) {
double accu = 0;
for(int ii = 0; ii < 4; ii++) {
accu += m[indx(i1,ii)]*invOut[indx(ii,i2)];
}
unit[indx(i1,i2)] = accu;
}
}
#undef indx
#endif
return true;
}
#endif //PPF_MULTIHIST
......@@ -178,6 +178,9 @@ public:
void calcRiseTime ();
float FindRiseTime (float *data);
void CoreEnergyFromXtalk (int segMult);
bool InvertMatrix2 (const double m[ 4], double invOut[ 4]);
bool InvertMatrix3 (const double m[ 9], double invOut[ 9]);
bool InvertMatrix4 (const double m[16], double invOut[16]);
#endif //PPF_MULTIHIST
};
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment