Ignore:
Timestamp:
Jan 8, 2010, 3:02:48 PM (14 years ago)
Author:
garnier
Message:

update to geant4.9.3

Location:
trunk/examples/advanced/xray_fluorescence/src
Files:
10 edited

Legend:

Unmodified
Added
Removed

  • trunk/examples/advanced/xray_fluorescence/src/XrayFluoAnalysisManager.cc

    r807 r1230  
    3232// History:
    3333// -----------
     34// 20 Jul 2007 A.Mantero, code cleaning and update. Replaced histos with tuple
    3435// 11 Jul 2003 A.Mantero, code cleaning / Plotter-XML addiction
    3536//    Sep 2002 A.Mantero, AIDA3.0 Migration
    … …  
    4647#include "XrayFluoAnalysisManager.hh"
    4748#include "G4Step.hh"
     49#include "XrayFluoDetectorConstruction.hh"
     50#include "G4VPhysicalVolume.hh"
    4851
    4952XrayFluoAnalysisManager* XrayFluoAnalysisManager::instance = 0;
    … …  
    5356XrayFluoAnalysisManager::XrayFluoAnalysisManager()
    5457  :outputFileName("xrayfluo"), visPlotter(false), phaseSpaceFlag(false), physicFlag (false), persistencyType("xml"),
    55    deletePersistencyFile(true), gunParticleEnergies(0), gunParticleTypes(0), analysisFactory(0), tree(0),histogramFactory(0), plotter(0)
     58   deletePersistencyFile(true), gunParticleEnergies(0), gunParticleTypes(0), analysisFactory(0), tree(0), treeDet(0), histogramFactory(0), plotter(0)
    5659{
    5760  //creating the messenger
    … …  
    197200    // Book tuple
    198201    std::vector<std::string> columnNames;
    199     columnNames.push_back("Particle");
    200     columnNames.push_back("Energies");
     202    columnNames.push_back("particle");
     203    columnNames.push_back("energies");
    201204    columnNames.push_back("momentumTheta");
    202205    columnNames.push_back("momentumPhi");
    203     columnNames.push_back("Processes");
     206    columnNames.push_back("processes");
     207    columnNames.push_back("material");
     208    columnNames.push_back("origin");
     209    columnNames.push_back("depth");
    204210
    205211    std::vector<std::string> columnTypes;
    … …  
    209215    columnTypes.push_back("double");
    210216    columnTypes.push_back("int"); // useful for hbk
    211 
     217    columnTypes.push_back("int"); // useful for hbk
     218    columnTypes.push_back("int");
     219    columnTypes.push_back("double");
    212220
    213221    tupleFluo = tupleFactory->create("10", "Total Tuple", columnNames, columnTypes, "");
    … …  
    262270  gunParticleTypes = new std::vector<G4String>;
    263271
     272  // dal punto di vista della programmazione e' un po' una porcata.....
     273
    264274  AIDA::ITreeFactory* treeDataFactory = analysisFactory->createTreeFactory();
    265   AIDA::ITree* treeData = treeDataFactory->create(fileName,persistencyType,true,false); // input file
     275  AIDA::ITree* treeData = treeDataFactory->create(fileName,persistencyType,false, false); // input file
    266276  AIDA::IManagedObject* mo = treeData->find("10");
     277  assert(mo);
     278  G4cout <<"mo Name:" << mo->name()<< G4endl;
    267279  AIDA::ITuple* tupleData = dynamic_cast<AIDA::ITuple*>(mo);
     280
     281  //AIDA::ITuple* tupleData = (AIDA::ITuple*) treeData->find("10");
     282
     283  assert(tupleData);
     284  // da riscrivere un pochino melgio. Del resto serve per pochissmo, poi non serve piu'
     285
     286
    268287  tupleData->start();
     288  G4int processesIndex = tupleData->findColumn("processes");
     289  G4int energyIndex = tupleData->findColumn("energies");
     290  G4int particleIndex = tupleData->findColumn("particle");
    269291
    270292  while (tupleData->next()) {
    271     if (raileighFlag ^ (!raileighFlag && (tupleData->getInt(4)) ) ) {
    272       gunParticleEnergies->push_back(tupleData->getDouble(1));
    273       if (tupleData->getInt(0) == 1 ) gunParticleTypes->push_back("gamma");
    274       if (tupleData->getInt(0) == 0 ) gunParticleTypes->push_back("e-");
     293    if (raileighFlag ^ (!raileighFlag && (tupleData->getInt(processesIndex) == 1 ||
     294                                          tupleData->getInt(processesIndex) == 2)) )  {
     295      gunParticleEnergies->push_back(tupleData->getDouble(energyIndex)*MeV);
     296      if (tupleData->getInt(particleIndex) == 1 ) gunParticleTypes->push_back("gamma");
     297      if (tupleData->getInt(particleIndex) == 0 ) gunParticleTypes->push_back("e-");
    275298    }
    276299
    … …  
    344367  if(plotter) {
    345368    if (phaseSpaceFlag){
     369
     370
     371      // altra porcata di programmazione e cmq mi serve un istogramma. Da rivedere il tutto.
    346372      // Plotting the spectrum of Gamma exiting the sample - cloud1
    347373      AIDA::ICloud1D& cloud = *cloud_1;
    … …  
    349375                                              " Particle == std::string(\"gamma\") ");
    350376      AIDA::IEvaluator* evaluatorEnergy = tupleFactory->createEvaluator("Energies");
     377
     378      // di nuovo, del resto serve solo per un attimo.. da correggere cmq.
     379
    351380      filterGamma->initialize(*tupleFluo);
    352381      evaluatorEnergy->initialize(*tupleFluo);
    … …  
    396425    G4ThreeVector momentum=0;
    397426    G4double particleEnergy=0;
    398    
    399     // Select volume from wich the step starts
     427    G4String sampleMaterial="";
     428    G4double particleDepth=0;
     429    G4int isBornInTheSample=0;
     430    XrayFluoDetectorConstruction* detector = XrayFluoDetectorConstruction::GetInstance();   
     431
     432    // Select volume from which the step starts
    400433    if(aStep->GetPreStepPoint()->GetPhysicalVolume()->GetName()=="Sample"){
    401       // Select volume in wich the step ends
    402       if (physicFlag ^ (!physicFlag && (aStep->GetTrack()->GetNextVolume()->GetName() == "World" ))) {
     434      G4ThreeVector creationPos = aStep->GetTrack()->GetVertexPosition();
     435      // qua serve ancora una selezione: deve essere interno al sample
     436      //codice "a mano" per controllare il volume di orgine della traccia
     437      /*
     438      G4Double sampleXplus  = detector->GetSampleXplusFaceCoord();
     439      G4Double sampleXminus = detector->GetSampleXminusFaceCoord();
     440      G4Double sampleYplus  = detector->GetSampleYplusFaceCoord();
     441      G4Double sampleYminus = detector->GetSampleYminusFaceCoord();
     442      G4Double sampleZplus  = detector->GetSampleYplusFaceCoord();
     443      G4Double sampleZminus = detector->GetSampleYminusFaceCoord();
     444
     445      G4bool ZsampleCheck = (creationPos.z()<=sampleZminus) && (creationPos.z() >= sampleZplus);
     446      G4bool XsampleCheck = (creationPos.x()<=sampleZminus) && (creationPos.x() >= sampleXplus);
     447      G4bool YsampleCheck = (creationPos.y()<=sampleZminus) && (creationPos.y() >= sampleYplus);
     448      */
     449      G4VPhysicalVolume* creationPosVolume = detector->GetGeometryNavigator()->LocateGlobalPointAndSetup(creationPos);
     450
     451      // if physicflag is on, I record all the data of all the particle born in the sample.
     452      // If it is off (but phase space production is on) I record data
     453      // only for particles creted inside the sample and exiting it. 
     454      if (physicFlag ^ (!physicFlag &&
     455                        (aStep->GetPostStepPoint()->GetStepStatus() == fGeomBoundary)
     456                        //ZsampleCheck && XsampleCheck && YsampleCheck)) {
     457                        ) ) {
    403458        //
    404459        //
    … …  
    410465        momentum = aStep->GetTrack()->GetDynamicParticle()->GetMomentum();
    411466        particleEnergy = aStep->GetPreStepPoint()->GetKineticEnergy();
     467        if (creationPosVolume->GetName() == "Sample" ) {
     468          isBornInTheSample = 1;
     469          particleDepth = creationPosVolume->GetLogicalVolume()->GetSolid()->DistanceToOut(creationPos, G4ThreeVector(0,0,-1));
     470        }
     471        else {
     472          particleDepth = -1;
     473        }
     474        // metodo per ottenere la profondita' "a mano"
     475        //      particleDepth = std::abs(creationPos.z() - detector->GetSampleIlluminatedFaceCoord());
     476
    412477        G4int parent;
    413478        if(aStep->GetTrack()->GetCreatorProcess()){
    414479          parentProcess = aStep->GetTrack()->GetCreatorProcess()->GetProcessName();
    415           parent = 1;
     480          parent = 5;
     481          // hack for HBK
     482          if (parentProcess == "") parent = 0;
     483          if (parentProcess == "IONI") parent = 1;
     484          if (parentProcess == "LowEnPhotoElec") parent = 2;
     485          if (parentProcess == "Transportation") parent = 3;// should never happen
     486          if (parentProcess == "initStep") parent = 4;// should never happen
    416487        }       
    417488        else {
    418           parentProcess = "Not Known";
    419           parent = 0;
     489          parentProcess = "Not Known -- (primary generator + Rayleigh)";
     490          parent = 5;
    420491        }
     492        G4int sampleMat=0;
     493        if(aStep->GetTrack()){
     494          sampleMaterial = aStep->GetTrack()->GetMaterial()->GetName();
     495          if (sampleMaterial == ("Dolorite" || "Anorthosite" || "Mars1" || "IceBasalt" || "HPGe")) sampleMat=1;
     496            }
    421497        // filling tuple
    422498       
    423499        //      G4cout<< particleType << G4endl;
    424         G4int part = 2 ;
     500        G4int part = -1 ;
    425501        if (particleType == "gamma") part =1;
    426502        if (particleType == "e-") part = 0;
    427        
     503        if (particleType == "proton") part = 2;
     504       
     505
    428506        tupleFluo->fill(0,part);
    429507        tupleFluo->fill(1,particleEnergy);
    430508        tupleFluo->fill(2,momentum.theta());
    431509        tupleFluo->fill(3,momentum.phi());
    432         tupleFluo->fill(4,parent); //hacked to be useful for hbk
     510        tupleFluo->fill(4,parent); //hacked to be used with hbk
     511        tupleFluo->fill(5,sampleMat); //hacked to be used with hbk
     512        tupleFluo->fill(6,isBornInTheSample);
     513        tupleFluo->fill(7,particleDepth);
    433514       
    434515        tupleFluo->addRow();
    … …  
    678759   
    679760   
    680     AIDA::IFilter* filterAngle = tupleFactory->createFilter("(momentumPhi   >= (220. * (3.1415926/180.) )) &&
    681                                                              (momentumPhi   <= (230. * (3.1415926/180.) )) &&
    682                                                              (momentumTheta >= (130. * (3.1415926/180.) )) &&
    683                                                              (momentumTheta <= (140. * (3.1415926/180.) )) " );
     761    AIDA::IFilter* filterAngle = tupleFactory->createFilter(
     762
     763"(momentumPhi   >= (220. * (3.1415926/180.) )) && (momentumPhi   <= (230. * (3.1415926/180.) )) && (momentumTheta >= (130. * (3.1415926/180.) )) && (momentumTheta <= (140. * (3.1415926/180.) ))" );
    684764   
    685765   
    … …  
    743823
    744824#endif
    745 
    746 
    747 
    748 
    749 
    750 
    751 
    752 

  • trunk/examples/advanced/xray_fluorescence/src/XrayFluoDataSet.cc

    r807 r1230  
    131131
    132132
    133 void XrayFluoDataSet::LoadData(const G4String& fileName)
    134 {
     133G4bool XrayFluoDataSet::LoadData(const G4String& fileName) {
     134
    135135  // Build the complete string identifying the file with the data set
    136136 
    … …  
    199199 
    200200  file.close();
     201  return true;
    201202}
    202203void XrayFluoDataSet::PrintData() const
    … …  
    217218
    218219
    219 
    220 
    221 
    222 
    223 
    224 
    225 
    226 
    227 
    228 
    229 
    230 

  • trunk/examples/advanced/xray_fluorescence/src/XrayFluoDetectorConstruction.cc

    r807 r1230  
    5959
    6060
    61 #include "G4Region.hh"
    62 #include "G4RegionStore.hh"
     61//#include "G4Region.hh"
     62//#include "G4RegionStore.hh"
    6363
    6464//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    6666
    6767XrayFluoDetectorConstruction::XrayFluoDetectorConstruction()
    68   : detectorType(0),sampleGranularity(false), phaseSpaceFlag(false),
     68  : aNavigator(0), detectorType(0),sampleGranularity(false), phaseSpaceFlag(false),
    6969    DeviceSizeX(0), DeviceSizeY(0),DeviceThickness(0),
    7070    solidWorld(0),logicWorld(0),physiWorld(0),
    … …  
    8383{
    8484  materials = XrayFluoNistMaterials::GetInstance();
     85
     86  aNavigator = new G4Navigator();
    8587 
    8688  DefineDefaultMaterials();
    … …  
    8991  NbOfPixelColumns  =  1; // should be 1
    9092  NbOfPixels        =  NbOfPixelRows*NbOfPixelColumns;
    91   PixelSizeXY       =  std::sqrt(40.) * mm; // should be std::sqrt(40) * mm
     93  PixelSizeXY       =  7 * cm;// should be std::sqrt(40) * mm
    9294  PixelThickness = 3.5 * mm; //should be 3.5 mm
    9395
    … …  
    9597  G4cout << "PixelSizeXY(cm): "<< PixelSizeXY/cm << G4endl;
    9698
    97   ContactSizeXY     = std::sqrt(40.) * mm; //std::sqrt(40) * mm; //should be the same as PixelSizeXY
     99  ContactSizeXY     = PixelSizeXY; //std::sqrt(40) * mm; //should be the same as PixelSizeXY
    98100  SampleThickness = 4 * mm;
    99101  SampleSizeXY = 3. * cm;
    … …  
    132134  ComputeApparateParameters();
    133135
    134   G4String regName = "SampleRegion";
    135   sampleRegion = new G4Region(regName); 
     136  // G4String regName = "SampleRegion";
     137  //sampleRegion = new G4Region(regName); 
    136138
    137139  if (!phaseSpaceFlag) SetDetectorType(defaultDetectorType);
    … …  
    205207  //define materials of the apparate
    206208
    207   sampleMaterial = materials->GetMaterial("Mars1");
     209  sampleMaterial = materials->GetMaterial("Dolorite");
    208210  Dia1Material = materials->GetMaterial("G4_Pb");
    209211  Dia3Material = materials->GetMaterial("Galactic");
    210   pixelMaterial = materials->GetMaterial("G4_Si");
    211   OhmicPosMaterial = materials->GetMaterial("G4_Cu");
     212  pixelMaterial = materials->GetMaterial("SiLi");
     213  //pixelMaterial = materials->GetMaterial(detectorType->GetDetectorMaterial());
     214  OhmicPosMaterial = materials->GetMaterial("Cu");
    212215  OhmicNegMaterial = materials->GetMaterial("G4_Pb");
    213216  defaultMaterial = materials->GetMaterial("Galactic");
    … …  
    246249  //ComputeApparateParameters();
    247250 
    248   //world
     251  //world and associated navigator
    249252 
    250253  solidWorld = new G4Box("World",                               //its name
    … …  
    261264                                 false,                 //no boolean operation
    262265                                 0);                    //copy number
     266
     267  aNavigator->SetWorldVolume(physiWorld);
     268
    263269 
    264270  //HPGeDetector
    … …  
    378384          PixelCopyNb += PixelCopyNb;
    379385          G4cout << "PixelCopyNb: " << PixelCopyNb << G4endl;
    380       }
     386        }
    381387     
    382388      }
    … …  
    611617  // cut per region
    612618 
    613   logicSample->SetRegion(sampleRegion);
    614   sampleRegion->AddRootLogicalVolume(logicSample);
     619  //logicSample->SetRegion(sampleRegion);
     620  //sampleRegion->AddRootLogicalVolume(logicSample);
    615621 
    616622
    … …  
    758764
    759765
    760     G4cout << "Material Change in Progress" << newMaterial << G4endl;
     766    G4cout << "Material Change in Progress " << newMaterial << G4endl;
    761767    sampleMaterial = materials->GetMaterial(newMaterial);
    762768    logicSample->SetMaterial(sampleMaterial);

  • trunk/examples/advanced/xray_fluorescence/src/XrayFluoDetectorMessenger.cc

    r807 r1230  
    6060
    6161  sampleCmd = new G4UIcmdWithAString("/apparate/sampleMaterial",this);
    62   sampleCmd->SetGuidance("select a diferent material for the sample: materials can be: Dolorite, Anorthosite, Mars1, IceBasalt, HPGe OR choosen from Nist database (see /material/nist/listMaterials for details");
     62  sampleCmd->SetGuidance("select a diferent material for the sample: materials can be: Dolorite, Anorthosite, Mars1, HawaiianWD, HawaiianRF, IceBasalt, IcelandicWD, IcelandicRF, Gabbro, GabbroWD, GabbroRF, HPGe OR choosen from Nist database (see /material/nist/listMaterials for details");
    6363  sampleCmd->SetParameterName("material",true);
    6464  sampleCmd->SetDefaultValue("mars1");

  • trunk/examples/advanced/xray_fluorescence/src/XrayFluoEventAction.cc

    r807 r1230  
    254254    return   EdepDetect;
    255255   
    256 };
    257 
    258 
     256}
     257
     258

  • trunk/examples/advanced/xray_fluorescence/src/XrayFluoNistMaterials.cc

    r807 r1230  
    4343XrayFluoNistMaterials::~XrayFluoNistMaterials()
    4444{
    45   delete    dolorite;
    46   delete    HPGe;
    47   delete    mars1;
    48   delete    galactic;
    49   delete    madaBasalt;
    50   delete    icelandicBasalt;
    51   delete    GaAs;
     45  delete    dolorite;       
     46  delete    HPGe;           
     47  delete    SiLi;           
     48  delete    mars1;           
     49  delete    anorthosite;     
     50  delete    basalt;         
     51  delete    gabbro;
     52  delete    gabbroWD;
     53  delete    gabbroRF;
     54  delete    Air;             
     55  delete    Sci;             
     56  delete    Vacuum;         
     57  delete    madaBasalt;     
     58  delete    icelandicBasalt;
     59  delete    icelandicWD;
     60  delete    icelandicRF;
     61  delete    GaAs;           
     62  delete    galactic;       
     63  delete    copper;
     64  delete    hawaiianRF;
     65  delete    hawaiianWD;         
     66
     67
     68
    5269}
    5370XrayFluoNistMaterials* XrayFluoNistMaterials::instance = 0;
    … …  
    165182
    166183
    167   ///////////////////////
    168   // Iceland    Basalt //
    169   ///////////////////////
     184  ///////////////////////////////////////////
     185  // Iceland    Basalt 0029.PP.0035 sample //
     186  ///////////////////////////////////////////
    170187
    171188  elements.push_back("Si");  fractionMass.push_back(0.2313);
    … …  
    337354  mars1->AddMaterial(mars1Main, 0.9955036837);
    338355
     356  /////////////////////////////////
     357  //  Hawaiian -- WD coposition  //
     358  /////////////////////////////////
     359
     360  density = 3*g/cm3;
     361
     362  elements.push_back("Fe");    fractionMass.push_back(1.1819860E-01); 
     363  elements.push_back("Ti");    fractionMass.push_back(2.2781000E-02);
     364  elements.push_back("Ca");    fractionMass.push_back(4.5026100E-02);
     365  elements.push_back("Si");    fractionMass.push_back(2.0518860E-01); 
     366  elements.push_back("Al");    fractionMass.push_back(1.3285430E-01); 
     367  elements.push_back("Mg");    fractionMass.push_back(2.4120000E-03);
     368  elements.push_back("Na");    fractionMass.push_back(2.2257000E-02);
     369  elements.push_back("K");     fractionMass.push_back(4.9812000E-03);
     370  elements.push_back("O");     fractionMass.push_back(4.4630120E-01);
     371
     372  hawaiianWD = nistMan->ConstructNewMaterial("HawaiianWD", elements, fractionMass, density);
     373
     374  elements.clear();
     375  fractionMass.clear();
     376
     377  //////////////////////////////////
     378  //  Hawaiian -- RF composition  //
     379  //////////////////////////////////
     380
     381  density = 3*g/cm3;
     382
     383
     384  elements.push_back("Fe");    fractionMass.push_back(1.1120460E-01); 
     385  elements.push_back("Ti");    fractionMass.push_back(2.1582000E-02);
     386  elements.push_back("Ca");    fractionMass.push_back(4.3596700E-02);
     387  elements.push_back("Si");    fractionMass.push_back(2.1313440E-01); 
     388  elements.push_back("Al");    fractionMass.push_back(1.0374280E-01); 
     389  elements.push_back("Mg");    fractionMass.push_back(1.9296000E-02);
     390  elements.push_back("Na");    fractionMass.push_back(2.8192200E-02);
     391  elements.push_back("K");     fractionMass.push_back(5.8114000E-03);
     392  elements.push_back("P");     fractionMass.push_back(4.8004000E-03);
     393  elements.push_back("Mn");    fractionMass.push_back(2.3235000E-03);
     394  elements.push_back("O");     fractionMass.push_back(4.4531600E-01);
     395
     396  hawaiianRF = nistMan->ConstructNewMaterial("HawaiianRF", elements, fractionMass, density);
     397
     398  elements.clear();
     399  fractionMass.clear();
     400
     401  //////////////////////////////////
     402  //  Icelandic -- WD composition  //
     403  //////////////////////////////////
     404
     405  density = 3*g/cm3;
     406
     407
     408  elements.push_back("Si");    fractionMass.push_back(2.2949340E-01);   
     409  elements.push_back("Ti");    fractionMass.push_back(1.1990000E-02);
     410  elements.push_back("Al");    fractionMass.push_back(7.0396900E-02); 
     411  elements.push_back("Fe");    fractionMass.push_back(1.1330280E-01); 
     412  elements.push_back("Mg");    fractionMass.push_back(3.4974000E-02);
     413  elements.push_back("Ca");    fractionMass.push_back(7.5758200E-02);
     414  elements.push_back("Na");    fractionMass.push_back(1.8547500E-02); 
     415  elements.push_back("K");     fractionMass.push_back(3.3208000E-03);
     416  elements.push_back("O");     fractionMass.push_back(4.4121640E-01);
     417 
     418  icelandicWD = nistMan->ConstructNewMaterial("IcelandicWD", elements, fractionMass, density);
     419 
     420  elements.clear();
     421  fractionMass.clear();
     422
     423
     424  //////////////////////////////////
     425  //  Icelandic -- RF composition  //
     426  //////////////////////////////////
     427
     428  density = 3*g/cm3;
     429
     430
     431  elements.push_back("Si");    fractionMass.push_back(2.4304800E-01);   
     432  elements.push_back("Ti");    fractionMass.push_back(1.3788500E-02);
     433  elements.push_back("Al");    fractionMass.push_back(6.5103900E-02); 
     434  elements.push_back("Fe");    fractionMass.push_back(1.1819860E-01); 
     435  elements.push_back("Mn");    fractionMass.push_back(2.3235000E-03);
     436  elements.push_back("Mg");    fractionMass.push_back(2.3517000E-02);
     437  elements.push_back("Ca");    fractionMass.push_back(8.2190500E-02);
     438  elements.push_back("K");     fractionMass.push_back(3.3208000E-03);
     439  elements.push_back("P");     fractionMass.push_back(1.3092000E-03);
     440  elements.push_back("O");     fractionMass.push_back(4.4620000E-01);
     441 
     442  icelandicRF = nistMan->ConstructNewMaterial("IcelandicRF", elements, fractionMass, density);
     443 
     444  elements.clear();
     445  fractionMass.clear();
     446
     447  //////////////////////////////////
     448  //  Gabbro -- WD composition  //
     449  //////////////////////////////////
     450
     451  density = 3*g/cm3;
     452
     453  elements.push_back("Si");    fractionMass.push_back(1.8696000E-01);   
     454  elements.push_back("Ti");    fractionMass.push_back(2.3380500E-02);
     455  elements.push_back("Al");    fractionMass.push_back(4.6049100E-02); 
     456  elements.push_back("Fe");    fractionMass.push_back(1.2239500E-01); 
     457  elements.push_back("Mg");    fractionMass.push_back(8.3817000E-02);
     458  elements.push_back("Ca");    fractionMass.push_back(1.0720500E-01);
     459  elements.push_back("Na");    fractionMass.push_back(5.9352000E-03); 
     460  elements.push_back("K");     fractionMass.push_back(1.6604000E-03);
     461  elements.push_back("O");     fractionMass.push_back(4.2259780E-01);
     462 
     463  gabbroWD = nistMan->ConstructNewMaterial("GabbroWD", elements, fractionMass, density);
     464 
     465  elements.clear();
     466  fractionMass.clear();
     467
     468  //////////////////////////////////
     469  //  Gabbro -- RF composition  //
     470  //////////////////////////////////
     471
     472  density = 3*g/cm3;
     473
     474
     475  elements.push_back("Si");    fractionMass.push_back(1.6826400E-01);   
     476  elements.push_back("Ti");    fractionMass.push_back(2.2781000E-02);
     477  elements.push_back("Al");    fractionMass.push_back(5.8223000E-02); 
     478  elements.push_back("Fe");    fractionMass.push_back(1.2729080E-01); 
     479  elements.push_back("Mn");    fractionMass.push_back(1.5490000E-03);
     480  elements.push_back("Mg");    fractionMass.push_back(8.3817000E-02);
     481  elements.push_back("Ca");    fractionMass.push_back(1.1721080E-01);
     482  elements.push_back("Na");    fractionMass.push_back(0.0000000E+00); 
     483  elements.push_back("K");     fractionMass.push_back(1.6604000E-03);
     484  elements.push_back("P");     fractionMass.push_back(1.7456000E-03);
     485  elements.push_back("O");     fractionMass.push_back(4.1845840E-01);
     486 
     487  gabbroRF = nistMan->ConstructNewMaterial("GabbroRF", elements, fractionMass, density);
     488 
     489  elements.clear();
     490  fractionMass.clear();
     491
     492
    339493  ///////////////////////
    340494  //     Anorthosite   //
    … …  
    371525  fractionMass.clear();
    372526
     527  ////////////////////////////////////////
     528  //     Gabbro       0059.PP.0048      //
     529  ////////////////////////////////////////
     530
     531
     532  density = 3.0*g/cm3;
     533
     534  elements.push_back("Si");    fractionMass.push_back(1.8284688E-01); 
     535  elements.push_back("Ti");    fractionMass.push_back(2.2601150E-02); 
     536  elements.push_back("Al");    fractionMass.push_back(4.4831710E-02); 
     537  elements.push_back("Fe");    fractionMass.push_back(1.2578402E-01); 
     538  elements.push_back("Mn");    fractionMass.push_back(1.3166500E-03); 
     539  elements.push_back("Mg");    fractionMass.push_back(8.1706500E-02); 
     540  elements.push_back("Ca");    fractionMass.push_back(1.0506090E-01); 
     541  elements.push_back("Na");    fractionMass.push_back(5.4900600E-03); 
     542  elements.push_back("K");     fractionMass.push_back(1.4943600E-03); 
     543  elements.push_back("P");     fractionMass.push_back(3.4912000E-04); 
     544  elements.push_back("O");     fractionMass.push_back(4.0651865E-01);
     545
     546  gabbro = nistMan->ConstructNewMaterial("Gabbro", elements, fractionMass, density);
     547
     548  elements.clear();
     549  fractionMass.clear();
     550
    373551  //define gallium arsenide
    374552
    … …  
    382560  natoms.clear();
    383561
    384 
     562  /*
    385563  // define germanium
    386564 
    … …  
    388566 
    389567  elements.push_back("Ge");     natoms.push_back(1);
    390   HPGe = nistMan->ConstructNewMaterial("HPGe",elements, natoms, density);
    391 
    392   elements.clear();
    393   natoms.clear();
    394  
     568
     569  G4cout << elements[1] <<", "<<natoms[1] <<", " << elements.size() << ", " << natoms.size() << G4endl;
     570 
     571
     572  HPGe = nistMan->ConstructNewMaterial("High Purity Germanium",elements, natoms, density);
     573
     574  elements.clear();
     575  natoms.clear();
     576  */
    395577  //define scintillator
    396578
    … …  
    411593  Vacuum = new G4Material("Galactic", 1., 1.01*g/mole, density,
    412594                                       kStateGas,temperature,pressure);
     595
     596  elements.clear();
     597  natoms.clear();
    413598
    414599  //define basalt
    … …  
    422607  elements.push_back("Mg");     fractionMass.push_back(0.0590);   
    423608  elements.push_back("O");      fractionMass.push_back(0.4430);
    424 
    425 
    426609 
    427610  basalt = nistMan->ConstructNewMaterial("Basalt", elements, fractionMass, density);
    428611
     612  elements.clear();
     613  fractionMass.clear();
     614
     615 
     616  // define silicon
     617
     618  density = 2330*kg/m3;
     619
     620  // workaround for a problem in nistMan: it doesn't like material with a single element.
     621
     622  elements.push_back("Si");       natoms.push_back(1);
     623  elements.push_back("Si");       natoms.push_back(1);
     624 
     625  SiLi = nistMan->ConstructNewMaterial("SiLi",elements, natoms, density);
     626
     627  elements.clear();
     628  natoms.clear();
     629 
     630
     631  // define copper
     632
     633  density = 8920*kg/m3;
     634
     635  // workaround for a problem in nistMan: it doesn't like material with a single element.
     636  elements.push_back("Cu");       natoms.push_back(1);
     637  elements.push_back("Cu");       natoms.push_back(1);
     638
     639  copper = nistMan->ConstructNewMaterial("Cu",elements, natoms, density);
     640
     641  elements.clear();
     642  natoms.clear();
     643
    429644  G4cout << *(G4Material::GetMaterialTable()) << G4endl;
    430645}

  • trunk/examples/advanced/xray_fluorescence/src/XrayFluoPhysicsList.cc

    r807 r1230  
    3636// -------------------------------------------------------------------
    3737
     38#include "G4ParticleDefinition.hh"
     39#include "G4ParticleTypes.hh"
     40#include "G4ProcessManager.hh"
    3841#include "XrayFluoPhysicsList.hh"
    3942#include "XrayFluoPhysicsListMessenger.hh"
    … …  
    4346
    4447/////////////////////////////////////////
    45 //#include "G4LeptonConstructor.hh"
    46 //#include "G4BosonConstructor.hh"
     48#include "G4LeptonConstructor.hh"
     49#include "G4BosonConstructor.hh"
    4750//#include "G4MesonConstructor.hh"
     51#include "G4IonConstructor.hh"
    4852#include "G4BaryonConstructor.hh"
    4953/////////////////////////////////////////
    5054
    51 
     55/*
    5256#include "G4ParticleDefinition.hh"
    5357#include "G4ParticleWithCuts.hh"
    … …  
    5660#include "G4ParticleTable.hh"
    5761#include "G4ios.hh"
    58 
    59 #include "G4Region.hh"
    60 #include "G4RegionStore.hh"
     62*/
     63
     64//#include "G4Region.hh"
     65//#include "G4RegionStore.hh"
    6166
    6267
    … …  
    6974  pDet = p;
    7075
    71   //  SetGELowLimit(250*eV);
     76  SetGELowLimit(250*eV);
    7277
    7378  defaultCutValue = 10e-6*mm;
    … …  
    7984  physicsListMessenger = new XrayFluoPhysicsListMessenger(this);
    8085
     86  /*
    8187  G4String regName = "SampleRegion";
    8288  G4double cutValue = 0.000001 * mm; 
    … …  
    8692  cuts->SetProductionCut(cutValue);
    8793  reg->SetProductionCuts(cuts);
    88 
     94  */
    8995
    9096
    … …  
    178184{
    179185//  Ions
    180   G4Alpha::AlphaDefinition();
     186  G4IonConstructor ions;
     187  ions.ConstructParticle();
     188  //  G4Alpha::AlphaDefinition();
    181189}
    182190
    … …  
    189197//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    190198
     199
     200#include "G4PhotoElectricEffect.hh"
     201#include "G4LivermorePhotoElectricModel.hh"
     202
     203#include "G4ComptonScattering.hh"
     204#include "G4LivermoreComptonModel.hh"
     205
     206#include "G4GammaConversion.hh"
     207#include "G4LivermoreGammaConversionModel.hh"
     208
     209#include "G4RayleighScattering.hh"
     210#include "G4LivermoreRayleighModel.hh"
     211/*
    191212#include "G4LowEnergyCompton.hh"
    192213#include "G4LowEnergyGammaConversion.hh"
    193214#include "G4LowEnergyPhotoElectric.hh"
    194215#include "G4LowEnergyRayleigh.hh"
    195 
    196 // e+
     216*/
     217
     218// e+-
     219/*
    197220#include "G4MultipleScattering.hh"
    198221#include "G4eIonisation.hh"
    199222#include "G4eBremsstrahlung.hh"
    200 #include "G4eplusAnnihilation.hh"
     223
    201224
    202225#include "G4LowEnergyIonisation.hh"
    203226#include "G4LowEnergyBremsstrahlung.hh"
     227*/
     228
     229#include "G4eMultipleScattering.hh"
     230#include "G4UniversalFluctuation.hh"
     231
     232#include "G4eIonisation.hh"
     233#include "G4LivermoreIonisationModel.hh"
     234
     235#include "G4eBremsstrahlung.hh"
     236#include "G4LivermoreBremsstrahlungModel.hh"
     237
     238#include "G4eplusAnnihilation.hh"
     239
    204240#include "G4hLowEnergyIonisation.hh"
     241#include "G4hMultipleScattering.hh"
     242
     243// options
     244
     245#include "G4LossTableManager.hh"
     246#include "G4EmProcessOptions.hh"
    205247
    206248//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    217259
    218260      // gamma         
     261      /*
    219262      pmanager->AddDiscreteProcess(new G4LowEnergyCompton);
    220263     
    221264      LePeprocess = new G4LowEnergyPhotoElectric();
    222265
    223       LePeprocess->ActivateAuger(true);
    224       LePeprocess->SetCutForLowEnSecPhotons(0.250 * keV);
    225       LePeprocess->SetCutForLowEnSecElectrons(0.250 * keV);
    226266
    227267      pmanager->AddDiscreteProcess(LePeprocess);
    … …  
    229269      pmanager->AddDiscreteProcess(new G4LowEnergyRayleigh("Rayleigh"));
    230270     
     271      */
     272
     273      G4PhotoElectricEffect* thePhotoElectricEffect = new G4PhotoElectricEffect();
     274      G4LivermorePhotoElectricModel* theLivermorePhotoElectricModel = new G4LivermorePhotoElectricModel();
     275
     276      theLivermorePhotoElectricModel->ActivateAuger(true);
     277      theLivermorePhotoElectricModel->SetCutForLowEnSecPhotons(0.010 * keV);
     278      theLivermorePhotoElectricModel->SetCutForLowEnSecElectrons(0.010 * keV);
     279
     280      thePhotoElectricEffect->AddEmModel(0, theLivermorePhotoElectricModel);
     281      pmanager->AddDiscreteProcess(thePhotoElectricEffect);
     282
     283
     284      G4ComptonScattering* theComptonScattering = new G4ComptonScattering();
     285      G4LivermoreComptonModel* theLivermoreComptonModel = new G4LivermoreComptonModel();
     286      theComptonScattering->AddEmModel(0, theLivermoreComptonModel);
     287      pmanager->AddDiscreteProcess(theComptonScattering);
     288
     289      G4GammaConversion* theGammaConversion = new G4GammaConversion();
     290      G4LivermoreGammaConversionModel* theLivermoreGammaConversionModel = new G4LivermoreGammaConversionModel();
     291      theGammaConversion->AddEmModel(0, theLivermoreGammaConversionModel);
     292      pmanager->AddDiscreteProcess(theGammaConversion);
     293
     294      G4RayleighScattering* theRayleigh = new G4RayleighScattering();
     295      G4LivermoreRayleighModel* theRayleighModel = new G4LivermoreRayleighModel();
     296      theRayleigh->AddEmModel(0, theRayleighModel);
     297      pmanager->AddDiscreteProcess(theRayleigh);
     298
     299
     300
    231301    } else if (particleName == "e-") {
    232302      //electron
    233       pmanager->AddProcess(new G4MultipleScattering,-1, 1,1);
    234      
     303      /*   
     304      pmanager->AddProcess(new G4MultipleScattering(),-1, 1,1);
     305
    235306      LeIoprocess = new G4LowEnergyIonisation("IONI");
    236307      LeIoprocess->ActivateAuger(true);
    237308      //eIoProcess = new G4eIonisation("stdIONI");
    238       LeIoprocess->SetCutForLowEnSecPhotons(0.1*keV);
    239       LeIoprocess->SetCutForLowEnSecElectrons(0.1*keV);
     309      LeIoprocess->SetCutForLowEnSecPhotons(0.01*keV);
     310      LeIoprocess->SetCutForLowEnSecElectrons(0.01*keV);
    240311      pmanager->AddProcess(LeIoprocess, -1,  2, 2);
    241312
    242313      LeBrprocess = new G4LowEnergyBremsstrahlung();
    243314      pmanager->AddProcess(LeBrprocess, -1, -1, 3);
     315      */     
     316     
     317      G4eMultipleScattering* msc = new G4eMultipleScattering();
     318      pmanager->AddProcess(msc,                   -1, 1, 1);
    244319     
     320      // Ionisation
     321      G4eIonisation* eIoni = new G4eIonisation();
     322      G4LivermoreIonisationModel* ioniModel = new G4LivermoreIonisationModel();
     323      eIoni->AddEmModel(0, ioniModel, new G4UniversalFluctuation() );
     324
     325      //      ioniModel->SetCutForLowEnSecPhotons(0.01*keV);
     326      //      ioniModel->SetCutForLowEnSecElectrons(0.01*keV);
     327
     328      eIoni->SetStepFunction(0.2, 100*um); //     
     329      pmanager->AddProcess(eIoni,                 -1, 2, 2);
     330     
     331      // Bremsstrahlung
     332      G4eBremsstrahlung* eBrem = new G4eBremsstrahlung();
     333      eBrem->AddEmModel(0, new G4LivermoreBremsstrahlungModel());
     334      pmanager->AddProcess(eBrem,                 -1,-3, 3);
     335
     336
     337
     338
    245339      } else if (particleName == "e+") {
    246340      //positron
    247       pmanager->AddProcess(new G4MultipleScattering,-1, 1,1);
    248       pmanager->AddProcess(new G4eIonisation,      -1, 2,2);
    249       pmanager->AddProcess(new G4eBremsstrahlung,   -1,-1,3);
    250       pmanager->AddProcess(new G4eplusAnnihilation,  0,-1,4);
     341      pmanager->AddProcess(new G4eMultipleScattering(),-1, 1,1);
     342      pmanager->AddProcess(new G4eIonisation(),      -1, 2,2);
     343      pmanager->AddProcess(new G4eBremsstrahlung(),   -1,-1,3);
     344      pmanager->AddProcess(new G4eplusAnnihilation(),  0,-1,4);
    251345     
    252346    } 
    253347    else if (particleName == "proton") {
    254348      //proton
    255       G4hLowEnergyIonisation* hIoni = new G4hLowEnergyIonisation;
     349      G4hLowEnergyIonisation* hIoni = new G4hLowEnergyIonisation();
    256350      hIoni->SetFluorescence(true);
    257       pmanager->AddProcess(new G4MultipleScattering,-1,1,1);
     351      pmanager->AddProcess(new G4hMultipleScattering(),-1,1,1);
    258352      pmanager->AddProcess(hIoni,-1, 2,2);
    259353    }
    … …  
    261355      {
    262356       
    263         pmanager->AddProcess(new G4MultipleScattering,-1,1,1);
     357        pmanager->AddProcess(new G4hMultipleScattering(),-1,1,1);
    264358        G4hLowEnergyIonisation* iIon = new G4hLowEnergyIonisation() ;
    265359        pmanager->AddProcess(iIon,-1,2,2);
    … …  
    332426void XrayFluoPhysicsList::SetCuts(){
    333427
    334    SetCutValue(cutForGamma,"gamma");
    335    SetCutValue(cutForElectron,"e-");
    336    SetCutValue(cutForElectron,"e+");
    337    SetCutValue(cutForProton, "proton");
    338    //SetCutValueForOthers(cutForProton);
    339    if (verboseLevel>0) DumpCutValuesTable();
     428  G4double lowlimit=20*eV;
     429  G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange(lowlimit, 100.*GeV);
     430  SetCutValue(cutForGamma,"gamma");
     431  SetCutValue(cutForElectron,"e-");
     432  SetCutValue(cutForElectron,"e+");
     433  SetCutValue(cutForProton, "proton");
     434  //SetCutValueForOthers(cutForProton);
     435  if (verboseLevel>0) DumpCutValuesTable();
    340436}
    341437
    … …  
    343439
    344440void XrayFluoPhysicsList::SetLowEnSecPhotCut(G4double cut){
    345  
    346   G4cout<<"Low energy secondary photons cut is now set to: "<<cut/MeV<<" (MeV)"<<G4endl;
    347   G4cout<<"for processes LowEnergyBremsstrahlung, LowEnergyPhotoElectric, LowEnergyIonisation"<<G4endl;
    348   LeBrprocess->SetCutForLowEnSecPhotons(cut);
    349   LePeprocess->SetCutForLowEnSecPhotons(cut);
    350   LeIoprocess->SetCutForLowEnSecPhotons(cut);
     441
     442  cut=0;
     443 
     444//  G4cout<<"Low energy secondary photons cut is now set to: "<<cut/MeV<<" (MeV)"<<G4endl;
     445//  G4cout<<"for processes LowEnergyBremsstrahlung, LowEnergyPhotoElectric, LowEnergyIonisation"<<G4endl;
     446//  LeBrprocess->SetCutForLowEnSecPhotons(cut);
     447//  LePeprocess->SetCutForLowEnSecPhotons(cut);
     448//  LeIoprocess->SetCutForLowEnSecPhotons(cut);
    351449}
    352450
    353451void XrayFluoPhysicsList::SetLowEnSecElecCut(G4double cut){
    354452 
    355   G4cout<<"Low energy secondary electrons cut is now set to: "<<cut/MeV<<" (MeV)"<<G4endl;
    356  
    357   G4cout<<"for processes LowEnergyIonisation"<<G4endl;
    358  
    359   LeIoprocess->SetCutForLowEnSecElectrons(cut);
     453  cut = 0;
     454//  G4cout<<"Low energy secondary electrons cut is now set to: "<<cut/MeV<<" (MeV)"<<G4endl;
     455// 
     456//  G4cout<<"for processes LowEnergyIonisation"<<G4endl;
     457// 
     458//  LeIoprocess->SetCutForLowEnSecElectrons(cut);
    360459}
    361460//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
    … …  
    397496  SetCutValue(cut, "proton");
    398497 
    399 //   part = theXrayFluoParticleTable->FindParticle("gamma");
    400 //   cut = G4EnergyLossTables::GetRange(part,val,currMat);
    401 //   SetCutValue(cut, "gamma");
    402 
    403 }
    404 
    405 
    406 
    407 
    408 
    409 
    410 
    411 
    412 
    413 
    414 
    415 
    416 
    417 
    418 
    419 
    420 
    421 
    422 
    423 
    424 
    425 
    426 
    427 
     498  part = theXrayFluoParticleTable->FindParticle("gamma");
     499  cut = G4EnergyLossTables::GetRange(part,val,currMat);
     500  SetCutValue(cut, "gamma");
     501
     502}
     503
     504
     505
     506
     507
     508
     509
     510
     511
     512
     513
     514
     515
     516
     517
     518
     519
     520
     521
     522
     523
     524
     525
     526

  • trunk/examples/advanced/xray_fluorescence/src/XrayFluoPrimaryGeneratorAction.cc

    r807 r1230  
    101101  particleTypes = analysis->GetEmittedParticleTypes();
    102102  detectorPosition = XrayFluoDetector->GetDetectorPosition();
    103   detectorPosition.setR(detectorPosition.r()-(5.*cm)); // 5 cm dietro
     103  detectorPosition.setR(detectorPosition.r()-(5.*cm)); // 5 cm before the detector, so in front of it.
    104104#endif
    105105
    … …  
    233233    }
    234234
     235  // using prevoiously genereated emissions from sample.....
    235236
    236237  if (phaseSpaceGunFlag){
    … …  
    254255
    255256    particleGun->SetParticleEnergy(energy);
    256 
    257 
    258 
    259257    particleGun->SetParticleDefinition(particle);
     258
     259
    260260  }
    261261

  • trunk/examples/advanced/xray_fluorescence/src/XrayFluoRunAction.cc

    r807 r1230  
    6363 
    6464 
    65   ReadData(MeV,"mercury_flx_solmin");
     65  ReadData(keV,"spec10");
    6666  //ReadResponse("SILIresponse");
    6767 
    6868  G4double minGamma = 0.*keV;
    6969  G4double maxGamma = 10. *keV;
    70   G4int nBinsGamma = 5;
     70  G4int nBinsGamma = 6;
    7171 
    7272
    7373  dataGammaSet = normalization.Normalize(minGamma, maxGamma, nBinsGamma,
    74                                   "FlatSpectrum");
     74                                  "M_flare");
    7575 
    7676

  • trunk/examples/advanced/xray_fluorescence/src/XrayFluoSiLiDetectorType.cc

    r807 r1230  
    4848
    4949XrayFluoSiLiDetectorType::XrayFluoSiLiDetectorType():
    50   detectorMaterial("Silicon"),efficiencySet(0)
     50  detectorMaterial("SiLi"),efficiencySet(0)
    5151{
    5252  LoadResponseData("SILIresponse");
Note: See TracChangeset for help on using the changeset viewer.