Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Validation/Geometry: Tracker Air Categorization #25231

Merged
merged 4 commits into from Nov 20, 2018
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
2 changes: 1 addition & 1 deletion Validation/Geometry/interface/MaterialBudgetCategorizer.h
Expand Up @@ -15,7 +15,7 @@
class MaterialBudgetCategorizer {

public:
MaterialBudgetCategorizer();
MaterialBudgetCategorizer(std::string mode);

int volume(std::string s){return theVolumeMap[s];}
int material(std::string s){return theMaterialMap[s];}
Expand Down
143 changes: 31 additions & 112 deletions Validation/Geometry/src/MaterialBudgetCategorizer.cc
Expand Up @@ -12,12 +12,7 @@
#include <fstream>
#include <vector>

MaterialBudgetCategorizer::MaterialBudgetCategorizer()
{
buildMaps();
}

void MaterialBudgetCategorizer::buildMaps()
MaterialBudgetCategorizer::MaterialBudgetCategorizer(std::string mode)
{
//----- Build map volume name - volume index
G4LogicalVolumeStore* lvs = G4LogicalVolumeStore::GetInstance();
Expand All @@ -26,110 +21,29 @@ void MaterialBudgetCategorizer::buildMaps()
for (ite = lvs->begin(); ite != lvs->end(); ite++) {
theVolumeMap[(*ite)->GetName()] = ii++;
}

//----- Build map material name - volume index
const G4MaterialTable* matTable = G4Material::GetMaterialTable();
G4int matSize = matTable->size();
for( ii = 0; ii < matSize; ii++ ) {
theMaterialMap[ (*matTable)[ii]->GetName()] = ii+1;
}

//----- Build map material name
for( int ii = 0; ii < matSize; ii++ ) {
float sup,sen,cab,col,ele,oth,air;
sup=sen=cab=col=ele=0.;
oth=1.;
air=0;
edm::LogInfo("MaterialBudget") << "MaterialBudgetCategorizer: Material " << (*matTable)[ii]->GetName() << " prepared";
if((*matTable)[ii]->GetName()=="Air"
||
(*matTable)[ii]->GetName()=="Vacuum"
) {
air=1.000;
oth=0.000;
}
// actually this class does not work if there are spaces into material names --> Recover material properties
if((*matTable)[ii]->GetName()=="Carbon fibre str.") {
sup=1.000;
oth=0.000;
}
// X0
theX0Map[ (*matTable)[ii]->GetName() ].push_back(sup); // sup
theX0Map[ (*matTable)[ii]->GetName() ].push_back(sen); // sen
theX0Map[ (*matTable)[ii]->GetName() ].push_back(cab); // cab
theX0Map[ (*matTable)[ii]->GetName() ].push_back(col); // col
theX0Map[ (*matTable)[ii]->GetName() ].push_back(ele); // ele
theX0Map[ (*matTable)[ii]->GetName() ].push_back(oth); // oth
theX0Map[ (*matTable)[ii]->GetName() ].push_back(air); // air
// L0
theL0Map[ (*matTable)[ii]->GetName() ].push_back(sup); // sup
theL0Map[ (*matTable)[ii]->GetName() ].push_back(sen); // sen
theL0Map[ (*matTable)[ii]->GetName() ].push_back(cab); // cab
theL0Map[ (*matTable)[ii]->GetName() ].push_back(col); // col
theL0Map[ (*matTable)[ii]->GetName() ].push_back(ele); // ele
theL0Map[ (*matTable)[ii]->GetName() ].push_back(oth); // oth
theL0Map[ (*matTable)[ii]->GetName() ].push_back(air); // air
}

//----- Build map material name - X0 contributes
edm::LogInfo("MaterialBudget") << "MaterialBudgetCategorizer: Fill X0 Map";
std::string theMaterialX0FileName = edm::FileInPath("Validation/Geometry/data/trackerMaterials.x0").fullPath();
buildCategoryMap(theMaterialX0FileName, theX0Map);
//For the HGCal
std::string thehgcalMaterialX0FileName = edm::FileInPath("Validation/Geometry/data/hgcalMaterials.x0").fullPath();
buildHGCalCategoryMap(thehgcalMaterialX0FileName, theHGCalX0Map);
//----- Build map material name - L0 contributes
edm::LogInfo("MaterialBudget") << "MaterialBudgetCategorizer: Fill L0 Map";
std::string theMaterialL0FileName = edm::FileInPath("Validation/Geometry/data/trackerMaterials.l0").fullPath();
buildCategoryMap(theMaterialL0FileName, theL0Map);
//For the HGCal
std::string thehgcalMaterialL0FileName = edm::FileInPath("Validation/Geometry/data/hgcalMaterials.l0").fullPath();
buildHGCalCategoryMap(thehgcalMaterialL0FileName, theHGCalL0Map);
// summary of all the materials loaded
edm::LogInfo("MaterialBudget") << "MaterialBudgetCategorizer: Material Summary Starts";
G4EmCalculator calc;
for( ii = 0; ii < matSize; ii++ ) {
edm::LogInfo("MaterialBudget")
<< "MaterialBudgetCateogorizer: Material " << (*matTable)[ii]->GetName()
<< std::endl
<< "\t density = " << G4BestUnit((*matTable)[ii]->GetDensity(),"Volumic Mass")
<< std::endl
<< "\t X0 = " << (*matTable)[ii]->GetRadlen() << " mm"
<< std::endl
<< "\t Energy threshold for photons for 100 mm range = "
<< G4BestUnit(calc.ComputeEnergyCutFromRangeCut(100, G4String("gamma"), (*matTable)[ii]->GetName()) , "Energy")
<< std::endl
<< " SUP " << theX0Map[ (*matTable)[ii]->GetName() ][0]
<< " SEN " << theX0Map[ (*matTable)[ii]->GetName() ][1]
<< " CAB " << theX0Map[ (*matTable)[ii]->GetName() ][2]
<< " COL " << theX0Map[ (*matTable)[ii]->GetName() ][3]
<< " ELE " << theX0Map[ (*matTable)[ii]->GetName() ][4]
<< " OTH " << theX0Map[ (*matTable)[ii]->GetName() ][5]
<< " AIR " << theX0Map[ (*matTable)[ii]->GetName() ][6]
<< std::endl
<< "\t Lambda0 = " << (*matTable)[ii]->GetNuclearInterLength() << " mm"
<< std::endl
<< " SUP " << theL0Map[ (*matTable)[ii]->GetName() ][0]
<< " SEN " << theL0Map[ (*matTable)[ii]->GetName() ][1]
<< " CAB " << theL0Map[ (*matTable)[ii]->GetName() ][2]
<< " COL " << theL0Map[ (*matTable)[ii]->GetName() ][3]
<< " ELE " << theL0Map[ (*matTable)[ii]->GetName() ][4]
<< " OTH " << theL0Map[ (*matTable)[ii]->GetName() ][5]
<< " AIR " << theL0Map[ (*matTable)[ii]->GetName() ][6]
<< std::endl;
if( theX0Map[ (*matTable)[ii]->GetName() ][5] == 1 || theL0Map[ (*matTable)[ii]->GetName() ][5] == 1 )
edm::LogWarning("MaterialBudget")
<< "MaterialBudgetCategorizer Material with no category: " << (*matTable)[ii]->GetName();

if ( mode.compare("Tracker") == 0 ) {
std::string theMaterialX0FileName = edm::FileInPath("Validation/Geometry/data/trackerMaterials.x0").fullPath();
buildCategoryMap(theMaterialX0FileName, theX0Map);
std::string theMaterialL0FileName = edm::FileInPath("Validation/Geometry/data/trackerMaterials.l0").fullPath();
buildCategoryMap(theMaterialL0FileName, theL0Map);
} else if ( mode.compare("HGCal") == 0 ){
std::string thehgcalMaterialX0FileName = edm::FileInPath("Validation/Geometry/data/hgcalMaterials.x0").fullPath();
buildHGCalCategoryMap(thehgcalMaterialX0FileName, theHGCalX0Map);
std::string thehgcalMaterialL0FileName = edm::FileInPath("Validation/Geometry/data/hgcalMaterials.l0").fullPath();
buildHGCalCategoryMap(thehgcalMaterialL0FileName, theHGCalL0Map);
}
}

void MaterialBudgetCategorizer::buildCategoryMap(std::string theMaterialFileName, std::map<std::string,std::vector<float> >& theMap) {

void MaterialBudgetCategorizer::buildCategoryMap(std::string theMaterialFileName,
std::map<std::string,std::vector<float> >& theMap)
{

std::ifstream theMaterialFile(theMaterialFileName);

if (!theMaterialFile)
cms::Exception("LogicError") << " File not found " << theMaterialFileName;

// fill everything as "other"
float sup,sen,cab,col,ele,oth,air;
sup=sen=cab=col=ele=0.;

Expand All @@ -138,6 +52,10 @@ void MaterialBudgetCategorizer::buildCategoryMap(std::string theMaterialFileName
while(theMaterialFile) {
theMaterialFile >> materialName;
theMaterialFile >> sup >> sen >> cab >> col >> ele;

if (materialName[0] == '#') //Ignore comments
continue;

oth = 0.000;
air = 0.000;
theMap[materialName].clear(); // clear before re-filling
Expand All @@ -148,15 +66,16 @@ void MaterialBudgetCategorizer::buildCategoryMap(std::string theMaterialFileName
theMap[materialName].push_back(ele); // ele
theMap[materialName].push_back(oth); // oth
theMap[materialName].push_back(air); // air
LogDebug("MaterialBudget")
<< "MaterialBudgetCategorizer: Material " << materialName << " filled "
<< " SUP " << sup
<< " SEN " << sen
<< " CAB " << cab
<< " COL " << col
<< " ELE " << ele
<< " OTH " << oth
<< " AIR " << air;
edm::LogInfo("MaterialBudget")
<< "MaterialBudgetCategorizer: Material " << materialName << " filled: "
<< "\n\tSUP " << sup
<< "\n\tSEN " << sen
<< "\n\tCAB " << cab
<< "\n\tCOL " << col
<< "\n\tELE " << ele
<< "\n\tOTH " << oth
<< "\n\tAIR " << air
<< "\n\tAdd up to: " << sup + sen + cab + col + ele + oth + air;
}
}

Expand Down
23 changes: 15 additions & 8 deletions Validation/Geometry/src/MaterialBudgetData.cc
Expand Up @@ -30,7 +30,11 @@ void MaterialBudgetData::dataStartTrack( const G4Track* aTrack )
const G4ThreeVector& dir = aTrack->GetMomentum() ;

if( myMaterialBudgetCategorizer == nullptr){
myMaterialBudgetCategorizer = std::make_unique<MaterialBudgetCategorizer>();
if(isHGCal){
myMaterialBudgetCategorizer = std::make_unique<MaterialBudgetCategorizer>("HGCal");
} else {
myMaterialBudgetCategorizer = std::make_unique<MaterialBudgetCategorizer>("Tracker");
}
}

theStepN=0;
Expand Down Expand Up @@ -226,13 +230,16 @@ void MaterialBudgetData::dataPerStep( const G4Step* aStep )

if(!isCtgOk)
{
theOtherFractionMB = 1;
theOtherFractionIL = 1;

edm::LogWarning("MaterialBudget")
<< "MaterialBudgetData: Material forced to 'Other': " << materialName
<< " in volume " << volumeName << ". Check Categorization.";

if(materialName.compare("Air") == 0){
theAirFractionMB = 1;
theAirFractionIL = 1;
} else {
theOtherFractionMB = 1;
theOtherFractionIL = 1;
edm::LogWarning("MaterialBudget")
<< "MaterialBudgetData: Material forced to 'Other': " << materialName
<< " in volume " << volumeName << ". Check Categorization.";
}
}
else
{
Expand Down