#include "hgeommdc.h"
#include "hgeommdchit.h"
#include "hgeombuilder.h"
#include "hgeommedia.h"
#include "hgeommedium.h"
#include "hgeomnode.h"
#include "hgeommdcwireplanes.h"
#include "hgeommdcwire.h"
#include "TFile.h"
#include <iostream>
#include <iomanip>
using namespace std;
ClassImp(HGeomMdc)
HGeomMdc::HGeomMdc() {
fName="mdc";
maxSectors=6;
maxModules=4;
pHit=new HGeomMdcHit(this);
wn0=33; wn1=0; wn2=0; wn3=-1;
memcpy(wnbuf,"0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ",36);
}
const Char_t* HGeomMdc::getModuleName(Int_t m) {
sprintf(modName,"DR%iM",m+1);
return modName;
}
const Char_t* HGeomMdc::getEleName(Int_t m) {
sprintf(eleName,"D%i",m+1);
return eleName;
}
Bool_t HGeomMdc::createAdditionalGeometry(HGeomBuilder* builder, const TString& paramFile, HGeomMedia* media) {
Bool_t rc = kTRUE;
if (builder && paramFile.Length() > 0) {
TFile* file = new TFile(paramFile.Data(),"READ");
if (file==NULL || file->IsOpen()== kFALSE) {
Error("createAdditionalGeometry", "Parameter file %s not found", paramFile.Data());
return kFALSE;
}
file->cd();
HGeomMdcWirePlanes *mdcWirePlanes = (HGeomMdcWirePlanes*)file->Get("HGeomMdcWirePlanes");
if (mdcWirePlanes) {
vector<HGeomMdcWirePlane>& pWirePlanes = mdcWirePlanes->getWirePlanes();
Int_t totNumWires=0;
for (UInt_t np=0; np<pWirePlanes.size(); np++) {
HGeomMdcWirePlane& plane = pWirePlanes[np];
HGeomNode* pMother=getVolume(plane.planeName.Data());
if (pMother==NULL || pMother->isActive()==kFALSE) continue;
Double_t planeX[4], planeY[4];
for (Int_t i=0; i<4; i++) {
HGeomVector* point = pMother->getPoint(i);
planeX[i] = point->getX();
planeY[i] = point->getY();
}
Int_t medId[2] = {-1,-1};
for (Int_t i = 0; i <= plane.planeType && rc; i++) {
HGeomMedium* medium = media->getMedium(plane.wireMedium[i]);
if (medium) {
medId[i] = medium->getMediumIndex();
if (medId[i] <= 0) medId[i] = builder->createMedium(medium);
if (medId[i] <= 0) {
Error("createAdditionalGeometry",
"Medium %s not created",plane.wireMedium[i].Data());
rc = kFALSE;
}
} else {
Error("createAdditionalGeometry",
"Medium %s not found in list of media",plane.wireMedium[i].Data());
rc = kFALSE;
}
}
if (rc && plane.planeType == 0) {
Float_t wireDist = plane.wireDist;
Float_t halfDist = wireDist * 0.5;
Float_t radius = plane.wireRadius[0];
Double_t posX = 0., posY = 0., posZ = 0.;
Double_t at = (planeY[1] - planeY[0]) / (planeX[1] - planeX[0]);
Double_t bt = planeY[0] - at * planeX[0];
Double_t halfLengthCenter = (planeY[1] - planeY[0]) * 0.5;
TString wireName;
generateWireName(wireName);
Int_t wireNum = 0;
Int_t copyNum = 1;
HGeomMdcWire* pWire = NULL;
Int_t arrIndex = -1;
do {
copyNum = wireNum + 1;
arrIndex = addWireObject(wireNum, wireName, copyNum, 0, radius, halfLengthCenter, posX, posY, posZ);
if (wireNum == 0) {
pWire = wireObjects[arrIndex];
} else {
wireObjects[arrIndex]->setCopyNode(pWire);
}
wireNum++;
posX += wireDist;
} while ((posX+radius) <= planeX[0]);
Int_t maxCopyNumLeft = copyNum;
copyNum = 1;
do {
generateWireName(wireName);
Double_t y = at * (posX + radius) + bt;
Double_t halfLength = (planeY[1] - y) / 2.;
posZ = halfLengthCenter - halfLength;
addWireObject(wireNum, wireName, copyNum, 0, radius, halfLength, posX, posY, posZ);
wireNum++;
posX += wireDist;
} while ((posX+halfDist) < planeX[1]);
Int_t maxInd = wireObjects.size();
Int_t rotInd = 0;
for (Int_t i=0; i<maxInd && rc; i++) {
pWire = wireObjects[i];
copyNum = pWire->getCopyNumber();
if (copyNum == 1) {
rc = builder->createVolume(pWire,medId[0]);
if (!rc) {
Error("createAdditionalGeometry",
"Plane %s wire %i not created",plane.planeName.Data(),i);
break;
}
}
rc = builder->positionNode(pWire,pMother,rotInd);
if (!rc) {
Error("createAdditionalGeometry",
"Plane %s wire %i not positioned",plane.planeName.Data(),i);
}
}
for (Int_t i=1; i<maxInd && rc; i++) {
pWire = wireObjects[i];
copyNum = pWire->getCopyNumber();
if (copyNum > 1) {
copyNum += maxCopyNumLeft -1;
} else {
copyNum += 1;
}
Float_t* wirePos = pWire->getPosition();
wirePos[0] *= -1.;
pWire->setWireNumber(wireNum);
pWire->setCopyNumber(copyNum);
rc = builder->positionNode(pWire,pMother,rotInd);
wireNum++;
}
totNumWires += wireNum;
clearWireObjects();
}
if (rc && plane.planeType == 1) {
Double_t wireOrient = -plane.wireOrient;
Double_t sinWireOr = TMath::Sin(wireOrient*TMath::DegToRad());
Double_t cosWireOr = TMath::Cos(wireOrient*TMath::DegToRad());
Double_t arr[] = {sinWireOr, 0., cosWireOr, 0., 1., 0., -cosWireOr, 0., sinWireOr};
HGeomRotation rotMatrix;
rotMatrix.setMatrix(arr);
Int_t rotInd = builder->createRotation(&rotMatrix);
Double_t planeCenter[2];
planeCenter[0] = (planeX[0] + planeX[1] + planeX[2] + planeX[3]) * .25;
planeCenter[1] = (planeY[0] + planeY[1] + planeY[2] + planeY[3]) * .25;
for (Int_t i=0; i<4; i++) {
planeX[i] -= planeCenter[0];
planeY[i] -= planeCenter[1];
}
Double_t wireOffset = (plane.centralWireNr - 1) * plane.wireDist;
Double_t at[4], bt[4];
Double_t radius = plane.wireRadius[0];
at[0] = (planeY[1] - planeY[0]) / (planeX[1] - planeX[0]);
Double_t dy = radius * cosWireOr;
Double_t dx1 = dy / at[0];
Double_t dx2 = 0;
if (wireOrient > 0.) {
dx2 = radius * (sinWireOr + cosWireOr / at[0]);
} else {
dx2 = -radius * (sinWireOr - cosWireOr / at[0]);
}
planeX[0] = planeX[0] + dx1 - dx2;
planeY[0] = planeY[0] + dy;
planeX[1] = planeX[1] - dx1 - dx2;
planeY[1] = planeY[1] - dy;
planeX[2] = -planeX[1];
planeY[2] = planeY[1];
planeX[3] = -planeX[0];
planeY[3] = planeY[0];
for (Int_t i1=0; i1<4; i1++) {
Int_t i2 = i1 + 1;
if (i2==4) i2=0;
at[i1] = (planeY[i2] - planeY[i1]) / (planeX[i2] - planeX[i1]);
bt[i1] = planeY[i1] - at[i1] * planeX[i1];
}
Double_t a = TMath::Tan(wireOrient*TMath::DegToRad());
Int_t numWiresPlane=0;
for(Int_t wireNum=0; wireNum<plane.numWires; wireNum++) {
Int_t wireType = wireNum%2;
Float_t radius = plane.wireRadius[wireType];
Double_t yWire = wireNum * plane.wireDist - wireOffset;
Double_t b = yWire/cosWireOr - planeCenter[1];
Double_t x1 = (bt[0]-b)/(a-at[0]);
Double_t x2 = (bt[2]-b)/(a-at[2]);
Double_t y1 = a*x1+b;
Double_t y2 = a*x2+b;
Int_t nLine1 = (x1>=planeX[0] && x1<=planeX[1] && y1>=planeY[0] && y1<=planeY[1]) ? 0:-1;
Int_t nLine2 = (x2<=planeX[3] && x2>=planeX[2] && y2>=planeY[3] && y2<=planeY[2]) ? 2:-1;
if(nLine1<0 && nLine2 <0) {
continue;
}
if(nLine1<0) {
x1 = (bt[1]-b)/(a-at[1]);
if(x1<planeX[2]) x1 = (bt[3]-b)/(a-at[3]);
y1 = a*x1+b;
} else if(nLine2<0) {
x2 = (bt[1]-b)/(a-at[1]);
if(x2>planeX[1]) x2 = (bt[3]-b)/(a-at[3]);
y2 = a*x2+b;
}
Double_t wireLength = TMath::Sqrt((x2-x1)*(x2-x1)+(y2-y1)*(y2-y1));
Double_t posX = (x2 + x1) * .5;
Double_t posY = 0.;
Double_t posZ = (y2 + y1) * .5;
TString wireName;
generateWireName(wireName);
HGeomMdcWire wire(wireNum, wireName, 1, wireType, radius, wireLength * 0.5, posX, posY, posZ);
rc = builder->createVolume(&wire, medId[wireType]);
if (!rc) {
Error("createAdditionalGeometry",
"Plane %s wire %i not created",plane.planeName.Data(),wireNum);
break;
}
rc = builder->positionNode(&wire,pMother,rotInd);
if (!rc) {
Error("createAdditionalGeometry",
"Plane %s wire %i not positioned",plane.planeName.Data(),wireNum);
}
numWiresPlane++;
}
totNumWires += numWiresPlane;
}
if (rc) builder->fillMdcCommonBlock(pMother);
}
cout<<" Wires created: "<<totNumWires<<endl;
delete mdcWirePlanes;
mdcWirePlanes = NULL;
} else {
Warning("createAdditionalGeometry", "HGeomMdcWires not found in parameter file %s", paramFile.Data());
}
file->Close();
}
return rc;
}
void HGeomMdc::generateWireName(TString& vName) {
if (wn3<35) wn3++;
else {
wn3 = 0;
if (wn2<35) wn2++;
else {
wn2 = 0;
if (wn1<35) wn1++;
else {
wn1 = 0;
wn0++;
}
}
}
vName.Form("%c%c%c%c",wnbuf[wn0],wnbuf[wn1],wnbuf[wn2],wnbuf[wn3]);
}
void HGeomMdc::clearWireObjects(){
for(UInt_t i = 0; i < wireObjects.size(); i ++) {
if (wireObjects[i]) {
delete wireObjects[i];
wireObjects[i] = 0;
}
}
wireObjects.clear();
}
Int_t HGeomMdc::addWireObject(Int_t wn, TString& wname, Int_t cn, Int_t wtype, Float_t radius,
Double_t hlen, Double_t xpos, Double_t ypos, Double_t zpos) {
HGeomMdcWire* wire = new HGeomMdcWire(wn, wname, cn, wtype, radius, hlen, xpos, ypos, zpos);
wireObjects.push_back(wire);
return (wireObjects.size()-1);
}