Merge branch 'master' of https://git.cern.ch/reps/AliRoot
[u/mrichter/AliRoot.git] / MUON / MUONChamberMaterialBudget.C
CommitLineData
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 40void 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.};
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};
69 Int_t iyMax[5] = {90, 120, 180, 250, 270};
4be12860 70
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],
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
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 (segmentationLevel*segmentationLevel*ixMax[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.000001*xIn, 1.000001*yIn, 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);
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[2*ist])) chId = 2*ist;
163 else if (currentNodePath.Contains(chamberName[2*ist+1])) chId = 2*ist+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[2*ist]->Fill(xIn,yIn,100.*localOneOverX0MeanCh[2*ist]);
177 hXOverX0[2*ist+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[2*ist] += localOneOverX0MeanCh[2*ist];
184 OneOverX0MeanCh[2*ist+1] += localOneOverX0MeanCh[2*ist+1];
185 totalLengthCh[2*ist] += localTotalLengthCh[2*ist];
186 totalLengthCh[2*ist+1] += localTotalLengthCh[2*ist+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=2*ist; i<=2*ist+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