In MUONChamberMaterialBudget.C
authorivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 6 Nov 2009 10:54:35 +0000 (10:54 +0000)
committerivana <ivana@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 6 Nov 2009 10:54:35 +0000 (10:54 +0000)
- Use the name of the mother volumes instead of the position to select
  material belonging to a given chamber
- The resolution can be changed by changing the sementation level:
  resolution = 1 cm / segmentationLevel
(Philippe P.)

MUON/MUONChamberMaterialBudget.C

index 060a71d..8a49d0a 100644 (file)
 
 #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.;
@@ -49,7 +51,7 @@ void MUONChamberMaterialBudget(Bool_t warn = kFALSE)
   
   // Import TGeo geometry
   if (!gGeoManager) {
-    TGeoManager::Import("geometry.root");
+    TGeoManager::Import(geoFilename);
     if (!gGeoManager) {
       cout<<"getting geometry from file geometry.root failed"<<endl;
       return;
@@ -57,19 +59,19 @@ void MUONChamberMaterialBudget(Bool_t warn = kFALSE)
   }
   
   // 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
@@ -80,9 +82,10 @@ void MUONChamberMaterialBudget(Bool_t warn = kFALSE)
   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);
@@ -97,20 +100,21 @@ void MUONChamberMaterialBudget(Bool_t warn = kFALSE)
     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]));
@@ -119,10 +123,7 @@ void MUONChamberMaterialBudget(Bool_t warn = kFALSE)
        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)
@@ -130,88 +131,54 @@ void MUONChamberMaterialBudget(Bool_t warn = kFALSE)
        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];
@@ -224,6 +191,7 @@ void MUONChamberMaterialBudget(Bool_t warn = kFALSE)
       }
       
     }
+    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++) {