TPC/AliTPCCalibGlobalMisalignment.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
  17 ////////////////////////////////////////////////////////////////////////////
  18 //                                                                        //
  19 // AliTPCCalibGlobalMisalignment class                                        //
  20 // The class calculates the space point distortions due to simple         //
  21 // misalignments like shifts in caresian coordinates or a rotation        //
  22 // of the TPC read out planes (A and C side)                              //
  23 // Optionaly possible to use it for visualization of the alignemnt form the Alignment OCDB //
  24 // fUseGeoManager has to be set to kTRUE to enable this option
  25 //                                                                        //
  26 // date: 06/05/2010                                                       //
  27 // Authors: Stefan Rossegger, Jim Thomas, Magnus Mager                    //
  28 ////////////////////////////////////////////////////////////////////////////
  29
  30 #include "AliTPCCalibGlobalMisalignment.h"
  31 #include "TMath.h"
  32 #include "TGeoMatrix.h"
  33 #include "AliTPCROC.h"
  34 #include "AliTPCcalibDB.h"
  35 #include "AliTPCParam.h"
  36 #include <TGeoPhysicalNode.h>
  37 //
  38 #include "AliAlignObjParams.h"
  39 #include "AliGeomManager.h"
  40 #include "AliCDBManager.h"
  41 #include "AliCDBEntry.h"
  42
  43 AliTPCCalibGlobalMisalignment::AliTPCCalibGlobalMisalignment()
  44   : AliTPCCorrection("mialign","Misalignment"),
  45     fXShift(0.),fYShift(0.),fZShift(0.),
  46     fRotPhiA(0.),fRotPhiC(0.),
  47     fdRPhiOffsetA(0.),
  48     fdRPhiOffsetC(0.),
  49     fQuadrantQ0(0),   //OROC medium pads -delta ly+ - ly - shift
  50     fQuadrantRQ0(0),  //OROC medium pads -delta ly+ - ly - rotation
  51     fQuadrantQ1(0),   //OROC long   pads -delta ly+ - ly - shift
  52     fQuadrantQ2(0),   //OROC long   pads -shift
  53     fQuadrantRQ1(0),  //OROC long   pads -delta ly+ - ly - rotation
  54     fQuadrantRQ2(0),  //OROC long   pads -rotation
  55     fMatrixGlobal(0), // global Alignment common
  56     fMatrixGlobalDelta(0), // global Alignment Delta A side-c side
  57     fArraySector(0)   // fArraySector
  58 {
  59   //
  60   // default constructor
  61   //
  62 }
  63
  64 AliTPCCalibGlobalMisalignment::~AliTPCCalibGlobalMisalignment() {
  65   //
  66   // destructor
  67   //
  68   delete fQuadrantQ0;   //OROC medium pads -delta ly+ - ly - shift
  69   delete fQuadrantRQ0;  //OROC medium pads -delta ly+ - ly - rotation
  70   delete fQuadrantQ1;   //OROC long   pads -delta ly+ - ly - shift
  71   delete fQuadrantQ2;   //OROC long   pads -shift
  72   delete fQuadrantRQ1;  //OROC long   pads -delta ly+ - ly - rotation
  73   delete fQuadrantRQ2;  //OROC long   pads -rotation
  74   delete fMatrixGlobal; // global matrix
  75   delete fMatrixGlobal; // global matrix
  76   delete fArraySector;  // sector matrices
  77 }
  78
  79 void AliTPCCalibGlobalMisalignment::SetQuadranAlign(const TVectorD *quadrantQ0, const TVectorD *quadrantRQ0, const TVectorD *quadrantQ1,const TVectorD *quadrantRQ1,  const TVectorD *quadrantQ2,  const TVectorD *quadrantRQ2){
  80   //
  81   // Set quadrant alignment
  82   // 6 vectors for 36 (super) sectors
  83   //
  84   if (quadrantQ0) fQuadrantQ0   = new TVectorD(*quadrantQ0);
  85   if (quadrantRQ0) fQuadrantRQ0 = new TVectorD(*quadrantRQ0);
  86   //
  87   if (quadrantQ1) fQuadrantQ1   = new TVectorD(*quadrantQ1);
  88   if (quadrantQ1) fQuadrantRQ1  = new TVectorD(*quadrantRQ1);
  89   if (quadrantQ2) fQuadrantQ2   = new TVectorD(*quadrantQ2);
  90   if (quadrantQ2) fQuadrantRQ2  = new TVectorD(*quadrantRQ2);
  91 }
  92
  93
  94
  95 void AliTPCCalibGlobalMisalignment::SetAlignGlobal(const TGeoMatrix * matrixGlobal){
  96   //
  97   // Set global misalignment
  98   // Object is OWNER
  99   //
 100   if (fMatrixGlobal) delete fMatrixGlobal;
 101   fMatrixGlobal=0;
 102   if (matrixGlobal) fMatrixGlobal = new TGeoHMatrix(*matrixGlobal);
 103 }
 104
 105 void AliTPCCalibGlobalMisalignment::SetAlignGlobalDelta(const TGeoMatrix * matrixGlobalDelta){
 106   //
 107   // Set global misalignment
 108   // Object is OWNER
 109   //
 110   if (fMatrixGlobalDelta) delete fMatrixGlobalDelta;
 111   fMatrixGlobalDelta=0;
 112   if (matrixGlobalDelta) fMatrixGlobalDelta = new TGeoHMatrix(*matrixGlobalDelta);
 113 }
 114
 115 void AliTPCCalibGlobalMisalignment::SetAlignSectors(const TObjArray *arraySector){
 116   //
 117   // Set misalignment TObjArray of TGeoMatrices  - for each sector
 118   // Object is OWNER
 119   //
 120   if (fArraySector) delete fArraySector;
 121   fArraySector=0;
 122   if (arraySector) fArraySector = (TObjArray*)arraySector->Clone();
 123 }
 124
 125
 126 //void AliTPCCalibGlobalMisalignment::Init() {
 127 //  //
 128 // // Initialization funtion
 129 //  //
 130
 131 //  // nothing to be initialized, results of this calibration class will go to the global aligment structure
 132
 133 //}
 134
 135 //void AliTPCCalibGlobalMisalignment::Update(const TTimeStamp &/*timeStamp*/) {
 136 //  //
 137 //  // Update function
 138 //  //
 139 //
 140 //  // nothing to be updated, results of this calibration class will go to the global aligment structure
 141 //
 142 //}
 143
 144
 145
 146 void AliTPCCalibGlobalMisalignment::GetCorrection(const Float_t x[],const Short_t roc,Float_t dx[]) {
 147   //
 148   // Calculates the simple correction due to a shift (in x,y,z) or an rotation of the TPC (around z)
 149   //
 150   static AliTPCROC *tpcRoc =AliTPCROC::Instance();
 151   Double_t xref  = ( tpcRoc->GetPadRowRadii(0,0)+tpcRoc->GetPadRowRadii(36,tpcRoc->GetNRows(36)-1))*0.5;
 152   Double_t xquadrant  = tpcRoc->GetPadRowRadii(36,53); // row 53 from uli
 153   Double_t xIO  = ( tpcRoc->GetPadRowRadii(0,tpcRoc->GetNRows(0)-1)+tpcRoc->GetPadRowRadii(36,0))*0.5;
 154   Double_t r=0, phi=0;
 155   r   = TMath::Sqrt( x[0]*x[0] + x[1]*x[1] );
 156   phi = TMath::ATan2(x[1],x[0]);
 157   // Getsector number
 158   Double_t sec=TMath::Nint(-0.5+(phi*9./TMath::Pi()));
 159   if (sec<0) sec+=18;
 160   Int_t isec = TMath::Nint(sec);
 161   if  (roc%36>=18) isec+=18;
 162   //
 163   // Get the point on the local coordiante frame
 164   //
 165   Double_t alpha=(sec+0.5)*TMath::Pi()/9;
 166   Double_t pos[3]={0,0,x[2]};
 167   pos[0]=  TMath::Cos(alpha)*x[0]+TMath::Sin(alpha)*x[1];
 168   pos[1]= -TMath::Sin(alpha)*x[0]+TMath::Cos(alpha)*x[1];
 169   if (pos[0]>tpcRoc->GetPadRowRadiiUp(0))  isec+=36;
 170
 171   //
 172   // apply quadrant alignment if available - in local coordinate frame
 173   //
 174   //
 175   Double_t posQG[3]={x[0],x[1],x[2]};
 176   {
 177     Double_t dly=0;
 178     Bool_t isQ0 = (pos[0]<xquadrant)&&(pos[0]>xIO);
 179     Bool_t isQ1 = (pos[0]>xquadrant);
 180     Double_t  sign = (pos[1]>0)? 1.: -1.;
 181     if (isQ0){
 182       if (fQuadrantQ0)  dly+=sign*(*fQuadrantQ0)[isec%36];  // shift in cm
 183       if (fQuadrantRQ0) dly+=sign*(*fQuadrantRQ0)[isec%36]*(pos[0]-xref);
 184     }
 185     if (isQ1){
 186       if (fQuadrantQ1)  dly+=sign*(*fQuadrantQ1)[isec%36];  // shift in cm
 187       if (fQuadrantRQ1) dly+=sign*(*fQuadrantRQ1)[isec%36]*(pos[0]-xref);
 188       if (fQuadrantQ2)  dly+=(*fQuadrantQ2)[isec%36];  // shift in cm
 189       if (fQuadrantRQ2) dly+=(*fQuadrantRQ2)[isec%36]*(pos[0]-xref);
 190     }
 191     // Tranform the corrected point to the global frame
 192     posQG[0]=  TMath::Cos(alpha)*pos[0]-TMath::Sin(alpha)*(pos[1]+dly);
 193     posQG[1]=  TMath::Sin(alpha)*pos[0]+TMath::Cos(alpha)*(pos[1]+dly);
 194   }
 195   //
 196   // rotation of the read-out planes
 197   if  (roc%36<18) // A side
 198     phi += fRotPhiA;
 199   else         // C side
 200     phi += fRotPhiC;
 201
 202   // Simply adding a constant dRPHi residual. PURELY FOR CALIBRATION PURPOSES
 203   if  (roc%36<18) // A side
 204     phi += fdRPhiOffsetA/r;
 205   else         // C side
 206     phi += fdRPhiOffsetC/r;
 207
 208   dx[0] = r * TMath::Cos(phi) - x[0];
 209   dx[1] = r * TMath::Sin(phi) - x[1];
 210   dx[2] = 0.;
 211
 212   // Simple shifts
 213   dx[0] -= fXShift;
 214   dx[1] -= fYShift;
 215   dx[2] -= fZShift;
 216   // quadrant shifts
 217   dx[0] += (posQG[0]-x[0]);
 218   dx[1] += (posQG[1]-x[1]);
 219   //
 220   //
 221   if (fMatrixGlobal){
 222     // apply global alignment matrix
 223     Double_t ppos[3]={x[0],x[1],x[2]};
 224     Double_t pposC[3]={x[0],x[1],x[2]};
 225     fMatrixGlobal->LocalToMaster(ppos,pposC);
 226     dx[0]+=pposC[0]-ppos[0];
 227     dx[1]+=pposC[1]-ppos[1];
 228     dx[2]+=pposC[2]-ppos[2];
 229   }
 230   if (fMatrixGlobalDelta){
 231     // apply global alignment matrix A-C Side side
 232     Double_t ppos[3]={x[0],x[1],x[2]};
 233     Double_t pposC[3]={x[0],x[1],x[2]};
 234     fMatrixGlobalDelta->LocalToMaster(ppos,pposC);
 235     Double_t ssign=(roc%36<18) ? 1.:-1.;
 236     dx[0]+=ssign*(pposC[0]-ppos[0]);
 237     dx[1]+=ssign*(pposC[1]-ppos[1]);
 238     dx[2]+=ssign*(pposC[2]-ppos[2]);
 239   }
 240
 241   if (fArraySector){
 242     // apply global alignment matrix
 243     TGeoMatrix  *mat = (TGeoMatrix*)fArraySector->At(isec);
 244     if (mat){
 245       Double_t ppos[3]={x[0],x[1],x[2]};
 246       Double_t pposC[3]={x[0],x[1],x[2]};
 247       mat->LocalToMaster(ppos,pposC);
 248       dx[0]+=pposC[0]-ppos[0];
 249       dx[1]+=pposC[1]-ppos[1];
 250       dx[2]+=pposC[2]-ppos[2];
 251     }
 252   }
 253 }
 254
 255 void AliTPCCalibGlobalMisalignment::Print(Option_t* /*option*/ ) const {
 256   //
 257   // Print function to check the settings
 258   //
 259   printf("%s",GetTitle());
 260   printf(" - Trivial Misalignments for calibration purposes: \n");
 261   printf(" - X-Shift: %1.3f cm, Y-Shift: %1.3f cm, Z-Shift: %1.3f cm \n",fXShift,fYShift,fZShift);
 262   printf(" - Phi-Rotations: A side: %1.5f rad, C side: %1.5f rad\n",fRotPhiA,fRotPhiC);
 263   printf(" - dRPhi offsets: A side: %1.5f cm, C side: %1.5f cm\n",fdRPhiOffsetA,fdRPhiOffsetC);
 264
 265
 266 }
 267
 268 void AliTPCCalibGlobalMisalignment::AddAlign(const  AliTPCCalibGlobalMisalignment & add){
 269   //
 270   // Add the alignmnet to current object
 271   //
 272   fXShift+=add.fXShift;               // Shift in global X [cm]
 273   fYShift+=add.fYShift;               // Shift in global Y [cm]
 274   fZShift+=add.fZShift;               // Shift in global Z [cm]
 275
 276   fRotPhiA+=add.fRotPhiA;      // simple rotation of A side read-out plane around the Z axis [rad]
 277   fRotPhiC+=add.fRotPhiC;      // simple rotation of C side read-out plane around the Z axis [rad]
 278   fdRPhiOffsetA+=add.fdRPhiOffsetA;  // add a constant offset of dRPhi (or local Y) in [cm]: purely for calibration purposes!
 279   fdRPhiOffsetC+=add.fdRPhiOffsetC;  // add a constant offset of dRPhi (or local Y) in [cm]: purely for calibration purposes!
 280   //
 281   // Quadrant alignment
 282   //
 283   if (add.fQuadrantQ0) {
 284     if (fQuadrantQ0) fQuadrantQ0->Add(*(add.fQuadrantQ0));
 285     if (!fQuadrantQ0) fQuadrantQ0 = (TVectorD*)(add.fQuadrantQ0->Clone());
 286   }
 287   if (add.fQuadrantRQ0) {
 288     if (fQuadrantRQ0) fQuadrantRQ0->Add(*(add.fQuadrantRQ0));
 289     if (!fQuadrantRQ0) fQuadrantRQ0 = (TVectorD*)(add.fQuadrantRQ0->Clone());
 290   }
 291   //
 292   if (add.fQuadrantQ1) {
 293     if (fQuadrantQ1) fQuadrantQ1->Add(*(add.fQuadrantQ1));
 294     if (!fQuadrantQ1) fQuadrantQ1 = (TVectorD*)(add.fQuadrantQ1->Clone());
 295   }
 296   if (add.fQuadrantRQ1) {
 297     if (fQuadrantRQ1) fQuadrantRQ1->Add(*(add.fQuadrantRQ1));
 298     if (!fQuadrantRQ1) fQuadrantRQ1 = (TVectorD*)(add.fQuadrantRQ1->Clone());
 299   }
 300   //
 301   if (add.fQuadrantQ2) {
 302     if (fQuadrantQ2) fQuadrantQ2->Add(*(add.fQuadrantQ2));
 303     if (!fQuadrantQ2) fQuadrantQ2 = (TVectorD*)(add.fQuadrantQ2->Clone());
 304   }
 305   if (add.fQuadrantRQ2) {
 306     if (fQuadrantRQ2) fQuadrantRQ2->Add(*(add.fQuadrantRQ2));
 307     if (!fQuadrantRQ2) fQuadrantRQ2 = (TVectorD*)(add.fQuadrantRQ2->Clone());
 308   }
 309   //
 310   // Global alignment - use native ROOT representation
 311   //
 312   if (add.fMatrixGlobal){
 313     if (!fMatrixGlobal)  fMatrixGlobal = new TGeoHMatrix(*(add.fMatrixGlobal));
 314     if (fMatrixGlobal)   ((TGeoHMatrix*)fMatrixGlobal)->Multiply(add.fMatrixGlobal);
 315   }
 316   if (add.fArraySector){
 317     if (!fArraySector) {SetAlignSectors(add.fArraySector);
 318     }else{
 319       for (Int_t isec=0; isec<72; isec++){
 320         TGeoHMatrix *mat0= (TGeoHMatrix*)fArraySector->At(isec);
 321         TGeoHMatrix *mat1= (TGeoHMatrix*)add.fArraySector->At(isec);
 322         if (mat0&&mat1) mat0->Multiply(mat1);
 323       }
 324     }
 325   }
 326 }
 327
 328
 329 AliTPCCalibGlobalMisalignment *  AliTPCCalibGlobalMisalignment::CreateOCDBAlign(){
 330   //
 331   // Create  AliTPCCalibGlobalMisalignment from OCDB Alignment entry
 332   // OCDB has to be initialized before in user code
 333   // All storages (defualt and specific)  and run number
 334   //
 335   AliCDBEntry * entry = AliCDBManager::Instance()->Get("TPC/Align/Data");
 336   if (!entry){
 337     printf("Missing alignmnet entry. OCDB not initialized?\n");
 338     return 0;
 339   }
 340   TClonesArray * array = (TClonesArray*)entry->GetObject();
 341   Int_t entries = array->GetEntries();
 342   TGeoHMatrix matrixGlobal;
 343   TObjArray *alignArrayOCDB= new TObjArray(73);  // sector misalignment + global misalignment
 344   //                                            // global is number 72
 345   //
 346   { for (Int_t i=0;i<entries; i++){
 347       //
 348       //
 349       TGeoHMatrix matrix;
 350       AliAlignObjParams *alignP = (AliAlignObjParams*)array->UncheckedAt(i);
 351       alignP->GetMatrix(matrix);
 352       Int_t imod;
 353       AliGeomManager::ELayerID ilayer;
 354       alignP->GetVolUID(ilayer, imod);
 355       if (ilayer==AliGeomManager::kInvalidLayer) {
 356         alignArrayOCDB->AddAt(matrix.Clone(),72);
 357         alignP->GetMatrix(matrixGlobal);
 358       }else{
 359         Int_t sector=imod;
 360         if (ilayer==AliGeomManager::kTPC2) sector+=36;
 361         alignArrayOCDB->AddAt(matrix.Clone(),sector);
 362       }
 363     }
 364   }
 365   AliTPCCalibGlobalMisalignment *align = new  AliTPCCalibGlobalMisalignment;
 366   align->SetAlignGlobal(&matrixGlobal);
 367   align->SetAlignSectors(alignArrayOCDB);
 368   return align;
 369 }
 370
 371
 372 AliTPCCalibGlobalMisalignment *  AliTPCCalibGlobalMisalignment::CreateMeanAlign(const AliTPCCalibGlobalMisalignment *alignIn){
 373   //
 374   // Create new object, disantangle common mean alignmnet and sector alignment
 375   //
 376   // 1. Try to get mean alignment
 377   // 2. Remove mean alignment from sector alignment
 378   // 3. Create new object
 379   //
 380   TObjArray * array = alignIn->GetAlignSectors();
 381   TObjArray * arrayNew = new TObjArray(72);
 382   //
 383   //Get mean transformation
 384   TGeoHMatrix matrix;
 385   {for (Int_t isec=0; isec<72; isec++){
 386       const TGeoMatrix* cmatrix=(TGeoMatrix*)array->At(isec);
 387       if (!cmatrix) continue;
 388       matrix.Multiply(cmatrix);
 389     }}
 390   TGeoHMatrix matrixMean(matrix);
 391   matrixMean.SetDx(matrix.GetTranslation()[0]/72.);
 392   matrixMean.SetDy(matrix.GetTranslation()[1]/72.);
 393   matrixMean.SetDz(matrix.GetTranslation()[2]/72.);
 394   Double_t rotation[12];
 395   {for (Int_t i=0; i<12; i++) {
 396       rotation[i]=1.0;
 397       if (TMath::Abs(matrix.GetRotationMatrix()[i]-1.)>0.1){
 398       rotation[i]=matrix.GetRotationMatrix()[i]/72.;
 399       }
 400     }}
 401   matrixMean.SetRotation(rotation);
 402   TGeoHMatrix matrixInv = matrixMean.Inverse();
 403   //
 404   {for (Int_t isec=0; isec<72; isec++){
 405       TGeoHMatrix* amatrix=(TGeoHMatrix*)(array->At(isec)->Clone());
 406       if (!amatrix) continue;
 407       amatrix->Multiply(&matrixInv);
 408       arrayNew->AddAt(amatrix,isec);
 409     }}
 410   if (alignIn->GetAlignGlobal()) matrixMean.Multiply((alignIn->GetAlignGlobal()));
 411   AliTPCCalibGlobalMisalignment *alignOut = new  AliTPCCalibGlobalMisalignment;
 412   alignOut->SetAlignGlobal(&matrixMean);
 413   alignOut->SetAlignSectors(arrayNew);
 414   /*
 415     Checks transformation:
 416     AliTPCCalibGlobalMisalignment *  alignIn =  AliTPCCalibGlobalMisalignment::CreateOCDBAlign()
 417     AliTPCCalibGlobalMisalignment * alignOut =  AliTPCCalibGlobalMisalignment::CreateMeanAlign(alignIn)
 418     alignOutM= (AliTPCCalibGlobalMisalignment*)alignOut->Clone();
 419     alignOutS= (AliTPCCalibGlobalMisalignment*)alignOut->Clone();
 420     alignOutS->SetAlignGlobal(0);
 421     alignOutM->SetAlignSectors(0);
 422     //
 423     AliTPCCorrection::AddVisualCorrection(alignOut,0);
 424     AliTPCCorrection::AddVisualCorrection(alignOutM,1);
 425     AliTPCCorrection::AddVisualCorrection(alignOutS,2);
 426     AliTPCCorrection::AddVisualCorrection(alignIn,3);
 427
 428     TF1 f0("f0","AliTPCCorrection::GetCorrSector(x,85,0.9,1,0)",0,18);
 429     TF1 f1("f1","AliTPCCorrection::GetCorrSector(x,85,0.9,1,1)",0,18);
 430     TF1 f2("f2","AliTPCCorrection::GetCorrSector(x,85,0.9,1,2)",0,18);
 431     TF1 f3("f3","AliTPCCorrection::GetCorrSector(x,85,0.9,1,3)",0,18);
 432     f0->SetLineColor(1);
 433     f1->SetLineColor(2);
 434     f2->SetLineColor(3);
 435     f3->SetLineColor(4);
 436     f0->Draw();
 437     f1->Draw("same");
 438     f2->Draw("same");
 439     f3->Draw("same");
 440
 441     TF2 f2D("f2D","AliTPCCorrection::GetCorrSector(x,y,0.9,1,0)-AliTPCCorrection::GetCorrSector(x,y,0.9,1,3)",0,18,85,245);
 442   */
 443   return alignOut;
 444 }
 445
 446
 447 void AliTPCCalibGlobalMisalignment::DumpAlignment( AliTPCCalibGlobalMisalignment* align, TTreeSRedirector *pcstream, const char *name){
 448   //
 449   // Dump alignment per sector into tree
 450   //
 451   TObjArray * array = align->GetAlignSectors();
 452   if (!array) return;
 453   //
 454   //Get mean transformation
 455   TGeoHMatrix matrix;
 456   {for (Int_t isec=0; isec<72; isec++){
 457       TGeoHMatrix* cmatrix=(TGeoHMatrix*)array->At(isec);
 458       TGeoHMatrix* cmatrixDown=(TGeoHMatrix*)array->At(isec%36);
 459       TGeoHMatrix* cmatrixUp=(TGeoHMatrix*)array->At(isec%36+36);
 460       TGeoHMatrix diff(*cmatrixDown);
 461       diff.Multiply(&(cmatrixUp->Inverse()));
 462       (*pcstream)<<name<<
 463         "isec="<<isec<<
 464         "m0.="<<cmatrix<<
 465         "diff.="<<&diff<<
 466         "\n";
 467     }
 468   }
 469
 470 }