--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+
+// ---
+// Macro to process survey and photogrammetry data of chamber 8L
+//
+// Macro loads the survey data from .txt file using AliSurveyObj.
+// Macro MUONSurveyUtil.C is then loaded.
+// The transformations of the slats are obatained in 2 steps:
+// 1. Fit a plane to the sticker targets -> psi, theta
+// 2. Using above psi in theta obtain xc, yc, zc and phi by solving
+// the equations from a local to global transformation of the
+// fixed button targets
+// Various histograms are filled and printed for monitoring.
+// MisAlignment object is then created.
+//
+// Author: Javier Castillo
+// ---
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "AliMUONGeometryTransformer.h"
+#include "AliMUONGeometryMisAligner.h"
+
+#include "AliSurveyObj.h"
+#include "AliSurveyPoint.h"
+#include "AliGeomManager.h"
+#include "AliCDBManager.h"
+#include "AliCDBMetaData.h"
+#include "AliCDBId.h"
+
+#include <TROOT.h>
+#include <TGeoManager.h>
+#include <TClonesArray.h>
+#include <TObjArray.h>
+#include <TObjString.h>
+#include <TMath.h>
+#include <TString.h>
+#include <Riostream.h>
+#include <TF2.h>
+#include <TH2.h>
+#include <TGraph2DErrors.h>
+#include <TGraph.h>
+#include <TVector3.h>
+#include <TCanvas.h>
+#include <TPad.h>
+#include <TPostScript.h>
+#include <TPaveLabel.h>
+#include <TStyle.h>
+
+#include <fstream>
+
+#endif
+
+Bool_t MatrixToAngles(const Double_t *rot, Double_t *angles);
+Double_t eqPlane(Double_t *x, Double_t *par);
+Double_t xpCenter(Double_t *x, Double_t *par);
+Double_t xnCenter(Double_t *x, Double_t *par);
+Double_t ypCenter(Double_t *x, Double_t *par);
+Double_t ynCenter(Double_t *x, Double_t *par);
+Double_t zpCenter(Double_t *x, Double_t *par);
+Double_t znCenter(Double_t *x, Double_t *par);
+Double_t phixpp(Double_t *x, Double_t *par);
+Double_t phixpn(Double_t *x, Double_t *par);
+Double_t phixnp(Double_t *x, Double_t *par);
+Double_t phixnn(Double_t *x, Double_t *par);
+Double_t phiypp(Double_t *x, Double_t *par);
+Double_t phiypn(Double_t *x, Double_t *par);
+Double_t phiynp(Double_t *x, Double_t *par);
+Double_t phiynn(Double_t *x, Double_t *par);
+AliMUONGeometryTransformer *ReAlign(const AliMUONGeometryTransformer * transformer,
+ int rMod, TGeoCombiTrans deltaDetElemTransf[], Bool_t verbose);
+void MUONSurveyCh8L() {
+
+ char str[100];
+ char opt[100];
+ char var[100];
+ char filename[100];
+ Char_t histoName[20];
+ Char_t histoTitle[50];
+
+ int saveps = 1;
+ const int font = 41; // Helvetica + precision 1
+ const int cWidth = (int)(700*(29./21.));
+ const int cHeight = 700;
+ const int lineColor = 1;
+ const int fillColor = 4;
+// const int filetype = 111; // portrait
+ const int filetype = 112; // landscape
+
+ sprintf(filename,"surveyChamber8L.ps");
+
+ Int_t nSlats = 13;
+
+ gROOT->LoadMacro("MUONSurveyUtil.C+");
+
+ AliSurveyObj *so = new AliSurveyObj();
+
+ Int_t size = so->GetEntries();
+ printf("-> %d\n", size);
+
+ so->FillFromLocalFile("Alice_MuonSystem_Chamber8LCavern_3561b.txt");
+ size = so->GetEntries();
+ printf("--> %d\n", size);
+
+ Printf("Title: \"%s\"", so->GetReportTitle().Data());
+ Printf("Date: \"%s\"", so->GetReportDate().Data());
+ Printf("Detector: \"%s\"", so->GetDetector().Data());
+ Printf("URL: \"%s\"", so->GetURL().Data());
+ Printf("Number: \"%d\"", so->GetReportNumber());
+ Printf("Version: \"%d\"", so->GetReportVersion());
+ Printf("Observations: \"%s\"", so->GetObservations().Data());
+ Printf("Coordinate System: \"%s\"", so->GetCoordSys().Data());
+ Printf("Measurement Units: \"%s\"", so->GetUnits().Data());
+ Printf("Nr Columns: \"%d\"", so->GetNrColumns());
+
+ TObjArray *colNames = so->GetColumnNames();
+ for (Int_t i = 0; i < colNames->GetEntries(); ++i)
+ Printf(" Column %d --> \"%s\"", i, ((TObjString *) colNames->At(i))->GetString().Data());
+
+ // Get Array of surveyed points
+ Printf("Points:");
+ TObjArray *points = so->GetData();
+
+ for (Int_t i = 0; i < points->GetEntries(); ++i)
+ Printf(" Point %d --> \"%s\" %s ", i, ((AliSurveyPoint *) points->At(i))->GetPointName().Data(), points->At(i)->GetName());
+
+
+ // Slats (#1 - #13) button targets local coordinates
+ Double_t lSBTLoc6[13][2][3] = {{{ -412.50, 0.0, -(11.75+ 8.20+20.00)},
+ { 412.50, 0.0, -(11.75+16.50+20.00)}},
+ {{ -612.50, 0.0, (11.75+ 8.20+20.00)},
+ { 612.50, 0.0, (11.75+16.50+20.00)}},
+ {{ - 812.50, 0.0, -(11.75+ 8.20+20.00)},
+ { 812.50, 0.0, -(11.75+16.50+20.00)}},
+ {{ -1012.50, 0.0, (11.75+ 8.20+20.00)},
+ { 1012.50, 0.0, (11.75+16.50+20.00)}},
+ {{ -1012.50, 0.0, -(11.75+ 8.20+20.00)},
+ { 1012.50, 0.0, -(11.75+16.50+20.00)}},
+ {{ -1212.50, 5.0, (11.75+ 8.20+20.00)},
+ { 1212.50, 0.0, (11.75+16.50+20.00)}},
+ {{ -1012.50, 0.0, -(11.75+ 8.20+20.00)},
+ { 1012.50, 0.0, -(11.75+16.50+20.00)}},
+ {{ -1212.50, 5.0, (11.75+ 8.20+20.00)},
+ { 1212.50, 0.0, (11.75+16.50+20.00)}},
+ {{ -1012.50, 0.0, -(11.75+ 8.20+20.00)},
+ { 1012.50, 0.0, -(11.75+16.50+20.00)}},
+ {{ -1012.50, 0.0, (11.75+ 8.20+20.00)},
+ { 1012.50, 0.0, (11.75+16.50+20.00)}},
+ {{ -812.50, 0.0, -(11.75+ 8.20+20.00)},
+ { 812.50, 0.0, -(11.75+16.50+20.00)}},
+ {{ -612.50, 0.0, (11.75+ 8.20+20.00)},
+ { 612.50, 0.0, (11.75+16.50+20.00)}},
+ {{ -412.50, 0.0, -(11.75+ 8.20+20.00)},
+ { 412.50, 0.0, -(11.75+16.50+20.00)}}};
+
+
+ AliSurveyPoint *pointSST = 0;
+ AliSurveyPoint *pointCST = 0;
+ AliSurveyPoint *pointCPST = 0;
+ AliSurveyPoint **pointSBT = new AliSurveyPoint*[2];
+
+ char sPointName[10] = "5000";
+
+ // Print length of slats
+ cout << "Slat lengths:" << endl;
+ TVector3 vTemp(0., 0., 0.);
+ TVector3 vSBT(0., 0., 0.);
+
+ for (Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ // Get button targets survey points
+ vTemp.SetXYZ(0., 0., 0.);
+ if (iSlat+1<10) {
+ sprintf(sPointName,"60%d%d",iSlat+1,1);
+ pointSBT[0] = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointSBT[0]) {
+ cout << "Error! No button targets ... " << endl;
+ break;
+ }
+ vSBT.SetXYZ(pointSBT[0]->GetX(),pointSBT[0]->GetY(),pointSBT[0]->GetZ());
+ vTemp+=vSBT;
+ sprintf(sPointName,"60%d%d",iSlat+1,2);
+ pointSBT[1] = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointSBT[1]) {
+ cout << "Error! No button targets ... " << endl;
+ break;
+ }
+ vSBT.SetXYZ(pointSBT[1]->GetX(),pointSBT[1]->GetY(),pointSBT[1]->GetZ());
+ vTemp-=vSBT;
+ }
+ else {
+ sprintf(sPointName,"6%d%d",iSlat+1,1);
+ pointSBT[0] = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointSBT[0]) {
+ cout << "Error! No button targets ... " << endl;
+ break;
+ }
+ vSBT.SetXYZ(pointSBT[0]->GetX(),pointSBT[0]->GetY(),pointSBT[0]->GetZ());
+ vTemp+=vSBT;
+ sprintf(sPointName,"6%d%d",iSlat+1,2);
+ pointSBT[1] = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointSBT[1]) {
+ cout << "Error! No button targets ... " << endl;
+ break;
+ }
+ vSBT.SetXYZ(pointSBT[1]->GetX(),pointSBT[1]->GetY(),pointSBT[1]->GetZ());
+ vTemp-=vSBT;
+ }
+ cout << "Slat " << iSlat+1 << ": " << vTemp.Mag() << endl;
+ }
+
+ // Histograms for monitoring
+ TH2F *hCPSTry = new TH2F("hCPSTry","hCPSTry",28,-600,2200,52,0,5200);
+ TH2F *hCPSTly = new TH2F("hCPSTly","hCPSTly",28,-600,2200,52,0,5200);
+ TH2F *hSSTry = new TH2F("hSSTry","hSSTry",70,-600,2200,130,0,5200);
+ TH2F *hSSTly = new TH2F("hSSTly","hSSTly",70,-600,2200,130,0,5200);
+ TH2F *hCSTy = new TH2F("hCSTy","hCSTy",70,-600,2200,130,0,5200);
+ TH2F *hSSTrpy = new TH2F("hSSTrpy","hSSTrpy",70,-600,2200,130,0,5200);
+ TH2F *hSSTlpy = new TH2F("hSSTlpy","hSSTlpy",70,-600,2200,130,0,5200);
+
+ // Chamber plane sticker targets
+ for (int iPoint=0; iPoint<9; iPoint++) {
+ sprintf(sPointName,"700%d",iPoint+1);
+ pointCPST = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointCPST) {
+ printf("Point %s is missing ...\n",sPointName);
+ break;
+ }
+ hCPSTry->Fill(pointCPST->GetX(),pointCPST->GetZ(),pointCPST->GetY());
+ }
+ for (int iPoint=9; iPoint<18; iPoint++) {
+ sprintf(sPointName,"70%d",iPoint+1);
+ pointCPST = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointCPST) {
+ printf("Point %s is missing ...\n",sPointName);
+ break;
+ }
+ hCPSTly->Fill(pointCPST->GetX(),pointCPST->GetZ(),pointCPST->GetY());
+ }
+
+ // Chamber Side Targets
+ for (int iPoint=0; iPoint<25; iPoint++) {
+ if (iPoint+1<10) {
+ sprintf(sPointName,"800%d",iPoint+1);
+ }
+ else {
+ sprintf(sPointName,"80%d",iPoint+1);
+ }
+ pointCST = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointCPST) {
+ printf("Point %s is missing ...\n",sPointName);
+ break;
+ }
+ hCSTy->Fill(pointCST->GetX(),pointCST->GetZ(),pointCST->GetY());
+ }
+
+
+ // Graphs used for plane fitting
+ TGraph2DErrors **gSST5 = new TGraph2DErrors*[nSlats];
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ gSST5[iSlat] = new TGraph2DErrors();
+ }
+
+ // Keep the id of slat sticker targets next to slat button targets
+ Int_t iSSBT[13][2] = {0};
+
+ // Fill graph with sticker target positions
+ for (int iSlat=0; iSlat<nSlats; iSlat++){
+ for (int iPoint=0; iPoint<9; iPoint++) {
+ // Second sticker target is next to first button target
+ // Previous to last sticker target is next to second button target
+ iSSBT[iSlat][0] = 2;
+ iSSBT[iSlat][1] = 8;
+ if (iSlat+1<10) {
+ sprintf(sPointName,"50%d%d",iSlat+1,iPoint+1);
+ }
+ else {
+ sprintf(sPointName,"5%d%d",iSlat+1,iPoint+1);
+ }
+ pointSST = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointSST) {
+ printf("%s\n",sPointName);
+ cout << iSlat << " " << iPoint << " " << pointSST << endl;
+ // Previous to last sticker target is next to second button target
+ iSSBT[iSlat][1] = iPoint+1-2;
+ break;
+ }
+
+ gSST5[iSlat]->SetPoint(iPoint,-pointSST->GetX(),pointSST->GetZ(),pointSST->GetY());
+ gSST5[iSlat]->SetPointError(iPoint,pointSST->GetPrecisionX(),pointSST->GetPrecisionZ(),pointSST->GetPrecisionY());
+
+ // Fill histograms of sticker targets. For monitoring purposes.
+ if((iSlat+1)%2==0){
+ hSSTly->Fill(pointSST->GetX(),pointSST->GetZ(),pointSST->GetY());
+ }
+ else {
+ hSSTry->Fill(pointSST->GetX(),pointSST->GetZ(),pointSST->GetY());
+ }
+ }
+ }
+
+ Float_t xMin = -2200.;
+ Float_t xMax = 2200.;
+ Float_t yMin = -5200.;
+ Float_t yMax = 5200.;
+ Float_t zMin = -200.;
+ Float_t zMax = 200.;
+
+ Double_t xMinSlat, xMaxSlat;
+ Double_t yMinSlat, yMaxSlat;
+
+ // Slat plane function
+ char fsName[100] = "fSlat00";
+ TF2 **fSlat = new TF2*[nSlats];
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ sprintf(fsName,"fSlat%d",iSlat+1);
+ fSlat[iSlat] = new TF2(fsName,eqPlane,xMin,xMax,yMin,yMax,3);
+ }
+
+ // Xcenter functions
+ char fxcName[100] = "fXcnSlat00";
+ TF2 ***fXcSlat = new TF2**[nSlats];
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ fXcSlat[iSlat] = new TF2*[2];
+ sprintf(fxcName,"fXcnSlat%d",iSlat+1);
+ fXcSlat[iSlat][0] = new TF2(fxcName,xnCenter,xMin,xMax,yMin,yMax,7);
+ sprintf(fxcName,"fXcpSlat%d",iSlat+1);
+ fXcSlat[iSlat][1] = new TF2(fxcName,xpCenter,xMin,xMax,yMin,yMax,7);
+ }
+
+ // Ycenter functions
+ char fycName[100] = "fYcnSlat00";
+ TF2 ***fYcSlat = new TF2**[nSlats];
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ fYcSlat[iSlat] = new TF2*[2];
+ sprintf(fycName,"fYcnSlat%d",iSlat+1);
+ fYcSlat[iSlat][0] = new TF2(fycName,ynCenter,yMin,yMax,yMin,yMax,8);
+ sprintf(fycName,"fYcpSlat%d",iSlat+1);
+ fYcSlat[iSlat][1] = new TF2(fycName,ypCenter,yMin,yMax,yMin,yMax,8);
+ }
+
+ // Zcenter functions
+ char fzcName[100] = "fZcnSlat00";
+ TF2 ***fZcSlat = new TF2**[nSlats];
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ fZcSlat[iSlat] = new TF2*[2];
+ sprintf(fzcName,"fZcnSlat%d",iSlat+1);
+ fZcSlat[iSlat][0] = new TF2(fzcName,znCenter,zMin,zMax,zMin,zMax,8);
+ sprintf(fzcName,"fZcpSlat%d",iSlat+1);
+ fZcSlat[iSlat][1] = new TF2(fzcName,zpCenter,zMin,zMax,zMin,zMax,8);
+ }
+
+ // Phi rotation using xglobal coords functions
+ char fphixName[100] = "fPhiXnnSlat00";
+ TF2 ****fPhiXSlat = new TF2***[nSlats];
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ fPhiXSlat[iSlat] = new TF2**[2];
+ for (Int_t iX =0; iX<2; iX++)
+ fPhiXSlat[iSlat][iX] = new TF2*[2];
+ sprintf(fphixName,"fPhiXnnSlat%d",iSlat+1);
+ fPhiXSlat[iSlat][0][0] = new TF2(fphixName,phixnn,xMin,xMax,xMin,xMax,7);
+ sprintf(fphixName,"fPhixnpSlat%d",iSlat+1);
+ fPhiXSlat[iSlat][0][1] = new TF2(fphixName,phixnp,xMin,xMax,xMin,xMax,7);
+ sprintf(fphixName,"fPhiXpnSlat%d",iSlat+1);
+ fPhiXSlat[iSlat][1][0] = new TF2(fphixName,phixpn,xMin,xMax,xMin,xMax,7);
+ sprintf(fphixName,"fPhixppSlat%d",iSlat+1);
+ fPhiXSlat[iSlat][1][1] = new TF2(fphixName,phixpp,xMin,xMax,xMin,xMax,7);
+ }
+
+ // Phi rotation using yglobal coords functions
+ char fphiyName[100] = "fPhiYnnSlat00";
+ TF2 ****fPhiYSlat = new TF2***[nSlats];
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ fPhiYSlat[iSlat] = new TF2**[2];
+ for (Int_t iY =0; iY<2; iY++)
+ fPhiYSlat[iSlat][iY] = new TF2*[2];
+ sprintf(fphiyName,"fPhiYnnSlat%d",iSlat+1);
+ fPhiYSlat[iSlat][0][0] = new TF2(fphiyName,phiynn,yMin,yMax,yMin,yMax,7);
+ sprintf(fphiyName,"fPhiYnpSlat%d",iSlat+1);
+ fPhiYSlat[iSlat][0][1] = new TF2(fphiyName,phiynp,yMin,yMax,yMin,yMax,7);
+ sprintf(fphiyName,"fPhiYpnSlat%d",iSlat+1);
+ fPhiYSlat[iSlat][1][0] = new TF2(fphiyName,phiypn,yMin,yMax,yMin,yMax,7);
+ sprintf(fphiyName,"fPhiYppSlat%d",iSlat+1);
+ fPhiYSlat[iSlat][1][1] = new TF2(fphiyName,phiypp,yMin,yMax,yMin,yMax,7);
+ }
+
+ Double_t *xce = new Double_t[nSlats];
+ Double_t *yce = new Double_t[nSlats];
+ Double_t *zce = new Double_t[nSlats];
+ Double_t *psi = new Double_t[nSlats];
+ Double_t *tht = new Double_t[nSlats];
+ Double_t *phi = new Double_t[nSlats];
+
+ Double_t **lCenSlat = new Double_t*[nSlats];
+ Double_t **lRotSlat = new Double_t*[nSlats];
+ Double_t **lDiffCenSlat0 = new Double_t*[nSlats];
+ Double_t **lDiffRotSlat0 = new Double_t*[nSlats];
+ Double_t **lDiffThCenSlat0 = new Double_t*[nSlats];
+ Double_t **lDiffThRotSlat0 = new Double_t*[nSlats];
+ Double_t **lDeltaDiffCenSlat0 = new Double_t*[nSlats];
+ Double_t **lDeltaDiffRotSlat0 = new Double_t*[nSlats];
+
+ for (int iSlat=0; iSlat<nSlats; iSlat++) {
+ lCenSlat[iSlat] = new Double_t[3];
+ lRotSlat[iSlat] = new Double_t[3];
+
+ lDiffCenSlat0[iSlat] = new Double_t[3];
+ lDiffRotSlat0[iSlat] = new Double_t[3];
+ lDiffThCenSlat0[iSlat] = new Double_t[3];
+ lDiffThRotSlat0[iSlat] = new Double_t[3];
+ lDeltaDiffCenSlat0[iSlat] = new Double_t[3];
+ lDeltaDiffRotSlat0[iSlat] = new Double_t[3];
+ }
+
+
+ TGeoTranslation transSlat[nSlats];
+ TGeoRotation rotSlat[nSlats];
+ TGeoCombiTrans trfSlat[nSlats];
+
+ TGeoTranslation dtransSlat[nSlats];
+ TGeoRotation drotSlat[nSlats];
+ TGeoCombiTrans dtrfSlat[nSlats];
+
+ TGeoTranslation transTemp;
+ TGeoRotation rotTemp;
+ TGeoCombiTrans trfTemp;
+
+
+ Double_t lCenTemp[3];
+ Double_t lRotTemp[3];
+
+ Double_t lDiffTemp[9];
+ Double_t lDiffMin[9];
+
+ double tempDiff = 0.;
+ double tempDiff1 = 0.;
+ double tempDiff2 = 0.;
+
+ AliSurveyPoint **pointSSBT = new AliSurveyPoint*[2];
+
+ //
+ // Get Slat transformation.
+ // Psi and Theta are obtained by fitting a plane to the sticker targets.
+ // Then Xc, Yc, Zc and Phi are obtained by solving the equations to the ref.
+ // syst. transformation of the button targets
+ //
+ for (int iSlat=0; iSlat<nSlats; iSlat++) {
+ sprintf(fsName,"fSlat%d",iSlat+1);
+ cout << "Fitting Slat" << iSlat+1 << " ..." << endl;
+
+ // Fit a plane to the sticker targets
+ gSST5[iSlat]->Fit(fsName,"","same");
+
+ psi[iSlat] = TMath::ATan(fSlat[iSlat]->GetParameter(1));
+ tht[iSlat] = TMath::ATan(fSlat[iSlat]->GetParameter(0));
+ if (iSlat==5)
+ psi[iSlat] += TMath::Pi(); // Rotated slat
+
+ lRotSlat[iSlat][0] = psi[iSlat];
+ lRotSlat[iSlat][1] = tht[iSlat];
+
+ for(Int_t iS=0; iS<2; iS++){
+ fXcSlat[iSlat][iS]->SetParameters(lSBTLoc6[iSlat][0][0],lSBTLoc6[iSlat][0][1],lSBTLoc6[iSlat][0][2],lSBTLoc6[iSlat][1][0],lSBTLoc6[iSlat][1][1],lSBTLoc6[iSlat][1][2],tht[iSlat]);
+ fYcSlat[iSlat][iS]->SetParameters(lSBTLoc6[iSlat][0][0],lSBTLoc6[iSlat][0][1],lSBTLoc6[iSlat][0][2],lSBTLoc6[iSlat][1][0],lSBTLoc6[iSlat][1][1],lSBTLoc6[iSlat][1][2],psi[iSlat],tht[iSlat]);
+ fZcSlat[iSlat][iS]->SetParameters(lSBTLoc6[iSlat][0][0],lSBTLoc6[iSlat][0][1],lSBTLoc6[iSlat][0][2],lSBTLoc6[iSlat][1][0],lSBTLoc6[iSlat][1][1],lSBTLoc6[iSlat][1][2],psi[iSlat],tht[iSlat]);
+ for(Int_t jS=0; jS<2; jS++){
+ fPhiXSlat[iSlat][iS][jS]->SetParameters(lSBTLoc6[iSlat][0][0],lSBTLoc6[iSlat][0][1],lSBTLoc6[iSlat][0][2],lSBTLoc6[iSlat][1][0],lSBTLoc6[iSlat][1][1],lSBTLoc6[iSlat][1][2],tht[iSlat]);
+ fPhiYSlat[iSlat][iS][jS]->SetParameters(lSBTLoc6[iSlat][0][0],lSBTLoc6[iSlat][0][1],lSBTLoc6[iSlat][0][2],lSBTLoc6[iSlat][1][0],lSBTLoc6[iSlat][1][1],lSBTLoc6[iSlat][1][2],psi[iSlat],tht[iSlat]);
+ }
+ }
+
+ //
+ // Calculate Slat Center from button targets
+ //
+
+ // Get button targets survey points
+ for (Int_t iPoint=0; iPoint<2; iPoint++) {
+ if (iSlat+1<10) {
+ sprintf(sPointName,"60%d%d",iSlat+1,iPoint+1);
+ pointSBT[iPoint] = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointSBT[iPoint]) {
+ cout << "Error! No button targets ... " << endl;
+ break;
+ }
+ sprintf(sPointName,"50%d%d",iSlat+1,iSSBT[iSlat][iPoint]);
+ pointSSBT[iPoint] = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointSSBT[iPoint]) {
+ cout << "Error! No sticker target ... " << sPointName << endl;
+ break;
+ }
+ }
+ else {
+ sprintf(sPointName,"6%d%d",iSlat+1,iPoint+1);
+ pointSBT[iPoint] = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointSBT[iPoint]) {
+ cout << "Error! No button targets ... " << endl;
+ break;
+ }
+ sprintf(sPointName,"5%d%d",iSlat+1,iSSBT[iSlat][iPoint]);
+ pointSSBT[iPoint] = (AliSurveyPoint *)points->FindObject(sPointName);
+ if(!pointSSBT[iPoint]) {
+ cout << "Error! No sticker targets ... " << sPointName << endl;
+ break;
+ }
+ }
+ }
+
+ tempDiff += TMath::Power(-1,iSlat)*((pointSBT[1]->GetY() - pointSSBT[1]->GetY())-(pointSBT[0]->GetY() - pointSSBT[0]->GetY()));
+ tempDiff1 += TMath::Abs(pointSBT[0]->GetY() - pointSSBT[0]->GetY())-20;
+ tempDiff2 += TMath::Abs(pointSBT[1]->GetY() - pointSSBT[1]->GetY())-20;
+ cout << "BSdiff: " << TMath::Abs(pointSBT[0]->GetY() - pointSSBT[0]->GetY()) << " " << TMath::Abs(pointSBT[1]->GetY() - pointSSBT[1]->GetY()) << " " << tempDiff1/(iSlat+1) << " " << tempDiff2/(iSlat+1) << " " << tempDiff/(iSlat+1) << endl;
+
+
+ Double_t p0l[3] = {0};
+ Double_t p1l[3] = {0};
+ Double_t p2l[3] = {0};
+ Double_t p0g[3] = {0};
+ Double_t p1g[3] = {0};
+ Double_t p2g[3] = {0};
+
+ p0l[2] = lSBTLoc6[iSlat][0][2];
+ // Button targets local coordinates
+ for(Int_t iCor=0; iCor<3; iCor++){
+ p1l[iCor]= lSBTLoc6[iSlat][0][iCor];
+ p2l[iCor]= lSBTLoc6[iSlat][1][iCor];
+ }
+ for(Int_t i=0; i<9; i++){
+ lDiffMin[i]=1000000.;
+ }
+ // Trying 2x*2y*2z*2phi possibilities
+ for(Int_t iX=0; iX<2; iX++){
+ for(Int_t iY=0; iY<2; iY++){
+ for(Int_t iZ=0; iZ<2; iZ++){
+ lCenTemp[0] = fXcSlat[iSlat][iX]->Eval(-pointSBT[0]->GetX(),-pointSBT[1]->GetX());
+ lCenTemp[1] = fYcSlat[iSlat][iY]->Eval(pointSBT[0]->GetZ(),pointSBT[1]->GetZ());
+ lCenTemp[2] = fZcSlat[iSlat][iZ]->Eval(pointSBT[0]->GetY(),pointSBT[1]->GetY());
+ // lCenTemp[2] = fZcSlat[iSlat][iZ]->Eval(pointSSBT[0]->GetY(),pointSSBT[1]->GetY());
+ lRotTemp[0] = psi[iSlat];
+ lRotTemp[1] = tht[iSlat];
+ for(Int_t iP=0; iP<2; iP++){
+ lRotTemp[2] = fPhiXSlat[iSlat][iX][iP]->Eval(-pointSBT[0]->GetX(),-pointSBT[1]->GetX());
+
+ trfTemp.SetTranslation(transTemp);
+ trfTemp.SetRotation(rotTemp);
+ trfTemp.Clear();
+ trfTemp.RotateZ(TMath::RadToDeg()*lRotTemp[2]);
+ trfTemp.RotateY(TMath::RadToDeg()*lRotTemp[1]);
+ trfTemp.RotateX(TMath::RadToDeg()*lRotTemp[0]);
+ trfTemp.SetTranslation(lCenTemp);
+
+ // trfTemp.LocalToMaster(p0l, p0g);
+ // lCenTemp[2]= fSlat[iSlat]->Eval(p0g[0],p0g[1]) - trfTemp.GetRotationMatrix()[8]*p0l[2];
+ // trfTemp.SetTranslation(lCenTemp);
+
+ trfTemp.LocalToMaster(p1l, p1g);
+ trfTemp.LocalToMaster(p2l, p2g);
+
+ lDiffTemp[0] = (-pointSBT[0]->GetX()-p1g[0]);
+ lDiffTemp[1] = (pointSBT[0]->GetZ()-p1g[1]);
+ lDiffTemp[2] = (pointSBT[0]->GetY()-p1g[2]);
+ // lDiffTemp[2] = (pointSSBT[0]->GetY()-p1g[2]);
+ lDiffTemp[3] = TMath::Sqrt(lDiffTemp[0]*lDiffTemp[0]+lDiffTemp[1]*lDiffTemp[1]+lDiffTemp[2]*lDiffTemp[2]);
+ lDiffTemp[4] = (-pointSBT[1]->GetX()-p2g[0]);
+ lDiffTemp[5] = (pointSBT[1]->GetZ()-p2g[1]);
+ lDiffTemp[6] = (pointSBT[1]->GetY()-p2g[2]);
+ // lDiffTemp[6] = (pointSSBT[1]->GetY()-p2g[2]);
+ lDiffTemp[7] = TMath::Sqrt(lDiffTemp[4]*lDiffTemp[4]+lDiffTemp[5]*lDiffTemp[5]+lDiffTemp[6]*lDiffTemp[6]);
+ lDiffTemp[8] = TMath::Sqrt(lDiffTemp[3]*lDiffTemp[3]+lDiffTemp[7]*lDiffTemp[7]);
+
+ if(lDiffTemp[8]<lDiffMin[8]){
+ cout << "Diffs" ;
+ for(Int_t i=0; i<9; i++){
+ cout << " " << lDiffTemp[i];
+ }
+ cout << endl;
+ cout << "Slat" << iSlat+1 << " : mycenX" << iX << "Y" << iY << "Z" << iZ << "(" << lCenTemp[0] << "," << lCenTemp[1] << "," << lCenTemp[2] << "); rotx" << iP << "(" << lRotTemp[0] << "," << lRotTemp[1] << "," << lRotTemp[2] << ")" << endl;
+ cout << p1g[0] << " " << p1g[1] << " " << p1g[2] << " " << p2g[0] << " " << p2g[1] << " " << p2g[2] << endl;
+ cout << "Transformation improved ..." << endl;
+ lCenSlat[iSlat][0] = lCenTemp[0]; lCenSlat[iSlat][1] = lCenTemp[1]; lCenSlat[iSlat][2] = lCenTemp[2];
+ lRotSlat[iSlat][2] = lRotTemp[2];
+ for(Int_t i=0; i<9; i++){
+ lDiffMin[i]=lDiffTemp[i];
+ }
+ if((lDiffMin[3]*lDiffMin[3]<0.1*0.1+0.1*0.1+0.1*0.1)&&
+ (lDiffMin[7]*lDiffMin[7]<0.1*0.1+0.1*0.1+0.1*0.1)){
+ cout << "Correct Transformation found X " << iX << " Y " << iY << " Z " << iZ << " xP " << iP << endl;
+ lCenSlat[iSlat][0] = lCenTemp[0]; lCenSlat[iSlat][1] = lCenTemp[1]; lCenSlat[iSlat][2] = lCenTemp[2];
+ lRotSlat[iSlat][2] = lRotTemp[2];
+ }
+ }
+ }
+ for(Int_t iP=0; iP<2; iP++){
+ lRotTemp[2] = fPhiYSlat[iSlat][iY][iP]->Eval(pointSBT[0]->GetZ(),pointSBT[1]->GetZ());
+ Double_t lPhi = TMath::ATan2((pointSBT[1]->GetZ()-pointSBT[0]->GetZ()),-(pointSBT[1]->GetX()-pointSBT[0]->GetX()));
+
+ trfTemp.Clear();
+ trfTemp.RotateZ(TMath::RadToDeg()*lRotTemp[2]);
+ trfTemp.RotateY(TMath::RadToDeg()*lRotTemp[1]);
+ trfTemp.RotateX(TMath::RadToDeg()*lRotTemp[0]);
+ trfTemp.SetTranslation(lCenTemp);
+
+ // trfTemp.LocalToMaster(p0l, p0g);
+ // lCenTemp[2]= fSlat[iSlat]->Eval(p0g[0],p0g[1]) - trfTemp.GetRotationMatrix()[8]*p0l[2];
+ // trfTemp.SetTranslation(lCenTemp);
+
+ trfTemp.LocalToMaster(p1l, p1g);
+ trfTemp.LocalToMaster(p2l, p2g);
+
+ lDiffTemp[0] = (-pointSBT[0]->GetX()-p1g[0]);
+ lDiffTemp[1] = (pointSBT[0]->GetZ()-p1g[1]);
+ lDiffTemp[2] = (pointSBT[0]->GetY()-p1g[2]);
+ // lDiffTemp[2] = (pointSSBT[0]->GetY()-p1g[2]);
+ lDiffTemp[3] = TMath::Sqrt(lDiffTemp[0]*lDiffTemp[0]+lDiffTemp[1]*lDiffTemp[1]+lDiffTemp[2]*lDiffTemp[2]);
+ lDiffTemp[4] = (-pointSBT[1]->GetX()-p2g[0]);
+ lDiffTemp[5] = (pointSBT[1]->GetZ()-p2g[1]);
+ lDiffTemp[6] = (pointSBT[1]->GetY()-p2g[2]);
+ // lDiffTemp[6] = (pointSSBT[1]->GetY()-p2g[2]);
+ lDiffTemp[7] = TMath::Sqrt(lDiffTemp[4]*lDiffTemp[4]+lDiffTemp[5]*lDiffTemp[5]+lDiffTemp[6]*lDiffTemp[6]);
+ lDiffTemp[8] = TMath::Sqrt(lDiffTemp[3]*lDiffTemp[3]+lDiffTemp[7]*lDiffTemp[7]);
+
+ if(lDiffTemp[8]<lDiffMin[8]){
+ cout << "Diffs" ;
+ for(Int_t i=0; i<9; i++){
+ cout << " " << lDiffTemp[i];
+ }
+ cout << endl;
+ cout << "Slat" << iSlat+1 << " : mycenX" << iX << "Y" << iY << "Z" << iZ << "(" << lCenTemp[0] << "," << lCenTemp[1] << "," << lCenTemp[2] << "); roty" << iP << "(" << lRotTemp[0] << "," << lRotTemp[1] << "," << lRotTemp[2] << "(" << lPhi << "))" << endl;
+ cout << p1g[0] << " " << p1g[1] << " " << p1g[2] << " " << p2g[0] << " " << p2g[1] << " " << p2g[2] << endl;
+ cout << "Transformation improved ..." << endl;
+ lCenSlat[iSlat][0] = lCenTemp[0]; lCenSlat[iSlat][1] = lCenTemp[1]; lCenSlat[iSlat][2] = lCenTemp[2];
+ lRotSlat[iSlat][2] = lRotTemp[2];
+
+ for(Int_t i=0; i<9; i++){
+ lDiffMin[i]=lDiffTemp[i];
+ }
+ if((lDiffMin[3]*lDiffMin[3]<0.1*0.1+0.1*0.1+0.1*0.1)&&
+ (lDiffMin[7]*lDiffMin[7]<0.1*0.1+0.1*0.1+0.1*0.1)){
+ cout << "Correct Transformation found X " << iX << " Y " << iY << " Z " << iZ << " yP " << iP << endl;
+ lCenSlat[iSlat][0] = lCenTemp[0]; lCenSlat[iSlat][1] = lCenTemp[1]; lCenSlat[iSlat][2] = lCenTemp[2];
+ lRotSlat[iSlat][2] = lRotTemp[2];
+ }
+ }
+ }
+ }
+ }
+ }
+
+ // Fill slat plane for fit monitor.
+ xMinSlat = hSSTrpy->GetXaxis()->FindBin(pointSBT[0]->GetX());
+ xMaxSlat = hSSTrpy->GetXaxis()->FindBin(pointSBT[1]->GetX());
+ yMinSlat = hSSTrpy->GetYaxis()->FindBin(pointSBT[0]->GetZ()-200.);
+ yMaxSlat = hSSTrpy->GetYaxis()->FindBin(pointSBT[0]->GetZ()+200.);
+
+ for (int i=(int)xMinSlat; i<=(int)xMaxSlat; i++) {
+ for (int j=(int)yMinSlat; j<=(int)yMaxSlat; j++) {
+ Double_t zSlat = fSlat[iSlat]->Eval(-hSSTrpy->GetXaxis()->GetBinCenter(i),
+ hSSTrpy->GetYaxis()->GetBinCenter(j));
+ if((iSlat+1)%2==0){
+ hSSTlpy->SetBinContent(i,j,zSlat);
+ }
+ else {
+ hSSTrpy->SetBinContent(i,j,zSlat);
+ }
+ }
+ }
+ }
+
+ //
+ // Compare transformations to expected ones
+ //
+ Int_t iSlatToPos[13] = {0, 11, 2, 9, 4, 7, 6, 5, 8, 3, 10, 1, 12};
+
+ // Theoretical differences with respect to Slat 513 which is here Slat07
+ lDiffThCenSlat0[0][0] = (2-5)*400./2. +12.5 -375. -25.;
+ lDiffThCenSlat0[1][0] = (3-5)*400./2. +12.5 -375. -25.;
+ lDiffThCenSlat0[2][0] = (4-5)*400./2. +12.5 -375. -25.;
+ lDiffThCenSlat0[3][0] = (5-5)*400./2. +12.5 -375. -25.;
+ lDiffThCenSlat0[4][0] = (5-5)*400./2. +12.5 -375. -25.;
+ lDiffThCenSlat0[5][0] = (6-5)*400./2. +12.5 -375. -25.;
+ lDiffThCenSlat0[6][0] = (5-5)*400./2. +12.5 -375. -25. +375. +25. -12.5;
+ lDiffThCenSlat0[7][0] = (6-5)*400./2. +12.5 -375. -25.;
+ lDiffThCenSlat0[8][0] = (5-5)*400./2. +12.5 -375. -25.;
+ lDiffThCenSlat0[9][0] = (5-5)*400./2. +12.5 -375. -25.;
+ lDiffThCenSlat0[10][0] = (4-5)*400./2. +12.5 -375. -25.;
+ lDiffThCenSlat0[11][0] = (3-5)*400./2. +12.5 -375. -25.;
+ lDiffThCenSlat0[12][0] = (2-5)*400./2. +12.5 -375. -25.;
+
+ lDiffThCenSlat0[12][1] = 382.0 +378.5 +375.5 +294.0 +370.0 +286.0;
+ lDiffThCenSlat0[1][1] = 382.0 +378.5 +375.5 +294.0 +370.0;
+ lDiffThCenSlat0[10][1] = 382.0 +378.5 +375.5 +294.0;
+ lDiffThCenSlat0[3][1] = 382.0 +378.5 +375.5;
+ lDiffThCenSlat0[8][1] = 382.0 +378.5;
+ lDiffThCenSlat0[5][1] = 382.0;
+ lDiffThCenSlat0[6][1] = 0.0;
+ lDiffThCenSlat0[7][1] = -382.0;
+ lDiffThCenSlat0[4][1] = -382.0 -378.5;
+ lDiffThCenSlat0[9][1] = -382.0 -378.5 -375.5;
+ lDiffThCenSlat0[2][1] = -382.0 -378.5 -375.5 -294.0;
+ lDiffThCenSlat0[11][1] = -382.0 -378.5 -375.5 -294.0 -370.0;
+ lDiffThCenSlat0[0][1] = -382.0 -378.5 -375.5 -294.0 -370.0 -286.0;
+
+ lDiffThCenSlat0[12][2] = (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[1][2] = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[10][2] = (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[3][2] = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[8][2] = (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[5][2] = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[6][2] = (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[7][2] = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[4][2] = (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[9][2] = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[2][2] = (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[11][2] = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
+ lDiffThCenSlat0[0][2] = (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
+
+ TGraph *gDeltaDiffCenXSlat0 = new TGraph(nSlats);
+ TGraph *gDeltaDiffCenYSlat0 = new TGraph(nSlats);
+ TGraph *gDeltaDiffCenZSlat0 = new TGraph(nSlats);
+ TGraph *gDeltaDiffPsiSlat0 = new TGraph(nSlats);
+ TGraph *gDeltaDiffThtSlat0 = new TGraph(nSlats);
+ TGraph *gDeltaDiffPhiSlat0 = new TGraph(nSlats);
+
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ trfSlat[iSlat].SetTranslation(transSlat[iSlat]);
+ trfSlat[iSlat].SetRotation(rotSlat[iSlat]);
+ trfSlat[iSlat].Clear();
+ trfSlat[iSlat].RotateZ(TMath::RadToDeg()*lRotSlat[iSlat][2]);
+ trfSlat[iSlat].RotateY(TMath::RadToDeg()*lRotSlat[iSlat][1]);
+ trfSlat[iSlat].RotateX(TMath::RadToDeg()*lRotSlat[iSlat][0]);
+ trfSlat[iSlat].SetTranslation(lCenSlat[iSlat]);
+ }
+
+ Double_t *myTrans = 0;
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ dtrfSlat[iSlat].SetTranslation(dtransSlat[iSlat]);
+ dtrfSlat[iSlat].SetRotation(drotSlat[iSlat]);
+ dtrfSlat[iSlat].Clear();
+ dtrfSlat[iSlat] = trfSlat[6].Inverse()*trfSlat[iSlat];
+ dtrfSlat[iSlat].Print();
+ lDiffCenSlat0[iSlat] = (Double_t*)dtrfSlat[iSlat].GetTranslation();
+ MatrixToAngles(dtrfSlat[iSlat].GetRotationMatrix(),lDiffRotSlat0[iSlat]);
+
+ lDeltaDiffCenSlat0[iSlat][0] = lDiffCenSlat0[iSlat][0]-lDiffThCenSlat0[iSlat][0];
+ lDeltaDiffCenSlat0[iSlat][1] = -lDiffCenSlat0[iSlat][1]-lDiffThCenSlat0[iSlat][1];
+ lDeltaDiffCenSlat0[iSlat][2] = -lDiffCenSlat0[iSlat][2]-lDiffThCenSlat0[iSlat][2];
+ lDeltaDiffRotSlat0[iSlat][0] = lDiffRotSlat0[iSlat][0];
+ lDeltaDiffRotSlat0[iSlat][1] = lDiffRotSlat0[iSlat][1];
+ lDeltaDiffRotSlat0[iSlat][2] = lDiffRotSlat0[iSlat][2];
+ gDeltaDiffCenXSlat0->SetPoint(iSlat,lDeltaDiffCenSlat0[iSlat][0],iSlatToPos[iSlat]+1);
+ gDeltaDiffCenYSlat0->SetPoint(iSlat,lDeltaDiffCenSlat0[iSlat][1],iSlatToPos[iSlat]+1);
+ gDeltaDiffCenZSlat0->SetPoint(iSlat,lDeltaDiffCenSlat0[iSlat][2],iSlatToPos[iSlat]+1);
+ gDeltaDiffPsiSlat0->SetPoint(iSlat,1e3*lDeltaDiffRotSlat0[iSlat][0],iSlatToPos[iSlat]+1);
+ gDeltaDiffThtSlat0->SetPoint(iSlat,1e3*lDeltaDiffRotSlat0[iSlat][1],iSlatToPos[iSlat]+1);
+ gDeltaDiffPhiSlat0->SetPoint(iSlat,1e3*lDeltaDiffRotSlat0[iSlat][2],iSlatToPos[iSlat]+1);
+ }
+
+ // Import TGeo geometry
+ char* geoFilename = "geometry.root";
+ cout << "geometry imported" << endl;
+ if ( ! AliGeomManager::GetGeometry() ) {
+ AliGeomManager::LoadGeometry(geoFilename);
+ if (! AliGeomManager::GetGeometry() ) {
+ printf("MUONSurveyCh8L: getting geometry from file %s failed\n", geoFilename);
+ return;
+ }
+ }
+
+ AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer();
+// transform->ReadGeometryData("volpath.dat", gGeoManager);
+ transform->LoadGeometryData();
+ cout << "geometry data read" << endl;
+ AliMUONGeometryTransformer *newTransform = ReAlign(transform,11,dtrfSlat,true);
+ newTransform->WriteTransformations("transform2ReAlign.dat");
+ cout << "newtransform read" << endl;
+ // Generate realigned data in local cdb
+ const TClonesArray* array = newTransform->GetMisAlignmentData();
+
+ // CDB manager
+ AliCDBManager* cdbManager = AliCDBManager::Instance();
+ cdbManager->SetDefaultStorage("local://ReAlignCDB");
+
+ AliCDBMetaData* cdbData = new AliCDBMetaData();
+ cdbData->SetResponsible("Dimuon Offline project");
+ cdbData->SetComment("MUON alignment objects with residual misalignment");
+ AliCDBId id("MUON/Align/Data", 0, 0);
+ cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
+
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ cout << lDeltaDiffCenSlat0[iSlat][0] << " " << lDiffCenSlat0[iSlat][0] << " " << lDiffThCenSlat0[iSlat][0] << endl;
+ }
+ cout << endl;
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ // lDiffCenSlat0[iSlat][1] = lCenSlat[iSlat][1]-lCenSlat[6][1];
+ // lDeltaDiffCenSlat0[iSlat][1] = lDiffCenSlat0[iSlat][1]-lDiffThCenSlat0[iSlat][1];
+ cout << lDeltaDiffCenSlat0[iSlat][1] << " " << lDiffCenSlat0[iSlat][1] << " " << lDiffThCenSlat0[iSlat][1] << endl;
+ }
+ cout << endl;
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ cout << lDeltaDiffCenSlat0[iSlat][2] << " " << lDiffCenSlat0[iSlat][2] << " " << lDiffThCenSlat0[iSlat][2] << endl;
+ }
+ cout << endl;
+ cout << endl;
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ cout << lDeltaDiffRotSlat0[iSlat][0] << " " << lDiffRotSlat0[iSlat][0] << " " << lDiffThRotSlat0[iSlat][0] << endl;
+ }
+ cout << endl;
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ cout << lDeltaDiffRotSlat0[iSlat][1] << " " << lDiffRotSlat0[iSlat][1] << " " << lDiffThRotSlat0[iSlat][1] << endl;
+ }
+ cout << endl;
+ for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
+ cout << lDeltaDiffRotSlat0[iSlat][2] << " " << lDiffRotSlat0[iSlat][2] << " " << lDiffThRotSlat0[iSlat][2] << endl;
+ }
+
+ TH1F *mySlatDeltaDiffCenX = new TH1F("mySlatDeltaDiffCenX","mySlatDeltaDiffCenX",100,-10,10);
+ mySlatDeltaDiffCenX->SetMaximum(15);
+ mySlatDeltaDiffCenX->SetMinimum(0);
+ TH1F *mySlatDeltaDiffCenY = new TH1F("mySlatDeltaDiffCenY","mySlatDeltaDiffCenY",100,-10,10);
+ mySlatDeltaDiffCenY->SetMaximum(15);
+ mySlatDeltaDiffCenY->SetMinimum(0);
+ TH1F *mySlatDeltaDiffCenZ = new TH1F("mySlatDeltaDiffCenZ","mySlatDeltaDiffCenZ",100,-20,20);
+ mySlatDeltaDiffCenZ->SetMaximum(15);
+ mySlatDeltaDiffCenZ->SetMinimum(0);
+
+ TH1F *mySlatDeltaDiffRotX = new TH1F("mySlatDeltaDiffRotX","mySlatDeltaDiffRotX",100,-10,10);
+ mySlatDeltaDiffRotX->SetMaximum(15);
+ mySlatDeltaDiffRotX->SetMinimum(0);
+ TH1F *mySlatDeltaDiffRotY = new TH1F("mySlatDeltaDiffRotY","mySlatDeltaDiffRotY",100,-10,10);
+ mySlatDeltaDiffRotY->SetMaximum(15);
+ mySlatDeltaDiffRotY->SetMinimum(0);
+ TH1F *mySlatDeltaDiffRotZ = new TH1F("mySlatDeltaDiffRotZ","mySlatDeltaDiffRotZ",100,-5,5);
+ mySlatDeltaDiffRotZ->SetMaximum(15);
+ mySlatDeltaDiffRotZ->SetMinimum(0);
+ //
+ // ******** Starting plots
+ //
+ TCanvas *canvas;
+ TPad *pad;
+ TPaveLabel *theTitle;
+
+ TPostScript *ps = 0;
+
+ if( saveps ){
+ ps = new TPostScript(filename,filetype);
+ ps->NewPage();
+ }
+
+ // Inv Mass, Multiplicity
+ sprintf(str,"Chamber 8L");
+ TCanvas *cvn0 = new TCanvas("cvn0",str,cWidth,cHeight);
+ canvas = cvn0;
+ canvas->Range(0,0,21,29);
+
+ TPaveLabel *theTitle0 = new TPaveLabel(3,27.0,18,28.5," Deformations of chamber 8L ","br");
+ theTitle = theTitle0;
+ theTitle->SetFillColor(18);
+ theTitle->SetTextFont(32);
+ theTitle->SetTextSize(0.4);
+ theTitle->SetTextColor(1);
+ theTitle->Draw();
+
+ TPad *pad0 = new TPad("pad0","pad0",0.01,0.01,0.98,0.91,0);
+ pad = pad0;
+ pad->Draw();
+ pad->Divide(2,2);
+
+ pad->cd(1);
+ gStyle->SetPalette(1);
+ hCPSTry->SetMinimum(100);
+ hCPSTry->SetMaximum(120);
+ hCPSTry->Draw("lego2z");
+
+ pad->cd(2);
+ gStyle->SetPalette(1);
+ hSSTry->SetMinimum(60);
+ hSSTry->SetMaximum(80);
+ hSSTry->Draw("lego2z");
+
+ pad->cd(3);
+ gStyle->SetPalette(1);
+ hCSTy->SetMinimum(110);
+ hCSTy->SetMaximum(130);
+ hCSTy->Draw("lego2z");
+
+ pad->cd(4);
+ gStyle->SetPalette(1);
+ hSSTly->SetMinimum(165);
+ hSSTly->SetMaximum(185);
+ hSSTly->Draw("lego2z");
+
+ if(saveps){
+ ps->NewPage();
+ }
+
+ // Inv Mass, Multiplicity
+ sprintf(str,"Chamber 8L");
+ TCanvas *cvn1 = new TCanvas("cvn1",str,cWidth,cHeight);
+ canvas = cvn1;
+ canvas->Range(0,0,21,29);
+
+ TPaveLabel *theTitle1 = new TPaveLabel(3,27.0,18,28.5," Deformations of chamber 8L ","br");
+ theTitle = theTitle1;
+ theTitle->SetFillColor(18);
+ theTitle->SetTextFont(32);
+ theTitle->SetTextSize(0.4);
+ theTitle->SetTextColor(1);
+ theTitle->Draw();
+
+ TPad *pad1 = new TPad("pad1","pad1",0.01,0.01,0.98,0.91,0);
+ pad = pad1;
+ pad->Draw();
+ pad->Divide(2,2);
+
+ pad->cd(1);
+ gStyle->SetPalette(1);
+ hCPSTly->SetMinimum(120);
+ hCPSTly->SetMaximum(140);
+ hCPSTly->Draw("lego2z");
+
+ pad->cd(2);
+ gStyle->SetPalette(1);
+ hSSTrpy->SetMinimum(60);
+ hSSTrpy->SetMaximum(80);
+ hSSTrpy->Draw("surf2z");
+
+ pad->cd(3);
+ gStyle->SetPalette(1);
+ hCSTy->SetMinimum(110);
+ hCSTy->SetMaximum(130);
+ hCSTy->Draw("lego2z");
+
+ pad->cd(4);
+ gStyle->SetPalette(1);
+ hSSTlpy->SetMinimum(165);
+ hSSTlpy->SetMaximum(185);
+ hSSTlpy->Draw("surf2z");
+
+ // Inv Mass, Multiplicity
+ sprintf(str,"Chamber 8L");
+ TCanvas *cvn2 = new TCanvas("cvn2",str,cWidth,cHeight);
+ canvas = cvn2;
+ canvas->Range(0,0,21,29);
+
+ TPaveLabel *theTitle2 = new TPaveLabel(3,27.0,18,28.5," Deformations of chamber 8L ","br");
+ theTitle = theTitle2;
+ theTitle->SetFillColor(18);
+ theTitle->SetTextFont(32);
+ theTitle->SetTextSize(0.4);
+ theTitle->SetTextColor(1);
+ theTitle->Draw();
+
+ TPad *pad2 = new TPad("pad2","pad2",0.01,0.01,0.98,0.91,0);
+ pad = pad2;
+ pad->Draw();
+ pad->Divide(3,2);
+
+ pad->cd(1);
+ mySlatDeltaDiffCenX->Draw();
+ mySlatDeltaDiffCenX->SetXTitle("#Delta[(xc_{i}^{m}-xc_{0}^{m})-(xc_{i}^{th}-xc_{0}^{th})] (mm)");
+ mySlatDeltaDiffCenX->SetYTitle("Slat ascending vertical ordering");
+ gDeltaDiffCenXSlat0->SetMarkerStyle(20);
+ gDeltaDiffCenXSlat0->Draw("P");
+
+ pad->cd(2);
+ mySlatDeltaDiffCenY->Draw();
+ mySlatDeltaDiffCenY->SetXTitle("#Delta[(yc_{i}^{m}-yc_{0}^{m})-(yc_{i}^{th}-yc_{0}^{th})] (mm)");
+ mySlatDeltaDiffCenY->SetYTitle("Slat ascending vertical ordering");
+ gDeltaDiffCenYSlat0->SetMarkerStyle(20);
+ gDeltaDiffCenYSlat0->Draw("P");
+
+ pad->cd(3);
+ mySlatDeltaDiffCenZ->Draw();
+ mySlatDeltaDiffCenZ->SetXTitle("#Delta[(zc_{i}^{m}-zc_{0}^{m})-(zc_{i}^{th}-zc_{0}^{th})] (mm)");
+ mySlatDeltaDiffCenZ->SetYTitle("Slat ascending vertical ordering");
+ gDeltaDiffCenZSlat0->SetMarkerStyle(20);
+ gDeltaDiffCenZSlat0->Draw("P");
+
+ pad->cd(4);
+ mySlatDeltaDiffRotX->Draw();
+ mySlatDeltaDiffRotX->SetXTitle("#Delta[(#psi_{i}^{m}-#psi_{0}^{m})-(#psi_{i}^{th}-#psi_{0}^{th})] (mrad)");
+ mySlatDeltaDiffRotX->SetYTitle("Slat ascending vertical ordering");
+ gDeltaDiffPsiSlat0->SetMarkerStyle(20);
+ gDeltaDiffPsiSlat0->Draw("P");
+
+ pad->cd(5);
+ mySlatDeltaDiffRotY->Draw();
+ mySlatDeltaDiffRotY->SetXTitle("#Delta[(#theta_{i}^{m}-#theta_{0}^{m})-(#theta_{i}^{th}-#theta_{0}^{th})] (mrad)");
+ mySlatDeltaDiffRotY->SetYTitle("Slat ascending vertical ordering");
+ gDeltaDiffThtSlat0->SetMarkerStyle(20);
+ gDeltaDiffThtSlat0->Draw("P");
+
+ pad->cd(6);
+ mySlatDeltaDiffRotZ->Draw();
+ mySlatDeltaDiffRotZ->SetXTitle("#Delta[(#phi_{i}^{m}-#phi_{0}^{m})-(#phi_{i}^{th}-#phi_{0}^{th})] (mrad)");
+ mySlatDeltaDiffRotZ->SetYTitle("Slat ascending vertical ordering");
+ gDeltaDiffPhiSlat0->SetMarkerStyle(20);
+ gDeltaDiffPhiSlat0->Draw("P");
+
+ if( saveps ){
+ ps->Close();
+ }
+}
--- /dev/null
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * *
+ * Author: The ALICE Off-line Project. *
+ * Contributors are mentioned in the code where appropriate. *
+ * *
+ * Permission to use, copy, modify and distribute this software and its *
+ * documentation strictly for non-commercial purposes is hereby granted *
+ * without fee, provided that the above copyright notice appears in all *
+ * copies and that both the copyright notice and this permission notice *
+ * appear in the supporting documentation. The authors make no claims *
+ * about the suitability of this software for any purpose. It is *
+ * provided "as is" without express or implied warranty. *
+ **************************************************************************/
+
+/* $Id$ */
+
+// ---
+// Utility macro for survey data to alignment transformation.
+//
+// Macro contains various functions to calculate misalignement parameters
+// from survey data and designed positions of survey targets.
+// Macro also includes a method to get the new AliMUONGeometryTranformer.
+// It is intended to be loaded by chamber specific macros.
+//
+// Author: Javier Castillo
+// ---
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+
+#include "AliMUONGeometryTransformer.h"
+#include "AliMUONGeometryMisAligner.h"
+#include "AliMUONGeometryModuleTransformer.h"
+#include "AliMUONGeometryDetElement.h"
+#include "AliMUONGeometryBuilder.h"
+#include "AliMpExMap.h"
+
+#include "AliGeomManager.h"
+#include "AliCDBManager.h"
+#include "AliCDBMetaData.h"
+#include "AliCDBId.h"
+
+#include <TGeoManager.h>
+#include <TClonesArray.h>
+#include <TMath.h>
+#include <TString.h>
+#include <Riostream.h>
+
+#include <fstream>
+
+#endif
+
+static int fgNDetElemCh[10] = {4,4,4,4,18,18,26,26,26,26};
+
+Bool_t MatrixToAngles(const Double_t *rot, Double_t *angles)
+{
+ // Calculates the Euler angles in "x y z" notation
+ // using the rotation matrix
+ // Returns false in case the rotation angles can not be
+
+ // extracted from the matrix
+ //
+ if(TMath::Abs(rot[0])<1e-7 || TMath::Abs(rot[8])<1e-7) {
+ printf("Failed to extract roll-pitch-yall angles!");
+ return kFALSE;
+ }
+ // Double_t raddeg = TMath::RadToDeg();
+ angles[0]=TMath::ATan2(-rot[5],rot[8]);
+ angles[1]=TMath::ASin(rot[2]);
+ angles[2]=TMath::ATan2(-rot[1],rot[0]);
+ return kTRUE;
+}
+
+Double_t eqPlane(Double_t *x, Double_t *par){
+ return (-par[0]*x[0] -par[1]*x[1] -par[2]);
+}
+
+Double_t xpCenter(Double_t *x, Double_t *par){
+
+ Double_t lCos2Tht = TMath::Cos(2*par[6]);
+ Double_t lSinTht = TMath::Sin(par[6]);
+
+ Double_t inSqrt = TMath::Abs((par[0] - par[3])*(par[0] - par[3])
+ -2*(x[0] -x[1])*(x[0] -x[1])
+ +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5])
+ +((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4])
+ +(par[2] - par[5])*(par[2] - par[5]))*lCos2Tht
+ +4*(x[0] - x[1])*(par[2] - par[5])*lSinTht);
+
+ Double_t xD = ((2*(par[0]*par[0]*x[1]
+ -par[0]*par[3]*(x[0] + x[1])
+ +x[1]*par[1]*(par[1] - par[4])
+ +x[0]*(par[3]*par[3] - par[1]*par[4] + par[4]*par[4]))
+ -2*(par[3]*par[3]*par[2]
+ +par[0]*par[0]*par[5]
+ -par[0]*par[3]*(par[2] + par[5])
+ +(par[1] - par[4])*(-par[4]*par[2] +par[1]*par[5]))*lSinTht
+ +TMath::Sqrt(2)*(-par[3]*par[1] + par[0]*par[4])
+ *TMath::Sqrt(inSqrt))
+ /(2*((par[0] - par[3])*(par[0] - par[3]) + (par[1] - par[4])*(par[1] - par[4]))));
+
+ return xD;
+}
+
+Double_t xnCenter(Double_t *x, Double_t *par){
+
+ Double_t lCos2Tht = TMath::Cos(2*par[6]);
+ Double_t lSinTht = TMath::Sin(par[6]);
+
+ Double_t inSqrt = TMath::Abs((par[0] - par[3])*(par[0] - par[3])
+ -2*(x[0] - x[1])*(x[0] - x[1])
+ +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5])
+ +((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4])
+ +(par[2] - par[5])*(par[2] - par[5]))*lCos2Tht
+ +4*(x[0] - x[1])*(par[2] - par[5])*lSinTht);
+
+ Double_t xD = ((2*(par[0]*par[0]*x[1]
+ -par[0]*par[3]*(x[0] + x[1])
+ +x[1]*par[1]*(par[1] - par[4])
+ +x[0]*(par[3]*par[3] - par[1]*par[4] + par[4]*par[4]))
+ -2*(par[3]*par[3]*par[2] + par[0]*par[0]*par[5]
+ -par[0]*par[3]*(par[2] + par[5])
+ +(par[1] - par[4])*(-par[4]*par[2] + par[1]*par[5]))*lSinTht
+ +TMath::Sqrt(2)*(par[3]*par[1] - par[0]*par[4])
+ *TMath::Sqrt(inSqrt))
+ /(2*((par[0] - par[3])*(par[0] - par[3]) + (par[1] - par[4])*(par[1] - par[4]))));
+
+ return xD;
+}
+
+Double_t phixpn(Double_t *x, Double_t *par){
+
+ Double_t inSqrt = TMath::Abs(((par[0] - par[3])*(par[0] - par[3])
+ -2*(x[0] - x[1])*(x[0] - x[1])
+ +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5])
+ +(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4])
+ +(par[2] - par[5])*(par[2] - par[5]))*TMath::Cos(2*par[6])
+ +4*(x[0] - x[1])*(par[2] - par[5])*TMath::Sin(par[6])));
+
+ Double_t phix = ((+2*(par[0] - par[3])*(x[0] - x[1])
+ -2*(par[0] - par[3])*(par[2] - par[5])*TMath::Sin(par[6])
+ +TMath::Sqrt(2)*(par[1] - par[4])
+ *TMath::Sqrt(inSqrt))
+ /(2*(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*TMath::Cos(par[6])));
+
+ phix = -TMath::ACos(phix);
+
+ return phix;
+}
+
+Double_t phixpp(Double_t *x, Double_t *par){
+
+ Double_t inSqrt = TMath::Abs(+(par[0] - par[3])*(par[0] - par[3])
+ -2*(x[0] - x[1])*(x[0] - x[1])
+ +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5])
+ +(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4])
+ +(par[2] - par[5])*(par[2] - par[5]))*TMath::Cos(2*par[6])
+ +4*(x[0] - x[1])*(par[2] - par[5])*TMath::Sin(par[6]));
+
+ Double_t phix = ((+2*(par[0] - par[3])*(x[0] - x[1])
+ -2*(par[0] - par[3])*(par[2] - par[5])*TMath::Sin(par[6])
+ +TMath::Sqrt(2)*(par[1] - par[4])
+ *TMath::Sqrt(inSqrt))
+ /(2*(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*TMath::Cos(par[6])));
+
+ phix = TMath::ACos(phix);
+
+ return phix;
+}
+
+Double_t phixnn(Double_t *x, Double_t *par){
+
+ Double_t inSqrt = TMath::Abs(+(par[0] - par[3])*(par[0] - par[3])
+ -2*(x[0] - x[1])*(x[0] - x[1])
+ +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5])
+ +(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4])
+ +(par[2] - par[5])*(par[2] - par[5]))*TMath::Cos(2*par[6])
+ + 4*(x[0] - x[1])*(par[2] - par[5])*TMath::Sin(par[6]));
+
+ Double_t phix = (+(+2*(par[0] - par[3])*(x[0] - x[1])
+ -2*(par[0] - par[3])*(par[2] - par[5])*TMath::Sin(par[6])
+ +TMath::Sqrt(2)*(-par[1] + par[4])
+ *TMath::Sqrt(inSqrt))
+ /(2*(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*TMath::Cos(par[6])));
+
+ phix = -TMath::ACos(phix);
+
+ return phix;
+}
+
+Double_t phixnp(Double_t *x, Double_t *par){
+
+ Double_t inSqrt = TMath::Abs(+(par[0] - par[3])*(par[0] - par[3])
+ -2*(x[0] - x[1])*(x[0] - x[1])
+ +(par[1] - par[4] + par[2] - par[5])*(par[1] - par[4] - par[2] + par[5])
+ +(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4])
+ +(par[2] - par[5])*(par[2] - par[5]))*TMath::Cos(2*par[6])
+ +4*(x[0] - x[1])*(par[2] - par[5])*TMath::Sin(par[6]));
+
+ Double_t phix = (+(+2*(par[0] - par[3])*(x[0] - x[1])
+ -2*(par[0] - par[3])*(par[2] - par[5])*TMath::Sin(par[6])
+ +TMath::Sqrt(2)*(-par[1] + par[4])
+ *TMath::Sqrt(inSqrt))
+ /(2*(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*TMath::Cos(par[6])));
+
+ phix = TMath::ACos(phix);
+
+ return phix;
+}
+
+Double_t ypCenter(Double_t *x, Double_t *par){
+ // par : x1l, y1l, z1l, x2l, y2l, z2, lpsi, tht,
+
+ Double_t lCosPsi = TMath::Cos(par[6]);
+ Double_t lSinPsi = TMath::Sin(par[6]);
+ Double_t lCosTht = TMath::Cos(par[7]);
+ Double_t lSinTht = TMath::Sin(par[7]);
+
+ Double_t yD = ((1./((par[0] - par[3])*(par[0] - par[3]) + (par[1] - par[4])*(par[1] - par[4])))
+ *(+par[3]*par[3]*x[0]
+ +par[0]*par[0]*x[1]
+ -par[0]*par[3]*(x[0] + x[1])
+ +(par[1] - par[4])*(-x[0]*par[4] + par[1]*x[1])
+ +(par[3]*par[3]*par[2]
+ +par[0]*par[0]*par[5]
+ -par[0]*par[3]*(par[2] + par[5])
+ +(par[1] - par[4])*(-par[4]*par[2] + par[1]*par[5]))*lCosTht*lSinPsi
+ +(-par[3]*par[1] + par[0]*par[4])
+ *TMath::Sqrt(-(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ *(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ + ((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht))));
+
+ return yD;
+}
+
+Double_t phiypn(Double_t *x, Double_t *par){
+
+ Double_t lCosPsi = TMath::Cos(par[6]);
+ Double_t lSinPsi = TMath::Sin(par[6]);
+ Double_t lCosTht = TMath::Cos(par[7]);
+ Double_t lSinTht = TMath::Sin(par[7]);
+
+ Double_t phiy = ((lCosPsi*((par[1] - par[4])*(x[0] - x[1])
+ +(par[1] - par[4])*(par[2] - par[5])*lCosTht*lSinPsi
+ +(-par[0] + par[3])
+ *TMath::Sqrt(-(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ *(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ +(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht)))
+ +lSinPsi*lSinTht*((par[0] - par[3])*(x[0] - x[1])
+ +(par[0] - par[3])*(par[2] - par[5])*lCosTht*lSinPsi
+ +(par[1] - par[4])
+ *TMath::Sqrt(-(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ *(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ + ((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht))))
+ /((+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht)));
+
+ phiy = -TMath::ACos(phiy);
+
+
+ return phiy;
+}
+
+Double_t phiypp(Double_t *x, Double_t *par){
+
+ Double_t lCosPsi = TMath::Cos(par[6]);
+ Double_t lSinPsi = TMath::Sin(par[6]);
+ Double_t lCosTht = TMath::Cos(par[7]);
+ Double_t lSinTht = TMath::Sin(par[7]);
+
+ Double_t phiy = ((lCosPsi*((par[1] - par[4])*(x[0] - x[1])
+ +(par[1] - par[4])*(par[2] - par[5])*lCosTht*lSinPsi
+ +(-par[0] + par[3])
+ *TMath::Sqrt(-(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ *(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ +((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht)))
+ +lSinPsi*lSinTht*((par[0] - par[3])*(x[0] - x[1])
+ +(par[0] - par[3])*(par[2] - par[5])*lCosTht*lSinPsi
+ +(par[1] - par[4])*TMath::Sqrt(-(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ *(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ +((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht))))
+ /(((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht)));
+
+ phiy = TMath::ACos(phiy);
+
+ return phiy;
+}
+
+Double_t ynCenter(Double_t *x, Double_t *par){
+
+ Double_t lCosPsi = TMath::Cos(par[6]);
+ Double_t lSinPsi = TMath::Sin(par[6]);
+ Double_t lCosTht = TMath::Cos(par[7]);
+ Double_t lSinTht = TMath::Sin(par[7]);
+
+ Double_t yD = ((1./(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4])))
+ *(+par[3]*par[3]*x[0]
+ +par[0]*par[0]*x[1]
+ -par[0]*par[3]*(x[0] + x[1])
+ +(par[1] - par[4])*(-x[0]*par[4] + par[1]*x[1])
+ +(+par[3]*par[3]*par[2]
+ +par[0]*par[0]*par[5]
+ -par[0]*par[3]*(par[2] + par[5])
+ +(par[1] - par[4])*(-par[4]*par[2] + par[1]*par[5]))*lCosTht*lSinPsi
+ +(par[3]*par[1] - par[0]*par[4])
+ *TMath::Sqrt(-(+x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ *(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ +((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht))));
+
+ return yD;
+}
+
+
+Double_t phiynn(Double_t *x, Double_t *par){
+
+ Double_t lCosPsi = TMath::Cos(par[6]);
+ Double_t lSinPsi = TMath::Sin(par[6]);
+ Double_t lCosTht = TMath::Cos(par[7]);
+ Double_t lSinTht = TMath::Sin(par[7]);
+
+ Double_t phiy = ((lCosPsi*(+(par[1] - par[4])*(x[0] - x[1])
+ +(par[1] - par[4])*(par[2] - par[5])*lCosTht*lSinPsi
+ +(par[0] - par[3])
+ *TMath::Sqrt(-(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ *(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ +(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht)))
+ +lSinPsi*lSinTht*(+(par[0] - par[3])*(x[0] - x[1])
+ +(par[0] - par[3])*(par[2] - par[5])*lCosTht*lSinPsi
+ +(-par[1] + par[4])
+ *TMath::Sqrt(-(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ *(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ +(+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht))))
+ /((+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht)));
+
+ phiy = -TMath::ACos(phiy);
+
+ return phiy;
+}
+
+
+Double_t phiynp(Double_t *x, Double_t *par){
+
+ Double_t lCosPsi = TMath::Cos(par[6]);
+ Double_t lSinPsi = TMath::Sin(par[6]);
+ Double_t lCosTht = TMath::Cos(par[7]);
+ Double_t lSinTht = TMath::Sin(par[7]);
+
+ Double_t phiy = ((lCosPsi*(+(par[1] - par[4])*(x[0] - x[1])
+ +(par[1] - par[4])*(par[2] - par[5])*lCosTht*lSinPsi
+ +(par[0] - par[3])
+ *TMath::Sqrt(-(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ *(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ +((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht)))
+ +lSinPsi*lSinTht*(+(par[0] - par[3])*(x[0] - x[1])
+ +(par[0] - par[3])*(par[2] - par[5])*lCosTht*lSinPsi
+ +(-par[1] + par[4])
+ *TMath::Sqrt(-(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ *(x[0] - x[1]
+ +(par[2] - par[5])*lCosTht*lSinPsi)
+ +((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht))))
+ /((+(par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*(+lCosPsi*lCosPsi
+ +lSinPsi*lSinPsi*lSinTht*lSinTht)));
+
+ phiy = TMath::ACos(phiy);
+
+ return phiy;
+}
+
+Double_t znCenter(Double_t *x, Double_t *par){
+ // par : x1l, y1l, z1l, x2l, y2l, z2l, psi, tht
+
+ Double_t lCosPsi = TMath::Cos(par[6]);
+ Double_t lSinPsi = TMath::Sin(par[6]);
+ Double_t lCosTht = TMath::Cos(par[7]);
+ Double_t lSinTht = TMath::Sin(par[7]);
+
+ Double_t inSqrt = ((par[3]*par[1] - par[0]*par[4])*(par[3]*par[1] - par[0]*par[4])
+ *((-(x[0] - x[1])*(x[0] - x[1]))
+ +(((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4])))*lSinPsi*lSinPsi
+ +lCosPsi*((-(par[2] - par[5]))
+ *lCosTht*(-2*x[0]+2*x[1]
+ +(par[2] - par[5])*lCosPsi*lCosTht)
+ +((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*lCosPsi*lSinTht*lSinTht)));
+
+ Double_t zD = ((1./((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4])))
+ *(-par[1]*par[4]*x[0]
+ +par[4]*par[4]*x[0]
+ +par[0]*par[0]*x[1]
+ +par[1]*par[1]*x[1]
+ -par[1]*par[4]*x[1]
+ -par[0]*par[3]*(x[0] + x[1])
+ +par[3]*par[3]*x[0]
+ +(+par[1]*par[4]*par[2]
+ -par[4]*par[4]*par[2]
+ -par[0]*par[0]*par[5]
+ -par[1]*par[1]*par[5]
+ +par[1]*par[4]*par[5]
+ +par[0]*par[3]*(par[2] + par[5])
+ -par[3]*par[3]*par[2])*lCosPsi*lCosTht
+ -TMath::Sqrt(inSqrt)));
+
+ return zD;
+}
+
+Double_t zpCenter(Double_t *x, Double_t *par){
+ // par : x1l, y1l, z1l, x2l, y2l, z2l, psi, tht
+
+ Double_t lCosPsi = TMath::Cos(par[6]);
+ Double_t lSinPsi = TMath::Sin(par[6]);
+ Double_t lCosTht = TMath::Cos(par[7]);
+ Double_t lSinTht = TMath::Sin(par[7]);
+
+ Double_t inSqrt = ((par[3]*par[1] - par[0]*par[4])*(par[3]*par[1] - par[0]*par[4])
+ *((-(x[0] - x[1])*(x[0] - x[1]))
+ +(((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4])))*lSinPsi*lSinPsi
+ +lCosPsi*((-(par[2] - par[5]))
+ *lCosTht*(-2*x[0]+2*x[1]
+ +(par[2] - par[5])*lCosPsi*lCosTht)
+ +((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4]))*lCosPsi*lSinTht*lSinTht)));
+
+ Double_t zD = ((1./((par[0] - par[3])*(par[0] - par[3])
+ +(par[1] - par[4])*(par[1] - par[4])))
+ *(-par[1]*par[4]*x[0]
+ +par[4]*par[4]*x[0]
+ +par[0]*par[0]*x[1]
+ +par[1]*par[1]*x[1]
+ -par[1]*par[4]*x[1]
+ -par[0]*par[3]*(x[0] + x[1])
+ +par[3]*par[3]*x[0]
+ +(+par[1]*par[4]*par[2]
+ -par[4]*par[4]*par[2]
+ -par[0]*par[0]*par[5]
+ -par[1]*par[1]*par[5]
+ +par[1]*par[4]*par[5]
+ +par[0]*par[3]*(par[2] + par[5])
+ -par[3]*par[3]*par[2])*lCosPsi*lCosTht
+ +TMath::Sqrt(inSqrt)));
+
+ return zD;
+}
+
+//______________________________________________________________________
+AliMUONGeometryTransformer *ReAlign(const AliMUONGeometryTransformer * transformer,
+ int rMod, TGeoCombiTrans deltaDetElemTransf[], Bool_t verbose)
+{
+ /////////////////////////////////////////////////////////////////////
+ // Takes the internal geometry module transformers, copies them
+ // and gets the Detection Elements from them.
+ // Takes misalignment parameters and applies these
+ // to the local transform of the Detection Element
+ // Obtains the global transform by multiplying the module transformer
+ // transformation with the local transformation
+ // Applies the global transform to a new detection element
+ // Adds the new detection element to a new module transformer
+ // Adds the new module transformer to a new geometry transformer
+ // Returns the new geometry transformer
+
+
+ Int_t iDetElemId = 0;
+ Int_t iDetElemNumber = 0;
+ Int_t iDetElemIndex = 0;
+ Int_t iCh = 0;
+
+ AliMUONGeometryTransformer *newGeometryTransformer =
+ new AliMUONGeometryTransformer();
+ for (Int_t iMt = 0; iMt < transformer->GetNofModuleTransformers(); iMt++) {
+ // module transformers
+ const AliMUONGeometryModuleTransformer *kModuleTransformer =
+ transformer->GetModuleTransformer(iMt, true);
+
+ AliMUONGeometryModuleTransformer *newModuleTransformer =
+ new AliMUONGeometryModuleTransformer(iMt);
+ newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
+
+ TGeoCombiTrans moduleTransform =
+ TGeoCombiTrans(*kModuleTransformer->GetTransformation());
+ // New module transformation
+ TGeoCombiTrans *newModuleTransform;
+ if (iMt==rMod) {
+ newModuleTransform = new TGeoCombiTrans(moduleTransform);
+ } else {
+ newModuleTransform = new TGeoCombiTrans(moduleTransform);
+ }
+ newModuleTransformer->SetTransformation(*newModuleTransform);
+
+ // For the selected chamber add misalign module
+ if (iMt==rMod) {
+ // Get delta transformation:
+ // Tdelta = Tnew * Told.inverse
+ TGeoHMatrix deltaModuleTransform =
+ AliMUONGeometryBuilder::Multiply(*newModuleTransform,
+ kModuleTransformer->GetTransformation()->Inverse());
+ // Create module mis alignment matrix
+ newGeometryTransformer
+ ->AddMisAlignModule(kModuleTransformer->GetModuleId(), deltaModuleTransform);
+ }
+
+ AliMpExMap *detElements = kModuleTransformer->GetDetElementStore();
+
+ if (verbose)
+ printf("%i DEs in old GeometryStore %i\n",detElements->GetSize(), iMt);
+ TGeoCombiTrans *deltaLocalTransform;
+ for (Int_t iDe = 0; iDe < detElements->GetSize(); iDe++) {
+ // detection elements.
+ AliMUONGeometryDetElement *detElement =
+ (AliMUONGeometryDetElement *) detElements->GetObject(iDe);
+ if (!detElement) {
+ printf("Detection element not found.\n");
+ break;
+ }
+ /// make a new detection element
+ AliMUONGeometryDetElement *newDetElement =
+ new AliMUONGeometryDetElement(detElement->GetId(),
+ detElement->GetVolumePath());
+ TString lDetElemName(detElement->GetDEName());
+ lDetElemName.ReplaceAll("DE","");
+ iDetElemId = lDetElemName.Atoi();
+ iDetElemNumber = iDetElemId%100;
+ iCh = iDetElemId/100 -1;
+ if(iMt==rMod){
+ if (iCh<4) {
+ iDetElemIndex = iDetElemId;
+ } else {
+ if ((iDetElemNumber > (fgNDetElemCh[iCh]-2)/4) &&
+ (iDetElemNumber < fgNDetElemCh[iCh]-(fgNDetElemCh[iCh]-2)/4)) {
+ iDetElemIndex = (+fgNDetElemCh[iCh]
+ -(1+(fgNDetElemCh[iCh]-2)/4)
+ -iDetElemNumber);
+ } else {
+ iDetElemIndex = (+fgNDetElemCh[iCh]
+ -fgNDetElemCh[iCh]/2
+ -((1+(fgNDetElemCh[iCh]-2)/4)
+ -TMath::Min(iDetElemNumber,
+ TMath::Abs(iDetElemNumber-fgNDetElemCh[iCh]))));
+ }
+ }
+ deltaLocalTransform = new TGeoCombiTrans(deltaDetElemTransf[iDetElemIndex]);
+ } else {
+ deltaLocalTransform = new TGeoCombiTrans(*gGeoIdentity);
+ }
+
+ // local transformation of this detection element.
+ TGeoCombiTrans localTransform
+ = TGeoCombiTrans(*detElement->GetLocalTransformation());
+ // TGeoHMatrix newLocalMatrix = localTransform * (*deltaLocalTransform);
+ TGeoCombiTrans newLocalTransform
+ = TGeoCombiTrans(localTransform * (*deltaLocalTransform));
+ newDetElement->SetLocalTransformation(newLocalTransform);
+ // global transformation
+ TGeoHMatrix newGlobalTransform =
+ AliMUONGeometryBuilder::Multiply(*newModuleTransform,
+ newLocalTransform);
+ newDetElement->SetGlobalTransformation(newGlobalTransform);
+
+ // add this det element to module
+ newModuleTransformer->GetDetElementStore()->Add(newDetElement->GetId(),
+ newDetElement);
+
+ // In the Alice Alignment Framework misalignment objects store
+ // global delta transformation
+ // Get detection "intermediate" global transformation
+ TGeoHMatrix newOldGlobalTransform = (*newModuleTransform) * localTransform;
+ // Get detection element global delta transformation:
+ // Tdelta = Tnew * Told.inverse
+ TGeoHMatrix deltaGlobalTransform
+ = AliMUONGeometryBuilder::Multiply(newGlobalTransform,
+ newOldGlobalTransform.Inverse());
+
+ // Create mis alignment matrix
+ newGeometryTransformer
+ ->AddMisAlignDetElement(detElement->GetId(), deltaGlobalTransform);
+ }
+
+ if (verbose)
+ printf("Added module transformer %i to the transformer\n", iMt);
+ newGeometryTransformer->AddModuleTransformer(newModuleTransformer);
+ }
+ return newGeometryTransformer;
+}
git://git.uio.no / u / mrichter / AliRoot.git / commitdiff