MUON/MUONSurveyCh8L.C

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 /* $Id$ */
  17
  18 // ---
  19 // Macro to process survey and photogrammetry data of chamber 8L
  20 //
  21 // Macro loads the survey data from .txt file using AliSurveyObj.
  22 // Macro MUONSurveyUtil.C is then loaded.
  23 // The transformations of the slats are obatained in 2 steps:
  24 //   1. Fit a plane to the sticker targets -> psi, theta
  25 //   2. Using above psi in theta obtain xc, yc, zc and phi by solving
  26 //      the equations from a local to global transformation of the
  27 //      fixed button targets
  28 // Various histograms are filled and printed for monitoring.
  29 // MisAlignment object is then created.
  30 //
  31 // Author: Javier Castillo
  32 // ---
  33
  34 #if !defined(__CINT__) || defined(__MAKECINT__)
  35
  36 #include "AliMUONGeometryTransformer.h"
  37 #include "AliMUONGeometryMisAligner.h"
  38
  39 #include "AliSurveyObj.h"
  40 #include "AliSurveyPoint.h"
  41 #include "AliGeomManager.h"
  42 #include "AliCDBManager.h"
  43 #include "AliCDBMetaData.h"
  44 #include "AliCDBId.h"
  45
  46 #include <TROOT.h>
  47 #include <TGeoManager.h>
  48 #include <TClonesArray.h>
  49 #include <TObjArray.h>
  50 #include <TObjString.h>
  51 #include <TMath.h>
  52 #include <TString.h>
  53 #include <Riostream.h>
  54 #include <TF2.h>
  55 #include <TH2.h>
  56 #include <TGraph2DErrors.h>
  57 #include <TGraph.h>
  58 #include <TVector3.h>
  59 #include <TCanvas.h>
  60 #include <TPad.h>
  61 #include <TPostScript.h>
  62 #include <TPaveLabel.h>
  63 #include <TStyle.h>
  64
  65 #include <fstream>
  66
  67 #endif
  68
  69 Bool_t MatrixToAngles(const Double_t *rot, Double_t *angles);
  70 Double_t eqPlane(Double_t *x, Double_t *par);
  71 Double_t xpCenter(Double_t *x, Double_t *par);
  72 Double_t xnCenter(Double_t *x, Double_t *par);
  73 Double_t ypCenter(Double_t *x, Double_t *par);
  74 Double_t ynCenter(Double_t *x, Double_t *par);
  75 Double_t zpCenter(Double_t *x, Double_t *par);
  76 Double_t znCenter(Double_t *x, Double_t *par);
  77 Double_t phixpp(Double_t *x, Double_t *par);
  78 Double_t phixpn(Double_t *x, Double_t *par);
  79 Double_t phixnp(Double_t *x, Double_t *par);
  80 Double_t phixnn(Double_t *x, Double_t *par);
  81 Double_t phiypp(Double_t *x, Double_t *par);
  82 Double_t phiypn(Double_t *x, Double_t *par);
  83 Double_t phiynp(Double_t *x, Double_t *par);
  84 Double_t phiynn(Double_t *x, Double_t *par);
  85 AliMUONGeometryTransformer *ReAlign(const AliMUONGeometryTransformer * transformer,
  86                                     int rMod, TGeoCombiTrans deltaDetElemTransf[], Bool_t verbose);
  87 void MUONSurveyCh8L() {
  88
  89   char str[100];
  90   char opt[100];
  91   char var[100];
  92   char filename[100];
  93   Char_t histoName[20];
  94   Char_t histoTitle[50];
  95
  96   int saveps = 1;
  97   const int font = 41; // Helvetica + precision 1
  98   const int cWidth = (int)(700*(29./21.));
  99   const int cHeight = 700;
 100   const int   lineColor = 1;
 101   const int   fillColor = 4;
 102 //   const int   filetype  = 111; // portrait
 103   const int   filetype  = 112; // landscape
 104
 105   sprintf(filename,"surveyChamber8L.ps");
 106
 107   Int_t nSlats = 13;
 108
 109   gROOT->LoadMacro("MUONSurveyUtil.C+");
 110
 111   AliSurveyObj *so = new AliSurveyObj();
 112
 113   Int_t size = so->GetEntries();
 114   printf("-> %d\n", size);
 115
 116   so->FillFromLocalFile("Alice_MuonSystem_Chamber8LCavern_3561b.txt");
 117   size = so->GetEntries();
 118   printf("--> %d\n", size);
 119
 120   Printf("Title: \"%s\"", so->GetReportTitle().Data());
 121   Printf("Date: \"%s\"", so->GetReportDate().Data());
 122   Printf("Detector: \"%s\"", so->GetDetector().Data());
 123   Printf("URL: \"%s\"", so->GetURL().Data());
 124   Printf("Number: \"%d\"", so->GetReportNumber());
 125   Printf("Version: \"%d\"", so->GetReportVersion());
 126   Printf("Observations: \"%s\"", so->GetObservations().Data());
 127   Printf("Coordinate System: \"%s\"", so->GetCoordSys().Data());
 128   Printf("Measurement Units: \"%s\"", so->GetUnits().Data());
 129   Printf("Nr Columns: \"%d\"", so->GetNrColumns());
 130
 131   TObjArray *colNames = so->GetColumnNames();
 132   for (Int_t i = 0; i < colNames->GetEntries(); ++i)
 133     Printf("  Column %d --> \"%s\"", i, ((TObjString *) colNames->At(i))->GetString().Data());
 134
 135   // Get Array of surveyed points
 136   Printf("Points:");
 137   TObjArray *points = so->GetData();
 138
 139   for (Int_t i = 0; i < points->GetEntries(); ++i)
 140     Printf("  Point %d --> \"%s\"  %s ", i, ((AliSurveyPoint *) points->At(i))->GetPointName().Data(), points->At(i)->GetName());
 141
 142
 143   // Slats (#1 - #13) button targets local coordinates
 144   Double_t lSBTLoc6[13][2][3] = {{{  -412.50,  0.0, -(11.75+ 8.20+20.00)},
 145                                   {   412.50,  0.0, -(11.75+16.50+20.00)}},
 146                                  {{  -612.50,  0.0,  (11.75+ 8.20+20.00)},
 147                                   {   612.50,  0.0,  (11.75+16.50+20.00)}},
 148                                  {{ - 812.50,  0.0, -(11.75+ 8.20+20.00)},
 149                                   {   812.50,  0.0, -(11.75+16.50+20.00)}},
 150                                  {{ -1012.50,  0.0,  (11.75+ 8.20+20.00)},
 151                                   {  1012.50,  0.0,  (11.75+16.50+20.00)}},
 152                                  {{ -1012.50,  0.0, -(11.75+ 8.20+20.00)},
 153                                   {  1012.50,  0.0, -(11.75+16.50+20.00)}},
 154                                  {{ -1212.50,  5.0,  (11.75+ 8.20+20.00)},
 155                                   {  1212.50,  0.0,  (11.75+16.50+20.00)}},
 156                                  {{ -1012.50,  0.0, -(11.75+ 8.20+20.00)},
 157                                   {  1012.50,  0.0, -(11.75+16.50+20.00)}},
 158                                  {{ -1212.50,  5.0,  (11.75+ 8.20+20.00)},
 159                                   {  1212.50,  0.0,  (11.75+16.50+20.00)}},
 160                                  {{ -1012.50,  0.0, -(11.75+ 8.20+20.00)},
 161                                   {  1012.50,  0.0, -(11.75+16.50+20.00)}},
 162                                  {{ -1012.50,  0.0,  (11.75+ 8.20+20.00)},
 163                                   {  1012.50,  0.0,  (11.75+16.50+20.00)}},
 164                                  {{  -812.50,  0.0, -(11.75+ 8.20+20.00)},
 165                                   {   812.50,  0.0, -(11.75+16.50+20.00)}},
 166                                  {{  -612.50,  0.0,  (11.75+ 8.20+20.00)},
 167                                   {   612.50,  0.0,  (11.75+16.50+20.00)}},
 168                                  {{  -412.50,  0.0, -(11.75+ 8.20+20.00)},
 169                                   {   412.50,  0.0, -(11.75+16.50+20.00)}}};
 170
 171
 172   AliSurveyPoint *pointSST = 0;
 173   AliSurveyPoint *pointCST = 0;
 174   AliSurveyPoint *pointCPST = 0;
 175   AliSurveyPoint **pointSBT = new AliSurveyPoint*[2];
 176
 177   char sPointName[10] = "5000";
 178
 179   // Print length of slats
 180   cout << "Slat lengths:" << endl;
 181   TVector3 vTemp(0., 0., 0.);
 182   TVector3 vSBT(0., 0., 0.);
 183
 184   for (Int_t iSlat=0; iSlat<nSlats; iSlat++){
 185     // Get button targets survey points
 186     vTemp.SetXYZ(0., 0., 0.);
 187     if (iSlat+1<10) {
 188       sprintf(sPointName,"60%d%d",iSlat+1,1);
 189       pointSBT[0] = (AliSurveyPoint *)points->FindObject(sPointName);
 190       if(!pointSBT[0]) {
 191         cout << "Error! No button targets ... " << endl;
 192         break;
 193       }
 194       vSBT.SetXYZ(pointSBT[0]->GetX(),pointSBT[0]->GetY(),pointSBT[0]->GetZ());
 195       vTemp+=vSBT;
 196       sprintf(sPointName,"60%d%d",iSlat+1,2);
 197       pointSBT[1] = (AliSurveyPoint *)points->FindObject(sPointName);
 198       if(!pointSBT[1]) {
 199         cout << "Error! No button targets ... " << endl;
 200         break;
 201       }
 202       vSBT.SetXYZ(pointSBT[1]->GetX(),pointSBT[1]->GetY(),pointSBT[1]->GetZ());
 203       vTemp-=vSBT;
 204     }
 205     else {
 206       sprintf(sPointName,"6%d%d",iSlat+1,1);
 207       pointSBT[0] = (AliSurveyPoint *)points->FindObject(sPointName);
 208       if(!pointSBT[0]) {
 209         cout << "Error! No button targets ... " << endl;
 210         break;
 211       }
 212       vSBT.SetXYZ(pointSBT[0]->GetX(),pointSBT[0]->GetY(),pointSBT[0]->GetZ());
 213       vTemp+=vSBT;
 214       sprintf(sPointName,"6%d%d",iSlat+1,2);
 215       pointSBT[1] = (AliSurveyPoint *)points->FindObject(sPointName);
 216       if(!pointSBT[1]) {
 217         cout << "Error! No button targets ... " << endl;
 218         break;
 219       }
 220       vSBT.SetXYZ(pointSBT[1]->GetX(),pointSBT[1]->GetY(),pointSBT[1]->GetZ());
 221       vTemp-=vSBT;
 222     }
 223     cout << "Slat " << iSlat+1 << ": " << vTemp.Mag() << endl;
 224   }
 225
 226   // Histograms for monitoring
 227   TH2F *hCPSTry = new TH2F("hCPSTry","hCPSTry",28,-600,2200,52,0,5200);
 228   TH2F *hCPSTly = new TH2F("hCPSTly","hCPSTly",28,-600,2200,52,0,5200);
 229   TH2F *hSSTry = new TH2F("hSSTry","hSSTry",70,-600,2200,130,0,5200);
 230   TH2F *hSSTly = new TH2F("hSSTly","hSSTly",70,-600,2200,130,0,5200);
 231   TH2F *hCSTy = new TH2F("hCSTy","hCSTy",70,-600,2200,130,0,5200);
 232   TH2F *hSSTrpy = new TH2F("hSSTrpy","hSSTrpy",70,-600,2200,130,0,5200);
 233   TH2F *hSSTlpy = new TH2F("hSSTlpy","hSSTlpy",70,-600,2200,130,0,5200);
 234
 235   // Chamber plane sticker targets
 236   for (int iPoint=0; iPoint<9; iPoint++) {
 237     sprintf(sPointName,"700%d",iPoint+1);
 238     pointCPST = (AliSurveyPoint *)points->FindObject(sPointName);
 239     if(!pointCPST) {
 240       printf("Point %s is missing ...\n",sPointName);
 241       break;
 242     }
 243     hCPSTry->Fill(pointCPST->GetX(),pointCPST->GetZ(),pointCPST->GetY());
 244   }
 245   for (int iPoint=9; iPoint<18; iPoint++) {
 246     sprintf(sPointName,"70%d",iPoint+1);
 247     pointCPST = (AliSurveyPoint *)points->FindObject(sPointName);
 248     if(!pointCPST) {
 249       printf("Point %s is missing ...\n",sPointName);
 250       break;
 251     }
 252     hCPSTly->Fill(pointCPST->GetX(),pointCPST->GetZ(),pointCPST->GetY());
 253   }
 254
 255   // Chamber Side Targets
 256   for (int iPoint=0; iPoint<25; iPoint++) {
 257     if (iPoint+1<10) {
 258       sprintf(sPointName,"800%d",iPoint+1);
 259     }
 260     else {
 261       sprintf(sPointName,"80%d",iPoint+1);
 262     }
 263     pointCST = (AliSurveyPoint *)points->FindObject(sPointName);
 264     if(!pointCPST) {
 265       printf("Point %s is missing ...\n",sPointName);
 266       break;
 267     }
 268     hCSTy->Fill(pointCST->GetX(),pointCST->GetZ(),pointCST->GetY());
 269   }
 270
 271
 272   // Graphs used for plane fitting
 273   TGraph2DErrors **gSST5 = new TGraph2DErrors*[nSlats];
 274   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 275     gSST5[iSlat] = new TGraph2DErrors();
 276   }
 277
 278   // Keep the id of slat sticker targets next to slat button targets
 279   Int_t iSSBT[13][2] = {0};
 280
 281   // Fill graph with sticker target positions
 282   for (int iSlat=0; iSlat<nSlats; iSlat++){
 283     for (int iPoint=0; iPoint<9; iPoint++) {
 284       // Second sticker target is next to first button target
 285       // Previous to last sticker target is next to second button target
 286       iSSBT[iSlat][0] = 2;
 287       iSSBT[iSlat][1] = 8;
 288       if (iSlat+1<10) {
 289         sprintf(sPointName,"50%d%d",iSlat+1,iPoint+1);
 290       }
 291       else {
 292         sprintf(sPointName,"5%d%d",iSlat+1,iPoint+1);
 293       }
 294       pointSST = (AliSurveyPoint *)points->FindObject(sPointName);
 295       if(!pointSST) {
 296         printf("%s\n",sPointName);
 297         cout << iSlat << " " << iPoint  << " " << pointSST << endl;
 298         // Previous to last sticker target is next to second button target
 299         iSSBT[iSlat][1] = iPoint+1-2;
 300         break;
 301       }
 302
 303       gSST5[iSlat]->SetPoint(iPoint,-pointSST->GetX(),pointSST->GetZ(),pointSST->GetY());
 304       gSST5[iSlat]->SetPointError(iPoint,pointSST->GetPrecisionX(),pointSST->GetPrecisionZ(),pointSST->GetPrecisionY());
 305
 306       // Fill histograms of sticker targets. For monitoring purposes.
 307       if((iSlat+1)%2==0){
 308         hSSTly->Fill(pointSST->GetX(),pointSST->GetZ(),pointSST->GetY());
 309       }
 310       else {
 311         hSSTry->Fill(pointSST->GetX(),pointSST->GetZ(),pointSST->GetY());
 312       }
 313     }
 314   }
 315
 316   Float_t xMin = -2200.;
 317   Float_t xMax =  2200.;
 318   Float_t yMin = -5200.;
 319   Float_t yMax =  5200.;
 320   Float_t zMin = -200.;
 321   Float_t zMax =  200.;
 322
 323   Double_t xMinSlat, xMaxSlat;
 324   Double_t yMinSlat, yMaxSlat;
 325
 326   // Slat plane function
 327   char fsName[100] = "fSlat00";
 328   TF2 **fSlat = new TF2*[nSlats];
 329   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 330     sprintf(fsName,"fSlat%d",iSlat+1);
 331     fSlat[iSlat] = new TF2(fsName,eqPlane,xMin,xMax,yMin,yMax,3);
 332   }
 333
 334   // Xcenter functions
 335   char fxcName[100] = "fXcnSlat00";
 336   TF2 ***fXcSlat = new TF2**[nSlats];
 337   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 338     fXcSlat[iSlat] = new TF2*[2];
 339     sprintf(fxcName,"fXcnSlat%d",iSlat+1);
 340     fXcSlat[iSlat][0] = new TF2(fxcName,xnCenter,xMin,xMax,yMin,yMax,7);
 341     sprintf(fxcName,"fXcpSlat%d",iSlat+1);
 342     fXcSlat[iSlat][1] = new TF2(fxcName,xpCenter,xMin,xMax,yMin,yMax,7);
 343   }
 344
 345   // Ycenter functions
 346   char fycName[100] = "fYcnSlat00";
 347   TF2 ***fYcSlat = new TF2**[nSlats];
 348   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 349     fYcSlat[iSlat] = new TF2*[2];
 350     sprintf(fycName,"fYcnSlat%d",iSlat+1);
 351     fYcSlat[iSlat][0] = new TF2(fycName,ynCenter,yMin,yMax,yMin,yMax,8);
 352     sprintf(fycName,"fYcpSlat%d",iSlat+1);
 353     fYcSlat[iSlat][1] = new TF2(fycName,ypCenter,yMin,yMax,yMin,yMax,8);
 354   }
 355
 356   // Zcenter functions
 357   char fzcName[100] = "fZcnSlat00";
 358   TF2 ***fZcSlat = new TF2**[nSlats];
 359   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 360     fZcSlat[iSlat] = new TF2*[2];
 361     sprintf(fzcName,"fZcnSlat%d",iSlat+1);
 362     fZcSlat[iSlat][0] = new TF2(fzcName,znCenter,zMin,zMax,zMin,zMax,8);
 363     sprintf(fzcName,"fZcpSlat%d",iSlat+1);
 364     fZcSlat[iSlat][1] = new TF2(fzcName,zpCenter,zMin,zMax,zMin,zMax,8);
 365   }
 366
 367   // Phi rotation using xglobal coords functions
 368   char fphixName[100] = "fPhiXnnSlat00";
 369   TF2 ****fPhiXSlat = new TF2***[nSlats];
 370   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 371     fPhiXSlat[iSlat] = new TF2**[2];
 372     for (Int_t iX =0; iX<2; iX++)
 373       fPhiXSlat[iSlat][iX] = new TF2*[2];
 374     sprintf(fphixName,"fPhiXnnSlat%d",iSlat+1);
 375     fPhiXSlat[iSlat][0][0] = new TF2(fphixName,phixnn,xMin,xMax,xMin,xMax,7);
 376     sprintf(fphixName,"fPhixnpSlat%d",iSlat+1);
 377     fPhiXSlat[iSlat][0][1] = new TF2(fphixName,phixnp,xMin,xMax,xMin,xMax,7);
 378     sprintf(fphixName,"fPhiXpnSlat%d",iSlat+1);
 379     fPhiXSlat[iSlat][1][0] = new TF2(fphixName,phixpn,xMin,xMax,xMin,xMax,7);
 380     sprintf(fphixName,"fPhixppSlat%d",iSlat+1);
 381     fPhiXSlat[iSlat][1][1] = new TF2(fphixName,phixpp,xMin,xMax,xMin,xMax,7);
 382   }
 383
 384   // Phi rotation using yglobal coords functions
 385   char fphiyName[100] = "fPhiYnnSlat00";
 386   TF2 ****fPhiYSlat = new TF2***[nSlats];
 387   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 388     fPhiYSlat[iSlat] = new TF2**[2];
 389     for (Int_t iY =0; iY<2; iY++)
 390       fPhiYSlat[iSlat][iY] = new TF2*[2];
 391     sprintf(fphiyName,"fPhiYnnSlat%d",iSlat+1);
 392     fPhiYSlat[iSlat][0][0] = new TF2(fphiyName,phiynn,yMin,yMax,yMin,yMax,7);
 393     sprintf(fphiyName,"fPhiYnpSlat%d",iSlat+1);
 394     fPhiYSlat[iSlat][0][1] = new TF2(fphiyName,phiynp,yMin,yMax,yMin,yMax,7);
 395     sprintf(fphiyName,"fPhiYpnSlat%d",iSlat+1);
 396     fPhiYSlat[iSlat][1][0] = new TF2(fphiyName,phiypn,yMin,yMax,yMin,yMax,7);
 397     sprintf(fphiyName,"fPhiYppSlat%d",iSlat+1);
 398     fPhiYSlat[iSlat][1][1] = new TF2(fphiyName,phiypp,yMin,yMax,yMin,yMax,7);
 399   }
 400
 401   Double_t *xce = new Double_t[nSlats];
 402   Double_t *yce = new Double_t[nSlats];
 403   Double_t *zce = new Double_t[nSlats];
 404   Double_t *psi = new Double_t[nSlats];
 405   Double_t *tht = new Double_t[nSlats];
 406   Double_t *phi = new Double_t[nSlats];
 407
 408   Double_t **lCenSlat = new Double_t*[nSlats];
 409   Double_t **lRotSlat = new Double_t*[nSlats];
 410   Double_t **lDiffCenSlat0 = new Double_t*[nSlats];
 411   Double_t **lDiffRotSlat0 = new Double_t*[nSlats];
 412   Double_t **lDiffThCenSlat0 = new Double_t*[nSlats];
 413   Double_t **lDiffThRotSlat0 = new Double_t*[nSlats];
 414   Double_t **lDeltaDiffCenSlat0 = new Double_t*[nSlats];
 415   Double_t **lDeltaDiffRotSlat0 = new Double_t*[nSlats];
 416
 417   for (int iSlat=0; iSlat<nSlats; iSlat++) {
 418     lCenSlat[iSlat] = new Double_t[3];
 419     lRotSlat[iSlat] = new Double_t[3];
 420
 421     lDiffCenSlat0[iSlat] = new Double_t[3];
 422     lDiffRotSlat0[iSlat] = new Double_t[3];
 423     lDiffThCenSlat0[iSlat] = new Double_t[3];
 424     lDiffThRotSlat0[iSlat] = new Double_t[3];
 425     lDeltaDiffCenSlat0[iSlat] = new Double_t[3];
 426     lDeltaDiffRotSlat0[iSlat] = new Double_t[3];
 427   }
 428
 429
 430   TGeoTranslation transSlat[nSlats];
 431   TGeoRotation rotSlat[nSlats];
 432   TGeoCombiTrans trfSlat[nSlats];
 433
 434   TGeoTranslation dtransSlat[nSlats];
 435   TGeoRotation drotSlat[nSlats];
 436   TGeoCombiTrans dtrfSlat[nSlats];
 437
 438   TGeoTranslation transTemp;
 439   TGeoRotation rotTemp;
 440   TGeoCombiTrans trfTemp;
 441
 442
 443   Double_t lCenTemp[3];
 444   Double_t lRotTemp[3];
 445
 446   Double_t lDiffTemp[9];
 447   Double_t lDiffMin[9];
 448
 449   double tempDiff = 0.;
 450   double tempDiff1 = 0.;
 451   double tempDiff2 = 0.;
 452
 453   AliSurveyPoint **pointSSBT = new AliSurveyPoint*[2];
 454
 455   //
 456   // Get Slat transformation.
 457   // Psi and Theta are obtained by fitting a plane to the sticker targets.
 458   // Then Xc, Yc, Zc and Phi are obtained by solving the equations to the ref.
 459   // syst. transformation of the button targets
 460   //
 461   for (int iSlat=0; iSlat<nSlats; iSlat++) {
 462     sprintf(fsName,"fSlat%d",iSlat+1);
 463     cout << "Fitting Slat" << iSlat+1 << " ..." << endl;
 464
 465     // Fit a plane to the sticker targets
 466     gSST5[iSlat]->Fit(fsName,"","same");
 467
 468     psi[iSlat] = TMath::ATan(fSlat[iSlat]->GetParameter(1));
 469     tht[iSlat] = TMath::ATan(fSlat[iSlat]->GetParameter(0));
 470     if (iSlat==5)
 471       psi[iSlat] += TMath::Pi(); // Rotated slat
 472
 473     lRotSlat[iSlat][0] = psi[iSlat];
 474     lRotSlat[iSlat][1] = tht[iSlat];
 475
 476     for(Int_t iS=0; iS<2; iS++){
 477       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]);
 478       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]);
 479       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]);
 480       for(Int_t jS=0; jS<2; jS++){
 481         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]);
 482         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]);
 483       }
 484     }
 485
 486     //
 487     // Calculate Slat Center from button targets
 488     //
 489
 490     // Get button targets survey points
 491     for (Int_t iPoint=0; iPoint<2; iPoint++) {
 492       if (iSlat+1<10) {
 493         sprintf(sPointName,"60%d%d",iSlat+1,iPoint+1);
 494         pointSBT[iPoint] = (AliSurveyPoint *)points->FindObject(sPointName);
 495         if(!pointSBT[iPoint]) {
 496           cout << "Error! No button targets ... " << endl;
 497           break;
 498         }
 499         sprintf(sPointName,"50%d%d",iSlat+1,iSSBT[iSlat][iPoint]);
 500         pointSSBT[iPoint] = (AliSurveyPoint *)points->FindObject(sPointName);
 501         if(!pointSSBT[iPoint]) {
 502           cout << "Error! No sticker target ... " << sPointName << endl;
 503           break;
 504         }
 505       }
 506       else {
 507         sprintf(sPointName,"6%d%d",iSlat+1,iPoint+1);
 508         pointSBT[iPoint] = (AliSurveyPoint *)points->FindObject(sPointName);
 509         if(!pointSBT[iPoint]) {
 510           cout << "Error! No button targets ... " << endl;
 511           break;
 512         }
 513         sprintf(sPointName,"5%d%d",iSlat+1,iSSBT[iSlat][iPoint]);
 514         pointSSBT[iPoint] = (AliSurveyPoint *)points->FindObject(sPointName);
 515         if(!pointSSBT[iPoint]) {
 516           cout << "Error! No sticker targets ... " << sPointName << endl;
 517           break;
 518         }
 519       }
 520     }
 521
 522     tempDiff += TMath::Power(-1,iSlat)*((pointSBT[1]->GetY() - pointSSBT[1]->GetY())-(pointSBT[0]->GetY() - pointSSBT[0]->GetY()));
 523     tempDiff1 += TMath::Abs(pointSBT[0]->GetY() - pointSSBT[0]->GetY())-20;
 524     tempDiff2 += TMath::Abs(pointSBT[1]->GetY() - pointSSBT[1]->GetY())-20;
 525     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;
 526
 527
 528     Double_t p0l[3] = {0};
 529     Double_t p1l[3] = {0};
 530     Double_t p2l[3] = {0};
 531     Double_t p0g[3] = {0};
 532     Double_t p1g[3] = {0};
 533     Double_t p2g[3] = {0};
 534
 535     p0l[2] = lSBTLoc6[iSlat][0][2];
 536     // Button targets local coordinates
 537     for(Int_t iCor=0; iCor<3; iCor++){
 538       p1l[iCor]= lSBTLoc6[iSlat][0][iCor];
 539       p2l[iCor]= lSBTLoc6[iSlat][1][iCor];
 540     }
 541     for(Int_t i=0; i<9; i++){
 542       lDiffMin[i]=1000000.;
 543     }
 544     // Trying 2x*2y*2z*2phi possibilities
 545     for(Int_t iX=0; iX<2; iX++){
 546       for(Int_t iY=0; iY<2; iY++){
 547         for(Int_t iZ=0; iZ<2; iZ++){
 548           lCenTemp[0] = fXcSlat[iSlat][iX]->Eval(-pointSBT[0]->GetX(),-pointSBT[1]->GetX());
 549           lCenTemp[1] = fYcSlat[iSlat][iY]->Eval(pointSBT[0]->GetZ(),pointSBT[1]->GetZ());
 550           lCenTemp[2] = fZcSlat[iSlat][iZ]->Eval(pointSBT[0]->GetY(),pointSBT[1]->GetY());
 551           //     lCenTemp[2] = fZcSlat[iSlat][iZ]->Eval(pointSSBT[0]->GetY(),pointSSBT[1]->GetY());
 552           lRotTemp[0] = psi[iSlat];
 553           lRotTemp[1] = tht[iSlat];
 554           for(Int_t iP=0; iP<2; iP++){
 555             lRotTemp[2] = fPhiXSlat[iSlat][iX][iP]->Eval(-pointSBT[0]->GetX(),-pointSBT[1]->GetX());
 556
 557             trfTemp.SetTranslation(transTemp);
 558             trfTemp.SetRotation(rotTemp);
 559             trfTemp.Clear();
 560             trfTemp.RotateZ(TMath::RadToDeg()*lRotTemp[2]);
 561             trfTemp.RotateY(TMath::RadToDeg()*lRotTemp[1]);
 562             trfTemp.RotateX(TMath::RadToDeg()*lRotTemp[0]);
 563             trfTemp.SetTranslation(lCenTemp);
 564
 565             //     trfTemp.LocalToMaster(p0l, p0g);
 566             //     lCenTemp[2]= fSlat[iSlat]->Eval(p0g[0],p0g[1]) - trfTemp.GetRotationMatrix()[8]*p0l[2];
 567             //     trfTemp.SetTranslation(lCenTemp);
 568
 569             trfTemp.LocalToMaster(p1l, p1g);
 570             trfTemp.LocalToMaster(p2l, p2g);
 571
 572             lDiffTemp[0] = (-pointSBT[0]->GetX()-p1g[0]);
 573             lDiffTemp[1] = (pointSBT[0]->GetZ()-p1g[1]);
 574             lDiffTemp[2] = (pointSBT[0]->GetY()-p1g[2]);
 575             //     lDiffTemp[2] = (pointSSBT[0]->GetY()-p1g[2]);
 576             lDiffTemp[3] = TMath::Sqrt(lDiffTemp[0]*lDiffTemp[0]+lDiffTemp[1]*lDiffTemp[1]+lDiffTemp[2]*lDiffTemp[2]);
 577             lDiffTemp[4] = (-pointSBT[1]->GetX()-p2g[0]);
 578             lDiffTemp[5] = (pointSBT[1]->GetZ()-p2g[1]);
 579             lDiffTemp[6] = (pointSBT[1]->GetY()-p2g[2]);
 580             //     lDiffTemp[6] = (pointSSBT[1]->GetY()-p2g[2]);
 581             lDiffTemp[7] = TMath::Sqrt(lDiffTemp[4]*lDiffTemp[4]+lDiffTemp[5]*lDiffTemp[5]+lDiffTemp[6]*lDiffTemp[6]);
 582             lDiffTemp[8] = TMath::Sqrt(lDiffTemp[3]*lDiffTemp[3]+lDiffTemp[7]*lDiffTemp[7]);
 583
 584             if(lDiffTemp[8]<lDiffMin[8]){
 585               cout << "Diffs" ;
 586               for(Int_t i=0; i<9; i++){
 587                 cout << " " << lDiffTemp[i];
 588               }
 589               cout << endl;
 590               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;
 591               cout << p1g[0] << " " << p1g[1] << " " << p1g[2] << " " << p2g[0] << " " << p2g[1] << " " << p2g[2] << endl;
 592               cout << "Transformation improved ..." << endl;
 593               lCenSlat[iSlat][0] = lCenTemp[0]; lCenSlat[iSlat][1] = lCenTemp[1]; lCenSlat[iSlat][2] = lCenTemp[2];
 594               lRotSlat[iSlat][2] = lRotTemp[2];
 595               for(Int_t i=0; i<9; i++){
 596                 lDiffMin[i]=lDiffTemp[i];
 597               }
 598               if((lDiffMin[3]*lDiffMin[3]<0.1*0.1+0.1*0.1+0.1*0.1)&&
 599                  (lDiffMin[7]*lDiffMin[7]<0.1*0.1+0.1*0.1+0.1*0.1)){
 600                 cout << "Correct Transformation found X " << iX  << " Y " << iY << " Z " << iZ << " xP " << iP << endl;
 601                 lCenSlat[iSlat][0] = lCenTemp[0]; lCenSlat[iSlat][1] = lCenTemp[1]; lCenSlat[iSlat][2] = lCenTemp[2];
 602                 lRotSlat[iSlat][2] = lRotTemp[2];
 603               }
 604             }
 605           }
 606           for(Int_t iP=0; iP<2; iP++){
 607             lRotTemp[2] = fPhiYSlat[iSlat][iY][iP]->Eval(pointSBT[0]->GetZ(),pointSBT[1]->GetZ());
 608             Double_t lPhi = TMath::ATan2((pointSBT[1]->GetZ()-pointSBT[0]->GetZ()),-(pointSBT[1]->GetX()-pointSBT[0]->GetX()));
 609
 610             trfTemp.Clear();
 611             trfTemp.RotateZ(TMath::RadToDeg()*lRotTemp[2]);
 612             trfTemp.RotateY(TMath::RadToDeg()*lRotTemp[1]);
 613             trfTemp.RotateX(TMath::RadToDeg()*lRotTemp[0]);
 614             trfTemp.SetTranslation(lCenTemp);
 615
 616             //     trfTemp.LocalToMaster(p0l, p0g);
 617             //     lCenTemp[2]= fSlat[iSlat]->Eval(p0g[0],p0g[1]) - trfTemp.GetRotationMatrix()[8]*p0l[2];
 618             //     trfTemp.SetTranslation(lCenTemp);
 619
 620             trfTemp.LocalToMaster(p1l, p1g);
 621             trfTemp.LocalToMaster(p2l, p2g);
 622
 623             lDiffTemp[0] = (-pointSBT[0]->GetX()-p1g[0]);
 624             lDiffTemp[1] = (pointSBT[0]->GetZ()-p1g[1]);
 625             lDiffTemp[2] = (pointSBT[0]->GetY()-p1g[2]);
 626             //     lDiffTemp[2] = (pointSSBT[0]->GetY()-p1g[2]);
 627             lDiffTemp[3] = TMath::Sqrt(lDiffTemp[0]*lDiffTemp[0]+lDiffTemp[1]*lDiffTemp[1]+lDiffTemp[2]*lDiffTemp[2]);
 628             lDiffTemp[4] = (-pointSBT[1]->GetX()-p2g[0]);
 629             lDiffTemp[5] = (pointSBT[1]->GetZ()-p2g[1]);
 630             lDiffTemp[6] = (pointSBT[1]->GetY()-p2g[2]);
 631             //     lDiffTemp[6] = (pointSSBT[1]->GetY()-p2g[2]);
 632             lDiffTemp[7] = TMath::Sqrt(lDiffTemp[4]*lDiffTemp[4]+lDiffTemp[5]*lDiffTemp[5]+lDiffTemp[6]*lDiffTemp[6]);
 633             lDiffTemp[8] = TMath::Sqrt(lDiffTemp[3]*lDiffTemp[3]+lDiffTemp[7]*lDiffTemp[7]);
 634
 635             if(lDiffTemp[8]<lDiffMin[8]){
 636               cout << "Diffs" ;
 637               for(Int_t i=0; i<9; i++){
 638                 cout << " " << lDiffTemp[i];
 639               }
 640               cout << endl;
 641               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;
 642               cout << p1g[0] << " " << p1g[1] << " " << p1g[2] << " " << p2g[0] << " " << p2g[1] << " " << p2g[2] << endl;
 643               cout << "Transformation improved ..." << endl;
 644               lCenSlat[iSlat][0] = lCenTemp[0]; lCenSlat[iSlat][1] = lCenTemp[1]; lCenSlat[iSlat][2] = lCenTemp[2];
 645               lRotSlat[iSlat][2] = lRotTemp[2];
 646
 647               for(Int_t i=0; i<9; i++){
 648                 lDiffMin[i]=lDiffTemp[i];
 649               }
 650               if((lDiffMin[3]*lDiffMin[3]<0.1*0.1+0.1*0.1+0.1*0.1)&&
 651                  (lDiffMin[7]*lDiffMin[7]<0.1*0.1+0.1*0.1+0.1*0.1)){
 652                 cout << "Correct Transformation found X " << iX  << " Y " << iY << " Z " << iZ << " yP " << iP << endl;
 653                 lCenSlat[iSlat][0] = lCenTemp[0]; lCenSlat[iSlat][1] = lCenTemp[1]; lCenSlat[iSlat][2] = lCenTemp[2];
 654                 lRotSlat[iSlat][2] = lRotTemp[2];
 655               }
 656             }
 657           }
 658         }
 659       }
 660     }
 661
 662     // Fill slat plane for fit monitor.
 663     xMinSlat = hSSTrpy->GetXaxis()->FindBin(pointSBT[0]->GetX());
 664     xMaxSlat = hSSTrpy->GetXaxis()->FindBin(pointSBT[1]->GetX());
 665     yMinSlat = hSSTrpy->GetYaxis()->FindBin(pointSBT[0]->GetZ()-200.);
 666     yMaxSlat = hSSTrpy->GetYaxis()->FindBin(pointSBT[0]->GetZ()+200.);
 667
 668     for (int i=(int)xMinSlat; i<=(int)xMaxSlat; i++) {
 669       for (int j=(int)yMinSlat; j<=(int)yMaxSlat; j++) {
 670         Double_t zSlat = fSlat[iSlat]->Eval(-hSSTrpy->GetXaxis()->GetBinCenter(i),
 671                                             hSSTrpy->GetYaxis()->GetBinCenter(j));
 672         if((iSlat+1)%2==0){
 673           hSSTlpy->SetBinContent(i,j,zSlat);
 674         }
 675         else {
 676           hSSTrpy->SetBinContent(i,j,zSlat);
 677         }
 678       }
 679     }
 680   }
 681
 682   //
 683   // Compare transformations to expected ones
 684   //
 685   Int_t iSlatToPos[13] = {0, 11, 2, 9, 4, 7, 6, 5, 8, 3, 10, 1, 12};
 686
 687   // Theoretical differences with respect to Slat 513 which is here Slat07
 688   lDiffThCenSlat0[0][0]  = (2-5)*400./2. +12.5 -375. -25.;
 689   lDiffThCenSlat0[1][0]  = (3-5)*400./2. +12.5 -375. -25.;
 690   lDiffThCenSlat0[2][0]  = (4-5)*400./2. +12.5 -375. -25.;
 691   lDiffThCenSlat0[3][0]  = (5-5)*400./2. +12.5 -375. -25.;
 692   lDiffThCenSlat0[4][0]  = (5-5)*400./2. +12.5 -375. -25.;
 693   lDiffThCenSlat0[5][0]  = (6-5)*400./2. +12.5 -375. -25.;
 694   lDiffThCenSlat0[6][0]  = (5-5)*400./2. +12.5 -375. -25. +375. +25. -12.5;
 695   lDiffThCenSlat0[7][0]  = (6-5)*400./2. +12.5 -375. -25.;
 696   lDiffThCenSlat0[8][0]  = (5-5)*400./2. +12.5 -375. -25.;
 697   lDiffThCenSlat0[9][0]  = (5-5)*400./2. +12.5 -375. -25.;
 698   lDiffThCenSlat0[10][0] = (4-5)*400./2. +12.5 -375. -25.;
 699   lDiffThCenSlat0[11][0] = (3-5)*400./2. +12.5 -375. -25.;
 700   lDiffThCenSlat0[12][0] = (2-5)*400./2. +12.5 -375. -25.;
 701
 702   lDiffThCenSlat0[12][1] = 382.0 +378.5 +375.5 +294.0 +370.0 +286.0;
 703   lDiffThCenSlat0[1][1]  = 382.0 +378.5 +375.5 +294.0 +370.0;
 704   lDiffThCenSlat0[10][1] = 382.0 +378.5 +375.5 +294.0;
 705   lDiffThCenSlat0[3][1]  = 382.0 +378.5 +375.5;
 706   lDiffThCenSlat0[8][1]  = 382.0 +378.5;
 707   lDiffThCenSlat0[5][1]  = 382.0;
 708   lDiffThCenSlat0[6][1]  = 0.0;
 709   lDiffThCenSlat0[7][1]  = -382.0;
 710   lDiffThCenSlat0[4][1]  = -382.0 -378.5;
 711   lDiffThCenSlat0[9][1]  = -382.0 -378.5 -375.5;
 712   lDiffThCenSlat0[2][1]  = -382.0 -378.5 -375.5 -294.0;
 713   lDiffThCenSlat0[11][1] = -382.0 -378.5 -375.5 -294.0 -370.0;
 714   lDiffThCenSlat0[0][1]  = -382.0 -378.5 -375.5 -294.0 -370.0 -286.0;
 715
 716   lDiffThCenSlat0[12][2] =  (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
 717   lDiffThCenSlat0[1][2]  = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
 718   lDiffThCenSlat0[10][2] =  (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
 719   lDiffThCenSlat0[3][2]  = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
 720   lDiffThCenSlat0[8][2]  =  (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
 721   lDiffThCenSlat0[5][2]  = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
 722   lDiffThCenSlat0[6][2]  =  (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
 723   lDiffThCenSlat0[7][2]  = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
 724   lDiffThCenSlat0[4][2]  =  (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
 725   lDiffThCenSlat0[9][2]  = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
 726   lDiffThCenSlat0[2][2]  =  (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
 727   lDiffThCenSlat0[11][2] = -(42.5)-(42.5); //-(42.5-1.175)-(42.5-1.175);
 728   lDiffThCenSlat0[0][2]  =  (42.5)-(42.5); // (42.5-1.175)-(42.5-1.175);
 729
 730   TGraph *gDeltaDiffCenXSlat0 = new TGraph(nSlats);
 731   TGraph *gDeltaDiffCenYSlat0 = new TGraph(nSlats);
 732   TGraph *gDeltaDiffCenZSlat0 = new TGraph(nSlats);
 733   TGraph *gDeltaDiffPsiSlat0 = new TGraph(nSlats);
 734   TGraph *gDeltaDiffThtSlat0 = new TGraph(nSlats);
 735   TGraph *gDeltaDiffPhiSlat0 = new TGraph(nSlats);
 736
 737   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 738     trfSlat[iSlat].SetTranslation(transSlat[iSlat]);
 739     trfSlat[iSlat].SetRotation(rotSlat[iSlat]);
 740     trfSlat[iSlat].Clear();
 741     trfSlat[iSlat].RotateZ(TMath::RadToDeg()*lRotSlat[iSlat][2]);
 742     trfSlat[iSlat].RotateY(TMath::RadToDeg()*lRotSlat[iSlat][1]);
 743     trfSlat[iSlat].RotateX(TMath::RadToDeg()*lRotSlat[iSlat][0]);
 744     trfSlat[iSlat].SetTranslation(lCenSlat[iSlat]);
 745   }
 746
 747   Double_t *myTrans = 0;
 748   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 749     dtrfSlat[iSlat].SetTranslation(dtransSlat[iSlat]);
 750     dtrfSlat[iSlat].SetRotation(drotSlat[iSlat]);
 751     dtrfSlat[iSlat].Clear();
 752     dtrfSlat[iSlat] = trfSlat[6].Inverse()*trfSlat[iSlat];
 753     dtrfSlat[iSlat].Print();
 754     lDiffCenSlat0[iSlat] = (Double_t*)dtrfSlat[iSlat].GetTranslation();
 755     MatrixToAngles(dtrfSlat[iSlat].GetRotationMatrix(),lDiffRotSlat0[iSlat]);
 756
 757     lDeltaDiffCenSlat0[iSlat][0] = lDiffCenSlat0[iSlat][0]-lDiffThCenSlat0[iSlat][0];
 758     lDeltaDiffCenSlat0[iSlat][1] =  -lDiffCenSlat0[iSlat][1]-lDiffThCenSlat0[iSlat][1];
 759     lDeltaDiffCenSlat0[iSlat][2] =  -lDiffCenSlat0[iSlat][2]-lDiffThCenSlat0[iSlat][2];
 760     lDeltaDiffRotSlat0[iSlat][0] = lDiffRotSlat0[iSlat][0];
 761     lDeltaDiffRotSlat0[iSlat][1] = lDiffRotSlat0[iSlat][1];
 762     lDeltaDiffRotSlat0[iSlat][2] = lDiffRotSlat0[iSlat][2];
 763     gDeltaDiffCenXSlat0->SetPoint(iSlat,lDeltaDiffCenSlat0[iSlat][0],iSlatToPos[iSlat]+1);
 764     gDeltaDiffCenYSlat0->SetPoint(iSlat,lDeltaDiffCenSlat0[iSlat][1],iSlatToPos[iSlat]+1);
 765     gDeltaDiffCenZSlat0->SetPoint(iSlat,lDeltaDiffCenSlat0[iSlat][2],iSlatToPos[iSlat]+1);
 766     gDeltaDiffPsiSlat0->SetPoint(iSlat,1e3*lDeltaDiffRotSlat0[iSlat][0],iSlatToPos[iSlat]+1);
 767     gDeltaDiffThtSlat0->SetPoint(iSlat,1e3*lDeltaDiffRotSlat0[iSlat][1],iSlatToPos[iSlat]+1);
 768     gDeltaDiffPhiSlat0->SetPoint(iSlat,1e3*lDeltaDiffRotSlat0[iSlat][2],iSlatToPos[iSlat]+1);
 769   }
 770
 771   // Import TGeo geometry
 772   char* geoFilename = "geometry.root";
 773   cout << "geometry imported" << endl;
 774   if ( ! AliGeomManager::GetGeometry() ) {
 775     AliGeomManager::LoadGeometry(geoFilename);
 776     if (! AliGeomManager::GetGeometry() ) {
 777       printf("MUONSurveyCh8L: getting geometry from file %s failed\n", geoFilename);
 778       return;
 779     }
 780   }
 781
 782   AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer();
 783 //   transform->ReadGeometryData("volpath.dat", gGeoManager);
 784   transform->LoadGeometryData();
 785   cout << "geometry data read" << endl;
 786   AliMUONGeometryTransformer *newTransform = ReAlign(transform,11,dtrfSlat,true);
 787   newTransform->WriteTransformations("transform2ReAlign.dat");
 788   cout << "newtransform read" << endl;
 789   // Generate realigned data in local cdb
 790   const TClonesArray* array = newTransform->GetMisAlignmentData();
 791
 792   // CDB manager
 793   AliCDBManager* cdbManager = AliCDBManager::Instance();
 794   cdbManager->SetDefaultStorage("local://ReAlignCDB");
 795
 796   AliCDBMetaData* cdbData = new AliCDBMetaData();
 797   cdbData->SetResponsible("Dimuon Offline project");
 798   cdbData->SetComment("MUON alignment objects with residual misalignment");
 799   AliCDBId id("MUON/Align/Data", 0, 0);
 800   cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
 801
 802   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 803     cout << lDeltaDiffCenSlat0[iSlat][0] << " " << lDiffCenSlat0[iSlat][0] << " " << lDiffThCenSlat0[iSlat][0] << endl;
 804   }
 805   cout << endl;
 806   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 807     //    lDiffCenSlat0[iSlat][1] = lCenSlat[iSlat][1]-lCenSlat[6][1];
 808     //    lDeltaDiffCenSlat0[iSlat][1] = lDiffCenSlat0[iSlat][1]-lDiffThCenSlat0[iSlat][1];
 809     cout << lDeltaDiffCenSlat0[iSlat][1] << " " << lDiffCenSlat0[iSlat][1] << " " << lDiffThCenSlat0[iSlat][1] << endl;
 810   }
 811   cout << endl;
 812   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 813     cout << lDeltaDiffCenSlat0[iSlat][2] << " " << lDiffCenSlat0[iSlat][2] << " " << lDiffThCenSlat0[iSlat][2] << endl;
 814   }
 815   cout << endl;
 816   cout << endl;
 817   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 818     cout << lDeltaDiffRotSlat0[iSlat][0] << " " << lDiffRotSlat0[iSlat][0] << " " << lDiffThRotSlat0[iSlat][0] << endl;
 819   }
 820   cout << endl;
 821   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 822     cout << lDeltaDiffRotSlat0[iSlat][1] << " " << lDiffRotSlat0[iSlat][1] << " " << lDiffThRotSlat0[iSlat][1] << endl;
 823   }
 824   cout << endl;
 825   for(Int_t iSlat=0; iSlat<nSlats; iSlat++){
 826     cout << lDeltaDiffRotSlat0[iSlat][2] << " " << lDiffRotSlat0[iSlat][2] << " " << lDiffThRotSlat0[iSlat][2] << endl;
 827   }
 828
 829   TH1F *mySlatDeltaDiffCenX = new TH1F("mySlatDeltaDiffCenX","mySlatDeltaDiffCenX",100,-10,10);
 830   mySlatDeltaDiffCenX->SetMaximum(15);
 831   mySlatDeltaDiffCenX->SetMinimum(0);
 832   TH1F *mySlatDeltaDiffCenY = new TH1F("mySlatDeltaDiffCenY","mySlatDeltaDiffCenY",100,-10,10);
 833   mySlatDeltaDiffCenY->SetMaximum(15);
 834   mySlatDeltaDiffCenY->SetMinimum(0);
 835   TH1F *mySlatDeltaDiffCenZ = new TH1F("mySlatDeltaDiffCenZ","mySlatDeltaDiffCenZ",100,-20,20);
 836   mySlatDeltaDiffCenZ->SetMaximum(15);
 837   mySlatDeltaDiffCenZ->SetMinimum(0);
 838
 839   TH1F *mySlatDeltaDiffRotX = new TH1F("mySlatDeltaDiffRotX","mySlatDeltaDiffRotX",100,-10,10);
 840   mySlatDeltaDiffRotX->SetMaximum(15);
 841   mySlatDeltaDiffRotX->SetMinimum(0);
 842   TH1F *mySlatDeltaDiffRotY = new TH1F("mySlatDeltaDiffRotY","mySlatDeltaDiffRotY",100,-10,10);
 843   mySlatDeltaDiffRotY->SetMaximum(15);
 844   mySlatDeltaDiffRotY->SetMinimum(0);
 845   TH1F *mySlatDeltaDiffRotZ = new TH1F("mySlatDeltaDiffRotZ","mySlatDeltaDiffRotZ",100,-5,5);
 846   mySlatDeltaDiffRotZ->SetMaximum(15);
 847   mySlatDeltaDiffRotZ->SetMinimum(0);
 848   //
 849   // ******** Starting plots
 850   //
 851   TCanvas *canvas;
 852   TPad *pad;
 853   TPaveLabel *theTitle;
 854
 855   TPostScript *ps = 0;
 856
 857   if( saveps ){
 858     ps = new TPostScript(filename,filetype);
 859     ps->NewPage();
 860   }
 861
 862   // Inv Mass, Multiplicity
 863   sprintf(str,"Chamber 8L");
 864   TCanvas *cvn0 = new TCanvas("cvn0",str,cWidth,cHeight);
 865   canvas = cvn0;
 866   canvas->Range(0,0,21,29);
 867
 868   TPaveLabel *theTitle0 = new TPaveLabel(3,27.0,18,28.5," Deformations of chamber 8L ","br");
 869   theTitle = theTitle0;
 870   theTitle->SetFillColor(18);
 871   theTitle->SetTextFont(32);
 872   theTitle->SetTextSize(0.4);
 873   theTitle->SetTextColor(1);
 874   theTitle->Draw();
 875
 876   TPad *pad0 = new TPad("pad0","pad0",0.01,0.01,0.98,0.91,0);
 877   pad = pad0;
 878   pad->Draw();
 879   pad->Divide(2,2);
 880
 881   pad->cd(1);
 882   gStyle->SetPalette(1);
 883   hCPSTry->SetMinimum(100);
 884   hCPSTry->SetMaximum(120);
 885   hCPSTry->Draw("lego2z");
 886
 887   pad->cd(2);
 888   gStyle->SetPalette(1);
 889   hSSTry->SetMinimum(60);
 890   hSSTry->SetMaximum(80);
 891   hSSTry->Draw("lego2z");
 892
 893   pad->cd(3);
 894   gStyle->SetPalette(1);
 895   hCSTy->SetMinimum(110);
 896   hCSTy->SetMaximum(130);
 897   hCSTy->Draw("lego2z");
 898
 899   pad->cd(4);
 900   gStyle->SetPalette(1);
 901   hSSTly->SetMinimum(165);
 902   hSSTly->SetMaximum(185);
 903   hSSTly->Draw("lego2z");
 904
 905   if(saveps){
 906     ps->NewPage();
 907   }
 908
 909   // Inv Mass, Multiplicity
 910   sprintf(str,"Chamber 8L");
 911   TCanvas *cvn1 = new TCanvas("cvn1",str,cWidth,cHeight);
 912   canvas = cvn1;
 913   canvas->Range(0,0,21,29);
 914
 915   TPaveLabel *theTitle1 = new TPaveLabel(3,27.0,18,28.5," Deformations of chamber 8L ","br");
 916   theTitle = theTitle1;
 917   theTitle->SetFillColor(18);
 918   theTitle->SetTextFont(32);
 919   theTitle->SetTextSize(0.4);
 920   theTitle->SetTextColor(1);
 921   theTitle->Draw();
 922
 923   TPad *pad1 = new TPad("pad1","pad1",0.01,0.01,0.98,0.91,0);
 924   pad = pad1;
 925   pad->Draw();
 926   pad->Divide(2,2);
 927
 928   pad->cd(1);
 929   gStyle->SetPalette(1);
 930   hCPSTly->SetMinimum(120);
 931   hCPSTly->SetMaximum(140);
 932   hCPSTly->Draw("lego2z");
 933
 934   pad->cd(2);
 935   gStyle->SetPalette(1);
 936   hSSTrpy->SetMinimum(60);
 937   hSSTrpy->SetMaximum(80);
 938   hSSTrpy->Draw("surf2z");
 939
 940   pad->cd(3);
 941   gStyle->SetPalette(1);
 942   hCSTy->SetMinimum(110);
 943   hCSTy->SetMaximum(130);
 944   hCSTy->Draw("lego2z");
 945
 946   pad->cd(4);
 947   gStyle->SetPalette(1);
 948   hSSTlpy->SetMinimum(165);
 949   hSSTlpy->SetMaximum(185);
 950   hSSTlpy->Draw("surf2z");
 951
 952   // Inv Mass, Multiplicity
 953   sprintf(str,"Chamber 8L");
 954   TCanvas *cvn2 = new TCanvas("cvn2",str,cWidth,cHeight);
 955   canvas = cvn2;
 956   canvas->Range(0,0,21,29);
 957
 958   TPaveLabel *theTitle2 = new TPaveLabel(3,27.0,18,28.5," Deformations of chamber 8L ","br");
 959   theTitle = theTitle2;
 960   theTitle->SetFillColor(18);
 961   theTitle->SetTextFont(32);
 962   theTitle->SetTextSize(0.4);
 963   theTitle->SetTextColor(1);
 964   theTitle->Draw();
 965
 966   TPad *pad2 = new TPad("pad2","pad2",0.01,0.01,0.98,0.91,0);
 967   pad = pad2;
 968   pad->Draw();
 969   pad->Divide(3,2);
 970
 971   pad->cd(1);
 972   mySlatDeltaDiffCenX->Draw();
 973   mySlatDeltaDiffCenX->SetXTitle("#Delta[(xc_{i}^{m}-xc_{0}^{m})-(xc_{i}^{th}-xc_{0}^{th})] (mm)");
 974   mySlatDeltaDiffCenX->SetYTitle("Slat ascending vertical ordering");
 975   gDeltaDiffCenXSlat0->SetMarkerStyle(20);
 976   gDeltaDiffCenXSlat0->Draw("P");
 977
 978   pad->cd(2);
 979   mySlatDeltaDiffCenY->Draw();
 980   mySlatDeltaDiffCenY->SetXTitle("#Delta[(yc_{i}^{m}-yc_{0}^{m})-(yc_{i}^{th}-yc_{0}^{th})] (mm)");
 981   mySlatDeltaDiffCenY->SetYTitle("Slat ascending vertical ordering");
 982   gDeltaDiffCenYSlat0->SetMarkerStyle(20);
 983   gDeltaDiffCenYSlat0->Draw("P");
 984
 985   pad->cd(3);
 986   mySlatDeltaDiffCenZ->Draw();
 987   mySlatDeltaDiffCenZ->SetXTitle("#Delta[(zc_{i}^{m}-zc_{0}^{m})-(zc_{i}^{th}-zc_{0}^{th})] (mm)");
 988   mySlatDeltaDiffCenZ->SetYTitle("Slat ascending vertical ordering");
 989   gDeltaDiffCenZSlat0->SetMarkerStyle(20);
 990   gDeltaDiffCenZSlat0->Draw("P");
 991
 992   pad->cd(4);
 993   mySlatDeltaDiffRotX->Draw();
 994   mySlatDeltaDiffRotX->SetXTitle("#Delta[(#psi_{i}^{m}-#psi_{0}^{m})-(#psi_{i}^{th}-#psi_{0}^{th})] (mrad)");
 995   mySlatDeltaDiffRotX->SetYTitle("Slat ascending vertical ordering");
 996   gDeltaDiffPsiSlat0->SetMarkerStyle(20);
 997   gDeltaDiffPsiSlat0->Draw("P");
 998
 999   pad->cd(5);
1000   mySlatDeltaDiffRotY->Draw();
1001   mySlatDeltaDiffRotY->SetXTitle("#Delta[(#theta_{i}^{m}-#theta_{0}^{m})-(#theta_{i}^{th}-#theta_{0}^{th})] (mrad)");
1002   mySlatDeltaDiffRotY->SetYTitle("Slat ascending vertical ordering");
1003   gDeltaDiffThtSlat0->SetMarkerStyle(20);
1004   gDeltaDiffThtSlat0->Draw("P");
1005
1006   pad->cd(6);
1007   mySlatDeltaDiffRotZ->Draw();
1008   mySlatDeltaDiffRotZ->SetXTitle("#Delta[(#phi_{i}^{m}-#phi_{0}^{m})-(#phi_{i}^{th}-#phi_{0}^{th})] (mrad)");
1009   mySlatDeltaDiffRotZ->SetYTitle("Slat ascending vertical ordering");
1010   gDeltaDiffPhiSlat0->SetMarkerStyle(20);
1011   gDeltaDiffPhiSlat0->Draw("P");
1012
1013   if( saveps ){
1014     ps->Close();
1015   }
1016 }