ITS/UPGRADE/AliITSUGeomTGeo.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 //    AliITSUGeomTGeo is a simple interface class to TGeoManager       //
  18 //    It is used in the simulation and reconstruction in order to        //
  19 //    query the TGeo ITS geometry                                        //
  20 //                                                                       //
  21 //    author - cvetan.cheshkov@cern.ch                                   //
  22 //    15/02/2007                                                         //
  23 //    adapted to ITSupg 18/07/2012 - ruben.shahoyan@cern.ch              //
  24 //                                                                       //
  25 //    ATTENTION: In opposite to ols AliITSgeomTGeo, all indices start    //
  26 //    from 0, not from 1!!!                                              //
  27 //                                                                       //
  28 ///////////////////////////////////////////////////////////////////////////
  29
  30 #include <TClass.h>
  31 #include <TString.h>
  32 #include <TGeoManager.h>
  33 #include <TGeoPhysicalNode.h>
  34 #include <TGeoShape.h>
  35 #include <TGeoBBox.h>
  36 #include <TDatime.h>
  37 #include <TMath.h>
  38 #include <TSystem.h>
  39
  40 #include "AliITSUGeomTGeo.h"
  41 #include "AliLog.h"
  42 #include "AliAlignObj.h"
  43 #include "AliITSsegmentation.h"
  44 #include "AliITSUSegmentationPix.h"
  45 using namespace TMath;
  46
  47 ClassImp(AliITSUGeomTGeo)
  48
  49 UInt_t AliITSUGeomTGeo::fgUIDShift = 16;                // bit shift to go from mod.id to modUUID for TGeo
  50 TString AliITSUGeomTGeo::fgITSVolName = "ITSV";
  51 TString AliITSUGeomTGeo::fgITSLrName  = "ITSULayer";
  52 TString AliITSUGeomTGeo::fgITSStaveName = "ITSUStave";
  53 TString AliITSUGeomTGeo::fgITSHalfStaveName = "ITSUHalfStave";
  54 TString AliITSUGeomTGeo::fgITSModuleName = "ITSUModule";
  55 TString AliITSUGeomTGeo::fgITSChipName = "ITSUChip";
  56 TString AliITSUGeomTGeo::fgITSSensName = "ITSUSensor";
  57 TString AliITSUGeomTGeo::fgITSWrapVolName = "ITSUWrapVol";
  58 TString AliITSUGeomTGeo::fgITSChipTypeName[AliITSUGeomTGeo::kNChipTypes] = {"Pix"};
  59 //
  60 TString AliITSUGeomTGeo::fgITSsegmFileName = "itsSegmentations.root";
  61
  62 //______________________________________________________________________
  63 AliITSUGeomTGeo::AliITSUGeomTGeo(Bool_t build, Bool_t loadSegm)
  64   :fVersion(kITSVNA)
  65   ,fNLayers(0)
  66   ,fNChips(0)
  67   ,fNStaves(0)
  68   ,fNHalfStaves(0)
  69   ,fNModules(0)
  70   ,fNChipsPerModule(0)
  71   ,fNChipRowsPerModule(0)
  72   ,fNChipsPerHalfStave(0)
  73   ,fNChipsPerStave(0)
  74   ,fNChipsPerLayer(0)
  75   ,fLrChipType(0)
  76   ,fLastChipIndex(0)
  77   ,fMatSens(0)
  78   ,fMatT2L(0)
  79   ,fSegm(0)
  80 {
  81   // default c-tor
  82   for (int i=AliITSUAux::kMaxLayers;i--;) fLr2Wrapper[i] = -1;
  83   if (build) BuildITS(loadSegm);
  84 }
  85
  86 //______________________________________________________________________
  87 AliITSUGeomTGeo::AliITSUGeomTGeo(const AliITSUGeomTGeo &src)
  88   :TObject(src)
  89   ,fVersion(src.fVersion)
  90   ,fNLayers(src.fNLayers)
  91   ,fNChips(src.fNChips)
  92   ,fNStaves(0)
  93   ,fNHalfStaves(0)
  94   ,fNModules(0)
  95   ,fNChipsPerModule(0)
  96   ,fNChipRowsPerModule(0)
  97   ,fNChipsPerHalfStave(0)
  98   ,fNChipsPerStave(0)
  99   ,fNChipsPerLayer(0)
 100   ,fLrChipType(0)
 101   ,fLastChipIndex(0)
 102   ,fMatSens(0)
 103   ,fMatT2L(0)
 104   ,fSegm(0)
 105 {
 106   // copy c-tor
 107   if (fNLayers) {
 108     fNStaves   = new Int_t[fNLayers];
 109     fNChipsPerModule = new Int_t[fNLayers];
 110     fNChipRowsPerModule = new Int_t[fNLayers];
 111     fLrChipType  = new Int_t[fNLayers];
 112     fLastChipIndex   = new Int_t[fNLayers];
 113     fNChipsPerHalfStave = new Int_t[fNLayers];
 114     fNChipsPerStave = new Int_t[fNLayers];
 115     fNChipsPerLayer = new Int_t[fNLayers];
 116     //
 117     for (int i=fNLayers;i--;) {
 118       fNStaves[i] = src.fNStaves[i];
 119       fNHalfStaves[i] = src.fNHalfStaves[i];
 120       fNModules[i] = src.fNModules[i];
 121       fNChipsPerModule[i] = src.fNChipsPerModule[i];
 122       fNChipRowsPerModule[i] = src.fNChipRowsPerModule[i];
 123       fNChipsPerHalfStave[i] = src.fNChipsPerHalfStave[i];
 124       fNChipsPerStave[i] = src.fNChipsPerStave[i];
 125       fNChipsPerLayer[i] = src.fNChipsPerLayer[i];
 126       fLrChipType[i]  = src.fLrChipType[i];
 127       fLastChipIndex[i] = src.fLastChipIndex[i];
 128     }
 129     if (src.fMatSens) {
 130       fMatSens = new TObjArray(fNChips);
 131       fMatSens->SetOwner(kTRUE);
 132       for (int i=0;i<fNChips;i++) {
 133         const TGeoHMatrix* mat = (TGeoHMatrix*)src.fMatSens->At(i);
 134         fMatSens->AddAt(new TGeoHMatrix(*mat),i);
 135       }
 136     }
 137     if (src.fMatT2L) {
 138       fMatT2L = new TObjArray(fNChips);
 139       fMatT2L->SetOwner(kTRUE);
 140       for (int i=0;i<fNChips;i++) {
 141         const TGeoHMatrix* mat =(TGeoHMatrix*) src.fMatT2L->At(i);
 142         fMatT2L->AddAt(new TGeoHMatrix(*mat),i);
 143       }
 144     }
 145     if (src.fSegm) {
 146       int sz = src.fSegm->GetEntriesFast();
 147       fSegm = new TObjArray(sz);
 148       fSegm->SetOwner(kTRUE);
 149       for (int i=0;i<sz;i++) {
 150         AliITSsegmentation* sg = (AliITSsegmentation*)src.fSegm->UncheckedAt(i);
 151         if (!sg) continue;
 152         fSegm->AddAt(sg->Clone(),i);
 153       }
 154     }
 155   }
 156   for (int i=AliITSUAux::kMaxLayers;i--;) fLr2Wrapper[i] = src.fLr2Wrapper[i];
 157 }
 158
 159 //______________________________________________________________________
 160 AliITSUGeomTGeo::~AliITSUGeomTGeo()
 161 {
 162   //d-tor
 163   delete[] fNStaves;
 164   delete[] fNHalfStaves;
 165   delete[] fNModules;
 166   delete[] fLrChipType;
 167   delete[] fNChipsPerModule;
 168   delete[] fNChipRowsPerModule;
 169   delete[] fNChipsPerHalfStave;
 170   delete[] fNChipsPerStave;
 171   delete[] fNChipsPerLayer;
 172   delete[] fLastChipIndex;
 173   delete fMatT2L;
 174   delete fMatSens;
 175   delete fSegm;
 176 }
 177
 178
 179 //______________________________________________________________________
 180 AliITSUGeomTGeo& AliITSUGeomTGeo::operator=(const AliITSUGeomTGeo &src)
 181 {
 182   // cp op.
 183   if (this!=&src) {
 184     delete[] fNStaves;
 185     delete[] fNHalfStaves;
 186     delete[] fNModules;
 187     delete[] fLrChipType;
 188     delete[] fNChipsPerModule;
 189     delete[] fNChipRowsPerModule;
 190     delete[] fNChipsPerHalfStave;
 191     delete[] fNChipsPerStave;
 192     delete[] fNChipsPerLayer;
 193     delete[] fLastChipIndex;
 194     fNStaves = fNHalfStaves = fNModules = fLrChipType = fNChipsPerModule = fLastChipIndex = 0;
 195     fVersion = src.fVersion;
 196     fNLayers = src.fNLayers;
 197     fNChips = src.fNChips;
 198     if (src.fMatSens) {
 199       delete fMatSens;
 200       fMatSens = new TObjArray(fNChips);
 201       fMatSens->SetOwner(kTRUE);
 202       for (int i=0;i<fNChips;i++) {
 203         const TGeoHMatrix* mat = (TGeoHMatrix*) src.fMatSens->At(i);
 204         fMatSens->AddAt(new TGeoHMatrix(*mat),i);
 205       }
 206     }
 207     if (src.fMatT2L) {
 208       delete fMatT2L;
 209       fMatT2L = new TObjArray(fNChips);
 210       fMatT2L->SetOwner(kTRUE);
 211       for (int i=0;i<fNChips;i++) {
 212         const TGeoHMatrix* mat = (TGeoHMatrix*) src.fMatT2L->At(i);
 213         fMatT2L->AddAt(new TGeoHMatrix(*mat),i);
 214       }
 215     }
 216     if (src.fSegm) {
 217       int sz = src.fSegm->GetEntriesFast();
 218       fSegm = new TObjArray(sz);
 219       fSegm->SetOwner(kTRUE);
 220       for (int i=0;i<sz;i++) {
 221         AliITSsegmentation* sg = (AliITSsegmentation*)src.fSegm->UncheckedAt(i);
 222         if (!sg) continue;
 223         fSegm->AddAt(sg->Clone(),i);
 224       }
 225     }
 226     //
 227     if (fNLayers) {
 228       fNStaves   = new Int_t[fNLayers];
 229       fNHalfStaves   = new Int_t[fNLayers];
 230       fNModules     = new Int_t[fNLayers];
 231       fNChipsPerModule = new Int_t[fNLayers];
 232       fNChipRowsPerModule = new Int_t[fNLayers];
 233       fNChipsPerHalfStave = new Int_t[fNLayers];
 234       fNChipsPerStave = new Int_t[fNLayers];
 235       fNChipsPerLayer = new Int_t[fNLayers];
 236       fLrChipType  = new Int_t[fNLayers];
 237       fLastChipIndex   = new Int_t[fNLayers];
 238       for (int i=fNLayers;i--;) {
 239         fNStaves[i] = src.fNStaves[i];
 240         fNHalfStaves[i] = src.fNHalfStaves[i];
 241         fNModules[i]   = src.fNModules[i];
 242         fNChipsPerModule[i] = src.fNChipsPerModule[i];
 243         fNChipRowsPerModule[i] = src.fNChipRowsPerModule[i];
 244         fNChipsPerHalfStave[i] = src.fNChipsPerHalfStave[i];
 245         fNChipsPerStave[i] = src.fNChipsPerStave[i];
 246         fNChipsPerLayer[i] = src.fNChipsPerLayer[i];
 247         fLrChipType[i]  = src.fLrChipType[i];
 248         fLastChipIndex[i] = src.fLastChipIndex[i];
 249       }
 250     }
 251   }
 252   return *this;
 253 }
 254
 255 //______________________________________________________________________
 256 Int_t AliITSUGeomTGeo::GetChipIndex(Int_t lay,Int_t sta,Int_t chipInStave) const
 257 {
 258   // This routine computes the chip index number from the layer,
 259   // stave, and chip number in stave.
 260   // Inputs:
 261   //    Int_t lay  The layer number. Starting from 0.
 262   //    Int_t sta  The stave number. Starting from 0
 263   //    Int_t chipInStave  The chip number in the stave. Starting from 0
 264   //
 265   return GetFirstChipIndex(lay) + fNChipsPerStave[lay]*sta + chipInStave;
 266 }
 267
 268 //______________________________________________________________________
 269 Int_t AliITSUGeomTGeo::GetChipIndex(Int_t lay,Int_t sta, Int_t substa, Int_t chipInSStave) const
 270 {
 271   // This routine computes the chip index number from the layer,
 272   // stave, substave and chip number in substave.
 273   // Inputs:
 274   //    Int_t lay  The layer number. Starting from 0.
 275   //    Int_t sta  The stave number. Starting from 0
 276   //    Int_t substa  The substave number. Starting from 0
 277   //    Int_t chipInSStave  The chip number in the sub stave. Starting from 0
 278   //
 279   int n = GetFirstChipIndex(lay) + fNChipsPerStave[lay]*sta + chipInSStave;
 280   if (fNHalfStaves[lay] && substa>0) n += fNChipsPerHalfStave[lay]*substa;
 281   return n;
 282 }
 283
 284 //______________________________________________________________________
 285 Int_t AliITSUGeomTGeo::GetChipIndex(Int_t lay,Int_t sta, Int_t substa, Int_t md, Int_t chipInMod) const
 286 {
 287   // This routine computes the chip index number from the layer,
 288   // stave, substave module and chip number in module.
 289   // Inputs:
 290   //    Int_t lay  The layer number. Starting from 0.
 291   //    Int_t sta  The stave number. Starting from 0
 292   //    Int_t substa  The substave number. Starting from 0
 293   //    Int_t module  The module number ...
 294   //    Int_t chipInSStave  The chip number in the module. Starting from 0
 295   //
 296   int n = GetFirstChipIndex(lay) + fNChipsPerStave[lay]*sta + chipInMod;
 297   if (fNHalfStaves[lay] && substa>0) n += fNChipsPerHalfStave[lay]*substa;
 298   if (fNModules[lay] && md>0)       n += fNChipsPerModule[lay]*md;
 299   return n;
 300 }
 301
 302 //______________________________________________________________________
 303 Bool_t AliITSUGeomTGeo::GetLayer(Int_t index,Int_t &lay,Int_t &indexInLr)  const
 304 {
 305   // This routine computes the layer number a
 306   // given the chip index. The
 307   // Inputs:
 308   //     Int_t index  The chip index number, starting from zero.
 309   // Outputs:
 310   //     Int_t indexInLr The chip index inside a layer, starting from zero.
 311   //     Int_t lay    The layer number. Starting from 0.
 312   //
 313   lay = GetLayer(index);
 314   indexInLr = index - GetFirstChipIndex(lay);
 315   return kTRUE;
 316   //
 317 }
 318
 319 //______________________________________________________________________
 320 Int_t AliITSUGeomTGeo::GetLayer(Int_t index) const
 321 {
 322   // Get chip layer, from 0
 323   //
 324   int lay = 0;
 325   while(index>fLastChipIndex[lay]) lay++;
 326   return lay;
 327 }
 328
 329 //______________________________________________________________________
 330 Int_t AliITSUGeomTGeo::GetStave(Int_t index) const
 331 {
 332   // Get chip stave, from 0
 333   //
 334   int lay = 0;
 335   while(index>fLastChipIndex[lay]) lay++;
 336   index -= GetFirstChipIndex(lay);
 337   return index/fNChipsPerStave[lay];
 338 }
 339
 340 //______________________________________________________________________
 341 Int_t AliITSUGeomTGeo::GetHalfStave(Int_t index) const
 342 {
 343   // Get chip substave id in stave, from 0
 344   //
 345   int lay = 0;
 346   while(index>fLastChipIndex[lay]) lay++;
 347   if (fNHalfStaves[lay]<0) return -1;
 348   index -= GetFirstChipIndex(lay);
 349   index %= fNChipsPerStave[lay];
 350   return index/fNChipsPerHalfStave[lay];
 351 }
 352
 353 //______________________________________________________________________
 354 Int_t AliITSUGeomTGeo::GetModule(Int_t index) const
 355 {
 356   // Get chip module id in substave, from 0
 357   //
 358   int lay = 0;
 359   while(index>fLastChipIndex[lay]) lay++;
 360   if (fNModules[lay]<0) return 0;
 361   index -= GetFirstChipIndex(lay);
 362   index %= fNChipsPerStave[lay];
 363   if (fNHalfStaves[lay]) index %= fNChipsPerHalfStave[lay];
 364   return index/fNChipsPerModule[lay];
 365 }
 366
 367 //______________________________________________________________________
 368 Int_t AliITSUGeomTGeo::GetChipIdInLayer(Int_t index) const
 369 {
 370   // Get chip number within layer, from 0
 371   //
 372   int lay = 0;
 373   while(index>fLastChipIndex[lay]) lay++;
 374   index -= GetFirstChipIndex(lay);
 375   return index;
 376 }
 377
 378 //______________________________________________________________________
 379 Int_t AliITSUGeomTGeo::GetChipIdInStave(Int_t index) const
 380 {
 381   // Get chip number within stave, from 0
 382   //
 383   int lay = 0;
 384   while(index>fLastChipIndex[lay]) lay++;
 385   index -= GetFirstChipIndex(lay);
 386   return index%fNChipsPerStave[lay];
 387 }
 388
 389 //______________________________________________________________________
 390 Int_t AliITSUGeomTGeo::GetChipIdInHalfStave(Int_t index) const
 391 {
 392   // Get chip number within stave, from 0
 393   //
 394   int lay = 0;
 395   while(index>fLastChipIndex[lay]) lay++;
 396   index -= GetFirstChipIndex(lay);
 397   return index%fNChipsPerHalfStave[lay];
 398 }
 399
 400 //______________________________________________________________________
 401 Int_t AliITSUGeomTGeo::GetChipIdInModule(Int_t index) const
 402 {
 403   // Get chip number within module, from 0
 404   //
 405   int lay = 0;
 406   while(index>fLastChipIndex[lay]) lay++;
 407   index -= GetFirstChipIndex(lay);
 408   return index%fNChipsPerModule[lay];
 409 }
 410
 411 //______________________________________________________________________
 412 Bool_t AliITSUGeomTGeo::GetChipId(Int_t index,Int_t &lay,Int_t &sta,Int_t &hsta, Int_t &mod, Int_t &chip)  const
 413 {
 414   //
 415   // This routine computes the layer, stave, substave, module and chip number
 416   // given the chip index number.
 417   // Inputs:
 418   //     Int_t index  The chip index number, starting from zero.
 419   // Outputs:
 420   //     Int_t lay    The layer number. Starting from 0
 421   //     Int_t sta    The stave number. Starting from 0
 422   //     Int_t ssta   The halfstave number. Starting from 0
 423   //     Int_t mod    The module number. Starting from 0
 424   //     Int_t chip   The detector number. Starting from 0
 425   //
 426   lay  = GetLayer(index);
 427   index -= GetFirstChipIndex(lay);
 428   sta  = index/fNChipsPerStave[lay];
 429   index %= fNChipsPerStave[lay];
 430   hsta = fNHalfStaves[lay]>0 ? index/fNChipsPerHalfStave[lay] : -1;
 431   index %= fNChipsPerHalfStave[lay];
 432   mod  = fNModules[lay]>0 ? index/fNChipsPerModule[lay] : -1;
 433   chip = index%fNChipsPerModule[lay];
 434   //
 435   return kTRUE;
 436 }
 437
 438 //______________________________________________________________________
 439 const char* AliITSUGeomTGeo::GetSymName(Int_t index)  const
 440 {
 441   // Get the TGeoPNEntry symbolic name
 442   // for a given chip identified by 'index'
 443   //
 444   Int_t lay, index2;
 445   if (!GetLayer(index,lay,index2)) return NULL;
 446   // return AliGeomManager::SymName((AliGeomManager::ELayerID)((lay-1)+AliGeomManager::kSPD1),index2);
 447   // RS: this is not optimal, but we cannod access directly AliGeomManager, since the latter has hardwired layers
 448   //  TGeoPNEntry* pne = gGeoManager->GetAlignableEntryByUID( AliGeomManager::LayerToVolUID(lay+1,index2) );
 449   TGeoPNEntry* pne = gGeoManager->GetAlignableEntryByUID( ChipVolUID(index) );
 450   if (!pne) {
 451     AliError(Form("Failed to find alignable entry with index %d: (Lr%d Chip:%d) !",index,lay,index2));
 452     return NULL;
 453   }
 454   return pne->GetName();
 455 }
 456
 457 //______________________________________________________________________
 458 const char* AliITSUGeomTGeo::ComposeSymNameITS()
 459 {
 460   // sym name of the layer
 461   return "ITS";
 462 }
 463
 464 //______________________________________________________________________
 465 const char* AliITSUGeomTGeo::ComposeSymNameLayer(Int_t lr)
 466 {
 467   // sym name of the layer
 468   return Form("%s/%s%d",ComposeSymNameITS(),GetITSLayerPattern(),lr);
 469 }
 470
 471 //______________________________________________________________________
 472 const char* AliITSUGeomTGeo::ComposeSymNameStave(Int_t lr, Int_t stave)
 473 {
 474   // sym name of the stave at given layer
 475   return Form("%s/%s%d",ComposeSymNameLayer(lr),GetITSStavePattern(),stave);
 476 }
 477
 478 //______________________________________________________________________
 479 const char* AliITSUGeomTGeo::ComposeSymNameHalfStave(Int_t lr, Int_t stave, Int_t substave)
 480 {
 481   // sym name of the stave at given layer
 482   return substave>=0 ?
 483     Form("%s/%s%d",ComposeSymNameStave(lr,stave),GetITSHalfStavePattern(),substave) :
 484     ComposeSymNameStave(lr,stave);
 485 }
 486
 487 //______________________________________________________________________
 488 const char* AliITSUGeomTGeo::ComposeSymNameModule(Int_t lr, Int_t stave, Int_t substave, Int_t mod)
 489 {
 490   // sym name of the substave at given layer/stave
 491   return mod>=0 ?
 492     Form("%s/%s%d",ComposeSymNameHalfStave(lr,stave,substave),GetITSModulePattern(),mod) :
 493     ComposeSymNameHalfStave(lr,stave,substave);
 494 }
 495
 496 //______________________________________________________________________
 497 const char* AliITSUGeomTGeo::ComposeSymNameChip(Int_t lr, Int_t sta, Int_t substave, Int_t mod, Int_t chip)
 498 {
 499   // sym name of the chip in the given layer/stave/substave/module
 500   return Form("%s/%s%d",ComposeSymNameModule(lr,sta,substave,mod),GetITSChipPattern(),chip);
 501 }
 502
 503 //______________________________________________________________________
 504 TGeoHMatrix* AliITSUGeomTGeo::GetMatrix(Int_t index)  const
 505 {
 506   // Get the transformation matrix for a given chip 'index'
 507   // by quering the TGeoManager
 508   static TGeoHMatrix matTmp;
 509   TGeoPNEntry *pne = GetPNEntry(index);
 510   if (!pne) return NULL;
 511
 512   TGeoPhysicalNode *pnode = pne->GetPhysicalNode();
 513   if (pnode) return pnode->GetMatrix();
 514
 515   const char* path = pne->GetTitle();
 516   gGeoManager->PushPath(); // Preserve the modeler state.
 517   if (!gGeoManager->cd(path)) {
 518     gGeoManager->PopPath();
 519     AliError(Form("Volume path %s not valid!",path));
 520     return NULL;
 521   }
 522   matTmp = *gGeoManager->GetCurrentMatrix();
 523   gGeoManager->PopPath();
 524   return &matTmp;
 525 }
 526
 527 //______________________________________________________________________
 528 Bool_t AliITSUGeomTGeo::GetTranslation(Int_t index, Double_t t[3])  const
 529 {
 530   // Get the translation vector for a given chip 'index'
 531   // by quering the TGeoManager
 532   TGeoHMatrix *m = GetMatrix(index);
 533   if (!m) return kFALSE;
 534
 535   Double_t *trans = m->GetTranslation();
 536   for (Int_t i = 0; i < 3; i++) t[i] = trans[i];
 537
 538   return kTRUE;
 539 }
 540
 541 //______________________________________________________________________
 542 Bool_t AliITSUGeomTGeo::GetRotation(Int_t index, Double_t r[9])  const
 543 {
 544   // Get the rotation matrix for a given chip 'index'
 545   // by quering the TGeoManager
 546   TGeoHMatrix *m = GetMatrix(index);
 547   if (!m) return kFALSE;
 548
 549   Double_t *rot = m->GetRotationMatrix();
 550   for (Int_t i = 0; i < 9; i++) r[i] = rot[i];
 551
 552   return kTRUE;
 553 }
 554
 555 //______________________________________________________________________
 556 Bool_t AliITSUGeomTGeo::GetOrigMatrix(Int_t index, TGeoHMatrix &m) const
 557 {
 558   // Get the original (ideal geometry) TGeo matrix for
 559   // a given chip identified by 'index'.
 560   // The method is slow, so it should be used
 561   // with great care.
 562   m.Clear();
 563
 564   const char *symname = GetSymName(index);
 565   if (!symname) return kFALSE;
 566
 567   return AliGeomManager::GetOrigGlobalMatrix(symname,m);
 568 }
 569
 570 //______________________________________________________________________
 571 Bool_t AliITSUGeomTGeo::GetOrigTranslation(Int_t index, Double_t t[3])  const
 572 {
 573   // Get the original translation vector (ideal geometry)
 574   // for a given chip 'index' by quering the TGeoManager
 575   TGeoHMatrix m;
 576   if (!GetOrigMatrix(index,m)) return kFALSE;
 577
 578   Double_t *trans = m.GetTranslation();
 579   for (Int_t i = 0; i < 3; i++) t[i] = trans[i];
 580
 581   return kTRUE;
 582 }
 583
 584 //______________________________________________________________________
 585 Bool_t AliITSUGeomTGeo::GetOrigRotation(Int_t index, Double_t r[9])  const
 586 {
 587   // Get the original rotation matrix (ideal geometry)
 588   // for a given chip 'index' by quering the TGeoManager
 589   TGeoHMatrix m;
 590   if (!GetOrigMatrix(index,m)) return kFALSE;
 591
 592   Double_t *rot = m.GetRotationMatrix();
 593   for (Int_t i = 0; i < 9; i++) r[i] = rot[i];
 594
 595   return kTRUE;
 596 }
 597
 598 //______________________________________________________________________
 599 TGeoHMatrix* AliITSUGeomTGeo::ExtractMatrixT2L(Int_t index) const
 600 {
 601   // Get the matrix which transforms from the tracking to local r.s.
 602   // The method queries directly the TGeoPNEntry
 603   TGeoPNEntry *pne = GetPNEntry(index);
 604   if (!pne) return NULL;
 605
 606   TGeoHMatrix *m = (TGeoHMatrix*) pne->GetMatrix();
 607   if (!m) AliError(Form("TGeoPNEntry (%s) contains no matrix !",pne->GetName()));
 608
 609   return m;
 610 }
 611
 612 //______________________________________________________________________
 613 Bool_t AliITSUGeomTGeo::GetTrackingMatrix(Int_t index, TGeoHMatrix &m)
 614 {
 615   // Get the matrix which transforms from the tracking r.s. to
 616   // the global one.
 617   // Returns kFALSE in case of error.
 618   m.Clear();
 619
 620   TGeoHMatrix *m1 = GetMatrix(index);
 621   if (!m1) return kFALSE;
 622
 623   const TGeoHMatrix *m2 = GetMatrixT2L(index);
 624   if (!m2) return kFALSE;
 625
 626   m = *m1;
 627   m.Multiply(m2);
 628
 629   return kTRUE;
 630 }
 631
 632 //______________________________________________________________________
 633 TGeoHMatrix* AliITSUGeomTGeo::ExtractMatrixSens(Int_t index) const
 634 {
 635   // Get the transformation matrix of the SENSOR (not ncessary the same as the chip)
 636   // for a given chip 'index' by quering the TGeoManager
 637   Int_t lay,stav,sstav,mod,chipInMod;
 638   GetChipId(index,lay,stav,sstav,mod,chipInMod);
 639   int wrID = fLr2Wrapper[lay];
 640   TString path = Form("/ALIC_1/%s_2/",AliITSUGeomTGeo::GetITSVolPattern());
 641   if (wrID>=0) path += Form("%s%d_1/",GetITSWrapVolPattern(),wrID);
 642   path += Form("%s%d_1/%s%d_%d/",AliITSUGeomTGeo::GetITSLayerPattern(),lay,AliITSUGeomTGeo::GetITSStavePattern(),lay,stav);
 643   if (fNHalfStaves[lay]>0) path += Form("%s%d_%d/",AliITSUGeomTGeo::GetITSHalfStavePattern(),lay,sstav);
 644   if (fNModules[lay]>0)   path += Form("%s%d_%d/",AliITSUGeomTGeo::GetITSModulePattern(),lay,mod);
 645   path += Form("%s%d_%d/%s%d_1",AliITSUGeomTGeo::GetITSChipPattern(),lay,chipInMod,AliITSUGeomTGeo::GetITSSensorPattern(),lay);
 646   static TGeoHMatrix matTmp;
 647   gGeoManager->PushPath();
 648   if (!gGeoManager->cd(path.Data())) {
 649     gGeoManager->PopPath();
 650     AliError(Form("Error in cd-ing to %s",path.Data()));
 651     return 0;
 652   } // end if !gGeoManager
 653   matTmp = *gGeoManager->GetCurrentMatrix(); // matrix may change after cd
 654   //RSS
 655   //  printf("%d/%d/%d %s\n",lay,stav,detInSta,path.Data());
 656   //  mat->Print();
 657   // Retstore the modeler state.
 658   gGeoManager->PopPath();
 659   return &matTmp;
 660 }
 661
 662
 663 //______________________________________________________________________
 664 TGeoPNEntry* AliITSUGeomTGeo::GetPNEntry(Int_t index) const
 665 {
 666   // Get a pointer to the TGeoPNEntry of a chip
 667   // identified by 'index'
 668   // Returns NULL in case of invalid index,
 669   // missing TGeoManager or invalid symbolic name
 670   //
 671   if (index >= fNChips) {
 672     AliError(Form("Invalid ITS chip index: %d (0 -> %d) !",index,fNChips));
 673     return NULL;
 674   }
 675
 676   if (!gGeoManager || !gGeoManager->IsClosed()) {
 677     AliError("Can't get the matrix! gGeoManager doesn't exist or it is still opened!");
 678     return NULL;
 679   }
 680   TGeoPNEntry* pne = gGeoManager->GetAlignableEntryByUID( ChipVolUID(index) );
 681   //  TGeoPNEntry* pne = gGeoManager->GetAlignableEntry(GetSymName(index));
 682   if (!pne) AliError(Form("The index %d does not correspond to a physical entry!",index));
 683   //
 684   return pne;
 685 }
 686
 687 //______________________________________________________________________
 688 void AliITSUGeomTGeo::BuildITS(Bool_t loadSegm)
 689 {
 690   // exract upg ITS parameters from TGeo
 691   if (fVersion!=kITSVNA) {AliWarning("Already built"); return; // already initialized}
 692     if (!gGeoManager) AliFatal("Geometry is not loaded");
 693   }
 694   fNLayers    = ExtractNumberOfLayers();
 695   if (!fNLayers) return;
 696   //
 697   fNStaves         = new Int_t[fNLayers];
 698   fNHalfStaves     = new Int_t[fNLayers];
 699   fNModules        = new Int_t[fNLayers];
 700   fNChipsPerModule = new Int_t[fNLayers];
 701   fNChipRowsPerModule = new Int_t[fNLayers];
 702   fNChipsPerHalfStave = new Int_t[fNLayers];
 703   fNChipsPerStave  = new Int_t[fNLayers];
 704   fNChipsPerLayer  = new Int_t[fNLayers];
 705   fLrChipType      = new Int_t[fNLayers];
 706   fLastChipIndex   = new Int_t[fNLayers];
 707   fNChips = 0;
 708
 709   for (int i=0;i<fNLayers;i++) {
 710     fLrChipType[i]      = ExtractLayerChipType(i);
 711     fNStaves[i]         = ExtractNumberOfStaves(i);
 712     fNHalfStaves[i]     = ExtractNumberOfHalfStaves(i);
 713     fNModules[i]        = ExtractNumberOfModules(i);
 714     fNChipsPerModule[i] = ExtractNChipsPerModule(i,fNChipRowsPerModule[i]);
 715     fNChipsPerHalfStave[i]= fNChipsPerModule[i]*Max(1,fNModules[i]);
 716     fNChipsPerStave[i]    = fNChipsPerHalfStave[i]*Max(1,fNHalfStaves[i]);
 717     fNChipsPerLayer[i]    = fNChipsPerStave[i]*fNStaves[i];
 718     fNChips              += fNChipsPerLayer[i];
 719     fLastChipIndex[i]     = fNChips-1;
 720   }
 721   //
 722   FetchMatrices();
 723   fVersion = kITSVUpg;
 724   //
 725   if (loadSegm) {  // fetch segmentations
 726     fSegm = new TObjArray();
 727     AliITSUSegmentationPix::LoadSegmentations(fSegm,GetITSsegmentationFileName());
 728   }
 729   //
 730 }
 731
 732 //______________________________________________________________________
 733 Int_t AliITSUGeomTGeo::ExtractNumberOfLayers()
 734 {
 735   // Determines the number of layers in the Upgrade Geometry
 736   //
 737   Int_t numberOfLayers = 0;
 738   //
 739   TGeoVolume *itsV = gGeoManager->GetVolume(GetITSVolPattern());
 740   if (!itsV) AliFatal(Form("ITS volume %s is not in the geometry",GetITSVolPattern()));
 741   SetUIDShift(itsV->GetUniqueID());
 742   //
 743   // Loop on all ITSV nodes, count Layer volumes by checking names
 744   // Build on the fly layer - wrapper correspondence
 745   TObjArray* nodes = itsV->GetNodes();
 746   Int_t nNodes = nodes->GetEntriesFast();
 747   //
 748   int nWrp = 0;
 749   for (Int_t j=0; j<nNodes; j++) {
 750     TGeoNode* nd = (TGeoNode*)nodes->At(j);
 751     const char* name = nd->GetName();
 752     if      (strstr(name,GetITSLayerPattern())) numberOfLayers++;
 753     else if (strstr(name,GetITSWrapVolPattern())) { // this is a wrapper volume, may cointain layers
 754       TObjArray* nodesW = nd->GetNodes();
 755       int nNodesW = nodesW->GetEntriesFast();
 756       for (Int_t jw=0; jw<nNodesW; jw++) {
 757         TGeoNode* ndW = (TGeoNode*)nodesW->At(jw);
 758         if (strstr(ndW->GetName(),GetITSLayerPattern())) fLr2Wrapper[numberOfLayers++] = nWrp;
 759       }
 760       nWrp++;
 761     }
 762   }
 763   //
 764   return numberOfLayers;
 765 }
 766
 767 //______________________________________________________________________
 768 Int_t AliITSUGeomTGeo::ExtractNumberOfStaves(Int_t lay) const
 769 {
 770   // Determines the number of layers in the Upgrade Geometry
 771   //
 772   // Inputs:
 773   //   lay: layer number, starting from 0
 774   //
 775   // MS
 776   Int_t numberOfStaves = 0;
 777   char laynam[30];
 778   snprintf(laynam, 30, "%s%d",GetITSLayerPattern(),lay);
 779   TGeoVolume* volLr = gGeoManager->GetVolume(laynam);
 780   if (!volLr) AliFatal(Form("can't find %s volume",laynam));
 781   //
 782   // Loop on all layer nodes, count Stave volumes by checking names
 783   Int_t nNodes = volLr->GetNodes()->GetEntries();
 784   for (Int_t j=0; j<nNodes; j++) {
 785     //    AliInfo(Form("L%d %d of %d %s %s -> %d",lay,j,nNodes,volLr->GetNodes()->At(j)->GetName(),GetITSStavePattern(),numberOfStaves));
 786     if (strstr(volLr->GetNodes()->At(j)->GetName(),GetITSStavePattern())) numberOfStaves++;
 787   }
 788   //
 789   return numberOfStaves;
 790   //
 791 }
 792
 793 //______________________________________________________________________
 794 Int_t AliITSUGeomTGeo::ExtractNumberOfHalfStaves(Int_t lay) const
 795 {
 796   // Determines the number of substaves in the stave of the layer
 797   //
 798   // Inputs:
 799   //   lay: layer number, starting from 0
 800   //
 801   // MS
 802   if (fgITSHalfStaveName.IsNull()) return 0; // for the setup w/o substave defined the stave and the substave is the same thing
 803   Int_t nSS = 0;
 804   char stavnam[30];
 805   snprintf(stavnam, 30, "%s%d", GetITSStavePattern(),lay);
 806   TGeoVolume* volLd = gGeoManager->GetVolume(stavnam);
 807   if (!volLd) AliFatal(Form("can't find %s volume",stavnam));
 808   //
 809   // Loop on all stave nodes, count Chip volumes by checking names
 810   Int_t nNodes = volLd->GetNodes()->GetEntries();
 811   for (Int_t j=0; j<nNodes; j++) if (strstr(volLd->GetNodes()->At(j)->GetName(),GetITSHalfStavePattern())) nSS++;
 812   //
 813   return nSS;
 814   //
 815 }
 816
 817 //______________________________________________________________________
 818 Int_t AliITSUGeomTGeo::ExtractNumberOfModules(Int_t lay) const
 819 {
 820   // Determines the number of modules in substave in the stave of the layer
 821   //
 822   // Inputs:
 823   //   lay: layer number, starting from 0
 824   //
 825   // for the setup w/o modules defined the module and the stave or the substave is the same thing
 826   if (fgITSModuleName.IsNull()) return 0;
 827   char stavnam[30];
 828   TGeoVolume* volLd = 0;
 829   if (!fgITSHalfStaveName.IsNull()) {
 830     snprintf(stavnam, 30, "%s%d", GetITSHalfStavePattern(),lay);
 831     volLd = gGeoManager->GetVolume(stavnam);
 832   }
 833   if (!volLd) { // no substaves, check staves
 834     snprintf(stavnam, 30, "%s%d", GetITSStavePattern(),lay);
 835     volLd = gGeoManager->GetVolume(stavnam);
 836   }
 837   if (!volLd) return 0;
 838   Int_t nMod = 0;
 839   //
 840   // Loop on all substave nodes, count module volumes by checking names
 841   Int_t nNodes = volLd->GetNodes()->GetEntries();
 842   for (Int_t j=0; j<nNodes; j++) if (strstr(volLd->GetNodes()->At(j)->GetName(),GetITSModulePattern())) nMod++;
 843   //
 844   return nMod;
 845   //
 846 }
 847
 848 //______________________________________________________________________
 849 Int_t AliITSUGeomTGeo::ExtractNChipsPerModule(Int_t lay, int &nrow)  const
 850 {
 851   // Determines the number of chips per module on the (sub)stave in the Upgrade Geometry
 852   // Also extract the layout: span of module centers in Z and X
 853   // Inputs:
 854   //   lay: layer number from 0
 855   // MS
 856   Int_t numberOfChips = 0;
 857   char stavnam[30];
 858   TGeoVolume* volLd = 0;
 859   if (!fgITSModuleName.IsNull()) {
 860     snprintf(stavnam, 30, "%s%d", GetITSModulePattern(),lay);
 861     volLd = gGeoManager->GetVolume(stavnam);
 862   }
 863   if (!volLd) { // no modules on this layer, check substaves
 864     if (!fgITSHalfStaveName.IsNull()) {
 865       snprintf(stavnam, 30, "%s%d", GetITSHalfStavePattern(),lay);
 866       volLd = gGeoManager->GetVolume(stavnam);
 867     }
 868   }
 869   if (!volLd) { // no substaves on this layer, check staves
 870     snprintf(stavnam, 30, "%s%d", GetITSStavePattern(),lay);
 871     volLd = gGeoManager->GetVolume(stavnam);
 872   }
 873   if (!volLd) AliFatal(Form("can't find volume containing chips on layer %d",lay));
 874   //
 875   // Loop on all stave nodes, count Chip volumes by checking names
 876   Int_t nNodes = volLd->GetNodes()->GetEntries();
 877   //
 878   double xmin=1e9,xmax=-1e9, zmin=1e9,zmax=-1e9;
 879   double lab[3],loc[3]={0,0,0};
 880   double dx=-1,dz=-1;
 881   for (Int_t j=0; j<nNodes; j++) {
 882     //    AliInfo(Form("L%d %d of %d %s %s -> %d",lay,j,nNodes,volLd->GetNodes()->At(j)->GetName(),GetITSChipPattern(),numberOfChips));
 883     TGeoNodeMatrix* node = (TGeoNodeMatrix*)volLd->GetNodes()->At(j);
 884     if (!strstr(node->GetName(),GetITSChipPattern())) continue;
 885     node->LocalToMaster(loc,lab);
 886     if (lab[0]>xmax) xmax=lab[0];
 887     if (lab[0]<xmin) xmin=lab[0];
 888     if (lab[2]>zmax) zmax=lab[2];
 889     if (lab[2]<zmin) zmin=lab[2];
 890     //
 891     numberOfChips++;
 892     //
 893     if (dx<0) {
 894       TGeoShape* chShape = node->GetVolume()->GetShape();
 895       TGeoBBox* bbox = dynamic_cast<TGeoBBox*>(chShape);
 896       if (!bbox) {
 897         AliFatal(Form("Chip %s volume is of unprocessed shape %s",node->GetName(),chShape->IsA()->GetName()));
 898       }
 899       else {
 900         dx = 2*bbox->GetDX();
 901         dz = 2*bbox->GetDZ();
 902       }
 903     }
 904   }
 905   //
 906   double spanX = xmax-xmin;
 907   double spanZ = zmax-zmin;
 908   nrow = TMath::Nint(spanX/dx + 1);
 909   int ncol = TMath::Nint(spanZ/dz + 1);
 910   if (nrow*ncol != numberOfChips)
 911     AliError(Form("Inconsistency between Nchips=%d and Nrow*Ncol=%d*%d->%d\n"
 912                   "Extracted chip dimensions (x,z): %.4f %.4f, Module Span: %.4f %.4f",
 913                   numberOfChips,nrow,ncol,nrow*ncol,
 914                   dx,dz,spanX,spanZ));
 915   return numberOfChips;
 916   //
 917 }
 918
 919 //______________________________________________________________________
 920 Int_t AliITSUGeomTGeo::ExtractLayerChipType(Int_t lay)  const
 921 {
 922   // Determines the layer detector type the Upgrade Geometry
 923   //
 924   // Inputs:
 925   //   lay: layer number from 0
 926   // Outputs:
 927   //   none
 928   // Return:
 929   //   detector type id for the layer
 930   // MS
 931   char stavnam[30];
 932   snprintf(stavnam, 30, "%s%d", GetITSLayerPattern(),lay);
 933   TGeoVolume* volLd = gGeoManager->GetVolume(stavnam);
 934   if (!volLd) {AliFatal(Form("can't find %s volume",stavnam)); return -1;}
 935   //
 936   return volLd->GetUniqueID();
 937   //
 938 }
 939
 940 //______________________________________________________________________
 941 UInt_t AliITSUGeomTGeo::ComposeChipTypeID(UInt_t segmId)
 942 {
 943   if (segmId>=kMaxSegmPerChipType) AliFatalClass(Form("Id=%d is >= max.allowed %d",segmId,kMaxSegmPerChipType));
 944   return segmId + kChipTypePix*kMaxSegmPerChipType;
 945 }
 946
 947 //______________________________________________________________________
 948 void AliITSUGeomTGeo::Print(Option_t *) const
 949 {
 950   // print
 951   printf("Geometry version %d, NLayers:%d NChips:%d\n",fVersion,fNLayers,fNChips);
 952   if (fVersion==kITSVNA) return;
 953   for (int i=0;i<fNLayers;i++) {
 954     printf("Lr%2d\tNStav:%2d\tNChips:%2d (%dx%-2d)\tNMod:%d\tNSubSt:%d\tNSt:%3d\tChipType:%3d\tChip#:%5d:%-5d\tWrapVol:%d\n",
 955            i,fNStaves[i],fNChipsPerModule[i],fNChipRowsPerModule[i],
 956            fNChipRowsPerModule[i] ? fNChipsPerModule[i]/fNChipRowsPerModule[i] : 0,
 957            fNModules[i],fNHalfStaves[i],fNStaves[i],
 958            fLrChipType[i],GetFirstChipIndex(i),GetLastChipIndex(i),fLr2Wrapper[i]);
 959   }
 960 }
 961
 962 //______________________________________________________________________
 963 void AliITSUGeomTGeo::FetchMatrices()
 964 {
 965   // store pointer on often used matrices for faster access
 966   if (!gGeoManager) AliFatal("Geometry is not loaded");
 967   fMatSens = new TObjArray(fNChips);
 968   fMatSens->SetOwner(kTRUE);
 969   for (int i=0;i<fNChips;i++) fMatSens->AddAt(new TGeoHMatrix(*ExtractMatrixSens(i)),i);
 970   CreateT2LMatrices();
 971 }
 972
 973 //______________________________________________________________________
 974 void AliITSUGeomTGeo::CreateT2LMatrices()
 975 {
 976   // create tracking to local (Sensor!) matrices
 977   fMatT2L  = new TObjArray(fNChips);
 978   fMatT2L->SetOwner(kTRUE);
 979   TGeoHMatrix matLtoT;
 980   double loc[3]={0,0,0},glo[3];
 981   const double *rotm;
 982   for (int isn=0;isn<fNChips;isn++) {
 983     const TGeoHMatrix* matSens = GetMatrixSens(isn);
 984     if (!matSens) {AliFatal(Form("Failed to get matrix for sensor %d",isn)); return;}
 985     matSens->LocalToMaster(loc,glo);
 986     rotm = matSens->GetRotationMatrix();
 987     Double_t al = -ATan2(rotm[1],rotm[0]);
 988     double sn=Sin(al), cs=Cos(al), r=glo[0]*sn-glo[1]*cs, x=r*sn, y=-r*cs; // sensor plane PCA to origin
 989     TGeoHMatrix* t2l = new TGeoHMatrix();
 990     t2l->RotateZ(ATan2(y,x)*RadToDeg()); // rotate in direction of normal to the sensor plane
 991     t2l->SetDx(x);
 992     t2l->SetDy(y);
 993     t2l->MultiplyLeft(&matSens->Inverse());
 994     fMatT2L->AddAt(t2l,isn);
 995     /*
 996     const double *gtrans = matSens->GetTranslation();
 997     memcpy(&rotMatrix[0], matSens->GetRotationMatrix(), 9*sizeof(Double_t));
 998     Double_t al = -ATan2(rotMatrix[1],rotMatrix[0]);
 999     Double_t rSens = Sqrt(gtrans[0]*gtrans[0] + gtrans[1]*gtrans[1]);
1000     Double_t tanAl = ATan2(gtrans[1],gtrans[0]) - Pi()/2; //angle of tangent
1001     Double_t alTr = tanAl - al;
1002     //
1003     // The X axis of tracking frame must always look outward
1004     loc[1] = rSens/2;
1005     matSens->LocalToMaster(loc,glo);
1006     double rPos = Sqrt(glo[0]*glo[0] + glo[1]*glo[1]);
1007     Bool_t rotOutward = rPos>rSens ? kFALSE : kTRUE;
1008     //
1009     // Transformation matrix
1010     matLtoT.Clear();
1011     matLtoT.SetDx(-rSens*Sin(alTr)); // translation
1012     matLtoT.SetDy(0.);
1013     matLtoT.SetDz(gtrans[2]);
1014     // Rotation matrix
1015     rotMatrix[0]= 0;  rotMatrix[1]= 1;  rotMatrix[2]= 0; // + rotation
1016     rotMatrix[3]=-1;  rotMatrix[4]= 0;  rotMatrix[5]= 0;
1017     rotMatrix[6]= 0;  rotMatrix[7]= 0;  rotMatrix[8]= 1;
1018     //
1019     TGeoRotation rot;
1020     rot.SetMatrix(rotMatrix);
1021     matLtoT.MultiplyLeft(&rot);
1022     if (rotOutward) matLtoT.RotateZ(180.);
1023     // Inverse transformation Matrix
1024     fMatT2L->AddAt(new TGeoHMatrix(matLtoT.Inverse()),isn);
1025     */
1026   }
1027   //
1028 }
1029