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       }
 291
 292       AliITSdigit* d = new ((*digits)[digits->GetEntriesFast()]) AliITSdigit();
 293       d->SetCoord1(column);
 294       d->SetCoord2(row);
 295       d->SetSignal(1);
 296
 297       // printf("SPD: %d %d %d\n",module,column,row);
 298     }
 299     raw->Reset();
 300   }
 301
 302   if ((mode & 4) || (mode & 8))
 303   {
 304     AliITSRawStream* inputSDD = AliITSRawStreamSDD::CreateRawStreamSDD(raw);
 305     inputSDD->SetDDLModuleMap(fgDDLMapSDD);
 306     TClonesArray* digits = 0;
 307     while (inputSDD->Next())
 308     {
 309       Int_t module = inputSDD->GetModuleID();
 310
 311       digits = fSDDmap[module];
 312       if (digits == 0)
 313         fSDDmap[module] = digits = new TClonesArray("AliITSdigit", 0);
 314
 315       if (inputSDD->IsCompletedModule()==kFALSE && inputSDD->IsCompletedDDL()==kFALSE)
 316       {
 317         Int_t anode  = inputSDD->GetCoord1()+inputSDD->GetChannel()*AliITSsegmentationSDD::GetNAnodesPerHybrid();
 318         Int_t time   = inputSDD->GetCoord2();
 319         Int_t signal = inputSDD->GetSignal();
 320         AliITSdigit* d = new ((*digits)[digits->GetEntriesFast()]) AliITSdigit();
 321         d->SetCoord1(anode);
 322         d->SetCoord2(time);
 323         d->SetSignal(signal);
 324       }
 325       // printf("SDD: %d %d %d %d\n",module,anode,time,signal);
 326     }
 327     delete inputSDD;
 328     raw->Reset();
 329   }
 330
 331   if ((mode & 16) || (mode & 32))
 332   {
 333     AliITSRawStreamSSD input(raw);
 334     TClonesArray* digits = 0;
 335     while (input.Next())
 336     {
 337       Int_t module  = input.GetModuleID();
 338       Int_t side    = input.GetSideFlag();
 339       Int_t strip   = input.GetStrip();
 340       Int_t signal  = input.GetSignal();
 341
 342       // !!!! IsNewModule() is false in the beginning of the stream, so we also
 343       // !!!! check for digits == 0. Should be fixed in SSD stream.
 344       if (input.IsNewModule() || digits == 0)
 345       {
 346         digits = fSSDmap[module];
 347         if (digits == 0)
 348           fSSDmap[module] = digits = new TClonesArray("AliITSdigit", 0);
 349       }
 350
 351       AliITSdigit* d = new ((*digits)[digits->GetEntriesFast()]) AliITSdigit();
 352       d->SetCoord1(side);
 353       d->SetCoord2(strip);
 354       d->SetSignal(signal);
 355
 356       // printf("SSD: %d %d %d %d\n",module,side,strip,signal);
 357     }
 358     raw->Reset();
 359   }
 360 }
 361
 362 /******************************************************************************/
 363
 364 void AliEveITSDigitsInfo::SetITSSegmentation()
 365 {
 366   // Create the segmentation objects and fill internal
 367   // data-structures.
 368
 369   // SPD
 370   fSegSPD = new AliITSsegmentationSPD("TGeo");
 371
 372   Float_t fNzSPD = 160;
 373   Float_t fHlSPD = 3.48;
 374   Float_t fZ1pitchSPD = 0.0425, fZ2pitchSPD = 0.0625;
 375
 376   fSPDZCoord[0] = fZ1pitchSPD - fHlSPD;
 377   for (Int_t m=1; m<fNzSPD; m++)
 378   {
 379     Float_t dz = fZ1pitchSPD;
 380     if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || m==127 || m==128)
 381     {
 382       dz = fZ2pitchSPD;
 383     }
 384     fSPDZCoord[m] = fSPDZCoord[m-1] + dz;
 385   }
 386
 387   for (Int_t m=0; m<fNzSPD; m++)
 388   {
 389     Float_t dz = fZ1pitchSPD;
 390     if (m==31 || m==32 || m==63 || m==64 || m==95 || m==96 || m==127 || m==128)
 391     {
 392       dz = fZ2pitchSPD;
 393     }
 394     fSPDZCoord[m] -= dz;
 395   }
 396
 397   // SDD
 398   fSegSDD = new AliITSsegmentationSDD("TGeo");
 399   // !!!! Set default drift speed, eventually need to get it from CDB.
 400   fSegSDD->SetDriftSpeed(7.3);
 401
 402   // SSD
 403   fSegSSD = new AliITSsegmentationSSD("TGeo");
 404 }
 405
 406 /******************************************************************************/
 407
 408 TClonesArray* AliEveITSDigitsInfo::GetDigits(Int_t mod, Int_t subdet) const
 409 {
 410   // Return TClonesArray of digits for specified module and sub-detector-id.
 411
 412   switch(subdet)
 413   {
 414     case 0:
 415     {
 416       std::map<Int_t, TClonesArray*>::iterator i = fSPDmap.find(mod);
 417       if (i != fSPDmap.end())
 418       {
 419         return i->second;
 420       }
 421       else if (fTree)
 422       {
 423         TClonesArray *digitsSPD = 0;
 424         TBranch *br = fTree->GetBranch("ITSDigitsSPD");
 425         br->SetAddress(&digitsSPD);
 426         br->GetEntry(mod);
 427         fSPDmap[mod] = digitsSPD;
 428         return digitsSPD;
 429       }
 430       break;
 431     }
 432
 433     case 1:
 434     {
 435       std::map<Int_t, TClonesArray*>::iterator i = fSDDmap.find(mod);
 436       if (i != fSDDmap.end())
 437       {
 438         return i->second;
 439       }
 440       else if (fTree)
 441       {
 442         TClonesArray *digitsSDD = 0;
 443         TBranch *br = fTree->GetBranch("ITSDigitsSDD");
 444         br->SetAddress(&digitsSDD);
 445         br->GetEntry(mod);
 446         fSDDmap[mod] = digitsSDD;
 447         return digitsSDD;
 448       }
 449       break;
 450     }
 451
 452     case 2:
 453     {
 454       std::map<Int_t, TClonesArray*>::iterator i = fSSDmap.find(mod);
 455       if (i != fSSDmap.end())
 456       {
 457         return i->second;
 458       }
 459       else if (fTree)
 460       {
 461         TClonesArray *digitsSSD = 0;
 462         TBranch *br = fTree->GetBranch("ITSDigitsSSD");
 463         br->SetAddress(&digitsSSD);
 464         br->GetEntry(mod);
 465         fSSDmap[mod] = digitsSSD;
 466         return digitsSSD;
 467       }
 468       break;
 469     }
 470   }
 471
 472   return 0;
 473 }
 474
 475 /******************************************************************************/
 476
 477 void AliEveITSDigitsInfo::GetModuleIDs(AliEveITSModuleSelection* sel,
 478                                        std::vector<UInt_t>& ids) const
 479 {
 480   // Fill the id-vector with ids of modules that satisfy conditions
 481   // given by the AliEveITSModuleSelection object.
 482
 483   Int_t idx0 = 0, idx1 = 0;
 484   switch(sel->GetType())
 485   {
 486     case 0:
 487       idx0 = 0;
 488       idx1 = AliITSgeomTGeo::GetModuleIndex(3, 1, 1) - 1;
 489       break;
 490     case 1:
 491       idx0 = AliITSgeomTGeo::GetModuleIndex(3, 1, 1);
 492       idx1 = AliITSgeomTGeo::GetModuleIndex(5, 1, 1) - 1;
 493       break;
 494     case 2:
 495       idx0 = AliITSgeomTGeo::GetModuleIndex(5, 1, 1);
 496       idx1 = AliITSgeomTGeo::GetNModules() - 1;
 497       break;
 498     default:
 499       idx1 = 0;
 500       idx1 = AliITSgeomTGeo::GetNModules() - 1;
 501       break;
 502   }
 503
 504   TVector3  v;
 505   TEveTrans mx;
 506   Int_t     lay, lad, det;
 507   for (Int_t id = idx0; id < idx1; ++id)
 508   {
 509     AliITSgeomTGeo::GetModuleId(id, lay, lad, det);
 510     if (sel->GetLayer() == lay || sel->GetLayer() == -1)
 511     {
 512       // check data from matrix
 513       mx.SetFrom(*AliITSgeomTGeo::GetMatrix(id));
 514       mx.GetPos(v);
 515       if (v.Phi()   <= sel->GetMaxPhi()   && v.Phi()   >= sel->GetMinPhi()   &&
 516           v.Theta() <= sel->GetMaxTheta() && v.Theta() >= sel->GetMinTheta())
 517       {
 518         ids.push_back(id);
 519       }
 520     }
 521   }
 522 }
 523
 524 /******************************************************************************/
 525
 526 Bool_t AliEveITSDigitsInfo::HasData(Int_t module, Int_t det_id) const
 527 {
 528   // Return true if given module has data.
 529
 530   TClonesArray *digits = GetDigits(module, det_id);
 531   return digits && digits->GetEntriesFast() > 0;
 532 }
 533
 534 Bool_t AliEveITSDigitsInfo::IsDead (Int_t module, Int_t det_id) const
 535 {
 536   // Return true if given module is dead.
 537   // Only implemented for SPD.
 538
 539   if (det_id == 0 && fgDeadModSPD)
 540     return ((AliITSCalibrationSPD*) fgDeadModSPD->At(module))->IsBad();
 541   if (det_id == 1 && fgDeadModSDD)
 542     return ((AliITSCalibrationSDD*) fgDeadModSDD->At(module))->IsBad();
 543   if (det_id == 2 && fgDeadModSSD)
 544     return ((AliITSCalibrationSSD*) fgDeadModSSD->At(module))->IsBad();
 545   return kFALSE;
 546 }
 547
 548 /******************************************************************************/
 549
 550 void AliEveITSDigitsInfo::Print(Option_t* ) const
 551 {
 552   // Print information about stored geometry and segmentation.
 553
 554   printf("*********************************************************\n");
 555   printf("SPD module dimension (%f,%f)\n",           fSegSPD->Dx()*0.0001, fSegSPD->Dz()*0.0001);
 556   printf("SPD first,last module: %d, %d\n",
 557          AliITSgeomTGeo::GetModuleIndex(1,1,1),
 558          AliITSgeomTGeo::GetModuleIndex(3,1,1) - 1);
 559   printf("SPD num cells per module (x::%d,z::%d)\n", fSegSPD->Npx(), fSegSPD->Npz());
 560   Int_t iz = 0, ix = 0;
 561   printf("SPD dimesion of (%d,%d) in pixel(%f,%f)\n",   ix, iz, fSegSPD->Dpx(ix), fSegSPD->Dpz(iz));
 562   iz = 32;
 563   printf("SPD dimesion of pixel (%d,%d) are (%f,%f)\n", ix, iz, fSegSPD->Dpx(ix)*0.001, fSegSPD->Dpz(iz)*0.001);
 564
 565   printf("*********************************************************\n");
 566   printf("SDD module dimension (%f,%f)\n",           fSegSDD->Dx()*0.0001, fSegSDD->Dz()*0.0001);
 567   printf("SDD first,last module: %d, %d\n",
 568          AliITSgeomTGeo::GetModuleIndex(3,1,1),
 569          AliITSgeomTGeo::GetModuleIndex(5,1,1) - 1);
 570   printf("SDD num cells per module (x::%d,z::%d)\n", fSegSDD->Npx(), fSegSDD->Npz());
 571   printf("SDD dimesion of pixel are (%f,%f)\n",      fSegSDD->Dpx(1)*0.001, fSegSDD->Dpz(1)*0.001);
 572
 573   Float_t ap, an;
 574   printf("*********************************************************\n");
 575   printf("SSD module dimension (%f,%f)\n",  fSegSSD->Dx()*0.0001, fSegSSD->Dz()*0.0001);
 576   printf("SSD first, last module: %d, %d\n",
 577          AliITSgeomTGeo::GetModuleIndex(5,1,1),
 578          AliITSgeomTGeo::GetNModules() - 1);
 579   printf("SSD strips in module %d\n",       fSegSSD->Npx());
 580   printf("SSD strip sizes are (%f,%f)\n",   fSegSSD->Dpx(1), fSegSSD->Dpz(1));
 581   fSegSSD->SetLayer(5);  fSegSSD->Angles(ap,an);
 582   printf("SSD layer 5 stereoP %f stereoN %f angle\n", ap, an);
 583   fSegSSD->SetLayer(6);  fSegSSD->Angles(ap,an);
 584   printf("SSD layer 6 stereoP %f stereoN %f angle\n", ap, an);
 585 }
 586
 587
 588 /*
 589   printf("num cells %d,%d scaled %d,%d \n",fSegSPD->Npz(),fSegSPD->Npx(),Nz,Nx);
 590   printf("%d digits in AliEveITSModule %d\n",ne, module);
 591   Float_t zn = i*(3.48*2)/Nz - 3.48 ;
 592   Float_t xo =  -fSegSPD->Dx()*0.00005 + fSegSPD->Dpx(0)*od->GetCoord2()*0.0001;
 593   Float_t xn =  -fSegSPD->Dx()*0.00005 + j*0.0001*fSegSPD->Dx()/Nx;
 594   Float_t dpx = 0.0001*fSegSPD->Dx()/Nx;
 595   Float_t dpz = 3.48*2/Nz;
 596   printf("Z::original (%3f) scaled (%3f, %3f) \n", zo, zn-dpz/2, zn+dpz/2);
 597   printf("X::original (%3f) scaled (%3f, %3f) \n", xo, xn-dpx/2, xn+dpx/2);
 598   printf("%d,%d maped to %d,%d \n", od->GetCoord1(), od->GetCoord2(), i,j );
 599 */