G4LossTableBuilder.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id$
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:     G4LossTableBuilder
00034 //
00035 // Author:        Vladimir Ivanchenko
00036 //
00037 // Creation date: 03.01.2002
00038 //
00039 // Modifications:
00040 //
00041 // 23-01-03 V.Ivanchenko Cut per region
00042 // 21-07-04 V.Ivanchenko Fix problem of range for dedx=0
00043 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko)
00044 // 07-12-04 Fix of BuildDEDX table (V.Ivanchenko)
00045 // 27-03-06 Add bool options isIonisation (V.Ivanchenko)
00046 // 16-01-07 Fill new (not old) DEDX table (V.Ivanchenko)
00047 // 12-02-07 Use G4LPhysicsFreeVector for the inverse range table (V.Ivanchenko)
00048 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
00049 //
00050 // Class Description:
00051 //
00052 // -------------------------------------------------------------------
00053 //
00054 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00056 
00057 #include "G4LossTableBuilder.hh"
00058 #include "G4SystemOfUnits.hh"
00059 #include "G4PhysicsTable.hh"
00060 #include "G4PhysicsLogVector.hh"
00061 #include "G4PhysicsTableHelper.hh"
00062 #include "G4LPhysicsFreeVector.hh"
00063 #include "G4ProductionCutsTable.hh"
00064 #include "G4MaterialCutsCouple.hh"
00065 #include "G4Material.hh"
00066 #include "G4VEmModel.hh"
00067 #include "G4ParticleDefinition.hh"
00068 #include "G4LossTableManager.hh"
00069 
00070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00071 
00072 G4LossTableBuilder::G4LossTableBuilder() 
00073 {
00074   splineFlag = true;
00075   isInitialized = false;
00076 
00077   theDensityFactor = new std::vector<G4double>;
00078   theDensityIdx = new std::vector<G4int>;
00079   theFlag = new std::vector<G4bool>;
00080 }
00081 
00082 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00083 
00084 G4LossTableBuilder::~G4LossTableBuilder() 
00085 {
00086   delete theDensityFactor;
00087   delete theDensityIdx;
00088   delete theFlag;
00089 }
00090 
00091 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00092 
00093 void 
00094 G4LossTableBuilder::BuildDEDXTable(G4PhysicsTable* dedxTable,
00095                                    const std::vector<G4PhysicsTable*>& list)
00096 {
00097   size_t n_processes = list.size();
00098   //G4cout << "Nproc= " << n_processes << " Ncoup= " << dedxTable->size() << G4endl;
00099   if(1 >= n_processes) { return; }
00100 
00101   size_t nCouples = dedxTable->size();
00102   if(0 >= nCouples) { return; }
00103 
00104   for (size_t i=0; i<nCouples; ++i) {
00105     //    if ((*theFlag)[i]) {
00106     G4PhysicsLogVector* pv0 = static_cast<G4PhysicsLogVector*>((*(list[0]))[i]);
00107     if(pv0) {
00108       size_t npoints = pv0->GetVectorLength();
00109       G4PhysicsLogVector* pv = new G4PhysicsLogVector(*pv0);
00110       //    pv = new G4PhysicsLogVector(elow, ehigh, npoints-1);
00111       pv->SetSpline(splineFlag);
00112       for (size_t j=0; j<npoints; ++j) {
00113         G4double dedx = 0.0;
00114         for (size_t k=0; k<n_processes; ++k) {
00115           G4PhysicsVector* pv1   = (*(list[k]))[i];
00116           dedx += (*pv1)[j];
00117         }
00118         pv->PutValue(j, dedx);
00119       }
00120       if(splineFlag) { pv->FillSecondDerivatives(); }
00121       G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv);
00122     }
00123   }
00124 }
00125 
00126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00127 
00128 void G4LossTableBuilder::BuildRangeTable(const G4PhysicsTable* dedxTable,
00129                                          G4PhysicsTable* rangeTable,
00130                                          G4bool isIonisation)
00131 // Build range table from the energy loss table
00132 {
00133   size_t nCouples = dedxTable->size();
00134   if(0 >= nCouples) { return; }
00135 
00136   size_t n = 100;
00137   G4double del = 1.0/(G4double)n;
00138 
00139   for (size_t i=0; i<nCouples; ++i) {
00140     if(isIonisation) {
00141       if( !(*theFlag)[i] ) { continue; }
00142     }
00143     G4PhysicsLogVector* pv = static_cast<G4PhysicsLogVector*>((*dedxTable)[i]);
00144     size_t npoints = pv->GetVectorLength();
00145     size_t bin0    = 0;
00146     G4double elow  = pv->Energy(0);
00147     G4double ehigh = pv->Energy(npoints-1);
00148     G4double dedx1 = (*pv)[0];
00149 
00150     //G4cout << "i= " << i << "npoints= " << npoints << " dedx1= " << dedx1 << G4endl;
00151 
00152     // protection for specific cases dedx=0
00153     if(dedx1 == 0.0) {
00154       for (size_t k=1; k<npoints; ++k) {
00155         bin0++;
00156         elow  = pv->Energy(k);
00157         dedx1 = (*pv)[k];
00158         if(dedx1 > 0.0) { break; }
00159       }
00160       npoints -= bin0;
00161     }
00162     //G4cout<<"New Range vector" << G4endl;
00163     //G4cout<<"nbins= "<<npoints-1<<" elow= "<<elow<<" ehigh= "<<ehigh
00164     //      <<" bin0= " << bin0 <<G4endl;
00165 
00166     // initialisation of a new vector
00167     if(npoints < 2) { npoints = 2; }
00168 
00169     delete (*rangeTable)[i];
00170     G4PhysicsLogVector* v;
00171     if(0 == bin0) { v = new G4PhysicsLogVector(*pv); }
00172     else { v = new G4PhysicsLogVector(elow, ehigh, npoints-1); }
00173 
00174     // dedx is exact zero cannot build range table
00175     if(2 == npoints) {
00176       v->PutValue(0,1000.);
00177       v->PutValue(1,2000.);
00178       G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
00179       return;
00180     }
00181     v->SetSpline(splineFlag);
00182 
00183     // assumed dedx proportional to beta
00184     G4double energy1 = v->Energy(0);
00185     G4double range   = 2.*energy1/dedx1;
00186     //G4cout << "range0= " << range << G4endl;
00187     v->PutValue(0,range);
00188 
00189     for (size_t j=1; j<npoints; ++j) {
00190 
00191       G4double energy2 = v->Energy(j);
00192       G4double de      = (energy2 - energy1) * del;
00193       G4double energy  = energy2 + de*0.5;
00194       G4double sum = 0.0;
00195       //G4cout << "j= " << j << " e1= " << energy1 << " e2= " << energy2 
00196       //       << " n= " << n << G4endl;
00197       for (size_t k=0; k<n; ++k) {
00198         energy -= de;
00199         dedx1 = pv->Value(energy);
00200         if(dedx1 > 0.0) { sum += de/dedx1; }
00201       }
00202       range += sum;
00203       v->PutValue(j,range);
00204       energy1 = energy2;
00205     }
00206     if(splineFlag) { v->FillSecondDerivatives(); }
00207     G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
00208   }
00209 }
00210 
00211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00212 
00213 void G4LossTableBuilder::BuildInverseRangeTable(const G4PhysicsTable* rangeTable,
00214                                                 G4PhysicsTable* invRangeTable,
00215                                                 G4bool isIonisation)
00216 // Build inverse range table from the energy loss table
00217 {
00218   size_t nCouples = rangeTable->size();
00219   if(0 >= nCouples) { return; }
00220 
00221   for (size_t i=0; i<nCouples; ++i) {
00222 
00223     if(isIonisation) {
00224       if( !(*theFlag)[i] ) { continue; }
00225     }
00226     G4PhysicsVector* pv = (*rangeTable)[i];
00227     size_t npoints = pv->GetVectorLength();
00228     G4double rlow  = (*pv)[0];
00229     G4double rhigh = (*pv)[npoints-1];
00230       
00231     delete (*invRangeTable)[i];
00232     G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(npoints,rlow,rhigh);
00233     v->SetSpline(splineFlag);
00234 
00235     for (size_t j=0; j<npoints; ++j) {
00236       G4double e  = pv->Energy(j);
00237       G4double r  = (*pv)[j];
00238       v->PutValues(j,r,e);
00239     }
00240     if(splineFlag) { v->FillSecondDerivatives(); }
00241 
00242     G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v);
00243   }
00244 }
00245 
00246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00247 
00248 void 
00249 G4LossTableBuilder::InitialiseBaseMaterials(G4PhysicsTable* table)
00250 {
00251   size_t nCouples = table->size();
00252   size_t nFlags = theFlag->size();
00253 
00254   if(nCouples == nFlags && isInitialized) { return; }
00255 
00256   isInitialized = true;
00257 
00258   //G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials Ncouples= " 
00259   //     << nCouples << "  FlagSize= " << nFlags << G4endl; 
00260 
00261   // variable density check
00262   const G4ProductionCutsTable* theCoupleTable=
00263     G4ProductionCutsTable::GetProductionCutsTable();
00264 
00265   /*
00266   for(size_t i=0; i<nFlags; ++i) {
00267     G4cout << "CoupleIdx= " << i << "  Flag= " <<  (*theFlag)[i] 
00268            << " tableFlag= " << table->GetFlag(i) << "  "
00269            << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
00270            << G4endl;
00271   }
00272   */
00273 
00274   // expand vectors
00275   if(nFlags < nCouples) {
00276     for(size_t i=nFlags; i<nCouples; ++i) { theDensityFactor->push_back(1.0); }
00277     for(size_t i=nFlags; i<nCouples; ++i) { theDensityIdx->push_back(-1); }
00278     for(size_t i=nFlags; i<nCouples; ++i) { theFlag->push_back(true); }
00279   }
00280   for(size_t i=0; i<nCouples; ++i) {
00281 
00282     // base material is needed only for a couple which is not
00283     // initialised and for which tables will be computed
00284     (*theFlag)[i] = table->GetFlag(i);
00285     if ((*theDensityIdx)[i] < 0) {
00286       (*theDensityIdx)[i] = i;
00287       const G4MaterialCutsCouple* couple =
00288         theCoupleTable->GetMaterialCutsCouple(i);
00289       const G4ProductionCuts* pcuts = couple->GetProductionCuts();
00290       const G4Material* mat  = couple->GetMaterial();
00291       const G4Material* bmat = mat->GetBaseMaterial();
00292 
00293       // base material exists - find it and check if it can be reused
00294       if(bmat) {
00295         for(size_t j=0; j<nCouples; ++j) {
00296 
00297           if(j == i) { continue; }
00298           const G4MaterialCutsCouple* bcouple =
00299               theCoupleTable->GetMaterialCutsCouple(j);
00300 
00301           if(bcouple->GetMaterial() == bmat && 
00302              bcouple->GetProductionCuts() == pcuts) {
00303 
00304             // based couple exist in the same region
00305             (*theDensityIdx)[i] = j;
00306             (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity();
00307             (*theFlag)[i] = false;
00308 
00309             // ensure that there will no double initialisation
00310             (*theDensityIdx)[j] = j;
00311             (*theDensityFactor)[j] = 1.0;
00312             (*theFlag)[j] = true;
00313             break;
00314           }
00315         }
00316       }
00317     }
00318   }
00319   /*
00320   for(size_t i=0; i<nCouples; ++i) {
00321     G4cout << "CoupleIdx= " << i << "  Flag= " <<  (*theFlag)[i] 
00322            << "  TableFlag= " << table->GetFlag(i) << "  "
00323            << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
00324            << G4endl;
00325   }
00326   G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials end" 
00327          << G4endl; 
00328   */
00329 }
00330 
00331 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00332 
00333 void G4LossTableBuilder::InitialiseCouples()
00334 {
00335   isInitialized = true;
00336 
00337   // variable density initialisation for the cas without tables 
00338   const G4ProductionCutsTable* theCoupleTable=
00339     G4ProductionCutsTable::GetProductionCutsTable();
00340   size_t nCouples = theCoupleTable->GetTableSize();
00341   //G4cout << "%%%%%% G4LossTableBuilder::InitialiseCouples() nCouples= " 
00342   //     << nCouples << G4endl; 
00343 
00344   theDensityFactor->resize(nCouples, 1.0);
00345   theDensityIdx->resize(nCouples, -1);
00346   theFlag->resize(nCouples, true);
00347 
00348   for(size_t i=0; i<nCouples; ++i) {
00349 
00350     // base material is needed only for a couple which is not
00351     // initialised and for which tables will be computed
00352     if ((*theDensityIdx)[i] < 0) {
00353       (*theDensityIdx)[i] = i;
00354       const G4MaterialCutsCouple* couple =
00355         theCoupleTable->GetMaterialCutsCouple(i);
00356       const G4ProductionCuts* pcuts = couple->GetProductionCuts();
00357       const G4Material* mat  = couple->GetMaterial();
00358       const G4Material* bmat = mat->GetBaseMaterial();
00359 
00360       // base material exists - find it and check if it can be reused
00361       if(bmat) {
00362         for(size_t j=0; j<nCouples; ++j) {
00363 
00364           if(j == i) { continue; }
00365           const G4MaterialCutsCouple* bcouple =
00366             theCoupleTable->GetMaterialCutsCouple(j);
00367 
00368           if(bcouple->GetMaterial() == bmat && 
00369              bcouple->GetProductionCuts() == pcuts) {
00370 
00371             // based couple exist in the same region
00372             (*theDensityIdx)[i] = j;
00373             (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity();
00374             (*theFlag)[i] = false;
00375 
00376             // ensure that there will no double initialisation
00377             (*theDensityIdx)[j] = j;
00378             (*theDensityFactor)[j] = 1.0;
00379             (*theFlag)[j] = true;
00380             break;
00381           }
00382         }
00383       }
00384     }
00385   }
00386   /*  
00387   for(size_t i=0; i<nCouples; ++i) {
00388     G4cout << "CoupleIdx= " << i << "  Flag= " <<  (*theFlag)[i] << "  "
00389            << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName()
00390            << "  DensityFactor= " << (*theDensityFactor)[i]
00391            << G4endl;
00392   }
00393   G4cout << "%%%%%% G4LossTableBuilder::InitialiseCouples() end" << G4endl; 
00394   */
00395 }
00396 
00397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00398 
00399 G4PhysicsTable* 
00400 G4LossTableBuilder::BuildTableForModel(G4PhysicsTable* aTable, 
00401                                        G4VEmModel* model, 
00402                                        const G4ParticleDefinition* part,
00403                                        G4double emin, G4double emax,
00404                                        G4bool spline)
00405 {
00406   // check input
00407   G4PhysicsTable* table = G4PhysicsTableHelper::PreparePhysicsTable(aTable);
00408   if(!table) { return table; }
00409   if(emin >= emax) { 
00410     table->clearAndDestroy();
00411     delete table;
00412     table = 0;
00413     return table; 
00414   }
00415   InitialiseBaseMaterials(table);
00416 
00417   G4int nbins = G4int(std::log10(emax/emin) + 0.5)
00418     *G4LossTableManager::Instance()->GetNumberOfBinsPerDecade();
00419   if(nbins < 3) { nbins = 3; }
00420 
00421   // Access to materials
00422   const G4ProductionCutsTable* theCoupleTable=
00423         G4ProductionCutsTable::GetProductionCutsTable();
00424   size_t numOfCouples = theCoupleTable->GetTableSize();
00425 
00426   G4PhysicsLogVector* aVector = 0;
00427   G4PhysicsLogVector* bVector = 0;
00428 
00429   for(size_t i=0; i<numOfCouples; ++i) {
00430 
00431     //G4cout<< "i= " << i << " Flag=  " << GetFlag(i) << G4endl;
00432 
00433     if (GetFlag(i)) {
00434 
00435       // create physics vector and fill it
00436       const G4MaterialCutsCouple* couple = 
00437         theCoupleTable->GetMaterialCutsCouple(i);
00438       delete (*table)[i];
00439 
00440       // if start from zero then change the scale
00441 
00442       const G4Material* mat = couple->GetMaterial();
00443 
00444       G4double tmin = std::max(emin,model->MinPrimaryEnergy(mat,part));
00445       if(0.0 >= tmin) { tmin = eV; }
00446       G4int n = nbins + 1;
00447 
00448       if(tmin >= emax) {
00449         aVector = 0;
00450       } else if(tmin > emin) {
00451         G4int bin = nbins*G4int(std::log10(emax/tmin) + 0.5);
00452         if(bin < 3) { bin = 3; }
00453         n = bin + 1;
00454         aVector = new G4PhysicsLogVector(tmin, emax, bin);
00455 
00456       } else if(!bVector) {
00457         aVector = new G4PhysicsLogVector(emin, emax, nbins);
00458         bVector = aVector;
00459 
00460       } else {
00461         aVector = new G4PhysicsLogVector(*bVector);
00462       }
00463 
00464       if(aVector) {
00465         aVector->SetSpline(spline);
00466         for(G4int j=0; j<n; ++j) {
00467           aVector->PutValue(j, model->Value(couple, part, aVector->Energy(j)));
00468         }
00469         if(spline) { aVector->FillSecondDerivatives(); }
00470       }
00471       G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
00472     }
00473   }
00474   /*
00475   G4cout << "G4LossTableBuilder::BuildTableForModel done for "
00476          << part->GetParticleName() << " and "<< model->GetName()
00477          << "  " << table << G4endl;
00478   */
00479   return table; 
00480 }
00481 
00482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
00483 

Generated on Mon May 27 17:48:49 2013 for Geant4 by  doxygen 1.4.7