updates to read survey parameters from EDMS info
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALSurvey.cxx
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 // Objects of this class read txt file with survey data
19 // and convert the data into AliAlignObjParams of alignable EMCAL volumes.
20 // AliEMCALSurvey inherits TObject only to use AliLog "functions".
21 //
22 // Dummy functions originally written before EMCAL installation and
23 // survey are kept for backward compatibility, but now they are not
24 // used.
25 // Surveyed points on the face of each module are read in from file
26 // and converted to position of the center and roll, pitch, yaw angles
27 // of each surveyed SM.
28 //
29 // J.L. Klay - Cal Poly
30 // 07-Apr-2010
31 //
32
33 #include <fstream>
34
35 #include <TClonesArray.h>
36 #include <TGeoManager.h>
37 #include <TString.h>
38 #include <TMath.h>
39
40 #include "AliSurveyObj.h"
41 #include "AliSurveyPoint.h"
42
43 #include "AliAlignObjParams.h"
44 #include "AliEMCALGeometry.h"
45 #include "AliEMCALSurvey.h"
46 #include "AliLog.h"
47
48 ClassImp(AliEMCALSurvey)
49
50 //____________________________________________________________________________
51 AliEMCALSurvey::AliEMCALSurvey()
52                 : fNSuperModule(0),
53                   fSuperModuleData(0)
54 {
55   //Default constructor.
56 }
57
58 namespace {
59
60   //Coordinates for each SM described in survey reports
61
62   struct AliEMCALSuperModuleCoords {
63     Double_t fX1; //x coordinate of the center of supermodule
64     Double_t fY1; //y coordinate of the center of supermodule
65     Double_t fZ1; //z coordinate of the center of supermodule
66     Double_t fPhi;  //roll angle (phi) of supermodule
67     Double_t fTheta; //pitch (theta) of supermodule
68     Double_t fPsi;   //yaw (psi) of supermodule
69   };
70
71 }
72
73 //____________________________________________________________________________
74 AliEMCALSurvey::AliEMCALSurvey(const TString &txtFileName)
75                 : fNSuperModule(0),
76                   fSuperModuleData(0)
77 {
78   //Read survey data from txt file.
79   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
80   if (!geom) {
81     AliError("Cannot obtain AliEMCALGeometry instance.");
82     return;
83   }
84
85   AliSurveyObj *s1 = new AliSurveyObj();
86   s1->FillFromLocalFile(txtFileName);
87
88   TObjArray* points = s1->GetData();
89
90   fNSuperModule = geom->GetNumberOfSuperModules();
91
92   InitSuperModuleData(points);
93
94   //////////////////////////////////
95   //Old way with dummy survey file
96   //////////////////////////////////
97   //std::ifstream inputFile(txtFileName.Data());
98   //if (!inputFile) {
99   //  AliError(("Cannot open the survey file " + txtFileName).Data());
100   //  return;
101   //}
102   //
103   //Int_t dummyInt = 0;
104   //Double_t *xReal = new Double_t[fNSuperModule];
105   //Double_t *yReal = new Double_t[fNSuperModule];
106   //Double_t *zReal = new Double_t[fNSuperModule];
107   //
108   //for (Int_t i = 0; i < fNSuperModule; ++i) {
109   //  if (!inputFile) {
110   //    AliError("Error while reading input file.");
111   //    delete [] xReal;
112   //    delete [] yReal;
113   //    delete [] zReal;
114   //    return;
115   //  }
116   //  inputFile>>dummyInt>>xReal[i]>>yReal[i]>>zReal[i];
117   //}
118   //
119   //InitSuperModuleData(xReal, yReal, zReal);
120   //
121   //delete [] xReal;
122   //delete [] yReal;
123   //delete [] zReal;
124
125 }
126
127 //____________________________________________________________________________
128 AliEMCALSurvey::~AliEMCALSurvey()
129 {
130   //destructor
131   delete [] fSuperModuleData;
132 }
133
134 //____________________________________________________________________________
135 void AliEMCALSurvey::CreateAliAlignObjParams(TClonesArray &array)
136 {
137   //Create AliAlignObjParams.
138   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
139   if (!geom) {
140     AliError("Cannot obtain AliEMCALGeometry instance.");
141     return;
142   }
143
144   if (!gGeoManager) {
145     AliWarning("Cannot create local transformations for supermodules - gGeoManager does not exist.");
146     AliInfo("Null shifts and rotations will be created instead.");
147     return CreateNullObjects(array, geom);
148   }
149
150   Int_t arrayInd = array.GetEntries(), iIndex = 0;
151   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
152   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iIndex);
153
154   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
155     TString smodName(TString::Format("EMCAL/FullSupermodule%d", smodnum+1));
156     if(geom->GetKey110DEG() && smodnum >= 10) {
157       smodName = "EMCAL/HalfSupermodule";
158       smodName += (smodnum-10+1);
159     }    
160     AliEMCALSuperModuleDelta t(GetSuperModuleTransformation(smodnum));
161     new(array[arrayInd])
162       AliAlignObjParams(
163                         smodName.Data(), volid, 
164                         t.fXShift, t.fYShift, t.fZShift, 
165                         -t.fPsi, -t.fTheta, -t.fPhi, 
166                         false
167                         );
168     ++arrayInd;
169
170   }
171
172 }
173
174 //____________________________________________________________________________
175 void AliEMCALSurvey::CreateNullObjects(TClonesArray &array, const AliEMCALGeometry *geom)const
176 {
177   //Create null shifts and rotations.
178   Int_t arrayInd = array.GetEntries(), iIndex = 0;
179   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
180   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iIndex);
181
182   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
183     TString smodName(TString::Format("EMCAL/FullSupermodule%d", smodnum+1));
184     if(geom->GetKey110DEG() && smodnum >= 10) {
185       smodName = "EMCAL/HalfSupermodule";
186       smodName += (smodnum-10+1);
187     }
188     new(array[arrayInd]) AliAlignObjParams(smodName.Data(), volid, 0., 0., 0., 0., 0., 0., true);
189     ++arrayInd;
190   }
191 }
192
193 //____________________________________________________________________________
194 AliEMCALSurvey::AliEMCALSuperModuleDelta AliEMCALSurvey::GetSuperModuleTransformation(Int_t supModIndex)const
195 {
196   //Supermodule transformation.
197   AliEMCALSuperModuleDelta t = {0., 0., 0., 0., 0., 0.};
198   if (!fSuperModuleData)
199     return t;
200
201   return fSuperModuleData[supModIndex];
202 }
203
204 //____________________________________________________________________________
205 void AliEMCALSurvey::InitSuperModuleData(const TObjArray *svypts)
206 {
207   //This method uses the data points from the actual EMCAL survey to
208   //create the alignment matrices.  Only valid for measured(installed)
209   //SM, others will have null objects
210   
211   
212   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
213   //Center of supermodules
214   Float_t *pars = geom->GetSuperModulesPars();
215   Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
216   Double_t phi, phiRad, xpos, ypos, zpos;
217
218   AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
219   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
220     AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
221     phiRad = geom->GetPhiCenterOfSM(smodnum); //comes in radians
222     phi = phiRad*180./TMath::Pi(); //need degrees for AliAlignObjParams
223     xpos = rpos * TMath::Cos(phiRad);
224     ypos = rpos * TMath::Sin(phiRad);
225     zpos = pars[2];
226     if(geom->GetKey110DEG() && smodnum >= 10) {
227       xpos += (pars[1]/2. * TMath::Sin(phiRad));
228       ypos -= (pars[1]/2. * TMath::Cos(phiRad));
229     }
230     smc.fX1 = xpos;
231     smc.fY1 = ypos;
232     smc.fPhi = phi; //degrees
233     smc.fTheta = 0.; //degrees
234     smc.fPsi = 0.; //degrees
235     if(smodnum%2==0) {
236       smc.fZ1 = zpos;
237     } else {
238       smc.fZ1 = -zpos;
239     }
240
241     printf("PHI OF IDEAL SM = %.2f\n",smc.fPhi);
242
243   }
244
245   //Real coordinates of center and rotation angles need to be computed
246   //from survey points
247
248   char substr[100];
249   AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
250   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
251     AliEMCALSuperModuleCoords &smc = realSM[smodnum];
252     zpos = pars[2]; //center of SM in z from geometry
253
254     sprintf(substr,"4096%d",smodnum);
255     //retrieve components of four face points and determine average position of center
256     //in x,y,z
257
258     std::vector<Double_t> xval;
259     std::vector<Double_t> yval;
260     std::vector<Double_t> zval;
261     
262     for(Int_t i = 0; i < svypts->GetEntries(); i++) {
263       AliSurveyPoint* pt = (AliSurveyPoint*)svypts->At(i);
264       TString ptname = pt->GetPointName();
265       if(ptname.Contains(substr)) {
266         xval.push_back(pt->GetX()*100.); //convert m to cm
267         yval.push_back(pt->GetY()*100.); 
268         zval.push_back(pt->GetZ()*100.); 
269       }
270     }
271
272     //Take average of all relevant pairs of points
273     Double_t xReal = ((xval[1]+xval[7])/2.      //4X02 and 4X08
274                     + (xval[2]+xval[6])/2.      //4X03 and 4X07
275                     + (xval[3]+xval[5])/2.      //4X04 and 4X06
276                     + (xval[2]+xval[4])/2.)/4.; //4X03 and 4X05
277     Double_t yReal = ((yval[1]+yval[7])/2.      //4X02 and 4X08
278                     + (yval[2]+yval[6])/2.      //4X03 and 4X07
279                     + (yval[3]+yval[5])/2.      //4X04 and 4X06
280                     + (yval[2]+yval[4])/2.)/4.; //4X03 and 4X05
281     smc.fX1 = xReal;
282     smc.fY1 = yReal;
283
284     //Find average value of z for front face of SM
285     Double_t zReal = 0.;
286     Int_t nPoints = zval.size();
287     for(Int_t iz = 0; iz < nPoints; iz++) {
288       zReal += zval[iz];
289     }
290     if(nPoints > 0) zReal = zReal/nPoints;
291
292     if(smodnum%2==0) {
293       smc.fZ1 = zReal-zpos;  //z measured is along end,
294     } else {                 //convert to middle of SM
295       smc.fZ1 = zReal+zpos;
296     }
297
298     Double_t roll = ( TMath::ATan((yval[5]-yval[3])/(xval[5]-xval[3])) //4X04 and 4X06
299                       + TMath::ATan((yval[4]-yval[2])/(xval[4]-xval[2])) )/2.; //4X05 and 4X03
300     smc.fPhi = 90. + roll*TMath::RadToDeg();
301
302     //Note pitch calc only uses first 8 values, even though 10 are
303     //measured on the topmost modules
304     Double_t pitch = ( TMath::ATan((zval[0]-zval[1])/(yval[1]-yval[0])) //4X01 and 4X02
305                        + TMath::ATan((zval[2]-zval[3])/(yval[3]-yval[2])) //4X03 and 4X04
306                        + TMath::ATan((zval[4]-zval[5])/(yval[5]-yval[4])) //4X05 and 4X06
307                        + TMath::ATan((zval[6]-zval[7])/(yval[7]-yval[6])) )/4.; //4X07 and 4X08
308     smc.fTheta = 0. + pitch*TMath::RadToDeg();
309
310     Double_t yaw = ( TMath::ATan((zval[3]-zval[5])/(xval[5]-xval[3])) //4X04 and 4X06
311                      + TMath::ATan((zval[2]-zval[4])/(xval[4]-xval[2])) //4X03 and 4X05
312                      + TMath::ATan((zval[1]-zval[7])/(xval[7]-xval[1])) //4X02 and 4X08
313                      + TMath::ATan((zval[0]-zval[6])/(xval[6]-xval[0])) )/4.; //4X01 and 4X07
314     smc.fPsi = 0. + yaw*TMath::RadToDeg();
315
316   }//loop over supermodules
317   
318   fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
319   
320   for (Int_t i = 0; i < fNSuperModule; ++i) {
321     const AliEMCALSuperModuleCoords &real = realSM[i];
322     const AliEMCALSuperModuleCoords &ideal = idealSM[i];
323     AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
324     t.fXShift = real.fX1 - ideal.fX1;
325     t.fYShift = real.fY1 - ideal.fY1;
326     t.fZShift = real.fZ1 - ideal.fZ1;
327     t.fPhi = real.fPhi - ideal.fPhi;
328     t.fTheta = real.fTheta - ideal.fTheta;
329     t.fPsi = real.fPsi - ideal.fPsi;
330
331     printf("===================== SM %d =======================\n",i);
332     printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
333     printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
334     printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
335     printf("real theta (%.2f) - ideal theta (%.2f) = shift in theta %.2f\n",real.fTheta,ideal.fTheta,t.fTheta);
336     printf("real psi (%.2f) - ideal psi (%.2f) = shift in psi %.2f\n",real.fPsi,ideal.fPsi,t.fPsi);
337     printf("real phi (%.2f) - ideal phi (%.2f) = shift in phi %.2f\n",real.fPhi,ideal.fPhi,t.fPhi);
338     printf("===================================================\n");    
339   }
340
341   delete [] realSM;
342   delete [] idealSM;
343 }
344
345
346 //____________________________________________________________________________
347 void AliEMCALSurvey::InitSuperModuleData(const Double_t *xReal, const Double_t *yReal, const Double_t *zReal)
348 {
349   ///////////////////////////////////////
350   //Old dummy file way of doing it
351   //////////////////////////////////////
352   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
353   //Center of supermodules
354   Float_t *pars = geom->GetSuperModulesPars();
355   Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
356   Double_t phi, phiRad, xpos, ypos, zpos;
357
358   AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
359   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
360     AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
361     phiRad = geom->GetPhiCenterOfSM(smodnum); //comes in radians
362     phi = phiRad*180./TMath::Pi(); //need degrees for AliAlignObjParams
363     xpos = rpos * TMath::Cos(phiRad);
364     ypos = rpos * TMath::Sin(phiRad);
365     zpos = pars[2];
366     if(geom->GetKey110DEG() && smodnum >= 10) {
367       xpos += (pars[1]/2. * TMath::Sin(phiRad));
368       ypos -= (pars[1]/2. * TMath::Cos(phiRad));
369     }
370     smc.fX1 = xpos;
371     smc.fY1 = ypos;
372     if(smodnum%2==0) {
373       smc.fZ1 = zpos;
374     } else {
375       smc.fZ1 = -zpos;
376     }
377
378   }
379
380   AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
381   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
382     AliEMCALSuperModuleCoords &smc = realSM[smodnum];
383     zpos = pars[2];
384     smc.fX1 = xReal[smodnum];  //x and y match
385     smc.fY1 = yReal[smodnum];  //x and y match
386     if(smodnum%2==0) {
387       smc.fZ1 = zReal[smodnum]-zpos;  //z measured is along end,
388     } else {                          //convert to middle of SM
389       smc.fZ1 = zReal[smodnum]+zpos;
390     }
391   }
392   
393   fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
394   
395   for (Int_t i = 0; i < fNSuperModule; ++i) {
396     const AliEMCALSuperModuleCoords &real = realSM[i];
397     const AliEMCALSuperModuleCoords &ideal = idealSM[i];
398     AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
399     t.fTheta = TMath::ATan(real.fZ1  / real.fX1) - 
400                TMath::ATan(ideal.fZ1 / ideal.fX1);
401     t.fTheta *= TMath::RadToDeg();
402     t.fPsi = 0.;
403     t.fPhi = TMath::ATan(real.fY1 / real.fX1) - TMath::ATan(ideal.fY1 / ideal.fX1);
404     t.fPhi *= TMath::RadToDeg();
405     t.fXShift = real.fX1 - ideal.fX1;
406     t.fYShift = real.fY1 - ideal.fY1;
407     t.fZShift = real.fZ1 - ideal.fZ1;
408
409     printf("===================== SM %d =======================\n",i);
410     printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
411     printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
412     printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
413     printf("theta %.2f\n",t.fTheta);
414     printf("psi %.2f\n",t.fPsi);
415     printf("phi %.2f\n",t.fPhi);
416     printf("===================================================\n");    
417   }
418
419   delete [] realSM;
420   delete [] idealSM;
421 }
422