// ESAF : Euso Simulation and Analysis Framework
// $Id: DatabaseOpticalSystem.cc,v 1.16 2005/05/17 21:48:30 thea Exp $
// A.Thea, J.Watts created Mar, 23 2004

#include "DatabaseOpticalSystem.hh"
#include "EsafRandom.hh"

ClassImp(DatabaseOpticalSystem)

// ctor
 DatabaseOpticalSystem::DatabaseOpticalSystem() {
    fKeyLength = (size_t)Conf()->GetNum("DatabaseOpticalSystem.fKeyLength");

    string name = Conf()->GetStr("DatabaseOpticalSystem.fBinFilePath");
    fBinFile = new ifstream(name.c_str());

    if (!fBinFile->is_open())
        throw runtime_error("DatabaseOpticalSystem: unable to open "+name);

    Init();
}

// dtor
 DatabaseOpticalSystem::~DatabaseOpticalSystem() {
    if ( !fBinFile) {
        fBinFile->close();
        delete fBinFile;
    }
}

// transport a photon through the optics
 Photon* DatabaseOpticalSystem::Transport( Photon *p ) const {
    // check the photon to be on the first lens
    if ( (FirstLensTop()-p->pos[Z]) < -kTolerance )
       return 0;
   
    // track oly photons going up 
    if ( p->dir.Unit()[Z] < kTolerance )
        return 0;

    if ( p->dir.Theta() < fAngles.front() || p->dir.Theta() > fAngles.back() ) {
        cout << "out of angle" << endl;
        return 0;
    }
    
    // find the nearest wl
    Double_t dWl = TMath::Abs(p->wl-fWavelengths.front());
    Int_t iWl(0);
    for (size_t k(0); k<fWavelengths.size(); k++)
        if ( TMath::Abs(p->wl-fWavelengths[k]) < dWl ){
            dWl = TMath::Abs(p->wl-fWavelengths[k]);
            iWl = k;
        }    
     
    // find the two nearest angles
    Double_t dAng = TMath::Abs(p->dir.Theta()-fAngles.front());
    Int_t iAngA(0), iAngB(0);
    for (size_t k(0); k<fAngles.size(); k++)
        if ( TMath::Abs(p->dir.Theta()-fAngles[k]) < dAng ){
            dAng = TMath::Abs(p->dir.Theta()-fAngles[k]);
            iAngA = k;
        }

    if ( TMath::Abs(fAngles[iAngA+1]-p->dir.Theta()) < TMath::Abs(fAngles[iAngA-1]-p->dir.Theta()) )
        iAngB = iAngA+1;
    else
        iAngB = iAngA-1;

    // trace photon back to the entrance disc
    EVector olddir = p->dir.Unit();
    EVector dummy = p->pos-olddir*(p->pos[Z]-fEntranceZ);
     
    if ( dummy.Perp() > fEntranceDiscRadius )
       return 0; 

    // rotate dummy in the database reference system
    dummy.SetPhi(dummy.Phi()+TMath::PiOver2()-olddir.Phi());
    Bool_t isMirror = dummy[X] < 0;
    if ( isMirror ) dummy[X] = -dummy[X];
    
    Double_t yOffset = fXYStep*(fNumCols.size()-1)/2;
    
    // find the nearest spot inside EntranceDiscRadius
    Int_t iCol, jRow;
  
    iCol = (Int_t) TMath::Floor(dummy[X]/fXYStep); 
    if ( dummy[X]>(iCol+0.5)*fXYStep ) iCol++;

    jRow = (Int_t)(TMath::Floor((yOffset-dummy[Y])/fXYStep)); 
    if ( yOffset - dummy[Y]>(jRow+0.5)*fXYStep ) jRow++;
    
    Double_t xCol = iCol*fXYStep;
    Double_t yRow = yOffset-jRow*fXYStep;
    if ( TMath::Sqrt(xCol*xCol+yRow*yRow) > fEntranceDiscRadius ){
        cout << "out of radius" << endl;
        cout << TMath::Sqrt(xCol*xCol+yRow*yRow) << " xCol  = " << xCol << " yRow = " << yRow << endl;
        return 0;
    }
  
    Float_t photonA[5], photonB[5], photon[5];
    streampos sA, sB;
    sA = fKeyMap[iWl][iAngA]+(fRowsOffset[jRow]+(streampos)iCol)*5*sizeof(Float_t);
    sB = fKeyMap[iWl][iAngB]+(fRowsOffset[jRow]+(streampos)iCol)*5*sizeof(Float_t);

    fBinFile->seekg(sA, ios::beg
    fBinFile->read((char*)photonA, sizeof(Float_t)*5);
    fBinFile->seekg(sB, ios::beg
    fBinFile->read((char*)photonB, sizeof(Float_t)*5);
    
    // merge the two db entries
    Float_t weightA, weightB;
    weightA = 1-(p->dir.Theta()-fAngles[iAngA])/(fAngles[iAngB]-fAngles[iAngA]); // between 0 and 1
    weightB = 1-weightA;

    if ( photonA[0] < 0 && photonB[0] < 0 )
        return 0;
    else if ( photonA[0] < 0 ) {
        weightA = 0;
    } else if ( photonB[0] < 0 ) {
        weightB = 0;
    } else {
    }
    if ( photonA[0] < 0 || photonB[0] < 0 )
        return 0;

    // merge the two db entries
    photon[0]=weightA*photonA[0]+weightB*photonB[0];
    photon[1]=weightA*photonA[1]+weightB*photonB[1];
    photon[2]=weightA*photonA[2]+weightB*photonB[2];
    photon[3]=weightA*photonA[3]+weightB*photonB[3];
    photon[4]=weightA*photonA[4]+weightB*photonB[4];


    
    // check if the ph gets through the optics
    TRandom *rndm = EsafRandom::Get();
    if ( photon[0] < 0 || rndm->Rndm() > photon[0] )
        return 0;

    p->pos[X] = isMirror ? -photon[1] : photon[1];
    p->pos[Y] = photon[2];
    p->pos[Z] = fExitZ;
    p->dir.SetTheta(photon[3]);
    p->dir.SetPhi(photon[4]);
    p->dir[X] = isMirror ? -p->dir[X] : p->dir[X];
    // inverse rotation
    p->pos.SetPhi(p->pos.Phi()-TMath::PiOver2()+olddir.Phi());
    p->dir.SetPhi(p->dir.Phi()-TMath::PiOver2()+olddir.Phi());
    return p;
}

// initialize the optics module 
 void DatabaseOpticalSystem::Init() {
    cout << "DatabaseOpticalSystem: Initializing..." << flush;
    fBinFile->seekg(0,ios::beg
    streampos cur = fBinFile->tellg();
    cur = LoadHeader( cur );
    cur = LoadKeyMap( cur );
    cout << " done." << endl;
}

 streampos DatabaseOpticalSystem::LoadHeader( streampos start ) {

    fBinFile->seekg(start);
    streampos cur = fBinFile->tellg();

    Float_t *values(NULL);
    char buffer[fKeyLength];
    
    // Open header
    fBinFile->read(buffer, fKeyLength);
    if ( (KeyWord(kBegin)+KeyWord(kHeader)) != buffer) {
        cout << buffer << endl;
        throw runtime_error("Header missing or misplaced. Exiting...");
}
    
    size_t headSize;
    fBinFile->read((char*)&headSize, sizeof(size_t));
   
    // check the end of the header to be where it should
    cur = fBinFile->tellg();
    fBinFile->seekg(headSize, ios::cur

    fBinFile->read(buffer, fKeyLength);
    if ( (KeyWord(kEnd)+KeyWord(kHeader)) != buffer) {
        cout << buffer << endl;
        throw runtime_error("End Header missing or misplaced. Exiting...");
}
    
    fBinFile->seekg(cur);

    // read comment
    size_t commentSize;
    char* comment;

    fBinFile->read((char*)&commentSize, sizeof(size_t));
    comment = new char[commentSize];
    fBinFile->read((char*)comment, commentSize);
    delete [] comment;

    Int_t numData;
    fBinFile->read((char*)&numData, sizeof(Int_t));

    values = new Float_t[numData];
    fBinFile->read((char*)values, sizeof(Float_t)*numData);
    /* values now contains the Physical Parameters of the detector in the
     * following order:
     *
     * XYStep
     * OpticsRadius
     * EntranceDiscRadius
     * FirstLensBottom
     * FirstLensTop
     * SecondLensBottom
     * SecondLensTop
     * EntranceZ
     * ExitZ
     * PosZ
     */

    fXYStep = values[0];
    fR = values[1];
    fEntranceDiscRadius = values[2];
    fFirstLensDZ = values[4]-values[3];
    fSecondLensDZ = values[6]-values[5];
    fEntranceZ = values[7];
    fExitZ = values[8];
    fPos[Z] = values[9];

    fDZdown = 0;
    fDZup = values[6]-values[3];

    delete [] values;

    // read wls
    Int_t numWls;
    fBinFile->read((char*)&numWls, sizeof(Int_t));

    values = new Float_t[numWls];
    fBinFile->read((char*)values, sizeof(Float_t)*numWls);
    for(Int_t i(0); i<numWls; i++) {
        fWavelengths.push_back(values[i]);
    }
    delete [] values;
    
    // read angles
    Int_t numAngles;
    fBinFile->read((char*)&numAngles, sizeof(Int_t));

    values = new Float_t[numAngles];
    fBinFile->read((char*)values, sizeof(Float_t)*numAngles);
    for(Int_t i(0); i<numAngles; i++) { 
        fAngles.push_back(values[i]);
    }
    delete [] values;

    // positions
    fBinFile->read((char*)&fNumPos, sizeof(Int_t));

    Int_t rows;
    Int_t *cols(NULL);
    fBinFile->read((char*)&rows, sizeof(Int_t));

    Int_t countPos(0);
    cols = new Int_t[rows];
    fBinFile->read((char*)cols, sizeof(Float_t)*rows);
    for(Int_t i(0); i<rows; i++) {
        fRowsOffset.push_back(countPos);
        fNumCols.push_back(cols[i]);
        countPos+=fNumCols[i];
    }
    delete [] cols;

    if ( fNumPos != countPos ) {
        string dummy;
        dummy = "Mismatch between number of columns and total "; 
        dummy += "number of positions in the header. Exiting...";
        throw runtime_error(dummy);
    }

    // Header closed
    fBinFile->read(buffer, fKeyLength);
    if ( (KeyWord(kEnd)+KeyWord(kHeader)) != buffer) {
        throw runtime_error("End Header not found. Something's gone wrong. Exiting...");
    }

    return fBinFile->tellg();
}

 streampos DatabaseOpticalSystem::LoadKeyMap( streampos start) {
    
    fBinFile->seekg(start);
    streampos cur = fBinFile->tellg();

    // Photon list found
    char buffer[fKeyLength];
    fBinFile->read(buffer, fKeyLength);
    if ( (KeyWord(kBegin)+KeyWord(kPhotons)) != buffer) {
        cout << buffer << endl;
        throw runtime_error("Photon missing or misplaced. Exiting");
    }

    // calculate the length of the photons' list
    size_t angleSize = sizeof(char)*(2*fKeyLength)+sizeof(Int_t)+sizeof(Float_t)*(5*fNumPos);
    size_t wlsize = sizeof(char)*(2*fKeyLength)+sizeof(Int_t)+angleSize*fAngles.size();
    size_t phSize = wlsize*fWavelengths.size();

    cout << "Calculated list size " << phSize << endl;
    
    // check the end of the list to be where it should
    cur = fBinFile->tellg();
    fBinFile->seekg(phSize, ios::cur

    // End Photons found
    fBinFile->read(buffer, fKeyLength);
    if ( (KeyWord(kEnd)+KeyWord(kPhotons)) != buffer) {
        cout << buffer << endl;
        throw runtime_error("End Photons missing or misplaced. Exiting...");
    }

    fBinFile->seekg(cur);
    
    // checking all the key positions
    fKeyMap.clear();
    for ( Int_t iWl(0); iWl < (Int_t)fWavelengths.size(); iWl++) {
        // check the wl key
        fBinFile->read(buffer, fKeyLength);
        if ( (KeyWord(kBegin)+KeyWord(kWavelength)) != buffer) {
            cout << buffer << endl;
            throw runtime_error(Form("Wavelength keyword %d not found. Exiting", iWl));
        }
        
        // check wl id
        Int_t wlId;
        fBinFile->read((char*)&wlId, sizeof(Int_t));
        if ( wlId != iWl)
            throw runtime_error(Form("Wavelength id mismatch: wlId = %d   iWl = %d", wlId, iWl));
       
        vector<streampos> dummy;
        for ( Int_t iAng(0); iAng < (Int_t)fAngles.size(); iAng++) {
            // check the angle key
            fBinFile->read(buffer, fKeyLength);
            if ( (KeyWord(kBegin)+KeyWord(kAngle)) != buffer) {
                cout << buffer << endl;
                throw runtime_error(Form("Angle keyword %d not found. Exiting", iAng));
            }

            // check angle id 
            Int_t angleId;
            fBinFile->read((char*)&angleId, sizeof(Int_t));
            if ( angleId != iAng)
                throw runtime_error(Form("Angle id mismatch: angleId = %d   iAng = %d", wlId, iWl));
            
            // keep the pos of the first ph
            cur = fBinFile->tellg();
            
            // move at the end
            fBinFile->seekg(sizeof(Float_t)*5*fNumPos,ios::cur
            
            // check the end key
            fBinFile->read(buffer, fKeyLength);
            if ( (KeyWord(kEnd)+KeyWord(kAngle)) != buffer) {
                cout << buffer << endl;
                throw runtime_error(Form("End Angle keyword %d not found. Exiting...", iAng));
            }

            // save the streamer position of the first ph
            dummy.push_back(cur);
        }
        
        // check the end key
        fBinFile->read(buffer, fKeyLength);
        if ( (KeyWord(kEnd)+KeyWord(kWavelength)) != buffer) {
            cout << buffer << endl;
            throw runtime_error(Form("End Wavelength keyword %d not found. Exiting...",iWl));
        }
        
        fKeyMap.push_back(dummy); 
    }

    return fBinFile->tellg();
}