Updated version.
[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 //
26 // Surveyed points on the EMCAL support rails were used with the CATIA
27 // 3D graphics program to determine the positions of the bottom
28 // corners of the active area for each supermodule.  These numbers are
29 // read in from file and converted to position of the center and roll, 
30 // pitch, yaw angles of each installed SM.
31 //
32 // J.L. Klay - Cal Poly
33 // 21-May-2010
34 //
35
36 #include <fstream>
37
38 #include <TClonesArray.h>
39 #include <TGeoManager.h>
40 #include <TString.h>
41 #include <TMath.h>
42
43 #include "AliSurveyObj.h"
44 #include "AliSurveyPoint.h"
45
46 #include "AliAlignObjParams.h"
47 #include "AliEMCALGeometry.h"
48 #include "AliEMCALSurvey.h"
49 #include "AliLog.h"
50
51 ClassImp(AliEMCALSurvey)
52
53 //____________________________________________________________________________
54 AliEMCALSurvey::AliEMCALSurvey()
55  : fNSuperModule(0),
56    fSuperModuleData(0),
57    fDataType(kSurvey)
58 {
59   //Default constructor.
60 }
61
62 namespace {
63
64   //Coordinates for each SM described in survey reports
65
66   struct AliEMCALSuperModuleCoords {
67     Double_t fX1; //x coordinate of the center of supermodule
68     Double_t fY1; //y coordinate of the center of supermodule
69     Double_t fZ1; //z coordinate of the center of supermodule
70     Double_t fPsi;   //yaw (psi) of supermodule
71     Double_t fTheta; //pitch (theta) of supermodule
72     Double_t fPhi;  //roll angle (phi) of supermodule
73
74   };
75
76 }
77
78 //____________________________________________________________________________
79 AliEMCALSurvey::AliEMCALSurvey(const TString &txtFileName,const SurveyDataType_t type)
80   : fNSuperModule(0),
81     fSuperModuleData(0),
82     fDataType(type)
83 {
84   //Get the geometry object and then attempt to
85   //read survey data from a file, depending on which
86   //method (kSurvey or kDummy) is selected.
87
88   const AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
89   if (!geom) {
90     AliError("Cannot obtain AliEMCALGeometry instance.");
91     return;
92   }
93
94   fNSuperModule = geom->GetNumberOfSuperModules();
95
96   if(fDataType == kSurvey) {
97
98     AliSurveyObj *s1 = new AliSurveyObj();
99     s1->FillFromLocalFile(txtFileName);
100     TObjArray* points = s1->GetData();
101     InitSuperModuleData(points);
102
103   } else {
104
105     //Use a dummy file that stores x,y,z of the center of each SM
106     //useful for testing...
107     std::ifstream inputFile(txtFileName.Data());
108     if (!inputFile) {
109       AliError(("Cannot open the survey file " + txtFileName).Data());
110       return;
111     }
112     Int_t dummyInt = 0;
113     Double_t *xReal = new Double_t[fNSuperModule];
114     Double_t *yReal = new Double_t[fNSuperModule];
115     Double_t *zReal = new Double_t[fNSuperModule];
116     Double_t *psiReal = new Double_t[fNSuperModule];
117     Double_t *thetaReal = new Double_t[fNSuperModule];
118     Double_t *phiReal = new Double_t[fNSuperModule];
119
120     for (Int_t i = 0; i < fNSuperModule; ++i) {
121       if (!inputFile) {
122         AliError("Error while reading input file.");
123         delete [] xReal;
124         delete [] yReal;
125         delete [] zReal;
126         delete [] psiReal;
127         delete [] thetaReal;
128         delete [] phiReal;
129         return;
130       }
131       inputFile>>dummyInt>>xReal[i]>>yReal[i]>>zReal[i]>>psiReal[i]>>thetaReal[i]>>phiReal[i];
132     }
133
134     InitSuperModuleData(xReal, yReal, zReal, psiReal, thetaReal, phiReal);
135
136     delete [] xReal;
137     delete [] yReal;
138     delete [] zReal;
139     delete [] psiReal;
140     delete [] thetaReal;
141     delete [] phiReal;
142
143   } //kDummy way of doing it
144
145 }
146
147 //____________________________________________________________________________
148 AliEMCALSurvey::~AliEMCALSurvey()
149 {
150   //destructor
151   delete [] fSuperModuleData;
152 }
153
154 //____________________________________________________________________________
155 void AliEMCALSurvey::CreateAliAlignObjParams(TClonesArray &array)
156 {
157   //Create AliAlignObjParams.
158   const AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance();
159   if (!geom) {
160     AliError("Cannot obtain AliEMCALGeometry instance.");
161     return;
162   }
163
164   if (!gGeoManager) {
165     AliWarning("Cannot create local transformations for supermodules - gGeoManager does not exist.");
166     AliInfo("Null shifts and rotations will be created instead.");
167     return CreateNullObjects(array, geom);
168   }
169
170   Int_t arrayInd = array.GetEntries(), iIndex = 0;
171   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
172   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iIndex);
173   AliAlignObjParams* myobj = 0x0;
174
175   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
176     TString smodName(TString::Format("EMCAL/FullSupermodule%d", smodnum+1));
177     if(geom->GetKey110DEG() && smodnum >= 10) {
178       smodName = "EMCAL/HalfSupermodule";
179       smodName += (smodnum-10+1);
180     }    
181     AliEMCALSuperModuleDelta t(GetSuperModuleTransformation(smodnum));
182
183     ///////////////////////////////
184     // JLK 13-July-2010
185     //
186     // VERY IMPORTANT!!!!
187     //
188     // All numbers were calculated in ALICE global c.s., which means
189     // that the last argument in the creation of AliAlignObjParams
190     // MUST BE set to true
191     //////////////////////////////
192     new(array[arrayInd])
193       AliAlignObjParams(
194                         smodName.Data(), volid, 
195                         t.fXShift, t.fYShift, t.fZShift, 
196                         -t.fPsi, -t.fTheta, -t.fPhi, 
197                         true
198                         );
199     ++arrayInd;
200     myobj = (AliAlignObjParams*)array.UncheckedAt(smodnum);
201     printf("==== AliAlignObjParams for SM %d ====\n",smodnum);
202     myobj->Print("");
203
204   }
205
206 }
207
208 //____________________________________________________________________________
209 void AliEMCALSurvey::CreateNullObjects(TClonesArray &array, const AliEMCALGeometry *geom)const
210 {
211   //Create null shifts and rotations.
212   Int_t arrayInd = array.GetEntries(), iIndex = 0;
213   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
214   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iIndex);
215
216   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
217     TString smodName(TString::Format("EMCAL/FullSupermodule%d", smodnum+1));
218     if(geom->GetKey110DEG() && smodnum >= 10) {
219       smodName = "EMCAL/HalfSupermodule";
220       smodName += (smodnum-10+1);
221     }
222     new(array[arrayInd]) AliAlignObjParams(smodName.Data(), volid, 0., 0., 0., 0., 0., 0., true);
223     ++arrayInd;
224   }
225 }
226
227 //____________________________________________________________________________
228 AliEMCALSurvey::AliEMCALSuperModuleDelta AliEMCALSurvey::GetSuperModuleTransformation(Int_t supModIndex)const
229 {
230   //Supermodule transformation.
231   AliEMCALSuperModuleDelta t = {0., 0., 0., 0., 0., 0.};
232   if (!fSuperModuleData)
233     return t;
234
235   return fSuperModuleData[supModIndex];
236 }
237
238 //____________________________________________________________________________
239 void AliEMCALSurvey::InitSuperModuleData(const TObjArray *svypts)
240 {
241   //This method uses the data points from the EMCAL survey and CATIA program to
242   //create the alignment matrices.  Only valid for (installed)
243   //SM, others will have null objects
244
245   /*--------------------------------------
246     The bottom edges of the strip modules
247     define the active area of the EMCAL, but
248     in software we have a box to hold them which
249     is longer than that.  We need to convert
250     info about the position of the corners of the
251     bottom of the active area to the center of
252     the software box that contains the strip
253     modules.
254
255     View from beam axis up to EMCAL
256               Ai                Ci
257
258         0,1         0,0 1,0          1,1
259     xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
260     x   x             x x              x   x
261     x   x    % *      x x      * %     x   x
262     x   x             x x              x   x
263     xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
264         1,1         1,0 0,0          0,1
265     <--> = added length                 <--> = added length
266
267     * represents the center of the active area
268     % represents the center of the full box (with added length)
269
270          View from side of topmost SM
271
272               Ai                Ci
273
274         0,1         0,0 1,0          1,1
275     xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
276     x   x    % *      x x      % *     x   x
277     xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
278         1,1         1,0 0,0          0,1
279     <--> = added length                 <--> = added length
280
281     * represents the center of the active area
282     % represents the center of the full box (with added length)
283
284     -------------------------------------*/
285
286   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
287   //Center of supermodules
288   Float_t *pars = geom->GetSuperModulesPars();
289   Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
290   Double_t phi, phiRad, xpos, ypos, zpos;
291
292   AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
293   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
294     AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
295     phiRad = geom->GetPhiCenterOfSM(smodnum); //comes in radians
296     phi = phiRad*180./TMath::Pi(); //need degrees for AliAlignObjParams
297     xpos = rpos * TMath::Cos(phiRad);
298     ypos = rpos * TMath::Sin(phiRad);
299     zpos = pars[2];
300     if(geom->GetKey110DEG() && smodnum >= 10) {
301       xpos += (pars[1]/2. * TMath::Sin(phiRad));
302       ypos -= (pars[1]/2. * TMath::Cos(phiRad));
303     }
304     smc.fX1 = xpos;
305     smc.fY1 = ypos;
306     smc.fPhi = phi; //degrees
307     smc.fTheta = 0.; //degrees
308     smc.fPsi = 0.; //degrees
309     if(smodnum%2==0) {
310       smc.fZ1 = zpos;
311     } else {
312       smc.fZ1 = -zpos;
313     }
314
315     //printf("PHI OF IDEAL SM = %.2f\n",smc.fPhi);
316
317   }
318
319   //Real coordinates of center and rotation angles need to be computed
320   //from the survey/CATIA points
321
322   char substr[100];
323   AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
324   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
325     AliEMCALSuperModuleCoords &smc = realSM[smodnum];
326     Double_t zLength = pars[2]*2.; //length of SM in z from software
327     Double_t halfHeight = pars[0]; //half the height of the SM in y direction
328
329     sprintf(substr,"4096%d",smodnum);
330     //retrieve components of four face points and determine average position of center
331     //in x,y,z
332
333     std::vector<Double_t> xval;
334     std::vector<Double_t> yval;
335     std::vector<Double_t> zval;
336     
337     for(Int_t i = 0; i < svypts->GetEntries(); i++) {
338       AliSurveyPoint* pt = (AliSurveyPoint*)svypts->At(i);
339       TString ptname = pt->GetPointName();
340       if(ptname.Contains(substr)) {
341         //Note: order of values is 00, 01, 10, 11
342         xval.push_back(pt->GetX()*100.); //convert m to cm
343         yval.push_back(pt->GetY()*100.); 
344         zval.push_back(pt->GetZ()*100.); 
345       }
346     }
347
348     //compute center of active area of each SM on bottome face from survey points
349     Double_t activeX = ((xval[0] + (xval[2] - xval[0])/2.)        //x00 and x10
350                         +(xval[1] + (xval[3] - xval[1])/2.) ) /2.; //x01 and x11
351     
352 //    Double_t activeY = ((yval[0] + (yval[2] - yval[0])/2.)
353 //                      +(yval[1] + (yval[3] - yval[1])/2.) ) /2.;
354 //    
355 //    Double_t activeZ = ((zval[0] + (zval[2] - zval[0])/2.)
356 //                      +(zval[1] + (zval[3] - zval[1])/2.) ) /2.;
357     
358     //printf("Bottom Center of active area of SM %s: %.2f, %.2f, %.2f\n",substr,activeX,activeY,activeZ);
359     
360     //compute angles for each SM
361     //rotation about each axis
362     //phi = angle in x-y plane
363     
364     Double_t realphi = 0.;
365     //Note: this is phi wrt y axis.  To get phi wrt to x, add pi/2
366     if(smodnum%2 == 0) {
367       realphi = (TMath::ATan((yval[2] - yval[0])/(xval[2] - xval[0])) 
368                  +TMath::ATan((yval[3] - yval[1])/(xval[3] - xval[1])) )/2.;
369     } else {
370       realphi = (TMath::ATan((yval[0] - yval[2])/(xval[0] - xval[2]))
371                  +TMath::ATan((yval[1] - yval[3])/(xval[1] - xval[3])) )/2.;
372     }
373
374     //NOTE: Psi angle is always zero because the two z values being
375     //subtracted are exactly the same, but just in case that could change...
376     //psi = angle in x-z plane
377     Double_t realpsi = (TMath::ATan((zval[0] - zval[2])/(xval[2] - xval[0]))
378                         +TMath::ATan((zval[1] - zval[3])/(xval[3] - xval[1])) )/2.;
379     
380     //theta = angle in y-z plane
381     Double_t realtheta = TMath::Pi()/2. 
382                        - (TMath::ATan((zval[2] - zval[3])/(yval[3] - yval[2]))
383                           +TMath::ATan((zval[0] - zval[1])/(yval[1] - yval[0])) )/2.;
384
385     //printf("Old edge of %s 01: %.2f, %.2f, %.2f\n",substr,xval[1],yval[1],zval[1]);
386     //printf("Old edge of %s 11: %.2f, %.2f, %.2f\n",substr,xval[3],yval[3],zval[3]);
387
388     //Now calculate the center of the box in z with length added to the 01
389     //and 11 corners, corrected by the theta angle
390     Double_t activeLength = TMath::Abs(((zval[1] - zval[0]) + (zval[3] - zval[2]))/2.);
391     //printf("ACTIVE LENGTH = %.2f\n",activeLength);
392     if(smodnum%2 == 0) {
393       yval[1] += (zLength - activeLength)*sin(realtheta);
394       yval[3] += (zLength - activeLength)*sin(realtheta);
395       zval[1] += (zLength - activeLength)*cos(realtheta);
396       zval[3] += (zLength - activeLength)*cos(realtheta);
397     } else {
398       yval[1] -= (zLength - activeLength)*sin(realtheta);
399       yval[3] -= (zLength - activeLength)*sin(realtheta);
400       zval[1] -= (zLength - activeLength)*cos(realtheta);
401       zval[3] -= (zLength - activeLength)*cos(realtheta);
402     }
403
404     //printf("New extended edge of %s 01: %.2f, %.2f, %.2f\n",substr,xval[1],yval[1],zval[1]);            
405     //printf("New extended edge of %s 11: %.2f, %.2f, %.2f\n",substr,xval[3],yval[3],zval[3]);
406
407     //Compute the center of the bottom of the box in x,y,z
408     Double_t realX = activeX;    
409     Double_t realY = ((yval[0] + (yval[2] - yval[0])/2.)
410                         +(yval[1] + (yval[3] - yval[1])/2.) ) /2.;    
411     Double_t realZ = ((zval[0] + (zval[2] - zval[0])/2.)
412                         +(zval[1] + (zval[3] - zval[1])/2.) ) /2.;
413
414
415     //printf("Bottom Center of SM %s Box: %.2f, %.2f, %.2f\n",substr,realX,realY,realZ);
416
417     //correct the SM centers so that we have the center of the box in
418     //x,y using the phi,theta angles                   
419     realX += halfHeight*TMath::Cos(TMath::Pi()/2+realphi);
420     realY += halfHeight*(TMath::Sin(TMath::Pi()/2+realphi) + TMath::Sin(realtheta));
421     realZ += halfHeight*TMath::Cos(TMath::Pi()/2-realtheta);
422
423     //printf("Rotation angles of SM %s: %.4f, %.4f, %.4f\n",substr,realphi*TMath::RadToDeg(),realpsi*TMath::RadToDeg(),realtheta*TMath::RadToDeg());
424     //printf("Middle of SM %s: %.2f, %.2f, %.2f\n\n",substr,realX,realY,realZ);
425
426     smc.fX1 = realX;
427     smc.fY1 = realY;
428     smc.fZ1 = realZ;
429
430     smc.fPhi = 90. + realphi*TMath::RadToDeg();
431     smc.fTheta = 0. + realtheta*TMath::RadToDeg();
432     smc.fPsi = 0. + realpsi*TMath::RadToDeg();
433
434   }//loop over supermodules
435   
436   fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
437   
438   for (Int_t i = 0; i < fNSuperModule; ++i) {
439     const AliEMCALSuperModuleCoords &real = realSM[i];
440     const AliEMCALSuperModuleCoords &ideal = idealSM[i];
441     AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
442     t.fXShift = real.fX1 - ideal.fX1;
443     t.fYShift = real.fY1 - ideal.fY1;
444     t.fZShift = real.fZ1 - ideal.fZ1;
445     t.fPhi = real.fPhi - ideal.fPhi;
446     t.fTheta = real.fTheta - ideal.fTheta;
447     t.fPsi = real.fPsi - ideal.fPsi;
448
449     printf("===================== SM %d =======================\n",i);
450     printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
451     printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
452     printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
453     printf("real theta (%.2f) - ideal theta (%.2f) = shift in theta %.2f\n",real.fTheta,ideal.fTheta,t.fTheta);
454     printf("real psi (%.2f) - ideal psi (%.2f) = shift in psi %.2f\n",real.fPsi,ideal.fPsi,t.fPsi);
455     printf("real phi (%.2f) - ideal phi (%.2f) = shift in phi %.2f\n",real.fPhi,ideal.fPhi,t.fPhi);
456     printf("===================================================\n");    
457   }
458
459   delete [] realSM;
460   delete [] idealSM;
461 }
462
463
464 //____________________________________________________________________________
465 void AliEMCALSurvey::InitSuperModuleData(const Double_t *xReal, const Double_t *yReal, 
466                                          const Double_t *zReal, const Double_t *psiReal,
467                                          const Double_t *thetaReal, const Double_t *phiReal)
468 {
469   ///////////////////////////////////////
470   //Dummy method just takes the inputted values and applies them
471   //
472   //Useful for testing small changes
473   //////////////////////////////////////
474   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
475   //Center of supermodules
476   Float_t *pars = geom->GetSuperModulesPars();
477   Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
478   Double_t phi, phiRad, xpos, ypos, zpos;
479
480   zpos = pars[2];
481
482   AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
483   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
484     AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
485     phiRad = geom->GetPhiCenterOfSM(smodnum); //comes in radians
486     phi = phiRad*180./TMath::Pi(); //need degrees for AliAlignObjParams
487     xpos = rpos * TMath::Cos(phiRad);
488     ypos = rpos * TMath::Sin(phiRad);
489     if(geom->GetKey110DEG() && smodnum >= 10) {
490       xpos += (pars[1]/2. * TMath::Sin(phiRad));
491       ypos -= (pars[1]/2. * TMath::Cos(phiRad));
492     }
493     smc.fX1 = xpos;
494     smc.fY1 = ypos;
495
496     smc.fPhi = phi; //degrees
497     smc.fTheta = 0.; //degrees
498     smc.fPsi = 0.; //degrees
499
500     if(smodnum%2==0) {
501       smc.fZ1 = zpos;
502     } else {
503       smc.fZ1 = -zpos;
504     }
505
506   }
507
508   AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
509   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
510     AliEMCALSuperModuleCoords &smc = realSM[smodnum];
511     smc.fX1    = xReal[smodnum];  
512     smc.fY1    = yReal[smodnum];  
513     smc.fZ1    = zReal[smodnum];  
514     smc.fTheta = thetaReal[smodnum];
515     smc.fPsi   = psiReal[smodnum];
516     smc.fPhi   = phiReal[smodnum];
517   }
518   
519   fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
520   
521   for (Int_t i = 0; i < fNSuperModule; ++i) {
522     const AliEMCALSuperModuleCoords &real = realSM[i];
523     const AliEMCALSuperModuleCoords &ideal = idealSM[i];
524     AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
525     t.fTheta = real.fTheta - ideal.fTheta;
526     t.fPsi = 0.;
527     t.fPhi = real.fPhi - ideal.fPhi;
528     t.fXShift = real.fX1 - ideal.fX1;
529     t.fYShift = real.fY1 - ideal.fY1;
530     t.fZShift = real.fZ1 - ideal.fZ1;
531
532     printf("===================== SM %d =======================\n",i);
533     printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
534     printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
535     printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
536     printf("real theta (%.2f) - ideal theta (%.2f) = shift in theta %.2f\n",real.fTheta,ideal.fTheta,t.fTheta);
537     printf("real psi (%.2f) - ideal psi (%.2f) = shift in psi %.2f\n",real.fPsi,ideal.fPsi,t.fPsi);
538     printf("real phi (%.2f) - ideal phi (%.2f) = shift in phi %.2f\n",real.fPhi,ideal.fPhi,t.fPhi);
539     printf("===================================================\n");    
540   }
541
542   delete [] realSM;
543   delete [] idealSM;
544 }
545