//*--- Author : Vladimir Pechenov
//*--- Modified: 16.08.05 V.Pechenov
//*--- Modified: 11.05.02 by V. Pechenov
//*--- Modified: 08.05.01 by V. Pechenov
//*--- Modified: 07.03.01 by V. Pechenov
//*--- Modified: 27.05.00 by V. Pechenov
//*--- Modified: 24.05.00 by V. Pechenov
//*--- Modified: 07.03.00 by R. Holzmann
//*--- Modified: 31.08.99 by V. Pechenov
using namespace std;
#include "hmdcsizescells.h"
#include "hspectrometer.h"
#include "hmdcdetector.h"
#include "hmdcgeomstruct.h"
#include <iostream>
#include <iomanip>
#include "math.h"
#include "hgeomvolume.h"
#include "hgeomcompositevolume.h"
#include "hmdcgetcontainers.h"
#include "hmdclayergeompar.h"
#include "hmdcgeompar.h"
#include "hmdclookupraw.h"
#include "hmdcrawstruct.h"
#include "hruntimedb.h"
#include "hspecgeompar.h"
#include <stdlib.h>
//_HADES_CLASS_DESCRIPTION
//////////////////////////////////////////////////////////////////////////////
//
// HMdcSizesCells
// HMdcSizesCellsSec
// HMdcSizesCellsMod
// HMdcSizesCellsLayer
// HMdcSizesCellsCell
//
// Creating the container and geting the pointer to cont.:
// HMdcSizesCells* fSizesCells=HMdcSizesCells::getObject();
//
// Geting the pointer to cont. (without creating):
// HMdcSizesCells* fSizesCells=HMdcSizesCells::getExObject();
//
// Deleting of container:
// HMdcSizesCells::deleteCont();
//
// Container class keep volumes sizes of wires (2-points/wire)
// in coordinate system of sector and transformations layer<->mod,
// layer<->sec, ...
// It has same useful functions:
// HMdcSizesCellsCell
// Float_t getWireLength(void)
//
// HMdcSizesCellsLayer:
// HGeomTransform* getSecTrans();
// HGeomTransform* getModTrans();
// void transSecToRotLayer(HGeomVector& point);
// Double_t getAlpha(const HGeomVector& p1, const HGeomVector& p2);
// Double_t getDist(const HGeomVector& p1,
// const HGeomVector& p2, Int_t cell);
// void getImpact(const HGeomVector& p1, const HGeomVector& p2,
// Int_t cell, Double_t &alpha, Double_t &per);
// Double_t getImpact(Double_t x1,Double_t y1,Double_t z1,
// Double_t x2,Double_t y2,Double_t z2, Int_t cell, Double_t &alphaPerp);
// Bool_t getImpact(Double_t x1,Double_t y1,Double_t z1,
// Double_t x2, Double_t y2, Double_t z2,
// Int_t cell, Double_t &alphaDrDist, Double_t &driftDist);
// Int_t getCellNum(Double_t x, Double_t y):
// calculation the the cell number for the point x,y
// on the wire plane. (x,y - in coor. sys. of layer, z=0)
// Double_t getXSign(Double_t x1, Double_t y1, Double_t z1,
// Double_t x2, Double_t y2, Double_t z2, Int_t cell)
// Float_t getSignPath(Double_t x1,Double_t y1,Double_t z1,
// Double_t x2, Double_t y2, Double_t z2, Int_t cell);
// Float_t getSignPath(Double_t x1,Double_t y1,Double_t z1,
// Double_t x2, Double_t y2, Double_t z2, Int_t cell, Float_t& outside);
//
// HMdcSizesCellsMod:
// HGeomTransform* getSecTrans();
//
// HMdcSizesCellsSec:
// HGeomTransform* getLabTrans();
//
// HMdcSizesCells:
// Bool_t getCellVol(Int_t sec, Int_t mod, Int_t lay,
// Int_t cell, HMdcTrap& volCell);
// Char_t findReadOutSide(Int_t s,Int_t m,Int_t l,Int_t c)
//
//////////////////////////////////////////////////////////////////////////////
ClassImp(HMdcSizesCellsCell)
ClassImp(HMdcSizesCellsLayer)
ClassImp(HMdcSizesCellsMod)
ClassImp(HMdcSizesCellsSec)
ClassImp(HMdcSizesCells)
HMdcSizesCells* HMdcSizesCellsLayer::getSizesCells(void) {
return pToMod->getSizesCells();
}
HMdcSizesCells* HMdcSizesCellsMod::getSizesCells(void) {
return pToSec->getSizesCells();
}
HMdcSizesCells* HMdcSizesCellsSec::getSizesCells(void) {
return pToSC;
}
void HMdcSizesCellsCell::clear(void) {
wirePnt1.clear();
wirePnt2.clear();
status=kFALSE;
xReadout=0.;
readOutSide='0';
}
HMdcSizesCellsLayer::HMdcSizesCellsLayer(HMdcSizesCellsMod* pMod, Int_t lay) {
// Constructor creates an array of pointers of type HMdcSizesCellsCell
sector=pMod->getSec();
module=pMod->getMod();
layer=lay;
pToMod=pMod;
HMdcGetContainers* fGetCont=HMdcGetContainers::getObject();
HMdcGeomStruct* fMdcGeomStruct=fGetCont->getMdcGeomStruct();
nCells=((*fMdcGeomStruct)[sector][module])[layer];
if( nCells > 0 ) {
array = new TClonesArray("HMdcSizesCellsCell",nCells);
for (Int_t cell = 0; cell < nCells; cell++)
new((*array)[cell]) HMdcSizesCellsCell(cell);
}
}
HMdcSizesCellsLayer::~HMdcSizesCellsLayer(void) {
// destructor
if(array) delete array;
array=0;
}
Int_t HMdcSizesCellsLayer::getSize(void) const {
// return the size of the pointer array
return array->GetLast()+1;
}
Bool_t HMdcSizesCellsLayer::setSecTrans(Double_t corZ) {
// alignTrans: coord. system of layer will transformed to "alignTrans" sys.
HMdcGetContainers* fGetCont=HMdcGetContainers::getObject();
if(!fGetCont) return kFALSE;
if(!fGetCont->getModTransLayerZ0(sysRMod,module,layer)) return kFALSE;
HGeomVector vc(sysRMod.getTransVector());
if(corZ!=0) { // For alignment.
vc.setZ(vc.getZ()+corZ);
sysRMod.setTransVector(vc);
}
sysRSec=sysRMod;
sysRSec.transFrom(pToMod->sysRSec);
HMdcSizesCells::copyTransfToArr(sysRSec,tSysRSec);
setPlanePar(sysRSec);
rotLaySysRSec.clear();
HGeomRotation rotLay;
rotLay.setUnitMatrix();
rotLay.setElement( cosWireOr,0);
rotLay.setElement(-sinWireOr,1);
rotLay.setElement( sinWireOr,3);
rotLay.setElement( cosWireOr,4);
rotLaySysRSec.setRotMatrix(rotLay);
rotLaySysRSec.transFrom(sysRSec);
HMdcSizesCells::copyTransfToArr(rotLaySysRSec,tRLaySysRSec);
return kTRUE;
}
Int_t HMdcSizesCellsLayer::calcCellNum(Double_t x, Double_t y) const{
// calculation the the cell number for the point x,y
// on the wire plane. (x,y - in coor. sys. of layer, z=0)
Double_t wire=(wireOffSet + y*cosWireOr-x*sinWireOr)/pitch+0.5;
Int_t nWire=(Int_t)wire;
return (wire<0. || nWire>=nCells) ? -1:nWire;
}
Bool_t HMdcSizesCellsLayer::calcCrossedCells(Double_t x1, Double_t y1,
Double_t z1, Double_t x2, Double_t y2, Double_t z2,
Float_t& cell1, Float_t& cell2) const {
// calculation of cells crossed by line x1,y1,z1-x2,y2,z2
// x1,y1,z1, x2,y2,z2 - sector coordinate system
transTo(x1,y1,z1);
transTo(x2,y2,z2);
Double_t dX=x2-x1;
Double_t dY=y2-y1;
Double_t dZ=z2-z1;
if(dZ==0.) return kFALSE;
Double_t tmp=(halfCatDist-z1)/dZ;
Double_t x=tmp*dX+x1;
Double_t y=tmp*dY+y1;
cell1=(wireOffSet + y*cosWireOr-x*sinWireOr)/pitch+0.5;
tmp=(-halfCatDist-z1)/dZ;
x=tmp*dX+x1;
y=tmp*dY+y1;
cell2=(wireOffSet + y*cosWireOr-x*sinWireOr)/pitch+0.5;
if(cell1>cell2) {
Float_t c=cell1;
cell1=cell2;
cell2=c;
}
if(cell1<0.) {
if(cell2<0.) return kFALSE;
cell1=0.;
}
if(cell2>=nCells) {
if(cell1>=nCells) return kFALSE;
cell2=nCells-0.0001;//nCells-1;
}
return kTRUE;
}
void HMdcSizesCellsLayer::transSecToRotLayer(HGeomVector& point) const {
// transformation to coor. sys. of rotated layer from sector sys.
point=rotLaySysRSec.transTo(point);
}
Double_t HMdcSizesCellsLayer::getAlpha(const HGeomVector& p1,
const HGeomVector& p2) const {
// return the angle alpha (in degree) between projection of linep1-p2
// on the Y-Z plane and Y axis in coor.sys. of wires in degree.
// p1,p2 in coordinate system of sector.
HGeomVector tranP1(p1);
HGeomVector tranP2(p2);
tranP1=rotLaySysRSec.transTo(tranP1);
tranP2=rotLaySysRSec.transTo(tranP2);
return atan2( fabs(tranP2.getZ()-tranP1.getZ()),
fabs(tranP2.getY()-tranP1.getY()) )*TMath::RadToDeg();
}
Double_t HMdcSizesCellsLayer::getDist(const HGeomVector& p1,
const HGeomVector& p2,Int_t cell){
// return the minimal distance from line (p1-p2) to wire.
// p1,p2 are in coordinate system of sector.
HGeomVector tranP1(p1);
HGeomVector tranP2(p2);
tranP1=rotLaySysRSec.transTo(tranP1);
tranP2=rotLaySysRSec.transTo(tranP2);
Double_t yWire=cell*pitch-wireOffSet;
Double_t y1=tranP1.getY() - yWire;
Double_t y2=tranP2.getY() - yWire;
Double_t dz=tranP1.getZ()-tranP2.getZ();
Double_t dy=y1-y2;
return fabs((tranP1.getZ()*y2-tranP2.getZ()*y1)/sqrt(dz*dz+dy*dy));
}
void HMdcSizesCellsLayer::getImpact(const HGeomVector& p1,
const HGeomVector& p2, Int_t cell, Double_t &alpha, Double_t &minDist) const {
// calc. the angle alpha (in degree) between projection of line p1-p2
// on the Y-Z plane and Y axis in coor.sys. of wire (y=z=0) and
// the minimal distance from line (p1-p2) to wire.
// p1,p2 - in sector coor.sys.
HGeomVector tranP1(p1);
HGeomVector tranP2(p2);
tranP1=rotLaySysRSec.transTo(tranP1);
tranP2=rotLaySysRSec.transTo(tranP2);
Double_t yWire=cell*pitch-wireOffSet;
Double_t y1=tranP1.getY() - yWire;
Double_t y2=tranP2.getY() - yWire;
Double_t dz=tranP2.getZ()-tranP1.getZ();
Double_t dy=y2-y1;
Double_t lng=1./sqrt(dz*dz+dy*dy);
alpha=asin(fabs(dz)*lng)*TMath::RadToDeg();
minDist = fabs((tranP1.getZ()*y2-tranP2.getZ()*y1)*lng);
}
Double_t HMdcSizesCellsLayer::HMdcSizesCellsLayer::getImpact(
Double_t x1,Double_t y1,Double_t z1,
Double_t x2,Double_t y2,Double_t z2, Int_t cell,Double_t &alpha) const {
// calc. the angle alpha (in degree) between projection of line
// x1,y1,z1-x2,y2,z2 on the Y-Z plane and Y axis in coor.sys.
// of wires and the minimal distance from line x1,y1,z1-x2,y2,z2
// to wire.
// x1,y1,z1,x2,y2,z2 - in sector coor.sys.
Double_t y,z,dy,dz;
rotYZtoRotLay(x1-tRLaySysRSec[9],y1-tRLaySysRSec[10],z1-tRLaySysRSec[11],y,z);
y += wireOffSet - cell*pitch;
rotYZtoRotLay(x2-x1,y2-y1,z2-z1,dy,dz);
alpha=atan2(fabs(dz),fabs(dy))*TMath::RadToDeg();
return fabs((y*dz-z*dy)/sqrt(dz*dz+dy*dy));
}
Bool_t HMdcSizesCellsLayer::getImpact(Double_t x1, Double_t y1, Double_t z1,
Double_t x2, Double_t y2, Double_t z2,
Int_t cell, Double_t &alpha, Double_t &minDist) const {
// return kTRUE if line x1,y1,z1-x2,y2,z2 intersect cell or kFALSE
// if not intersect and calc. the angle alpha (in degree) between projection
// of the line x1,y1,z1-x2,y2,z2 on the Y-Z plane and Y axis in coor.
// sys. of wires and the minimal distance from line x1,y1,z1-x2,y2,z2
// to wire.
// x1,y1,z1,x2,y2,z2 - in sector coor.sys.
Double_t y,z,dy,dz;
rotYZtoRotLay(x1-tRLaySysRSec[9],y1-tRLaySysRSec[10],z1-tRLaySysRSec[11],y,z);
y += wireOffSet - cell*pitch;
rotYZtoRotLay(x2-x1,y2-y1,z2-z1,dy,dz);
Double_t lng=1./sqrt(dz*dz+dy*dy);
Double_t yDist = fabs(y*dz-z*dy); // abs(Ywire-Ycross_track)=yDist/dz
minDist = yDist*lng;
dz=fabs(dz);
dy=fabs(dy);
alpha=atan2(dz,dy)*TMath::RadToDeg();
if(minDist*dz*lng > halfPitch) {
if(dy==0.0) return kFALSE;
if((yDist-halfPitch*dz)/dy > halfCatDist) return kFALSE;
} else if(minDist*dy*lng > halfCatDist && // dy*lng = cos(alpha)
(yDist-halfCatDist*dy)/dz > halfPitch) return kFALSE;
return kTRUE;
}
Double_t HMdcSizesCellsLayer::getXSign(Double_t x1, Double_t y1,
Double_t z1, Double_t x2, Double_t y2, Double_t z2, Int_t cell) const {
// return X position of the signal on the wire in coor. sys. of this wire
// (Ywire=0, Zwire=0), x1,y1,z1 & x2,y2,z2 - track in sector coor.system
Double_t x,y,z;
rotVectToRotLay(x1-tRLaySysRSec[9],y1-tRLaySysRSec[10],z1-tRLaySysRSec[11],
x,y,z);
y += wireOffSet - cell*pitch;
Double_t dx,dy,dz;
rotVectToRotLay(x2-x1,y2-y1,z2-z1,dx,dy,dz);
return x-(z*dz+y*dy)/(dz*dz+dy*dy)*dx;
}
Float_t HMdcSizesCellsLayer::getSignPath(Double_t x1, Double_t y1, Double_t z1,
Double_t x2, Double_t y2, Double_t z2, Int_t cell) const {
// return the path of signal along wire number "cell"
// x1,y1,z1 & x2,y2,z2 - track in sector coor.system
// Path can be <0. or > wire_length !!!
// For GEANT data and don't connected wires return 0.
HMdcSizesCellsCell& fCell=(*this)[cell];
if(!&fCell || fCell.readOutSide=='0') return 0.;
Float_t x=getXSign(x1, y1,z1, x2, y2, z2,cell);
return (fCell.readOutSide=='1') ? fCell.xReadout-x : x-fCell.xReadout;
}
Double_t HMdcSizesCellsLayer::getXSign(const HMdcLineParam& ln,
Int_t cell) const {
// return X position of the signal on the wire in coor. sys. of this wire
// (Ywire=0, Zwire=0), "ln" - track in sector coor.system
Double_t x,y,z,dx,dy,dz;
rotVectToRotLay(ln.x1()-tRLaySysRSec[9],ln.y1()-tRLaySysRSec[10],
ln.z1()-tRLaySysRSec[11],x,y,z);
y += wireOffSet - cell*pitch;
rotVectToRotLay(ln.getDir(),dx,dy,dz);
return x-(z*dz+y*dy)/(dz*dz+dy*dy)*dx;
}
Double_t HMdcSizesCellsLayer::getSignPath(const HMdcLineParam& ln,
Int_t cell) const {
// return the path of signal along wire number "cell"
// "ln" - track in sector coor.system
// Path can be <0. or > wire_length !!!
// For GEANT data and don't connected wires return 0.
HMdcSizesCellsCell& fCell=(*this)[cell];
if(!&fCell || fCell.readOutSide=='0') return 0.;
Double_t x=getXSign(ln,cell);
return (fCell.readOutSide=='1') ? fCell.xReadout-x : x-fCell.xReadout;
}
Float_t HMdcSizesCellsLayer::getSignPath(Double_t x1, Double_t y1, Double_t z1,
Double_t x2, Double_t y2, Double_t z2, Int_t cell, Float_t& outside) const {
// return the path of signal along wire number "cell"
// x1,y1,z1 & x2,y2,z2 - track in sector coor.system
// For GEANT data and don't connected wires return 0.
//
// outside:
// if(path<0) outside=-path;
// else if(path>wire_length) outside=path-wire_length;
// else outside=0;
//
// path-outside give new path which >=0. and <=wire_length)
HMdcSizesCellsCell& fCell=(*this)[cell];
if(!&fCell || fCell.readOutSide=='0') {
outside=0.;
return 0.;
}
Float_t x=getXSign(x1, y1,z1, x2, y2, z2,cell);
Float_t path=(fCell.readOutSide=='1') ? fCell.xReadout-x : x-fCell.xReadout;
if(path<0.) outside=path;
else if(path>fCell.length) outside=path-fCell.length;
else outside=0.;
return path;
}
Bool_t HMdcSizesCellsLayer::getDriftDist(Double_t xP1,Double_t yP1,Double_t zP1,
Double_t xP2, Double_t yP2, Double_t zP2,
Int_t cell, Double_t &alphaDrDist, Double_t &driftDist) const {
// ---- esli ponadobitsya, nado proverit' rasschet uglov dlya treka vne cella
//
// return kTRUE if line p1-p2 intersect cell or kFALSE if not intersect
// calc. the angle driftAlpha (in degree) and drift distance in cell
// from line (p1-p2) to wire.
// p1,p2 - in sector coor.sys.
// 1)
// Z^
// | |------------*----|
// | | cell * |
// | | / * |
// | | /90 * |
// | | driftDist/ *|
// | | / *
// --|--|--------*^-----/-|*--------->
// | | / | * Y
// | | alphaDrDist| *
// | | *
// | | *
// |-----------------| * track
//
// 2) *** 90
// Z^ ***
// | |----------:--***-|
// | | cell : / ***
// | | : / | *** track
// | | : / | ***
// | | driftDist/ |
// | | :/ |
// --|--|--------*^-------|---------->
// | | | | Y
// | | | alphaDrDist ???
// | |
// | |
// |-----------------|
//
//
// 3) *
// *
// Z^ / * <-- no eg. 90 deg. !
// | |-----------------/ *
// | | cell /| *
// | | / | *
// | | / | *
// | | driftDist / | *
// | | / | *track
// | | / | *
// | | / | *
// | | / | Alpha*
// --|--|--------*^-------|------------->
// | | | | Y
// | | alphaDrDist|
// | |
// | |
// | |
// | |
// | |
// | |
// |-----------------|
Double_t yWire=cell*pitch-wireOffSet;
Double_t y1,y2,z1,z2;
rotYZtoRotLay(xP1-tRLaySysRSec[9],yP1-tRLaySysRSec[10],zP1-tRLaySysRSec[11],
y1,z1);
y1 -= yWire;
rotYZtoRotLay(xP2-tRLaySysRSec[9],yP2-tRLaySysRSec[10],zP2-tRLaySysRSec[11],
y2,z2);
y2 -= yWire;
Double_t dz=fabs(z1-z2);
Double_t dy=fabs(y1-y2);
Double_t length=sqrt(dz*dz+dy*dy);
Double_t ctgAlpha=dy/dz;
Double_t sinAlpha = dz/length;
Double_t yDist = fabs((z1*y2-z2*y1)/dz);
driftDist = yDist*sinAlpha;
alphaDrDist=atan2(dz,dy)*TMath::RadToDeg();
if(driftDist*sinAlpha > halfPitch) {
Double_t z = (dy != 0.0) ? ((yDist-halfPitch)/ctgAlpha) : 1.E+10;
if(z <= halfCatDist) {
driftDist = sqrt(halfPitch*halfPitch+z*z);
alphaDrDist=atan2((Double_t)halfPitch,z)*TMath::RadToDeg();
} else {
Double_t tgADrDist=halfCatDist/halfPitch;
alphaDrDist=90.-atan(tgADrDist)*TMath::RadToDeg();
driftDist=yDist/(fabs(ctgAlpha*tgADrDist)+1.)*sqrt(1+tgADrDist*tgADrDist);
return kFALSE;
}
} else if(driftDist*dy/length > halfCatDist) { // dy/length = cos(alpha)
Double_t y = yDist-halfCatDist*ctgAlpha;
if(y <= halfPitch) {
driftDist = sqrt(y*y+halfCatDist*halfCatDist);
alphaDrDist=atan2(y,(Double_t)halfCatDist)*TMath::RadToDeg();
} else {
Double_t tgADrDist=halfCatDist/halfPitch;
alphaDrDist=90.-atan(tgADrDist)*TMath::RadToDeg();
driftDist=yDist/(fabs(ctgAlpha*tgADrDist)+1.)*sqrt(1+tgADrDist*tgADrDist);
return kFALSE;
}
}
return kTRUE;
}
void HMdcSizesCellsLayer::print(void) const {
printf("HMSCLayer: %iS %iM %iL, %i cells, cat.dist=%5.2f, pitch=%f, ang.=%4.1fn",
sector+1,module+1,layer+1,nCells,halfCatDist*2,pitch,
asin(sinWireOr)*TMath::RadToDeg());
}
//----------Sector:--------------------------------
HMdcSizesCellsMod::HMdcSizesCellsMod(HMdcSizesCellsSec* pSec, Int_t mod) {
// constructor creates an array of pointers of type
// HMdcSizesCellsLayer
sector=pSec->sector;
module=mod;
pToSec=pSec;
array = new TObjArray(6);
for (Int_t layer = 0; layer < 6; layer++)
(*array)[layer] = new HMdcSizesCellsLayer(this,layer);
}
HMdcSizesCellsMod::~HMdcSizesCellsMod(void) {
// destructor
if(array) {
array->Delete();
delete array;
}
}
Int_t HMdcSizesCellsMod::getSize(void) const {
// return the size of the pointer array
return array->GetEntries();
}
Bool_t HMdcSizesCellsMod::setSecTrans(const HGeomTransform * alignTrans,
Int_t sysOfAlignTrans) {
if(!HMdcGetContainers::getObject()->getSecTransMod(sysRSec,sector,module))
return kFALSE;
if(alignTrans) {
if(sysOfAlignTrans==0) {
sysRSec.transTo(*alignTrans);
HMdcSizesCells::getExObject()->geomModified = kTRUE;
} else if(sysOfAlignTrans==1) {
HGeomTransform alSys=(*alignTrans);
alSys.transFrom(sysRSec);
sysRSec=alSys;
HMdcSizesCells::getExObject()->geomModified = kTRUE;
} else if(sysOfAlignTrans==2) {
sysRSec.transFrom(pToSec->sysRLab);
sysRSec.transTo(*alignTrans);
sysRSec.transTo(pToSec->sysRLab);
HMdcSizesCells::getExObject()->geomModified = kTRUE;
}
}
HMdcSizesCells::copyTransfToArr(sysRSec,tSysRSec);
setPlanePar(sysRSec);
ct[0]=(tSysRSec[0]-tSysRSec[6]*parA);
ct[1]=(tSysRSec[3]-tSysRSec[6]*parB);
ct[2]=tSysRSec[6]*(parD-parB*tSysRSec[10]-tSysRSec[11]-parA*tSysRSec[9]);
ct[3]=(tSysRSec[1]-tSysRSec[7]*parA);
ct[4]=(tSysRSec[4]-tSysRSec[7]*parB);
ct[5]=tSysRSec[7]*(parD-parB*tSysRSec[10]-tSysRSec[11]-parA*tSysRSec[9]);
return kTRUE;
}
void HMdcSizesCellsMod::calcMdcHit(Double_t x1, Double_t y1, Double_t z1,
Double_t x2, Double_t y2, Double_t z2,
Double_t eX1, Double_t eY1, Double_t eZ1, Double_t eX2, Double_t eY2,
Double_t dZ1dX1, Double_t dZ1dY1, Double_t dZ2dX2, Double_t dZ2dY2,
Double_t& x, Double_t& eX, Double_t& y, Double_t& eY,
Double_t& xDir, Double_t& eXDir, Double_t& yDir, Double_t& eYDir) const {
// Calcul. a cross of the line with plane MDC (parA*x+parB*y+c*z=parD),
// transform. this point to mdc coor.sys. (z=0) and calc.
// hit direction in mdc coor.sys. and errors
//
// Input param:
// x1,y1,z1 - target or point on plane
// x2,y2,z2 - point on plane
// eX1,eY1,eZ1, eX2,eY2 - errors (eZ2=0!)
// dZ2dX2 = -A2 where A2 - plane param. (z2=D2-A2*x2+B2*y2)
// dZ2dY2 = -B2 where B2 - plane param. (z2=D2-A2*x2+B2*y2)
// If P1-on plane
// eZ1=0;
// dZ1dX1 = -A1 where A1 - plane param. (z1=D1-A1*x1+B1*y1)
// dZ1dY1 = -B1 where B1 - plane param. (z1=D1-A1*x1+B1*y1)
// If P1-target:
// dZ1dX1 = dZ1/dY1 = 0;
//---Cross-----------------------------------------------
Double_t dX=x2-x1;
Double_t dY=y2-y1;
Double_t dZ=z2-z1;
Double_t cxz1=parD-z1-parA*x1;
Double_t cxz2=parD-z2-parA*x2;
Double_t cyz1=parD-z1-parB*y1;
Double_t cyz2=parD-z2-parB*y2;
Double_t del=1/(parA*dX+parB*dY+dZ);
Double_t xt0=(x2*cyz1 - x1*cyz2)*del;
Double_t yt0=(y2*cxz1 - y1*cxz2)*del;
Double_t xt=xt0-tSysRSec[ 9];
Double_t yt=yt0-tSysRSec[10];
x=ct[0]*xt+ct[1]*yt+ct[2];
y=ct[3]*xt+ct[4]*yt+ct[5];
// - Errors - - - - - - - - - - - - - - - - -
Double_t x1t=x1-xt0;
Double_t y1t=y1-yt0;
Double_t xt2=xt0-x2;
Double_t yt2=yt0-y2;
Double_t xtA=xt0*parA;
Double_t ytB=yt0*parB;
Double_t dXtdX2=x1t*dZ2dX2+cyz1-xtA;
Double_t dYtdY2=y1t*dZ2dY2+cxz1-ytB;
Double_t dXtdX1=xt2*dZ1dX1-cyz2+xtA;
Double_t dYtdY1=yt2*dZ1dY1-cxz2+ytB;
Double_t dXtdY2=x1t*(dZ2dY2+parB);
Double_t dYtdX2=y1t*(dZ2dX2+parA);
Double_t dXtdY1=xt2*(dZ1dY1+parB);
Double_t dYtdX1=yt2*(dZ1dX1+parA);
Double_t dXtdZ1=xt2;
Double_t dYtdZ1=yt2;
Double_t dXdX1=(ct[0]*dXtdX1+ct[1]*dYtdX1)*eX1;
Double_t dYdX1=(ct[3]*dXtdX1+ct[4]*dYtdX1)*eX1;
Double_t dXdY1=(ct[0]*dXtdY1+ct[1]*dYtdY1)*eY1;
Double_t dYdY1=(ct[3]*dXtdY1+ct[4]*dYtdY1)*eY1;
Double_t dXdX2=(ct[0]*dXtdX2+ct[1]*dYtdX2)*eX2;
Double_t dYdX2=(ct[3]*dXtdX2+ct[4]*dYtdX2)*eX2;
Double_t dXdY2=(ct[0]*dXtdY2+ct[1]*dYtdY2)*eY2;
Double_t dYdY2=(ct[3]*dXtdY2+ct[4]*dYtdY2)*eY2;
Double_t dXdZ1=(ct[0]*dXtdZ1+ct[1]*dYtdZ1)*eZ1;
Double_t dYdZ1=(ct[3]*dXtdZ1+ct[4]*dYtdZ1)*eZ1;
eX=sqrt(dXdX1*dXdX1+dXdY1*dXdY1+dXdZ1*dXdZ1+dXdX2*dXdX2+dXdY2*dXdY2)*del;
eY=sqrt(dYdX1*dYdX1+dYdY1*dYdY1+dYdZ1*dYdZ1+dYdX2*dYdX2+dYdY2*dYdY2)*del;
//---Hit direction----------------------------------------------------
Double_t length=1/sqrt(dX*dX+dY*dY+dZ*dZ);
dX *= length; // unit vector
dY *= length;
dZ *= length;
xDir=tSysRSec[0]*dX+tSysRSec[3]*dY+tSysRSec[6]*dZ;
yDir=tSysRSec[1]*dX+tSysRSec[4]*dY+tSysRSec[7]*dZ;
// - Errors - - - - - - - - - - - - - - - - -
Double_t c0xx=(tSysRSec[0]-xDir*dX);
Double_t c3xy=(tSysRSec[3]-xDir*dY);
Double_t c6xz=(tSysRSec[6]-xDir*dZ);
Double_t c1yx=(tSysRSec[1]-yDir*dX);
Double_t c4yy=(tSysRSec[4]-yDir*dY);
Double_t c7yz=(tSysRSec[7]-yDir*dZ);
Double_t dXDdX2 = (c0xx+c6xz*dZ2dX2)*eX2;
Double_t dXDdY2 = (c3xy+c6xz*dZ2dY2)*eY2;
Double_t dXDdX1 = (c0xx+c6xz*dZ1dX1)*eX1; // *(-1)
Double_t dXDdY1 = (c3xy+c6xz*dZ1dY1)*eY1; // *(-1)
Double_t dXDdZ1 = c6xz *eZ1; // *(-1)
Double_t dYDdX2 = (c1yx+c7yz*dZ2dX2)*eX2;
Double_t dYDdY2 = (c4yy+c7yz*dZ2dY2)*eY2;
Double_t dYDdX1 = (c1yx+c7yz*dZ1dX1)*eX1; // *(-1)
Double_t dYDdY1 = (c4yy+c7yz*dZ1dY1)*eY1; // *(-1)
Double_t dYDdZ1 = c7yz *eZ1; // *(-1)
eXDir=sqrt(dXDdX2*dXDdX2+dXDdY2*dXDdY2+
dXDdX1*dXDdX1+dXDdY1*dXDdY1+dXDdZ1*dXDdZ1)*length;
eYDir=sqrt(dYDdX2*dYDdX2+dYDdY2*dYDdY2+
dYDdX1*dYDdX1+dYDdY1*dYDdY1+dYDdZ1*dYDdZ1)*length;
}
//----------Sector:--------------------------------
HMdcSizesCellsSec::HMdcSizesCellsSec(HMdcSizesCells* pSC, Int_t sec) {
// constructor creates an array of pointers of type
// HMdcSizesCellsMod
sector=sec;
array = new TObjArray(4);
mdcStatSec=pSC->mdcStatus[sector];
pToSC=pSC;
for (Int_t mod=0; mod<4; mod++) if(mdcStatSec[mod])
(*array)[mod] = new HMdcSizesCellsMod(this,mod);
targets=0;
}
HMdcSizesCellsSec::~HMdcSizesCellsSec(void) {
// destructor
if(array) {
array->Delete();
delete array;
}
if(targets) delete [] targets;
targets=0;
}
Int_t HMdcSizesCellsSec::getSize(void) const {
// return the size of the pointer array
return array->GetEntries();
}
//==================================================
HMdcSizesCells* HMdcSizesCells::fMdcSizesCells=0;
HMdcSizesCells::HMdcSizesCells(void) {
// constructor creates an array of pointers of type HMdcSizesCellsSec
fGetCont=HMdcGetContainers::getObject();
fLayerGeomPar = fGetCont->getMdcLayerGeomPar();
fGeomPar = fGetCont->getMdcGeomPar();
fSpecGeomPar = fGetCont->getSpecGeomPar();
fMdcDet = fGetCont->getMdcDetector();
fMdcGeomStruct= fGetCont->getMdcGeomStruct();
fMdcRawStruct = fGetCont->getMdcRawStruct();
if (!fMdcRawStruct->isStatic() || !HMdcGetContainers::isInited(fMdcRawStruct))
((HParSet*)fMdcRawStruct)->init(); //Needs explicit initialization
fMdcLookupRaw = (HMdcLookupRaw*)gHades->getRuntimeDb()->
getContainer("MdcLookupRaw");
changed=kFALSE;
geomModified=kFALSE;
array = new TObjArray(6);
for(Int_t i=0;i<3;i++)
verLayerGeomPar[i]=verGeomPar[i]=verLookupRaw[i]=-111111;
targets=0;
numTargets=0;
return;
}
HMdcSizesCells::~HMdcSizesCells(void) {
// destructor
if(array) {
array->Delete();
delete array;
}
HMdcGetContainers::deleteCont();
if(targets) delete [] targets;
targets=0;
}
Bool_t HMdcSizesCells::initContainer(void) {
if( !fMdcDet || !fMdcGeomStruct || !HMdcGetContainers::isInited(fLayerGeomPar)
|| !HMdcGetContainers::isInited(fGeomPar)
|| !HMdcGetContainers::isInited(fMdcLookupRaw)
|| !HMdcGetContainers::isInited(fMdcRawStruct) ) return kFALSE;
if( fLayerGeomPar->hasChanged() || fGeomPar->hasChanged() ||
fMdcLookupRaw->hasChanged()) {
for(Int_t i=1; i<3; i++) {
if(fLayerGeomPar->getInputVersion(i) == verLayerGeomPar[i] &&
fGeomPar->getInputVersion(i) == verGeomPar[i] &&
fMdcLookupRaw->getInputVersion(i) == verLookupRaw[i]) continue;
for(Int_t j=i; j<3; j++) {
verLayerGeomPar[j]=fLayerGeomPar->getInputVersion(j);
verGeomPar[j]=fGeomPar->getInputVersion(j);
verLookupRaw[j]=fMdcLookupRaw->getInputVersion(j);
}
changed=kTRUE;
for (Int_t sec = 0; sec < 6; sec++) {
nModsSeg[sec][0]=nModsSeg[sec][1]=0;
for(Int_t mod=0;mod<4;mod++) {
mdcStatus[sec][mod]=(fMdcDet->getModule(sec,mod)) ? kTRUE:kFALSE;
if(mdcStatus[sec][mod]) nModsSeg[sec][mod/2]++;
}
if( !fMdcDet->isSectorActive(sec) ) continue;
if( !(*array)[sec] ) (*array)[sec]=new HMdcSizesCellsSec(this,sec);
if( !fillCont(sec) ) return kFALSE;
}
fillTargetPos();
break;
}
} else changed=kFALSE;
if(geomModified) changed=kTRUE;
return kTRUE;
}
HMdcSizesCells* HMdcSizesCells::getObject(void) {
if(!fMdcSizesCells) fMdcSizesCells=new HMdcSizesCells();
return fMdcSizesCells;
}
void HMdcSizesCells::deleteCont(void) {
if(fMdcSizesCells) {
delete fMdcSizesCells;
fMdcSizesCells=0;
}
}
HMdcSizesCells* HMdcSizesCells::getExObject(void) {
return fMdcSizesCells;
}
Int_t HMdcSizesCells::getSize(void) const {
// return the size of the pointer array
return array->GetEntries();
}
void HMdcSizesCells::clear(void) {
// clears the container
for(Int_t s=0;s<6;s++) {
HMdcSizesCellsSec* sec=&(*this)[s];
if( !sec ) continue;
for(Int_t m=0;m<4;m++) {
HMdcSizesCellsMod* mod=&(*sec)[m];
if( !mod ) continue;
Int_t nl=mod->getSize();
for(Int_t l=0;l<nl;l++) {
HMdcSizesCellsLayer* lay=&(*mod)[l];
if( !lay ) continue;
Int_t nc=lay->getSize();
for(Int_t c=0;c<nc;c++) (*lay)[c].clear();
}
}
}
changed=kFALSE;
}
Bool_t HMdcSizesCells::fillCont(Int_t sec) {
// calculation of the container HMdcSizesCells
// with wire's sizes (2-pint) in coordinate system of sector.
HMdcSizesCellsSec& fSizesCellsSec=(*this)[sec];
if(!fGetCont->getLabTransSec(fSizesCellsSec.sysRLab,sec)) return kFALSE;
for(Int_t mod=0; mod<4; mod++) if(!fillModCont(sec,mod)) return kFALSE;
return kTRUE;
}
Bool_t HMdcSizesCells::fillModCont(Int_t sec, Int_t mod,
HGeomTransform * alignTrans, Int_t sysOfAlignTrans) {
// filling HMdcSizesCellsMod, HMdcSizesCellsLayer, HMdcSizesCellsCell
// objects for MDC module "mod" in sector "sec".
// if alignTrans>0 - correction of mdc position (must be done after
// filling HMdcSizesCells!)
// sysOfAlignTrans =0 transform mdc coor.sys. by "alignTrans" respect
// coor.sys.of sector
// =1 transform mdc coor.sys. by "alignTrans" respect
// coor.sys.of mdc module
// =2 transform mdc coor.sys. by "alignTrans" respect
// lab.coor.sys.
if(sec<0 || sec>5 || mod<0 || mod>3) return kFALSE;
HMdcSizesCellsSec& fSizesCellsSec=(*this)[sec];
if(&fSizesCellsSec==0) return kFALSE;
if( !mdcStatus[sec][mod] ) return kTRUE;
HMdcSizesCellsMod& fSCellsMod=fSizesCellsSec[mod];
if(!fSCellsMod.setSecTrans(alignTrans,sysOfAlignTrans)) return kFALSE;
for(Int_t layer=0; layer<6; layer++)
if(!fSCellsMod[layer].fillLayerCont()) return kFALSE;
return kTRUE;
}
Bool_t HMdcSizesCellsLayer::getParamCont(void) {
HMdcGetContainers* pGetCont=HMdcGetContainers::getObject();
HGeomCompositeVolume* pCV =pGetCont->getGeomCompositeVolume(module);
if(!pCV) {
Error("setParamCont","GeomCompositeVolume for MDC%i in sec.%i is absent.",
module+1,sector+1);
return kFALSE;
}
pGeomVol=pCV->getComponent(layer);
if(pGeomVol==0) {
Error("setParamCont",
"HGeomVolume for lay. %i in sec.%i mod.%i doesn't exist!",layer+1,
sector+1,module+1);
return kFALSE;
}
pLayerGeomParLay=&(*(pGetCont->getMdcLayerGeomPar()))[sector][module][layer];
if(pLayerGeomParLay==0) {
Error("setParamCont",
"HLayerGeomParLay for lay. %i in sec.%i mod.%i doesn't exist!",layer+1,
sector+1,module+1);
return kFALSE;
}
return kTRUE;
}
Bool_t HMdcSizesCellsLayer::fillLayerCont(const Double_t* corr) {
if(!getParamCont()) return kFALSE;
if(nCells<0) return kTRUE;
pitch=pLayerGeomParLay->getPitch();
wireOffSet=(pLayerGeomParLay->getCentralWireNr()-1.)*pitch;
halfPitch=pitch*0.5;
halfCatDist=pLayerGeomParLay->getCatDist()*0.5;
Double_t wOr=pLayerGeomParLay->getWireOrient();
// Layer geometry correction:
if(corr) {
wOr += corr[layer];
wireOffSet += corr[layer + 6];
}
wOr *= TMath::DegToRad();
cosWireOr=cos(wOr);
sinWireOr=sin(wOr);
if(!setSecTrans(corr ? corr[layer+12]:0.)) return kFALSE;
Int_t nPoint=pGeomVol->getNumPoints();
if(nPoint != 8) {
Error("fillCont","NumPoints for layer volume not eg. 8.");
return kFALSE;
}
Double_t x[4],y[4];
for(nPoint=0; nPoint<4;nPoint++) {
const HGeomVector* pnt=pGeomVol->getPoint(nPoint); // mm!
x[nPoint]=pnt->getX();
y[nPoint]=pnt->getY();
}
Double_t a=tan(wOr); // y=a*x+ln[1]
// ----- calcBLines ----------------
Double_t at[4],bt[4];
for(Int_t i1=0; i1<4; i1++) {
Int_t i2=i1+1;
if(i2==4) i2=0;
at[i1]=(y[i2]-y[i1])/(x[i2]-x[i1]);
bt[i1]=y[i1]-at[i1]*x[i1];
}
//----------------------------------
Double_t ymax=-1.e+30;
Double_t ymin=1.e+30;
for(Int_t i=0; i<4; i++) {
Double_t yp=y[i]*cosWireOr-x[i]*sinWireOr;
if(yp>ymax) ymax=yp;
if(yp<ymin) ymin=yp;
}
for(Int_t cell=0; cell<nCells; cell++) {
Double_t b=(cell*pitch-wireOffSet)/cosWireOr;
HMdcSizesCellsCell& fSizesCellsCell=(*this)[cell];
HGeomVector& lineP1=fSizesCellsCell.wirePnt1;
HGeomVector& lineP2=fSizesCellsCell.wirePnt2;
//---------- Calculation of wire size ------------------
//
//
// Line 1
// ---------- ^ (y)
// Layre / |
// Line 0 -> MDC / <- Line 2 |
// ____/ |
// Line 3 |
// (x) <--------|----------------------
// 0,0
//
Double_t x1=(bt[0]-b)/(a-at[0]); //Xb line 0
Double_t x2=(bt[2]-b)/(a-at[2]); //Xe line 2
Double_t y1=a*x1+b; // Yb
Double_t y2=a*x2+b; // Ye
Int_t nLine1=(x1>=x[0] && x1<=x[1] && y1>=y[0] && y1<=y[1]) ? 0:-1;
Int_t nLine2=(x2<=x[3] && x2>=x[2] && y2>=y[3] && y2<=y[2]) ? 2:-1;
if(nLine1<0 && nLine2 <0) {
if(cell<100) {
Warning("fillCont","S%i M%i L%i: Cell %i is out of MDC.",
sector+1,module+1,layer+1,cell+1);
continue;
} else {
Warning("fillCont","S%i M%i L%i: Cells %i-%i are out of MDC.",
sector+1,module+1,layer+1,cell+1,nCells);
break;
}
}
if(nLine1<0) {
x1=(bt[1]-b)/(a-at[1]); // Xb line 1
if(x1<x[2]) x1=(bt[3]-b)/(a-at[3]); // Xb line 3
y1=a*x1+b; // Yb
} else if(nLine2<0) {
x2=(bt[1]-b)/(a-at[1]); // Xe line 1
if(x2>x[1]) x2=(bt[3]-b)/(a-at[3]); // Xe line 3
y2=a*x2+b; // Ye
}
HMdcSizesCells* pSC=HMdcSizesCells::getExObject();
Char_t readOutSide=pSC->findReadOutSide(sector,module,layer,cell);
if( readOutSide=='0'){
if( cell>3 && cell<(((*(pSC->fMdcGeomStruct))[sector][module][layer])-4)) {
if(corr==0) Warning("fillCont",
" %iS%iM%iL: Wire %i is not connected.",sector+1,module+1,layer+1,cell+1);
}
} else {
if( (readOutSide=='1' && x1<x2) || (readOutSide=='2' && x1>x2) ||
(readOutSide=='3' && y1<y2) ) {
Double_t tmp;
tmp=x1; x1=x2; x2=tmp;
tmp=y1; y1=y2; y2=tmp;
}
if(readOutSide=='3') readOutSide=(x1>x2) ? '1':'2';
}
fSizesCellsCell.readOutSide=readOutSide;
fSizesCellsCell.xReadout=x1*cosWireOr + y1*sinWireOr;
fSizesCellsCell.length=fabs(x2*cosWireOr + y2*sinWireOr-
fSizesCellsCell.xReadout);
lineP1.setXYZ(x1,y1,0.);
lineP2.setXYZ(x2,y2,0.);
lineP1=sysRSec.transFrom(lineP1);
lineP2=sysRSec.transFrom(lineP2);
fSizesCellsCell.setStatus(kTRUE);
}
return kTRUE;
}
Bool_t HMdcSizesCells::fillModCont(Int_t sec, Int_t mod, Double_t * corr) {
// filling HMdcSizesCellsMod, HMdcSizesCellsLayer, HMdcSizesCellsCell
// objects for MDC module "mod" in sector "sec"
if(sec<0 || sec>5 || mod<0 || mod>3) return kFALSE;
HMdcSizesCellsSec& fSizesCellsSec=(*this)[sec];
if(&fSizesCellsSec==0) return kFALSE;
if( !mdcStatus[sec][mod] ) return kTRUE;
HMdcSizesCellsMod& fSCellsMod=fSizesCellsSec[mod];
if(!fSCellsMod.setSecTrans()) return kFALSE;
for(Int_t layer=0; layer<6; layer++)
if(!fSCellsMod[layer].fillLayerCont(corr)) return kFALSE;
return kTRUE;
}
Bool_t HMdcSizesCells::getCellVol(Int_t sec,Int_t mod,Int_t lay,Int_t cell,
HMdcTrap& volCell, Double_t sizeFactor) const {
// Return the cell volume as HMdcTrap (8 points) in sector(!) coordinate sys..
// Use "sizeFactor" param. for drawing purpose only.
// This parameter increase Y and Z cell size on this factor.
HGeomVector layerVol[8]; // Drift gas vol. (8-points)
if( &((*this)[sec]) == 0 ) return kFALSE;
HMdcSizesCellsMod& fSCellsMod=(*this)[sec][mod];
if(!&fSCellsMod) return kFALSE;
HGeomCompositeVolume* fComVol=fGetCont->getGeomCompositeVolume(mod);
if(!fComVol) return kFALSE;
HMdcSizesCellsLayer& fSCellsLay=fSCellsMod[lay];
if(cell<0 || cell>=fSCellsLay.getSize() || !fSCellsLay[cell].getStatus())
return kFALSE;
HGeomVolume* fGeomVol=fComVol->getComponent(lay);
for(Int_t point=0; point<4; point++) {
layerVol[point]=*(fGeomVol->getPoint(point)); // mm!
layerVol[point+4]=layerVol[point];
layerVol[point].setZ(-fSCellsLay.halfCatDist);
layerVol[point+4].setZ(+fSCellsLay.halfCatDist);
}
Double_t a,b;
a=fSCellsLay.sinWireOr/fSCellsLay.cosWireOr; // y=a*x+ln[1]
// ----- calcBLines ----------------
Double_t at[4],bt[4];
for(Int_t i1=0; i1<4; i1++) {
Int_t i2=i1+1;
if(i2==4) i2=0;
at[i1]=(layerVol[i2].getY()-layerVol[i1].getY())/
(layerVol[i2].getX()-layerVol[i1].getX());
bt[i1]=layerVol[i1].getY()-at[i1]*layerVol[i1].getX();
}
Double_t ymax, ymin;
ymax=-1.e+30;
ymin=1.e+30;
for(Int_t i=0; i<4; i++) {
Double_t y=layerVol[i].getY()*fSCellsLay.cosWireOr-
layerVol[i].getX()*fSCellsLay.sinWireOr;
if(y>ymax) ymax=y;
if(y<ymin) ymin=y;
}
b=(cell*fSCellsLay.pitch-fSCellsLay.wireOffSet)/fSCellsLay.cosWireOr;
//---------- Calculation of wire size ------------------
//
//
// Line 1
// ---------- ^ (y)
// Layre / |
// Line 0 -> MDC / <- Line 2 |
// ____/ |
// Line 3 |
// (x) <--------|----------------------
// 0,0
//
Double_t x1=(bt[0]-b)/(a-at[0]); //Xb line 0
Double_t x2=(bt[2]-b)/(a-at[2]); //Xe line 2
Double_t y1=a*x1+b; // Yb
Double_t y2=a*x2+b; // Ye
Int_t nLine1=(x1>=layerVol[0].getX() && x1<=layerVol[1].getX() &&
y1>=layerVol[0].getY() && y1<=layerVol[1].getY()) ? 0:-1;
Int_t nLine2=(x2<=layerVol[3].getX() && x2>=layerVol[2].getX() &&
y2>=layerVol[3].getY() && y2<=layerVol[2].getY()) ? 2:-1;
if(nLine1<0 && nLine2 <0) return kFALSE;
if(nLine1<0) nLine1=((bt[1]-b)/(a-at[1]) >= layerVol[2].getX()) ? 1:3;
else if(nLine2<0) nLine2=((bt[1]-b)/(a-at[1]) <= layerVol[1].getX()) ? 1:3;
Double_t tLine,bLine;
Double_t halfPitch=fSCellsLay.halfPitch;
if(sizeFactor>0.) halfPitch*=sizeFactor;
bLine=cell*fSCellsLay.pitch-fSCellsLay.wireOffSet-halfPitch;
tLine=cell*fSCellsLay.pitch-fSCellsLay.wireOffSet+halfPitch;
if(tLine > ymax) tLine=ymax;
if(bLine < ymin) bLine=ymin;
bLine=bLine/fSCellsLay.cosWireOr;
tLine=tLine/fSCellsLay.cosWireOr;
Double_t x,z;
z=layerVol[0](2);
if(sizeFactor>0.) z*=sizeFactor;
x=(bt[nLine1]-bLine)/(a-at[nLine1]);
volCell[0].setXYZ(x,a*x+bLine,z);
x=(bt[nLine1]-tLine)/(a-at[nLine1]);
volCell[1].setXYZ(x,a*x+tLine,z);
x=(bt[nLine2]-tLine)/(a-at[nLine2]);
volCell[2].setXYZ(x,a*x+tLine,z);
x=(bt[nLine2]-bLine)/(a-at[nLine2]);
volCell[3].setXYZ(x,a*x+bLine,z);
z=layerVol[4](2);
if(sizeFactor>0.) z*=sizeFactor;
for(Int_t i=0; i<4; i++) {
volCell[i+4]=volCell[i];
volCell[i+4].setZ(z);
}
volCell.transFrom(fSCellsLay.sysRSec);
return kTRUE;
}
Char_t HMdcSizesCells::findReadOutSide(Int_t s,Int_t m,Int_t l,Int_t c) const {
// Function to get the side from which the corresponding wire is read out.
// Function return values are '1' for left, '2' for right, '3' for top
// and '0' for GEANT data and not connected wires.
if(!fMdcLookupRaw || !fMdcRawStruct) return '0';
HMdcLookupCell &cell=(*fMdcLookupRaw)[s][m][l][c];
if(!&cell) return '0';
Int_t mbo=cell.getNMoth();
if(mbo < 0) return '0'; // don't conected wires have mbo==-1
HMdcRawMothStru& rMbo=(*fMdcRawStruct)[s][m][mbo];
if(!&rMbo) return '0';
Char_t mbn=rMbo.GetName()[0];
return (mbn=='1' || mbn=='2' || mbn=='3') ? mbn:'0';
}
Bool_t HMdcSizesCells::fillTargetPos(HGeomVector* targetShift) {
// Filling target parameters
if(targets && numTargets<fSpecGeomPar->getNumTargets()) {
delete [] targets;
targets=0;
}
numTargets=fSpecGeomPar->getNumTargets();
if( numTargets < 1 ) {
Error("calcTarget","Number of targets = %i!",numTargets);
return kFALSE;
}
if(targets==0) targets=new HGeomVector [numTargets];
HGeomVolume* tr1=0;
HGeomVolume* tr2=0;
for(Int_t t=0;t<numTargets;t++) {
HGeomVolume* tr=fSpecGeomPar->getTarget(t);
if(!tr) {
Error("calcTarget","GeomVolume for target absent!");
return kFALSE;
}
targets[t]=tr->getTransform().getTransVector();
if(tr1==0) tr1=tr2=tr;
else if( tr ->getTransform().getTransVector().getZ() <
tr1->getTransform().getTransVector().getZ() ) tr1=tr;
else if( tr ->getTransform().getTransVector().getZ() >
tr2->getTransform().getTransVector().getZ() ) tr2=tr;
}
targ3p[0]=tr1->getTransform().getTransVector();
Double_t dz1=tr1->getPoint(0)->getZ();
if(dz1 > tr1->getPoint(2)->getZ()) dz1=tr1->getPoint(2)->getZ();
targ3p[0].setZ(targ3p[0].getZ() + dz1);
targ3p[2]=tr2->getTransform().getTransVector();
Double_t dz2=tr2->getPoint(2)->getZ();
if(dz2 < tr2->getPoint(0)->getZ()) dz2=tr2->getPoint(0)->getZ();
targ3p[2].setZ(targ3p[2].getZ() + dz2);
targ3p[1]=targ3p[0]+targ3p[2];
targ3p[1] *= 0.5;
if(targetShift) {
for(Int_t i=0;i<3;i++) targ3p[i] += *targetShift;
for(Int_t i=0;i<numTargets;i++) targets[i] += *targetShift;
}
for(Int_t sec=0;sec<6;sec++) {
if((*array)[sec] == 0) continue;
HMdcSizesCellsSec* fSec=(HMdcSizesCellsSec*)((*array)[sec]);
const HGeomTransform* trans=fSec->getLabTrans();
for(Int_t i=0;i<3;i++) fSec->targ3p[i]=trans->transTo(targ3p[i]);
if(fSec->targets && fSec->numTargets<numTargets) {
delete [] fSec->targets;
fSec->targets=0;
}
if(fSec->targets==0) fSec->targets=new HGeomVector [numTargets];
for(Int_t i=0;i<numTargets;i++) fSec->targets[i]=
trans->transTo(targets[i]);
fSec->numTargets=numTargets;
}
return kTRUE;
}
void HMdcSizesCells::calcMdcSeg(Double_t x1, Double_t y1, Double_t z1,
Double_t x2, Double_t y2, Double_t z2,
Double_t eX1, Double_t eY1, Double_t eZ1, Double_t eX2, Double_t eY2,
Double_t dZ1dX1, Double_t dZ1dY1, Double_t dZ2dX2, Double_t dZ2dY2,
Double_t& zm, Double_t& eZm, Double_t& r0, Double_t& eR0,
Double_t& theta, Double_t& eTheta, Double_t& phi, Double_t& ePhi) {
// Calcul. a HMdcSeg hit.
// Input and output are in sector coor.sys.
// theta=atan(sqrt(dX*dX+dY*dY)/dZ); phi=atan(dY/dX)
// zm= z1 - cos(theta)/sin(theta) * (x1*cos(phi)+y1*sin(phi))
// r0=y1*cos(phi)-x1*sin(phi)
//
// If (x1,y1,z1)=(x2,y2,z2) dZ will eq.1.!
// Input param:
// x1,y1,z1 - target or point on plane
// x2,y2,z2 - point on plane
// eX1,eY1,eZ1, eX2,eY2 - errors (eZ2=0!)
// dZ2dX2 = -A2 where A2 - plane param. (z2=D2-A2*x2+B2*y2)
// dZ2dY2 = -B2 where B2 - plane param. (z2=D2-A2*x2+B2*y2)
// If P1-on plane
// eZ1=0;
// dZ1dX1 = -A1 where A1 - plane param. (z1=D1-A1*x1+B1*y1)
// dZ1dY1 = -B1 where B1 - plane param. (z1=D1-A1*x1+B1*y1)
// If P1-target:
// dZ1dX1 = dZ1/dY1 = 0;
Double_t dX=x2-x1;
Double_t dY=y2-y1;
Double_t dZ=z2-z1;
Double_t lenXY=sqrt(dX*dX+dY*dY);
if(lenXY<2.E-5) {
dX =x2 * 1.E-5/sqrt(x2*x2+y2*y2);
dY =y2 * 1.E-5/sqrt(x2*x2+y2*y2);
lenXY=1.E-5; // = dX*dX+dY*dY;
dZ=1.;
}
dX /= lenXY;
dY /= lenXY;
dZ /= lenXY;
phi=atan2(dY,dX);
theta=atan2(1.,dZ);
Double_t xy2=(x2*dX+y2*dY);
Double_t xy1=(x1*dX+y1*dY);
zm=(z1*xy2-z2*xy1)/lenXY; //= z1-dZ*(x1*dX+y1*dY);
r0=y1*dX-x1*dY;
Double_t dR0dX1=(y2-r0*dX)*eX1; // *(-1)
Double_t dR0dX2=(y1-r0*dX)*eX2;
Double_t dR0dY1=(x2+r0*dY)*eY1;
Double_t dR0dY2=(x1+r0*dY)*eY2; // *(-1)
eR0=sqrt(dR0dX1*dR0dX1+dR0dY1*dR0dY1+dR0dX2*dR0dX2+dR0dY2*dR0dY2)/lenXY;
Double_t zx12=(z1*dX-x1*dZ)-2*zm*dX;
Double_t zy12=(z1*dY-y1*dZ)-2*zm*dY;
Double_t dZmdZ1= xy2*eZ1;
Double_t dZmdX1=(xy2*dZ1dX1-zx12-z2*dX)*eX1;
Double_t dZmdX2=(xy1*dZ2dX2-zx12-z1*dX)*eX2; // *(-1)
Double_t dZmdY1=(xy2*dZ1dY1-zy12-z2*dY)*eY1;
Double_t dZmdY2=(xy1*dZ2dY2-zy12-z1*dY)*eY2; // *(-1)
eZm=sqrt(dZmdZ1*dZmdZ1+dZmdX1*dZmdX1+dZmdX2*dZmdX2+dZmdY1*dZmdY1+
dZmdY2*dZmdY2)/lenXY;
Double_t drvTh=1./(1.+dZ*dZ);
Double_t dThdX1=(dZ*dX-dZ1dX1)*eX1; // *(-1)
Double_t dThdY1=(dZ*dY-dZ1dY1)*eY1; // *(-1)
Double_t dThdX2=(dZ*dX-dZ2dX2)*eX2;
Double_t dThdY2=(dZ*dY-dZ2dY2)*eY2;
// Double_t dThdZ1=eZ1;
eTheta=sqrt(dThdX1*dThdX1+dThdY1*dThdY1+eZ1*eZ1+dThdX2*dThdX2+dThdY2*dThdY2)*
drvTh/lenXY;
Double_t dPhdX1=dY*eX1;
Double_t dPhdX2=dY*eX2; // *(-1)
Double_t dPhdY1=dX*eY1; // *(-1)
Double_t dPhdY2=dX*eY2;
ePhi=sqrt(dPhdX1*dPhdX1+dPhdX2*dPhdX2+dPhdY1*dPhdY1+dPhdY2*dPhdY2)/lenXY;
}
void HMdcSizesCells::calcRZToTargLine(Double_t x1, Double_t y1, Double_t z1,
Double_t x2, Double_t y2, Double_t z2,
Double_t &zm, Double_t &r0) const {
// Input and output are in lab.sys.
calcRZToLineXY(x1,y1,z1,x2,y2,z2,targ3p[2].getX(),targ3p[2].getY(),zm,r0);
}
void HMdcSizesCellsSec::calcRZToTargLine(Double_t x1, Double_t y1, Double_t z1,
Double_t x2, Double_t y2, Double_t z2,
Double_t &zm, Double_t &r0) const {
// Input and output are in sector coor.sys.
HMdcSizesCells::calcRZToLineXY(x1,y1,z1,x2,y2,z2,
targ3p[2].getX(),targ3p[2].getY(),zm,r0);
}
void HMdcSizesCells::calcRZToLineXY(Double_t x1, Double_t y1, Double_t z1,
Double_t x2, Double_t y2, Double_t z2,
Double_t xBeam, Double_t yBeam, Double_t &zm, Double_t &r0) {
// zm= z1 - cos(theta)/sin(theta) * ((x1-xBeam)*cos(phi)+(y1-yBeam)*sin(phi))
// r0=(y1-yBeam)*cos(phi)-(x1-xBeam)*sin(phi)
//
// If (x1,y1,z1)=(x2,y2,z2) dZ will eq.1.!
Double_t dX=x2-x1;
Double_t dY=y2-y1;
Double_t dZ=z2-z1;
Double_t lenXY2=dX*dX+dY*dY;
Double_t lenXY=sqrt(lenXY2);
if(lenXY<2.E-5) {
dX =x2 * 1.E-5/sqrt(x2*x2+y2*y2);
dY =y2 * 1.E-5/sqrt(x2*x2+y2*y2);
lenXY=1.E-5; //dX*dX+dY*dY;
lenXY2=lenXY*lenXY;
dZ=1.;
}
x1-=xBeam; // Shifting to beam line
y1-=yBeam; // Shifting to beam line
zm=z1-dZ*(x1*dX+y1*dY)/lenXY2;
r0=(y1*dX-x1*dY)/lenXY;
}
void HMdcSizesCells::setTransform(Double_t* par, HGeomTransform& trans) {
// Filling HGeomTransform object by 6 parameters
// par[0] - shift along X
// par[1] - shift along Y
// par[2] - shift along Z
// par[3] - rotation around X (in degree)
// par[4] - rotation around Y (in degree)
// par[5] - rotation around Z (in degree)
if(par==0) return;
Double_t alpha=par[3]*TMath::DegToRad();
Double_t beta=par[4]*TMath::DegToRad();
Double_t gamma=par[5]*TMath::DegToRad();
Double_t cosAlpha=TMath::Cos(alpha);
Double_t sinAlpha=TMath::Sin(alpha);
Double_t cosBeta=TMath::Cos(beta);
Double_t sinBeta=TMath::Sin(beta);
Double_t cosGamma=TMath::Cos(gamma);
Double_t sinGamma=TMath::Sin(gamma);
Double_t rot[9];
rot[0]=cosBeta*cosGamma;
rot[1]=-cosBeta*sinGamma;
rot[2]=sinBeta;
rot[3]=sinAlpha*sinBeta*cosGamma+cosAlpha*sinGamma;
rot[4]=-sinAlpha*sinBeta*sinGamma+cosAlpha*cosGamma;
rot[5]=-sinAlpha*cosBeta;
rot[6]=-cosAlpha*sinBeta*cosGamma+sinAlpha*sinGamma;
rot[7]=cosAlpha*sinBeta*sinGamma+sinAlpha*cosGamma;
rot[8]=cosAlpha*cosBeta;
trans.setRotMatrix(rot);
trans.setTransVector(par);
}
void HMdcSizesCells::copyTransfToArr(const HGeomTransform& trans,
Double_t* arr) {
// Copy trans to array
// Array size must be eq. or more 12!
const HGeomRotation& rot=trans.getRotMatrix();
const HGeomVector& trn=trans.getTransVector();
for(Int_t i=0; i<9; i++) arr[i]=rot(i);
arr[ 9]=trn(0);
arr[10]=trn(1);
arr[11]=trn(2);
}
ROOT page - Class index - Class Hierarchy - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.