#include <fstream>
#include <TClonesArray.h>
#include <TGeoManager.h>
#include <TString.h>
#include <TMath.h>
#include "AliSurveyObj.h"
#include "AliSurveyPoint.h"
#include "AliAlignObjParams.h"
#include "AliEMCALGeometry.h"
#include "AliEMCALSurvey.h"
#include "AliLog.h"
ClassImp(AliEMCALSurvey)
AliEMCALSurvey::AliEMCALSurvey()
: fNSuperModule(0),
fSuperModuleData(0),
fDataType(kSurvey)
{
}
namespace {
struct AliEMCALSuperModuleCoords {
Double_t fX1;
Double_t fY1;
Double_t fZ1;
Double_t fPsi;
Double_t fTheta;
Double_t fPhi;
};
}
AliEMCALSurvey::AliEMCALSurvey(const TString &txtFileName,const SurveyDataType_t type)
: fNSuperModule(0),
fSuperModuleData(0),
fDataType(type)
{
const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
if (!geom) {
AliError("Cannot obtain AliEMCALGeometry instance.");
return;
}
fNSuperModule = geom->GetNumberOfSuperModules();
if(fDataType == kSurvey) {
AliSurveyObj *s1 = new AliSurveyObj();
s1->FillFromLocalFile(txtFileName);
TObjArray* points = s1->GetData();
InitSuperModuleData(points);
} else {
std::ifstream inputFile(txtFileName.Data());
if (!inputFile) {
AliError(("Cannot open the survey file " + txtFileName).Data());
return;
}
Int_t dummyInt = 0;
Double_t *xReal = new Double_t[fNSuperModule];
Double_t *yReal = new Double_t[fNSuperModule];
Double_t *zReal = new Double_t[fNSuperModule];
Double_t *psiReal = new Double_t[fNSuperModule];
Double_t *thetaReal = new Double_t[fNSuperModule];
Double_t *phiReal = new Double_t[fNSuperModule];
memset(xReal, 0,sizeof(Int_t)*fNSuperModule);
memset(yReal, 0,sizeof(Int_t)*fNSuperModule);
memset(zReal, 0,sizeof(Int_t)*fNSuperModule);
memset(psiReal, 0,sizeof(Int_t)*fNSuperModule);
memset(thetaReal, 0,sizeof(Int_t)*fNSuperModule);
memset(phiReal, 0,sizeof(Int_t)*fNSuperModule);
for (Int_t i = 0; i < fNSuperModule; ++i) {
if (!inputFile) {
AliError("Error while reading input file.");
delete [] xReal;
delete [] yReal;
delete [] zReal;
delete [] psiReal;
delete [] thetaReal;
delete [] phiReal;
return;
}
inputFile>>dummyInt>>xReal[i]>>yReal[i]>>zReal[i]>>psiReal[i]>>thetaReal[i]>>phiReal[i];
}
InitSuperModuleData(xReal, yReal, zReal, psiReal, thetaReal, phiReal);
delete [] xReal;
delete [] yReal;
delete [] zReal;
delete [] psiReal;
delete [] thetaReal;
delete [] phiReal;
}
}
AliEMCALSurvey::~AliEMCALSurvey()
{
delete [] fSuperModuleData;
}
void AliEMCALSurvey::CreateAliAlignObjParams(TClonesArray &array)
{
const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
if (!geom) {
AliError("Cannot obtain AliEMCALGeometry instance.");
return;
}
if (!gGeoManager) {
AliWarning("Cannot create local transformations for supermodules - gGeoManager does not exist.");
AliInfo("Null shifts and rotations will be created instead.");
return CreateNullObjects(array, geom);
}
Int_t arrayInd = array.GetEntries(), iIndex = 0;
AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iIndex);
AliAlignObjParams* myobj = 0x0;
TString SMName;
Int_t tmpType = -1;
Int_t SMOrder = 0;
for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
if(geom->GetSMType(smodnum) == AliEMCALGeometry::kEMCAL_Standard ) SMName = "FullSupermodule";
else if(geom->GetSMType(smodnum) == AliEMCALGeometry::kEMCAL_Half ) SMName = "HalfSupermodule";
else if(geom->GetSMType(smodnum) == AliEMCALGeometry::kEMCAL_3rd ) SMName = "OneThrdSupermodule";
else if( geom->GetSMType(smodnum) == AliEMCALGeometry::kDCAL_Standard ) SMName = "DCALSupermodule";
else if( geom->GetSMType(smodnum) == AliEMCALGeometry::kDCAL_Ext ) SMName = "DCALExtensionSM";
else AliError("Unkown SM Type!!");
if(geom->GetSMType(smodnum) == tmpType) {
SMOrder++;
} else {
tmpType = geom->GetSMType(smodnum);
SMOrder = 1;
}
TString smodName(TString::Format("EMCAL/%s%d", SMName.Data(), SMOrder));
AliEMCALSuperModuleDelta t(GetSuperModuleTransformation(smodnum));
new(array[arrayInd])
AliAlignObjParams(
smodName.Data(), volid,
t.fXShift, t.fYShift, t.fZShift,
-t.fPsi, -t.fTheta, -t.fPhi,
true
);
++arrayInd;
myobj = (AliAlignObjParams*)array.UncheckedAt(smodnum);
printf("==== AliAlignObjParams for SM %d ====\n",smodnum);
myobj->Print("");
}
}
void AliEMCALSurvey::CreateNullObjects(TClonesArray &array, const AliEMCALGeometry *geom)const
{
Int_t arrayInd = array.GetEntries(), iIndex = 0;
AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iIndex);
TString SMName;
Int_t tmpType = -1;
Int_t SMOrder = 0;
for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
if(geom->GetSMType(smodnum) == AliEMCALGeometry::kEMCAL_Standard ) SMName = "FullSupermodule";
else if(geom->GetSMType(smodnum) == AliEMCALGeometry::kEMCAL_Half ) SMName = "HalfSupermodule";
else if(geom->GetSMType(smodnum) == AliEMCALGeometry::kEMCAL_3rd ) SMName = "OneThrdSupermodule";
else if( geom->GetSMType(smodnum) == AliEMCALGeometry::kDCAL_Standard ) SMName = "DCALSupermodule";
else if( geom->GetSMType(smodnum) == AliEMCALGeometry::kDCAL_Ext ) SMName = "DCALExtensionSM";
else AliError("Unkown SM Type!!");
if(geom->GetSMType(smodnum) == tmpType) {
SMOrder++;
} else {
tmpType = geom->GetSMType(smodnum);
SMOrder = 1;
}
TString smodName(TString::Format("EMCAL/%s%d", SMName.Data(), SMOrder));
new(array[arrayInd]) AliAlignObjParams(smodName.Data(), volid, 0., 0., 0., 0., 0., 0., true);
++arrayInd;
}
}
AliEMCALSurvey::AliEMCALSuperModuleDelta AliEMCALSurvey::GetSuperModuleTransformation(Int_t supModIndex)const
{
AliEMCALSuperModuleDelta t = {0., 0., 0., 0., 0., 0.};
if (!fSuperModuleData)
return t;
return fSuperModuleData[supModIndex];
}
void AliEMCALSurvey::InitSuperModuleData(const TObjArray *svypts)
{
AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
Float_t pars[] = {geom->GetSuperModulesPar(0),geom->GetSuperModulesPar(1),geom->GetSuperModulesPar(2)};
Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
Float_t fInnerEdge = geom->GetDCALInnerEdge();
Double_t phi=0, phiRad=0, xpos=0, ypos=0, zpos=0;
AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
phiRad = geom->GetPhiCenterOfSMSec(smodnum);
phi = phiRad*180./TMath::Pi();
xpos = rpos * TMath::Cos(phiRad);
ypos = rpos * TMath::Sin(phiRad);
zpos = pars[2];
if( geom->GetSMType(smodnum) == AliEMCALGeometry::kEMCAL_Half
|| geom->GetSMType(smodnum) == AliEMCALGeometry::kEMCAL_3rd
|| geom->GetSMType(smodnum) == AliEMCALGeometry::kDCAL_Ext ) {
xpos += (pars[1]/2. * TMath::Sin(phiRad));
ypos -= (pars[1]/2. * TMath::Cos(phiRad));
} else if( geom->GetSMType(smodnum) == AliEMCALGeometry::kDCAL_Standard) {
zpos = pars[2] + fInnerEdge/2.;
}
smc.fX1 = xpos;
smc.fY1 = ypos;
smc.fPhi = phi;
smc.fTheta = 0.;
smc.fPsi = 0.;
if(smodnum%2==0) {
smc.fZ1 = zpos;
} else {
smc.fZ1 = -zpos;
}
printf("<SM %d> IDEAL x,y,z positions: %.2f,%.2f,%.2f, IDEAL phi,theta,psi angles: %.2f,%.2f,%.2f\n",smodnum,smc.fX1,smc.fY1,smc.fZ1,smc.fPhi,smc.fTheta,smc.fPsi);
}
const Int_t buffersize = 255;
char substr[buffersize];
AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
AliEMCALSuperModuleCoords &smc = realSM[smodnum];
Double_t zLength = pars[2]*2.;
Double_t halfHeight = pars[0];
Double_t halfWidth = pars[1];
printf("AliEMCALGeometry says zlength = %.2f, halfheight = %.2f, halfwidth = %.2f\n",zLength,halfHeight,halfWidth);
snprintf(substr,buffersize,"4096%d",smodnum);
std::vector<Double_t> xval;
std::vector<Double_t> yval;
std::vector<Double_t> zval;
for(Int_t i = 0; i < svypts->GetEntries(); i++) {
AliSurveyPoint* pt = (AliSurveyPoint*)svypts->At(i);
TString ptname = pt->GetPointName();
if(ptname.Contains(substr)) {
xval.push_back(pt->GetX()*100.);
yval.push_back(pt->GetY()*100.);
zval.push_back(pt->GetZ()*100.);
}
}
Double_t activeX = ((xval[0] + (xval[2] - xval[0])/2.)
+(xval[1] + (xval[3] - xval[1])/2.) ) /2.;
Double_t activeY = ((yval[0] + (yval[2] - yval[0])/2.)
+(yval[1] + (yval[3] - yval[1])/2.) ) /2.;
Double_t activeZ = ((zval[0] + (zval[2] - zval[0])/2.)
+(zval[1] + (zval[3] - zval[1])/2.) ) /2.;
printf("Bottom Center of active area of SM %s: %.2f, %.2f, %.2f\n",substr,activeX,activeY,activeZ);
Double_t realphi = 0.;
if(smodnum%2 == 0) {
realphi = (TMath::ATan((yval[2] - yval[0])/(xval[2] - xval[0]))
+TMath::ATan((yval[3] - yval[1])/(xval[3] - xval[1])) )/2.;
} else {
realphi = (TMath::ATan((yval[0] - yval[2])/(xval[0] - xval[2]))
+TMath::ATan((yval[1] - yval[3])/(xval[1] - xval[3])) )/2.;
}
Double_t realpsi = (TMath::ATan((zval[0] - zval[2])/(xval[2] - xval[0]))
+TMath::ATan((zval[1] - zval[3])/(xval[3] - xval[1])) )/2.;
Double_t realtheta = TMath::Pi()/2.
+ (TMath::ATan((zval[2] - zval[3])/(yval[3] - yval[2]))
+TMath::ATan((zval[0] - zval[1])/(yval[1] - yval[0])) )/2.;
printf("Old edge of %s 01: %.2f, %.2f, %.2f\n",substr,xval[1],yval[1],zval[1]);
printf("Old edge of %s 11: %.2f, %.2f, %.2f\n",substr,xval[3],yval[3],zval[3]);
printf("Real theta angle (degrees) = %.2f\n",realtheta*TMath::RadToDeg());
Double_t activeLength = TMath::Abs(((zval[1] - zval[0]) + (zval[3] - zval[2]))/2.);
printf("ACTIVE LENGTH = %.2f\n",activeLength);
if(smodnum%2 == 0) {
yval[1] += (zLength - activeLength)*sin(realtheta);
yval[3] += (zLength - activeLength)*sin(realtheta);
zval[1] += (zLength - activeLength)*cos(realtheta);
zval[3] += (zLength - activeLength)*cos(realtheta);
} else {
yval[1] -= (zLength - activeLength)*sin(realtheta);
yval[3] -= (zLength - activeLength)*sin(realtheta);
zval[1] -= (zLength - activeLength)*cos(realtheta);
zval[3] -= (zLength - activeLength)*cos(realtheta);
}
printf("New extended edge of %s 01: %.2f, %.2f, %.2f\n",substr,xval[1],yval[1],zval[1]);
printf("New extended edge of %s 11: %.2f, %.2f, %.2f\n",substr,xval[3],yval[3],zval[3]);
Double_t realX = activeX;
Double_t realY = ((yval[0] + (yval[2] - yval[0])/2.)
+(yval[1] + (yval[3] - yval[1])/2.) ) /2.;
Double_t realZ = ((zval[0] + (zval[2] - zval[0])/2.)
+(zval[1] + (zval[3] - zval[1])/2.) ) /2.;
printf("Bottom Center of SM %s Box: %.2f, %.2f, %.2f\n",substr,realX,realY,realZ);
realX += halfHeight*TMath::Cos(TMath::Pi()/2+realphi);
realY += halfHeight*(TMath::Sin(TMath::Pi()/2+realphi) + TMath::Sin(realtheta));
realZ += halfHeight*TMath::Cos(TMath::Pi()/2-realtheta);
printf("Rotation angles of SM %s (phi,psi,theta) in degrees: %.4f, %.4f, %.4f\n",substr,realphi*TMath::RadToDeg(),realpsi*TMath::RadToDeg(),realtheta*TMath::RadToDeg());
printf("Middle of SM %s: %.2f, %.2f, %.2f\n\n",substr,realX,realY,realZ);
smc.fX1 = realX;
smc.fY1 = realY;
smc.fZ1 = realZ;
smc.fPhi = 90. + realphi*TMath::RadToDeg();
smc.fTheta = 0. + realtheta*TMath::RadToDeg();
smc.fPsi = 0. + realpsi*TMath::RadToDeg();
}
for (Int_t i = 0; i < fNSuperModule; i++) {
if(i%2==0) {
AliEMCALSuperModuleCoords &realA = realSM[i];
AliEMCALSuperModuleCoords &realC = realSM[i+1];
Double_t avgx = (realA.fX1 + realC.fX1)/2.;
Double_t avgy = (realA.fY1 + realC.fY1)/2.;
Double_t avgphi = (realA.fPhi + realC.fPhi)/2.;
Double_t avgtheta = (realA.fTheta + realC.fTheta)/2.;
Double_t avgpsi = (realA.fPsi + realC.fPsi)/2.;
printf("AVERAGE VALUES: %.2f,%.2f,%.2f,%.2f,%.2f\n",avgx,avgy,avgphi,avgtheta,avgpsi);
realA.fX1 = avgx; realC.fX1 = avgx;
realA.fY1 = avgy; realC.fY1 = avgy;
realA.fPhi = avgphi; realC.fPhi = avgphi;
realA.fTheta = avgtheta; realC.fTheta = avgtheta;
realA.fPsi = avgpsi; realC.fPhi = avgphi;
}
}
fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
for (Int_t i = 0; i < fNSuperModule; ++i) {
const AliEMCALSuperModuleCoords &real = realSM[i];
const AliEMCALSuperModuleCoords &ideal = idealSM[i];
AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
t.fXShift = real.fX1 - ideal.fX1;
t.fYShift = real.fY1 - ideal.fY1;
t.fZShift = real.fZ1 - ideal.fZ1;
t.fPhi = real.fPhi - ideal.fPhi;
t.fTheta = real.fTheta - ideal.fTheta;
t.fPsi = real.fPsi - ideal.fPsi;
printf("===================== SM %d =======================\n",i);
printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
printf("real theta (%.2f) - ideal theta (%.2f) = shift in theta %.2f\n",real.fTheta,ideal.fTheta,t.fTheta);
printf("real psi (%.2f) - ideal psi (%.2f) = shift in psi %.2f\n",real.fPsi,ideal.fPsi,t.fPsi);
printf("real phi (%.2f) - ideal phi (%.2f) = shift in phi %.2f\n",real.fPhi,ideal.fPhi,t.fPhi);
printf("===================================================\n");
}
delete [] realSM;
delete [] idealSM;
}
void AliEMCALSurvey::InitSuperModuleData(const Double_t *xReal, const Double_t *yReal,
const Double_t *zReal, const Double_t *psiReal,
const Double_t *thetaReal, const Double_t *phiReal)
{
AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
Float_t pars[] = {geom->GetSuperModulesPar(0),geom->GetSuperModulesPar(1),geom->GetSuperModulesPar(2)};
Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
Float_t fInnerEdge = geom->GetDCALInnerEdge();
Double_t phi=0, phiRad=0, xpos=0, ypos=0, zpos=0;
zpos = pars[2];
AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
phiRad = geom->GetPhiCenterOfSMSec(smodnum);
phi = phiRad*180./TMath::Pi();
xpos = rpos * TMath::Cos(phiRad);
ypos = rpos * TMath::Sin(phiRad);
if( geom->GetSMType(smodnum) == AliEMCALGeometry::kEMCAL_Half
|| geom->GetSMType(smodnum) == AliEMCALGeometry::kEMCAL_3rd
|| geom->GetSMType(smodnum) == AliEMCALGeometry::kDCAL_Ext ) {
xpos += (pars[1]/2. * TMath::Sin(phiRad));
ypos -= (pars[1]/2. * TMath::Cos(phiRad));
} else if( geom->GetSMType(smodnum) == AliEMCALGeometry::kDCAL_Standard) {
zpos = pars[2] + fInnerEdge/2.;
}
smc.fX1 = xpos;
smc.fY1 = ypos;
smc.fPhi = phi;
smc.fTheta = 0.;
smc.fPsi = 0.;
if(smodnum%2==0) {
smc.fZ1 = zpos;
} else {
smc.fZ1 = -zpos;
}
}
AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
AliEMCALSuperModuleCoords &smc = realSM[smodnum];
smc.fX1 = xReal[smodnum];
smc.fY1 = yReal[smodnum];
smc.fZ1 = zReal[smodnum];
smc.fTheta = thetaReal[smodnum];
smc.fPsi = psiReal[smodnum];
smc.fPhi = phiReal[smodnum];
}
fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
for (Int_t i = 0; i < fNSuperModule; ++i) {
const AliEMCALSuperModuleCoords &real = realSM[i];
const AliEMCALSuperModuleCoords &ideal = idealSM[i];
AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
t.fTheta = real.fTheta - ideal.fTheta;
t.fPsi = 0.;
t.fPhi = real.fPhi - ideal.fPhi;
t.fXShift = real.fX1 - ideal.fX1;
t.fYShift = real.fY1 - ideal.fY1;
t.fZShift = real.fZ1 - ideal.fZ1;
printf("===================== SM %d =======================\n",i);
printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
printf("real theta (%.2f) - ideal theta (%.2f) = shift in theta %.2f\n",real.fTheta,ideal.fTheta,t.fTheta);
printf("real psi (%.2f) - ideal psi (%.2f) = shift in psi %.2f\n",real.fPsi,ideal.fPsi,t.fPsi);
printf("real phi (%.2f) - ideal phi (%.2f) = shift in phi %.2f\n",real.fPhi,ideal.fPhi,t.fPhi);
printf("===================================================\n");
}
delete [] realSM;
delete [] idealSM;
}