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