[u/mrichter/AliRoot.git] / MUON / MUONChamberMaterialBudget.C

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

// $Id$

/// \ingroup macros
/// \file MUONChamberMaterialBudget.C
/// \brief Utility macro to check tracking chamber material properties. 
///
/// \author Philippe Pillot, Subatech, Oct. 2009  

#if !defined(__CINT__) || defined(__MAKECINT__)

#include <Riostream.h>
#include <TH2F.h>
#include <TFile.h>
#include <TStyle.h>
#include <TBox.h>
#include <TString.h>
#include <TMath.h>
#include <TGeoManager.h>
#include <TGeoNode.h>
#include <TGeoMaterial.h>
#include <TGeoShape.h>

#endif

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.
  /// 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.;
  Double_t totalLength = 0.;
  
  // Import TGeo geometry
  if (!gGeoManager) {
    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.,  -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, 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.,  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
  gStyle->SetPalette(1);
  TFile *f = TFile::Open("MaterialBudget.root","RECREATE");
  TH2F* hXOverX0[10];
  TBox* bXOverX0[10];
  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),
			   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);
    bXOverX0[i]->SetLineWidth(2);
    bXOverX0[i]->SetFillStyle(0);
    hXOverX0[i]->GetListOfFunctions()->Add(bXOverX0[i]);
  }
  
  // loop over stations
  for (Int_t ist=0; ist<5; ist++) {
    
    Int_t nPoints = 0;
    
    // loop over position in non bending direction (by step of 1cm)
    for (Int_t ix=0; ix<segmentationLevel*ixMax[ist]; 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<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] + ((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 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]));
	Double_t b[3];
	b[0] = (trackXYZOut[0] - trackXYZIn[0]) / pathLength;
	b[1] = (trackXYZOut[1] - trackXYZIn[1]) / pathLength;
	b[2] = (trackXYZOut[2] - trackXYZIn[2]) / pathLength;
	TGeoNode *currentnode = gGeoManager->InitTrack(trackXYZIn, b);
	if (!currentnode) break;
	
	Bool_t OK = kTRUE;
	Double_t x0 = 0.;  // radiation-length (cm-1)
	Double_t localOneOverX0MeanCh[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
	Double_t localTotalLengthCh[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
	Double_t localPathLength = 0.;
	Double_t remainingPathLength = pathLength;
	do {
	  // Get material properties
	  TGeoMaterial *material = currentnode->GetVolume()->GetMedium()->GetMaterial();
	  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
	  if (localPathLength >= remainingPathLength) localPathLength = remainingPathLength;
	  else {
	    currentnode = gGeoManager->Step();
	    if (!currentnode) { OK = kFALSE; break; }
	    if (!gGeoManager->IsEntering()) {
	      // make another small step to try to enter in new slice
	      gGeoManager->SetStep(0.001);
	      currentnode = gGeoManager->Step();
	      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;
	  
	  // add current material budget
	  localOneOverX0MeanCh[chId] += localPathLength / x0;
	  localTotalLengthCh[chId] += localPathLength;
	  
	} while (remainingPathLength > TGeoShape::Tolerance());
	
	// account for the local material characteristic if computed successfully
	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]);
	  
	  // 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];
	    totalLengthCh[2*ist] += localTotalLengthCh[2*ist];
	    totalLengthCh[2*ist+1] += localTotalLengthCh[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++) {
      OneOverX0MeanCh[i] /= (Double_t) nPoints;
      totalLengthCh[i] /= (Double_t) nPoints;
    }
    
  }
  
  // print results
  cout<<endl<<endl;
  cout<<"chamber   thickness (cm)    x/x0 (%)"<<endl;
  cout<<"------------------------------------"<<endl;
  for (Int_t i=0; i<10; i++) {
    printf("  %2d          %4.2f            %4.2f\n",i+1,totalLengthCh[i],100.*OneOverX0MeanCh[i]);
    totalLength += totalLengthCh[i];
    OneOverX0Mean += OneOverX0MeanCh[i];
  }
  cout<<"------------------------------------"<<endl;
  printf("  tot         %4.1f            %4.1f\n",totalLength,100.*OneOverX0Mean);
  cout<<endl;
  
  // save histograms
  f->Write();
  f->Close();
  
}
Commit	Line	Data
4be12860	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	// $Id$
	17
	18	/// \ingroup macros
	19	/// \file MUONChamberMaterialBudget.C
	20	/// \brief Utility macro to check tracking chamber material properties.
	21	///
	22	/// \author Philippe Pillot, Subatech, Oct. 2009
	23
	24	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	25
	26	#include <Riostream.h>
	27	#include <TH2F.h>
	28	#include <TFile.h>
	29	#include <TStyle.h>
	30	#include <TBox.h>
	31	#include <TString.h>
	32	#include <TMath.h>
	33	#include <TGeoManager.h>
	34	#include <TGeoNode.h>
	35	#include <TGeoMaterial.h>
	36	#include <TGeoShape.h>
	37
	38	#endif
	39
fb2d9d49	40	void MUONChamberMaterialBudget(const char* geoFilename = "geometry.root", Int_t segmentationLevel = 1)
4be12860	41	{
fb2d9d49	42	/// Draw the local chamber thickness over x0 (x/x0) used in the computation of Multiple Coulomb Scattering effets.
4be12860	43	/// Compute <x> and <x/x0> in a limited area (displayed on the histograms) to avoid edge effets.
fb2d9d49	44	/// The resolution can be changed by changing the sementation level: resolution = 1 cm / segmentationLevel.
4be12860	45
fb2d9d49	46	const char* chamberName[10] = {"SC01", "SC02", "SC03", "SC04", "SC05", "SC06", "SC07", "SC08", "SC09", "SC10"};
4be12860	47	Double_t OneOverX0MeanCh[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
	48	Double_t totalLengthCh[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
	49	Double_t OneOverX0Mean = 0.;
	50	Double_t totalLength = 0.;
	51
	52	// Import TGeo geometry
	53	if (!gGeoManager) {
fb2d9d49	54	TGeoManager::Import(geoFilename);
4be12860	55	if (!gGeoManager) {
	56	cout<<"getting geometry from file geometry.root failed"<<endl;
	57	return;
	58	}
	59	}
	60
	61	// z intervals where to find the stations
fb2d9d49	62	Double_t zIn[5] = {-510., -600., -800., -1150., -1350.};
fb2d9d49	63	Double_t zOut[5] = {-600., -800., -1150., -1350., -1470.};
4be12860	64
	65	// transverse area where to compute locally x and x/x0
	66	Double_t xIn0[5] = {0., 0., 0., 0., 0.};
	67	Double_t yIn0[5] = {0., 0., 0., 0., 0.};
fb2d9d49	68	Int_t ixMax[5] = {90, 120, 165, 250, 260};
fb2d9d49	69	Int_t iyMax[5] = {90, 120, 180, 250, 270};
4be12860	70
4be12860	71	// transverse area where to compute <x> and <x/x0> for each chamber
fb2d9d49	72	Double_t xIn0ForMean[5] = { 5., 5., 35., 40., 40.};
	73	Double_t yIn0ForMean[5] = {20., 25., 0., 0., 0.};
	74	Int_t ixMaxForMean[5] = { 50, 65, 85, 120, 160 };
4be12860	75	Int_t iyMaxForMean[5] = { 60, 70, 110, 150, 150 };
	76
	77	// define output histograms
	78	gStyle->SetPalette(1);
	79	TFile *f = TFile::Open("MaterialBudget.root","RECREATE");
	80	TH2F* hXOverX0[10];
	81	TBox* bXOverX0[10];
	82	for (Int_t i=0; i<10; i++) {
	83	Int_t st = i/2;
	84	hXOverX0[i] = new TH2F(Form("hXOverX0_%d",i+1), Form("x/x0 on ch %d (%%)",i+1),
fb2d9d49	85	segmentationLevel*ixMax[st], xIn0[st], xIn0[st]+ixMax[st],
fb2d9d49	86	segmentationLevel*iyMax[st], yIn0[st], yIn0[st]+iyMax[st]);
4be12860	87	hXOverX0[i]->SetOption("COLZ");
fb2d9d49	88	hXOverX0[i]->SetStats(kFALSE);
4be12860	89	bXOverX0[i] = new TBox(xIn0ForMean[st], yIn0ForMean[st],
	90	xIn0ForMean[st]+ixMaxForMean[st], yIn0ForMean[st]+iyMaxForMean[st]);
	91	bXOverX0[i]->SetLineStyle(2);
	92	bXOverX0[i]->SetLineWidth(2);
	93	bXOverX0[i]->SetFillStyle(0);
	94	hXOverX0[i]->GetListOfFunctions()->Add(bXOverX0[i]);
	95	}
	96
	97	// loop over stations
	98	for (Int_t ist=0; ist<5; ist++) {
	99
	100	Int_t nPoints = 0;
	101
	102	// loop over position in non bending direction (by step of 1cm)
fb2d9d49	103	for (Int_t ix=0; ix<segmentationLevel*ixMax[ist]; ix++) {
4be12860	104
fb2d9d49	105	Double_t xIn = xIn0[ist] + ((Double_t)ix+0.5) / ((Double_t)segmentationLevel);
4be12860	106
4be12860	107	// loop over position in bending direction (by step of 1cm)
fb2d9d49	108	for (Int_t iy=0; iy<segmentationLevel*iyMax[ist]; iy++) {
	109	Int_t permilDone = 1000 * (ix * segmentationLevel*iyMax[ist] + iy + 1) /
	110	(segmentationLevelsegmentationLevelixMax[ist]*iyMax[ist]);
	111	if (permilDone%10 == 0) cout<<"\rStation "<<ist+1<<": processing... "<<permilDone/10<<"%"<<flush;
4be12860	112
fb2d9d49	113	Double_t yIn = yIn0[ist] + ((Double_t)iy+0.5) / ((Double_t)segmentationLevel);
4be12860	114
	115	// Initialize starting point and direction
	116	Double_t trackXYZIn[3] = {xIn, yIn, zIn[ist]};
	117	Double_t trackXYZOut[3] = {1.000001xIn, 1.000001yIn, zOut[ist]};
4be12860	118	Double_t pathLength = TMath::Sqrt((trackXYZOut[0] - trackXYZIn[0])*(trackXYZOut[0] - trackXYZIn[0])+
	119	(trackXYZOut[1] - trackXYZIn[1])*(trackXYZOut[1] - trackXYZIn[1])+
	120	(trackXYZOut[2] - trackXYZIn[2])*(trackXYZOut[2] - trackXYZIn[2]));
	121	Double_t b[3];
	122	b[0] = (trackXYZOut[0] - trackXYZIn[0]) / pathLength;
	123	b[1] = (trackXYZOut[1] - trackXYZIn[1]) / pathLength;
	124	b[2] = (trackXYZOut[2] - trackXYZIn[2]) / pathLength;
	125	TGeoNode *currentnode = gGeoManager->InitTrack(trackXYZIn, b);
fb2d9d49	126	if (!currentnode) break;
4be12860	127
	128	Bool_t OK = kTRUE;
	129	Double_t x0 = 0.; // radiation-length (cm-1)
	130	Double_t localOneOverX0MeanCh[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
	131	Double_t localTotalLengthCh[10] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.};
	132	Double_t localPathLength = 0.;
	133	Double_t remainingPathLength = pathLength;
4be12860	134	do {
	135	// Get material properties
	136	TGeoMaterial *material = currentnode->GetVolume()->GetMedium()->GetMaterial();
fb2d9d49	137	TString currentNodePath(gGeoManager->GetPath());
4be12860	138	x0 = material->GetRadLen();
	139
	140	// Get path length within this material
	141	gGeoManager->FindNextBoundary(remainingPathLength);
	142	localPathLength = gGeoManager->GetStep() + 1.e-6;
fb2d9d49	143
	144	// Check if boundary within remaining path length.
	145	// If so, make sure to cross the boundary to prepare the next step
4be12860	146	if (localPathLength >= remainingPathLength) localPathLength = remainingPathLength;
	147	else {
	148	currentnode = gGeoManager->Step();
fb2d9d49	149	if (!currentnode) { OK = kFALSE; break; }
4be12860	150	if (!gGeoManager->IsEntering()) {
fb2d9d49	151	// make another small step to try to enter in new slice
4be12860	152	gGeoManager->SetStep(0.001);
4be12860	153	currentnode = gGeoManager->Step();
fb2d9d49	154	if (!gGeoManager->IsEntering() \|\| !currentnode) { OK = kFALSE; break; }
4be12860	155	localPathLength += 0.001;
	156	}
	157	}
fb2d9d49	158	remainingPathLength -= localPathLength;
	159
	160	// check if entering a chamber of the current station or go to next step
	161	Int_t chId;
	162	if (currentNodePath.Contains(chamberName[2ist])) chId = 2ist;
	163	else if (currentNodePath.Contains(chamberName[2ist+1])) chId = 2ist+1;
	164	else continue;
4be12860	165
fb2d9d49	166	// add current material budget
	167	localOneOverX0MeanCh[chId] += localPathLength / x0;
	168	localTotalLengthCh[chId] += localPathLength;
4be12860	169
4be12860	170	} while (remainingPathLength > TGeoShape::Tolerance());
	171
	172	// account for the local material characteristic if computed successfully
fb2d9d49	173	if (OK) {
4be12860	174
	175	// fill histograms in the full space
	176	hXOverX0[2ist]->Fill(xIn,yIn,100.localOneOverX0MeanCh[2*ist]);
	177	hXOverX0[2ist+1]->Fill(xIn,yIn,100.localOneOverX0MeanCh[2*ist+1]);
	178
fb2d9d49	179	// computation of <x> and <x/x0> in a limited chamber region
	180	if (xIn > xIn0ForMean[ist] && xIn < xIn0ForMean[ist]+1.*ixMaxForMean[ist] &&
	181	yIn > yIn0ForMean[ist] && yIn < yIn0ForMean[ist]+1.*iyMaxForMean[ist]) {
4be12860	182	nPoints++;
	183	OneOverX0MeanCh[2ist] += localOneOverX0MeanCh[2ist];
	184	OneOverX0MeanCh[2ist+1] += localOneOverX0MeanCh[2ist+1];
	185	totalLengthCh[2ist] += localTotalLengthCh[2ist];
	186	totalLengthCh[2ist+1] += localTotalLengthCh[2ist+1];
	187	}
	188
	189	}
	190
	191	}
	192
	193	}
fb2d9d49	194	cout<<"\rStation "<<ist+1<<": processing... 100%"<<endl;
4be12860	195
	196	// normalize <x> and <x/x0> to the number of data points
	197	for (Int_t i=2ist; i<=2ist+1; i++) {
	198	OneOverX0MeanCh[i] /= (Double_t) nPoints;
	199	totalLengthCh[i] /= (Double_t) nPoints;
	200	}
	201
	202	}
	203
	204	// print results
	205	cout<<endl<<endl;
	206	cout<<"chamber thickness (cm) x/x0 (%)"<<endl;
	207	cout<<"------------------------------------"<<endl;
	208	for (Int_t i=0; i<10; i++) {
	209	printf(" %2d %4.2f %4.2f\n",i+1,totalLengthCh[i],100.*OneOverX0MeanCh[i]);
	210	totalLength += totalLengthCh[i];
	211	OneOverX0Mean += OneOverX0MeanCh[i];
	212	}
	213	cout<<"------------------------------------"<<endl;
	214	printf(" tot %4.1f %4.1f\n",totalLength,100.*OneOverX0Mean);
	215	cout<<endl;
	216
	217	// save histograms
	218	f->Write();
	219	f->Close();
	220
	221	}
	222