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