// ESAF : Euso Simulation and Analysis Framework
// $Id: DatabaseOpticalSystem.cc,v 1.15 2005/05/13 15:23:55 thea Exp $
// A.Thea, J.Watts created Mar, 23 2004
#include "DatabaseOpticalSystem.hh"
#include "EsafRandom.hh"
// SLC gcc 3.2.3 workaround for ifstream/TString conflict
char* operator+( std::streampos&, char* );
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();
}