#include "CbmRichRonchiAna.h"

#include <boost/gil/gil_all.hpp>
#include <boost/gil/extension/io/tiff_dynamic_io.hpp>
//#include <boost/gil/extension/io/png_dynamic_io.hpp>

#include <iostream>
#include "TH2D.h"
#include "TCanvas.h"
#include "CbmDrawHist.h"
#include "TEllipse.h"
#include <math.h>

using namespace boost::gil;
using namespace std;

// For Ubuntu one needs dev version of libtiff
// sudo apt-get install libtiff-dev

CbmRichRonchiAna::CbmRichRonchiAna():       // constructor

// constant values
fPi(3.141592654),
fRadiusMirror(3000000),                 // in microns   !!! MIGHT HAVE TO BE CHANGED TO DISTANCE_MIRROR-POINTSOURCE IF IT IS NOT IN THE CENTER OF CURVATURE
fEdgeLengthCCD(13300),                  // in microns
fEdgeLengthPixel(13),                   // in microns
fPitchGrid(200),                        // in microns
fImageWidth(1024),                      // in pixels
 
// values to be measured first
fDistRulingCCD(18200),                  // in microns
fDistMirrorRuling(2900000),             // in microns
fEffectiveEdgeLengthGridPlate(10000),   // in microns
fEffectiveEdgeLengthMirror(380000),     // in microns
fEffectiveEdgeLengthCCD(38*100),        // in microns
fLineDistance(20),                      // in pixels
fCenterX(500),                          // in pixels
fCenterY(480)                           // in pixels
{

}

CbmRichRonchiAna::~CbmRichRonchiAna()
{

}

void CbmRichRonchiAna::Run()
{
    // Initialization
    vector<vector<int> > dataV;
    vector<vector<int> > dataH;
        
    if ( fTiffFileNameV == "" || fTiffFileNameH == "" ) {
        Fatal("CbmRichRonchiAna::Run:", "No FileNameV or FileNameH!");
    }
    else { cout << "FileNameV: " << fTiffFileNameV << " and FileNameH: " << fTiffFileNameH << endl;
        dataV = ReadTiffFile(fTiffFileNameV);
        dataH = ReadTiffFile(fTiffFileNameH);
    }

    int width = dataV.size();
    int height = dataV[0].size();
    
    // initialisierung der histogramme: name, groesse usw
    TH2D* hInitH = new TH2D("hInitH", "hInitH;X [pixel];Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hMeanIntensityH = new TH2D("hMeanIntensityH", "hMeanIntensityH;X [pixel];Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hPeakH = new TH2D("hPeakH", "hPeakH;X [pixel];Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hMeanYH = new TH2D("hMeanYH", "hMeanYH;X [pixel];Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hEjectingH = new TH2D("hEjectingH", "hEjectingH; X [pixel]; Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hGapsFilledH = new TH2D("hGapsFilledH", "hGapsFilledH; X [pixel]; Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);

    TH2D* hInitV = new TH2D("hInitV", "hInitV;X [pixel];Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hMeanIntensityV = new TH2D("hMeanIntensityV", "hMeanIntensityV;X [pixel];Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hPeakV = new TH2D("hPeakV", "hPeakV;X [pixel];Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hMeanYV = new TH2D("hMeanYV", "hMeanYV;X [pixel];Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    
    TH2D* hSuperpose = new TH2D("hSuperpose", "hSuperpose;X [pixel];Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    
    TH2D* hIntersections = new TH2D("hIntersections", "hIntersections;X [pixel];Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hMatrix = new TH2D("hMatrix", "hIntersections;X [pixel];Y [pixel];Intensity", width, -.5, width - 0.5, height, -0.5, height - 0.5);
    
    TH2D* hNormals = new TH2D("normals", "single heights from normals; X [pixel]; Y [pixel];Single Heights", width, -0.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hHeightY = new TH2D("heightY", "heightY; X [pixel];Y [pixel]; summed Heights in Y", width, -0.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hHeights = new TH2D("heights", "summed heights; X [pixel];Y [pixel]; summed Heights", width, -0.5, width - 0.5, height, -0.5, height - 0.5);
    TH2D* hAvrdHeights = new TH2D("averaged heights", "averaged heights; X [pixel]; Y [pixel]; averaged Heights", width, -0.5, width - 0.5, height, -0.5, height - 0.5);


    // vertical image
    FillH2WithVector(hInitV, dataV);
    DoRotation(dataV);
    DoMeanIntensityY(dataV);
    FillH2WithVector(hMeanIntensityV, dataV);
    DoPeakFinderY(dataV);
    FillH2WithVector(hPeakV, dataV);
    DoMeanY(dataV);
    DoEjectSingle(dataV);
    DoFillGaps(dataV);
    DoRotation(dataV);
    FillH2WithVector(hMeanYV, dataV);
    DoShadowCenter(dataV);

    // horizontal image
    FillH2WithVector(hInitH, dataH);
    DoMeanIntensityY(dataH);
    FillH2WithVector(hMeanIntensityH, dataH);
    DoPeakFinderY(dataH);
    FillH2WithVector(hPeakH, dataH);
    DoMeanY(dataH);
    DoEjectSingle(dataH);
    FillH2WithVector(hMeanYH, dataH);
    DoFillGaps(dataH);
    FillH2WithVector(hGapsFilledH, dataH);
    DoShadowCenter(dataH);
    
    // finding intersections
    vector<vector<int> > intersectionXY = DoIntersection(dataH, dataV);
    vector<vector<int> > dataSup = DoSuperpose(dataH, dataV);
    FillH2WithVector(hSuperpose, dataSup);
    
    // numbering the intersections and calculating 
    vector<vector<int> >    matrix = DoNumberIntersections(intersectionXY);
    vector<vector<double> > normal = DoLocalNormal(matrix);
    vector<vector<double> > sphere = DoHeight(normal);
    DoSphere(sphere);
    DoIntegrate(sphere);
    
    {
        TCanvas* c = new TCanvas("ronchi_2d_normals", "ronchi_2d_normals", 1000, 1000);
                
        int nH = sphere.size();
        int radius = 10;
        
        for (int l = 0; l < nH; l++){
            for (int x = sphere[l][0]-radius; x <= sphere[l][0]+radius; x++) {
                for (int y = sphere[l][1]-radius; y <= sphere[l][1]+radius; y++) {
                    hNormals->SetBinContent(x, y, sphere[l][10]);
                }
            }
            
        }
        
        DrawH2(hNormals); 
    }
    
    {
        TCanvas* c = new TCanvas("ronchi_2d_heightY", "ronchi_2d_heightY", 1000, 1000);
                
        int nH = sphere.size();
        int radius = 8;
        
        for (int l = 0; l < nH; l++){
            for (int x = sphere[l][0]-radius; x <= sphere[l][0]+radius; x++) {
                for (int y = sphere[l][1]-radius; y <= sphere[l][1]+radius; y++) {
                    hHeightY->SetBinContent(x, y, sphere[l][11]);
                }
            }
            
        }
        
        DrawH2(hHeightY); 
    }
     
    {
        TCanvas* c = new TCanvas("ronchi_2d_summed_heights", "ronchi_2d_summed_heights", 1000, 1000);
                
        int nS = sphere.size();
        int radius = 10;    
        
        for (int l = 0; l < nS; l++){
            for (int x = sphere[l][0]-radius; x <= sphere[l][0]+radius; x++) {
                for (int y = sphere[l][1]-radius; y <= sphere[l][1]+radius; y++) {
                    hHeights->SetBinContent(x, y, sphere[l][14]);
                }
            }
            
        }
        
        DrawH2(hHeights); 
                 
    } 
    
    {
        TCanvas* c = new TCanvas("ronchi_2d_avrd_heights", "ronchi_2d_avrd_heights", 1000, 1000);
                
        int nS = sphere.size();
        int radius = 10;    
        
        for (int l = 0; l < nS; l++){
            for (int x = sphere[l][0]-radius; x <= sphere[l][0]+radius; x++) {
                for (int y = sphere[l][1]-radius; y <= sphere[l][1]+radius; y++) {
                    hAvrdHeights->SetBinContent(x, y, sphere[l][19]);
                }
            }
            
        }
        
        DrawH2(hAvrdHeights); 
                 
    } 
    
    {
        TCanvas* c = new TCanvas("ronchi_2d_horizontal", "ronchi_2d_horizontal", 2500, 500);
        c->Divide(5,1);
        c->cd(1);
        DrawH2(hInitH);
        c->cd(2);
        DrawH2(hMeanIntensityH);
        c->cd(3);
        DrawH2(hPeakH);
        c->cd(4);
        DrawH2(hMeanYH);
        c->cd(5);
        DrawH2(hGapsFilledH);
    }

    {
        TCanvas* c = new TCanvas("ronchi_2d_vertical", "ronchi_2d_vertical", 2000, 500);
        c->Divide(4,1);
        c->cd(1);
        DrawH2(hInitV);
        c->cd(2);
        DrawH2(hMeanIntensityV);
        c->cd(3);
        DrawH2(hPeakV);
        c->cd(4);
        DrawH2(hMeanYV);
    }
    
    {
        TCanvas* c = new TCanvas("Superpose", "Superpose", 1000, 1000);
        DrawH2(hSuperpose);
    }
    
    /*{
        TH1D* h1 = hInitH->ProjectionY("_py1", 100, 100);
        TH1D* h2 = hInitH->ProjectionY("_py2", 200, 200);
        TH1D* h3 = hInitH->ProjectionY("_py3", 300, 300);
        TH1D* h4 = hInitH->ProjectionY("_py4", 500, 500);

        TH1D* hM1 = hMeanIntensityH->ProjectionY("_pyM1", 100, 100);
        TH1D* hM2 = hMeanIntensityH->ProjectionY("_pyM2", 200, 200);
        TH1D* hM3 = hMeanIntensityH->ProjectionY("_pyM3", 300, 300);
        TH1D* hM4 = hMeanIntensityH->ProjectionY("_pyM4", 400, 400);

        TH1D* hP1 = hPeakH->ProjectionY("_pyP1", 100, 100);
        TH1D* hP2 = hPeakH->ProjectionY("_pyP2", 200, 200);
        TH1D* hP3 = hPeakH->ProjectionY("_pyP3", 300, 300);
        TH1D* hP4 = hPeakH->ProjectionY("_pyP4", 400, 400);

        TCanvas* c2 = new TCanvas("ronchi_1d_slices_horizontal", "ronchi_1d_slices_horizontal", 1000, 1000);
        c2->Divide(2,2);
        c2->cd(1);
        DrawH1({h1,hM1,hP1}, {"Init", "Mean", "Peak"});
        c2->cd(2);
        DrawH1({h2,hM2,hP2}, {"Init", "Mean", "Peak"});
        c2->cd(3);
        DrawH1({h3,hM3,hP3}, {"Init", "Mean", "Peak"});
        c2->cd(4);
        DrawH1({h4,hM4,hP4}, {"Init", "Mean", "Peak"});
    }
    
    {
        TCanvas* c = new TCanvas("ronchi_2d_intersection", "ronchi_2d_intersection", 1000, 1000);
        DrawH2(hSuperpose);
        
        int nI = intersectionXY.size();
        
        for (int i = 0; i < nI; i++) {
            if (intersectionXY[i][0]==0) continue;
            TEllipse* center = new TEllipse(intersectionXY[i][0], intersectionXY[i][1], 3);
            center->Draw();
        }
    }*/

   
   {
        TCanvas* c = new TCanvas("ronchi_2d_matrix", "ronchi_2d_matrixX", 1000, 1000);
        DrawH2(hMatrix);
        
        int nM = matrix.size();
        
        for (int i = 0; i < nM; i++) {
            TEllipse* center = new TEllipse(matrix[i][0], matrix[i][1], 3);
            center->Draw();
        }
    }
    
    {
        TCanvas* c = new TCanvas("ronchi_ejecting", "ronchi_ejecting", 1000, 1000);
                
        int nD = dataH.size();
        
        for (int x = 0; x < nD; x++) {
            for (int y = 0; y < nD; y++) {
                hEjectingH->SetBinContent(x, y, dataH[x][y]);
            }
        }
            
        DrawH2(hEjectingH); 
    }
}

// VECTORS:
// intXY:  [x,y,index1,index2] 
// matrix: [x,y,i,j,x_init,y_init,i_init,j_init]
// normal: [x,y,i,j,normalX,normalY,xM,yM]
// sphere: [x,y,i,j,normalX,normalY,upEdge,lowEdge,lEdge,rEdge,hCenterSingle, heightY, totHeight,index1, upper, lower, left, right, averageDev, hCenterSummed, index2, index3]
//         [0 1 2 3   4        5      6       7      8     9        10          11         12       13     14     15    16     17       18            19         20      21  ]
   

void CbmRichRonchiAna::DoIntegrate(vector<vector<double> >& sphere)
{
    int nS = sphere.size();
    int iFirst;             // data of segment that is closest to the center
    int jFirst; 
    int lFirst;
    int xMin, xMax, yMin, yMax;
    int radius    = 1;     
    int delta     = 0;      // integrated height differences of neighboured segments
    int threshold = 100;    // threshold for delta
    double avDiff;
    
    // to find the intersection that is closest to the center
    for (int l = 0; l < nS; l++) {      
        xMin = fCenterX - radius;
        xMax = fCenterX + radius;
        yMin = fCenterY - radius;
        yMax = fCenterY + radius;
        if (sphere[l][0] >= xMin && sphere[l][0] <= xMax && sphere[l][1] >= yMin && sphere[l][1] <= yMax) {
            iFirst = sphere[l][2];
            jFirst = sphere[l][3];
            lFirst = l;
            sphere[l][21] = 1;
            cout << endl << "Found closest: sphere[" << l << "] = [" << sphere[l][0] << "," << sphere[l][1] << "," << sphere[l][2] << "," << sphere[l][3] << ", ...], radius = " << radius << endl;
            break;
        }
        if (l == nS-1) {
            radius++;
            l = -1;
        }
    }
    
    sphere[lFirst][14] = sphere[lFirst][6];
    sphere[lFirst][15] = sphere[lFirst][7];
    sphere[lFirst][16] = sphere[lFirst][8];
    sphere[lFirst][17] = sphere[lFirst][9];
    sphere[lFirst][19] = sphere[lFirst][10];
        
    // to know range of i and j to terminate loops
    int iMin = 1000;
    int iMax = -1000;
    int jMin = 1000;
    int jMax = -1000;
    for (int l = 0; l < nS; l++) {  
        if (sphere[l][2] < iMin) iMin = sphere[l][2];
        if (sphere[l][2] > iMax) iMax = sphere[l][2];
        if (sphere[l][3] < jMin) jMin = sphere[l][3];
        if (sphere[l][3] > jMax) jMax = sphere[l][3];
    }
    
    bool last   = false;
    bool abort  = false;
    int counter = 0;
    radius = 1;
    int currentI, currentJ;
    
    // main integrating function
    while (counter < 100) {
        counter++;
        while (abort == false) {    
            // scanning left edge
            currentI = iFirst-radius;      
            currentJ = jFirst-radius;
            cout << "currentI/J = " << currentI << "," << currentJ << endl;
            for (int l = 0; l < nS; l++) {
                if (sphere[l][2] == currentI && sphere[l][3] == currentJ) {
                    sphere[l][20] = 1;
                    DoAverageSurroundings(sphere);
                    currentJ++;
                    l = -1;
                    if (currentJ > jFirst+radius or currentJ > jMax) break;    
                }
            }
            // scanning right edge
            currentI = iFirst+radius;
            currentJ = jFirst-radius;
            for (int l = 0; l < nS; l++) {
                if (sphere[l][2] == currentI && sphere[l][3] == currentJ) {
                    sphere[l][20] = 1;
                    DoAverageSurroundings(sphere);
                    currentJ++;
                    l = -1;
                    if (currentJ > jFirst+radius or currentJ > jMax) break;
                }
            } 
            // scanning upper edge
            currentI = iFirst-(radius-1);
            currentJ = jFirst+radius;
            for (int l = 0; l < nS; l++) {
                if (sphere[l][2] == currentI && sphere[l][3] == currentJ) {
                    sphere[l][20] = 1;
                    DoAverageSurroundings(sphere);
                    currentI++;
                    l = -1;
                    if (currentI > iFirst+(radius-1) or currentI > iMax) break;
                }
            }
            // scanning lower edge
            currentI = iFirst-(radius-1);
            currentJ = jFirst-radius;
            for (int l = 0; l < nS; l++) {
                if (sphere[l][2] == currentI && sphere[l][3] == currentJ) {
                    sphere[l][20] = 1;
                    DoAverageSurroundings(sphere);
                    currentI++;
                    l = -1;
                    if (currentI > iFirst+(radius-1) or currentI > iMax) break;
                }
            }
            radius++;
            if (radius > 1.2*((iMax-iMin)/2)) {     // breaking the loop if it is sure that all segments have been considered
                cout << "value 'radius' at limit" << endl;
                break;    
            }
        }
        for (int l = 0; l < nS; l++) {  //summing up all average differences
            avDiff += sphere[l][18];
        }
        cout << "-------------------------------------------------------avDiff = " << avDiff << " -----------------------------------------------------" << endl;            
        if (abs(avDiff) <= threshold) {
            cout << "---------------------------------------------------avDiff < Threshold !! -------------------------------------------------" << endl;
            break;
        }
        for (int l = 0; l < nS; l++) {
            sphere[l][18] = 0;
        }
        if (counter > 99) {
            cout << "------------------------------------------------------- COUNTER > 99 -----------------------------------------------------------" << endl;
        }
        radius = 1;
    } 
    for (int l = 0; l < nS; l++) {
        cout << "sphere[" << l << "] = " << "[" << sphere[l][0] << "," << sphere[l][1] << "," << sphere[l][2] << "," << sphere[l][3] << "," << sphere[l][4] << "," << sphere[l][5] << "," << sphere[l][6] << "," << sphere[l][7] << "," << sphere[l][8] << "," << sphere[l][9] << "," << sphere[l][10] << "," << sphere[l][11] << "," << sphere[l][12] << "," << sphere[l][13] << "," << sphere[l][14] << "," << sphere[l][15] << "," << sphere[l][16] << "," << sphere[l][17] << "," << sphere[l][18] << "," << sphere[l][19] << "," << sphere[l][20] << "]" << endl;
    }
    cout << endl << "sphere: [x,y,i,j,normalX,normalY,upEdge,lowEdge,lEdge,rEdge,hCenterSingle, heightY, totHeight,index1, upper, lower, left, right, averageDev, hCenterSummed, index2, index3]"<<endl;
    cout <<         "        [0 1 2 3   4        5      6       7      8     9        10          11         12       13     14     15    16     17       18           19          20      21  ]" << endl << endl;     
}


// calculating average height difference between the four edges of current segment with edges of neighboured segments
void CbmRichRonchiAna::DoAverageSurroundings(vector<vector<double> >& sphere)
{
    int nS = sphere.size();
    int initI, initJ, initL;
    int counter = 0;
    
    double diffTop   = 0; 
    double diffBot   = 0;
    double diffLeft  = 0;
    double diffRight = 0; 
    double avDiff    = 0;
    
    // to get dates of initial segment
    for (int l = 0; l < nS; l++) {
        if (sphere[l][20]==1) {
            initI = sphere[l][2];
            initJ = sphere[l][3];
            initL = l;
            sphere[l][20]=0;
            break;
        }
    }
    
    // calculating height differences of all four adjacent edges
    for (int l = 0; l < nS; l++) {      // doppelt belegte index-kombinationen lassen counter auf ueber 4 steigen !!
        if (sphere[l][2] == initI && sphere[l][3] == initJ+1) {
            diffTop = sphere[l][15]-sphere[initL][14];
            counter++;
        }
        if (sphere[l][2] == initI && sphere[l][3] == initJ-1) {
            diffBot = sphere[l][14]-sphere[initL][15];
            counter++;
        }
        if (sphere[l][2] == initI-1 && sphere[l][3] == initJ) {
            diffLeft = sphere[l][17]-sphere[initL][16];
            counter++;
        }
        if (sphere[l][2] == initI+1 && sphere[l][3] == initJ) {
            diffRight = sphere[l][16]-sphere[initL][17];
            counter++;
        }
    }
    
    double oldHeight;
    double newHeight;
    if (counter != 0) {
        avDiff = (diffTop+diffBot+diffLeft+diffRight)/counter;
        sphere[initL][18] = avDiff;
        oldHeight         = sphere[initL][19];
        sphere[initL][14]+= avDiff;
        sphere[initL][15]+= avDiff;
        sphere[initL][16]+= avDiff;
        sphere[initL][17]+= avDiff;
        sphere[initL][19]+= avDiff;
        newHeight  = sphere[initL][19];
        //cout << "At sphere[" << initL << "] " << "difference up/down/left/right/av | counter = " << diffTop << ", " << diffBot << ", " << diffLeft << ", " << diffRight << ", " << avDiff << " | " << counter << endl;
    }
    else {cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!NO SEGMENTS TO AVERAGE OVER!" << endl;}   
    //cout << "oldHeight = " << oldHeight << ", newHeight = " << newHeight << endl;
}


void CbmRichRonchiAna::DoSphere(vector<vector<double> >& sphere)
{
    int nS = sphere.size();    
    int xMin, xMax, yMin, yMax; // to find intersection that is closest to the center
    int iMin = 1000;
    int iMax = -1000;
    int jMin = 1000;
    int jMax = -1000;
    int iFirst, jFirst, lFirst; // dates of intersection that is closest to the center
    int currentJ;
    int radius1 = 1;
    double radius2;    // distance of first found intersection to the center (if function with radius1 would describe a circle (its a square), this would be redundant and radius1 coudl be used instead
        
    double totHeightY;
    
    for (int l = 0; l < nS; l++) {  // to know range of i and j to terminate loops
        if (sphere[l][2] < iMin) iMin = sphere[l][2];
        if (sphere[l][2] > iMax) iMax = sphere[l][2];
        if (sphere[l][3] < jMin) jMin = sphere[l][3];
        if (sphere[l][3] > jMax) jMax = sphere[l][3];
    }
    
    cout << "i/j Min/Max = " << iMin << "," << iMax << "," << jMin << "," << jMax << endl;
    
    // to find the intersection that is closest to the center
    for (int l = 0; l < nS; l++) {      
        xMin = fCenterX - radius1;
        xMax = fCenterX + radius1;
        yMin = fCenterY - radius1;
        yMax = fCenterY + radius1;
        if (sphere[l][0] >= xMin && sphere[l][0] <= xMax && sphere[l][1] >= yMin && sphere[l][1] <= yMax) {
            iFirst = sphere[l][2];
            jFirst = sphere[l][3];
            lFirst = l;
            sphere[l][21] = 1;
            cout << endl << "Found center: sphere[" << l << "] = [" << sphere[l][0] << "," << sphere[l][1] << "," << sphere[l][2] << "," << sphere[l][3] << ", ...], radius1 = " << radius1 << endl;
            break;
        }
        if (l == nS-1) {
            radius1++;
            l = -1;
        }
    }
    
    // summing up the center heights of segments
    // scanning upper part
    
    currentJ   = jFirst;
    radius2    = sqrt(pow((sphere[lFirst][0]-fCenterX)*fEdgeLengthPixel,2) + pow((sphere[lFirst][1]-fCenterY)*fEdgeLengthPixel,2));
    totHeightY = fRadiusMirror-sqrt(pow(fRadiusMirror,2) - pow(radius2,2)) -sphere[lFirst][10];     // rising height (Pfeilhoehe)
    
    cout << "radius2 = " << radius2 << endl;    
    cout << "totHeightY (h) - heightFirst = " << totHeightY << endl; 
    
    for (int l = 0; l < nS; l++) {
        if (sphere[l][2] == iFirst && sphere[l][3] == currentJ) {
            totHeightY   += sphere[l][10];
            sphere[l][11] = totHeightY;
            sphere[l][12] = totHeightY;
            sphere[l][14] = totHeightY + sphere[l][6];
            sphere[l][15] = totHeightY + sphere[l][7];
            sphere[l][16] = totHeightY + sphere[l][8];
            sphere[l][17] = totHeightY + sphere[l][9];
            sphere[l][19] = totHeightY;
            sphere[l][13] = 1;
            sphere[l][21] = 1;
            DoScanLineHeight(sphere);
            currentJ++;
            if (currentJ > jMax) break;
            l = -1;
        }
    }
       
    // scanning lower part 
    totHeightY = -sphere[lFirst][10];
    currentJ   = jFirst-1; 
    for (int l = 0; l < nS; l++) {
        if (sphere[l][2] == iFirst && sphere[l][3] == currentJ) {
            totHeightY   += sphere[l][10];
            sphere[l][11] = totHeightY;
            sphere[l][12] = totHeightY;
            sphere[l][14] = totHeightY + sphere[l][6];
            sphere[l][15] = totHeightY + sphere[l][7];
            sphere[l][16] = totHeightY + sphere[l][8];
            sphere[l][17] = totHeightY + sphere[l][9];
            sphere[l][19] = totHeightY;
            sphere[l][13] = 1;
            sphere[l][21] = 1;
            DoScanLineHeight(sphere);
            currentJ--;
            if (currentJ < jMin) break;
            l = -1;
        }
    }   
    
    // calculating remaining segments
    int counter = 0;    // to terminate loop
    for (int l = 0; l < nS; l++) {
        if (sphere[l][21] == 0) {
            sphere[l][21] = 2;
            DoCalculateRemaining(sphere);
            counter++;
            cout << "counter_remaining = " << counter << " at l = " << l << ": [" << sphere[l][0] << "," << sphere[l][1] << "]" << endl;
            l = -1;          
        }
        if (counter>1000) break;
    }
}


// searching for remainig, not calculated segments
void CbmRichRonchiAna::DoCalculateRemaining(vector<vector<double> >& sphere)
{
    int nS = sphere.size();
    int currentI, currentJ, initL;
    
    for (int l = 0; l < nS; l++) {
        if (sphere[l][21] == 2) {
            currentI = sphere[l][2];
            currentJ = sphere[l][3];
            initL = l;
            sphere[l][13] = 1;
            sphere[l][21] = 1;
            break;
        }
    }
    // !!! WARUM BEKOMMEN DIE SEGMENTE IN DEN ECKEN FALSCHE HEIGHTY-WERTE ZUGETEILT??? !!!
    int heightY;
    for (int l = 0; l < nS; l++) {
        if (sphere[l][2] == currentI && sphere[l][3] == currentJ-1) {
            if (sphere[initL][1] <= fCenterY) {heightY = sphere[l][11]+sphere[initL][10];}
            if (sphere[initL][1] >  fCenterY) {heightY = sphere[l][11]+sphere[initL][10];}
            sphere[initL][11] = heightY;
            sphere[initL][12] = heightY;
            sphere[initL][14] = heightY + sphere[l][6];
            sphere[initL][15] = heightY + sphere[l][7];
            sphere[initL][16] = heightY + sphere[l][8];
            sphere[initL][17] = heightY + sphere[l][9];
            sphere[initL][19] = heightY;
            DoScanLineHeight(sphere);
            break;
        }
    } 
    
    for (int l = 0; l < nS; l++) {
        if (sphere[l][2] == currentI && sphere[l][3] == currentJ+1) {
            if (sphere[initL][1] <= fCenterY) heightY = sphere[l][11]+sphere[initL][10];
            if (sphere[initL][1] >  fCenterY) heightY = sphere[l][11]+sphere[initL][10];
            sphere[initL][11] = heightY;
            sphere[initL][12] = heightY;
            sphere[initL][14] = heightY + sphere[l][6];
            sphere[initL][15] = heightY + sphere[l][7];
            sphere[initL][16] = heightY + sphere[l][8];
            sphere[initL][17] = heightY + sphere[l][9];
            sphere[initL][19] = heightY;
            DoScanLineHeight(sphere);
            break;
        }
    }  
}


void CbmRichRonchiAna::DoScanLineHeight(vector<vector<double> >& sphere)
{   
    int    nS = sphere.size();
    int    iInit, currentI, currentJ;
    int    iMin = 1000;
    int    iMax = -1000;
    double heightY, totHeight;    // heights in x and y, measured at center of segments
    
    for (int l = 0; l < nS; l++) {
        if (sphere[l][2] < iMin) iMin = sphere[l][2];
        if (sphere[l][2] > iMax) iMax = sphere[l][2];
    }
    
    for (int l = 0; l < nS; l++) {
        if (sphere[l][13] == 1) {
            heightY       = sphere[l][11];
            totHeight     = sphere[l][11];
            iInit         = sphere[l][2];
            currentI      = sphere[l][2];
            currentJ      = sphere[l][3];
            sphere[l][13] = 0;
            break;
        }
    }
    
    // scanning line to the left
    for (int l = 0; l < nS; l++) {
        if (sphere[l][2] == currentI-1 && sphere[l][3] == currentJ) {
            if (sphere[l][21] == 1) break;
            totHeight    += sphere[l][10];
            sphere[l][11] = heightY;
            sphere[l][12] = totHeight;
            sphere[l][14] = totHeight + sphere[l][6];
            sphere[l][15] = totHeight + sphere[l][7];
            sphere[l][16] = totHeight + sphere[l][8];
            sphere[l][17] = totHeight + sphere[l][9];
            sphere[l][19] = totHeight;
            sphere[l][21] = 1;
            currentI--;
            if (currentI < iMin) break;
            l = -1;
        }
    }
            
    // scanning line to the right
    totHeight = heightY;
    currentI  = iInit;
    for (int l = 0; l < nS; l++) {
        if (sphere[l][2] == currentI+1 && sphere[l][3] == currentJ) {
            if (sphere[l][21] == 1) break;
            totHeight += sphere[l][10];
            sphere[l][11] = heightY;
            sphere[l][12] = totHeight;
            sphere[l][14] = totHeight + sphere[l][6];
            sphere[l][15] = totHeight + sphere[l][7];
            sphere[l][16] = totHeight + sphere[l][8];
            sphere[l][17] = totHeight + sphere[l][9];
            sphere[l][19] = totHeight;
            sphere[l][21] = 1;
            currentI++;
            if (currentI > iMax) break;
            l = -1;
        }
    }   
}


vector<vector<double>> CbmRichRonchiAna::DoHeight(vector<vector<double> >& normal)  // here heights of single segments are being calculated
{
    vector<vector<double>> sphere;
    vector<double> here;
    
    int    nC = normal.size();
    
    double angleRadX, angleRadY;
    double upperEdge, lowerEdge, leftEdge, rightEdge;
    double heightCenter;
    
    for (int l = 0; l < nC; l++) {
        //angleRadX = (normal[l][4]*fPi)/180.0;   // transform angle from DEG to RAD
        //angleRadY = (normal[l][5]*fPi)/180.0;
        angleRadX = (normal[l][4]);   
        angleRadY = (normal[l][5]); 
        rightEdge = (fLineDistance*fEdgeLengthPixel*(fEffectiveEdgeLengthMirror/fEffectiveEdgeLengthCCD)/2.0)*sin(atan(angleRadX));
        leftEdge  = -rightEdge; 
        upperEdge = (fLineDistance*fEdgeLengthPixel*(fEffectiveEdgeLengthMirror/fEffectiveEdgeLengthCCD)/2.0)*sin(atan(angleRadY)); 
        lowerEdge = -upperEdge;
        heightCenter = abs(upperEdge) + abs(rightEdge); //fLineDistance*fEdgeLengthPixel*(sin(angleRadX) + sin(angleRadY))/2.0; // !!! stimmt wahrscheinlich noch nicht ganz; sin*sin anstatt sin+sin ??
        here.push_back(normal[l][0]);
        here.push_back(normal[l][1]);
        here.push_back(normal[l][2]);
        here.push_back(normal[l][3]);
        here.push_back(normal[l][4]);  
        here.push_back(normal[l][5]);
        here.push_back(upperEdge);
        here.push_back(lowerEdge);
        here.push_back(leftEdge);
        here.push_back(rightEdge);
        here.push_back(heightCenter);
        here.push_back(0);     // (11) heightY; this and next will be filled later
        here.push_back(0);     // (12) total height
        here.push_back(0);     // (13) index for intersection to calculate
        here.push_back(0);     // (14) indices for: upper
        here.push_back(0);     // (15) lower
        here.push_back(0);     // (16) left and
        here.push_back(0);     // (17) right Edge after integration
        here.push_back(0);     // (18) average height difference with neighboured segments
        here.push_back(0);     // (19) height of center after integration                   !! scheint sich nicht von (12) zu unterscheiden !?
        here.push_back(0);     // (20) index2: for integration function, indicates integrated segments
        here.push_back(0);     // (21) index3: to indicate if corresponding segment had been assigned to height value (in 'heightY' and 'summed height')
        sphere.push_back(here);
        here.clear();
    }
    return sphere;
}

vector<vector<double>> CbmRichRonchiAna::DoLocalNormal(vector<vector<int> >& matrix)
{
    vector<vector<double>> normal;
    vector<double> norm;
    
    int nM = matrix.size(); 
    
    double xR, yR, xM, yM;             // positions in x and y on ronchi ruling and on mirror in microns
    double sIncX, sIncY, sRefX, sRefY; // slopes in x and y direction of incident and reflected beam
    double normalX, normalY;           // local normal
    
    int iMin =  1000;
    int iMax = -1000;
    int jMin =  1000;
    int jMax = -1000;
    int iCenter, jCenter;
    
    for (int l = 0; l < nM; l++) {  // to know range of i and j to terminate loops
        if (matrix[l][2] < iMin) iMin = matrix[l][2];
        if (matrix[l][2] > iMax) iMax = matrix[l][2];
        if (matrix[l][3] < jMin) jMin = matrix[l][3];
        if (matrix[l][3] > jMax) jMax = matrix[l][3];
    }
    
    iCenter = (iMax-iMin)/2;
    jCenter = (jMax-jMin)/2;
    
    cout << "i/jCenter = " << iCenter << ", " << jCenter << endl;
    
    //calculating slope and filling vector
    for (int l = 0; l < nM; l++) {
        //xR      = (matrix[l][0]-fImageWidth/2.0)*fEdgeLengthPixel;
        //yR      = (matrix[l][1]-fImageWidth/2.0)*fEdgeLengthPixel;   
        xR      = (matrix[l][2]-iCenter)*fPitchGrid;   
        yR      = (matrix[l][3]-jCenter)*fPitchGrid;   
        sRefX   = (double)(xR/fDistRulingCCD)*(1-fEffectiveEdgeLengthCCD/fEffectiveEdgeLengthGridPlate)/2;    // dividing by TWO would give much more exact results for slopes !!!
        sRefY   = (double)(yR/fDistRulingCCD)*(1-fEffectiveEdgeLengthCCD/fEffectiveEdgeLengthGridPlate)/2;
        //sRefX   = (double)((matrix[l][2]-iCenter)*fPitchGrid-(matrix[l][0]-fCenterX)*fEdgeLengthPixel)/fDistRulingCCD;    // these are not working yet; grid and CCD are obviously behing CoC, so...
        //sRefY   = (double)((matrix[l][3]-jCenter)*fPitchGrid-(matrix[l][1]-fCenterY)*fEdgeLengthPixel)/fDistRulingCCD;    // ... that image on CCD is bigger than image on ruling !!
        xM      = (sRefX*fDistMirrorRuling);
        yM      = (sRefY*fDistMirrorRuling);
        sIncX   = (xM/fRadiusMirror);     // if point source is in the center of curvature !!!
        sIncY   = (yM/fRadiusMirror);
        normalX = (sRefX+sIncX)/2.0;
        normalY = (sRefY+sIncY)/2.0;   
        //cout << "yR = " << yR << ", sRefY = " << sRefY << ", yM = " << yM << ", NormalY = " << normalY << ", sIncY = " << sIncY << endl;
        norm.push_back(matrix[l][0]);
        norm.push_back(matrix[l][1]);
        norm.push_back(matrix[l][2]);
        norm.push_back(matrix[l][3]);
        norm.push_back(normalX);
        norm.push_back(normalY);
        norm.push_back(xR); // !! these next 4 values can be deleted afterwards, only for controlling of values
        norm.push_back(yR);
        norm.push_back(xM);
        norm.push_back(yM);
        normal.push_back(norm);
        norm.clear();
        //cout << "normal[" << l <<  "][4.5] = [" << normal[l][4] << "," << normal[l][5] << "]" << endl;
        
        cout << "normal[" << l << "] = " << "[" << normal[l][0] << "," << normal[l][1] << "," << normal[l][2] << "," << normal[l][3] << "," << normal[l][4] << "," << normal[l][5] << "," << normal[l][6] << "," << normal[l][7] << "," << normal[l][8] << "," << normal[l][9] << "]"  << endl;
   }
      
   return normal;
}


/*
EXPLANATION TO VECTORS:

1)
Entries of vector<vector<int>> intersectionXY: [x,y,indicator_assigned,indicator_for_DoSearchBelow/Above_function];
indicator_assigned: is ZERO in the beginning, turns to ONE if according intersection has been assigned to 'i' and 'j' and loaded into 'matrix'.
indicator_for_DoSearchBelow/Above_function: is usually ZERO, only if DoNumberInt finds this point to be assigned and to hand over to DoSearchBelow/Above function, it will be set to ONE, but directly deleted inbetween the latter named function

2)
Entries of vector<vector<int>> matrix: [x,y,i,j,xInit,yInit,iInit,jInit]; 
The ...Inits are for orientation when DoScanLine has finished with searching on the left side for further entries and continues at xInit with iInit to search on the right side of initial point of according line (jInit actually is not necessary, because inbetween one line it does not change)
*/


// numbering intersections 
vector<vector<int>> CbmRichRonchiAna::DoNumberIntersections(vector<vector<int> >& intersectionXY)    // includes Searching and Scanning Lines
{
    const int nI    = intersectionXY.size(); 
    
    vector<vector<int> > matrix;
    vector<int> mix;
    
    const int originX = intersectionXY[0][0];
    const int originY = intersectionXY[0][1];
    
    int currentX = originX;     // this and next line to enable including first line into search function as well (and not run it seperately in the beginning)
    int currentY = originY - fLineDistance/2;
    int xInit = originX;
    int yInit = originY;
    int iInit = 0;
    int jInit = -1;
              
    int i = 0;  // column
    int j = 0;  // line
    int l = 0;  // l_th argument of 'matrix'
      
    for (int l1 = 0; l1 < nI; l1++) {
        if ( (intersectionXY[l1][0] - currentX > -fLineDistance/2) && (intersectionXY[l1][0] - currentX < fLineDistance/2) && (intersectionXY[l1][1] - currentY < 1.5*fLineDistance) && (intersectionXY[l1][1] - currentY > 0) ) {
            intersectionXY[l1][2] = 1;
            currentX = intersectionXY[l1][0];
            currentY = intersectionXY[l1][1];
            xInit = intersectionXY[l1][0];
            yInit = intersectionXY[l1][1];
            i = iInit;
            jInit++;
            mix.push_back(currentX);
            mix.push_back(currentY);
            mix.push_back(i);
            mix.push_back(j);
            mix.push_back(xInit);
            mix.push_back(yInit);
            mix.push_back(iInit);
            mix.push_back(jInit);
            matrix.push_back(mix);
            mix.clear();
            //cout << "matrix[" << l << "] (main function) = [" << matrix[l][0] << "," << matrix[l][1] << "," << matrix[l][2] << "," << matrix[l][3] << "][" << matrix[l][4] << "," << matrix[l][5] << "," << matrix[l][6] << "," << matrix[l][7] << "]" << endl;
            DoScanLine(intersectionXY, matrix);
            l = matrix.size();
            currentX = matrix[l-1][4];
            currentY = matrix[l-1][5]+fLineDistance/2;
            break;
        }
    }
    
    //numbering remaining intersections
    for (int l1 = 0; l1 < nI; l1++) {
        if (intersectionXY[l1][2] == 0) {
            intersectionXY[l1][3] = 1;
            DoSearchBelow(intersectionXY, matrix);
        }
    }
    
    for (int l1 = nI-1; l1 >=0; l1--) {
        if (intersectionXY[l1][2] == 0) {
            intersectionXY[l1][3] = 1;
            DoSearchAbove(intersectionXY, matrix);
        }
    }
    
    l = matrix.size();    
    cout << "After assigning remaining Ints: l = " << l << endl;
     
    return matrix;
}


// scanning line 'j' and number intersections
void CbmRichRonchiAna::DoScanLine(vector<vector<int> >& intersectionXY, vector<vector<int> >& matrix)
{
    vector<int> mix;
    
    int nI       = intersectionXY.size(); 
    int l        = matrix.size();
    int currentX = matrix[l-1][4];
    int currentY = matrix[l-1][5];    
    int i        = matrix[l-1][6];  // column
    int j        = matrix[l-1][7];  // line
    int xInit    = matrix[l-1][4];  // position of initial point
    int yInit    = matrix[l-1][5];
    int iInit    = matrix[l-1][6];  
    int jInit    = matrix[l-1][7];  // is actually redundant
       
    // scanning intersections on the left side
    int counter = 0;
    while (counter < nI && l < nI) {
        for (int l1 = 0; l1 < nI; l1++) {
            counter++;
            if (intersectionXY[l1][2] == 1) {
            //    cout << "already counted left" << endl;
                counter = nI;
                continue;
            }
            if ( (intersectionXY[l1][0] - currentX > -1.3*fLineDistance) && (intersectionXY[l1][0] - currentX < 0) && (intersectionXY[l1][1] - currentY < fLineDistance/2) && (intersectionXY[l1][1] - currentY > -fLineDistance/2) ) {
                intersectionXY[l1][2] = 1;
                i--;
                currentX = intersectionXY[l1][0];
                currentY = intersectionXY[l1][1];
                mix.push_back(currentX);
                mix.push_back(currentY);
                mix.push_back(i);
                mix.push_back(j);
                mix.push_back(xInit);
                mix.push_back(yInit);
                mix.push_back(iInit);
                mix.push_back(jInit);
                matrix.push_back(mix);
                mix.clear();
                l = matrix.size()-1;
                //cout << "matrix[" << l << "] (ScanLine_left) = [" << matrix[l][0] << "," << matrix[l][1] << "," << matrix[l][2] << "," << matrix[l][3] << "][" << matrix[l][4] << "," << matrix[l][5] << "," << matrix[l][6] << "," << matrix[l][7] << "]" << endl;
                l1 = -1; 
                counter = 0; 
            } 
        }
    }
    
    // scanning intersections on the right side
    i = iInit;
    currentX = xInit;
    currentY = yInit;
    
    counter = 0;
    while (counter < nI && l < nI) {
        for (int l1 = 0; l1 < nI; l1++) {
            counter++;
            if (intersectionXY[l1][2] == 1) {
                counter = nI;
                continue;
            }
            if ( (intersectionXY[l1][0] - currentX < 1.3*fLineDistance) && (intersectionXY[l1][0] - currentX > 0) && (intersectionXY[l1][1] - currentY < fLineDistance/1.5) && (intersectionXY[l1][1] - currentY > -fLineDistance/1.5) ) {
                intersectionXY[l1][2] = 1;
                i++;
                currentX = intersectionXY[l1][0];
                currentY = intersectionXY[l1][1];
                mix.push_back(currentX);
                mix.push_back(currentY);
                mix.push_back(i);
                mix.push_back(j);
                mix.push_back(xInit);
                mix.push_back(yInit);
                mix.push_back(iInit);
                mix.push_back(jInit);
                matrix.push_back(mix);
                mix.clear();
                l = matrix.size()-1;
                //cout << "matrix[" << l << "] (ScanLine_right) = [" << matrix[l][0] << "," << matrix[l][1] << "," << matrix[l][2] << "," << matrix[l][3] << "][" << matrix[l][4] << "," << matrix[l][5] << "," << matrix[l][6] << "," << matrix[l][7] << "]" << endl;
                l1 = -1;
                counter = 0;
            }       
        }
    }
}


// searching for intersections below first found intersection in 'intersectionXY' that had not been assigned to 'i' and 'j'
void CbmRichRonchiAna::DoSearchBelow(vector<vector<int> >& intersectionXY, vector<vector<int> >& matrix) 
{
    vector<int> mix;
    
    int nI = intersectionXY.size();
    int l  = matrix.size();
    
    int currentX, currentY;
        
    for (int l1 = 0; l1 < nI; l1++) {   // searching for ONE in the 4th column of intersectionXY (intXY[l1][3]==1)
        if (intersectionXY[l1][3] == 1) {
            currentX = intersectionXY[l1][0];
            currentY = intersectionXY[l1][1];
            intersectionXY[l1][3] = 0;
            break;
        }
    }
    
    for (int l1 = 0; l1 < l; l1++) {
        if ( (matrix[l1][0] - currentX > -fLineDistance/2) && (matrix[l1][0] - currentX < fLineDistance/2) && (matrix[l1][1] - currentY > -1.5*fLineDistance) && (matrix[l1][1] - currentY < 0) ) {
            int currentI = matrix[l1][2];
            int currentJ = matrix[l1][3]+1;
            mix.push_back(currentX);
            mix.push_back(currentY);
            mix.push_back(currentI);
            mix.push_back(currentJ);
            mix.push_back(currentX);
            mix.push_back(currentY);
            mix.push_back(currentI);
            mix.push_back(currentJ);
            matrix.push_back(mix);
            l  = matrix.size();
            mix.clear();
            //cout << "matrix[" << l-1 << "] = [" << matrix[l-1][0] << "," << matrix[l-1][1] << "," << matrix[l-1][2] << "," << matrix[l-1][3] << "][" << matrix[l-1][4] << "," << matrix[l-1][5] << "," << matrix[l-1][6] << "," << matrix[l-1][7] << "]" << endl;
            for (int l2 = 0; l2 < nI; l2++) {       // indicates that this intersection had been assigned to line and column
                if (intersectionXY[l2][0] == currentX && intersectionXY[l2][1] == currentY) {
                    intersectionXY[l2][2] = 1;  
                    break; 
                }
            }
            DoScanLine(intersectionXY, matrix);
            break;
        }
    }   
}


// searching for intersections above first found intersection in 'intersectionXY' that had not been assigned to 'i' and 'j'
void CbmRichRonchiAna::DoSearchAbove(vector<vector<int> >& intersectionXY, vector<vector<int> >& matrix) 
{
    //cout << "   SEARCHING ABOVE" << endl;
    vector<int> mix;
    
    int nI = intersectionXY.size();
    int l  = matrix.size();
    int currentX;
    int currentY;
    
    for (int l1 = 0; l1 < nI; l1++) {   // searching for ONE in the 4th column of intersectionXY (intXY[l1][3]==1)
        if (intersectionXY[l1][3] == 1) {
            currentX = intersectionXY[l1][0];
            currentY = intersectionXY[l1][1];
            intersectionXY[l1][3] = 0;
            break;
        }
    }
    
    for (int l1 = 0; l1 < l; l1++) {
        if ( (matrix[l1][0] - currentX > -fLineDistance/2) && (matrix[l1][0] - currentX < fLineDistance/2) && (matrix[l1][1] - currentY < 1.5*fLineDistance) && (matrix[l1][1] - currentY > 0) ) {
            int currentI = matrix[l1][2];
            int currentJ = matrix[l1][3]-1;
            mix.push_back(currentX);
            mix.push_back(currentY);
            mix.push_back(currentI);
            mix.push_back(currentJ);
            mix.push_back(currentX);
            mix.push_back(currentY);
            mix.push_back(currentI);
            mix.push_back(currentJ);
            matrix.push_back(mix);
            l  = matrix.size();
            mix.clear();
            //cout << "matrix[" << l-1 << "] = [" << matrix[l-1][0] << "," << matrix[l-1][1] << "," << matrix[l-1][2] << "," << matrix[l-1][3] << "][" << matrix[l-1][4] << "," << matrix[l-1][5] << "," << matrix[l-1][6] << "," << matrix[l-1][7] << "]" << endl;
            for (int l2 = 0; l2 < nI; l2++) {
                if (intersectionXY[l2][0] == currentX && intersectionXY[l2][1] == currentY) {
                    intersectionXY[l2][2] = 1;  // indicates that this intersection had been assigned to line and column
                    break; 
                }
            }
            DoScanLine(intersectionXY, matrix);
            break;
        }
    }   
}


// fading out center 
void CbmRichRonchiAna::DoShadowCenter(vector<vector<int> >& data)
{
    int nX = data.size();
    int radius = 50;
   
    for (int x = 0; x < nX; x++) {
        for (int y = 0; y < nX; y++) {
            if (x > fCenterX-sqrt(pow(radius,2)-pow(y-fCenterY,2)) && x < fCenterX+sqrt(pow(radius,2)-pow(y-fCenterY,2)) && y > fCenterY-sqrt(pow(radius,2)-pow(x-fCenterX,2)) && y < fCenterY+sqrt(pow(radius,2)-pow(x-fCenterX,2)) ) {
                data[x][y] = 0;
            }
        }
    }
}


// erasing single pixels that dont belong to a line
void CbmRichRonchiAna::DoEjectSingle(vector<vector<int> >& data)
{   
    int height = data.size();   
    int windowX = 5;
    int windowY = 3;
    int threshold = 6;
    
    for (int x = windowX; x < height-windowX; x++) {
        for (int y = windowY; y < height-windowY; y++) {
            int counter = 0;
            if (data[x][y] > 0) {
                for (int x1 = x-windowX; x1 <= x+windowX; x1++) {
                    if (x1 < 0 or x1 >= height) continue;
                    for (int y1 = y-windowY; y1 <= y+windowY; y1++) {
                        //if (y1 < 0 or y1 >= height) continue;
                        if (data[x1][y1] > 0) counter++;
                    }
                }
                if (counter < threshold) {
                    for (int x1 = x-windowX; x1 <= x+windowX; x1++) {
                        for (int y1 = y-windowY; y1 <= y+windowY; y1++) {
                            data[x1][y1] = 0;
                        }
                    }
                }
            }
        }
    }
}


// filling gaps that are not too big with linear function
void CbmRichRonchiAna::DoFillGaps(vector<vector<int> >& data)
{    
    int height = data.size();
    int winHeight = 5;
    int winWidth  = 60;
    int xLeft, xRight, yLeft, yRight;
    int gap = 1;
    int noGap = 0;
    double slope;
    bool big = false;
    
    for (int x = fLineDistance; x < height - fLineDistance; x++) {
        for (int y = fLineDistance; y < height - fLineDistance; y++) {
            if (data[x][y] > 0) {   // if a pixel is detected, look for a pixel of this line in the next column
                xLeft = x;      // memorize the position of the last detected line pixel
                yLeft = y;
                for (int y1 = y-winHeight; y1 <= y+winHeight; y1++) {   // scanning the next column to see if there is a gap or not
                    if (data[x+1][y1] > 0) noGap++;
                }
                if (noGap == 0) {   // i.e. gap detected
                    for (int x1 = x+2; x1 <= height-winWidth; x1++) {
                        for (int y2 = y-winHeight; y2 <= y+winHeight; y2++) {
                            if (data[x1][y2] > 0) {
                                noGap++;
                                xRight = x1;
                                yRight = y2;
                            }
                        }
                        if (noGap == 0) gap++;  
                        else if (noGap == 1) break;
                        if (gap > winWidth) {
                            //cout << " --------------------- GAP AT POSITION [" << x << "," << y << "] IS TOO BIG -------------------------------"   << endl;
                            gap = 0;
                            big = true;
                            break;
                        }
                    }
                    
                    if (big == false && x < 900) {  // only fill gaps if they are not too big; the condition 'x < 900'prevents continuing the lines at the end
                        slope = (double)(yRight-yLeft)/(xRight-xLeft);
                        int mult = 1;
                    
                        for (int x2 = x+1; x2 <= x+gap+1; x2++) {   // filling gaps
                            data[x2][y+(int)(mult*slope)] = 10000;
                            mult++;
                        }
                    }
                }
                else if (noGap == 1) {  // i.e. no gap detected
                    noGap = 0;
                    xLeft = 0;
                    yLeft = 0;
                    continue;
                }
                
                xLeft  = 0;
                xRight = 0;
                yLeft  = 0;
                yRight = 0;
                gap    = 0;
                noGap  = 0;
                slope  = 0;
                big = false;
            }
        }
    }
}


// calculating mean position in Y
void CbmRichRonchiAna::DoMeanY(vector<vector<int> >& data)
{
    int meanHalfLength = 10;
    int meanHalfHeight = 3;
    int nX = data.size();
    int nY = data[0].size();
    
    for (int x = meanHalfLength; x < nX-meanHalfLength; x++) {
        for (int y = meanHalfHeight; y < nY-meanHalfHeight; y++) {
            if (data[x][y] == 0) continue;  // to consider only filled pixels
            int curData = data[x][y];
            int sumY = 0;
            int divider = 0;
            for (int x2 = -meanHalfLength; x2 <= meanHalfLength; x2++) {    // summing up positions in y of current and neighboured pixels
                for (int y2 = -meanHalfHeight; y2 <= meanHalfHeight; y2++) {
                    if (data[x+x2][y+y2] > 0) {
                        sumY += y+y2;
                        divider++;
                    }
                }
            }
            data[x][y] = 0;
            y=(int) sumY/divider;   // average y-value of pixel
            data[x][y] = curData;
            curData = 0;
            sumY = 0;
        }
    }
}


// rotating the image (from vertical to horizontal)
void CbmRichRonchiAna::DoRotation(vector<vector<int> >& data)
{
    int nX = data.size();
    int nY = data[0].size();
    for (int x = 0; x < nX; x++) {
        for (int y = x+1; y < nY; y++) {
            swap(data[x][y], data[y][x]);
        }
    }
}


// getting a one dimensional line
void CbmRichRonchiAna::DoPeakFinderY(vector<vector<int> >& data)
{
    int nX = data.size();
    int nY = data[0].size();
    int halfWindow = 10;

    for (int x = 0; x < nX; x++) {
        for (int y = halfWindow; y < nY - halfWindow; y++){
            bool isPeak = (data[x][y] >= data[x][y - 1]) && (data[x][y] >= data[x][y + 1]); // comparing current pixel with neighbours in y direction
            if (!isPeak) {          // deleting this pixel, if it's not a peak
                data[x][y] = 0;
                continue;
            }
            bool isBiggest = true;
            for (int iW = -halfWindow; iW <= halfWindow; iW++) {    // comparing current pixel with 'halfWindow' next neighbours
                if (iW == 0) continue;
                if (data[x][y + iW] >= data[x][y]) {    
                    isBiggest = false;
                    break;
                }
            }
            data[x][y] = (isBiggest && isPeak)? data[x][y] : 0;     // if in both cases 'true', this pixel will be kept (the others are deleted) 
        }
    }
}


// averaging intensity with neighboured values
void CbmRichRonchiAna::DoMeanIntensityY(vector<vector<int> >& data)
{
    int nX = data.size();
    int nY = data[0].size();

    int halfAvWindow = 6;   // number of neighbours to each side to average over
    int threshold = 25;     // threshold for average intensity to delete areas that are too dark (and so don't belong to mirror image)
    vector<int> weights = {160, 80, 40, 20, 10, 5};

    for (int x = 0; x < nX; x++) {
        for (int y = 0; y < nY; y++){
            int total = 0;
            int weightSum = 0;
            for (int iW = -halfAvWindow; iW <= halfAvWindow; iW++) {    // muliplying values of neighbours with their weights
                int iWAbs = std::abs(iW);
                int weight = (iWAbs < weights.size() )? weights[iWAbs] : weights[weights.size() - 1];   // specifying 'weights' in dependence of distance to current pixel
                weightSum += weight;
                int ind = y + iW;
                if (ind < 0) ind = 0;
                if (ind >= nY) ind = nY - 1;
                total += data[x][ind] * weight;
            }
            data[x][y] = total / weightSum;                 // filling averaged intensity 
            if (data[x][y] <= threshold) data[x][y] = 0;    // deleting pixels that are too dark to yield only data for mirror and not background
        }
    }
}

void CbmRichRonchiAna::FillH2WithVector(TH2* hist, const vector<vector<int> >& data)
{
    int nX = data.size();
    int nY = data[0].size();
    
    for (int x = 0; x < nX; x++) {
        for (int y = 0; y < nY; y++){
            if (data[x][y] != 0) {
                hist->SetBinContent(x, y, data[x][y]);
            }
        }
    }
}

/*void CbmRichRonchiAna::H2FillPlane(TH2* hist, const vector<vector<int> >& data)
{
    int nX = data.size();
    int nY = data[0].size();
    
    for (int x = 0; x < nX; x++) {
        for (int y = 0; y < nY; y++){
            if (data[x][y] != 0) {
                hist->SetBinContent(x, y, data[x][y]);
            }
        }
    }
}*/


// finding intersections
vector<vector<int>> CbmRichRonchiAna::DoIntersection(vector<vector<int> >& dataH, const vector<vector<int> >& dataV)
{
    int nX = dataV.size();
    int nY = dataV[0].size();
    
    vector<vector<int>> intersectionXYPrelim;   // preliminary vector, stores all intersection points; will be deleted afterwards
    vector<vector<int>> intersectionXY;         // contains all intersections that are in 'intXYPrelim', except that ones that are not on regularly line (too close to each other)
    vector<int> intXY;                          // vector to fill vector 'intXYPrelim'
    
    // filling 'intXYPrelim' with x,y-data of intersections
    for (int y = fLineDistance; y < nY-fLineDistance; y++) {
        for (int x = fLineDistance; x < nX-fLineDistance; x++) {
            if (dataH[x][y] > 0 && dataV[x][y] > 0) {     // if horizontal and vertical image have data > ZERO at point [x][y]
                for (int x1 = x-1; x1<=x+1; x1++) {       // to prevent double counting of intersection points
                    for (int y1 = y-1; y1<=y+1; y1++) {
                        if (x1==x && y1==y) continue;
                        if (dataH[x1][y1] > 0 && dataV[x1][y1] > 0) dataH[x1][y1] = 0;  // delete one of the double counted intersection
                    }
                }
                intXY.push_back(x);
                intXY.push_back(y);
                intXY.push_back(0);   // this value will turn to ONE if associated intersection will be assigned to 'i' (column) and 'j' (line)
                intXY.push_back(0);   // to hand over intersection that is in question from 'DoNumberInt'-function to 'DoSearchAbove/Below'-function
                intersectionXYPrelim.push_back(intXY);
                intXY.clear();
            }
        }
    }
    
    int nI = intersectionXYPrelim.size();
    
    cout << "Number of intersections before sorting out: nI1 = " << nI << endl;
    
    // sorting out intersections that are between two regularly intersections
    double tol = 0.6*fLineDistance;     // tolerance value for how close two intersections can be to still being considered
    for (int l = 0; l < nI; l++) {   
        int x = intersectionXYPrelim[l][0];
        int y = intersectionXYPrelim[l][1];
        int counter = 0;
        for (int l1 = 0; l1 < nI; l1++) {   // looking around each intersection, if there are at least two, that are too close to current one
            if (intersectionXYPrelim[l1][0]-x < tol && intersectionXYPrelim[l1][0]-x > -tol && intersectionXYPrelim[l1][1]-y < tol && intersectionXYPrelim[l1][1]-y > -tol ) {
                counter++;
            }
        }
        if (counter > 1) continue;      // if an intersection point has at least two neighbours, that are too close, it's presumed, that its not on a regular line and so not being counted
        else intersectionXY.push_back(intersectionXYPrelim[l]);
    }
    
    nI = intersectionXY.size();
    cout << "Number of intersections after sorting out: nI2 = " << nI << endl;
    
    return intersectionXY;
}


// superposing horizontal and vertical image
vector<vector<int>> CbmRichRonchiAna::DoSuperpose(const vector<vector<int> >& dataH, const vector<vector<int> >& dataV)
{
    int nX = dataV.size();
    int nY = dataV[0].size();
    
    vector<vector<int> > dataSup;
    dataSup.resize(nX);
    for (int x = 0; x < nX; x++) {
        dataSup[x].resize(nY);
        for (int y = 0; y < nY; y++) {
            dataSup[x][y] = dataH[x][y] + dataV[x][y];  // data of horizontal and vertical image are being summed
        }
    }
    return dataSup;
}


vector<vector<int> > CbmRichRonchiAna::ReadTiffFile(const string& fileName)
{
    vector<vector<int> > data;
    rgba8_image_t img;
    tiff_read_and_convert_image(fileName,img);

    int height = img.height();
    int width = img.width();

    data.resize(width);

    auto view = const_view(img);
    for (int x = 0; x < width; ++x) 
    {
        auto it = view.col_begin(x);
        data[x].resize(height);
        for (int y = 0; y < height; ++y) 
        {
            int r = boost::gil::at_c<0>(it[y]);
            //int g = boost::gil::at_c<1>(it[y]);
            //int b = boost::gil::at_c<2>(it[y]);
            //int a = boost::gil::at_c<3>(it[y]);
            data[x][y] = r;
        }
    }
    return data;
}


ClassImp(CbmRichRonchiAna)

/*
// FINDING BASE INTERSECTION
vector<vector<int>> CbmRichRonchiAna::DoFindBasePoint(const vector<vector<int>>& dataH, const vector<vector<int>>& dataV, const vector<pair<int,int>>& intersectionXY)
{
    vector<vector<int>> intNumberXY; 
    pair<int,int> basePointXY;

   
    int nX = dataH[0].size();
    int nY = dataH.size();
    int nI = intersectionXY.size();
    int l = 0;
    
    int currentY = 0;
    int startx = 0;
    
    //bool line = false;
    //bool intersection = false;
    
    
    // if first continuous line is considered as base line   
    /*
    for (int y0 = fLineDistance; y0 < nY; y0++) {    // finding base line
        if (dataH[511][y0] > 0) {
            currentY = y0;
            break;
        }
    }
            
    for (int x = 510; x > fLineDistance; x--) {  // finding beginning of base line
        line = false;
        for (int ySearch = currentY-1; ySearch <= currentY+1; ySearch++) {
            if (dataH[x][ySearch] > 0) {
                currentY = ySearch;
                line = true;
                break;
            }
        }
        if (line == false) {
            startx = x+1;
            break;
        }
    }
    
    intNumberXY.resize(nI);
    for (int x = startx; x < nX; x++) {  // finding first intersection on base line
        intNumberXY[0].resize(4);
        for (int ySearch = currentY-1; ySearch <= currentY+1; ySearch++) {
            if (dataH[x][ySearch] > 0 && dataV[x][ySearch] > 0) {
                intersection = true;
                //intNumberXY[0].push_back(x);
                //intNumberXY[0].push_back(ySearch);
                //intNumberXY[0].push_back(0);
                //intNumberXY[0].push_back(0);
                intNumberXY[0][0] = x;
                intNumberXY[0][1] = ySearch;
                intNumberXY[0][2] = 0;
                intNumberXY[0][3] = 0;
            }
        }
        if (intersection == true) break;
    }
    
    
    
    // if first line bottom left is considered as base line
    bool foundLine = false;
    
    intNumberXY.resize(nI);
    for (int x = 1.5*fLineDistance; x < nX; x++) {       
        for (int y = fLineDistance; y < 3*fLineDistance; y++) {
            if (dataH[x][y] > 0) {
                intNumberXY[0].push_back(x);
                intNumberXY[0].push_back(y);
                intNumberXY[0].push_back(0);
                intNumberXY[0].push_back(0);
                //intNumberXY[0][0] = x;
                //intNumberXY[0][1] = y;
                foundLine = true;
                break;
            }
        }
        if (foundLine == true) break;
    }
    cout << "Base Point at [x,y;i,j] = [" << intNumberXY[0][0] << "," << intNumberXY[0][1] << ";" << intNumberXY[0][2] << "," << intNumberXY[0][3] << "]" << endl; 
    return intNumberXY;  
}*/

/*
// searching for next line
int CbmRichRonchiAna::DoSearchNextLine(const vector<vector<int>>& dataH, const vector<vector<int>>& dataV, vector<vector<int>>& intNumberXY)
{
    int nX = dataH[0].size();
    int nY = dataH.size();
    
    int i = 0;  // line index of intersection
    int j = 0;  // column index of intersection
    int l = 0;  // l-th argument of vector intNumberXY   
    int xInit = intNumberXY[0][0];     
    int yInit = intNumberXY[0][1];        
                
    int values[] = {i,j,l,xInit,yInit,0};   // last value indicates if last line had been reached, to cancel for further searching
    
   cout << "Intersection 0: ["  << intNumberXY[0][0] << "," << intNumberXY[0][1] << "," << intNumberXY[0][2] << "," << intNumberXY[0][3] << "]" << endl; 
        
    DoScanLine(dataH, dataV, intNumberXY, values);  // scanning and numbering intersections of first line
    
    int currentX = xInit;
    int currentY = yInit;
    int counterX = 0;
    int counterY = 0;
    
    for (int y = yInit+1; y <= nY-2*fLineDistance; y++) {  // search next intersection above
        counterY++;
        //if (values[5] == 1) break;
        if (y > 0.9*nY) break;          // can be deleted if bool 'lastLine' in 'DoScanLine' is active 
        for (int x = currentX-1; x <= currentX+1;x++) {
            if (dataV[x][y] > 0 && dataH[x][y] > 0) {
                values[1]++;    // j
                values[2]++;    // l
                values[3] = x;  // xInit
                values[4] = y;  // yInit
                j = values[1];
                l = values[2];
               cout << "l = " << l << "    j = " << j << endl;
                //xInit = x;    not needed
                yInit = y; 
                intNumberXY[l].push_back(x);
                intNumberXY[l].push_back(y);
                intNumberXY[l].push_back(i);
                intNumberXY[l].push_back(j);
                counterY = 0;    
               cout << "IntersectionBase  " << l << ": [" << intNumberXY[l][0] << "," << intNumberXY[l][1] << "," << intNumberXY[l][2] << "," << intNumberXY[l][3] << "]" << endl;
                DoScanLine(dataH, dataV, intNumberXY, values);
                break;
            }
            else if (dataV[x][y] > 0 && dataH[x][y] == 0) {
                currentX = x;
                break;
            }
            if (counterY > 2*fLineDistance) {  // if no intersection in range of 2*fLineDistance, go to next intersection on the right and search from there
               cout << "Skipping" << endl;
                xInit = values[3];
                yInit = values[4];
               cout << "xInit before = " << xInit << "    yInit before = " << yInit << endl;
                currentY = yInit;
                bool foundNext = false;
                for (int x2 = xInit+1; x2 < nX-2*fLineDistance; x2++) {
                    counterX++;
                    for (int y2 = currentY-1; y2 <= currentY+1; y2++) {
                        if (dataV[x2][y2] > 0 && dataH[x2][y2] > 0) {
                            currentX = x2;
                            cout << "y = " << y << endl;
                            xInit = x2;
                            yInit = y2;
                            y = yInit;  // y = yInit+1? why does it work yet?
                            i++;
                            values[0] = i;
                            foundNext = true;
                            counterX = 0;
                            counterY = 0;
                           cout << "xInit after = " << xInit << "    yInit after = " << yInit << endl;
                            break;                           
                        }
                        else if (dataH[x2][y2] > 0) {
                            currentY = y2;
                            break;
                        }
                        if (counterX > 2*fLineDistance) {
                           cout << "Gap or end of line!" << endl;
                            counterX = 0;
                            break;
                        }
                    }
                    if (foundNext == true) break;
                }
                if (foundNext == true) break;
            }            
        }        
    }
}  


// scanning the line 'j' and numbering intersections
void CbmRichRonchiAna::DoScanLine(const vector<vector<int>>& dataH, const vector<vector<int>>& dataV, vector<vector<int>>& intNumberXY, int values[])     // marks intersections on current line and puts data into vector intNumberXY
{
    int nX = intNumberXY.size();
    int initialI = values[0];
    int currentY = values[4];
        
    int i = values[0];
    int j = values[1];
    int l = values[2];
    int xInit = values[3];
    int yInit = values[4];
           
    bool line = false;
    bool lineAbove = false; // both this and next bool to detect last line and end the program (has yet to be installed)
    bool lastLine = false;
    
    for (int x = xInit-1; x > 30; x--) {    // searching for intersections on the left side of initial intersection
        line = false;
        for (int ySearch = currentY-1; ySearch <= currentY+1; ySearch++) {
            if (dataH[x][ySearch] > 0 && dataV[x][ySearch] > 0) {
                currentY = ySearch;
                i--;
                l++;
                intNumberXY[l].push_back(x);
                intNumberXY[l].push_back(ySearch);
                intNumberXY[l].push_back(i);
                intNumberXY[l].push_back(j);
                line = true;
                cout << "IntersectionLeft " << l << ": [" << intNumberXY[l][0] << "," << intNumberXY[l][1] << "," << intNumberXY[l][2] << "," << intNumberXY[l][3] << "]" << endl;
                break;
            }
            else if (dataH[x][ySearch] > 0) {
                currentY = ySearch;
                line = true;
                break;
            }
        }
        if (line == false) break;   // end of line or a gap had been reached
    }
    
    currentY = yInit;
    i = initialI;
    
    for (int x = xInit+1; x < nX-30; x++) {      // searching for intersections on the right side of initial intersection
        line = false;        
        if (x == 511) {     // looking, if this is the last continuous line
            lineAbove = false;
            for (int y2 = currentY; y2 <= currentY+2*fLineDistance; y2++) {
                if (dataH[x][y2] > 0) {
                    lineAbove = true;
                }
            }
            if (lineAbove == false) values[5] = 1;
           cout << "                                        values[5] = " << values[5] << endl;
        }
        
        for (int ySearch = currentY-1; ySearch <= currentY+1; ySearch++) {
            
            if (dataH[x][ySearch] > 0 && dataV[x][ySearch] > 0) {
                currentY = ySearch;
                i++;
                l++;                
                intNumberXY[l].push_back(x);
                intNumberXY[l].push_back(ySearch);
                intNumberXY[l].push_back(i);
                intNumberXY[l].push_back(j);
                line = true;
                //cout << "IntersectionRight " << l << ": [" << intNumberXY[l][0] << "," << intNumberXY[l][1] << "," << intNumberXY[l][2] << "," << intNumberXY[l][3] << "]" << endl;
                break;
            }
            else if (dataH[x][ySearch] > 0) {
                currentY = ySearch;
                line = true;
                break;
            }
        }
        if (line == false) break;
    }
    
    values[2] = l;
    
    cout << "values (i="<< i << ", j=" << j << ", l=" << l << ", xInit=" << xInit << ", yInit=" << yInit << ")"  << endl;
} */