{
gROOT->Reset();

// open the files

gStyle-> SetOptStat(1);


// TFile* file = new TFile("Probe1.root");
TFile* file = new TFile("Probe1.root");

float crystalWidthY = 53.41;
float crystalWidthX = 53.41;

file->cd();

TTree* singlestree = (TTree*)gDirectory->Get("Singles");
TTree* hitstree = (TTree*)gDirectory->Get("Hits");
TTree* opticalsinglestree = (TTree*)gDirectory->Get("OpticalSingles");

// **** Array

int     maxNumberOfSingles = 100000;

int     APD[25][25][maxNumberOfSingles];
int   truePosition[3][maxNumberOfSingles];

for(int i = 0; i < 25 ; i++)
{
    for(int j = 0; j < 25 ; j++)
    {
        for(int k = 0; k < maxNumberOfSingles ; k++)
        {
            APD[i][j][maxNumberOfSingles] = 0;
        }
    }
}

int y , x ;

// ****

// *****************************
// *** hits of optical Photons***
// *****************************

Float_t  localPosX;
Float_t  localPosY;
Float_t  localPosZ;
Float_t  posZ;
Float_t  posX;
Float_t  posY;
Float_t  sourcePosZ_hits;
Float_t  sourcePosX_hits;
Float_t  sourcePosY_hits;
Int_t runID_hits;
Int_t eventID_hits;
Int_t sourceID_hits;
Int_t trackID;
Int_t parentID;



// *****************************
// *** hard Photons ***
// *****************************


Float_t  energy;
Float_t  globalPosX;
Float_t  globalPosY;
Float_t  globalPosZ;
Float_t  sourcePosZ;
Float_t  sourcePosX;
Float_t  sourcePosY;
Int_t runID;
Int_t eventID;
Int_t sourceID;


// *****************************
// *** Optical Photons ***
// *****************************

Float_t energy_optical;

Float_t  globalPosX_optical;
Float_t  globalPosY_optical;
Float_t  globalPosZ_optical;
Float_t  sourcePosZ_optical;
Float_t  sourcePosX_optical;
Float_t  sourcePosY_optical;
Int_t runID_optical;
Int_t eventID_optical;
Int_t sourceID_optical;

hitstree->SetBranchAddress("localPosX",&localPosX);
hitstree->SetBranchAddress("localPosY",&localPosY);
hitstree->SetBranchAddress("localPosZ",&localPosZ);
hitstree->SetBranchAddress("posX",&posX);
hitstree->SetBranchAddress("posY",&posY);
hitstree->SetBranchAddress("posZ",&posZ);
hitstree->SetBranchAddress("sourcePosX",&sourcePosX_hits);
hitstree->SetBranchAddress("sourcePosY",&sourcePosY_hits);
hitstree->SetBranchAddress("sourcePosZ",&sourcePosZ_hits);
hitstree->SetBranchAddress("runID", &runID_hits);
hitstree->SetBranchAddress("eventID", &eventID_hits);
hitstree->SetBranchAddress("sourceID", &sourceID_hits);
hitstree->SetBranchAddress("trackID", &trackID);
hitstree->SetBranchAddress("parentID", &parentID);


singlestree->SetBranchAddress("globalPosX",&globalPosX);
singlestree->SetBranchAddress("globalPosY",&globalPosY);
singlestree->SetBranchAddress("globalPosZ",&globalPosZ);
singlestree->SetBranchAddress("sourcePosX",&sourcePosX);
singlestree->SetBranchAddress("sourcePosY",&sourcePosY);
singlestree->SetBranchAddress("sourcePosZ",&sourcePosZ);
singlestree->SetBranchAddress("runID", &runID);
singlestree->SetBranchAddress("eventID", &eventID);
singlestree->SetBranchAddress("sourceID", &sourceID);
singlestree->SetBranchAddress("energy", &energy);


opticalsinglestree->SetBranchAddress("globalPosX",&globalPosX_optical);
opticalsinglestree->SetBranchAddress("globalPosY",&globalPosY_optical);
opticalsinglestree->SetBranchAddress("globalPosZ",&globalPosZ_optical);
opticalsinglestree->SetBranchAddress("sourcePosX",&sourcePosX_optical);
opticalsinglestree->SetBranchAddress("sourcePosY",&sourcePosY_optical);
opticalsinglestree->SetBranchAddress("sourcePosZ",&sourcePosZ_optical);
opticalsinglestree->SetBranchAddress("runID", &runID_optical);
opticalsinglestree->SetBranchAddress("eventID", &eventID_optical);
opticalsinglestree->SetBranchAddress("sourceID", &sourceID_optical);
opticalsinglestree->SetBranchAddress("energy", &energy_optical);


Int_t  nhits = hitstree-> GetEntries();
Int_t  nsingles  = singlestree-> GetEntries();
Int_t  nopticalsingles  = opticalsinglestree-> GetEntries();

cout << "Number of hard Photons: " << nsingles << endl;
cout << "Number of optical Photons: " << nhits << endl;

//----------------------------------------------------------

int i = 0;

// load data of first gamma
singlestree->GetEntry(i);
truePosition[0][0] = globalPosX;
truePosition[1][0] = globalPosY;
truePosition[2][0] = globalPosZ; // DOI

//run though all optical photons
for (k=0; k<nhits; k++)
{
    // get next data of next photon
    hitstree->GetEntry(k);

    // if this photon belongs to the actual gamma
    if (runID_hits == runID  &&   eventID_hits == eventID)
    {
        if (energy>0.5)
        {
            x= (posX+crystalWidthX/2)/(crystalWidthX/25);
            y= (posY+crystalWidthY/2)/(crystalWidthY/25); if(y > 24) y = 24;

            if(x<25 && y<25) APD[ x ][ y ][ i ] ++;
        }
    }
    else  // if photon belongs to next gamma
    {
        // get data of next gamma and fill strings again with the new gamma data
        i++;
        singlestree->GetEntry(i);

        truePosition[0][i] = globalPosX;
        truePosition[1][i] = globalPosY;
        truePosition[2][i] = globalPosZ;
    }
}

// Output

// Monitor

for(int k = 0; k <= i ; k++)
{
    std:: cout << "\t" << truePosition[0][k] << "\t" << truePosition[1][k] << "\t" << truePosition[2][k] << std:: endl;

    for(int x = 0; x < 25 ; x++)
    {
        for(int y = 0; y < 25 ; y++)
        {
            std:: cout << "\t" << APD[x][y][k];
        }
        std:: cout << std::endl;
    }
    std:: cout << std::endl << std::endl << std::endl;
}

// END Monitor

// ASCII

std::ofstream outfile;
outfile.open("Probe2---.dat");

for(int k = 0; k <= i ; k++)
{
    outfile << truePosition[0][k] << "\t" << truePosition[1][k] << "\t" << truePosition[2][k];

    for(int x = 0; x < 25 ; x++)
    {
        for(int y = 0; y < 25 ; y++)
        {
            outfile << "\t" << APD[x][y][k];
        }
    }
    outfile << std::endl;
}
outfile.close();


// ENDE ASCII
// ENDE: Output */

}