]> git.uio.no Git - u/mrichter/AliRoot.git/blob - EMCAL/AliEMCALSurvey.cxx
Remove dependency on AliRecPoint from AliEMCALGeometry
[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
174   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
175     TString smodName(TString::Format("EMCAL/FullSupermodule%d", smodnum+1));
176     if(geom->GetKey110DEG() && smodnum >= 10) {
177       smodName = "EMCAL/HalfSupermodule";
178       smodName += (smodnum-10+1);
179     }    
180     AliEMCALSuperModuleDelta t(GetSuperModuleTransformation(smodnum));
181     new(array[arrayInd])
182       AliAlignObjParams(
183                         smodName.Data(), volid, 
184                         t.fXShift, t.fYShift, t.fZShift, 
185                         -t.fPsi, -t.fTheta, -t.fPhi, 
186                         false
187                         );
188     ++arrayInd;
189
190   }
191
192 }
193
194 //____________________________________________________________________________
195 void AliEMCALSurvey::CreateNullObjects(TClonesArray &array, const AliEMCALGeometry *geom)const
196 {
197   //Create null shifts and rotations.
198   Int_t arrayInd = array.GetEntries(), iIndex = 0;
199   AliGeomManager::ELayerID iLayer = AliGeomManager::kInvalidLayer;
200   UShort_t volid = AliGeomManager::LayerToVolUID(iLayer,iIndex);
201
202   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
203     TString smodName(TString::Format("EMCAL/FullSupermodule%d", smodnum+1));
204     if(geom->GetKey110DEG() && smodnum >= 10) {
205       smodName = "EMCAL/HalfSupermodule";
206       smodName += (smodnum-10+1);
207     }
208     new(array[arrayInd]) AliAlignObjParams(smodName.Data(), volid, 0., 0., 0., 0., 0., 0., true);
209     ++arrayInd;
210   }
211 }
212
213 //____________________________________________________________________________
214 AliEMCALSurvey::AliEMCALSuperModuleDelta AliEMCALSurvey::GetSuperModuleTransformation(Int_t supModIndex)const
215 {
216   //Supermodule transformation.
217   AliEMCALSuperModuleDelta t = {0., 0., 0., 0., 0., 0.};
218   if (!fSuperModuleData)
219     return t;
220
221   return fSuperModuleData[supModIndex];
222 }
223
224 //____________________________________________________________________________
225 void AliEMCALSurvey::InitSuperModuleData(const TObjArray *svypts)
226 {
227   //This method uses the data points from the EMCAL survey and CATIA program to
228   //create the alignment matrices.  Only valid for (installed)
229   //SM, others will have null objects
230
231   /*--------------------------------------
232     The bottom edges of the strip modules
233     define the active area of the EMCAL, but
234     in software we have a box to hold them which
235     is longer than that.  We need to convert
236     info about the position of the corners of the
237     bottom of the active area to the center of
238     the software box that contains the strip
239     modules.
240
241     View from beam axis up to EMCAL
242               Ai                Ci
243
244         0,1         0,0 1,0          1,1
245     xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
246     x   x             x x              x   x
247     x   x    % *      x x      * %     x   x
248     x   x             x x              x   x
249     xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
250         1,1         1,0 0,0          0,1
251     <--> = added length                 <--> = added length
252
253     * represents the center of the active area
254     % represents the center of the full box (with added length)
255
256          View from side of topmost SM
257
258               Ai                Ci
259
260         0,1         0,0 1,0          1,1
261     xxxxxxxxxxxxxxxxxxx xxxxxxxxxxxxxxxxxxxx
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     -------------------------------------*/
271
272   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
273   //Center of supermodules
274   Float_t *pars = geom->GetSuperModulesPars();
275   Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
276   Double_t phi, phiRad, xpos, ypos, zpos;
277
278   AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
279   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
280     AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
281     phiRad = geom->GetPhiCenterOfSM(smodnum); //comes in radians
282     phi = phiRad*180./TMath::Pi(); //need degrees for AliAlignObjParams
283     xpos = rpos * TMath::Cos(phiRad);
284     ypos = rpos * TMath::Sin(phiRad);
285     zpos = pars[2];
286     if(geom->GetKey110DEG() && smodnum >= 10) {
287       xpos += (pars[1]/2. * TMath::Sin(phiRad));
288       ypos -= (pars[1]/2. * TMath::Cos(phiRad));
289     }
290     smc.fX1 = xpos;
291     smc.fY1 = ypos;
292     smc.fPhi = phi; //degrees
293     smc.fTheta = 0.; //degrees
294     smc.fPsi = 0.; //degrees
295     if(smodnum%2==0) {
296       smc.fZ1 = zpos;
297     } else {
298       smc.fZ1 = -zpos;
299     }
300
301     //printf("PHI OF IDEAL SM = %.2f\n",smc.fPhi);
302
303   }
304
305   //Real coordinates of center and rotation angles need to be computed
306   //from the survey/CATIA points
307
308   char substr[100];
309   AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
310   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
311     AliEMCALSuperModuleCoords &smc = realSM[smodnum];
312     Double_t zLength = pars[2]*2.; //length of SM in z from software
313     Double_t halfHeight = pars[0]; //half the height of the SM in y direction
314
315     sprintf(substr,"4096%d",smodnum);
316     //retrieve components of four face points and determine average position of center
317     //in x,y,z
318
319     std::vector<Double_t> xval;
320     std::vector<Double_t> yval;
321     std::vector<Double_t> zval;
322     
323     for(Int_t i = 0; i < svypts->GetEntries(); i++) {
324       AliSurveyPoint* pt = (AliSurveyPoint*)svypts->At(i);
325       TString ptname = pt->GetPointName();
326       if(ptname.Contains(substr)) {
327         //Note: order of values is 00, 01, 10, 11
328         xval.push_back(pt->GetX()*100.); //convert m to cm
329         yval.push_back(pt->GetY()*100.); 
330         zval.push_back(pt->GetZ()*100.); 
331       }
332     }
333
334     //compute center of active area of each SM on bottome face from survey points
335     Double_t activeX = ((xval[0] + (xval[2] - xval[0])/2.)        //x00 and x10
336                         +(xval[1] + (xval[3] - xval[1])/2.) ) /2.; //x01 and x11
337     
338     Double_t activeY = ((yval[0] + (yval[2] - yval[0])/2.)
339                         +(yval[1] + (yval[3] - yval[1])/2.) ) /2.;
340     
341     Double_t activeZ = ((zval[0] + (zval[2] - zval[0])/2.)
342                         +(zval[1] + (zval[3] - zval[1])/2.) ) /2.;
343     
344     //printf("Bottom Center of active area of SM %s: %.2f, %.2f, %.2f\n",substr,activeX,activeY,activeZ);
345     
346     //compute angles for each SM
347     //rotation about each axis
348     //phi = angle in x-y plane
349     
350     Double_t realphi = 0.;
351     //Note: this is phi wrt y axis.  To get phi wrt to x, add pi/2
352     if(smodnum%2 == 0) {
353       realphi = (TMath::ATan((yval[2] - yval[0])/(xval[2] - xval[0])) 
354                  +TMath::ATan((yval[3] - yval[1])/(xval[3] - xval[1])) )/2.;
355     } else {
356       realphi = (TMath::ATan((yval[0] - yval[2])/(xval[0] - xval[2]))
357                  +TMath::ATan((yval[1] - yval[3])/(xval[1] - xval[3])) )/2.;
358     }
359
360     //NOTE: Psi angle is always zero because the two z values being
361     //subtracted are exactly the same, but just in case that could change...
362     //psi = angle in x-z plane
363     Double_t realpsi = (TMath::ATan((zval[0] - zval[2])/(xval[2] - xval[0]))
364                         +TMath::ATan((zval[1] - zval[3])/(xval[3] - xval[1])) )/2.;
365     
366     //theta = angle in y-z plane
367     Double_t realtheta = TMath::Pi()/2. 
368                        - (TMath::ATan((zval[2] - zval[3])/(yval[3] - yval[2]))
369                           +TMath::ATan((zval[0] - zval[1])/(yval[1] - yval[0])) )/2.;
370
371     //printf("Old edge of %s 01: %.2f, %.2f, %.2f\n",substr,xval[1],yval[1],zval[1]);
372     //printf("Old edge of %s 11: %.2f, %.2f, %.2f\n",substr,xval[3],yval[3],zval[3]);
373
374     //Now calculate the center of the box in z with length added to the 01
375     //and 11 corners, corrected by the theta angle
376     Double_t activeLength = TMath::Abs(((zval[1] - zval[0]) + (zval[3] - zval[2]))/2.);
377     //printf("ACTIVE LENGTH = %.2f\n",activeLength);
378     if(smodnum%2 == 0) {
379       yval[1] += (zLength - activeLength)*sin(realtheta);
380       yval[3] += (zLength - activeLength)*sin(realtheta);
381       zval[1] += (zLength - activeLength)*cos(realtheta);
382       zval[3] += (zLength - activeLength)*cos(realtheta);
383     } else {
384       yval[1] -= (zLength - activeLength)*sin(realtheta);
385       yval[3] -= (zLength - activeLength)*sin(realtheta);
386       zval[1] -= (zLength - activeLength)*cos(realtheta);
387       zval[3] -= (zLength - activeLength)*cos(realtheta);
388     }
389
390     //printf("New extended edge of %s 01: %.2f, %.2f, %.2f\n",substr,xval[1],yval[1],zval[1]);            
391     //printf("New extended edge of %s 11: %.2f, %.2f, %.2f\n",substr,xval[3],yval[3],zval[3]);
392
393     //Compute the center of the bottom of the box in x,y,z
394     Double_t realX = activeX;    
395     Double_t realY = ((yval[0] + (yval[2] - yval[0])/2.)
396                         +(yval[1] + (yval[3] - yval[1])/2.) ) /2.;    
397     Double_t realZ = ((zval[0] + (zval[2] - zval[0])/2.)
398                         +(zval[1] + (zval[3] - zval[1])/2.) ) /2.;
399
400
401     //printf("Bottom Center of SM %s Box: %.2f, %.2f, %.2f\n",substr,realX,realY,realZ);
402
403     //correct the SM centers so that we have the center of the box in
404     //x,y using the phi,theta angles                   
405     realX += halfHeight*TMath::Cos(TMath::Pi()/2+realphi);
406     realY += halfHeight*(TMath::Sin(TMath::Pi()/2+realphi) + TMath::Sin(realtheta));
407     realZ += halfHeight*TMath::Cos(TMath::Pi()/2-realtheta);
408
409     //printf("Rotation angles of SM %s: %.4f, %.4f, %.4f\n",substr,realphi*TMath::RadToDeg(),realpsi*TMath::RadToDeg(),realtheta*TMath::RadToDeg());
410     //printf("Middle of SM %s: %.2f, %.2f, %.2f\n\n",substr,realX,realY,realZ);
411
412     smc.fX1 = realX;
413     smc.fY1 = realY;
414     smc.fZ1 = realZ;
415
416     smc.fPhi = 90. + realphi*TMath::RadToDeg();
417     smc.fTheta = 0. + realtheta*TMath::RadToDeg();
418     smc.fPsi = 0. + realpsi*TMath::RadToDeg();
419
420   }//loop over supermodules
421   
422   fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
423   
424   for (Int_t i = 0; i < fNSuperModule; ++i) {
425     const AliEMCALSuperModuleCoords &real = realSM[i];
426     const AliEMCALSuperModuleCoords &ideal = idealSM[i];
427     AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
428     t.fXShift = real.fX1 - ideal.fX1;
429     t.fYShift = real.fY1 - ideal.fY1;
430     t.fZShift = ideal.fZ1 - real.fZ1; //due to z flip for C side
431     if(i%2==0) {
432       t.fZShift *= -1.0;  //correct shift for C side
433     }
434     t.fPhi = real.fPhi - ideal.fPhi;
435     t.fTheta = real.fTheta - ideal.fTheta;
436     t.fPsi = real.fPsi - ideal.fPsi;
437
438     printf("===================== SM %d =======================\n",i);
439     printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
440     printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
441     printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
442     printf("real theta (%.2f) - ideal theta (%.2f) = shift in theta %.2f\n",real.fTheta,ideal.fTheta,t.fTheta);
443     printf("real psi (%.2f) - ideal psi (%.2f) = shift in psi %.2f\n",real.fPsi,ideal.fPsi,t.fPsi);
444     printf("real phi (%.2f) - ideal phi (%.2f) = shift in phi %.2f\n",real.fPhi,ideal.fPhi,t.fPhi);
445     printf("===================================================\n");    
446   }
447
448   delete [] realSM;
449   delete [] idealSM;
450 }
451
452
453 //____________________________________________________________________________
454 void AliEMCALSurvey::InitSuperModuleData(const Double_t *xReal, const Double_t *yReal, 
455                                          const Double_t *zReal, const Double_t *psiReal,
456                                          const Double_t *thetaReal, const Double_t *phiReal)
457 {
458   ///////////////////////////////////////
459   //Dummy method just takes the inputted values and applies them
460   //
461   //Useful for testing small changes
462   //////////////////////////////////////
463   AliEMCALGeometry *geom = AliEMCALGeometry::GetInstance();
464   //Center of supermodules
465   Float_t *pars = geom->GetSuperModulesPars();
466   Double_t rpos = (geom->GetEnvelop(0) + geom->GetEnvelop(1))/2.;
467   Double_t phi, phiRad, xpos, ypos, zpos;
468
469   zpos = pars[2];
470
471   AliEMCALSuperModuleCoords *idealSM = new AliEMCALSuperModuleCoords[fNSuperModule];
472   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
473     AliEMCALSuperModuleCoords &smc = idealSM[smodnum];
474     phiRad = geom->GetPhiCenterOfSM(smodnum); //comes in radians
475     phi = phiRad*180./TMath::Pi(); //need degrees for AliAlignObjParams
476     xpos = rpos * TMath::Cos(phiRad);
477     ypos = rpos * TMath::Sin(phiRad);
478     if(geom->GetKey110DEG() && smodnum >= 10) {
479       xpos += (pars[1]/2. * TMath::Sin(phiRad));
480       ypos -= (pars[1]/2. * TMath::Cos(phiRad));
481     }
482     smc.fX1 = xpos;
483     smc.fY1 = ypos;
484
485     smc.fPhi = phi; //degrees
486     smc.fTheta = 0.; //degrees
487     smc.fPsi = 0.; //degrees
488
489     if(smodnum%2==0) {
490       smc.fZ1 = zpos;
491     } else {
492       smc.fZ1 = -zpos;
493     }
494
495   }
496
497   AliEMCALSuperModuleCoords *realSM = new AliEMCALSuperModuleCoords[fNSuperModule];
498   for (Int_t smodnum = 0; smodnum < geom->GetNumberOfSuperModules(); ++smodnum) {
499     AliEMCALSuperModuleCoords &smc = realSM[smodnum];
500     smc.fX1    = xReal[smodnum];  
501     smc.fY1    = yReal[smodnum];  
502     smc.fZ1    = zReal[smodnum];  
503     smc.fTheta = thetaReal[smodnum];
504     smc.fPsi   = psiReal[smodnum];
505     smc.fPhi   = phiReal[smodnum];
506   }
507   
508   fSuperModuleData = new AliEMCALSuperModuleDelta[fNSuperModule];
509   
510   for (Int_t i = 0; i < fNSuperModule; ++i) {
511     const AliEMCALSuperModuleCoords &real = realSM[i];
512     const AliEMCALSuperModuleCoords &ideal = idealSM[i];
513     AliEMCALSuperModuleDelta &t = fSuperModuleData[i];
514     t.fTheta = real.fTheta - ideal.fTheta;
515     t.fPsi = 0.;
516     t.fPhi = real.fPhi - ideal.fPhi;
517     t.fXShift = real.fX1 - ideal.fX1;
518     t.fYShift = real.fY1 - ideal.fY1;
519     t.fZShift = real.fZ1 - ideal.fZ1;
520
521     printf("===================== SM %d =======================\n",i);
522     printf("real x (%.2f) - ideal x (%.2f) = shift in x (%.2f)\n",real.fX1,ideal.fX1,t.fXShift);
523     printf("real y (%.2f) - ideal y (%.2f) = shift in y (%.2f)\n",real.fY1,ideal.fY1,t.fYShift);
524     printf("real z (%.2f) - ideal z (%.2f) = shift in z (%.2f)\n",real.fZ1,ideal.fZ1,t.fZShift);
525     printf("real theta (%.2f) - ideal theta (%.2f) = shift in theta %.2f\n",real.fTheta,ideal.fTheta,t.fTheta);
526     printf("real psi (%.2f) - ideal psi (%.2f) = shift in psi %.2f\n",real.fPsi,ideal.fPsi,t.fPsi);
527     printf("real phi (%.2f) - ideal phi (%.2f) = shift in phi %.2f\n",real.fPhi,ideal.fPhi,t.fPhi);
528     printf("===================================================\n");    
529   }
530
531   delete [] realSM;
532   delete [] idealSM;
533 }
534