// ----------------------------------------------------------------------------
// $Id: LCDGetParameters.cxx,v 1.2 2001/06/23 16:14:14 toshi Exp $
// ----------------------------------------------------------------------------
//
// $Log: LCDGetParameters.cxx,v $
// Revision 1.2  2001/06/23 16:14:14  toshi
// Char_t* -> const Char_t*
//
// Revision 1.1  2001/05/10 20:39:55  toshi
// C++ file name convention has been changed from *.C to *.cxx .
//
//
#include "LCDGetParameters.h"
#include <stdlib.h>

//______________________________________________________________________
// LCDGetParameters
//
ClassImp(LCDGetParameters)
  
  LCDGetParameters::LCDGetParameters() {
  m_Field=0;

  m_ipSigmaX=0;
  m_ipSigmaY=0;
  m_ipSigmaZ=0;

  m_minHit=0;
  m_PtMin=0;
  m_PolarInner=0;
  m_PolarOuter=0;

  m_CalEnergyEMA=0;
  m_CalEnergyEMB=0;
  m_CalEnergyHADA=0;
  m_CalEnergyHADB=0;
  m_CalTransEMA=0;
  m_CalTransEMB=0;
  m_CalTransHADA=0;
  m_CalTransHADB=0;

  m_EMEnergyScale=1;
  m_HADEnergyScale=1;
  m_MUEnergyScale=1;
  m_LUMEnergyScale=1;

  m_EMThetaSeg=0;
  m_HADThetaSeg=0;
  m_MUThetaSeg=0;
  m_LUMThetaSeg=0;

  m_EMphiSeg=0;
  m_HADphiSeg=0;
  m_MUphiSeg=0;
  m_LUMphiSeg=0;

  m_EMBarrelInnerR=0;
  m_HADBarrelInnerR=0;
  m_MUBarrelInnerR=0;
  m_EMBarrelOuterR=0;
  m_HADBarrelOuterR=0;
  m_MUBarrelOuterR=0;
  m_EMBarrelOuterZ=0;
  m_HADBarrelOuterZ=0;
  m_MUBarrelOuterZ=0;
  m_EMBarrelLayer=0;
  m_HADBarrelLayer=0;
  m_MUBarrelLayer=0;
  m_EMBarrelThick=0;
  m_HADBarrelThick=0;
  m_MUBarrelThick=0;
  m_EMBarrelAbsThick=0;
  m_HADBarrelAbsThick=0;
  m_MUBarrelAbsThick=0;

  m_EMEndcapInnerZ=0;
  m_HADEndcapInnerZ=0;
  m_MUEndcapInnerZ=0;
  m_EMEndcapLayer=0;
  m_HADEndcapLayer=0;
  m_MUEndcapLayer=0;
  m_EMEndcapThick=0;
  m_HADEndcapThick=0;
  m_MUEndcapThick=0;

  m_CoilInnerR=0;
  m_CoilOuterZ=0;
  m_CoilThick=0;
  m_CoilIntLeng=0;
  m_PipeOuter=0;
  
  m_DetectorName="";
  m_CoilPlace="";
  m_EMBarrelAbsMat="";
  m_HADBarrelAbsMat="";
  m_MUBarrelAbsMat="";

  m_Volumes=0;
  m_InnerVolume=0;
  m_nVolumes=0;
}

 LCDGetParameters::LCDGetParameters(FILE* parFile):
m_Volumes(0),m_InnerVolume(0),m_nVolumes(0)
{
  SetDetectorGeomFile(parFile);
}

 void LCDGetParameters::SetDetectorGeomFile(FILE* parFile) {
  
  TString input_line, input_string;
  
  Char_t temp[100];
  Char_t name[100];
  Char_t equals[100];
  
  Double_t read_value = 0.;
  
  Int_t pos;
  
  // rewind the file
  fseek(parFile, 0L, 0);
  
  Int_t inpos; 
  while ( (pos = fscanf(parFile,"%s",&temp[0])) > 0 ) {    
    input_line = temp;    
    if(input_line.Contains("Detector:")) {      
      // form of parameters is: name = read_value      
      while( (inpos = fscanf(parFile,"%s",&name[0]) ) > 0) {	
	input_string = name; 	
	if ( input_string.Contains("end")) break;
	inpos = fscanf(parFile,"%s %s",&equals[0],&temp[0]);	
	if ( input_string.Contains("DetectorName" )) {
	  m_DetectorName = temp;
	}
      }
    }
        
    if(input_line.Contains("IP:")) {      
      // form of parameters is: name = read_value      
      while( (inpos =fscanf(parFile,"%s",&name[0]) ) > 0){	
	input_string = name; 
	if ( input_string.Contains("end")) break;
	inpos = fscanf(parFile,"%s %s",&equals[0],&temp[0]);
	read_value = strtod(temp,(Char_t**)NULL);	
	
	if ( input_string.Contains("IPsigmaX" )) {
	  m_ipSigmaX = read_value;
	} else if ( input_string.Contains("IPsigmaY" )) {
	  m_ipSigmaY = read_value;
	} else if ( input_string.Contains("IPsigmaZ" )) {
	  m_ipSigmaZ = read_value;
	}
      }
    }
        
    if(input_line.Contains("Tracking:")) {      
      // form of parameters is: name = read_value      
      while( (inpos =fscanf(parFile,"%s",&name[0]) ) > 0){	
	input_string = name; 	
	if ( input_string.Contains("end")) break;
	inpos = fscanf(parFile,"%s %s",&equals[0],&temp[0]);
	read_value = strtod(temp,(Char_t**)NULL);
	
	if ( input_string.Contains("minHit" )) {
	  m_minHit = (Int_t) read_value;
	} else if ( input_string.Contains("PtMin" )) {
	  m_PtMin = read_value;
	} else if ( input_string.Contains("PolarInner" )) {
	  m_PolarInner = read_value;
	} else if ( input_string.Contains("PolarOuter" )) {
	  m_PolarOuter = read_value;
	} else if ( input_string.Contains("Magnetic_Field" )) {
	  m_Field = read_value;
	}
      }
    }

    if(input_line.Contains("CalSmearParam:")) {      
      // form of parameters is: name = read_value      
      while( (inpos =fscanf(parFile,"%s",&name[0]) ) > 0){	
	input_string = name; 
	if ( input_string.Contains("end")) break;
	inpos = fscanf(parFile,"%s %s",&equals[0],&temp[0]);
	read_value = strtod(temp,(Char_t**)NULL);	
	
	if ( input_string.Contains("CalEnergyEM_A" )) {
	  m_CalEnergyEMA  = read_value;
	} else if ( input_string.Contains("CalEnergyEM_B" )) {
	  m_CalEnergyEMB  = read_value;
	} else if ( input_string.Contains("CalEnergyHAD_A" )) {
	  m_CalEnergyHADA = read_value;
	} else if ( input_string.Contains("CalEnergyHAD_B" )) {
	  m_CalEnergyHADB = read_value;
	} else if ( input_string.Contains("CalTransEM_A" )) {
	  m_CalTransEMA   = read_value;
	} else if ( input_string.Contains("CalTransEM_B" )) {
	  m_CalTransEMB   = read_value;
	} else if ( input_string.Contains("CalTransHAD_A" )) {
	  m_CalTransHADA  = read_value;
	} else if ( input_string.Contains("CalTransHAD_B" )) {
	  m_CalTransHADB  = read_value;
	}
      }
    }

    if (input_line.Contains("CalMergeParam:")) {      
      // form of parameters is: name = read_value      
      while( (inpos =fscanf(parFile,"%s",&name[0]) ) > 0){	
	input_string = name; 
	if ( input_string.Contains("end")) break;
	inpos = fscanf(parFile,"%s %s",&equals[0],&temp[0]);
	read_value = strtod(temp,(Char_t**)NULL);	
	
	if      ( input_string.Contains("CalMergeSize_EM" )) {
	  m_EMMergeSize   = read_value;
	} else if ( input_string.Contains("CalMergeAngle_EM" )) {
	  m_EMMergeAngle  = read_value;
	} else if ( input_string.Contains("CalMergeSize_HAD" )) {
	  m_HADMergeSize  = read_value;
	} else if ( input_string.Contains("CalMergeAngle_HAD" )) {
	  m_HADMergeAngle = read_value;
	}
      }
    }

    if(input_line.Contains("Calorimeter:")) {      
      // form of parameters is: name = read_value
      while( (inpos =fscanf(parFile,"%s",&name[0]) ) > 0){      
	input_string = name;       
	if ( input_string.Contains("end")) break;
	inpos = fscanf(parFile,"%s %s",&equals[0],&temp[0]);
	
	if ( input_string.Contains("EMBarrelAbsMat" )) {
	  m_EMBarrelAbsMat = temp;
	} else if ( input_string.Contains("HADBarrelAbsMat" )) {
	  m_HADBarrelAbsMat = temp;
	} else if ( input_string.Contains("MUBarrelAbsMat" )) {
	  m_MUBarrelAbsMat = temp;
	} else{ 	
	  read_value = strtod(temp,(Char_t**)NULL);	
	  if ( input_string.Contains("EMEnergyScale" )) {
	    m_EMEnergyScale  = read_value;
	  } else if ( input_string.Contains("HADEnergyScale" )) {
	    m_HADEnergyScale = read_value;
	  } else if ( input_string.Contains("MUEnergyScale" )) {
	    m_MUEnergyScale  = read_value;
	  } else if ( input_string.Contains("LUMEnergyScale" )) {
	    m_LUMEnergyScale = read_value;
	    
	  } else if ( input_string.Contains("EMThetaSeg" )) {
	    m_EMThetaSeg  = read_value;
	  } else if ( input_string.Contains("HADThetaSeg" )) {
	    m_HADThetaSeg = read_value;
	  } else if ( input_string.Contains("MUThetaSeg" )) {
	    m_MUThetaSeg  = read_value;
	    
	  } else if ( input_string.Contains("LUMphiSeg" )) {
	    m_LUMphiSeg = read_value;
	  } else if ( input_string.Contains("EMphiSeg" )) {
	    m_EMphiSeg  = read_value;
	  } else if ( input_string.Contains("HADphiSeg" )) {
	    m_HADphiSeg = read_value;
	  } else if ( input_string.Contains("MUphiSeg" )) {
	    m_MUphiSeg  = read_value;
	  } else if ( input_string.Contains("LUMphiSeg" )) {
	    m_LUMphiSeg = read_value;
	    
	  } else if ( input_string.Contains("EMBarrelInnerR" )) {
	    m_EMBarrelInnerR  = read_value;
	  } else if ( input_string.Contains("HADBarrelInnerR" )) {
	    m_HADBarrelInnerR = read_value;
	  } else if ( input_string.Contains("MUBarrelInnerR" )) {
	    m_MUBarrelInnerR  = read_value;
	    
	  } else if ( input_string.Contains("EMBarrelOuterR" )) {
	    m_EMBarrelOuterR  = read_value;
	  } else if ( input_string.Contains("HADBarrelOuterR" )) {
	    m_HADBarrelOuterR = read_value;
	  } else if ( input_string.Contains("MUBarrelOuterR" )) {
	    m_MUBarrelOuterR  = read_value;
	    
	  } else if ( input_string.Contains("EMBarrelOuterZ" )) {
	    m_EMBarrelOuterZ  = read_value;
	  } else if ( input_string.Contains("HADBarrelOuterZ" )) {
	    m_HADBarrelOuterZ = read_value;
	  } else if ( input_string.Contains("MUBarrelOuterZ" )) {
	    m_MUBarrelOuterZ  = read_value;
	    
	  } else if ( input_string.Contains("EMBarrelLayer" )) {
	    m_EMBarrelLayer  = (Int_t) read_value;
	  } else if ( input_string.Contains("HADBarrelLayer" )) {
	    m_HADBarrelLayer = (Int_t) read_value;
	  } else if ( input_string.Contains("MUBarrelLayer" )) {
	    m_MUBarrelLayer  = (Int_t) read_value;
	    
	  } else if ( input_string.Contains("EMBarrelThick" )) {
	    m_EMBarrelThick  = read_value;
	  } else if ( input_string.Contains("HADBarrelThick" )) {
	    m_HADBarrelThick = read_value;
	  } else if ( input_string.Contains("MUBarrelThick" )) {
	    m_MUBarrelThick  = read_value;
	    
	  } else if ( input_string.Contains("EMBarrelAbsThick" )) {
	    m_EMBarrelAbsThick  = read_value;
	  } else if ( input_string.Contains("HADBarrelAbsThick" )) {
	    m_HADBarrelAbsThick = read_value;
	  } else if ( input_string.Contains("MUBarrelAbsThick" )) {
	    m_MUBarrelAbsThick  = read_value;
	    
	  } else if ( input_string.Contains("EMEndcapInnerZ" )) {
	    m_EMEndcapInnerZ  = read_value;
	  } else if ( input_string.Contains("HADEndcapInnerZ" )) {
	    m_HADEndcapInnerZ = read_value;
	  } else if ( input_string.Contains("MUEndcapInnerZ" )) {
	    m_MUEndcapInnerZ  = read_value;
	  } else if ( input_string.Contains("LUMEndcapInnerZ" )) {
	    m_LUMEndcapInnerZ = read_value;
	    
	  } else if ( input_string.Contains("EMEndcapLayer" )) {
	    m_EMEndcapLayer  = (Int_t) read_value;
	  } else if ( input_string.Contains("HADEndcapLayer" )) {
	    m_HADEndcapLayer = (Int_t) read_value;
	  } else if ( input_string.Contains("MUEndcapLayer" )) {
	    m_MUEndcapLayer  = (Int_t) read_value;
	  } else if ( input_string.Contains("LUMEndcapLayer" )) {
	    m_LUMEndcapLayer = (Int_t) read_value;
	    
	  } else if ( input_string.Contains("EMEndcapThick" )) {
	    m_EMEndcapThick  = read_value;
	  } else if ( input_string.Contains("HADEndcapThick" )) {
	    m_HADEndcapThick = read_value;
	  } else if ( input_string.Contains("MUEndcapThick" )) {
	    m_MUEndcapThick  = read_value;
	  } else if ( input_string.Contains("LUMEndcapThick" )) {
	    m_LUMEndcapThick = read_value; 
	  }
	}
      }
    }

    if (input_line.Contains("Coil_Pipe:")) {
      while( (inpos =fscanf(parFile,"%s",&name[0]) ) > 0){
	input_string = name;
	if ( input_string.Contains("end")) break;
	inpos = fscanf(parFile,"%s %s",&equals[0],&temp[0]);
	if ( input_string.Contains("CoilPlace" )) {
	  m_CoilPlace = temp;
	}
	read_value = strtod(temp,(Char_t**)NULL);	
	if ( input_string.Contains("CoilInnerR" )) {
	  m_CoilInnerR  = read_value;
	} else if ( input_string.Contains("CoilOuterZ" )) {
	  m_CoilOuterZ  = read_value;
	} else if ( input_string.Contains("CoilThick" )) {
	  m_CoilThick   = read_value;
	} else if ( input_string.Contains("CoilIntLeng" )) {
	  m_CoilIntLeng = read_value;
	} else if ( input_string.Contains("PipeOuter" )) {
	  m_PipeOuter   = read_value;
	}
	
      }
    }
  }
  printf("Tracking Parameters Used: n n");
  printf("PtMin= %f minHit= %i n",m_PtMin,m_minHit);
  printf("PolarInner= %f  n",m_PolarInner);
  printf("PolarOuter= %f n", m_PolarOuter);
  printf("Magnetic Field = %f n", m_Field);
  printf("DetectorName= %s n", m_DetectorName.Data());
  
  // rewind the file again
  fseek(parFile, 0L, 0);

};

 Int_t LCDGetParameters::SetDetectorVolumes(){

  TString* coil_place = GetCoilPlace();
  if       (coil_place->CompareTo("CoilinHAD") == 0) {
    printf("HAD cal is out of Coiln");
  }else if (coil_place->CompareTo("CoiloutHAD") == 0) {
    printf("HAD cal is inner of Coiln");
  }else {
    printf("Unsupported detectorn"); return 0;
  }

  static const Double_t FIELD = GetMagneticField();

  static const Double_t EM_ENDCAP_INNER_Z   = GetEMEndcapInnerZ();
  static const Double_t HAD_ENDCAP_INNER_Z  = GetHADEndcapInnerZ();
  static const Double_t MUON_ENDCAP_INNER_Z = GetMUEndcapInnerZ();
  static const Double_t EM_BARREL_INNER_R   = GetEMBarrelInnerR();
  static const Double_t HAD_BARREL_INNER_R  = GetHADBarrelInnerR();
  static const Double_t MUON_BARREL_INNER_R = GetMUBarrelInnerR();

  static const Double_t THICK_EM  = 
    (GetEMBarrelThick()) *(GetEMBarrelLayer());
  static const Double_t THICK_HAD = 
    (GetHADBarrelThick())*(GetHADBarrelLayer());
  static const Double_t THICK_MUON = 
    (GetMUBarrelThick()) *(GetMUBarrelLayer());

  static const Double_t COIL_INNER_R   = GetCoilInnerR();
  static const Double_t THICK_COIL     = GetCoilThick();

  Double_t Lrad=0.0;
  TString* EMmaterial = GetEMBarrelAbsMat();
  if      (EMmaterial->CompareTo("W") == 0)  { Lrad = 0.35; }
  else if (EMmaterial->CompareTo("Pb") == 0) { Lrad = 0.56; }

  Double_t EM_radlength = 
    (GetEMBarrelAbsThick()) *(GetEMBarrelLayer())/Lrad;
  static const Double_t RADLENPERCM_EM = EM_radlength / THICK_EM;

  Double_t Lint=0.0;
  if      (EMmaterial->CompareTo("W") == 0)  { Lint =  9.6; }
  else if (EMmaterial->CompareTo("Pb") == 0) { Lint = 17.0; }

  Double_t EM_intlength = 
    (GetEMBarrelAbsThick()) *(GetEMBarrelLayer())/Lint;
  static const Double_t ILENPERCM_EM   = EM_intlength / THICK_EM;

  TString* HADmaterial = GetHADBarrelAbsMat();
  if      (HADmaterial->CompareTo("Cu") == 0) { Lint = 14.8; }
  else if (HADmaterial->CompareTo("Pb") == 0) { Lint = 17.0; }
  Double_t HAD_intlength = 
    (GetHADBarrelAbsThick()) *(GetHADBarrelLayer())/Lint;
  static const Double_t ILENPERCM_HAD  = HAD_intlength / THICK_HAD;

  static const Double_t ILENPERCM_COIL =(GetCoilIntLeng())/ THICK_COIL;

  static const Double_t ILENPERCM_MUON = 6.0 / THICK_MUON; // Large
//static const Double_t ILENPERCM_MUON = 7.1 / THICK_MUON; // Small

  static const Double_t OUTER_PIPE =  GetPipeOuter();

  Int_t emBrlIx;
  Int_t coilIx=0;
  Int_t hadBrlIx;
  Int_t emEcIx;
  //  Int_t hadEc1Ix;

  Int_t N_BARREL=0, N_ENDCAP=0;

  if       (coil_place->CompareTo("CoilinHAD") == 0) {
    N_BARREL = 7;           // includes muon
    N_ENDCAP = 6;           // no muon yet
  // nBarrel = 7;  // air, EM, air, coil, had, air, muon.
  // nEndcap = 6;  // air, EM at innermost Z,
  //                  air at next Z
  //                  air, had-1, had-2  at outermost Z
  //                  leaving out muon endcap for now

  }else if (coil_place->CompareTo("CoiloutHAD") == 0) {
    N_BARREL = 6;
    N_ENDCAP = 5;
  }

  m_nVolumes = N_BARREL + N_ENDCAP;

  LCDDetectorVolume* pCur;
  pCur = m_Volumes = new LCDDetectorVolume [m_nVolumes]; 

  Double_t Router[7]; Double_t Zouter[7];
  Router[0] = EM_BARREL_INNER_R;      
  Router[1] = Router[0] + THICK_EM;
  Zouter[0] = EM_ENDCAP_INNER_Z;
  if       (N_BARREL == 7) {
    Router[2] = COIL_INNER_R;
    Router[3] = Router[2] + THICK_COIL;
    Router[4] = Router[3] + THICK_HAD;
    Router[5] = MUON_BARREL_INNER_R;
    Router[6] = Router[5] + THICK_MUON;    
    Zouter[1] = Zouter[2] = Zouter[3] = HAD_ENDCAP_INNER_Z;
    Zouter[4] = Zouter[5] = Zouter[6] = MUON_ENDCAP_INNER_Z;
  }else if (N_BARREL == 6) {
    Router[2] = HAD_BARREL_INNER_R;
    Router[3] = Router[2] + THICK_HAD;
    Router[4] = MUON_BARREL_INNER_R;
    Router[5] = Router[4] + THICK_MUON;
    Zouter[1] = Zouter[2] = HAD_ENDCAP_INNER_Z;
    Zouter[3] = Zouter[4] = Zouter[5] = MUON_ENDCAP_INNER_Z;
  }

  int i=0;
  for (i = 0; i < m_nVolumes; i++) {
    pCur->barrel = (i < N_BARREL);
    pCur->iLenPerCm  = 0.0;  pCur->radLenPerCm = 0.0; pCur->pNextZ = 0;
    if (i < N_BARREL ){
      pCur->outerR = Router[i]; pCur->outerZ = Zouter[i];
      pCur->innerZ = 0.0;     

      if (i == 0){
	pCur->innerR = 0; 
	pCur->pPrevR = 0;
      }else      {
	pCur->innerR = (pCur - 1)->outerR; 
	pCur->pPrevR = pCur - 1;  
      }
    }
    pCur++;
  }

  pCur = m_Volumes;
  for (i = 0; i < N_BARREL-1; i++) { 
    pCur->pNextR = pCur + 1;
    pCur++;
  }
  pCur->pNextR = 0;

  // Tracker volume
  pCur = m_Volumes;

  Int_t  iVol = 0;
  pCur->volKind = VOL_AIR;
  pCur->pNextZ = m_Volumes + N_BARREL;  //adjacent air volume inside EM endcap
  m_InnerVolume = pCur;

  // EM cal
  pCur++;  iVol++;
  emBrlIx = iVol;
  pCur->volKind = VOL_EM;
  pCur->iLenPerCm = ILENPERCM_EM;
  pCur->radLenPerCm =  RADLENPERCM_EM;
  if       (coil_place->CompareTo("CoilinHAD") == 0) {    
    pCur->pNextZ = m_Volumes + N_BARREL + 4;  // adjacent to first Had ec vol.
  }else if (coil_place->CompareTo("CoiloutHAD") == 0) {
    pCur->pNextZ = m_Volumes + N_BARREL + 2;  // next to air following EC endcap
  }

  // Air
  pCur++;  iVol++;
  pCur->volKind = VOL_AIR;
  if       (coil_place->CompareTo("CoilinHAD") == 0) {    
    pCur->pNextZ = m_Volumes + N_BARREL + 4;
  }else if (coil_place->CompareTo("CoiloutHAD") == 0) {
    pCur->pNextZ = (pCur - 1)->pNextZ;    // same as adjacent EM section
  }    

  if       (coil_place->CompareTo("CoilinHAD") == 0) {    
    // Coil
    pCur++;  iVol++; coilIx = iVol;
    pCur->volKind = VOL_COIL;
    pCur->iLenPerCm = ILENPERCM_COIL;
    pCur->pNextZ = m_Volumes + N_BARREL + 5;  // adjacent to 2nd had ec vol.
    
    //HAD Cal 
    pCur++;  iVol++;  hadBrlIx = iVol; 
    pCur->volKind = VOL_HAD;
    pCur->iLenPerCm = ILENPERCM_HAD;
    pCur->pPrevR = pCur + coilIx;
    
    //Air
    pCur++;  iVol++;
    pCur->volKind = VOL_AIR;
    
    // MU Cal
    pCur++;  iVol++;
    pCur->volKind = VOL_MUON;
    pCur->iLenPerCm = ILENPERCM_MUON;
    
    //Now do something similar for endcaps
    pCur++;  iVol++; // next to beam pipe, at Z of EM endcap
    pCur->volKind = VOL_AIR;
    pCur->innerR = 0.0;
    pCur->outerR = OUTER_PIPE;
    pCur->innerZ = EM_ENDCAP_INNER_Z;
    pCur->outerZ = pCur->innerZ + THICK_EM;
    pCur->pPrevR = 0;
    pCur->pNextR = pCur + 1;
    pCur->pNextZ = pCur + 2;
    
    // EM endcap    
    pCur++;  iVol++; emEcIx = iVol;                        
    pCur->volKind = VOL_EM;
    pCur->innerR = OUTER_PIPE;
    pCur->outerR = EM_BARREL_INNER_R;
    pCur->innerZ = EM_ENDCAP_INNER_Z;
    pCur->outerZ = pCur->innerZ + THICK_EM;
    pCur->iLenPerCm = ILENPERCM_EM;
    pCur->radLenPerCm = RADLENPERCM_EM;
    pCur->pNextR = m_Volumes + emBrlIx;
    pCur->pPrevR = pCur - 1;        // air by beampipe
    pCur->pNextZ = pCur + 1;        // air between EM and had endcaps
        
    // just outside of EM endcap in z
    pCur++;  iVol++;              
    pCur->volKind = VOL_AIR;
    pCur->innerR = 0.0;
    pCur->outerR = (pCur - 1)->outerR;
    pCur->innerZ = (pCur - 1)->outerZ;
    pCur->outerZ = HAD_ENDCAP_INNER_Z;
    pCur->pNextR = m_Volumes + emBrlIx;
    pCur->pPrevR = 0;
    pCur->pNextZ = pCur + 1;   // adjacent to another air volume 
    
    // next to beam pipe, at Z of had endcap
    pCur++;  iVol++; 
    pCur->volKind = VOL_AIR;
    pCur->innerR = 0.0;
    pCur->outerR = OUTER_PIPE;
    pCur->innerZ = HAD_ENDCAP_INNER_Z;
    pCur->outerZ = pCur->innerZ + THICK_HAD;
    pCur->pPrevR = 0;
    pCur->pNextR = pCur + 1;
    
    // First EC had
    pCur++; iVol++;
    pCur->volKind = VOL_HAD;
    pCur->innerR = OUTER_PIPE;
    pCur->outerR = COIL_INNER_R;
    pCur->innerZ = HAD_ENDCAP_INNER_Z;
    pCur->outerZ = pCur->innerZ + THICK_HAD;         
    pCur->iLenPerCm = ILENPERCM_HAD;
    pCur->pNextR = pCur + 1;  
    pCur->pPrevR = pCur - 1;
    
    // Outer EC had
    pCur++; iVol++;
    pCur->volKind = VOL_HAD;
    pCur->innerR = (pCur - 1)->outerR;
    pCur->outerR = HAD_BARREL_INNER_R;  
    pCur->innerZ =  HAD_ENDCAP_INNER_Z;
    pCur->outerZ = pCur->innerZ + THICK_HAD;    
    pCur->iLenPerCm = ILENPERCM_HAD;
    pCur->pNextR = m_Volumes + hadBrlIx;
    pCur->pPrevR = pCur - 1;
    
  }else if (coil_place->CompareTo("CoiloutHAD") == 0) {
    pCur++;  iVol++;  hadBrlIx = iVol;       // had
    pCur->volKind = VOL_HAD;
    pCur->iLenPerCm = ILENPERCM_HAD;
    
    pCur++; iVol++;
    pCur->volKind = VOL_COIL;
    pCur->iLenPerCm = ILENPERCM_COIL;
    
    pCur++;  iVol++;
    coilIx = iVol;
    pCur->volKind = VOL_MUON;
    pCur->iLenPerCm = ILENPERCM_MUON;
    
    //Now do something similar for endcaps
    pCur++;  iVol++; // next to beam pipe, at Z of EM endcap
    pCur->volKind = VOL_AIR;
    pCur->innerR = 0.0;
    pCur->outerR = OUTER_PIPE;
    pCur->innerZ = EM_ENDCAP_INNER_Z;
    pCur->outerZ = pCur->innerZ + THICK_EM;
    pCur->pPrevR = 0;
    pCur->pNextR = pCur + 1;
    pCur->pNextZ = pCur + 2;
        
    pCur++;  iVol++; emEcIx = iVol;                        // EM endcap
    pCur->volKind = VOL_EM;
    pCur->innerR = OUTER_PIPE;
    pCur->outerR = EM_BARREL_INNER_R;
    pCur->innerZ = EM_ENDCAP_INNER_Z;
    pCur->outerZ = pCur->innerZ + THICK_EM;
    pCur->iLenPerCm = ILENPERCM_EM;
    pCur->radLenPerCm = RADLENPERCM_EM;
    pCur->pNextR = m_Volumes + emBrlIx;
    pCur->pPrevR = pCur - 1;        // air by beampipe
    pCur->pNextZ = pCur + 1;        // air between EM and had endcaps
        
    pCur++;  iVol++;              // just outside of EM endcap in z
    pCur->volKind = VOL_AIR;
    pCur->innerR = 0.0;
    pCur->outerR = EM_BARREL_INNER_R;
    pCur->innerZ = (pCur - 1)->outerZ;
    pCur->outerZ = HAD_ENDCAP_INNER_Z;         // goes to beginning of ec had
    pCur->pNextR = m_Volumes + emBrlIx;
    pCur->pPrevR = 0;
    pCur->pNextZ = pCur + 1;   // adjacent to air, had ec vol.
    
    pCur++;  iVol++; // next to beam pipe, at Z of had endcap
    pCur->volKind = VOL_AIR;
    pCur->innerR = 0.0;
    pCur->outerR = OUTER_PIPE;
    pCur->innerZ = HAD_ENDCAP_INNER_Z;
    pCur->outerZ = pCur->innerZ + THICK_HAD;
    pCur->pPrevR = 0;
    pCur->pNextR = pCur + 1;
    
    // EC had
    pCur++; iVol++;
    pCur->volKind = VOL_HAD;
    pCur->innerR = OUTER_PIPE;
    pCur->outerR = HAD_BARREL_INNER_R;
    pCur->innerZ = HAD_ENDCAP_INNER_Z;
    pCur->outerZ = pCur->innerZ + THICK_HAD;         
    pCur->iLenPerCm = ILENPERCM_HAD;
    pCur->pNextR = m_Volumes + hadBrlIx;
    pCur->pPrevR = pCur - 1;    
  }

  // Fill in field based on radius
  pCur = m_Volumes;
  for (int j = 0; j < m_nVolumes; j++) {
    if (pCur->outerR <= (m_Volumes + coilIx)->innerR) {
      pCur->field = FIELD;
    } else if (pCur->innerR >= (m_Volumes + coilIx)->outerR) {
      pCur->field = -0.1 * FIELD;
    } else {
      pCur->field = FIELD * .45;   // average field in coil
    }
    pCur++;
  }
  // Take care of field stuff here while we're at it.
  // fieldInner = 6.0; 
  // fieldOuter = -0.1 * (fieldInner);
  // fieldInnerR = pCur->innerR;
  // fieldOuterR = pCur->outerR`;

  return 1;

}