// // ******************************************************************** // * 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 "globals.hh" #include #include "RadioTherapy/include/iaea_header.hh" #include "RadioTherapy/include/iaea_record.hh" //#include "RadioTherapy/include/iaeaUtilities.hh" #include "RadioTherapy/include/RTPhaseSpaceHistos.hh" #include "RadioTherapy/include/EGSPhspParticle.hh" #include "GamosCore/GamosUtils/include/GmGenUtils.hh" #include "GamosCore/GamosUtils/include/GmNumberOfEvent.hh" #include "GamosCore/GamosBase/Base/include/GmAnalysisMgr.hh" #include "GamosCore/GamosBase/Base/include/GmParameterMgr.hh" #include "G4UImanager.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void printStat( iaea_header_type *p_iaea_header1 ); void GetError( G4double& rat, G4double& err, G4double val1, G4double val2 ); //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... int main(int argc,char** argv) { if(argc<3) { G4Exception("analysePS_EGS", "Error", FatalErrorInArgument, "!!! FATAL ERROR: YOU MUST SUPPLY AT LEAST TWO ARGUMENTS: PHASESPACE_FILE PHASESPACE_Z"); } G4String theFileName; G4String hisEMax = "10.*MeV"; G4String hisRMax = "100.*mm"; G4String hisNbins = "100"; G4String theFileNameOut = "phaseSpace"; G4int nPartLoop = -1; G4double thePosZ = DBL_MAX; if( argc == 2 ) { theFileName = argv[1]; } else { if(argc%2 != 1) { G4Exception("analysePS_EGS", "Error", FatalErrorInArgument, "WRONG NUMBER OF ARGUMENTS: THEY MUST BE -XX1 VAL_XX1 -XX2 VAL_XX2 ... "); } for( G4int ii = 1; ii < argc; ii+=2 ){ G4String argvstr = argv[ii]; if( argvstr == "-f" ) { theFileName = argv[ii+1]; } else if( argvstr == "-posZ" ) { thePosZ = GmGenUtils::GetValue(argv[ii+1]); } else if( argvstr == "-EMax" ) { hisEMax = argv[ii+1]; } else if( argvstr == "-RMax" ) { hisRMax = argv[ii+1]; } else if( argvstr == "-NBins" ) { hisNbins = argv[ii+1]; } else if( argvstr == "-fOut" ) { theFileNameOut = argv[ii+1]; } else if( argvstr == "-NPart" ) { nPartLoop = GmGenUtils::GetInt(argv[ii+1]); } } } G4cout << " READING FILE " << theFileName << G4endl; FILE* phspfile = fopen(theFileName.c_str(), "rb"); if( thePosZ == DBL_MAX) { G4Exception("analysePS_EGS", "Error", FatalErrorInArgument, "YOU MUST SUPPLY ALWAYS -posZ ARGUMENT"); } char mode[6]; fread(mode,5*sizeof(char),1,phspfile); mode[5]='\0'; int mode_n; if(mode[4] == '0') mode_n = 0; else mode_n = 2; std::cout << " MODE = " << mode[4] << " => " << mode_n << std::endl; // int matrixInd1, matrixInd2; //these are the matrix indexes iaea_header_type *p_iaea_header; p_iaea_header = (iaea_header_type *) calloc(1, sizeof(iaea_header_type)); fread(&(p_iaea_header->nParticles),sizeof(int),1,phspfile); fread(&(p_iaea_header->particle_number[0]),sizeof(int),1,phspfile); std::cout << " N_PARTICLES " << p_iaea_header->nParticles << " N_PHOTONS " << p_iaea_header->particle_number[0] << std::endl; fread(&(p_iaea_header->maximumKineticEnergy[0]),sizeof(float),1,phspfile); fread(&(p_iaea_header->minimumKineticEnergy[0]),sizeof(float),1,phspfile); float p_incident; fread(&p_incident,sizeof(float),1,phspfile); std::cout << " MAX_KINE " << p_iaea_header->maximumKineticEnergy[0] << " MIN_KINE " << p_iaea_header->minimumKineticEnergy[0] << " P_INCIDENT " << p_incident << std::endl; /* fread(&n_of_particles,sizeof(int),1,phspfile); fread(&n_of_photons,sizeof(int),1,phspfile); fread(&energy_max,sizeof(float),1,phspfile); fread(&energy_min,sizeof(float),1,phspfile); fread(&p_incident,sizeof(float),1,phspfile); */ // long file_pointer = 5*sizeof(char)+2*sizeof(int)+3*sizeof(float)+3; int count_particles = 0; RTPhaseSpaceHistos* theHistos = 0; GmParameterMgr::GetInstance()->PrintParametersUsage(2); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo1Max PhaseSpace*Energy* ")+hisEMax); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo2MaxY PhaseSpace*Energy* ")+hisEMax); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo1Max PhaseSpace*X*at* ")+hisRMax); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo1Max PhaseSpace*Y*at* ")+hisRMax); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo1Max PhaseSpace*Z*at* ")+hisRMax); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo2MaxX PhaseSpace*XY*at* ")+hisRMax); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo2MaxY PhaseSpace*XY*at* ")+hisRMax); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo2MaxX PhaseSpace*R*vs* ")+hisRMax); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo2MaxX PhaseSpace*X*vs* ")+hisRMax); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo1Nbins PhaseSpace* ")+hisNbins); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo2NbinsX PhaseSpace* ")+hisNbins); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/analysis/histo2NbinsY PhaseSpace* ")+hisNbins); G4UImanager::GetUIpointer()->ApplyCommand(G4String("/gamos/setParam RTPhaseSpaceHistos:FileName ")+theFileNameOut); G4UImanager::GetUIpointer()->ApplyCommand("/gamos/setParam phaseSpace:NormalizeToNEvents 0"); EGSPhspParticle* particle = new EGSPhspParticle(phspfile,mode); int nParticles = p_iaea_header->nParticles; if( nPartLoop == 1 ) nPartLoop = nParticles; std::cout << " NPARTICLES TO BE READ " << nPartLoop << std::endl; for(long irecord=0; irecordRefreshParticle(phspfile,irecord+1); // std::cout << irecord << " PARTICLE " << particle->theLatch // << " " << particle->theEnergy << " " << particle->thePosX << " " << particle->thePosY << " " << particle->theDirU << " " << particle->theDirV << " " << particle->theWeight << std::endl; if(!iError){ /* float xcorr = (particle->thePosX - theXmin)/theDeltaX; float ycorr = (particle->thePosY - theYmin)/theDeltaY; if(particle->theCharge == 0){ count_photons++; }*/ count_particles++; } else { break; } if( !theHistos ) theHistos = new RTPhaseSpaceHistos(thePosZ ); theHistos->FillHistos( particle, thePosZ ); } G4cout << " NPARTICLES READ " << count_particles << G4endl; /* if(write==1){ fclose(phspfile); n_of_particles=count_particles; n_of_photons=count_photons; p_incident=n; } if(write==0){ mat -> save_matrix(firstname+"MATRIX"); mat -> colorscale("colors"); mat -> plot_matrix(firstname+"MATRIX","colors"); } */ // printStat( p_iaea_header ); GmAnalysisMgr::DeleteInstances(); GmParameterMgr::GetInstance()->CheckParametersUsage(); } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void printStat( iaea_header_type *p_iaea_header1 ) { long theNParticlesInFile = p_iaea_header1->orig_histories; G4cout << "ORIGINAL HISTORIES= " << theNParticlesInFile << G4endl; GmNumberOfEvent::SetNumberOfEvent(theNParticlesInFile); G4double val1, err1; GetError( val1, err1, p_iaea_header1->nParticles, theNParticlesInFile); G4cout << "PARTICLES PER EVENT= " << val1 << " +- " << err1 << G4endl; GetError( val1, err1, p_iaea_header1->particle_number[0], theNParticlesInFile); G4cout << "GAMMAS PER EVENT= " << val1 << " +- " << err1 << G4endl; GetError( val1, err1, p_iaea_header1->particle_number[1], theNParticlesInFile); G4cout << "ELECTRONS PER EVENT " << val1 << " +- " << err1 << G4endl; GetError( val1, err1, p_iaea_header1->particle_number[2], theNParticlesInFile); G4cout << "POSITRONS PER EVENT= " << val1 << " +- " << err1 << G4endl; } //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... void GetError( G4double& rat, G4double& err, G4double val1, G4double val2 ) { rat = val1/val2; if( val1 == 0 ) { err = 0.; } else { err = rat*sqrt(1./val1 + 1./val2 ); } // err *= sqrt( err * (1-err) ); }