/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Emc/EmcRecGeoSvc/EmcRecGeoSvc-01-01-07/src/EmcRecBarrelGeo.cxx

Go to the documentation of this file.

//
//   EmcRecBarrelGeo
//
//   Dec 18, 2003, Created by Wang.Zhe
//
//   unit: mm, radian
//
//#include "BesGeoEMC/DB2BesGeoEMC.h"
#include "EmcRecGeoSvc/EmcRecBarrelGeo.h"

EmcRecBarrelGeo::EmcRecBarrelGeo()
{
   ParameterInitialize();
   CalculateStandardCrystal();
   Transform2Column1();
   FillCCenterVector();
}

EmcRecBarrelGeo::~EmcRecBarrelGeo()
{
   fStandard.clear();
   fCCenter.clear();
   fCFrontCenter.clear();
}
  
void EmcRecBarrelGeo::ParameterInitialize()
{
   /*
   DB2BesGeoEMC obj1;
   vector<BesGeoEMC> VecA;
   
   obj1.Get_BesGeoEMC(VecA);
   fBarrelR=VecA[0].Get_BarrelR();
   fBarrelOffset1=VecA[0].Get_BarrelOffset1();
   fBarrelOffset2=VecA[0].Get_BarrelOffset2();
   fBarrelh1=VecA[0].Get_Barrelh1();
   fBarrelh2=VecA[0].Get_Barrelh2();
   fBarrelh3=VecA[0].Get_Barrelh3();
   fBarrelL=VecA[0].Get_BarrelL();
   fBarrelL2=VecA[0].Get_BarrelL2();
   fBarrelNPhiMax=(int)VecA[0].Get_BarrelNPhiMax();
   fBarrelNThetaMax=(int)VecA[0].Get_BarrelNThetaMax();
   */
    
   fBarrelR=942;
   fBarrelOffset1=25;
   fBarrelOffset2=50;
   fBarrelh1=51;
   fBarrelh2=52;
   fBarrelh3=52.466;
   fBarrelL=280;
   fBarrelL2=240; 
   fBarrelNPhiMax=120;
   fBarrelNThetaMax=22;
   
   fBarrelAlpha=twopi/fBarrelNPhiMax;

   fStandard.clear();
   fCCenter.clear();
   fCFrontCenter.clear();
}

void EmcRecBarrelGeo::CalculateStandardCrystal()
{
   double R1,R2,a;
   EmcRecCrystal pre,now;
   double dx,dy,dz;
   HepPoint3D t1,t2;
   double L,h;
   
   HepPoint3D O1(0,0,0);
   HepPoint3D O2(0,0,fBarrelOffset1);
   HepPoint3D O3(0,0,fBarrelOffset2);

   R1=2*fBarrelR*sin(fBarrelAlpha/2)/sin(fBarrelAlpha);
   R2=2*fBarrelR*sin(fBarrelAlpha/2)*(tan(fBarrelAlpha)+1/tan(fBarrelAlpha));
   a=2*fBarrelR*sin(fBarrelAlpha/2)/cos(fBarrelAlpha);
   
   HepPoint3D M4(0,R1,0);
   HepPoint3D A1(a,R1,0);
   
//   cout<<O1<<", "<<O2<<", "<<O3<<endl;
//   cout<<M4<<", "<<A1<<endl;
//   cout<<R1<<", "<<R2<<endl;
   
   // No. 1  theta: 21; phi: 30
   EmcRecCrystal Crystal1;

   Crystal1.Set(0,M4);
   Crystal1.Set(1,A1);
   Crystal1.Set(2,A1);   Crystal1.SetZ(2,A1.z()+fBarrelh1);
   Crystal1.Set(3,M4);   Crystal1.SetZ(3,M4.z()+fBarrelh1);
   
   Crystal1.Set(4,M4);   Crystal1.SetY(4,M4.y()+fBarrelL);
   dx=a*(R1+fBarrelL)/R1;
   Crystal1.Set(5,Crystal1.Get(4));   Crystal1.SetX(5,dx);
   dz=fBarrelOffset1+(fBarrelh1-fBarrelOffset1)*(R1+fBarrelL)/R1;
   Crystal1.Set(6,Crystal1.Get(5));   Crystal1.SetZ(6,Crystal1.Get(5).z()+dz);
   Crystal1.Set(7,Crystal1.Get(4));   Crystal1.SetZ(7,Crystal1.Get(4).z()+dz);
   
   fStandard.push_back(Crystal1);
   pre=fStandard[0];
   
   // No. 2  theta: 20; phi: 30
   EmcRecCrystal Crystal2; 
   double sin_gamma,cos_gamma,tan_gamma;
   
   tan_gamma=R1/(fBarrelh1-fBarrelOffset1);
   sin_gamma=tan_gamma/(sqrt(1+tan_gamma*tan_gamma));
   cos_gamma=1/(sqrt(1+tan_gamma*tan_gamma));
   
   double aa,bb,Rp;
   aa=fBarrelh1/sin_gamma;
   bb=fBarrelh1/tan_gamma; 
   
   Crystal2.Set(0,pre.Get(3));   Crystal2.SetZ(0,pre.Get(3).z()+bb*cos_gamma);
                                 Crystal2.SetY(0,pre.Get(3).y()+bb*sin_gamma);
   Rp=R1/sin_gamma;
   dx=a*(Rp+bb)/Rp;
   Crystal2.Set(1,Crystal2.Get(0));   Crystal2.SetX(1,dx);   
   Crystal2.Set(3,pre.Get(3));        Crystal2.SetZ(3,pre.Get(3).z()+aa);
   Crystal2.Set(2,Crystal2.Get(3));   Crystal2.SetX(2,a);
   
   dz=(fBarrelL+bb)*cos_gamma;
   dy=(fBarrelL+bb)*sin_gamma;
   Crystal2.Set(4,pre.Get(3));   Crystal2.SetZ(4,pre.Get(3).z()+dz);
                                 Crystal2.SetY(4,pre.Get(3).y()+dy); 
   dx=a*(Rp+bb+fBarrelL)/Rp; 
   Crystal2.Set(5,Crystal2.Get(4));     Crystal2.SetX(5,dx);
   double gam2,gam3,Lp,Rpp;
   t1=pre.Get(3)-Crystal2.Get(3);
   t2=O3-Crystal2.Get(3);
   gam2=acos(t1*t2/(t1.mag()*t2.mag()));
   gam3=acos(cos_gamma)-gam2;
   Lp=fBarrelL/cos(gam3);
   dz=Lp*cos(gam2);
   dy=Lp*sin(gam2);
   Crystal2.Set(7,Crystal2.Get(3));     Crystal2.SetZ(7,Crystal2.Get(3).z()+dz);
                                        Crystal2.SetY(7,Crystal2.Get(3).y()+dy);
   Rpp=R1/sin(gam2);
   dx=a*(Rpp+Lp)/Rpp;
   Crystal2.Set(6,Crystal2.Get(7));     Crystal2.SetX(6,dx);
   
   fStandard.push_back(Crystal2);

   // No. 3 -- No. 22  theta: 19-0; phi: 30
   for(int i=3; i<=fBarrelNThetaMax; i++) {
      pre=fStandard[i-2];
      
      if(i<fBarrelNThetaMax) {
         L=fBarrelL;
         h=fBarrelh2;
      }
      else {
         L=fBarrelL2;
         h=fBarrelh3;
      }
      
      tan_gamma=R1/(pre.Get(3).z()-fBarrelOffset2);
      sin_gamma=tan_gamma/(sqrt(1+tan_gamma*tan_gamma));
      cos_gamma=1/(sqrt(1+tan_gamma*tan_gamma));

      aa=h/sin_gamma;
      bb=h/tan_gamma;
      
      now.Set(3,pre.Get(3));   now.SetZ(3,pre.Get(3).z()+aa);
      now.Set(0,pre.Get(3));   now.SetZ(0,pre.Get(3).z()+bb*cos_gamma);
                               now.SetY(0,pre.Get(3).y()+bb*sin_gamma);
      
      now.Set(4,pre.Get(3));   now.SetZ(4,pre.Get(3).z()+(L+bb)*cos_gamma);
                               now.SetY(4,pre.Get(3).y()+(L+bb)*sin_gamma);
      
      gam2=atan(R1/(now.Get(3).z()-fBarrelOffset2));
      gam3=acos(cos_gamma)-gam2;
      
      Lp=L/cos(gam3);
      dz=Lp*cos(gam2);
      dy=Lp*sin(gam2);
   
      now.Set(7,now.Get(3));   now.SetZ(7,now.Get(3).z()+dz);
                               now.SetY(7,now.Get(3).y()+dy);

      Rp=R1/sin_gamma;
      dx=a*(Rp+bb)/Rp;
      now.Set(1,now.Get(0));   now.SetX(1,dx);
      dx=a*(Rp+bb+L)/Rp;
      now.Set(5,now.Get(4));   now.SetX(5,dx);
      
      now.Set(2,now.Get(3));   now.SetX(2,a);
      
      Rpp=R1/sin(gam2);
      dx=a*(Rpp+Lp)/Rpp;
      now.Set(6,now.Get(7));   now.SetX(6,dx);
      
      fStandard.push_back(now);
   }
   
//   for(int i=1; i<=fBarrelNThetaMax; i++) {
//      cout<<endl;
//      cout<<"No. "<<i<<endl;
//      now=fStandard[i-1];
//      now.BarrelCheckout();
//      now.Dump();
//   }  
}

void EmcRecBarrelGeo::Transform2Column1()  // to theta: 21-0; phi: 0
{
   HepPoint3D O1(fBarrelR*sin(fBarrelAlpha/2),
                 fBarrelR*cos(fBarrelAlpha/2)-2*fBarrelR*sin(fBarrelAlpha/2)/tan(fBarrelAlpha),
                 0);
   
   HepPoint3D R=fStandard[0].Get(5);
   double phi=(R.rotateZ(fBarrelAlpha)+O1).phi();

   for(int i=1; i<=fBarrelNThetaMax; i++) {
      for(int m=0;m<8;m++)
        {
           fStandard[i-1].Type(EmcRecCrystal::SixPlane());
           fStandard[i-1].Set(m,fStandard[i-1].Get(m).rotateZ(fBarrelAlpha));
           fStandard[i-1].Set(m,fStandard[i-1].Get(m)+O1);
           fStandard[i-1].Set(m,fStandard[i-1].Get(m).rotateZ(-phi));
        }

//      cout<<endl;
//      cout<<"No. "<<i<<endl;
//      fStandard[i-1].BarrelCheckout();
//      fStandard[i-1].Dump();
   }
}

void EmcRecBarrelGeo::FillCCenterVector()
{
   int phi,theta;
   EmcRecCrystal aCry;
   HepPoint3D aCenter;
   HepPoint3D aFrontCenter;
   
   for(theta=0;theta<2*fBarrelNThetaMax;++theta) {
      for(phi=0;phi<fBarrelNPhiMax;++phi) {
         aCry=GetCrystal(EmcID::crystal_id(EmcID::getBARREL(),theta,phi));
         aCenter=aCry.Center();
   aFrontCenter=aCry.FrontCenter();
         fCCenter.push_back(aCenter);
   fCFrontCenter.push_back(aFrontCenter);
//       aCry.Dump();
//       cout<<"("<<theta<<","<<phi<<")  ";
//       cout<<aCenter<<endl;
      }
   }
}


// Access
EmcRecCrystal EmcRecBarrelGeo::GetCrystal(const Identifier& id) const
{
   EmcRecCrystal cry;
   unsigned int theta=EmcID::theta_module(id);
   unsigned int phi=EmcID::phi_module(id);
   
   double dphi=phi*fBarrelAlpha;
   
   if((int)theta>=fBarrelNThetaMax) {
      cry=fStandard[theta-fBarrelNThetaMax];
      for(int m=0;m<8;++m)
        {
           cry.SetZ(m,-cry.Get(m).z());
        }
   }
   else {
      cry=fStandard[fBarrelNThetaMax-theta-1];
   }
   
   for(int m=0;m<8;++m)
     {
        cry.Set(m,cry.Get(m).rotateZ(dphi));
     }
   
   return cry;
}

HepPoint3D EmcRecBarrelGeo::GetCCenter(const Identifier& id) const
{
   unsigned int theta=EmcID::theta_module(id);
   unsigned int phi=EmcID::phi_module(id);
   
   return fCCenter[theta*fBarrelNPhiMax+phi];
}

HepPoint3D EmcRecBarrelGeo::GetCFrontCenter(const Identifier& id) const
{
  unsigned int theta=EmcID::theta_module(id);
  unsigned int phi=EmcID::phi_module(id);

  return fCFrontCenter[theta*fBarrelNPhiMax+phi];
}

---