Point5MaterialMap.cc

Go to the documentation of this file.

#ifndef Point5MaterialMap_cc
#define Point5MaterialMap_cc

#include <iostream>
#include <cmath>

#include "GeneratorInterface/CosmicMuonGenerator/interface/CosmicMuonParameters.h"



inline int inPlug(double vx, double vy, double vz,
          double PlugVx = PlugOnShaftVx, double PlugVz = PlugOnShaftVz) {
  if (vy > SurfaceOfEarth && vy < SurfaceOfEarth + PlugWidth) {
    if (vx > PlugVx - PlugXlength/2. && vx < PlugVx + PlugXlength/2. &&
    vz > PlugVz - PlugZlength/2. && vz < PlugVz + PlugZlength/2.) return Plug;
    if (vz >= PlugVz - PlugZlength/2. - PlugNoseZlength && vz < PlugVz - PlugZlength/2. &&
    vx > PlugVx - PlugNoseXlength/2. && vx < PlugVx + PlugNoseXlength/2.) return Plug;
  }
  return Unknown;
}


inline int inAirAfterPlug(double vx, double vy, double vz) {
  // particles above surface of earth
  if (vy >= SurfaceOfEarth) return Air;
  
  // CMS cavern (UXC 55)
  if (std::fabs(vz) < 26548. && sqrt((vx*1.1576)*(vx*1.1576) + vy*vy) < 15460. &&
      vy > -8762) return Air;
  
  // access shaft (PX 56)
  if (vy > 0. && vy < (SurfaceOfEarth-2250.) && 
      sqrt(vx*vx + (vz-Z_PX56)*(vz-Z_PX56)) < 10250.) return Air;
  
  //surface hall ground floor
  if (vy >= SurfaceOfEarth-2250. && vy < SurfaceOfEarth) {
    if (sqrt(vx*vx + (vz-Z_PX56)*(vz-Z_PX56)) < 10250.
    && vz-Z_PX56 > -7000. && vz-Z_PX56 < 7000.) return Air;
    if (vx > -2400. && vx < 2400. && vz-Z_PX56 >= -9800. && vz-Z_PX56 < -7000. )
      return Air;
  }
  
    // Shaft (PM 54)
  if (vy > 3233. && vy < (SurfaceOfEarth) && 
      sqrt((vx-26600.)*(vx-26600.) + (vz-30100.-Z_PX56)*(vz-30100.-Z_PX56)) < 6050.)
    //sqrt((vx-5000.)*(vx-5000.) + (vz-30100.-Z_PX56)*(vz-30100.-Z_PX56)) < 10050.)//@@@@@
    return Air;
  
  
  // Shaft (PM 56)
  if (vy > -3700. && vy < (SurfaceOfEarth) &&
      sqrt((vx-18220.)*(vx-18220.) + (vz+24227.-Z_PX56)*(vz+24227.-Z_PX56)) < 3550.)
      return Air;
  
  //Service cavern USC 55
  if (vz > -22050. && vz < 62050. &&
      sqrt((vx-29550.)*(vx-29550.) + (vy-3233.)*(vy-3233.)) < 9050. && 
      vy > -3650.) return Air;
  
  //UXC55 cavern endcap beam openings
  if (((vz >= -29278. && vz < -26548.) || (vz >= 26548 && vz <= 29278.)) &&
      sqrt(vy*vy + vx*vx) < 800.) return Air; 
  
  //Pillar access between CMS collision cavern and service cavern
  if (vx > 10000. && vx < 20500.+3050. && //TX54 galerie
      vz > 14460. && vz < 16260. &&//estimated wall thickness 1m
      vy > -8680. && vy < -1507.) return Air;
  
  if (vx > 10000. && vx < 16500. && //TX54 galerie 30.9 degree opening to UXC55
      vy > -8680. && vy < -1507. && //@vx=13300 delta(vz)=1865., 2vx=1650. delta(vz)=0
      vz > 14460.-1865.*(16500.-vx)/(3200.) && vz <= 14460.) return Air;
  
  if (vx > 13300. && vx < 20500.+3050. && //TX54 galerie
      vz > 14460. && vz < 16260. && //estimated wall thickness 1m
      vy > -8680. && vy < -1507.) return Air;
  
  if (vx > 26600-6050. && vx < 20500.+3050. && //TX54 going up in sewrvice cavern
      vz > 14460. && vz < 16260. && //estimated wall thickness 1m
      vy >= -1507. && vy < 1000.) return Air;
  
  //R56, LHC beam tunnel East
  if (vz > -85000. && vz <= -29278. && //UJ57 junction to UXC55 cavern
      sqrt(vy*vy + (vx-350.)*(vx-350.)) < 1900. &&
      vy > -1000.) return Air;
  
  //R54, LHC beam tunnel West
  if (vz >= 29278. && vz < 63000. && //UJ57 junction to UXC55 cavern
      sqrt(vy*vy + (vx-350.)*(vx-350.)) < 1900. &&
      vy > -1000.) return Air;
  
  //UJ56 cavern
  if (vz > -58875. && vz < -33927. &&
      sqrt(vy*vy + (vx-4450.)*(vx-4450.)) < 6750. && vy > -1000. && //and beam shielding
      !(vx > 2250. && vx < 4250. && (vz > -(33927.+18000.) || vz < -(33927.+19500.))) && 
      !(vx >= 4250. && vx < 6650. && vz > -(33927.+18000.) && vz < -(33927.+16000.)))
    return Air;
  
  //connection between PM56 shaft and UJ56 cavern
  if (vx > 9000. && vx < 18220. &&
      sqrt((vy-50.)*(vy-50.)+(vz+24227.-Z_PX56)*(vz+24227.-Z_PX56)) < 3550.
      && vy > -1000.) return Air;
  
  return Unknown;
  
}


inline int inWallAfterAir(double vx, double vy, double vz) {
  // phase II surface building
  if (vy < SurfaceOfEarth && vy >= (SurfaceOfEarth-2250.)) {
    if (std::fabs(vz-Z_PX56) < 30000. && std::fabs(vx) < 10950) return Wall;
    // foundation of crane
    if (std::fabs(vz-Z_PX56) < 9000. && std::fabs(vx) >= 10950 && std::fabs(vx) < 16950) 
      return Wall;
  }
  
  // CMS cavern (UXC 55)
  if (std::fabs(vz) < 29278. && sqrt((vx*1.1576)*(vx*1.1576) + vy*vy) < 16830. &&
      vy > -11762.) return Wall;
  
  // access shaft (PX 56)
  if (vy > 0. && vy < (SurfaceOfEarth-2250.) && //t(shaft wall)=2150.
      sqrt(vx*vx + (vz-Z_PX56)*(vz-Z_PX56)) < 12400.) return Wall;
  
  // Shaft (PM 54)
  if (vy > 3233. && vy < (SurfaceOfEarth-1000.) && //t~=t(PX56)/R(PX56)*R(PM54)
      sqrt((vx-26600.)*(vx-26600.) + (vz-30100.-Z_PX56)*(vz-30100.-Z_PX56)) 
      < 6050.+2150./10250.*6050.) return Wall;
  //sqrt((vx-5000.)*(vx-5000.) + (vz-30100.-Z_PX56)*(vz-30100.-Z_PX56))//@@@@@ 
  //< 10050.+2150./10250.*6050.) return Wall;//@@@@@
  else if (vy >= SurfaceOfEarth-1000. && vy < SurfaceOfEarth && //t~=t(PX56)/R(PX56)*R(PM54)
       sqrt((vx-26600.)*(vx-26600.) + (vz-30100.-Z_PX56)*(vz-30100.-Z_PX56)) 
       < 6050.+2150./10250.*6050. +1800.) return Wall;
  //sqrt((vx-5000.)*(vx-5000.) + (vz-30100.-Z_PX56)*(vz-30100.-Z_PX56))//@@@@@ 
  //< 10050.+2150./10250.*10050. +1800.) return Wall; //@@@@@@
  
  // Shaft (PM 56)
  if (vy > -5450. && vy < (SurfaceOfEarth-1000.) && //t~=t(PX56)/R(PX56)*R(PM56)
      sqrt((vx-18220.)*(vx-18220.) + (vz+24227.-Z_PX56)*(vz+24227.-Z_PX56)) 
      < 3550.+2150./10250.*3550.) return Wall;
  else if (vy > SurfaceOfEarth-1000. && vy < SurfaceOfEarth && //t~=t(PX56)/R(PX56)*R(PM56)
       sqrt((vx-18220.)*(vx-18220.) + (vz+24227.-Z_PX56)*(vz+24227.-Z_PX56)) 
       < 3550.+2150./10250.*3550. +1800.) return Wall;
  
  //Service cavern USC 55
  if (vz > -(22050.+1150.) && vz < (62050.+1150.) &&
      sqrt((vx-29550.)*(vx-29550.) + (vy-3233.)*(vy-3233.)) < 9050.+950. &&
      vy > -3650.-2000.) return Wall; //8762=estimate, to be checked 
  
  //Pillar between CMS collision cavern and service cavern
  if (vz > -29278.+1000. && vz < 29278.+1000.) {
    if (vy > -17985. && vy < 10410. && vx > 13300. && vx < 20500.) 
      return Wall;
    
    if (vy > 0. && vy < 10410. && vx > 10000. && vx <= 13300.)
      return Wall;
    
    if (vy > -3650.-2000. && vy < -3233. && // bottom edge between pillar and service cavern
    vx > 20000. && vx < 24000.) return Wall;
    if (vy > -11762. && vy < -5000. && // bottom edge between pillar and UXC55 cavern
    vx > 10500. && vx < 14000.) return Wall;
  }
  
  if (vy > -14000. && vy < -1450. && //TX54 galerie surrounding
      vz > 13460. && vz < 17260. && vx >= 20500. && vx < 24550.) 
    return Wall;
  
  //R56, LHC beam tunnel East
  if (vz > -85000. && vz < -28510.) { //UJ57 junction to UXC55 cavern
    if (sqrt(vy*vy + (vx-350.)*(vx-350.)) < 2250.) return Wall;
  }
  
  //R54, LHC beam tunnel West
  if (vz > 26550. && vz < 63000. && //UJ57 junction to UXC55 cavern
      sqrt(vy*vy + (vx-350.)*(vx-350.)) < 2250.) return Wall;
  
  //UJ56 cavern
  if (vz > -(58875.+500.) && vz < -(33927.-500.) &&
      sqrt(vy*vy + (vx-4450.)*(vx-4450.)) < (6750.+500.) && vy > -3650.) 
    return Wall;
  
  //connection between PM56 shaft and UJ56 cavern
  if (vx > 9000. && vx < 18220. &&
      sqrt((vy-50.)*(vy-50.)+(vz+24227.-Z_PX56)*(vz+24227.-Z_PX56)) < 3550.+500. 
      && vy > -3650.) return Wall;
  
  return Unknown;
}



inline int inClayOrRockAfterWall(double vx, double vy, double vz, double ClayWidth) {
  
  //So, it is not plug, air and wall, Check for clay
  if (vy >= SurfaceOfEarth - ClayWidth && vy < SurfaceOfEarth)
    return Clay;
  
  //So, it is not plug, air, wall and clay, Check for rock
  if (vy < SurfaceOfEarth - ClayWidth)
    return Rock;
  
  return Unknown;
  
}



inline int inClayAfterWall(double vx, double vy, double vz, double ClayWidth) {
  
  //So, it is not plug, air and wall, Check for clay
  if (vy >= SurfaceOfEarth - ClayWidth && vy < SurfaceOfEarth)
    return Clay;
  
  return Unknown;
  
}



inline int inRockAfterClay(double vx, double vy, double vz, double ClayWidth) {
  
  //So, it is not plug, air, wall and clay, Check for rock
  if (vy < SurfaceOfEarth - ClayWidth)
    return Rock;
  
  return Unknown;
  
}





inline int inMat(double vx, double vy, double vz,
         double PlugVx = PlugOnShaftVx, double PlugVz = PlugOnShaftVz,
         double ClayWidth = DefaultClayWidth) {
  
  //check for Plug
  if (inPlug(vx, vy, vz, PlugVx, PlugVz)) return Plug;
  
  //So, it is not plug, Check for air
  if (inAirAfterPlug(vx, vy, vz)) return Air;
  
  //So, it is not plug and air, Check for wall
  if (inWallAfterAir(vx, vy, vz)) return Wall;
  
  //So, it is not plug, air and wall, Check for clay
  if (vy >= SurfaceOfEarth - ClayWidth && vy < SurfaceOfEarth)
    return Clay;
  
  //So, it is not plug, air, wall and clay, Check for rock
  if (vy < SurfaceOfEarth - ClayWidth)
    return Rock;
  
  
  std::cout << "Point5MaterialMap.h: Warning! No Material recognised for point: vx="
        << vx << " vy=" << vy << " vz=" << vz << std::endl;
  //Something went wrong
  return Unknown;
  
}




#endif