+//_____________________________________________________________________________
+void AliMUONClusterFinderAZ::AddVirtualPad()
+{
+/// Add virtual pad (with small charge) to improve fit for some
+/// clusters (when pad with max charge is at the extreme of the cluster)
+
+ // Get number of pads in X and Y-directions
+ Int_t nInX = -1, nInY;
+ PadsInXandY(nInX, nInY);
+ //return;
+
+ // Add virtual pad only if number of pads per direction == 2
+ if (nInX != 2 && nInY != 2) return;
+
+ // Find pads with max charge
+ Int_t maxpad[2][2] = {{-1, -1}, {-1, -1}}, cath;
+ Double_t sigmax[2] = {0}, aamax[2] = {0};
+ for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
+ if (fPadIJ[1][j] != 0) continue;
+ cath = fPadIJ[0][j];
+ if (fXyq[2][j] > sigmax[cath]) {
+ maxpad[cath][1] = maxpad[cath][0];
+ aamax[cath] = sigmax[cath];
+ sigmax[cath] = fXyq[2][j];
+ maxpad[cath][0] = j;
+ }
+ }
+ if (maxpad[0][0] >= 0 && maxpad[0][1] < 0 || maxpad[1][0] >= 0 && maxpad[1][1] < 0) {
+ for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
+ if (fPadIJ[1][j] != 0) continue;
+ cath = fPadIJ[0][j];
+ if (j == maxpad[cath][0] || j == maxpad[cath][1]) continue;
+ if (fXyq[2][j] > aamax[cath]) {
+ aamax[cath] = fXyq[2][j];
+ maxpad[cath][1] = j;
+ }
+ }
+ }
+ // Check for mirrors (side X on cathode 0)
+ Bool_t mirror = kFALSE;
+ if (maxpad[0][0] >= 0 && maxpad[1][0] >= 0) {
+ mirror = fXyq[3][maxpad[0][0]] < fXyq[4][maxpad[0][0]];
+ if (!mirror && TMath::Abs(fXyq[3][maxpad[0][0]]-fXyq[3][maxpad[1][0]]) < 0.001) {
+ // Special case when pads on both cathodes have the same size
+ Int_t yud[2] = {0};
+ for (Int_t j = 0; j < fnPads[0]+fnPads[1]; j++) {
+ cath = fPadIJ[0][j];
+ if (j == maxpad[cath][0]) continue;
+ if (fPadIJ[2][j] != fPadIJ[2][maxpad[cath][0]]) continue;
+ if (fPadIJ[3][j] + 1 == fPadIJ[3][maxpad[cath][0]] ||
+ fPadIJ[3][j] - 1 == fPadIJ[3][maxpad[cath][0]]) yud[cath]++;
+ }
+ if (!yud[0]) mirror = kTRUE; // take the other cathode
+ } // if (!mirror &&...
+ } // if (maxpad[0][0] >= 0 && maxpad[1][0] >= 0)
+
+ // Find neughbours of pads with max charges
+ Int_t nn, xList[10], yList[10], ix0, iy0, ix, iy, neighb;
+ for (cath=0; cath<2; cath++) {
+ if (!cath && maxpad[0][0] < 0) continue; // one-sided cluster - cathode 1
+ if (cath && maxpad[1][0] < 0) break; // one-sided cluster - cathode 0
+ if (maxpad[1][0] >= 0) {
+ if (!mirror) {
+ if (!cath && nInY != 2) continue;
+ if (cath && nInX != 2 && (maxpad[0][0] >= 0 || nInY != 2)) continue;
+ } else {
+ if (!cath && nInX != 2) continue;
+ if (cath && nInY != 2 && (maxpad[0][0] >= 0 || nInX != 2)) continue;
+ }
+ }
+
+ Int_t iAddX = 0, iAddY = 0, ix1 = 0, iy1 = 0, iPad = 0;
+ if (maxpad[0][0] < 0) iPad = 1;
+
+ for (iPad=0; iPad<2; iPad++) {
+ if (maxpad[cath][iPad] < 0) continue;
+ if (iPad && !iAddX && !iAddY) break;
+ if (iPad && fXyq[2][maxpad[cath][1]] / sigmax[cath] < 0.5) break;
+
+ Int_t neighbx = 0, neighby = 0;
+ ix0 = fPadIJ[2][maxpad[cath][iPad]];
+ iy0 = fPadIJ[3][maxpad[cath][iPad]];
+ fSegmentation[cath]->Neighbours(ix0, iy0, &nn, xList, yList);
+ Float_t zpad;
+ for (Int_t j=0; j<nn; j++) {
+ if (TMath::Abs(xList[j]-ix0) == 1 || xList[j]*ix0 == -1) neighbx++;
+ if (TMath::Abs(yList[j]-iy0) == 1 || yList[j]*iy0 == -1) neighby++;
+ }
+ if (!mirror) {
+ if (cath) neighb = neighbx;
+ else neighb = neighby;
+ if (maxpad[0][0] < 0) neighb += neighby;
+ else if (maxpad[1][0] < 0) neighb += neighbx;
+ } else {
+ if (!cath) neighb = neighbx;
+ else neighb = neighby;
+ if (maxpad[0][0] < 0) neighb += neighbx;
+ else if (maxpad[1][0] < 0) neighb += neighby;
+ }
+
+ for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
+ if (fPadIJ[0][j] != cath) continue;
+ ix = fPadIJ[2][j];
+ iy = fPadIJ[3][j];
+ if (iy == iy0 && ix == ix0) continue;
+ for (Int_t k=0; k<nn; k++) {
+ if (xList[k] != ix || yList[k] != iy) continue;
+ if (!mirror) {
+ if ((!cath || maxpad[0][0] < 0) &&
+ (TMath::Abs(iy-iy0) == 1 || iy*iy0 == -1)) {
+ if (!iPad && TMath::Abs(ix-ix0) == 1 || ix*ix0 == -1) ix1 = xList[k]; //19-12-05
+ xList[k] = yList[k] = 0;
+ neighb--;
+ break;
+ }
+ if ((cath || maxpad[1][0] < 0) &&
+ (TMath::Abs(ix-ix0) == 1 || ix*ix0 == -1)) {
+ if (!iPad) ix1 = xList[k]; //19-12-05
+ xList[k] = yList[k] = 0;
+ neighb--;
+ }
+ } else {
+ if ((!cath || maxpad[0][0] < 0) &&
+ (TMath::Abs(ix-ix0) == 1 || ix*ix0 == -1)) {
+ if (!iPad) ix1 = xList[k]; //19-12-05
+ xList[k] = yList[k] = 0;
+ neighb--;
+ break;
+ }
+ if ((cath || maxpad[1][0] < 0) &&
+ (TMath::Abs(iy-iy0) == 1 || iy*iy0 == -1)) {
+ xList[k] = yList[k] = 0;
+ neighb--;
+ }
+ }
+ break;
+ } // for (Int_t k=0; k<nn;
+ if (!neighb) break;
+ } // for (Int_t j=0; j<fnPads[0]+fnPads[1];
+ if (!neighb) continue;
+
+ // Add virtual pad
+ Int_t npads, isec;
+ isec = 0;
+ for (Int_t j=0; j<nn; j++) {
+ if (xList[j] == 0 && yList[j] == 0) continue;
+ npads = fnPads[0] + fnPads[1];
+ fPadIJ[0][npads] = cath;
+ fPadIJ[1][npads] = 0;
+ ix = xList[j];
+ iy = yList[j];
+ if (TMath::Abs(ix-ix0) == 1 || ix*ix0 == -1) {
+ if (iy != iy0) continue; // new segmentation - check
+ if (nInX != 2) continue; // new
+ if (!mirror) {
+ if (!cath && maxpad[1][0] >= 0) continue;
+ } else {
+ if (cath && maxpad[0][0] >= 0) continue;
+ }
+ if (iPad && !iAddX) continue;
+ fSegmentation[cath]->GetPadC(ix, iy, fXyq[0][npads], fXyq[1][npads], zpad);
+ if (fXyq[0][npads] > 1.e+5) continue; // temporary fix
+ if (ix == ix1) continue; //19-12-05
+ if (ix1 == ix0) continue;
+ if (maxpad[1][0] < 0 || mirror && maxpad[0][0] >= 0) {
+ if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[0]/100, 5.);
+ else fXyq[2][npads] = TMath::Min (aamax[0]/100, 5.);
+ }
+ else {
+ if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[1]/100, 5.);
+ else fXyq[2][npads] = TMath::Min (aamax[1]/100, 5.);
+ }
+ fXyq[2][npads] = TMath::Max (fXyq[2][npads], (float)1);
+ fXyq[3][npads] = -2; // flag
+ fPadIJ[2][npads] = ix;
+ fPadIJ[3][npads] = iy;
+ fnPads[1]++;
+ iAddX = npads;
+ if (fDebug) printf(" ***** Add virtual pad in X ***** %f %f %f %3d %3d \n", fXyq[2][npads],
+ fXyq[0][npads], fXyq[1][npads], ix, iy);
+ ix1 = ix0;
+ continue;
+ }
+ if (nInY != 2) continue;
+ if (!mirror && cath && maxpad[0][0] >= 0) continue;
+ if (mirror && !cath && maxpad[1][0] >= 0) continue;
+ if (TMath::Abs(iy-iy0) == 1 || TMath::Abs(iy*iy0) == 1) {
+ if (ix != ix0) continue; // new segmentation - check
+ if (iPad && !iAddY) continue;
+ fSegmentation[cath]->GetPadC(ix, iy, fXyq[0][npads], fXyq[1][npads], zpad);
+ if (iy1 == iy0) continue;
+ //if (iPad && iy1 == iy0) continue;
+ if (maxpad[0][0] < 0 || mirror && maxpad[1][0] >= 0) {
+ if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[1]/15, fgkZeroSuppression);
+ else fXyq[2][npads] = TMath::Min (aamax[1]/15, fgkZeroSuppression);
+ }
+ else {
+ if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[0]/15, fgkZeroSuppression);
+ else fXyq[2][npads] = TMath::Min (aamax[0]/15, fgkZeroSuppression);
+ }
+ fXyq[2][npads] = TMath::Max (fXyq[2][npads], (float)1);
+ fXyq[3][npads] = -2; // flag
+ fPadIJ[2][npads] = ix;
+ fPadIJ[3][npads] = iy;
+ fnPads[1]++;
+ iAddY = npads;
+ if (fDebug) printf(" ***** Add virtual pad in Y ***** %f %f %f %3d %3d \n", fXyq[2][npads],
+ fXyq[0][npads], fXyq[1][npads], ix, iy);
+ iy1 = iy0;
+ }
+ } // for (Int_t j=0; j<nn;
+ } // for (Int_t iPad=0;
+ } // for (cath=0; cath<2;
+ return;
+}
+
+//_____________________________________________________________________________
+void AliMUONClusterFinderAZ::PadsInXandY(Int_t &nInX, Int_t &nInY)
+{
+/// Find number of pads in X and Y-directions (excluding virtual ones and
+/// overflows)
+
+ static Int_t nXsaved = 0, nYsaved = 0;
+ nXsaved = nYsaved = 0;
+ //if (nInX >= 0) {nInX = nXsaved; nInY = nYsaved; return; }
+ Float_t *xPad0 = NULL, *yPad0 = NULL, *xPad1 = NULL, *yPad1 = NULL;
+ Float_t wMinX[2] = {99, 99}, wMinY[2] = {99, 99};
+ Int_t *nPad0 = NULL, *nPad1 = NULL;
+ Int_t nPads = fnPads[0] + fnPads[1];
+ if (fnPads[0]) {
+ xPad0 = new Float_t[nPads];
+ yPad0 = new Float_t[nPads];
+ nPad0 = new Int_t[nPads];
+ }
+ if (fnPads[1]) {
+ xPad1 = new Float_t[nPads];
+ yPad1 = new Float_t[nPads];
+ nPad1 = new Int_t[nPads];
+ }
+ Int_t n0 = 0, n1 = 0, cath, npadx[2] = {1, 1}, npady[2] = {1, 1};
+ for (Int_t j = 0; j < nPads; j++) {
+ if (nInX < 0 && fPadIJ[1][j] != 0) continue; // before fit
+ else if (nInX == 0 && fPadIJ[1][j] != 1) continue; // fit - exclude overflows
+ else if (nInX > 0 && fPadIJ[1][j] != 1 && fPadIJ[1][j] != -9) continue; // exclude non-marked
+ if (nInX <= 0 && fXyq[2][j] > fgkSaturation-1) continue; // skip overflows
+ cath = fPadIJ[0][j];
+ if (fXyq[3][j] > 0) { // exclude virtual pads
+ wMinX[cath] = TMath::Min (wMinX[cath], fXyq[3][j]);
+ wMinY[cath] = TMath::Min (wMinY[cath], fXyq[4][j]);
+ //20-12-05 }
+ if (cath) { xPad1[n1] = fXyq[0][j]; yPad1[n1++] = fXyq[1][j]; }
+ else { xPad0[n0] = fXyq[0][j]; yPad0[n0++] = fXyq[1][j]; }
+ }
+ }
+
+ // Sort
+ if (n0) {
+ TMath::Sort (n0, xPad0, nPad0); // in X
+ for (Int_t i = 1; i < n0; i++)
+ if (xPad0[nPad0[i]] - xPad0[nPad0[i-1]] < -0.01) npadx[0]++;
+ TMath::Sort (n0, yPad0, nPad0); // in Y
+ for (Int_t i = 1; i < n0; i++)
+ if (yPad0[nPad0[i]] - yPad0[nPad0[i-1]] < -0.01) npady[0]++;
+ }
+
+ if (n1) {
+ TMath::Sort (n1, xPad1, nPad1); // in X
+ for (Int_t i = 1; i < n1; i++)
+ if (xPad1[nPad1[i]] - xPad1[nPad1[i-1]] < -0.01) npadx[1]++;
+ TMath::Sort (n1, yPad1, nPad1); // in Y
+ for (Int_t i = 1; i < n1; i++)
+ if (yPad1[nPad1[i]] - yPad1[nPad1[i-1]] < -0.01) npady[1]++;
+ }
+ if (fnPads[0]) { delete [] xPad0; delete [] yPad0; delete [] nPad0; }
+ if (fnPads[1]) { delete [] xPad1; delete [] yPad1; delete [] nPad1; }
+ if (TMath::Abs (wMinY[0] - wMinY[1]) < 1.e-3) nInY = TMath::Max (npady[0], npady[1]);
+ else nInY = wMinY[0] < wMinY[1] ? npady[0] : npady[1];
+ if (TMath::Abs (wMinX[0] - wMinX[1]) < 1.e-3) nInX = TMath::Max (npadx[0], npadx[1]);
+ else nInX = wMinX[0] < wMinX[1] ? npadx[0] : npadx[1];
+}
+
+//_____________________________________________________________________________
+void AliMUONClusterFinderAZ::Simple()
+{
+/// Process simple cluster (small number of pads) without EM-procedure
+
+ Int_t nForFit = 1, clustFit[1] = {0}, nfit;
+ Double_t parOk[3] = {0.};
+ TObjArray *clusters[1];
+ clusters[0] = fPixArray;
+ for (Int_t i = 0; i < fnPads[0]+fnPads[1]; i++) {
+ if (fXyq[2][i] > fgkSaturation-1) fPadIJ[1][i] = -9;
+ else fPadIJ[1][i] = 1;
+ }
+ nfit = Fit(1, nForFit, clustFit, clusters, parOk);
+}
+
+//_____________________________________________________________________________
+void AliMUONClusterFinderAZ::Errors(AliMUONRawCluster *clus)
+{
+/// Correct reconstructed coordinates for some clusters and evaluate errors
+
+ Double_t qTot = clus->GetCharge(0), fmin = clus->GetChi2(0);
+ Double_t xreco = clus->GetX(0), yreco = clus->GetY(0), zreco = clus->GetZ(0);
+ Double_t sigmax[2] = {0};
+
+ Int_t nInX = 1, nInY, maxdig[2] ={-1, -1}, digit, cath1, isec;
+ PadsInXandY(nInX, nInY);
+
+ // Find pad with maximum signal
+ for (Int_t cath = 0; cath < 2; cath++) {
+ for (Int_t j = 0; j < clus->GetMultiplicity(cath); j++) {
+ cath1 = cath;
+ digit = clus->GetIndex(j, cath);
+ if (digit < 0) { cath1 = TMath::Even(cath); digit = -digit - 1; } // from the other cathode
+
+ if (clus->GetContrib(j,cath) > sigmax[cath1]) {
+ sigmax[cath1] = clus->GetContrib(j,cath);
+ maxdig[cath1] = digit;
+ }
+ }
+ }
+
+ // Size of pad with maximum signal and reco coordinate distance from the pad center
+ AliMUONDigit *mdig = 0;
+ Double_t wx[2], wy[2], dxc[2], dyc[2];
+ Float_t xpad, ypad, zpad;
+ Int_t ix, iy;
+ for (Int_t cath = 0; cath < 2; cath++) {
+ if (maxdig[cath] < 0) continue;
+ mdig = fInput->Digit(cath,maxdig[cath]);
+ isec = fSegmentation[cath]->Sector(mdig->PadX(), mdig->PadY());
+ wx[cath] = fSegmentation[cath]->Dpx(isec);
+ wy[cath] = fSegmentation[cath]->Dpy(isec);
+ fSegmentation[cath]->GetPadI(xreco, yreco, zreco, ix, iy);
+ isec = fSegmentation[cath]->Sector(ix, iy);
+ if (isec > 0) {
+ fSegmentation[cath]->GetPadC(ix, iy, xpad, ypad, zpad);
+ dxc[cath] = xreco - xpad;
+ dyc[cath] = yreco - ypad;
+ }
+ }
+
+ // Check if pad with max charge at the edge (number of neughbours)
+ Int_t nn, xList[10], yList[10], neighbx[2][2] = {{0,0}, {0,0}}, neighby[2][2]= {{0,0}, {0,0}};
+ for (Int_t cath = 0; cath < 2; cath++) {
+ if (maxdig[cath] < 0) continue;
+ mdig = fInput->Digit(cath,maxdig[cath]);
+ fSegmentation[cath]->Neighbours(mdig->PadX(), mdig->PadY(), &nn, xList, yList);
+ isec = fSegmentation[cath]->Sector(mdig->PadX(), mdig->PadY());
+ /*??
+ Float_t sprX = fResponse->SigmaIntegration() * fResponse->ChargeSpreadX();
+ Float_t sprY = fResponse->SigmaIntegration() * fResponse->ChargeSpreadY();
+ //fSegmentation[cath]->FirstPad(fInput->DetElemId(),muons[ihit][1], muons[ihit][2], muons[ihit][3], sprX, sprY);
+ //fSegmentation[cath]->FirstPad(fInput->DetElemId(),xreco, yreco, zreco, sprX, sprY);
+ fSegmentation[cath]->FirstPad(xreco, yreco, zreco, sprX, sprY);
+ Int_t border = 0;
+ //if (fSegmentation[cath]->Sector(fInput->DetElemId(),fSegmentation[cath]->Ix(),fSegmentation[cath]->Iy()) <= 0) {
+ if (fSegmentation[cath]->Sector(fSegmentation[cath]->Ix(), fSegmentation[cath]->Iy()) <= 0) {
+ //fSegmentation[cath]->NextPad(fInput->DetElemId());
+ fSegmentation[cath]->NextPad();
+ border = 1;
+ }
+ */
+ for (Int_t j=0; j<nn; j++) {
+ //if (border && yList[j] < fSegmentation[cath]->Iy()) continue;
+ fSegmentation[cath]->GetPadC(xList[j], yList[j], xpad, ypad, zpad);
+ //cout << ch << " " << xList[j] << " " << yList[j] << " " << border << " " << x << " " << y << " " << xpad << " " << ypad << endl;
+ if (TMath::Abs(xpad) < 1 && TMath::Abs(ypad) < 1) continue;
+ if (xList[j] == mdig->PadX()-1 || mdig->PadX() == 1 &&
+ xList[j] == -1) neighbx[cath][0] = 1;
+ else if (xList[j] == mdig->PadX()+1 || mdig->PadX() == -1 &&
+ xList[j] == 1) neighbx[cath][1] = 1;
+ if (yList[j] == mdig->PadY()-1 || mdig->PadY() == 1 &&
+ yList[j] == -1) neighby[cath][0] = 1;
+ else if (yList[j] == mdig->PadY()+1 || mdig->PadY() == -1 &&
+ yList[j] == 1) neighby[cath][1] = 1;
+ } // for (Int_t j=0; j<nn;
+ if (neighbx[cath][0] && neighbx[cath][1]) neighbx[cath][0] = 0;
+ else if (neighbx[cath][1]) neighbx[cath][0] = -1;
+ else neighbx[cath][0] = 1;
+ if (neighby[cath][0] && neighby[cath][1]) neighby[cath][0] = 0;
+ else if (neighby[cath][1]) neighby[cath][0] = -1;
+ else neighby[cath][0] = 1;
+ }
+
+ Int_t iOver = clus->GetClusterType();
+ // One-sided cluster
+ if (!clus->GetMultiplicity(0)) {
+ neighby[0][0] = neighby[1][0];
+ wy[0] = wy[1];
+ if (iOver < 99) iOver += 100 * iOver;
+ dyc[0] = dyc[1];
+ } else if (!clus->GetMultiplicity(1)) {
+ neighbx[1][0] = neighbx[0][0];
+ wx[1] = wx[0];
+ if (iOver < 99) iOver += 100 * iOver;
+ dxc[1] = dxc[0];
+ }
+
+ // Apply corrections and evaluate errors
+ Double_t errY, errX;
+ Errors(nInY, nInX, neighby[0][0],neighbx[1][0], fmin, wy[0]*10, wx[1]*10, iOver,
+ dyc[0], dxc[1], qTot, yreco, xreco, errY, errX);
+ errY = TMath::Max (errY, 0.01);
+ //errY = 0.01;
+ //errX = TMath::Max (errX, 0.144);
+ clus->SetX(0, xreco); clus->SetY(0, yreco);
+ clus->SetErrX(errX); clus->SetErrY(errY);
+}
+
+//_____________________________________________________________________________
+void AliMUONClusterFinderAZ::Errors(Int_t ny, Int_t nx, Int_t iby, Int_t ibx, Double_t fmin,
+ Double_t wy, Double_t wx, Int_t iover,
+ Double_t dyc, Double_t /*dxc*/, Double_t qtot,
+ Double_t &yrec, Double_t &xrec, Double_t &erry, Double_t &errx)
+{
+/// Correct reconstructed coordinates for some clusters and evaluate errors
+
+ erry = 0.01;
+ errx = 0.144;
+ Int_t iovery = iover % 100;
+ Double_t corr = 0;
+
+/* ---> Ny = 1 */
+ if (ny == 1) {
+ if (iby != 0) {
+ // edge effect
+ yrec += iby * (0.1823+0.2008)/2;
+ erry = 0.04587;
+ } else {
+ // Find "effective pad width"
+ Double_t width = 0.218 / (1.31e-4 * TMath::Exp (2.688 * TMath::Log(qtot)) + 1) * 2;
+ width = TMath::Min (width, 0.4);
+ erry = width / TMath::Sqrt(12.);
+ erry = TMath::Max (erry, 0.01293);
+ }
+ goto x; //return;
+ }
+
+/* ---> "Bad" fit */
+ if (fmin > 0.4) {
+ erry = 0.1556;
+ if (ny == 5) erry = 0.06481;
+ goto x; //return;
+ }
+
+/* ---> By != 0 */
+ if (iby != 0) {
+ if (ny > 2) {
+ erry = 0.00417; //0.01010
+ } else {
+ // ny = 2
+ if (dyc * iby > -0.05) {
+ Double_t dyc2 = dyc * dyc;
+ if (iby < 0) {
+ corr = 0.019 - 0.602 * dyc + 8.739 * dyc2 - 44.209 * dyc2 * dyc;
+ corr = TMath::Min (corr, TMath::Abs(-0.25-dyc));
+ yrec -= corr;
+ //dyc -= corr;
+ erry = 0.00814;
+ } else {
+ corr = 0.006 + 0.300 * dyc + 6.147 * dyc2 + 42.039 * dyc2 * dyc;
+ corr = TMath::Min (corr, 0.25-dyc);
+ yrec += corr;
+ //dyc += corr;
+ erry = 0.01582;
+ }
+ } else {
+ erry = (0.00303 + 0.00296) / 2;
+ }
+ }
+ goto x; //return;
+ }
+
+/* ---> Overflows */
+ if (iovery != 0) {
+ if (qtot < 3000) {
+ erry = 0.0671;
+ } else {
+ if (iovery > 1) {
+ erry = 0.09214;
+ } else if (TMath::Abs(wy - 5) < 0.1) {
+ erry = 0.061; //0.06622
+ } else {
+ erry = 0.00812; // 0.01073
+ }
+ }
+ goto x; //return;
+ }
+
+/* ---> "Good" but very high signal */
+ if (qtot > 4000) {
+ if (TMath::Abs(wy - 4) < 0.1) {
+ erry = 0.00117;
+ } else if (fmin < 0.03 && qtot < 6000) {
+ erry = 0.01003;
+ } else {
+ erry = 0.1931;
+ }
+ goto x; //return;
+ }
+
+/* ---> "Good" clusters */
+ if (ny > 3) {
+ if (TMath::Abs(wy - 5) < 0.1) {
+ erry = 0.0011; //0.00304
+ } else if (qtot < 400.) {
+ erry = 0.0165;
+ } else {
+ erry = 0.00135; // 0.00358
+ }
+ } else if (ny == 3) {
+ if (TMath::Abs(wy - 4) < 0.1) {
+ erry = 35.407 / (1 + TMath::Exp(5.511*TMath::Log(qtot/265.51))) + 11.564;
+ //erry = 83.512 / (1 + TMath::Exp(3.344*TMath::Log(qtot/211.58))) + 12.260;
+ } else {
+ erry = 147.03 / (1 + TMath::Exp(1.713*TMath::Log(qtot/73.151))) + 9.575;
+ //erry = 91.743 / (1 + TMath::Exp(2.332*TMath::Log(qtot/151.67))) + 11.453;
+ }
+ erry *= 1.e-4;
+ } else {
+ // ny = 2
+ if (TMath::Abs(wy - 4) < 0.1) {
+ erry = 60.800 / (1 + TMath::Exp(3.305*TMath::Log(qtot/104.53))) + 11.702;
+ //erry = 73.128 / (1 + TMath::Exp(5.676*TMath::Log(qtot/120.93))) + 17.839;
+ } else {
+ erry = 117.98 / (1 + TMath::Exp(2.005*TMath::Log(qtot/37.649))) + 21.431;
+ //erry = 99.066 / (1 + TMath::Exp(4.900*TMath::Log(qtot/107.57))) + 25.315;
+ }
+ erry *= 1.e-4;
+ }
+ //return;
+
+ x:
+/* ---> X-coordinate */
+/* ---> Y-side */
+ if (wx > 11) {
+ errx = 0.0036;
+ xrec -= 0.1385;
+ return;
+ }
+/* ---> Nx = 1 */
+ if (nx == 1) {
+ if (TMath::Abs(wx - 6) < 0.1) {
+ if (qtot < 40) errx = 0.1693;
+ else errx = 0.06241;
+ } else if (TMath::Abs(wx - 7.5) < 0.1) {
+ if (qtot < 40) errx = 0.2173;
+ else errx = 0.07703;
+ } else if (TMath::Abs(wx - 10) < 0.1) {
+ if (ibx == 0) {
+ if (qtot < 40) errx = 0.2316;
+ else errx = 0.1426;
+ } else {
+ xrec += (0.2115 + 0.1942) / 2 * ibx;
+ errx = 0.1921;
+ }
+ }
+ return;
+ }
+/* ---> "Bad" fit */
+ if (fmin > 0.5) {
+ errx = 0.1591;
+ return;
+ }
+/* ---> Bx != 0 */
+ if (ibx != 0) {
+ if (ibx > 0) { errx = 0.06761; xrec -= 0.03832; }
+ else { errx = 0.06653; xrec += 0.02581; }
+ return;
+ }
+/* ---> Overflows */
+ if (iover != 0) {
+ if (TMath::Abs(wx - 6) < 0.1) errx = 0.06979;
+ else if (TMath::Abs(wx - 7.5) < 0.1) errx = 0.1089;
+ else if (TMath::Abs(wx - 10) < 0.1) errx = 0.09847;
+ return;
+ }
+/* ---> Good */
+ if (TMath::Abs(wx - 6) < 0.1) errx = 0.06022;
+ else if (TMath::Abs(wx - 7.5) < 0.1) errx = 0.07247;
+ else if (TMath::Abs(wx - 10) < 0.1) errx = 0.07359;
+}
+