]>
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.}; |
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 |