#endif
-void MUONChamberMaterialBudget(Bool_t warn = kFALSE)
+void MUONChamberMaterialBudget(const char* geoFilename = "geometry.root", Int_t segmentationLevel = 1)
{
- /// draw the local chamber thickness over x0 (x/x0) used in the computation of Multiple Coulomb Scattering effets.
+ /// Draw the local chamber thickness over x0 (x/x0) used in the computation of Multiple Coulomb Scattering effets.
/// Compute <x> and <x/x0> in a limited area (displayed on the histograms) to avoid edge effets.
+ /// The resolution can be changed by changing the sementation level: resolution = 1 cm / segmentationLevel.
+ const char* chamberName[10] = {"SC01", "SC02", "SC03", "SC04", "SC05", "SC06", "SC07", "SC08", "SC09", "SC10"};
Double_t OneOverX0MeanCh[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
Double_t totalLengthCh[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
Double_t OneOverX0Mean = 0.;
// Import TGeo geometry
if (!gGeoManager) {
- TGeoManager::Import("geometry.root");
+ TGeoManager::Import(geoFilename);
if (!gGeoManager) {
cout<<"getting geometry from file geometry.root failed"<<endl;
return;
}
// z intervals where to find the stations
- Double_t zIn[5] = {-510., -600., -850., -1250., -1350.};
- Double_t zOut[5] = {-600., -750., -1150., -1350., -1470.};
+ Double_t zIn[5] = {-510., -600., -800., -1150., -1350.};
+ Double_t zOut[5] = {-600., -800., -1150., -1350., -1470.};
// transverse area where to compute locally x and x/x0
Double_t xIn0[5] = {0., 0., 0., 0., 0.};
Double_t yIn0[5] = {0., 0., 0., 0., 0.};
- Int_t ixMax[5] = {90, 120, 180, 240, 270};
- Int_t iyMax[5] = {90, 120, 180, 240, 270};
+ Int_t ixMax[5] = {90, 120, 165, 250, 260};
+ Int_t iyMax[5] = {90, 120, 180, 250, 270};
// transverse area where to compute <x> and <x/x0> for each chamber
- Double_t xIn0ForMean[5] = { 5., 0., 40., 45., 50.};
- Double_t yIn0ForMean[5] = {21., 30., 0., 0., 0.};
- Int_t ixMaxForMean[5] = { 50, 70, 80, 115, 150 };
+ Double_t xIn0ForMean[5] = { 5., 5., 35., 40., 40.};
+ Double_t yIn0ForMean[5] = {20., 25., 0., 0., 0.};
+ Int_t ixMaxForMean[5] = { 50, 65, 85, 120, 160 };
Int_t iyMaxForMean[5] = { 60, 70, 110, 150, 150 };
// define output histograms
for (Int_t i=0; i<10; i++) {
Int_t st = i/2;
hXOverX0[i] = new TH2F(Form("hXOverX0_%d",i+1), Form("x/x0 on ch %d (%%)",i+1),
- ixMax[st], xIn0[st], xIn0[st]+ixMax[st],
- iyMax[st], yIn0[st], yIn0[st]+iyMax[st]);
+ segmentationLevel*ixMax[st], xIn0[st], xIn0[st]+ixMax[st],
+ segmentationLevel*iyMax[st], yIn0[st], yIn0[st]+iyMax[st]);
hXOverX0[i]->SetOption("COLZ");
+ hXOverX0[i]->SetStats(kFALSE);
bXOverX0[i] = new TBox(xIn0ForMean[st], yIn0ForMean[st],
xIn0ForMean[st]+ixMaxForMean[st], yIn0ForMean[st]+iyMaxForMean[st]);
bXOverX0[i]->SetLineStyle(2);
Int_t nPoints = 0;
// loop over position in non bending direction (by step of 1cm)
- for (Int_t ix=0; ix<ixMax[ist]; ix++) {
+ for (Int_t ix=0; ix<segmentationLevel*ixMax[ist]; ix++) {
- Double_t xIn = xIn0[ist] + 1.*ix;
+ Double_t xIn = xIn0[ist] + ((Double_t)ix+0.5) / ((Double_t)segmentationLevel);
// loop over position in bending direction (by step of 1cm)
- for (Int_t iy=0; iy<iyMax[ist]; iy++) {
+ for (Int_t iy=0; iy<segmentationLevel*iyMax[ist]; iy++) {
+ Int_t permilDone = 1000 * (ix * segmentationLevel*iyMax[ist] + iy + 1) /
+ (segmentationLevel*segmentationLevel*ixMax[ist]*iyMax[ist]);
+ if (permilDone%10 == 0) cout<<"\rStation "<<ist+1<<": processing... "<<permilDone/10<<"%"<<flush;
- Double_t yIn = yIn0[ist] + 1.*iy;
+ Double_t yIn = yIn0[ist] + ((Double_t)iy+0.5) / ((Double_t)segmentationLevel);
// Initialize starting point and direction
Double_t trackXYZIn[3] = {xIn, yIn, zIn[ist]};
Double_t trackXYZOut[3] = {1.000001*xIn, 1.000001*yIn, zOut[ist]};
- // Double_t trackXYZIn[3] = {35., -45., -510.};
- // Double_t trackXYZOut[3] = {100., -138., -1500.};
Double_t pathLength = TMath::Sqrt((trackXYZOut[0] - trackXYZIn[0])*(trackXYZOut[0] - trackXYZIn[0])+
(trackXYZOut[1] - trackXYZIn[1])*(trackXYZOut[1] - trackXYZIn[1])+
(trackXYZOut[2] - trackXYZIn[2])*(trackXYZOut[2] - trackXYZIn[2]));
b[1] = (trackXYZOut[1] - trackXYZIn[1]) / pathLength;
b[2] = (trackXYZOut[2] - trackXYZIn[2]) / pathLength;
TGeoNode *currentnode = gGeoManager->InitTrack(trackXYZIn, b);
- if (!currentnode) {
- if (warn) cout<<"start point out of geometry"<<endl;
- break;
- }
+ if (!currentnode) break;
Bool_t OK = kTRUE;
Double_t x0 = 0.; // radiation-length (cm-1)
Double_t localTotalLengthCh[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
Double_t localPathLength = 0.;
Double_t remainingPathLength = pathLength;
- Int_t chId = 2*ist;
- Bool_t enterCh = kFALSE;
do {
// Get material properties
TGeoMaterial *material = currentnode->GetVolume()->GetMedium()->GetMaterial();
- TString mName(material->GetName());
- mName.ToUpper();
+ TString currentNodePath(gGeoManager->GetPath());
x0 = material->GetRadLen();
// Get path length within this material
gGeoManager->FindNextBoundary(remainingPathLength);
localPathLength = gGeoManager->GetStep() + 1.e-6;
- // Check if boundary within remaining path length. If so, make sure to cross the boundary to prepare the next step
+
+ // Check if boundary within remaining path length.
+ // If so, make sure to cross the boundary to prepare the next step
if (localPathLength >= remainingPathLength) localPathLength = remainingPathLength;
else {
currentnode = gGeoManager->Step();
- if (!currentnode) {
- if (warn) cout<<"navigation failed(1)"<<endl;
- OK = kFALSE;
- break;
- }
+ if (!currentnode) { OK = kFALSE; break; }
if (!gGeoManager->IsEntering()) {
- // make another small step to try to enter in new absorber slice
+ // make another small step to try to enter in new slice
gGeoManager->SetStep(0.001);
currentnode = gGeoManager->Step();
- if (!gGeoManager->IsEntering() || !currentnode) {
- if (warn) cout<<"navigation failed(2)"<<endl;
- OK = kFALSE;
- break;
- }
+ if (!gGeoManager->IsEntering() || !currentnode) { OK = kFALSE; break; }
localPathLength += 0.001;
}
}
+ remainingPathLength -= localPathLength;
+
+ // check if entering a chamber of the current station or go to next step
+ Int_t chId;
+ if (currentNodePath.Contains(chamberName[2*ist])) chId = 2*ist;
+ else if (currentNodePath.Contains(chamberName[2*ist+1])) chId = 2*ist+1;
+ else continue;
- // check if entering a new chamber
- if (localPathLength > 10.) {
- if(!mName.Contains("AIR")) {
- if (warn) {
- material->Print();
- cout<<"changing chamber while crossing material"<<endl;
- }
- OK = kFALSE;
- break;
- }
- if (enterCh) chId++;
- //if (chId > 2*ist+1) break;
- enterCh = kFALSE;
- remainingPathLength -= localPathLength;
- continue;
- } else {
- enterCh = kTRUE;
- }
-
- // print parameters
- if(!mName.Contains("AIR")) {
- if(mName.Contains("DIPO")) {
- if (warn) cout<<"the track cross the dipole"<<endl;
- OK = kFALSE;
- break;
- }
- //material->Print();
- //cout<<"localPathLength = "<<localPathLength<<endl;
- if (chId <= 2*ist+1) {
- localOneOverX0MeanCh[chId] += localPathLength / x0;
- localTotalLengthCh[chId] += localPathLength;
- }
- }
+ // add current material budget
+ localOneOverX0MeanCh[chId] += localPathLength / x0;
+ localTotalLengthCh[chId] += localPathLength;
- // prepare next step
- remainingPathLength -= localPathLength;
} while (remainingPathLength > TGeoShape::Tolerance());
// account for the local material characteristic if computed successfully
- if (OK && (chId == 2*ist+2)) {
+ if (OK) {
// fill histograms in the full space
hXOverX0[2*ist]->Fill(xIn,yIn,100.*localOneOverX0MeanCh[2*ist]);
hXOverX0[2*ist+1]->Fill(xIn,yIn,100.*localOneOverX0MeanCh[2*ist+1]);
- // limit the chamber region for the computation of <x> and <x/x0>
- if (xIn > xIn0ForMean[ist]-0.5 && xIn < xIn0ForMean[ist]+1.*ixMaxForMean[ist]-0.5 &&
- yIn > yIn0ForMean[ist]-0.5 && yIn < yIn0ForMean[ist]+1.*iyMaxForMean[ist]-0.5) {
+ // computation of <x> and <x/x0> in a limited chamber region
+ if (xIn > xIn0ForMean[ist] && xIn < xIn0ForMean[ist]+1.*ixMaxForMean[ist] &&
+ yIn > yIn0ForMean[ist] && yIn < yIn0ForMean[ist]+1.*iyMaxForMean[ist]) {
nPoints++;
OneOverX0MeanCh[2*ist] += localOneOverX0MeanCh[2*ist];
OneOverX0MeanCh[2*ist+1] += localOneOverX0MeanCh[2*ist+1];
}
}
+ cout<<"\rStation "<<ist+1<<": processing... 100%"<<endl;
// normalize <x> and <x/x0> to the number of data points
for (Int_t i=2*ist; i<=2*ist+1; i++) {
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff