Skip to content
Snippets Groups Projects
Commit b6d194d6 authored by Thomas Latham's avatar Thomas Latham
Browse files

Add VELO-specific version of eta-phi scan

parent 1ef25c1b
No related tags found
No related merge requests found
Pipeline #6031668 failed with stages
in 4 minutes and 43 seconds
......@@ -21,6 +21,7 @@
#include "TH2F.h"
#include "TImage.h"
#include "TMath.h"
#include "TRandom3.h"
#include <map>
#include <string>
......@@ -474,7 +475,7 @@ static long tgeo_rad_scan_etaphi( dd4hep::Detector& /* description */, int argc,
auto nav = gGeoManager->GetCurrentNavigator();
std::vector<TH2F*> hrad_eta;
for ( size_t i = 0; i < zPoints.size() - 1; ++i ) {
for ( std::size_t i = 0; i < zPoints.size() - 1; ++i ) {
std::string name{"rad length eta vs. phi" + std::to_string( zPoints[i] )};
if ( i > 0 ) {
name = "rad length eta vs. phi from " + std::to_string( zPoints[i - 1] ) + " to " + std::to_string( zPoints[i] );
......@@ -494,11 +495,11 @@ static long tgeo_rad_scan_etaphi( dd4hep::Detector& /* description */, int argc,
double x{0.0};
double y{0.0};
for ( size_t k = 0; k < nPoints; ++k ) {
for ( std::size_t k = 0; k < nPoints; ++k ) {
eta = 1.5 + ( float( k ) + 0.5 ) * 3.5 / float( nPoints );
theta = 2.0 * atan( exp( -eta ) );
for ( size_t j = 0; j < nPoints; ++j ) {
for ( std::size_t j = 0; j < nPoints; ++j ) {
phi = -180. + ( float( j ) + 0.5 ) * 360. / float( nPoints );
phirad = phi / 180. * TMath::Pi();
......@@ -535,7 +536,7 @@ static long tgeo_rad_scan_etaphi( dd4hep::Detector& /* description */, int argc,
}
hrad_eta[zPoints.size() - 1]->Fill( phi, eta, totalRadLength );
for ( size_t i = 0; i < zPoints.size() - 1; ++i ) {
for ( std::size_t i = 0; i < zPoints.size() - 1; ++i ) {
const double r{zPoints[i] / cos( theta )};
x = r * cos( phirad ) * sin( theta );
y = r * sin( phirad ) * sin( theta );
......@@ -580,3 +581,331 @@ static long tgeo_rad_scan_etaphi( dd4hep::Detector& /* description */, int argc,
return 1;
}
DECLARE_APPLY( rad_scan_etaphi, tgeo_rad_scan_etaphi )
std::string tagMatName( const TString& matName, const TString& nodeName ) {
TString sanitisedName { matName };
sanitisedName.ReplaceAll(":","_");
if ( sanitisedName.Contains( "Aluminium" ) ) {
// Allow separate tagging of things made of Aluminium
if ( nodeName.Contains( "pvRFFoil" ) ) {
sanitisedName += "_RFFoil";
} else if ( nodeName.Contains( "RFBox" ) ) {
sanitisedName += "_RFBox";
} else if ( nodeName.Contains( "Foot" ) ) {
sanitisedName += "_Foot";
}
} else if ( sanitisedName.Contains( "Silicon" ) ) {
// Allow separate tagging of things made of Silicon
if ( nodeName.Contains( "pvDet" ) ) {
sanitisedName += "_Sensor";
} else if ( nodeName.Contains( "pvGuardRing" ) ) {
sanitisedName += "_GuardRing";
} else if ( nodeName.Contains( "pvChip" ) ) {
sanitisedName += "_ASIC";
} else if ( nodeName.Contains( "pvSubstrate" ) ) {
sanitisedName += "_Substrate";
}
}
return sanitisedName.Data();
}
void fillPerMaterialHistograms( std::map<std::string,TH2*>& hists, const TH2* histTemplate, const std::map<std::string,double>& radlengths, const double phi_deg, const double eta, const double ipSampleWeight )
{
for ( auto& [ matName, matRadLen ] : radlengths ) {
if ( hists.find( matName ) == hists.end() ) {
const std::string histName { std::string{histTemplate->GetName()} + "_" + matName };
const std::string histTitle { matName + " " + histTemplate->GetTitle() };
TH2* tmpHist = static_cast<TH2*>( histTemplate->Clone( histName.c_str() ) );
tmpHist->SetTitle( histTitle.c_str() );
tmpHist->Reset();
hists[ matName ] = tmpHist;
}
hists[ matName ]->Fill( phi_deg, eta, matRadLen * ipSampleWeight );
}
}
struct MaterialEntry {
// name
std::string name;
// index
std::size_t index;
// path length (in mm)
double thickness;
// percentage of X_0
double radLength;
};
struct HitEntry {
// index of corresponding entry in material list
std::size_t index;
// transverse radius sqrt(x^2 + y^2)
double rho;
// z position
double z;
};
static long tgeo_rad_scan_etaphi_velo( dd4hep::Detector& /* description */, int argc, char** argv ) {
unsigned int nPointsEta { 300 };
const double etaLow { 2.0 };
const double etaHigh { 5.0 };
const double etaRange { etaHigh - etaLow };
unsigned int nPointsPhi { 360 };
const double phiLow { -180.0 };
const double phiHigh { +180.0 };
const double phiRange { phiHigh - phiLow };
double maxRadius = 290.0;
double maxDistance = 835.0;
double meanIP = 0.0;
double sigmaIP = 53.0;
unsigned int nIPSamples = 100;
unsigned int seed = 0;
bool verbose = false;
bool help = false;
bool error = false;
for ( int i = 0; i < argc && argv[i]; ++i ) {
if ( 0 == ::strncmp( "-help", argv[i], 4 ) ) { help = true; }
}
const std::vector<std::string> cmdLine{ argv, argv+argc };
const std::size_t nArgs { cmdLine.size() };
//std::cout << "We got " << nArgs << " arguments : " << std::endl;
//for ( auto& arg : cmdLine ) {
//std::cout << arg << std::endl;
//}
if ( nArgs == 0 ) {
// This is ok but we don't need to do anything!
} else if ( nArgs == 7 ) {
nPointsEta = atoi( argv[0] );
nPointsPhi = atoi( argv[1] );
maxRadius = atof( argv[2] );
maxDistance = atof( argv[3] );
nIPSamples = atoi( argv[4] );
seed = atoi( argv[5] );
verbose = atoi( argv[6] );
} else {
error = true;
}
if ( help || error ) {
/// Help printout describing the basic command line interface
std::cout
<< "Usage: -plugin <name> npoints_eta npoints_phi max_radius max_distance n_ip_samples seed verbose\n"
" name: factory name rad_scan \n"
" -help Show this help. \n"
"\tArguments given: "
<< dd4hep::arguments( argc, argv ) << std::endl;
::exit( EINVAL );
}
auto nav = gGeoManager->GetCurrentNavigator();
TH2* histRadLen = new TH2F("RadLen", "Total radiation length [%]" , nPointsPhi, phiLow, phiHigh, nPointsEta, etaLow, etaHigh);
TH2* histBefore1stHitRadLen = new TH2F("RadLen_Before1stHit", "Radiation length [%] before 1st hit" , nPointsPhi, phiLow, phiHigh, nPointsEta, etaLow, etaHigh);
TH2* histBefore2ndHitRadLen = new TH2F("RadLen_Before2ndHit", "Radiation length [%] before 2nd hit" , nPointsPhi, phiLow, phiHigh, nPointsEta, etaLow, etaHigh);
TH2* histBeforeLastHitRadLen = new TH2F("RadLen_BeforeLastHit", "Radiation length [%] before last hit", nPointsPhi, phiLow, phiHigh, nPointsEta, etaLow, etaHigh);
TH2* histHitsPerTrack = new TH2F("HitsPerTrack", "Number of hits per track" , nPointsPhi, phiLow, phiHigh, nPointsEta, etaLow, etaHigh);
TH2* histRadius1stHit = new TH2F("Radius_1stHit", "Radius of the 1st hit [mm]" , nPointsPhi, phiLow, phiHigh, nPointsEta, etaLow, etaHigh);
TH2* histRadius2ndHit = new TH2F("Radius_2ndHit", "Radius of the 2nd hit [mm]" , nPointsPhi, phiLow, phiHigh, nPointsEta, etaLow, etaHigh);
TH2* histRadiusLastHit = new TH2F("Radius_LastHit", "Radius of the last hit [mm]" , nPointsPhi, phiLow, phiHigh, nPointsEta, etaLow, etaHigh);
TH2* histZ1stHit = new TH2F("Z_1stHit", "Z of the 1st hit [mm]" , nPointsPhi, phiLow, phiHigh, nPointsEta, etaLow, etaHigh);
TH2* histZ2ndHit = new TH2F("Z_2ndHit", "Z of the 2nd hit [mm]" , nPointsPhi, phiLow, phiHigh, nPointsEta, etaLow, etaHigh);
TH2* histZLastHit = new TH2F("Z_LastHit", "Z of the last hit [mm]" , nPointsPhi, phiLow, phiHigh, nPointsEta, etaLow, etaHigh);
histRadLen->SetMaximum(50.0);
histRadLen->SetMinimum(0.0);
histBefore1stHitRadLen->SetMaximum(20.0);
histBefore1stHitRadLen->SetMinimum(0.0);
histBefore2ndHitRadLen->SetMaximum(20.0);
histBefore2ndHitRadLen->SetMinimum(0.0);
histBeforeLastHitRadLen->SetMaximum(20.0);
histBeforeLastHitRadLen->SetMinimum(0.0);
histHitsPerTrack->SetMaximum(20.0);
histHitsPerTrack->SetMinimum( 0.0);
std::map<std::string,TH2*> matHist_Total;
std::map<std::string,TH2*> matHist_1stHit;
std::map<std::string,TH2*> matHist_2ndHit;
std::map<std::string,TH2*> matHist_LastHit;
ROOT::Math::XYZPoint ip{ 0.0, 0.0, meanIP };
const double ipSampleWeight { 1.0 / nIPSamples };
double eta{0.};
//double theta_deg{0.};
double theta_rad{0.};
double cosTheta{0.};
double sinTheta{0.};
double phi_deg{0.};
double phi_rad{0.};
double cosPhi{0.};
double sinPhi{0.};
// Store the different materials encountered along the path:
std::vector<MaterialEntry> materialList;
// Store the hits on each sensor volume
std::vector<HitEntry> hitList;
for ( std::size_t i = 0; i < nIPSamples; ++i ) {
if ( nIPSamples > 1 ) {
TRandom3 random{ seed };
ip.SetZ( random.Gaus( meanIP, sigmaIP ) );
}
for ( std::size_t j = 0; j < nPointsPhi; ++j ) {
phi_deg = phiLow + ( j + 0.5 ) * phiRange / nPointsPhi;
phi_rad = phi_deg / 180. * TMath::Pi();
cosPhi = cos(phi_rad);
sinPhi = sin(phi_rad);
for ( std::size_t k = 0; k < nPointsEta; ++k ) {
eta = etaLow + ( k + 0.5 ) * etaRange / nPointsEta;
theta_rad = 2.0 * atan( exp( -eta ) );
//theta_deg = theta_rad * 180. / TMath::Pi();
cosTheta = cos(theta_rad);
sinTheta = sin(theta_rad);
ROOT::Math::XYZVector direction { sinTheta*cosPhi, sinTheta*sinPhi, cosTheta };
nav->SetCurrentPoint( ip.x(), ip.y(), ip.z() );
nav->SetCurrentDirection( direction.x(), direction.y(), direction.z() );
hitList.clear();
materialList.clear();
bool navFinished = false;
nav->FindNode();
while ( !navFinished ) {
const double* ipos { nav->GetCurrentPoint() };
const std::string nodeName { nav->GetCurrentNode()->GetName() };
const std::string matName { nav->GetCurrentNode()->GetVolume()->GetMaterial()->GetName() };
const std::string taggedMatName { tagMatName( matName, nodeName ) };
const double nodeRadLength { nav->GetCurrentNode()->GetVolume()->GetMaterial()->GetRadLen() };
const bool isSensor { taggedMatName.find( "_Sensor" ) != std::string::npos };
// TODO - foil scaling
// Verbose printout
std::ostringstream os;
os << ipos[0] << ", " << ipos[1] << ", " << ipos[2] << " " << nodeName << " - " << matName << " - " << taggedMatName << " with radlen = " << nodeRadLength;
if ( isSensor ) {
os << " - A SENSOR!";
}
if ( verbose ) { dd4hep::printout( dd4hep::INFO, "rad_scan", "%s", os.str().c_str() ); }
const ROOT::Math::XYZPoint initPos{ipos[0], ipos[1], ipos[2]};
nav->FindNextBoundaryAndStep();
const double* npos = nav->GetCurrentPoint();
const ROOT::Math::XYZPoint newPos{npos[0], npos[1], npos[2]};
const ROOT::Math::XYZPoint avgPos { ( ROOT::Math::XYZVector{initPos} + ROOT::Math::XYZVector{newPos} ) / 2.0 };
const double currentThick { ( newPos - initPos ).r() };
const double currentRadPercent { 100.0 * currentThick / nodeRadLength };
const std::size_t index { materialList.size() };
if ( isSensor ) {
hitList.push_back( { index, avgPos.rho(), avgPos.z() } );
}
materialList.push_back( { taggedMatName, index, currentThick, currentRadPercent } );
navFinished = ( ( newPos.z() > maxDistance ) || ( newPos.rho() > maxRadius ) );
}
// Fill hit information histograms
const std::size_t nHits { hitList.size() };
histHitsPerTrack->Fill( phi_deg, eta, nHits * ipSampleWeight );
if ( nHits > 0 ) {
// fill 1st hit info
histRadius1stHit->Fill( phi_deg, eta, hitList.front().rho * ipSampleWeight );
histZ1stHit->Fill( phi_deg, eta, hitList.front().z * ipSampleWeight );
// fill last hit info
histRadiusLastHit->Fill( phi_deg, eta, hitList.back().rho * ipSampleWeight );
histZLastHit->Fill( phi_deg, eta, hitList.back().z * ipSampleWeight );
if ( nHits > 1 ) {
// fill 2nd hit info
histRadius2ndHit->Fill( phi_deg, eta, hitList[1].rho * ipSampleWeight );
histZ2ndHit->Fill( phi_deg, eta, hitList[1].z * ipSampleWeight );
}
}
const std::size_t index1stHit { nHits > 0 ? hitList.front().index : -1 };
const std::size_t index2ndHit { nHits > 1 ? hitList[1].index : -1 };
const std::size_t indexLastHit { nHits > 0 ? hitList.back().index : -1 };
// Fill radiation length histograms
double totalRadLen{0.0};
std::map<std::string,double> perMaterialRadLen;
for ( auto& [ matName, matIndex, matThick, matRadLen ] : materialList ) {
if ( matIndex == index1stHit ) {
histBefore1stHitRadLen->Fill( phi_deg, eta, totalRadLen * ipSampleWeight );
fillPerMaterialHistograms( matHist_1stHit, histBefore1stHitRadLen, perMaterialRadLen, phi_deg, eta, ipSampleWeight );
}
if ( matIndex == index2ndHit ) {
histBefore2ndHitRadLen->Fill( phi_deg, eta, totalRadLen * ipSampleWeight );
fillPerMaterialHistograms( matHist_2ndHit, histBefore2ndHitRadLen, perMaterialRadLen, phi_deg, eta, ipSampleWeight );
}
if ( matIndex == indexLastHit ) {
histBeforeLastHitRadLen->Fill( phi_deg, eta, totalRadLen * ipSampleWeight );
fillPerMaterialHistograms( matHist_LastHit, histBeforeLastHitRadLen, perMaterialRadLen, phi_deg, eta, ipSampleWeight );
}
if ( perMaterialRadLen.find( matName ) == perMaterialRadLen.end() ) {
perMaterialRadLen[ matName ] = matRadLen;
} else {
perMaterialRadLen[ matName ] += matRadLen;
}
totalRadLen += matRadLen;
}
// Fill total radlen histograms - total and per-material
histRadLen->Fill( phi_deg, eta, totalRadLen * ipSampleWeight );
fillPerMaterialHistograms( matHist_Total, histRadLen, perMaterialRadLen, phi_deg, eta, ipSampleWeight );
}
}
}
TFile* out_file = new TFile( "rad_scan_etaphi_velo.root", "recreate" );
histRadLen->Write();
histBefore1stHitRadLen->Write();
histBefore2ndHitRadLen->Write();
histBeforeLastHitRadLen->Write();
histHitsPerTrack->Write();
histRadius1stHit->Write();
histRadius2ndHit->Write();
histRadiusLastHit->Write();
histZ1stHit->Write();
histZ2ndHit->Write();
histZLastHit->Write();
for ( auto& [ _, hist ] : matHist_Total ) { hist->Write(); }
for ( auto& [ _, hist ] : matHist_1stHit ) { hist->Write(); }
for ( auto& [ _, hist ] : matHist_2ndHit ) { hist->Write(); }
for ( auto& [ _, hist ] : matHist_LastHit ) { hist->Write(); }
out_file->Close();
std::cout << "Saved histogram to rad_scan_etaphi_velo.root" << std::endl;
return 1;
}
DECLARE_APPLY( rad_scan_etaphi_velo, tgeo_rad_scan_etaphi_velo )
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment