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