Skip to content
Snippets Groups Projects
CRDEcal.cpp 6.07 KiB
Newer Older
//====================================================================
//  Detector description for CEPC Reference Detector ECal Barrel.
//--------------------------------------------------------------------
//
//  Author     : Fangyi Guo   email: guofangyi@ihep.ac.cn
//
//====================================================================
#include "DD4hep/DetFactoryHelper.h"
#include "XML/Layering.h"
#include "XML/Utilities.h"
#include "DDRec/DetectorData.h"
#include "DDSegmentation/Segmentation.h"

#define MYDEBUG(x) std::cout << __FILE__ << ":" << __LINE__ << ": " << x << std::endl;
#define MYDEBUGVAL(x) std::cout << __FILE__ << ":" << __LINE__ << ": " << #x << ": " << x << std::endl;

using dd4hep::rec::LayeredCalorimeterData;
using namespace dd4hep;
using namespace std;
static dd4hep::Ref_t create_detector(dd4hep::Detector& theDetector,
                                     xml_h e,
                                     dd4hep::SensitiveDetector sens) {

    xml_det_t x_det = e;

    std::string det_name = x_det.nameStr();
    std::string det_type = x_det.typeStr();
    MYDEBUGVAL(det_name);
    MYDEBUGVAL(det_type);
	 int detid = x_det.id();

	 //Global geometry
	 double R0 = theDetector.constant<double>("ecalbarrel_inner_radius");
	 double h0 = theDetector.constant<double>("ecalbarrel_thickness");
	 double Z0 = theDetector.constant<double>("ecalbarrel_zlength");

    double dim_x1 = R0*tan(22.5*degree) + sqrt(2)*h0/2.;
    double dim_x2 = dim_x1 - h0;
    double dim_y = Z0/2.;
    double dim_z = h0/2.;		                 
	 double dx = dim_x1 - R0*tan(22.5*degree);  //transport distance in x-axis
    double r0 = R0+h0/2.;							  //rotation radius 

	 //Crystal bar size
    double barx = 10*mm;             //Crystal size in R direction. 
    double bary = 10*mm;             //Crystal size in z/phi direction (z for odd layer, phi for even layer).
	 int Nlayers = (int)h0/barx;
	 int Nbarz_odd = (int)Z0/bary;
	 int Nbarphi_odd; 	  				 //Depends on layer. Layer1~5: 5. Layer7~23: 4. Layer25~27: 3. 
	 int Nbarz_even = 12; 
	 int Nbarphi_even;                //Depends on layer
	 double barz_even = Z0/Nbarz_even; //~38 cm
	 double barz_odd; 					  //Depends on layer

    //Define detector and motherVolume(world)
    dd4hep::DetElement ECAL(det_name, detid);
    dd4hep::Volume motherVol = theDetector.pickMotherVolume(ECAL);

	 // Create a Tube-like envelope representing the whole detector volume
	 dd4hep::Tube		envelope(R0, (R0+h0)/cos(22.5*degree), Z0/2, 0, 2*M_PI);
	 dd4hep::Material	air(theDetector.material("Air"));
	 dd4hep::Volume	envelopeVol(det_name, envelope, air);
	 dd4hep::PlacedVolume	envelopePlv = motherVol.placeVolume(envelopeVol, Position(0,0,0));
	 envelopeVol.setVisAttributes(theDetector, x_det.visStr() );
	 ECAL.setPlacement(envelopePlv);

//    dd4hep::Material det_pipe(theDetector.material("Steel235"));
//    dd4hep::Volume pipe_vol("pipe", dd4hep::Tube(10*mm, 20*mm, 2300*mm), det_pipe);
//    dd4hep::PlacedVolume pipePhv = motherVol.placeVolume(pipe_vol, Position(0.,0.,0.));

    //Define specific material and volumen for detElement
    dd4hep::Material mat_BGO(theDetector.material("G4_BGO"));
	 dd4hep::Trapezoid trap(dim_x1, dim_x2, dim_y, dim_y, dim_z);
    dd4hep::Volume det_vol("trap_vol", trap, mat_BGO);
//	 det_vol.setVisAttributes(theDetector, x_det.visStr());
	 dd4hep::DetElement stavedet(ECAL, "trap",detid);

	 //Loop to place crystalls in one part
	 //Outer loop: layer (odd layer). 
	 for(int ilayer=1; ilayer<=Nlayers; ilayer+=2){
		double lx = dim_x1 - ilayer*10*mm;
		//Loop in phi direction
		if(ilayer<=5) Nbarphi_odd=5; else if(ilayer<=23) Nbarphi_odd=4; else Nbarphi_odd=3;
		for(int iphi=1; iphi<=Nbarphi_odd; iphi++){
			barz_odd = 2*lx/Nbarphi_odd; 
			//Loop in Z direction
			dd4hep::Volume bar_odd("box_bar", dd4hep::Box(barz_odd/2, bary/2, barx/2), mat_BGO);
    		bar_odd.setSensitiveDetector(sens);
			for(int iz=1; iz<=Nbarz_odd;iz++){
				dd4hep::PlacedVolume plv = det_vol.placeVolume(bar_odd, Position((2*iphi-1)*barz_odd/2-lx, (2*iz-1)*barx/2-dim_y, (2*ilayer-1)*bary/2-dim_z));
				plv.addPhysVolID("layer", ilayer).addPhysVolID("piece", iphi).addPhysVolID("bar", iz);
				std::string barname = "CrystalBar_"+std::to_string(ilayer)+"_"+std::to_string(iphi)+"_"+std::to_string(iz);
				dd4hep::DetElement sd(stavedet, barname, detid);
				sd.setPlacement(plv);
			}
		}
	 }

	 //Loop in even layer
    dd4hep::Volume bar_even("box_bar", dd4hep::Box(bary/2, barz_even/2, barx/2), mat_BGO);
    bar_even.setSensitiveDetector(sens);
    for(int ilayer=2; ilayer<=Nlayers; ilayer+=2){
      double lx = dim_x1 - ilayer*10*mm;
      //Loop in phi direction
		Nbarphi_even = (int)floor(2*lx/bary);
      for(int iphi=1; iphi<=Nbarphi_even; iphi++){
         //Loop in Z direction
         for(int iz=1; iz<=Nbarz_even;iz++){
            dd4hep::PlacedVolume plv = det_vol.placeVolume(bar_even, Position(lx-(2*iphi-1)*bary/2, (2*iz-1)*barz_even/2-dim_y, (2*ilayer-1)*bary/2-dim_z));
				plv.addPhysVolID("layer", ilayer).addPhysVolID("piece", iphi).addPhysVolID("bar", iz);
            std::string barname = "CrystalBar_"+std::to_string(ilayer)+"_"+std::to_string(iphi)+"_"+std::to_string(iz);
            dd4hep::DetElement sd(stavedet, barname, detid);
            sd.setPlacement(plv);
         }
      }
    }

    //create a unique volume (PlacedVolume) for the detector
//    dd4hep::Transform3D trans0(dd4hep::RotationX(-90*degree),  dd4hep::Position(-dx, r0, 0.));
//    dd4hep::PlacedVolume plv0 = motherVol.placeVolume(det_vol, trans0);
//	 stavedet.setPlacement(plv0);
	 for(int i=0;i<8;i++){
		 double rotAngle = 45*i*degree;
		 double posx = -r0*sin(rotAngle) - dx*cos(rotAngle);
		 double posy = r0*cos(rotAngle) - dx*sin(rotAngle);
	    dd4hep::Transform3D transform(dd4hep::RotationZ(rotAngle)*dd4hep::RotationX(-90*degree),  dd4hep::Position(posx, posy, 0.));
   	 dd4hep::PlacedVolume plv = envelopeVol.placeVolume(det_vol, transform);
		 plv.addPhysVolID("block", i);
		 DetElement sd(ECAL, _toString(i,"trap%3d"), detid);
	 	 sd.setPlacement(plv);
	 }

	 sens.setType("calorimeter");

    MYDEBUG("create_detector DONE. ");
    return ECAL;
}

DECLARE_DETELEMENT(CRDEcalBarrel, create_detector)