// // ******************************************************************** // * License and Disclaimer * // * * // * The GAMOS software is copyright of the Copyright Holders of * // * the GAMOS Collaboration. It is provided under the terms and * // * conditions of the GAMOS Software License, included in the file * // * LICENSE and available at http://fismed.ciemat.es/GAMOS/license .* // * These include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GAMOS collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the GAMOS Software license. * // ******************************************************************** // #include "GmSqdose.hh" #include "GamosCore/GamosUtils/include/GmGenUtils.hh" #include "GamosCore/GamosBase/Base/include/GmParameterMgr.hh" //----------------------------------------------------------------------- GmSqdose::GmSqdose() { theHeader = 0; G4String stype = GmParameterMgr::GetInstance()->GetStringValue("GmSqdose:FileType","ALL"); if( stype == "ALL" || stype == "all") { theType = SqTALL; } else if( stype == "FILLED" || stype == "filled") { theType = SqTFILLED; } else { G4Exception("GmSqdose::GmSqdose", "Wrong FileType", FatalErrorInArgument, G4String("FileType can be ALL/all or FILLED/filled, it is " + stype ).c_str()); } } //----------------------------------------------------------------------- GmSqdose::~GmSqdose() { // delete theHeader; } //----------------------------------------------------------------------- void GmSqdose::Print( FILE* out ) { theHeader->Print(out); if(fwrite(&theType, sizeof(size_t),1,out)!=1){ G4Exception("GmSqdose::Print", "Error", FatalException, "Error writing type"); } if( theType == SqTALL ) { size_t siz = theHeader->GetNoVoxelX()* theHeader->GetNoVoxelY()* theHeader->GetNoVoxelZ(); //----- Loop to all voxels and print dose for( size_t ii = 0; ii < siz; ii++ ) { float dose = theDoses[ii]; if(fwrite((float *) &dose, sizeof(float),1,out)!=1) { G4Exception("GmSqdose::Print", "Error", FatalException, "Error writing dose"); } } for( size_t ii = 0; ii < siz; ii++ ) { float dose2 = theDoseSqs[ii]; if(fwrite((float *) &dose2, sizeof(float),1,out)!=1){ G4Exception("GmSqdose::Print", "Error", FatalException, "Error writing doseSq"); } } } else if( theType == SqTFILLED ) { size_t siz = theDoses.size(); for( size_t ii = 0; ii < siz; ii++ ) { float dose = theDoses[ii]; float dose2 = theDoseSqs[ii]; if( dose != 0. ) { if(fwrite(&ii, sizeof(size_t),1,out)!=1){ G4Exception("GmSqdose::Print", "Error", FatalException, "Error writing index"); } if(fwrite((float *) &dose, sizeof(float),1,out)!=1){ G4Exception("GmSqdose::Print", "Error", FatalException, "Error writing dose"); } if(fwrite((float *) &dose2, sizeof(float),1,out)!=1){ G4Exception("GmSqdose::Print", "Error", FatalException, "Error writing doseSq"); } // G4cout << " WRITING DOSE " << ii << " " << dose << " " << dose2 << G4endl; } } } } //----------------------------------------------------------------------- void GmSqdose::Read( G4String fileName ) { FILE * fin = fopen(fileName,"rb"); Read(fin); } //----------------------------------------------------------------------- void GmSqdose::Read( FILE* fin ) { if( theHeader != 0 ){ G4Exception("GmSqdose::Read", "Error", FatalException, "This dose does have a header already"); } theHeader = new GmSqdoseHeader; theHeader->Read( fin ); size_t nv = theHeader->GetNoVoxelX() * theHeader->GetNoVoxelY() * theHeader->GetNoVoxelZ(); if( fread(&theType, sizeof(size_t), 1, fin) != 1) { G4Exception("GmSqdose::Read", "Error", FatalException, "Problem reading type "); } G4cout << " GmSqdose::Read type " << theType << G4endl; if( theType == 1 ) { for( size_t iv = 0; iv < nv; iv++ ){ float ftmp; if( fread(&ftmp, sizeof(float), 1, fin) != 1) { G4Exception("GmSqdose::Read()", "Problem reading dose ",FatalErrorInArgument,G4String("Reading voxel number "+GmGenUtils::itoa(iv)).c_str()); } theDoses.push_back( ftmp ); // if( theDoses.size()%100000 == 1 ) G4cout << iv << " READ dose " << theDoses.size() << " = " << theDoses[theDoses.size()-1] << G4endl; } for( size_t iv = 0; iv < nv; iv++ ){ float ftmp; if( fread(&ftmp, sizeof(float), 1, fin) != 1) { G4Exception("GmSqdose::Read()","Problem reading dose ",FatalErrorInArgument,G4String("Reading voxel number "+GmGenUtils::itoa(iv)).c_str()); } theDoseSqs.push_back( ftmp ); // if( theDoseSqs.size()%100000 == 1 ) G4cout << iv << " READ dose2 " << theDoseSqs.size() << " = " << theDoseSqs[theDoseSqs.size()-1] << G4endl; } } else if( theType == 2 ){ G4int idPrev = -1; size_t iv; for( iv = 0; iv < nv; iv++ ){ size_t idNow; if( fread(&idNow, sizeof(size_t), 1, fin) != 1) { if(feof(fin)) { break; } else { G4Exception("GmSqdose::Read()","Problem reading voxel ID ",FatalErrorInArgument,G4String("Reading voxel number "+GmGenUtils::itoa(idNow)).c_str()); } } float ftmp; if( fread(&ftmp, sizeof(float), 1, fin) != 1) { G4Exception("GmSqdose::Read()","Problem reading dose ",FatalErrorInArgument,G4String("Reading voxel number "+GmGenUtils::itoa(idNow)).c_str()); } for( size_t ii = 0; ii < idNow-idPrev-1; ii++ ){ theDoses.push_back( 0. ); } theDoses.push_back( ftmp ); if( fread(&ftmp, sizeof(float), 1, fin) != 1) { G4Exception("GmSqdose::Read()","Problem reading dose ",FatalErrorInArgument,G4String("Reading voxel number "+GmGenUtils::itoa(idNow)).c_str()); } for( size_t ii = 0; ii < idNow-idPrev-1; ii++ ){ theDoseSqs.push_back( 0. ); } theDoseSqs.push_back( ftmp ); // G4cout << iv << " READ DOSE " << idNow << " prev " << idPrev << " dose " << ftmp1 << " dose2 " << ftmp << G4endl; idPrev = idNow; // if( theDoseSqs.size()%100000 == 1 ) G4cout << iv << " READ dose2 " << theDoseSqs.size() << " = " << theDoseSqs[theDoseSqs.size()-1] << G4endl; } G4cout << " GmSqdose Number of voxels read " << iv << " out of " << nv << G4endl; } } //----------------------------------------------------------------------- GmSqdose GmSqdose::operator+=( const GmSqdose& gmdose ) { theHeader->SetNumberOfEvents(theHeader->GetNumberOfEvents()+gmdose.GetHeader()->GetNumberOfEvents()); std::vector doses = gmdose.GetDoses(); std::vector doses2 = gmdose.GetDoseSqs(); size_t siz = theDoses.size(); if( doses.size() > siz ) { for( size_t ii = siz; ii < doses.size(); ii++ ) { theDoses[ii] = 0.; } } // G4double nevA = theHeader->GetNumberOfEvents(); // G4double nevB = gmdose.GetHeader()->GetNumberOfEvents(); // G4double nevSUM = nevA + nevB; // G4cout << " nevSUM " << nevSUM << " nevA " << nevA << " nevB " << nevB << G4endl; for( size_t ii = 0; ii < siz; ii++ ) { // if( ii == 0 ) G4cout << this << " dose1 " << theDoses[ii] << " " << &gmdose << " dose2 " << dose2[ii] << G4endl; // theDoses[ii] = (theDoses[ii]*nevA + doses[ii]*nevB) / nevSUM; // theDoseSqs[ii] = (theDoseSqs[ii]*nevA*nevA + doses2[ii]*nevB*nevB)/(nevSUM*nevSUM); if( ii < doses.size() ) { theDoses[ii] = theDoses[ii] + doses[ii]; theDoseSqs[ii] = theDoseSqs[ii] + doses2[ii]; } // if( ii == 0 ) G4cout << "summed dose1 " << theDoses[ii] << " dose2 " << doses[ii] << G4endl; // if( ii == 0) G4cout << "summed doseSq1 " << theDoseSqs[ii] << " dose2 " << doses2[ii] << G4endl; } // G4cout << "dose value sum " << theDoses[0] << G4endl; return *this; } //-------------------------------------------------------------------- void GmSqdose::CalculateErrors() { size_t nvox = theDoses.size(); G4double nEvents = theHeader->GetNumberOfEvents(); for( size_t ii = 0; ii < nvox; ii++ ){ G4double error = (theDoseSqs[ii]*nEvents - theDoses[ii]*theDoses[ii]) / (nEvents-1); if( error < 0. ) { if( error < -1.E-30 ) G4cerr << " !!WARNING GmSqdose::CalculateErrors negative error, set to 0. " << error << G4endl; error = 0.; } else { error = std::sqrt(error)/nEvents; } theErrors.push_back( error ); // G4cout << " GetError " << ii q<< " e= " << error << " S1 " << theDoseSqs[ii]*nEvents << " s2 " << theDoses[ii]*theDoses[ii] << " S3 " << (nEvents*nEvents*(nEvents-1)) << G4endl; } }