]> git.uio.no Git - u/mrichter/AliRoot.git/blame - MUON/MUONSurveyCh5.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coding rule violations fixed.
[u/mrichter/AliRoot.git] / MUON / MUONSurveyCh5.C
CommitLineData
ba8b0266 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 5
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
86void MUONSurveyCh5() {
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 = 4;
99 Int_t nChs = 2;
100 Int_t nDetElems = 18;
101 Int_t nDetElemsI = 9;
102 // Int_t nDetElemsO = 9;
103 Int_t iDetElemToDetElemId[18] = {514,506,516,508,500,510,502,512,504,513,503,511,501,509,517,507,515,505};
104 Int_t iDetElemPseudoIdToDetElem[18] = {4,12,6,10,8,17,1,15,3,13,5,11,7,9,0,16,2,14};
105 Int_t iDetElemsOfChamber[2][9] = {{0,16,2,14,4,12,6,10,8},
106 {9,7,11,5,13,3,15,1,17}};
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/AliceSt3_TC5_1381b.txt");
119 myChamberO = (AliMUONSurveyChamber*)myChamberArray->At(1);
120 myChamberO->GetSurveyObj()->FillFromLocalFile("../Reports/AliceSt3_TC5_1381b.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(Form("%dC_1000",chId+1),4);
131 myChamberI->AddStickerTargets(Form("%dA_702",chId+1),9);
132 myChamberI->AddStickerTargets(Form("%dA_703",chId+1),9);
133 myChamberI->AddStickerTargets(Form("%dC_700",chId+1),9);
134 myChamberI->AddStickerTargets(Form("%dC_701",chId+1),9);
135 myChamberI->AddLButtonTargets(lSO->GetData(),Form("%d_1000",chId+1),4);
136 myChamberO->AddGButtonTargets(Form("%dC_1010",chId+1),4);
137 myChamberO->AddStickerTargets(Form("%dA_700",chId+1),9);
138 myChamberO->AddStickerTargets(Form("%dA_701",chId+1),9);
139 myChamberO->AddStickerTargets(Form("%dC_702",chId+1),9);
140 myChamberO->AddStickerTargets(Form("%dC_703",chId+1),9);
141 myChamberO->AddLButtonTargets(lSO->GetData(),Form("%d_1010",chId+1),4);
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 if(iDetElem<9) baseName = Form("%dA_50%d",chId+1,iDetElem+1);
153 else baseName = Form("%dC_50%d",chId+1,iDetElem+1-9);
154 myDetElem->AddStickerTargets(baseName,9);
155 if(iDetElem<9) baseName = Form("%dA_60%d",chId+1,iDetElem+1);
156 else baseName = Form("%dC_60%d",chId+1,iDetElem+1-9);
157 myDetElem->AddGButtonTargets(baseName,2);
158 if(iDetElemToDetElemId[iDetElem]-500<10) baseName = Form("%d_50%d",chId+1,iDetElemToDetElemId[iDetElem]-500);
159 else baseName = Form("%d_5%d",chId+1,iDetElemToDetElemId[iDetElem]-500);
160 myDetElem->AddLButtonTargets(lSO->GetData(),baseName,2);
161 }
162 }
163 }
164 printf("All targets added! \n");
165
166 Double_t **lCenDetElem = new Double_t*[nDetElems+nChs];
167 Double_t **lRotDetElem = new Double_t*[nDetElems+nChs];
168 Double_t **lDiffCenDetElem0 = new Double_t*[nDetElems+nChs];
169 Double_t **lDiffRotDetElem0 = new Double_t*[nDetElems+nChs];
170 Double_t **lDiffThCenDetElem0 = new Double_t*[nDetElems+nChs];
171 Double_t **lDiffThRotDetElem0 = new Double_t*[nDetElems+nChs];
172 Double_t **lDeltaDiffCenDetElem0 = new Double_t*[nDetElems+nChs];
173 Double_t **lDeltaDiffRotDetElem0 = new Double_t*[nDetElems+nChs];
174
175 for (int iDetElem=0; iDetElem<nDetElems+nChs; iDetElem++) {
176 lCenDetElem[iDetElem] = new Double_t[3];
177 lRotDetElem[iDetElem] = new Double_t[3];
178 lDiffCenDetElem0[iDetElem] = new Double_t[3];
179 lDiffRotDetElem0[iDetElem] = new Double_t[3];
180 lDiffThCenDetElem0[iDetElem] = new Double_t[3];
181 lDiffThRotDetElem0[iDetElem] = new Double_t[3];
182 lDeltaDiffCenDetElem0[iDetElem] = new Double_t[3];
183 lDeltaDiffRotDetElem0[iDetElem] = new Double_t[3];
184 }
185
186 TGeoCombiTrans dtrfDetElem[nDetElems+nChs];
187 TGeoCombiTrans localTrfDetElem[nDetElems+nChs];
188 TGeoCombiTrans localTrfThDetElem[nDetElems+nChs];
189
190 // Import TGeo geometry
191 char* geoFilename = "geometry.root";
192 if ( ! AliGeomManager::GetGeometry() ) {
193 AliGeomManager::LoadGeometry(geoFilename);
194 if (! AliGeomManager::GetGeometry() ) {
195 printf("MUONSurveyCh%d: getting geometry from file %s failed\n", chId+1,geoFilename);
196 return;
197 }
198 cout << "geometry imported" << endl;
199 }
200
201 AliMUONGeometryTransformer *transform = new AliMUONGeometryTransformer();
202 transform->LoadGeometryData();
203
204 TGeoCombiTrans trfThChamber;
205 TGeoCombiTrans trfThDetElem;
206
207 for (int iCh =0; iCh<=1; iCh++) {
208 myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
209
210 trfThChamber = TGeoCombiTrans(*transform->GetModuleTransformerByDEId(iDetElemToDetElemId[iCh*nDetElemsI])->GetTransformation());
211 trfThChamber.SetTranslation(trfThChamber.GetTranslation()[0]*10,
212 trfThChamber.GetTranslation()[1]*10,
213 trfThChamber.GetTranslation()[2]*10);
214 trfThChamber.Print();
215 myChamber->SetLocalTransformation(new TGeoCombiTrans(trfThChamber),kTRUE);
216
217 // Set Chamber plane function
218 cout << "Setting plane for Chamber" << iCh+1 << " ..." << endl;
219 myChamber->SetPlane(Form("fChamber%d",iCh+1));
220 myChamber->SetPlaneParameters(0.,0.,0.);
221
222 // Fit a plane to sticker targets
223 Double_t lCChi2 = myChamber->FitPlane();
224 cout << "... chi2 = " << lCChi2 << " ..." << endl;
225
226 // Get best transformation from fitting method
227 // (local to global transformation)
228 cout << "Trying fitting method for chamber " << iCh << endl;
229 myChamber->SurveyToAlign();
230 // myChamber->SurveyToAlign(myChamber->GetPlane()->GetParameter(0),myChamber->GetPlane()->GetParameter(1),myChamber->GetPlane()->GetParError(0),myChamber->GetPlane()->GetParError(1));
231
232 for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
233 myDetElem = myChamber->GetDetElem(iDetElemI);
234 Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
235
236 trfThDetElem.Clear();
237 trfThDetElem = TGeoCombiTrans(*transform->GetDetElement(iDetElemToDetElemId[iDetElem])->GetLocalTransformation());
238 trfThDetElem.SetTranslation(trfThDetElem.GetTranslation()[0]*10,
239 trfThDetElem.GetTranslation()[1]*10,
240 trfThDetElem.GetTranslation()[2]*10);
241 trfThDetElem.Print();
242
243 myDetElem->SetLocalTransformation(new TGeoCombiTrans(trfThDetElem),kTRUE);
244
245 for (int iCor=0; iCor<3; iCor++){
246 lCenDetElem[iDetElem][iCor] = 0;
247 lRotDetElem[iDetElem][iCor] = 0;
248 }
249
250 if (bMonitor){
251 // MONITOR: Draw graph with sticker targets for plane fitting
252 myDetElem->DrawSTargets();
253 gPad->Update();
254 }
255
256 // Get DetElem transformation.
257 // Psi and Theta are obtained by fitting a plane to the sticker targets.
258 // Then Xc, Yc, Zc and Phi are obtained by solving the equations to the ref.
259 // syst. transformation of the button targets
260
261 // Set DetElem plane function
262 cout << "Setting plane for DetElem" << iDetElem+1 << " ..." << endl;
263 myDetElem->SetPlane(Form("fDetElem%d",iDetElem+1));
264 myDetElem->SetPlaneParameters(0.,0.,30.)
265;
266 // Fit a plane to sticker targets
267 Double_t lChi2 = myDetElem->FitPlane();
268 cout << "... chi2 = " << lChi2 << " ..." << endl;
269
270 lRotDetElem[iDetElem][0]=TMath::ATan(myDetElem->GetPlane()->GetParameter(0));
271 lRotDetElem[iDetElem][1]=TMath::ATan(myDetElem->GetPlane()->GetParameter(1));
272
273 // Calculate Mean transformation using previous plane fit
274 // and pairs of button targets
275 myDetElem->CalculateMeanTransf(lCenDetElem[iDetElem],lRotDetElem[iDetElem]);
276
277 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;
278
279
280 // Get best transformation from fitting method
281 // (local to global transformation)
282 cout << "Trying fitting method for DetElemId " << iDetElemToDetElemId[iDetElem] << endl;
283 // myDetElem->SurveyToAlign();
284 myDetElem->SurveyToAlign(lRotDetElem[iDetElem][0],lRotDetElem[iDetElem][1],myDetElem->GetPlane()->GetParError(0),myDetElem->GetPlane()->GetParError(1));
285 }
286 }
287
288 // Print found transformation of Detection Element (plane fit + loop)
289 for (int iCh =0; iCh<=1; iCh++) {
290 for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
291 Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
292 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;
293 }
294 }
295
296 // Print Theoretical transformation of Detection Element
297 for (int iCh =0; iCh<=1; iCh++) {
298 myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
299 for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
300 myChamber->GetDetElem(iDetElemI)->PrintLocalTrf();
301 }
302 }
303
304 // Print found delta transformation of Detection Element
305 for (int iCh =0; iCh<=1; iCh++) {
306 myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
307 for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
308 myChamber->GetDetElem(iDetElemI)->PrintAlignTrf();
309 }
310 }
311
312 //
313 // Compare transformations to expected ones
314 //
315 Int_t iDetElemToPos[18] = {0, 16, 2, 14, 4, 12, 6, 10, 8, 9, 7, 11, 5, 13, 3, 15, 1, 17};
316
317 TGraphErrors *gDeltaDiffCenXDetElem0 = new TGraphErrors(nDetElems);
318 TGraphErrors *gDeltaDiffCenYDetElem0 = new TGraphErrors(nDetElems);
319 TGraphErrors *gDeltaDiffCenZDetElem0 = new TGraphErrors(nDetElems);
320 TGraphErrors *gDeltaDiffPsiDetElem0 = new TGraphErrors(nDetElems);
321 TGraphErrors *gDeltaDiffThtDetElem0 = new TGraphErrors(nDetElems);
322 TGraphErrors *gDeltaDiffPhiDetElem0 = new TGraphErrors(nDetElems);
323
324 for (int iCh =0; iCh<=1; iCh++) {
325 myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
326 for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
327 myChamber->GetDetElem(iDetElemI)->GetAlignTrf()->Print();
328 }
329 }
330
331 for (int iCh =0; iCh<=1; iCh++) {
332 myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
333 // Store delta transformations to use for CDB entry creation
334 dtrfDetElem[nDetElems+iCh].Clear();
335 dtrfDetElem[nDetElems+iCh] = *(myChamber->GetAlignTrf());
336
337 // Those are for checks and visualizations
338 localTrfDetElem[nDetElems+iCh].Clear();
339 localTrfDetElem[nDetElems+iCh] = (*(myChamber->GetLocalTrf())*(*(myChamber->GetAlignTrf())));
340 localTrfDetElem[nDetElems+iCh].Print();
341 lDiffCenDetElem0[nDetElems+iCh] = (Double_t*)localTrfDetElem[nDetElems+iCh].GetTranslation();
342 AliMUONSurveyUtil::MatrixToAngles(localTrfDetElem[nDetElems+iCh].GetRotationMatrix(),lDiffRotDetElem0[nDetElems+iCh]);
343
344 localTrfThDetElem[nDetElems+iCh].Clear();
345 localTrfThDetElem[nDetElems+iCh] = (*(myChamber->GetLocalTrf()));
346 localTrfThDetElem[nDetElems+iCh].Print();
347 lDiffThCenDetElem0[nDetElems+iCh] = (Double_t*)localTrfThDetElem[nDetElems+iCh].GetTranslation();
348 AliMUONSurveyUtil::MatrixToAngles(localTrfThDetElem[nDetElems+iCh].GetRotationMatrix(),lDiffThRotDetElem0[nDetElems+iCh]);
349
350 for (int iCor=0; iCor<3; iCor++){
351 lDeltaDiffCenDetElem0[nDetElems+iCh][iCor] = lDiffCenDetElem0[nDetElems+iCh][iCor]-lDiffThCenDetElem0[nDetElems+iCh][iCor];
352 lDeltaDiffRotDetElem0[nDetElems+iCh][iCor] = lDiffRotDetElem0[nDetElems+iCh][iCor]-lDiffThRotDetElem0[nDetElems+iCh][iCor];
353 if (lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]>TMath::Pi()) lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]-=TMath::TwoPi();
354 if (lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]<-TMath::Pi()) lDeltaDiffRotDetElem0[nDetElems+iCh][iCor]+=TMath::TwoPi();
355 }
356
357 for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
358 Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
359 myDetElem = myChamber->GetDetElem(iDetElemI);
360 // Store delta transformations to use for CDB entry creation
361 dtrfDetElem[iDetElem].Clear();
362 dtrfDetElem[iDetElem] = *(myDetElem->GetAlignTrf());
363
364 // Those are for checks and visualizations
365 localTrfDetElem[iDetElem].Clear();
366 localTrfDetElem[iDetElem] = (*(myDetElem->GetLocalTrf())*(*(myDetElem->GetAlignTrf())));
367 // localTrfDetElem[iDetElem] = (*(myDetElem->GetBaseTrf())*(*(myDetElem->GetAlignTrf())));
368 localTrfDetElem[iDetElem].Print();
369 lDiffCenDetElem0[iDetElem] = (Double_t*)localTrfDetElem[iDetElem].GetTranslation();
370 AliMUONSurveyUtil::MatrixToAngles(localTrfDetElem[iDetElem].GetRotationMatrix(),lDiffRotDetElem0[iDetElem]);
371// lDiffCenDetElem0[iDetElem] = lCenDetElem[iDetElem];
372// lDiffRotDetElem0[iDetElem] = lRotDetElem[iDetElem];
373
374 localTrfThDetElem[iDetElem].Clear();
375 localTrfThDetElem[iDetElem] = (*(myDetElem->GetLocalTrf()));
376 // localTrfThDetElem[iDetElem] = (*(myDetElem->GetBaseTrf()));
377 localTrfThDetElem[iDetElem].Print();
378 lDiffThCenDetElem0[iDetElem] = (Double_t*)localTrfThDetElem[iDetElem].GetTranslation();
379 AliMUONSurveyUtil::MatrixToAngles(localTrfThDetElem[iDetElem].GetRotationMatrix(),lDiffThRotDetElem0[iDetElem]);
380
381 for (int iCor=0; iCor<3; iCor++){
382 lDeltaDiffCenDetElem0[iDetElem][iCor] = lDiffCenDetElem0[iDetElem][iCor]-lDiffThCenDetElem0[iDetElem][iCor];
383 lDeltaDiffRotDetElem0[iDetElem][iCor] = lDiffRotDetElem0[iDetElem][iCor]-lDiffThRotDetElem0[iDetElem][iCor];
384 if (lDeltaDiffRotDetElem0[iDetElem][iCor]>TMath::Pi()) lDeltaDiffRotDetElem0[iDetElem][iCor]-=TMath::TwoPi();
385 if (lDeltaDiffRotDetElem0[iDetElem][iCor]<-TMath::Pi()) lDeltaDiffRotDetElem0[iDetElem][iCor]+=TMath::TwoPi();
386 }
387
388 gDeltaDiffCenXDetElem0->SetPoint(iDetElem,lDeltaDiffCenDetElem0[iDetElem][0],iDetElemToPos[iDetElem]+1);
389 gDeltaDiffCenYDetElem0->SetPoint(iDetElem,lDeltaDiffCenDetElem0[iDetElem][1],iDetElemToPos[iDetElem]+1);
390 gDeltaDiffCenZDetElem0->SetPoint(iDetElem,lDeltaDiffCenDetElem0[iDetElem][2],iDetElemToPos[iDetElem]+1);
391 gDeltaDiffPsiDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][0],iDetElemToPos[iDetElem]+1);
392 gDeltaDiffThtDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][1],iDetElemToPos[iDetElem]+1);
393 gDeltaDiffPhiDetElem0->SetPoint(iDetElem,1e3*lDeltaDiffRotDetElem0[iDetElem][2],iDetElemToPos[iDetElem]+1);
394 gDeltaDiffCenXDetElem0->SetPointError(iDetElem,myDetElem->GetFitter()->GetParError(0),0.);
395 gDeltaDiffCenYDetElem0->SetPointError(iDetElem,myDetElem->GetFitter()->GetParError(1),0.);
396 gDeltaDiffCenZDetElem0->SetPointError(iDetElem,myDetElem->GetFitter()->GetParError(2),0.);
397 gDeltaDiffPsiDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(3),0.);
398 gDeltaDiffThtDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(4),0.);
399 gDeltaDiffPhiDetElem0->SetPointError(iDetElem,1e3*myDetElem->GetFitter()->GetParError(5),0.);
400 }
401 }
402
403 // Apply the found misalignments to the geometry
404 AliMUONGeometryTransformer *newTransform = AliMUONSurveyUtil::ReAlign(transform,chId,nDetElems,iDetElemPseudoIdToDetElem,dtrfDetElem,true);
405 newTransform->WriteTransformations(Form("transform2ReAlignSurveyCh%d.dat",chId+1));
406
407 if(bOCDB){
408 // Generate realigned data in local cdb
409 const TClonesArray* array = newTransform->GetMisAlignmentData();
410
411 // Set the alignment resolution in the align objects for this chamber
412 Double_t chResX = 0.1*myChamberI->GetAlignResX();
413 Double_t chResY = 0.1*myChamberI->GetAlignResY();
414 Double_t deResX = 0.1*(myChamberI->GetMeanDetElemAlignResX()+myChamberO->GetMeanDetElemAlignResX())/2.;
415 Double_t deResY = 0.1*(myChamberI->GetMeanDetElemAlignResY()+myChamberO->GetMeanDetElemAlignResY())/2.;
416 printf("Chamber alignment resolutions: resX=%f , resY=%f\n",chResX,chResY);
417 printf("Detection Elements alignment resolutions: resX=%f , resY=%f\n",deResX,deResY);
418 chResX = TMath::Sqrt(0.1*0.1+chResX*chResX);
419 chResY = TMath::Sqrt(0.1*0.1+chResY*chResY);
420 AliMUONSurveyUtil::SetAlignmentResolution(array,chId,chResX,chResY,deResX,deResY);
421
422 // CDB manager
423 AliCDBManager* cdbManager = AliCDBManager::Instance();
424 cdbManager->SetDefaultStorage(Form("local://ReAlignSurveyCh%dCDB",chId+1));
425
426 AliCDBMetaData* cdbData = new AliCDBMetaData();
427 cdbData->SetResponsible("Dimuon Offline project");
428 cdbData->SetComment("MUON alignment objects with survey realignment");
429 AliCDBId id("MUON/Align/Data", 0, AliCDBRunRange::Infinity());
430 cdbManager->Put(const_cast<TClonesArray*>(array), id, cdbData);
431 }
432
433 for(Int_t iCor=0; iCor<3; iCor++){
434 for(Int_t iDetElem=0; iDetElem<nDetElems; iDetElem++){
435 cout << lDeltaDiffCenDetElem0[iDetElem][iCor] << " " << lDiffCenDetElem0[iDetElem][iCor] << " " << lDiffThCenDetElem0[iDetElem][iCor] << endl;
436 }
437 cout << endl;
438 }
439 cout << endl;
440 for(Int_t iCor=0; iCor<3; iCor++){
441 for(Int_t iDetElem=0; iDetElem<nDetElems; iDetElem++){
442 cout << lDeltaDiffRotDetElem0[iDetElem][iCor] << " " << lDiffRotDetElem0[iDetElem][iCor] << " " << lDiffThRotDetElem0[iDetElem][iCor] << endl;
443 }
444 cout << endl;
445 }
446
447 TH1F *myDetElemDeltaDiffCenX = new TH1F("myDetElemDeltaDiffCenX","myDetElemDeltaDiffCenX",100,-10,10);
448 myDetElemDeltaDiffCenX->SetMaximum(nDetElems+1);
449 myDetElemDeltaDiffCenX->SetMinimum(0);
450 TH1F *myDetElemDeltaDiffCenY = new TH1F("myDetElemDeltaDiffCenY","myDetElemDeltaDiffCenY",100,-10,10);
451 myDetElemDeltaDiffCenY->SetMaximum(nDetElems+1);
452 myDetElemDeltaDiffCenY->SetMinimum(0);
453 TH1F *myDetElemDeltaDiffCenZ = new TH1F("myDetElemDeltaDiffCenZ","myDetElemDeltaDiffCenZ",100,-40,40);
454 myDetElemDeltaDiffCenZ->SetMaximum(nDetElems+1);
455 myDetElemDeltaDiffCenZ->SetMinimum(0);
456
457 TH1F *myDetElemDeltaDiffRotX = new TH1F("myDetElemDeltaDiffRotX","myDetElemDeltaDiffRotX",100,-40,40);
458 myDetElemDeltaDiffRotX->SetMaximum(nDetElems+1);
459 myDetElemDeltaDiffRotX->SetMinimum(0);
460 TH1F *myDetElemDeltaDiffRotY = new TH1F("myDetElemDeltaDiffRotY","myDetElemDeltaDiffRotY",100,-40,40);
461 myDetElemDeltaDiffRotY->SetMaximum(nDetElems+1);
462 myDetElemDeltaDiffRotY->SetMinimum(0);
463 TH1F *myDetElemDeltaDiffRotZ = new TH1F("myDetElemDeltaDiffRotZ","myDetElemDeltaDiffRotZ",100,-10,10);
464 myDetElemDeltaDiffRotZ->SetMaximum(nDetElems+1);
465 myDetElemDeltaDiffRotZ->SetMinimum(0);
466
467 //
468 // ******** Starting plots
469 //
470 TCanvas *canvas;
471 TPad *pad;
472 TPaveLabel *theTitle;
473 gStyle->SetPalette(1);
474
475 TPostScript *ps = 0;
476
477 if( saveps ){
478 ps = new TPostScript(Form("surveyChamber%d",chId+1),filetype);
479 ps->NewPage();
480 }
481
482 // Observed misalignments
483 str = Form("Chamber %d",chId+1);
484 TCanvas *cvn2 = new TCanvas("cvn2",str,cWidth,cHeight);
485 canvas = cvn2;
486 canvas->Range(0,0,21,29);
487
488 title = Form(" MisAlignments Chamber %d - PL2G - In Frame ",chId+1);
489 TPaveLabel *theTitle2 = new TPaveLabel(3,27.0,18,28.5,title,"br");
490 theTitle = theTitle2;
491 theTitle->SetFillColor(18);
492 theTitle->SetTextFont(32);
493 theTitle->SetTextSize(0.4);
494 theTitle->SetTextColor(1);
495 theTitle->Draw();
496
497 TPad *pad2 = new TPad("pad2","pad2",0.01,0.01,0.98,0.91,0);
498 pad = pad2;
499 pad->Draw();
500 TLine *ch0Line = new TLine(0,1,0,2);
501 TLine *ch1Line = new TLine(0,1,0,2);
502 ch1Line->SetLineStyle(2);
503 pad->Divide(3,2);
504
505 pad->cd(1);
506 myDetElemDeltaDiffCenX->Draw();
507 myDetElemDeltaDiffCenX->SetXTitle("#Delta[xc_{i}^{m}-xc_{i}^{th}] (mm)");
508 myDetElemDeltaDiffCenX->SetYTitle("DetElem arbitrary ordering");
509 gDeltaDiffCenXDetElem0->SetMarkerStyle(20);
510 gDeltaDiffCenXDetElem0->Draw("P");
511 ch0Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+0][0],0.5,lDeltaDiffCenDetElem0[nDetElems+0][0],9.5);
512 ch1Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+1][0],9.5,lDeltaDiffCenDetElem0[nDetElems+1][0],18.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],9.5);
521 ch1Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+1][1],9.5,lDeltaDiffCenDetElem0[nDetElems+1][1],18.5);
522
523 pad->cd(3);
524 myDetElemDeltaDiffCenZ->Draw();
525 myDetElemDeltaDiffCenZ->SetXTitle("#Delta[zc_{i}^{m}-zc_{i}^{th}] (mm)");
526 myDetElemDeltaDiffCenZ->SetYTitle("DetElem arbitrary ordering");
527 gDeltaDiffCenZDetElem0->SetMarkerStyle(20);
528 gDeltaDiffCenZDetElem0->Draw("P");
529 ch0Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+0][2],0.5,lDeltaDiffCenDetElem0[nDetElems+0][2],9.5);
530 ch1Line->DrawLine(lDeltaDiffCenDetElem0[nDetElems+1][2],9.5,lDeltaDiffCenDetElem0[nDetElems+1][2],18.5);
531
532 pad->cd(4);
533 myDetElemDeltaDiffRotX->Draw();
534 myDetElemDeltaDiffRotX->SetXTitle("#Delta[#psi_{i}^{m}-#psi_{i}^{th}] (mrad)");
535 myDetElemDeltaDiffRotX->SetYTitle("DetElem arbitrary ordering");
536 gDeltaDiffPsiDetElem0->SetMarkerStyle(20);
537 gDeltaDiffPsiDetElem0->Draw("P");
538 ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][0],9.5);
539 ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][0],9.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][0],18.5);
540
541 pad->cd(5);
542 myDetElemDeltaDiffRotY->Draw();
543 myDetElemDeltaDiffRotY->SetXTitle("#Delta[#theta_{i}^{m}-#theta_{i}^{th}] (mrad)");
544 myDetElemDeltaDiffRotY->SetYTitle("DetElem arbitrary ordering");
545 gDeltaDiffThtDetElem0->SetMarkerStyle(20);
546 gDeltaDiffThtDetElem0->Draw("P");
547 ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][1],9.5);
548 ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][1],9.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][1],18.5);
549
550 pad->cd(6);
551 myDetElemDeltaDiffRotZ->Draw();
552 myDetElemDeltaDiffRotZ->SetXTitle("#Delta[#phi_{i}^{m}-#phi_{i}^{th}] (mrad)");
553 myDetElemDeltaDiffRotZ->SetYTitle("DetElem arbitrary ordering");
554 gDeltaDiffPhiDetElem0->SetMarkerStyle(20);
555 gDeltaDiffPhiDetElem0->Draw("P");
556 ch0Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],0.5,1e3*lDeltaDiffRotDetElem0[nDetElems+0][2],9.5);
557 ch1Line->DrawLine(1e3*lDeltaDiffRotDetElem0[nDetElems+1][2],9.5,1e3*lDeltaDiffRotDetElem0[nDetElems+1][2],18.5);
558
559 pad->Update();
560
561 if(bMonitor){
562 // MONITOR: Histograms for monitoring
563 TH2F *hCPSTa = new TH2F("hCPSTa","hCPSTa",100,-2000,2000,100,-2000,2000);
564 TH2F *hCPSTc = new TH2F("hCPSTc","hCPSTc",100,-2000,2000,100,-2000,2000);
565 TH2F *hSSTa = new TH2F("hSSTa","hSSTa",100,-2000,2000,100,-2000,2000);
566 TH2F *hSSTc = new TH2F("hSSTc","hSSTc",100,-2000,2000,100,-2000,2000);
567 TH2F *hSSTap = new TH2F("hSSTap","hSSTap",800,-2000,2000,800,-2000,2000);
568 TH2F *hSSTcp = new TH2F("hSSTcp","hSSTcp",800,-2000,2000,800,-2000,2000);
569
570 // MONITOR: Fill histograms with chambers and slats sticker target positions
571 for (int iCh =0; iCh<=1; iCh++) {
572 myChamber = (AliMUONSurveyChamber*)myChamberArray->At(iCh);
573 cout << "Filling histograms of sticker target for chamber" << iCh+1 << " ..." << endl;
574 myChamber->FillCPSTHistograms(TString("C"),hCPSTc,TString("A"),hCPSTa);
575 myChamber->FillDESTHistograms(TString("C"),hSSTc,TString("A"),hSSTa);
576
577 for (int iDetElemI=0; iDetElemI<nDetElemsI; iDetElemI++){
578 // Int_t iDetElem = iDetElemsOfChamber[iCh][iDetElemI];
579 myDetElem = myChamber->GetDetElem(iDetElemI);
580 // MONITOR: Fill slat plane for fit monitor.
581 Double_t pl[3] = {0};
582 Double_t pg[3] = {0};
583 AliSurveyPoint *pointSBT0 = myDetElem->GetLButtonTarget(0);
584 AliSurveyPoint *pointSBT1 = myDetElem->GetLButtonTarget(1);
585 if(pointSBT0&&pointSBT1) {
586 if (pointSBT0->GetX()>pointSBT1->GetX()){
587 pointSBT0=myDetElem->GetLButtonTarget(1);
588 pointSBT1=myDetElem->GetLButtonTarget(0);
589 }
590 Double_t lX = pointSBT0->GetX();
591 while(lX<pointSBT1->GetX()){
592 Double_t lY = pointSBT0->GetY()-200;
593 while(lY<pointSBT0->GetY()+200){
594 pl[0] = lX; pl[1] = lY; pl[2] = 0.;
595 (TGeoCombiTrans((*(myDetElem->GetBaseTrf()))*(*(myDetElem->GetAlignTrf())))).LocalToMaster(pl,pg);
596 if(myDetElem->GetGButtonTarget(0)->GetPointName().Contains("A")){
597 hSSTap->Fill(pg[0],pg[1],-pg[2]);
598 }
599 else {
600 hSSTcp->Fill(pg[0],pg[1],-pg[2]);
601 }
602 lY+=hSSTap->GetYaxis()->GetBinWidth(1);
603 }
604 lX+=hSSTap->GetXaxis()->GetBinWidth(1);
605 }
606 }
607 }
608 }
609
610 if( saveps ){
611 ps->NewPage();
612 }
613
614 // View from side A
615 str = Form("Chamber %d - side A",chId+1);
616 TCanvas *cvn0 = new TCanvas("cvn0",str,cWidth,cHeight);
617 canvas = cvn0;
618 canvas->Range(0,0,21,29);
619
620 title = Form(" Deformations of chamber %d - View from side A ",chId+1);
621 TPaveLabel *theTitle0 = new TPaveLabel(3,27.0,18,28.5,title,"br");
622 theTitle = theTitle0;
623 theTitle->SetFillColor(18);
624 theTitle->SetTextFont(32);
625 theTitle->SetTextSize(0.4);
626 theTitle->SetTextColor(1);
627 theTitle->Draw();
628
629 TPad *pad0 = new TPad("pad0","pad0",0.01,0.01,0.98,0.91,0);
630 pad = pad0;
631 pad->Draw();
632 pad->Divide(2,2);
633
634 pad->cd(1);
635 hCPSTa->SetMinimum(9500);
636 hCPSTa->SetMaximum(9750);
637 hCPSTa->Draw("lego2z");
638
639 pad->cd(2);
640 hSSTa->SetMinimum(9500);
641 hSSTa->SetMaximum(9800);
642 hSSTa->Draw("lego2z");
643
644 pad->cd(3);
645
646 pad->cd(4);
647 hSSTap->SetMinimum(9500);
648 hSSTap->SetMaximum(9800);
649 hSSTap->Draw("surf2z");
650
651 pad->Update();
652 if(saveps){
653 ps->NewPage();
654 }
655
656 // Inv Mass, Multiplicity
657 str = Form("Chamber %d - side C",chId+1);
658 TCanvas *cvn1 = new TCanvas("cvn1",str,cWidth,cHeight);
659 canvas = cvn1;
660 canvas->Range(0,0,21,29);
661
662 title = Form(" Deformations of chamber %d - View from side C ",chId+1);
663 TPaveLabel *theTitle1 = new TPaveLabel(3,27.0,18,28.5,title,"br");
664 theTitle = theTitle1;
665 theTitle->SetFillColor(18);
666 theTitle->SetTextFont(32);
667 theTitle->SetTextSize(0.4);
668 theTitle->SetTextColor(1);
669 theTitle->Draw();
670
671 TPad *pad1 = new TPad("pad1","pad1",0.01,0.01,0.98,0.91,0);
672 pad = pad1;
673 pad->Draw();
674 pad->Divide(2,2);
675
676 pad->cd(1);
677 gStyle->SetPalette(1);
678 hCPSTc->SetMinimum(9550);
679 hCPSTc->SetMaximum(9850);
680 hCPSTc->Draw("lego2z");
681
682 pad->cd(2);
683 gStyle->SetPalette(1);
684 hSSTc->SetMinimum(9600);
685 hSSTc->SetMaximum(9850);
686 hSSTc->Draw("lego2z");
687
688 pad->cd(3);
689 gStyle->SetPalette(1);
690
691 pad->cd(4);
692 gStyle->SetPalette(1);
693 hSSTcp->SetMinimum(9600);
694 hSSTcp->SetMaximum(9850);
695 hSSTcp->Draw("surf2z");
696 }
697
698 pad->Update();
699 if( saveps ){
700 ps->Close();
701 }
702
703 TFile *hFile = new TFile(Form("spCH%d_PL2G_IF.root",chId+1),"RECREATE");
704 gDeltaDiffCenXDetElem0->Write("gDeltaDiffCenXDetElem0");
705 gDeltaDiffCenYDetElem0->Write("gDeltaDiffCenYDetElem0");
706 gDeltaDiffCenZDetElem0->Write("gDeltaDiffCenZDetElem0");
707 gDeltaDiffPsiDetElem0->Write("gDeltaDiffPsiDetElem0");
708 gDeltaDiffThtDetElem0->Write("gDeltaDiffThtDetElem0");
709 gDeltaDiffPhiDetElem0->Write("gDeltaDiffPhiDetElem0");
710 myDetElemDeltaDiffCenX->Write();
711 myDetElemDeltaDiffCenY->Write();
712 myDetElemDeltaDiffCenZ->Write();
713 myDetElemDeltaDiffRotX->Write();
714 myDetElemDeltaDiffRotY->Write();
715 myDetElemDeltaDiffRotZ->Write();
716 hFile->Close();
717 hFile->Delete();
718}