MUON/MUONChamberMaterialBudget.C

   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
  40 void MUONChamberMaterialBudget(const char* geoFilename = "geometry.root", Int_t segmentationLevel = 1)
  41 {
  42   /// Draw the local chamber thickness over x0 (x/x0) used in the computation of Multiple Coulomb Scattering effets.
  43   /// Compute <x> and <x/x0> in a limited area (displayed on the histograms) to avoid edge effets.
  44   /// The resolution can be changed by changing the sementation level: resolution = 1 cm / segmentationLevel.
  45
  46   const char* chamberName[10] = {"SC01", "SC02", "SC03", "SC04", "SC05", "SC06", "SC07", "SC08", "SC09", "SC10"};
  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) {
  54     TGeoManager::Import(geoFilename);
  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
  62   Double_t zIn[5] =  {-510., -600.,  -800., -1150., -1350.};
  63   Double_t zOut[5] = {-600., -800., -1150., -1350., -1470.};
  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.};
  68   Int_t ixMax[5] = {90, 120, 165, 250, 260};
  69   Int_t iyMax[5] = {90, 120, 180, 250, 270};
  70
  71   // transverse area where to compute <x> and <x/x0> for each chamber
  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 };
  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),
  85                            segmentationLevel*ixMax[st], xIn0[st], xIn0[st]+ixMax[st],
  86                            segmentationLevel*iyMax[st], yIn0[st], yIn0[st]+iyMax[st]);
  87     hXOverX0[i]->SetOption("COLZ");
  88     hXOverX0[i]->SetStats(kFALSE);
  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)
 103     for (Int_t ix=0; ix<segmentationLevel*ixMax[ist]; ix++) {
 104
 105       Double_t xIn = xIn0[ist] + ((Double_t)ix+0.5) / ((Double_t)segmentationLevel);
 106
 107       // loop over position in bending direction (by step of 1cm)
 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;
 112
 113         Double_t yIn = yIn0[ist] + ((Double_t)iy+0.5) / ((Double_t)segmentationLevel);
 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]};
 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);
 126         if (!currentnode) break;
 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;
 134         do {
 135           // Get material properties
 136           TGeoMaterial *material = currentnode->GetVolume()->GetMedium()->GetMaterial();
 137           TString currentNodePath(gGeoManager->GetPath());
 138           x0 = material->GetRadLen();
 139
 140           // Get path length within this material
 141           gGeoManager->FindNextBoundary(remainingPathLength);
 142           localPathLength = gGeoManager->GetStep() + 1.e-6;
 143
 144           // Check if boundary within remaining path length.
 145           // If so, make sure to cross the boundary to prepare the next step
 146           if (localPathLength >= remainingPathLength) localPathLength = remainingPathLength;
 147           else {
 148             currentnode = gGeoManager->Step();
 149             if (!currentnode) { OK = kFALSE; break; }
 150             if (!gGeoManager->IsEntering()) {
 151               // make another small step to try to enter in new slice
 152               gGeoManager->SetStep(0.001);
 153               currentnode = gGeoManager->Step();
 154               if (!gGeoManager->IsEntering() || !currentnode) { OK = kFALSE; break; }
 155               localPathLength += 0.001;
 156             }
 157           }
 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;
 165
 166           // add current material budget
 167           localOneOverX0MeanCh[chId] += localPathLength / x0;
 168           localTotalLengthCh[chId] += localPathLength;
 169
 170         } while (remainingPathLength > TGeoShape::Tolerance());
 171
 172         // account for the local material characteristic if computed successfully
 173         if (OK) {
 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
 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]) {
 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     }
 194     cout<<"\rStation "<<ist+1<<": processing... 100%"<<endl;
 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