// // ******************************************************************** // * 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. * // ******************************************************************** // // SHRadDecayScorerFromFlux //#define protected public //#define private public #include "SHRadDecayScorerFromFlux.hh" #include "SHUtils.hh" #include "ShieldingVerbosity.hh" #include "GamosCore/GamosBase/Base/include/GmParameterMgr.hh" #include "GamosCore/GamosUtils/include/GmNumberOfEvent.hh" #include "GamosCore/GamosUtils/include/GmG4Utils.hh" #include "GamosCore/GamosUtils/include/GmGenUtils.hh" #include "G4HadronicProcess.hh" #include "G4ParticleHPInelastic.hh" #include "G4ReactionProduct.hh" #include "G4Step.hh" #include "G4StepPoint.hh" #include "G4Track.hh" #include "G4VSolid.hh" #include "G4IonTable.hh" #include "G4ParticleTable.hh" //------------------------------------------------------------------- SHRadDecayScorerFromFlux::SHRadDecayScorerFromFlux() { bAdjustPHP = !getenv("G4PARTICLEHP_DO_NOT_ADJUST_FINAL_STATE"); bAdjustHP = !getenv("G4NEUTRONHP_DO_NOT_ADJUST_FINAL_STATE"); theEnergyInterpStep = GmParameterMgr::GetInstance()->GetNumericValue("SHRadDecayScorerFromFlux:EnergyInterpolationStep",DBL_MAX); theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable()); sumALL = 0.; bNewEvent = TRUE; theEvtMap = new G4THitsMap("SHRadDecayScorerFromFlux","SHRadDecayScorerFromFlux"); theRunMap = new G4THitsMap("SHRadDecayScorerFromFlux","SHRadDecayScorerFromFlux"); bInitialised = FALSE; } //------------------------------------------------------------------- SHRadDecayScorerFromFlux::~SHRadDecayScorerFromFlux() {;} //------------------------------------------------------------------- void SHRadDecayScorerFromFlux::ProcessHits(const G4Step* aStep) { if( !bInitialised ) { SHUtils::StoreHadronicProcesses(aStep->GetTrack()->GetDefinition()->GetParticleName()); bInitialised = TRUE; } G4cout << " SHRadDecayScorerFromFlux::ProcessHits " << G4endl;//GDEB G4double energyPre = aStep->GetPreStepPoint()->GetKineticEnergy(); G4double energyPost = aStep->GetPostStepPoint()->GetKineticEnergy(); G4LogicalVolume* lv = aStep->GetTrack()->GetVolume()->GetLogicalVolume(); G4Material* material = lv->GetMaterial(); G4double stepLength = aStep->GetStepLength(); G4double cubicVolume = lv->GetSolid()->GetCubicVolume(); size_t numberOfElem = material->GetNumberOfElements(); // G4double aTemp = material->GetTemperature(); G4ElementVector* elementVector = const_cast(material->GetElementVector()); const G4double* atomNumDensityVector = material->GetVecNbOfAtomsPerVolume(); //--- Loop to hadronic processes for( G4int iproc = 0; iproc < SHUtils::GetNHadronicProcesses(); iproc++ ){ G4HadronicProcess* hadProc = SHUtils::GetHadronicProcess(iproc); G4String hname = material->GetName() + "_" + hadProc->GetProcessName(); if( SHUtils::GetNChannelParticles() == 0 ) SHUtils::BuildChannelList(aStep->GetTrack()->GetDefinition(), "SHRadDecayScorerFromFlux"); //-- Inelastic G4HadronicInteraction* hadModel = hadProc->ChooseHadronicInteraction(energyPre, material, (*elementVector)[0]); G4cout << " hadModel " << hadModel->GetModelName() << G4endl;//GDEB // G4String modelName = hadModel->GetModelName(); if( dynamic_cast(hadModel) != 0 ) { G4ParticleHPInelastic* HPinelas = (G4ParticleHPInelastic*)(hadModel); G4cout << " SHRadDecayScorerFromFlux::ProcessHits G4ParticleHPInelastic nelem " << numberOfElem << G4endl;//GDEB //--- loop to elements for (size_t iele=0; ieleGetIndex(); G4ParticleHPChannel ** HPChannels = HPinelas->theInelastic[index]->theChannels; G4int nChannels = HPinelas->theInelastic[index]->nChannels; //--- loop to isotopes G4int numberOfIsos = HPChannels[0]->GetNiso(); // all channels have the same number of isotopes?? G4cout << " SHRadDecayScorerFromFlux::ProcessHits G4ParticleHPInelastic nchan " << nChannels << " nisot " << numberOfIsos << G4endl;//GDEB for( G4int ip = 0; ip < numberOfIsos; ip++ ) { // G4String isotName = elem->GetIsotope(ip)->GetName(); for(G4int ic=0; ictheFSType << G4endl; if( HPChannel->HasFSData(ip) ) {// one channel per isotope , one final state per isotope G4ParticleHPFinalState * FS = HPChannel->theFinalStates[ip]; //--- Get production for each ion G4ParticleHPEnAngCorrelation* energyAngData; G4ReactionProduct * aTarget = new G4ReactionProduct(); aTarget->SetMass(FS->theBaseA); //t if( dynamic_cast(FS) ) { // if( FSType == "F01/" || FSType == "F23/" || FSType == "F24/" || FSType == "F25/" || FSType == "F26/" || FSType == "F27/" ){ G4ParticleHPInelasticCompFS* FSC = (G4ParticleHPInelasticCompFS*)(FS); G4int it = FSC->SelectExitChannel(energyPre); energyAngData = FSC->theEnergyAngData[it]; if( !energyAngData ) continue; energyAngData->SetTarget(*aTarget); G4cout << FSC << " Use CompFS theEnergyAngData " << it << " " << FSC->theEnergyAngData[it] << G4endl; } else { G4ParticleHPInelasticBaseFS* FSB = (G4ParticleHPInelasticBaseFS*)(FS); energyAngData = FSB->theEnergyAngData; energyAngData->SetTarget(*aTarget); if( !energyAngData ) continue; G4cout << FSB << " Use BaseFS theEnergyAngData " << FSB->theEnergyAngData << G4endl; } if( !bAdjustPHP ) { G4int nProducts = energyAngData->nProducts; G4ParticleHPProduct * products = energyAngData->theProducts; for(G4int ip=0; ipGetMassCode() < 2004 ) continue; // only ions, bigger than alpha G4double ionProd = 0.; for( G4double ener = energyPre; ener < energyPost; ener += theEnergyInterpStep ) { // G4double dEdX = theEmCalculator->GetDEDX(ener, particleDef, material); G4double yield = product->GetMeanYield(ener); // Micro XS G4double FSXS = HPChannel->GetFSCrossSection(ener, ip); // Macro XS // G4double cross = FSXS*(nAtomsPermm3*isotFraction)*barn; G4double cross = FSXS*(nAtomsPermm3); // FSXS has into account the isotope fraction // G4double cross = FSXS*(nAtomsPermm3)*barn; // FSXS has into account the isotope fraction ionProd += yield * stepLength * cross / cubicVolume; #ifndef GAMOS_NO_VERBOSE if( SHVerb(debugVerb) && yield != 0 ) { G4String FSType = HPChannel->theFSType; G4String channelName = SHUtils::GetChannelName(int(FS->theBaseZ), int(FS->theBaseA), FSType.substr(0,3), GmG4Utils::GetParticleShortName(aStep->GetTrack()->GetDefinition()->GetParticleName()) ); G4String productName = SHUtils::GetProductName( product->theMassCode ); G4cout << material->GetName() << " " << hadProc->GetProcessName() << " CHANNEL= " << FSType << " " << channelName << " PRODUCT= " << productName << " YIELD= " << yield << " ION_PRODUCTION= " << ionProd << " XS micro " << FSXS << " XS macro " << cross << " Ener " << ener << " cross " << cross << " nAtomsPermm3 " << nAtomsPermm3 << G4endl; } #endif } // energy loop // Get index: G4int index = GetIndexFromProduct( product ); FillScorer( aStep, index, ionProd, aStep->GetPreStepPoint()->GetWeight()); // call directly the filling with index } // products loop } else { std::map theIonProds; std::map theAZs; for( G4double ener = energyPre; ener < energyPost; ener += theEnergyInterpStep ) { theAZs.clear(); // sample it N times and get all posible A-Z combinations G4int nSamples = 1000; for( G4int is = 0; is < nSamples; is++ ) { G4ReactionProductVector* prodVector = energyAngData->Sample(ener); G4ReactionProductVector::const_iterator itepv; G4int sumA = 0; G4int sumZ = 0; for( itepv = prodVector->begin(); itepv != prodVector->end(); itepv++ ){ G4ReactionProduct* product = (*itepv); G4cout << " PRODUCT GETMASS " << product->GetMass() << " > 4.? " << G4endl; if( product->GetMass() > 4. ) continue; // not ions bigger than alpha sumZ += product->GetDefinition()->GetAtomicNumber(); sumA += product->GetDefinition()->GetAtomicMass(); } G4int sumAZ = 10000 * sumZ + sumA; theAZs[sumAZ] ++; } std::map::const_iterator iteaz; for( iteaz = theAZs.begin(); iteaz != theAZs.end(); iteaz++ ) { G4double FSXS = HPChannel->GetFSCrossSection(ener, ip); // Macro XS // G4double cross = FSXS*(nAtomsPermm3*isotFraction)*barn; G4double cross = FSXS*(nAtomsPermm3); // FSXS has into account the isotope fraction // G4double cross = FSXS*(nAtomsPermm3)*barn; // FSXS has into account the isotope fraction G4double yield = (*iteaz).second / nSamples; theIonProds[(*iteaz).first] += yield * stepLength * cross / cubicVolume;// average ion production G4String FSType = HPChannel->theFSType; G4String channelName = SHUtils::GetChannelName(int(FS->theBaseZ), int(FS->theBaseA), FSType.substr(0,3), GmG4Utils::GetParticleShortName(aStep->GetTrack()->GetDefinition()->GetParticleName()) ); G4cout << material->GetName() << " " << hadProc->GetProcessName() << " CHANNEL= " << FSType << " " << channelName << " AZ= " << (*iteaz).first << " YIELD= " << yield << " ION_PRODUCTION= " << yield * stepLength * cross / cubicVolume << " XS micro " << FSXS << " XS macro " << cross << " Ener " << ener << " cross " << cross << " nAtomsPermm3 " << nAtomsPermm3 << G4endl; } } // energy loop // Get index: std::map::const_iterator iteip; const G4Isotope* isotope = elem->GetIsotope(ip); for( iteip = theIonProds.begin(); iteip != theIonProds.end(); iteip++ ){ G4int index = GetIndexAdjustingIon( (*iteip).first, isotope ); FillScorer( aStep, index, (*iteip).second, aStep->GetPreStepPoint()->GetWeight()); // call directly the filling with index } }// if bAdjust PHP } // if HasFSData } // numberOfChannels } // numberOfIsos } // numberOfElem /* }else if( modelName == "NeutronHPInelastic" ){ } else if( modelName == "NeutronHPCapture" ){ } else if( modelName == "ParticleHPnCapture" ){ */ } } } //------------------------------------------------------------------- G4int SHRadDecayScorerFromFlux::GetIndexFromProduct(G4ParticleHPProduct* product) { std::map::const_iterator ite = theHPProductIndices.find(product); if( ite == theHPProductIndices.end() ) { G4int massCode = G4int(product->GetMassCode()); G4Ions* ion = (G4Ions*)(theIonTable->GetIon(massCode/1000,massCode%1000,0.)); G4int index = theIonIndices.size(); theIonIndices[index] = ion; theHPProductIndices[product] = index; return index; } else { return (*ite).second; } } //------------------------------------------------------------------- G4int SHRadDecayScorerFromFlux::GetIndexAdjustingIon( G4int AZ, const G4Isotope* isotope ) { G4int newAZ = isotope->GetZ()*10000 + isotope->GetA() - AZ; if( newAZ / 10000 <= 0 || newAZ % 10000 <= 0 ) { G4Exception("SHRadDecayScorerFromFlux::GetIndexAdjustingIon", "Error", FatalException, ("WRONG AZ " +GmGenUtils::itoa(newAZ) + " THIS SHOULD NOT HAPPEN, CONTACT GAMOS AUTHORS").c_str()); } std::map::const_iterator ite = theAZindices.find(newAZ); if( ite == theAZindices.end() ) { G4Ions* ion = (G4Ions*)(theIonTable->GetIon(newAZ/10000, newAZ%10000,0.)); G4int index = theIonIndices.size(); theIonIndices[index] = ion; theAZindices[newAZ] = index; return index; } else { return (*ite).second; } } //-------------------------------------------------------------------- G4bool SHRadDecayScorerFromFlux::FillScorer(const G4Step*, G4int index, G4double val, G4double wei) { #ifndef GAMOS_NO_VERBOSE // if( SHVerb(debugVerb) ) G4cout << EvtMap << " EvtMap Number of entries " << EvtMap->entries() << " tmp " << theSumV_tmp.size() << " val " << val << " weight " << wei << G4endl; #endif if( bNewEvent ) { ScoreNewEvent(); } bNewEvent = false; G4double valwei = val*wei; theEvtMap->add(index,valwei); #ifndef GAMOS_NO_VERBOSE if( SHVerb(debugVerb) ) G4cout << "SHRadDecayScorerFromFlux" << " index " << index << " val " << val << " wei " << wei << G4endl; #endif theSumV_tmp[index] += val*wei; #ifndef GAMOS_NO_VERBOSE if( SHVerb(debugVerb) ) G4cout << "SHRadDecayScorerFromFlux chk theSumV_tmp " << index << " = " << theSumV_tmp[index] << G4endl; #endif sumALL += val*wei; #ifndef GAMOS_NO_VERBOSE if( SHVerb(debugVerb) ) G4cout << " EvtMap Number of entries " << theEvtMap->entries() << " tmp " << theSumV2.size() << " SUMALL " << sumALL << G4endl; #endif return true; } //-------------------------------------------------------------------- void SHRadDecayScorerFromFlux::ScoreNewEvent() { //--- Store X and X2 of last event // G4double eventSum = 0.; std::map::iterator ite; for( ite = theSumV_tmp.begin(); ite !=theSumV_tmp.end(); ite++ ){ G4int idx = (*ite).first; theSumV2[idx] += (*ite).second*(*ite).second; #ifndef GAMOS_NO_VERBOSE if( SHVerb(debugVerb) ) G4cout<< "SHRadDecayScorerFromFlux " << idx << " EVENT SumV= " << (*ite).second << " SumV2= " << theSumV2[idx] << G4endl; #endif /* if( theConvergenceTester ) { eventSum += (*ite).second; } */ } /* if( theConvergenceTester ) { theConvergenceTester->AddScore( eventSum ); } */ theSumV_tmp.clear(); *theRunMap += *theEvtMap; theEvtMap->clear(); } //-------------------------------------------------------------------- void SHRadDecayScorerFromFlux::CalculateErrors(G4THitsMap* RunMap) { G4double nev; nev = GmNumberOfEvent::GetNumberOfEvent(); std::map::iterator ite; for(ite = RunMap->GetMap()->begin(); ite != RunMap->GetMap()->end(); ite++){ G4double sumWX = (*(ite->second)); theError[ite->first] = GetError( ite->first, sumWX, nev ); #ifndef GAMOS_NO_VERBOSE if( SHVerb(debugVerb) ) G4cout << ite->first << " CERROR " << theError[ ite->first ] << G4endl; #endif } } //-------------------------------------------------------------------- G4double SHRadDecayScorerFromFlux::GetError( G4int index ) { return theError[index]; } //-------------------------------------------------------------------- G4double SHRadDecayScorerFromFlux::GetError( G4int index, G4double sumWX, G4double nEvents ) { if( theSumV2.find( index ) == theSumV2.end() ) return 0.; // G4double error = (theSumV2[index]*nEvents - sumWX*sumWX) / (nEvents*nEvents*(nEvents-1)); G4double error = (theSumV2[index]*nEvents - sumWX*sumWX) / (nEvents-1); #ifndef GAMOS_NO_VERBOSE if( SHVerb(debugVerb) ) G4cout << " GetError " << index << " e= " << error/(nEvents*nEvents) << " S1 " << theSumV2[index]*nEvents << " S2 " << sumWX*sumWX << " S3 " << (nEvents*nEvents*(nEvents-1)) << G4endl; #endif if( error <= 0. ) { if( error < -1.E-30 ) G4cerr << " !!WARNING Error squared in scorer SHRadDecayScorerFromFlux is negative: " << error << G4endl; error = 0.; } else { error = std::sqrt(error)/nEvents; } #ifndef GAMOS_NO_VERBOSE G4cout<< SHVerb(debugVerb) << " GetError " << index << " e= " << error << " e2= " << error*error << " nev " << nEvents << " sumWX " << sumWX << " sumW2 " << theSumV2[index] << G4endl; #endif // G4cout << " SHRadDecayScorerFromFlux::GetError " << index << " = " << error << " sumW " << sumWX << " sumW2 " << theSumV2[index] << " nEvents " << nEvents << G4endl; return error; } //-------------------------------------------------------------------- G4double SHRadDecayScorerFromFlux::GetErrorRelative( G4int index, G4double sumWX, G4double nEvents ) { G4double errorrel = GetError( index ); #ifndef GAMOS_NO_VERBOSE if( SHVerb(debugVerb) ) G4cout << " GetErrorRelative temp er " << errorrel << G4endl; #endif // divide by averageX if( sumWX != 0. ) errorrel /= (sumWX/nEvents); #ifndef GAMOS_NO_VERBOSE if( SHVerb(debugVerb) ) G4cout << " GetErrorRelative " << index << " er= " << errorrel << " nev " << nEvents << " sumWX " << sumWX << G4endl; #endif return errorrel; } G4Ions* SHRadDecayScorerFromFlux::GetIon(G4int index) const { std::map::const_iterator ite = theIonIndices.find(index); if( ite == theIonIndices.end() ) { G4Exception("SHRadDecayScorerFromFlux::GetIon", "Wron argument", FatalErrorInArgument, ("index not found " + GmGenUtils::itoa(index)).c_str()); } return (*ite).second; }