ESD file is added to the list of proof output files. It is then automatically merged...
[u/mrichter/AliRoot.git] / MUON / MUONSurveyCh3.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 3
21 /// 
22 /// Macro loads the survey data from .txt file using AliSurveyObj.
23 /// Macro uses AliMUONSurvey... classes.
24 /// The transformations of the detection elements or chambers can be obtained 
25 /// in three ways:
26 /// - A: Fit of local to global transformation using the fixed button targets.
27 /// - B: Fit a plane to the sticker targets -> psi, theta
28 ///      Use above psi and theta and fit remaining 4 parameters using the fixed 
29 ///      button targets
30 /// - C: Fit a plane to the sticker targets -> psi, theta
31 ///      Use above psi and theta to calculate xc, yc, zc and phi by solving 
32 ///      the equations from a local to global transformation of the fixed button 
33 ///      targets
34 ///
35 /// Methods A and B are prefered to C, and B is better if sticker targets are 
36 /// available and lie on a plane!
37 /// For slats only methods B and C can be used.
38 /// Various histograms are filled and printed for monitoring.
39 /// MisAlignment object is then created.
40 ///
41 /// \author Javier Castillo
42
43 #if !defined(__CINT__) || defined(__MAKECINT__)
44
45 #include "AliMUONGeometryTransformer.h"
46 #include "AliMUONGeometryModuleTransformer.h"
47 #include "AliMUONGeometryDetElement.h"
48 #include "AliMUONGeometryMisAligner.h"
49 #include "AliMUONSurveyChamber.h"
50 #include "AliMUONSurveyDetElem.h"
51 #include "AliMUONSurveyUtil.h"
52
53 #include "AliSurveyObj.h"
54 #include "AliSurveyPoint.h"
55 #include "AliGeomManager.h"
56 #include "AliCDBManager.h"
57 #include "AliCDBMetaData.h"
58 #include "AliCDBId.h"
59 #include "AliAlignObjMatrix.h"
60 #include "AliAlignObj.h"
61
62 #include <TROOT.h>
63 #include <TGeoManager.h>
64 #include <TClonesArray.h>
65 #include <TObjArray.h>
66 #include <TArrayD.h>
67 #include <TMath.h>
68 #include <TString.h>
69 #include <TFitter.h>
70 #include <TH2.h>
71 #include <TF2.h>
72 #include <TGraphErrors.h>
73 #include <TCanvas.h>
74 #include <TPad.h>
75 #include <TLine.h>
76 #include <TPostScript.h>
77 #include <TPaveLabel.h>
78 #include <TStyle.h>
79 #include <TFile.h>
80 #include <TMatrixDSym.h>
81
82 #include <fstream>
83
84 #endif
85
86 void MUONSurveyCh3() {
87   
88   TString str;
89   TString title;
90   
91   Bool_t bMonitor = kTRUE;
92   Bool_t bOCDB = kTRUE;
93   Bool_t saveps = kFALSE;
94   const int cWidth = (int)(700*(29./21.));
95   const int cHeight = 700;
96   const int filetype  = 112; // landscape  
97
98   Int_t chId = 2;
99   Int_t nChs = 2;
100   Int_t nDetElems = 4;
101   Int_t nDetElemsI = 2;
102   //  Int_t nDetElemsO = 2;
103   Int_t iDetElemToDetElemId[4] = {300,301,302,303};
104   Int_t iDetElemPseudoIdToDetElem[4] = {0,1,2,3};
105   Int_t iDetElemsOfChamber[2][2] = {{0,3},
106                                     {1,2}};
107   
108   TObjArray *myChamberArray = new TObjArray(); 
109   myChamberArray->Add(new AliMUONSurveyChamber(chId));  
110   myChamberArray->Add(new AliMUONSurveyChamber(chId));  
111
112   AliMUONSurveyChamber *myChamberI = 0x0;
113   AliMUONSurveyChamber *myChamberO = 0x0;
114   AliMUONSurveyChamber *myChamber = 0x0;
115   AliMUONSurveyDetElem *myDetElem = 0x0;
116
117   myChamberI = (AliMUONSurveyChamber*)myChamberArray->At(0);
118   myChamberI->GetSurveyObj()->FillFromLocalFile("../Reports/AliceSt2_TC3_2571e.txt");
119   myChamberO = (AliMUONSurveyChamber*)myChamberArray->At(1);
120   myChamberO->GetSurveyObj()->FillFromLocalFile("../Reports/AliceSt2_TC3_2571e.txt");
121
122   myChamber = myChamberI;
123   myChamber->PrintSurveyReport();
124
125   // Chamber & DetElems button targets local coordinates
126   AliSurveyObj *lSO = new AliSurveyObj();    
127   lSO->FillFromLocalFile("$ALICE_ROOT/MUON/data/MUONTrackingTargetsLocal.txt");
128
129   // Set survey targets for chambers 
130   myChamberI->AddGButtonTargets(TString("T1_3520"),3);
131   myChamberI->AddStickerTargets(TString("T1_3620"),6);
132   myChamberI->AddLButtonTargets(lSO->GetData(),Form("%d_520",chId+1),3);
133   myChamberI->AddGButtonTargets(TString("T1_3510"),3);
134   myChamberI->AddStickerTargets(TString("T1_3610"),6);
135   myChamberI->AddLButtonTargets(lSO->GetData(),Form("%d_510",chId+1),3);
136   myChamberO->AddGButtonTargets(TString("T1_3520"),3);
137   myChamberO->AddStickerTargets(TString("T1_3620"),6);
138   myChamberO->AddLButtonTargets(lSO->GetData(),Form("%d_520",chId+1),3);
139   myChamberO->AddGButtonTargets(TString("T1_3510"),3);
140   myChamberO->AddStickerTargets(TString("T1_3610"),6);
141   myChamberO->AddLButtonTargets(lSO->GetData(),Form("%d_510",chId+1),3);
142
143   // Set survey targets for detection elements 
144   for (int iCh =0; iCh<=1; iCh++) {
145     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
146     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
147       Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
148       myChamber->AddSurveyDetElem(iDetElemToDetElemId[iDetElem]);
149       TString baseName;
150       myDetElem =  myChamber->GetDetElem(iDetElemI);
151       if (myDetElem) {
152         baseName = Form("T1_353%d",iDetElem+1);
153         myDetElem->AddGButtonTargets(baseName,3);
154         baseName = Form("T1_363%d",iDetElem+1);
155         myDetElem->AddStickerTargets(baseName,6);
156         baseName = Form("T2_363%d",iDetElem+1);
157         myDetElem->AddStickerTargets(baseName,6);
158         baseName = Form("%d_53%d",chId+1,iDetElem+1);
159         myDetElem->AddLButtonTargets(lSO->GetData(),baseName,3);
160       }
161     }
162   }
163   printf("All targets added! \n");
164
165   Double_t **lCenDetElem = new Double_t*[nDetElems+nChs];
166   Double_t **lRotDetElem = new Double_t*[nDetElems+nChs];
167   Double_t **lDiffCenDetElem0 = new Double_t*[nDetElems+nChs];
168   Double_t **lDiffRotDetElem0 = new Double_t*[nDetElems+nChs];
169   Double_t **lDiffThCenDetElem0 = new Double_t*[nDetElems+nChs];
170   Double_t **lDiffThRotDetElem0 = new Double_t*[nDetElems+nChs];
171   Double_t **lDeltaDiffCenDetElem0 = new Double_t*[nDetElems+nChs];
172   Double_t **lDeltaDiffRotDetElem0 = new Double_t*[nDetElems+nChs];
173
174   for (int iDetElem=0; iDetElem<nDetElems+nChs; iDetElem++) {
175     lCenDetElem[iDetElem] = new Double_t[3];
176     lRotDetElem[iDetElem] = new Double_t[3];
177     lDiffCenDetElem0[iDetElem] = new Double_t[3];
178     lDiffRotDetElem0[iDetElem] = new Double_t[3];
179     lDiffThCenDetElem0[iDetElem] = new Double_t[3];
180     lDiffThRotDetElem0[iDetElem] = new Double_t[3];
181     lDeltaDiffCenDetElem0[iDetElem] = new Double_t[3];
182     lDeltaDiffRotDetElem0[iDetElem] = new Double_t[3];
183   }
184
185   TGeoCombiTrans dtrfDetElem[nDetElems+nChs];
186   TGeoCombiTrans localTrfDetElem[nDetElems+nChs];
187   TGeoCombiTrans localTrfThDetElem[nDetElems+nChs];
188
189   // Import TGeo geometry 
190   char* geoFilename = "geometry.root";
191   if ( ! AliGeomManager::GetGeometry() ) {
192     AliGeomManager::LoadGeometry(geoFilename);
193     if (! AliGeomManager::GetGeometry() ) {
194       printf("MUONSurveyCh%d: getting geometry from file %s failed\n", chId+1,geoFilename);
195       return;
196     }
197     cout << "geometry imported" << endl;
198   }
199
200   AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer();
201   transform->LoadGeometryData();
202
203   TGeoCombiTrans trfThChamber;
204   TGeoCombiTrans trfThDetElem;
205
206   for (int iCh =0; iCh<=1; iCh++) {
207     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
208
209     trfThChamber = TGeoCombiTrans(*transform->GetModuleTransformerByDEId(iDetElemToDetElemId[iCh*nDetElemsI])->GetTransformation());
210     trfThChamber.SetTranslation(trfThChamber.GetTranslation()[0]*10,
211                                 trfThChamber.GetTranslation()[1]*10,
212                                 trfThChamber.GetTranslation()[2]*10);
213     trfThChamber.Print();
214     myChamber->SetLocalTransformation(new TGeoCombiTrans(trfThChamber),kTRUE);
215
216     // Set Chamber plane function
217     cout << "Setting plane for Chamber" << iCh+1 << " ..." << endl;
218     myChamber->SetPlane(Form("fChamber%d",iCh+1));
219     myChamber->SetPlaneParameters(0.,0.,0.);
220
221     // Fit a plane to sticker targets
222     Double_t lCChi2 = myChamber->FitPlane();
223     cout << "... chi2 = " << lCChi2 << " ..." << endl; 
224     
225     // Get best transformation from fitting method 
226     // (local to global transformation)
227     cout << "Trying fitting method for chamber " << iCh << endl;
228     myChamber->SurveyToAlign();    
229     //    myChamber->SurveyToAlign(myChamber->GetPlane()->GetParameter(0),myChamber->GetPlane()->GetParameter(1),myChamber->GetPlane()->GetParError(0),myChamber->GetPlane()->GetParError(1));    
230     
231     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
232       myDetElem =  myChamber->GetDetElem(iDetElemI);
233       Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
234
235       trfThDetElem.Clear();
236       trfThDetElem = TGeoCombiTrans(*transform->GetDetElement(iDetElemToDetElemId[iDetElem])->GetLocalTransformation());
237       trfThDetElem.SetTranslation(trfThDetElem.GetTranslation()[0]*10,
238                                   trfThDetElem.GetTranslation()[1]*10,
239                                   trfThDetElem.GetTranslation()[2]*10);
240       trfThDetElem.Print();
241
242       myDetElem->SetLocalTransformation(new TGeoCombiTrans(trfThDetElem),kTRUE);
243
244       for (int iCor=0; iCor<3; iCor++){
245         lCenDetElem[iDetElem][iCor] = 0;
246         lRotDetElem[iDetElem][iCor] = 0;
247       }
248
249       if (bMonitor){
250         // MONITOR: Draw graph with sticker targets for plane fitting
251         myDetElem->DrawSTargets();
252         gPad->Update();
253       }
254
255       // Get DetElem transformation. 
256       // Psi and Theta are obtained by fitting a plane to the sticker targets.
257       // Then Xc, Yc, Zc and Phi are obtained by solving the equations to the ref.
258       // syst. transformation of the button targets
259
260       // Set DetElem plane function
261       cout << "Setting plane for DetElem" << iDetElem+1 << " ..." << endl;
262       myDetElem->SetPlane(Form("fDetElem%d",iDetElem+1));
263       myDetElem->SetPlaneParameters(0.,0.,30.)
264 ;
265       // Fit a plane to sticker targets
266       Double_t lChi2 = myDetElem->FitPlane();
267       cout << "... chi2 = " << lChi2 << " ..." << endl; 
268       
269       lRotDetElem[iDetElem][0]=TMath::ATan(myDetElem->GetPlane()->GetParameter(0));
270       lRotDetElem[iDetElem][1]=TMath::ATan(myDetElem->GetPlane()->GetParameter(1));
271       
272       // Calculate Mean transformation using previous plane fit 
273       // and pairs of button targets
274       myDetElem->CalculateMeanTransf(lCenDetElem[iDetElem],lRotDetElem[iDetElem]);
275     
276       cout << "DetElem" << iDetElem+1 << " : mycen(" << lCenDetElem[iDetElem][0] << "," << lCenDetElem[iDetElem][1] << "," << lCenDetElem[iDetElem][2] << "); rot(" << lRotDetElem[iDetElem][0] << "," << lRotDetElem[iDetElem][1] << "," << lRotDetElem[iDetElem][2] << ")"  << endl;          
277       
278      
279       // Get best transformation from fitting method 
280       // (local to global transformation)
281       cout << "Trying fitting method for DetElemId " << iDetElemToDetElemId[iDetElem] << endl;
282       //      myDetElem->SurveyToAlign();     
283       myDetElem->SurveyToAlign(lRotDetElem[iDetElem][0],lRotDetElem[iDetElem][1],myDetElem->GetPlane()->GetParError(0),myDetElem->GetPlane()->GetParError(1));    
284     }
285   }
286
287   // Print found transformation of Detection Element (plane fit + loop)        
288   for (int iCh =0; iCh<=1; iCh++) {
289     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
290       Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
291       cout << "DetElem" << iDetElem+1 << " : mycen(" << lCenDetElem[iDetElem][0] << "," << lCenDetElem[iDetElem][1] << "," << lCenDetElem[iDetElem][2] << "); rot(" << lRotDetElem[iDetElem][0] << "," << lRotDetElem[iDetElem][1] << "," << lRotDetElem[iDetElem][2] << ")"  << endl;          
292     }
293   }
294
295   // Print Theoretical transformation of Detection Element
296   for (int iCh =0; iCh<=1; iCh++) {
297     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
298     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
299       myChamber->GetDetElem(iDetElemI)->PrintLocalTrf();
300     }
301   }
302
303   // Print found delta transformation of Detection Element
304   for (int iCh =0; iCh<=1; iCh++) {
305     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);      
306     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
307       myChamber->GetDetElem(iDetElemI)->PrintAlignTrf();
308     }
309   }
310
311   //
312   // Compare transformations to expected ones 
313   //
314   Int_t iDetElemToPos[18] = {0, 1, 2, 3};
315  
316   TGraphErrors *gDeltaDiffCenXDetElem0 = new TGraphErrors(nDetElems);
317   TGraphErrors *gDeltaDiffCenYDetElem0 = new TGraphErrors(nDetElems);
318   TGraphErrors *gDeltaDiffCenZDetElem0 = new TGraphErrors(nDetElems);
319   TGraphErrors *gDeltaDiffPsiDetElem0 = new TGraphErrors(nDetElems);
320   TGraphErrors *gDeltaDiffThtDetElem0 = new TGraphErrors(nDetElems);
321   TGraphErrors *gDeltaDiffPhiDetElem0 = new TGraphErrors(nDetElems);
322
323   for (int iCh =0; iCh<=1; iCh++) {
324     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);      
325     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
326       myChamber->GetDetElem(iDetElemI)->GetAlignTrf()->Print();
327     }
328   }
329
330   for (int iCh =0; iCh<=1; iCh++) {
331     myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
332     // Store delta transformations to use for CDB entry creation
333     dtrfDetElem[nDetElems+iCh].Clear();
334     dtrfDetElem[nDetElems+iCh] = *(myChamber->GetAlignTrf());
335
336     // Those are for checks and visualizations
337     localTrfDetElem[nDetElems+iCh].Clear();
338     localTrfDetElem[nDetElems+iCh] = (*(myChamber->GetLocalTrf())*(*(myChamber->GetAlignTrf())));
339     localTrfDetElem[nDetElems+iCh].Print();
340     lDiffCenDetElem0[nDetElems+iCh] = (Double_t*)localTrfDetElem[nDetElems+iCh].GetTranslation();
341     AliMUONSurveyUtil::MatrixToAngles(localTrfDetElem[nDetElems+iCh].GetRotationMatrix(),lDiffRotDetElem0[nDetElems+iCh]);
342     
343     localTrfThDetElem[nDetElems+iCh].Clear();
344     localTrfThDetElem[nDetElems+iCh] = (*(myChamber->GetLocalTrf()));
345     localTrfThDetElem[nDetElems+iCh].Print();
346     lDiffThCenDetElem0[nDetElems+iCh] = (Double_t*)localTrfThDetElem[nDetElems+iCh].GetTranslation();
347     AliMUONSurveyUtil::MatrixToAngles(localTrfThDetElem[nDetElems+iCh].GetRotationMatrix(),lDiffThRotDetElem0[nDetElems+iCh]);
348
349     for (int iCor=0; iCor<3; iCor++){
350       lDeltaDiffCenDetElem0[nDetElems+iCh][iCor] = lDiffCenDetElem0[nDetElems+iCh][iCor]-lDiffThCenDetElem0[nDetElems+iCh][iCor];
351       lDeltaDiffRotDetElem0[nDetElems+iCh][iCor] = lDiffRotDetElem0[nDetElems+iCh][iCor]-lDiffThRotDetElem0[nDetElems+iCh][iCor];
352       if (lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]>TMath::Pi()) lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]-=TMath::TwoPi();
353       if (lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]<-TMath::Pi()) lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]+=TMath::TwoPi();
354     }      
355
356     for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
357       Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
358       myDetElem =  myChamber->GetDetElem(iDetElemI);
359       // Store delta transformations to use for CDB entry creation
360       dtrfDetElem[iDetElem].Clear();
361       dtrfDetElem[iDetElem] = *(myDetElem->GetAlignTrf());
362
363       // Those are for checks and visualizations
364       localTrfDetElem[iDetElem].Clear();
365       localTrfDetElem[iDetElem] = (*(myDetElem->GetLocalTrf())*(*(myDetElem->GetAlignTrf())));
366       //      localTrfDetElem[iDetElem] = (*(myDetElem->GetBaseTrf())*(*(myDetElem->GetAlignTrf())));
367       localTrfDetElem[iDetElem].Print();
368       lDiffCenDetElem0[iDetElem] = (Double_t*)localTrfDetElem[iDetElem].GetTranslation();
369       AliMUONSurveyUtil::MatrixToAngles(localTrfDetElem[iDetElem].GetRotationMatrix(),lDiffRotDetElem0[iDetElem]);
370 //       lDiffCenDetElem0[iDetElem] = lCenDetElem[iDetElem];
371 //       lDiffRotDetElem0[iDetElem] = lRotDetElem[iDetElem];
372
373       localTrfThDetElem[iDetElem].Clear();
374       localTrfThDetElem[iDetElem] = (*(myDetElem->GetLocalTrf()));
375       //      localTrfThDetElem[iDetElem] = (*(myDetElem->GetBaseTrf()));
376       localTrfThDetElem[iDetElem].Print();
377       lDiffThCenDetElem0[iDetElem] = (Double_t*)localTrfThDetElem[iDetElem].GetTranslation();
378       AliMUONSurveyUtil::MatrixToAngles(localTrfThDetElem[iDetElem].GetRotationMatrix(),lDiffThRotDetElem0[iDetElem]);
379       
380       for (int iCor=0; iCor<3; iCor++){
381         lDeltaDiffCenDetElem0[iDetElem][iCor] = lDiffCenDetElem0[iDetElem][iCor]-lDiffThCenDetElem0[iDetElem][iCor];
382         lDeltaDiffRotDetElem0[iDetElem][iCor] = lDiffRotDetElem0[iDetElem][iCor]-lDiffThRotDetElem0[iDetElem][iCor];
383         if (lDeltaDiffRotDetElem0[iDetElem][iCor]>TMath::Pi()) lDeltaDiffRotDetElem0[iDetElem][iCor]-=TMath::TwoPi();
384         if (lDeltaDiffRotDetElem0[iDetElem][iCor]<-TMath::Pi()) lDeltaDiffRotDetElem0[iDetElem][iCor]+=TMath::TwoPi();
385       }
386
387       gDeltaDiffCenXDetElem0->SetPoint(iDetElem,lDeltaDiffCenDetElem0[iDetElem][0],iDetElemToPos[iDetElem]+1);
388       gDeltaDiffCenYDetElem0->SetPoint(iDetElem,lDeltaDiffCenDetElem0[iDetElem][1],iDetElemToPos[iDetElem]+1);
389       gDeltaDiffCenZDetElem0->SetPoint(iDetElem,lDeltaDiffCenDetElem0[iDetElem][2],iDetElemToPos[iDetElem]+1);
390       gDeltaDiffPsiDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][0],iDetElemToPos[iDetElem]+1);
391       gDeltaDiffThtDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][1],iDetElemToPos[iDetElem]+1);
392       gDeltaDiffPhiDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][2],iDetElemToPos[iDetElem]+1);
393       gDeltaDiffCenXDetElem0->SetPointError(iDetElem,myDetElem->GetFitter()->GetParError(0),0.);
394       gDeltaDiffCenYDetElem0->SetPointError(iDetElem,myDetElem->GetFitter()->GetParError(1),0.);
395       gDeltaDiffCenZDetElem0->SetPointError(iDetElem,myDetElem->GetFitter()->GetParError(2),0.);
396       gDeltaDiffPsiDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(3),0.);
397       gDeltaDiffThtDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(4),0.);
398       gDeltaDiffPhiDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(5),0.);  
399     }
400   }
401
402   // Apply the found misalignments to the geometry
403   AliMUONGeometryTransformer *newTransform = AliMUONSurveyUtil::ReAlign(transform,chId,nDetElems,iDetElemPseudoIdToDetElem,dtrfDetElem,true); 
404   newTransform->WriteTransformations(Form("transform2ReAlignSurveyCh%d.dat",chId+1));
405
406   if(bOCDB){
407     // Generate realigned data in local cdb
408     const TClonesArray* array = newTransform->GetMisAlignmentData();
409     
410     // Set the alignment resolution in the align objects for this chamber
411     Double_t chResX = 0.1*myChamberI->GetAlignResX();
412     Double_t chResY = 0.1*myChamberI->GetAlignResY();
413     Double_t deResX = 0.1*(myChamberI->GetMeanDetElemAlignResX()+myChamberO->GetMeanDetElemAlignResX())/2.;
414     Double_t deResY = 0.1*(myChamberI->GetMeanDetElemAlignResY()+myChamberO->GetMeanDetElemAlignResY())/2.;
415     printf("Chamber alignment resolutions: resX=%f , resY=%f\n",chResX,chResY); 
416     printf("Detection Elements alignment resolutions: resX=%f , resY=%f\n",deResX,deResY); 
417     chResX = TMath::Sqrt(0.1*0.1+chResX*chResX);
418     chResY = TMath::Sqrt(0.1*0.1+chResY*chResY);
419     AliMUONSurveyUtil::SetAlignmentResolution(array,chId,chResX,chResY,deResX,deResY);
420     
421     // CDB manager
422     AliCDBManager* cdbManager = AliCDBManager::Instance();
423     cdbManager->SetDefaultStorage(Form("local://ReAlignSurveyCh%dCDB",chId+1));
424     
425     AliCDBMetaData* cdbData = new AliCDBMetaData();
426     cdbData->SetResponsible("Dimuon Offline project");
427     cdbData->SetComment("MUON alignment objects with survey realignment");
428     AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity()); 
429     cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
430   }
431
432   for(Int_t iCor=0; iCor<3; iCor++){
433     for(Int_t iDetElem=0; iDetElem<nDetElems; iDetElem++){
434       cout << lDeltaDiffCenDetElem0[iDetElem][iCor] << " " << lDiffCenDetElem0[iDetElem][iCor] << " " << lDiffThCenDetElem0[iDetElem][iCor] << endl;
435     }
436     cout << endl;
437   }
438   cout << endl;
439   for(Int_t iCor=0; iCor<3; iCor++){
440     for(Int_t iDetElem=0; iDetElem<nDetElems; iDetElem++){
441       cout << lDeltaDiffRotDetElem0[iDetElem][iCor] << " " << lDiffRotDetElem0[iDetElem][iCor] << " " << lDiffThRotDetElem0[iDetElem][iCor] << endl;
442     }
443     cout << endl;
444   }
445
446   TH1F *myDetElemDeltaDiffCenX = new TH1F("myDetElemDeltaDiffCenX","myDetElemDeltaDiffCenX",100,-10,10);
447   myDetElemDeltaDiffCenX->SetMaximum(nDetElems+1);
448   myDetElemDeltaDiffCenX->SetMinimum(0);
449   TH1F *myDetElemDeltaDiffCenY = new TH1F("myDetElemDeltaDiffCenY","myDetElemDeltaDiffCenY",100,-10,10);
450   myDetElemDeltaDiffCenY->SetMaximum(nDetElems+1);
451   myDetElemDeltaDiffCenY->SetMinimum(0);
452   TH1F *myDetElemDeltaDiffCenZ = new TH1F("myDetElemDeltaDiffCenZ","myDetElemDeltaDiffCenZ",100,-20,20);
453   myDetElemDeltaDiffCenZ->SetMaximum(nDetElems+1);
454   myDetElemDeltaDiffCenZ->SetMinimum(0);
455
456   TH1F *myDetElemDeltaDiffRotX = new TH1F("myDetElemDeltaDiffRotX","myDetElemDeltaDiffRotX",100,-15,15);
457   myDetElemDeltaDiffRotX->SetMaximum(nDetElems+1);
458   myDetElemDeltaDiffRotX->SetMinimum(0);
459   TH1F *myDetElemDeltaDiffRotY = new TH1F("myDetElemDeltaDiffRotY","myDetElemDeltaDiffRotY",100,-15,15);
460   myDetElemDeltaDiffRotY->SetMaximum(nDetElems+1);
461   myDetElemDeltaDiffRotY->SetMinimum(0);
462   TH1F *myDetElemDeltaDiffRotZ = new TH1F("myDetElemDeltaDiffRotZ","myDetElemDeltaDiffRotZ",100,-5,5);
463   myDetElemDeltaDiffRotZ->SetMaximum(nDetElems+1);
464   myDetElemDeltaDiffRotZ->SetMinimum(0);
465
466   //
467   // ******** Starting plots 
468   //
469   TCanvas *canvas;
470   TPad *pad;
471   TPaveLabel *theTitle;
472   gStyle->SetPalette(1);
473
474   TPostScript *ps = 0;
475   
476   if( saveps ){
477     ps = new TPostScript(Form("surveyChamber%d",chId+1),filetype); 
478     ps->NewPage();
479   }
480
481   // Observed misalignments
482   str = Form("Chamber %d",chId+1);
483   TCanvas *cvn2 = new TCanvas("cvn2",str,cWidth,cHeight);
484   canvas = cvn2;
485   canvas->Range(0,0,21,29);
486   
487   title = Form(" MisAlignments Chamber %d - PL2G - In Single Frame - DE2Ok - DE3Ok - DE4Ok ",chId+1);
488   TPaveLabel *theTitle2 = new TPaveLabel(3,27.0,18,28.5,title,"br");
489   theTitle = theTitle2;
490   theTitle->SetFillColor(18);
491   theTitle->SetTextFont(32);
492   theTitle->SetTextSize(0.4);
493   theTitle->SetTextColor(1);
494   theTitle->Draw();
495  
496   TPad *pad2 = new TPad("pad2","pad2",0.01,0.01,0.98,0.91,0);
497   pad = pad2;
498   pad->Draw();
499   TLine *ch0Line = new TLine(0,1,0,2);
500   TLine *ch1Line = new TLine(0,1,0,2);
501   ch1Line->SetLineStyle(2);
502   pad->Divide(3,2);
503   
504   pad->cd(1);
505   myDetElemDeltaDiffCenX->Draw();
506   myDetElemDeltaDiffCenX->SetXTitle("#Delta[xc_{i}^{m}-xc_{i}^{th}] (mm)");
507   myDetElemDeltaDiffCenX->SetYTitle("DetElem arbitrary ordering");
508   gDeltaDiffCenXDetElem0->SetMarkerStyle(20);
509   gDeltaDiffCenXDetElem0->Draw("P");
510   ch0Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+0][0],0.5,lDeltaDiffCenDetElem0[nDetElems+0][0],1.5);
511   ch0Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+0][0],3.5,lDeltaDiffCenDetElem0[nDetElems+0][0],4.5);
512   ch1Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+1][0],1.5,lDeltaDiffCenDetElem0[nDetElems+1][0],3.5);
513
514   pad->cd(2);
515   myDetElemDeltaDiffCenY->Draw();
516   myDetElemDeltaDiffCenY->SetXTitle("#Delta[yc_{i}^{m}-yc_{i}^{th}] (mm)");
517   myDetElemDeltaDiffCenY->SetYTitle("DetElem arbitrary ordering");
518   gDeltaDiffCenYDetElem0->SetMarkerStyle(20);
519   gDeltaDiffCenYDetElem0->Draw("P");
520   ch0Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+0][1],0.5,lDeltaDiffCenDetElem0[nDetElems+0][1],1.5);
521   ch0Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+0][1],3.5,lDeltaDiffCenDetElem0[nDetElems+0][1],4.5);
522   ch1Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+1][1],1.5,lDeltaDiffCenDetElem0[nDetElems+1][1],3.5);
523
524   pad->cd(3);
525   myDetElemDeltaDiffCenZ->Draw();
526   myDetElemDeltaDiffCenZ->SetXTitle("#Delta[zc_{i}^{m}-zc_{i}^{th}] (mm)");
527   myDetElemDeltaDiffCenZ->SetYTitle("DetElem arbitrary ordering");
528   gDeltaDiffCenZDetElem0->SetMarkerStyle(20);
529   gDeltaDiffCenZDetElem0->Draw("P");
530   ch0Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+0][2],0.5,lDeltaDiffCenDetElem0[nDetElems+0][2],1.5);
531   ch0Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+0][2],3.5,lDeltaDiffCenDetElem0[nDetElems+0][2],4.5);
532   ch1Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+1][2],1.5,lDeltaDiffCenDetElem0[nDetElems+1][2],3.5);
533
534   pad->cd(4);
535   myDetElemDeltaDiffRotX->Draw();
536   myDetElemDeltaDiffRotX->SetXTitle("#Delta[#psi_{i}^{m}-#psi_{i}^{th}] (mrad)");
537   myDetElemDeltaDiffRotX->SetYTitle("DetElem arbitrary ordering");
538   gDeltaDiffPsiDetElem0->SetMarkerStyle(20);
539   gDeltaDiffPsiDetElem0->Draw("P");
540   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],1.5);
541   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],3.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],4.5);
542   ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][0],1.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][0],3.5);
543
544   pad->cd(5);
545   myDetElemDeltaDiffRotY->Draw();
546   myDetElemDeltaDiffRotY->SetXTitle("#Delta[#theta_{i}^{m}-#theta_{i}^{th}] (mrad)");
547   myDetElemDeltaDiffRotY->SetYTitle("DetElem arbitrary ordering");
548   gDeltaDiffThtDetElem0->SetMarkerStyle(20);
549   gDeltaDiffThtDetElem0->Draw("P");
550   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],1.5);
551   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],3.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],4.5);
552   ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][1],1.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][1],3.5);
553
554   pad->cd(6);
555   myDetElemDeltaDiffRotZ->Draw();
556   myDetElemDeltaDiffRotZ->SetXTitle("#Delta[#phi_{i}^{m}-#phi_{i}^{th}] (mrad)");
557   myDetElemDeltaDiffRotZ->SetYTitle("DetElem arbitrary ordering");
558   gDeltaDiffPhiDetElem0->SetMarkerStyle(20);
559   gDeltaDiffPhiDetElem0->Draw("P");
560   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],1.5);
561   ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],3.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],4.5);
562   ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][2],1.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][2],3.5);
563
564   pad->Update();
565  
566   if(bMonitor){    
567     // MONITOR: Histograms for monitoring
568     TH2F *hCPSTa = new TH2F("hCPSTa","hCPSTa",100,-2000,2000,100,-2000,2000);
569     TH2F *hCPSTc = new TH2F("hCPSTc","hCPSTc",100,-2000,2000,100,-2000,2000);
570     TH2F *hSSTa = new TH2F("hSSTa","hSSTa",100,-2000,2000,100,-2000,2000);
571     TH2F *hSSTc = new TH2F("hSSTc","hSSTc",100,-2000,2000,100,-2000,2000);
572     TH2F *hSSTap = new TH2F("hSSTap","hSSTap",800,-2000,2000,800,-2000,2000);
573     TH2F *hSSTcp = new TH2F("hSSTcp","hSSTcp",800,-2000,2000,800,-2000,2000);
574     
575     // MONITOR: Fill histograms with chambers and slats sticker target positions
576     for (int iCh =0; iCh<=1; iCh++) {
577       myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);            
578       cout << "Filling histograms of sticker target for chamber" << iCh+1 << " ..." << endl;
579       myChamber->FillCPSTHistograms(Form("T1_36%d",2-iCh),hCPSTc);
580       myChamber->FillDESTHistograms(TString("T1_36"),hSSTc,TString("T2_36"),hSSTa);
581
582       for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
583         //      Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
584         myDetElem =  myChamber->GetDetElem(iDetElemI);
585         // MONITOR: Fill slat plane for fit monitor.
586         Double_t pl[3] = {0};
587         Double_t pg[3] = {0};
588         AliSurveyPoint *pointSBT0 = myDetElem->GetLButtonTarget(0);
589         AliSurveyPoint *pointSBT1 = myDetElem->GetLButtonTarget(1);
590         if(pointSBT0&&pointSBT1) {
591           if (pointSBT0->GetX()>pointSBT1->GetX()){
592             pointSBT0=myDetElem->GetLButtonTarget(1);
593             pointSBT1=myDetElem->GetLButtonTarget(0);
594           }
595           Double_t lX = pointSBT0->GetX();
596           while(lX<pointSBT1->GetX()){
597             Double_t lY = pointSBT0->GetY()-200;
598             while(lY<pointSBT0->GetY()+200){
599               pl[0] = lX;  pl[1] = lY;  pl[2] = 0.;
600               (TGeoCombiTrans((*(myDetElem->GetBaseTrf()))*(*(myDetElem->GetAlignTrf())))).LocalToMaster(pl,pg);
601               if(myDetElem->GetGButtonTarget(0)->GetPointName().Contains("A")){
602                 hSSTap->Fill(pg[0],pg[1],-pg[2]);
603               }
604               else {
605                 hSSTcp->Fill(pg[0],pg[1],-pg[2]);
606               }
607               lY+=hSSTap->GetYaxis()->GetBinWidth(1);
608             }
609             lX+=hSSTap->GetXaxis()->GetBinWidth(1);
610           }
611         }
612       }
613     }
614
615     if( saveps ){
616       ps->NewPage();
617     }
618
619     // View from side A
620     str = Form("Chamber %d - side A",chId+1);
621     TCanvas *cvn0 = new TCanvas("cvn0",str,cWidth,cHeight);
622     canvas = cvn0;
623     canvas->Range(0,0,21,29);
624   
625     title = Form(" Deformations of chamber %d - View from side A ",chId+1);
626     TPaveLabel *theTitle0 = new TPaveLabel(3,27.0,18,28.5,title,"br");
627     theTitle = theTitle0;
628     theTitle->SetFillColor(18);
629     theTitle->SetTextFont(32);
630     theTitle->SetTextSize(0.4);
631     theTitle->SetTextColor(1);
632     theTitle->Draw();
633  
634     TPad *pad0 = new TPad("pad0","pad0",0.01,0.01,0.98,0.91,0);
635     pad = pad0;
636     pad->Draw();
637     pad->Divide(2,2);
638
639     pad->cd(1);
640     hCPSTa->SetMinimum(9500);
641     hCPSTa->SetMaximum(9750);
642     hCPSTa->Draw("lego2z");
643
644     pad->cd(2);
645     hSSTa->SetMinimum(6680);
646     hSSTa->SetMaximum(6700);
647     hSSTa->Draw("lego2z");
648
649     pad->cd(3);
650
651     pad->cd(4);
652     hSSTap->SetMinimum(9500);
653     hSSTap->SetMaximum(9800);
654     hSSTap->Draw("surf2z");
655
656     pad->Update();
657     if(saveps){
658       ps->NewPage();
659     }
660
661     // Inv Mass, Multiplicity
662     str = Form("Chamber %d - side C",chId+1);
663     TCanvas *cvn1 = new TCanvas("cvn1",str,cWidth,cHeight);
664     canvas = cvn1;
665     canvas->Range(0,0,21,29);
666
667     title = Form(" Deformations of chamber %d - View from side C ",chId+1);
668     TPaveLabel *theTitle1 = new TPaveLabel(3,27.0,18,28.5,title,"br");
669     theTitle = theTitle1;
670     theTitle->SetFillColor(18);
671     theTitle->SetTextFont(32);
672     theTitle->SetTextSize(0.4);
673     theTitle->SetTextColor(1);
674     theTitle->Draw();
675  
676     TPad *pad1 = new TPad("pad1","pad1",0.01,0.01,0.98,0.91,0);
677     pad = pad1;
678     pad->Draw();
679     pad->Divide(2,2);
680
681     pad->cd(1);
682     hCPSTc->SetMinimum(6790);
683     hCPSTc->SetMaximum(6810);
684     hCPSTc->Draw("lego2z");
685
686     pad->cd(2);
687     hSSTc->SetMinimum(6820);
688     hSSTc->SetMaximum(6840);
689     hSSTc->Draw("lego2z");
690
691     pad->cd(3);
692
693     pad->cd(4);
694     hSSTcp->SetMinimum(9600);
695     hSSTcp->SetMaximum(9850);
696     hSSTcp->Draw("surf2z");
697   }
698
699   pad->Update();
700   if( saveps ){
701     ps->Close();
702   }
703
704   TFile *hFile = new TFile(Form("spCH%d_PL2G_ISF.root",chId+1),"RECREATE");
705   gDeltaDiffCenXDetElem0->Write("gDeltaDiffCenXDetElem0");
706   gDeltaDiffCenYDetElem0->Write("gDeltaDiffCenYDetElem0");
707   gDeltaDiffCenZDetElem0->Write("gDeltaDiffCenZDetElem0");
708   gDeltaDiffPsiDetElem0->Write("gDeltaDiffPsiDetElem0");
709   gDeltaDiffThtDetElem0->Write("gDeltaDiffThtDetElem0");
710   gDeltaDiffPhiDetElem0->Write("gDeltaDiffPhiDetElem0");
711   myDetElemDeltaDiffCenX->Write();
712   myDetElemDeltaDiffCenY->Write();
713   myDetElemDeltaDiffCenZ->Write();
714   myDetElemDeltaDiffRotX->Write();
715   myDetElemDeltaDiffRotY->Write();
716   myDetElemDeltaDiffRotZ->Write();
717   hFile->Close();
718   hFile->Delete();
719 }