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