EVE/EveDet/AliEveITSDigitsInfo.cxx

   1 // $Id$
   2 // Main authors: Matevz Tadel & Alja Mrak-Tadel: 2006, 2007
   3
   4 /**************************************************************************
   5  * Copyright(c) 1998-2008, ALICE Experiment at CERN, all rights reserved. *
   6  * See http://aliceinfo.cern.ch/Offline/AliRoot/License.html for          *
   7  * full copyright notice.                                                 *
   8  **************************************************************************/
   9
  10 #include "AliEveITSDigitsInfo.h"
  11 #include <EveBase/AliEveEventManager.h>
  12
  13 #include <AliITS.h>
  14 #include <AliITSInitGeometry.h>
  15 #include <AliITSgeomTGeo.h>
  16 #include <AliITSsegmentationSPD.h>
  17 #include <AliITSsegmentationSDD.h>
  18 #include <AliITSsegmentationSSD.h>
  19 #include <AliITSDDLModuleMapSDD.h>
  20
  21 #include <AliITSCalibrationSPD.h>
  22 #include <AliITSCalibrationSDD.h>
  23 #include <AliITSCalibrationSSD.h>
  24 #include <AliITSdigit.h>
  25 #include <AliITSdigitSPD.h>
  26
  27 #include <AliCDBEntry.h>
  28 #include <AliCDBManager.h>
  29
  30 #include <AliRawReader.h>
  31 #include <AliITSRawStreamSPD.h>
  32 #include <AliITSRawStreamSDD.h>
  33 #include <AliITSRawStreamSDDCompressed.h>
  34 #include <AliITSRawStreamSSD.h>
  35
  36 #include <TGeoMatrix.h>
  37 #include <TEveTrans.h>
  38
  39 #include <TMath.h>
  40 #include <TVector3.h>
  41
  42 //==============================================================================
  43 //==============================================================================
  44 // AliEveITSModuleSelection
  45 //==============================================================================
  46
  47 //______________________________________________________________________________
  48 //
  49 // Helper for selecting a range of ITS modules by type, layer, phi and
  50 // theta. Taken as an argument to AliEveITSDigitsInfo::GetModuleIDs().
  51
  52 ClassImp(AliEveITSModuleSelection)
  53
  54 AliEveITSModuleSelection::AliEveITSModuleSelection():
  55   fType(-1),
  56   fLayer(-1),
  57   fMinPhi(-TMath::Pi()),
  58   fMaxPhi(TMath::Pi()),
  59   fMinTheta(-TMath::Pi()),
  60   fMaxTheta(TMath::Pi())
  61 {
  62   // Constructor.
  63 }
  64
  65
  66 //==============================================================================
  67 //==============================================================================
  68 // AliEveITSDigitsInfo
  69 //==============================================================================
  70
  71 //______________________________________________________________________________
  72 //
  73 // Stores ITS geometry information and event-data in format suitable
  74 // for visualization.
  75
  76 ClassImp(AliEveITSDigitsInfo)
  77
  78 AliITSDDLModuleMapSDD *AliEveITSDigitsInfo::fgDDLMapSDD  = 0;
  79 TObjArray             *AliEveITSDigitsInfo::fgDeadModSPD = 0;
  80 TObjArray             *AliEveITSDigitsInfo::fgDeadModSDD = 0;
  81 TObjArray             *AliEveITSDigitsInfo::fgDeadModSSD = 0;
  82
  83 /******************************************************************************/
  84
  85 AliEveITSDigitsInfo::AliEveITSDigitsInfo() :
  86   TObject(),
  87   TEveRefCnt(),
  88   fTree       (0),
  89   fSegSPD     (0), fSegSDD     (0), fSegSSD     (0),
  90   fSPDMinVal  (0), fSSDMinVal  (0), fSDDMinVal  (0),
  91   fSPDMaxVal  (0), fSSDMaxVal  (0), fSDDMaxVal  (0),
  92   fSPDHighLim (0), fSDDHighLim (0), fSSDHighLim (0),
  93   fSPDmap     (),  fSDDmap     (),  fSSDmap     ()
  94 {
  95   // Default constructor.
  96
  97   InitInternals();
  98 }
  99
 100 void AliEveITSDigitsInfo::InitInternals()
 101 {
 102   // Initialize internal geometry structures, in particular the
 103   // module-id to transformation-matrix mapping and segmentation
 104   // classes and data-structures.
 105
 106   static const TEveException kEH("AliEveITSDigitsInfo::InitInternals ");
 107
 108   AliEveEventManager::AssertGeometry();
 109
 110   SetITSSegmentation();
 111
 112   // create tables for scaling
 113   fSPDMinVal = 0;
 114   fSDDMinVal = 5;
 115   fSSDMinVal = 2;
 116
 117   fSPDMaxVal = 1;
 118   fSDDMaxVal = 80;
 119   fSSDMaxVal = 100;
 120
 121   fSPDHighLim = 1;
 122   fSDDHighLim = 512;
 123   fSSDHighLim = 1024;
 124
 125   // lowest scale factor refers to unscaled ITS module
 126   fSPDScaleX[0] = 1;
 127   fSPDScaleZ[0] = 1;
 128   fSDDScaleX[0] = 1;
 129   fSDDScaleZ[0] = 1;
 130   fSSDScale [0] = 1;
 131
 132   // spd lowest resolution
 133   Int_t nx = 8; // fSegSPD->Npx()/8; // 32
 134   Int_t nz = 6; // fSegSPD->Npz()/2; // 128
 135   fSPDScaleX[1] = Int_t(nx);
 136   fSPDScaleZ[1] = Int_t(nz);
 137   fSPDScaleX[2] = Int_t(nx*2);
 138   fSPDScaleZ[2] = Int_t(nz*2);
 139   fSPDScaleX[3] = Int_t(nx*3);
 140   fSPDScaleZ[3] = Int_t(nz*3);
 141   fSPDScaleX[4] = Int_t(nx*4);
 142   fSPDScaleZ[4] = Int_t(nz*4);
 143
 144   fSDDScaleX[1] = 2;
 145   fSDDScaleZ[1] = 2;
 146   fSDDScaleX[2] = 8;
 147   fSDDScaleZ[2] = 8;
 148   fSDDScaleX[3] = 16;
 149   fSDDScaleZ[3] = 16;
 150   fSDDScaleX[4] = 25;
 151   fSDDScaleZ[4] = 25;
 152
 153   fSSDScale[1] = 3;
 154   fSSDScale[2] = 9;
 155   fSSDScale[3] = 20;
 156   fSSDScale[4] = 30;
 157
 158   if (fgDDLMapSDD == 0)
 159   {
 160     fgDDLMapSDD = new AliITSDDLModuleMapSDD();
 161     AliCDBManager *man = AliCDBManager::Instance();
 162     Bool_t cacheStatus = man->GetCacheFlag();
 163     AliCDBEntry   *ddlMapSDD = man->Get("ITS/Calib/DDLMapSDD");
 164
 165     if (!ddlMapSDD)
 166     {
 167       AliWarning("SDD DDL map file retrieval from OCDB failed! - Use default DDL map");
 168     }
 169     else
 170     {
 171       AliITSDDLModuleMapSDD *ddlsdd = (AliITSDDLModuleMapSDD*)ddlMapSDD->GetObject();
 172       if (!ddlsdd)
 173       {
 174         AliWarning("SDD DDL map object not found in OCDB file! - Use default DDL map");
 175       }
 176       else
 177       {
 178         if (!cacheStatus)
 179           ddlMapSDD->SetObject(0);
 180         ddlMapSDD->SetOwner(kTRUE);
 181         fgDDLMapSDD->SetDDLMap(ddlsdd);
 182       }
 183     }
 184
 185     if (!cacheStatus)
 186       delete ddlMapSDD;
 187   }
 188
 189   if (fgDeadModSPD == 0)
 190   {
 191     AliCDBManager *cdb = AliCDBManager::Instance();
 192
 193     AliCDBEntry *deadSPD = cdb->Get("ITS/Calib/SPDDead");
 194
 195     if (!deadSPD)
 196     {
 197       AliWarning("SPD Calibration object retrieval failed!");
 198     }
 199     else
 200     {
 201       fgDeadModSPD = (TObjArray*)deadSPD->GetObject();
 202       fgDeadModSPD->SetOwner(kTRUE);
 203     }
 204   }
 205
 206   // if (fgDeadModSDD == 0)
 207   // {
 208   //   AliCDBManager *cdb = AliCDBManager::Instance();
 209
 210   //   AliCDBEntry *deadSDD = cdb->Get("ITS/Calib/SDDDead");
 211
 212   //   if (!deadSDD)
 213   //   {
 214   //     AliWarning("SDD Calibration object retrieval failed!");
 215   //   }
 216   //   else
 217   //   {
 218   //     fgDeadModSDD = (TObjArray*)deadSDD->GetObject();
 219   //     fgDeadModSDD->SetOwner(kTRUE);
 220   //   }
 221   // }
 222
 223   // if (fgDeadModSSD == 0)
 224   // {
 225   //   AliCDBManager *cdb = AliCDBManager::Instance();
 226
 227   //   AliCDBEntry *deadSSD = cdb->Get("ITS/Calib/SSDDead");
 228
 229   //   if (!deadSSD)
 230   //   {
 231   //     AliWarning("SSD Calibration object retrieval failed!");
 232   //   }
 233   //   else
 234   //   {
 235   //     fgDeadModSSD = (TObjArray*)deadSSD->GetObject();
 236   //     fgDeadModSSD->SetOwner(kTRUE);
 237   //   }
 238   // }
 239 }
 240
 241 /******************************************************************************/
 242
 243 AliEveITSDigitsInfo::~AliEveITSDigitsInfo()
 244 {
 245   // Destructor.
 246   // Deletes the data-maps.
 247   // The tree, if it was passed into the object, is not touched.
 248
 249   std::map<Int_t, TClonesArray*>::iterator j;
 250   for (j = fSPDmap.begin(); j != fSPDmap.end(); ++j)
 251     delete j->second;
 252   for (j = fSDDmap.begin(); j != fSDDmap.end(); ++j)
 253     delete j->second;
 254   for (j = fSSDmap.begin(); j != fSSDmap.end(); ++j)
 255     delete j->second;
 256
 257   delete fSegSPD; delete fSegSDD; delete fSegSSD;
 258 }
 259
 260 /******************************************************************************/
 261
 262 void AliEveITSDigitsInfo::SetTree(TTree* tree)
 263 {
 264   // Set digit-tree to be used for digit retrieval. Data is loaded on
 265   // demand.
 266
 267   fTree = tree;
 268 }
 269
 270 void AliEveITSDigitsInfo::ReadRaw(AliRawReader* raw, Int_t mode)
 271 {
 272   // Read raw-data into internal structures. AliITSdigit is used to
 273   // store raw-adata for all sub-detectors.
 274
 275   if ((mode & 1) || (mode & 2))
 276   {
 277     AliITSRawStreamSPD inputSPD(raw);
 278     TClonesArray* digits = 0;
 279     while (inputSPD.Next())
 280     {
 281       Int_t module = inputSPD.GetModuleID();
 282       Int_t column = inputSPD.GetColumn();
 283       Int_t row    = inputSPD.GetRow();
 284
 285       if (inputSPD.IsNewModule())
 286       {
 287         digits = fSPDmap[module];
 288         if (digits == 0)
 289           fSPDmap[module] = digits = new TClonesArray("AliITSdigit", 16);
 290       } else if(!digits) {
 291                   AliFatal("No module found\n");
 292           }
 293
 294       AliITSdigit* d = new ((*digits)[digits->GetEntriesFast()]) AliITSdigit();
 295       d->SetCoord1(column);
 296       d->SetCoord2(row);
 297       d->SetSignal(1);
 298
 299       // printf("SPD: %d %d %d\n",module,column,row);
 300     }
 301     raw->Reset();
 302   }
 303
 304   if ((mode & 4) || (mode & 8))
 305   {
 306     AliITSRawStream* inputSDD = AliITSRawStreamSDD::CreateRawStreamSDD(raw);
 307     inputSDD->SetDDLModuleMap(fgDDLMapSDD);
 308     TClonesArray* digits = 0;
 309     while (inputSDD->Next())
 310     {
 311       Int_t module = inputSDD->GetModuleID();
 312
 313       digits = fSDDmap[module];
 314       if (digits == 0)
 315         fSDDmap[module] = digits = new TClonesArray("AliITSdigit", 0);
 316
 317       if (inputSDD->IsCompletedModule()==kFALSE && inputSDD->IsCompletedDDL()==kFALSE)
 318       {
 319         Int_t anode  = inputSDD->GetCoord1()+inputSDD->GetChannel()*AliITSsegmentationSDD::GetNAnodesPerHybrid();
 320         Int_t time   = inputSDD->GetCoord2();
 321         Int_t signal = inputSDD->GetSignal();
 322         AliITSdigit* d = new ((*digits)[digits->GetEntriesFast()]) AliITSdigit();
 323         d->SetCoord1(anode);
 324         d->SetCoord2(time);
 325         d->SetSignal(signal);
 326       }
 327       // printf("SDD: %d %d %d %d\n",module,anode,time,signal);
 328     }
 329     delete inputSDD;
 330     raw->Reset();
 331   }
 332
 333   if ((mode & 16) || (mode & 32))
 334   {
 335     AliITSRawStreamSSD input(raw);
 336     TClonesArray* digits = 0;
 337     while (input.Next())
 338     {
 339       Int_t module  = input.GetModuleID();
 340       Int_t side    = input.GetSideFlag();
 341       Int_t strip   = input.GetStrip();
 342       Int_t signal  = input.GetSignal();
 343
 344       // !!!! IsNewModule() is false in the beginning of the stream, so we also
 345       // !!!! check for digits == 0. Should be fixed in SSD stream.
 346       if (input.IsNewModule() || digits == 0)
 347       {
 348         digits = fSSDmap[module];
 349         if (digits == 0)
 350           fSSDmap[module] = digits = new TClonesArray("AliITSdigit", 0);
 351       }
 352
 353       AliITSdigit* d = new ((*digits)[digits->GetEntriesFast()]) AliITSdigit();
 354       d->SetCoord1(side);
 355       d->SetCoord2(strip);
 356       d->SetSignal(signal);
 357
 358       // printf("SSD: %d %d %d %d\n",module,side,strip,signal);
 359     }
 360     raw->Reset();
 361   }
 362 }
 363
 364 /******************************************************************************/
 365
 366 void AliEveITSDigitsInfo::SetITSSegmentation()
 367 {
 368   // Create the segmentation objects and fill internal
 369   // data-structures.
 370
 371   // SPD
 372   fSegSPD = new AliITSsegmentationSPD("TGeo");
 373
 374   Float_t fNzSPD = 160;
 375   Float_t fHlSPD = 3.48;
 376   Float_t fZ1pitchSPD = 0.0425, fZ2pitchSPD = 0.0625;
 377
 378   fSPDZCoord[0] = fZ1pitchSPD - fHlSPD;
 379   for (Int_t m=1; m<fNzSPD; m++)
 380   {
 381     Float_t dz = fZ1pitchSPD;
 382     if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || m==127 || m==128)
 383     {
 384       dz = fZ2pitchSPD;
 385     }
 386     fSPDZCoord[m] = fSPDZCoord[m-1] + dz;
 387   }
 388
 389   for (Int_t m=0; m<fNzSPD; m++)
 390   {
 391     Float_t dz = fZ1pitchSPD;
 392     if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || m==127 || m==128)
 393     {
 394       dz = fZ2pitchSPD;
 395     }
 396     fSPDZCoord[m] -= dz;
 397   }
 398
 399   // SDD
 400   fSegSDD = new AliITSsegmentationSDD("TGeo");
 401   // !!!! Set default drift speed, eventually need to get it from CDB.
 402   fSegSDD->SetDriftSpeed(7.3);
 403
 404   // SSD
 405   fSegSSD = new AliITSsegmentationSSD("TGeo");
 406 }
 407
 408 /******************************************************************************/
 409
 410 TClonesArray* AliEveITSDigitsInfo::GetDigits(Int_t mod, Int_t subdet) const
 411 {
 412   // Return TClonesArray of digits for specified module and sub-detector-id.
 413
 414   switch(subdet)
 415   {
 416     case 0:
 417     {
 418       std::map<Int_t, TClonesArray*>::iterator i = fSPDmap.find(mod);
 419       if (i != fSPDmap.end())
 420       {
 421         return i->second;
 422       }
 423       else if (fTree)
 424       {
 425         TClonesArray *digitsSPD = 0;
 426         TBranch *br = fTree->GetBranch("ITSDigitsSPD");
 427         br->SetAddress(&digitsSPD);
 428         br->GetEntry(mod);
 429         fSPDmap[mod] = digitsSPD;
 430         return digitsSPD;
 431       }
 432       break;
 433     }
 434
 435     case 1:
 436     {
 437       std::map<Int_t, TClonesArray*>::iterator i = fSDDmap.find(mod);
 438       if (i != fSDDmap.end())
 439       {
 440         return i->second;
 441       }
 442       else if (fTree)
 443       {
 444         TClonesArray *digitsSDD = 0;
 445         TBranch *br = fTree->GetBranch("ITSDigitsSDD");
 446         br->SetAddress(&digitsSDD);
 447         br->GetEntry(mod);
 448         fSDDmap[mod] = digitsSDD;
 449         return digitsSDD;
 450       }
 451       break;
 452     }
 453
 454     case 2:
 455     {
 456       std::map<Int_t, TClonesArray*>::iterator i = fSSDmap.find(mod);
 457       if (i != fSSDmap.end())
 458       {
 459         return i->second;
 460       }
 461       else if (fTree)
 462       {
 463         TClonesArray *digitsSSD = 0;
 464         TBranch *br = fTree->GetBranch("ITSDigitsSSD");
 465         br->SetAddress(&digitsSSD);
 466         br->GetEntry(mod);
 467         fSSDmap[mod] = digitsSSD;
 468         return digitsSSD;
 469       }
 470       break;
 471     }
 472   }
 473
 474   return 0;
 475 }
 476
 477 /******************************************************************************/
 478
 479 void AliEveITSDigitsInfo::GetModuleIDs(AliEveITSModuleSelection* sel,
 480                                        std::vector<UInt_t>& ids) const
 481 {
 482   // Fill the id-vector with ids of modules that satisfy conditions
 483   // given by the AliEveITSModuleSelection object.
 484
 485   Int_t idx0 = 0, idx1 = 0;
 486   switch(sel->GetType())
 487   {
 488     case 0:
 489       idx0 = 0;
 490       idx1 = AliITSgeomTGeo::GetModuleIndex(3, 1, 1) - 1;
 491       break;
 492     case 1:
 493       idx0 = AliITSgeomTGeo::GetModuleIndex(3, 1, 1);
 494       idx1 = AliITSgeomTGeo::GetModuleIndex(5, 1, 1) - 1;
 495       break;
 496     case 2:
 497       idx0 = AliITSgeomTGeo::GetModuleIndex(5, 1, 1);
 498       idx1 = AliITSgeomTGeo::GetNModules() - 1;
 499       break;
 500     default:
 501       idx1 = 0;
 502       idx1 = AliITSgeomTGeo::GetNModules() - 1;
 503       break;
 504   }
 505
 506   TVector3  v;
 507   TEveTrans mx;
 508   Int_t     lay, lad, det;
 509   for (Int_t id = idx0; id < idx1; ++id)
 510   {
 511     AliITSgeomTGeo::GetModuleId(id, lay, lad, det);
 512     if (sel->GetLayer() == lay || sel->GetLayer() == -1)
 513     {
 514       // check data from matrix
 515       mx.SetFrom(*AliITSgeomTGeo::GetMatrix(id));
 516       mx.GetPos(v);
 517       if (v.Phi()   <= sel->GetMaxPhi()   && v.Phi()   >= sel->GetMinPhi()   &&
 518           v.Theta() <= sel->GetMaxTheta() && v.Theta() >= sel->GetMinTheta())
 519       {
 520         ids.push_back(id);
 521       }
 522     }
 523   }
 524 }
 525
 526 /******************************************************************************/
 527
 528 Bool_t AliEveITSDigitsInfo::HasData(Int_t module, Int_t det_id) const
 529 {
 530   // Return true if given module has data.
 531
 532   TClonesArray *digits = GetDigits(module, det_id);
 533   return digits && digits->GetEntriesFast() > 0;
 534 }
 535
 536 Bool_t AliEveITSDigitsInfo::IsDead (Int_t module, Int_t det_id) const
 537 {
 538   // Return true if given module is dead.
 539   // Only implemented for SPD.
 540
 541   if (det_id == 0 && fgDeadModSPD)
 542     return ((AliITSCalibrationSPD*) fgDeadModSPD->At(module))->IsBad();
 543   if (det_id == 1 && fgDeadModSDD)
 544     return ((AliITSCalibrationSDD*) fgDeadModSDD->At(module))->IsBad();
 545   if (det_id == 2 && fgDeadModSSD)
 546     return ((AliITSCalibrationSSD*) fgDeadModSSD->At(module))->IsBad();
 547   return kFALSE;
 548 }
 549
 550 /******************************************************************************/
 551
 552 void AliEveITSDigitsInfo::Print(Option_t* ) const
 553 {
 554   // Print information about stored geometry and segmentation.
 555
 556   printf("*********************************************************\n");
 557   printf("SPD module dimension (%f,%f)\n",           fSegSPD->Dx()*0.0001, fSegSPD->Dz()*0.0001);
 558   printf("SPD first,last module: %d, %d\n",
 559          AliITSgeomTGeo::GetModuleIndex(1,1,1),
 560          AliITSgeomTGeo::GetModuleIndex(3,1,1) - 1);
 561   printf("SPD num cells per module (x::%d,z::%d)\n", fSegSPD->Npx(), fSegSPD->Npz());
 562   Int_t iz = 0, ix = 0;
 563   printf("SPD dimesion of (%d,%d) in pixel(%f,%f)\n",   ix, iz, fSegSPD->Dpx(ix), fSegSPD->Dpz(iz));
 564   iz = 32;
 565   printf("SPD dimesion of pixel (%d,%d) are (%f,%f)\n", ix, iz, fSegSPD->Dpx(ix)*0.001, fSegSPD->Dpz(iz)*0.001);
 566
 567   printf("*********************************************************\n");
 568   printf("SDD module dimension (%f,%f)\n",           fSegSDD->Dx()*0.0001, fSegSDD->Dz()*0.0001);
 569   printf("SDD first,last module: %d, %d\n",
 570          AliITSgeomTGeo::GetModuleIndex(3,1,1),
 571          AliITSgeomTGeo::GetModuleIndex(5,1,1) - 1);
 572   printf("SDD num cells per module (x::%d,z::%d)\n", fSegSDD->Npx(), fSegSDD->Npz());
 573   printf("SDD dimesion of pixel are (%f,%f)\n",      fSegSDD->Dpx(1)*0.001, fSegSDD->Dpz(1)*0.001);
 574
 575   Float_t ap, an;
 576   printf("*********************************************************\n");
 577   printf("SSD module dimension (%f,%f)\n",  fSegSSD->Dx()*0.0001, fSegSSD->Dz()*0.0001);
 578   printf("SSD first, last module: %d, %d\n",
 579          AliITSgeomTGeo::GetModuleIndex(5,1,1),
 580          AliITSgeomTGeo::GetNModules() - 1);
 581   printf("SSD strips in module %d\n",       fSegSSD->Npx());
 582   printf("SSD strip sizes are (%f,%f)\n",   fSegSSD->Dpx(1), fSegSSD->Dpz(1));
 583   fSegSSD->SetLayer(5);  fSegSSD->Angles(ap,an);
 584   printf("SSD layer 5 stereoP %f stereoN %f angle\n", ap, an);
 585   fSegSSD->SetLayer(6);  fSegSSD->Angles(ap,an);
 586   printf("SSD layer 6 stereoP %f stereoN %f angle\n", ap, an);
 587 }
 588
 589
 590 /*
 591   printf("num cells %d,%d scaled %d,%d \n",fSegSPD->Npz(),fSegSPD->Npx(),Nz,Nx);
 592   printf("%d digits in AliEveITSModule %d\n",ne, module);
 593   Float_t zn = i*(3.48*2)/Nz - 3.48 ;
 594   Float_t xo =  -fSegSPD->Dx()*0.00005 + fSegSPD->Dpx(0)*od->GetCoord2()*0.0001;
 595   Float_t xn =  -fSegSPD->Dx()*0.00005 + j*0.0001*fSegSPD->Dx()/Nx;
 596   Float_t dpx = 0.0001*fSegSPD->Dx()/Nx;
 597   Float_t dpz = 3.48*2/Nz;
 598   printf("Z::original (%3f) scaled (%3f, %3f) \n", zo, zn-dpz/2, zn+dpz/2);
 599   printf("X::original (%3f) scaled (%3f, %3f) \n", xo, xn-dpx/2, xn+dpx/2);
 600   printf("%d,%d maped to %d,%d \n", od->GetCoord1(), od->GetCoord2(), i,j );
 601 */