Fix for ITS dictionaries
[u/mrichter/AliRoot.git] / ITS / UPGRADE / ITSUpgradeBase / 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 necessary 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   fNLayers    = ExtractNumberOfLayers();
694   if (!fNLayers) return;
695   //
696   fNStaves         = new Int_t[fNLayers];
697   fNHalfStaves     = new Int_t[fNLayers];
698   fNModules        = new Int_t[fNLayers];
699   fNChipsPerModule = new Int_t[fNLayers];
700   fNChipRowsPerModule = new Int_t[fNLayers];
701   fNChipsPerHalfStave = new Int_t[fNLayers];
702   fNChipsPerStave  = new Int_t[fNLayers];
703   fNChipsPerLayer  = new Int_t[fNLayers];
704   fLrChipType      = new Int_t[fNLayers];
705   fLastChipIndex   = new Int_t[fNLayers];
706   fNChips = 0;
707   
708   for (int i=0;i<fNLayers;i++) {
709     fLrChipType[i]      = ExtractLayerChipType(i);
710     fNStaves[i]         = ExtractNumberOfStaves(i);
711     fNHalfStaves[i]     = ExtractNumberOfHalfStaves(i);
712     fNModules[i]        = ExtractNumberOfModules(i);
713     fNChipsPerModule[i] = ExtractNChipsPerModule(i,fNChipRowsPerModule[i]);
714     fNChipsPerHalfStave[i]= fNChipsPerModule[i]*Max(1,fNModules[i]);
715     fNChipsPerStave[i]    = fNChipsPerHalfStave[i]*Max(1,fNHalfStaves[i]);
716     fNChipsPerLayer[i]    = fNChipsPerStave[i]*fNStaves[i];
717     fNChips              += fNChipsPerLayer[i];
718     fLastChipIndex[i]     = fNChips-1;
719   }
720   //
721   FetchMatrices();
722   fVersion = kITSVUpg;
723   //
724   if (loadSegm) {  // fetch segmentations
725     fSegm = new TObjArray();
726     AliITSUSegmentationPix::LoadSegmentations(fSegm,GetITSsegmentationFileName());
727   }
728   //
729 }
730
731 //______________________________________________________________________
732 Int_t AliITSUGeomTGeo::ExtractNumberOfLayers()
733 {
734   // Determines the number of layers in the Upgrade Geometry
735   //
736   Int_t numberOfLayers = 0;
737   //
738   TGeoVolume *itsV = gGeoManager->GetVolume(GetITSVolPattern());
739   if (!itsV) AliFatal(Form("ITS volume %s is not in the geometry",GetITSVolPattern()));
740   SetUIDShift(itsV->GetUniqueID());
741   //
742   // Loop on all ITSV nodes, count Layer volumes by checking names
743   // Build on the fly layer - wrapper correspondence
744   TObjArray* nodes = itsV->GetNodes();
745   Int_t nNodes = nodes->GetEntriesFast();
746   //
747   for (Int_t j=0; j<nNodes; j++) {
748     int lrID = -1;
749     TGeoNode* nd = (TGeoNode*)nodes->At(j);
750     const char* name = nd->GetName();
751     if (strstr(name,GetITSLayerPattern())) {
752       numberOfLayers++;
753       if ( (lrID=ExtractVolumeCopy(name,AliITSUGeomTGeo::GetITSLayerPattern()))<0 ) {
754         AliFatal(Form("Failed to extract layer ID from the %s",name));
755         exit(1);
756       }
757       //
758       fLr2Wrapper[lrID] = -1; // not wrapped
759     }
760     else if (strstr(name,GetITSWrapVolPattern())) { // this is a wrapper volume, may cointain layers
761       int wrID = -1;
762       if ( (wrID=ExtractVolumeCopy(name,AliITSUGeomTGeo::GetITSWrapVolPattern()))<0 ) {
763         AliFatal(Form("Failed to extract wrapper ID from the %s",name));
764         exit(1);
765       }
766       //
767       TObjArray* nodesW = nd->GetNodes();
768       int nNodesW = nodesW->GetEntriesFast();
769       for (Int_t jw=0; jw<nNodesW; jw++) {
770         TGeoNode* ndW = (TGeoNode*)nodesW->At(jw);
771         if (strstr(ndW->GetName(),GetITSLayerPattern())) {
772           if ( (lrID=ExtractVolumeCopy(ndW->GetName(),AliITSUGeomTGeo::GetITSLayerPattern()))<0 ) {
773             AliFatal(Form("Failed to extract layer ID from the %s",name));
774             exit(1);
775           }
776           numberOfLayers++;
777           fLr2Wrapper[lrID] = wrID;
778         }
779       }
780     }
781   }
782   //  
783   return numberOfLayers;
784 }
785
786 //______________________________________________________________________
787 Int_t AliITSUGeomTGeo::ExtractNumberOfStaves(Int_t lay) const
788 {
789   // Determines the number of layers in the Upgrade Geometry
790   //
791   // Inputs:
792   //   lay: layer number, starting from 0
793   //
794   // MS
795   Int_t numberOfStaves = 0;
796   char laynam[30];
797   snprintf(laynam, 30, "%s%d",GetITSLayerPattern(),lay);
798   TGeoVolume* volLr = gGeoManager->GetVolume(laynam);
799   if (!volLr) { AliFatal(Form("can't find %s volume",laynam)); return -1; }
800   //
801   // Loop on all layer nodes, count Stave volumes by checking names
802   Int_t nNodes = volLr->GetNodes()->GetEntries();
803   for (Int_t j=0; j<nNodes; j++) {
804     //    AliInfo(Form("L%d %d of %d %s %s -> %d",lay,j,nNodes,volLr->GetNodes()->At(j)->GetName(),GetITSStavePattern(),numberOfStaves));
805     if (strstr(volLr->GetNodes()->At(j)->GetName(),GetITSStavePattern())) numberOfStaves++;
806   }
807   //
808   return numberOfStaves;
809   //
810 }
811
812 //______________________________________________________________________
813 Int_t AliITSUGeomTGeo::ExtractNumberOfHalfStaves(Int_t lay) const
814 {
815   // Determines the number of substaves in the stave of the layer
816   //
817   // Inputs:
818   //   lay: layer number, starting from 0
819   //
820   // MS
821   if (fgITSHalfStaveName.IsNull()) return 0; // for the setup w/o substave defined the stave and the substave is the same thing
822   Int_t nSS = 0;
823   char stavnam[30];
824   snprintf(stavnam, 30, "%s%d", GetITSStavePattern(),lay);
825   TGeoVolume* volLd = gGeoManager->GetVolume(stavnam);
826   if (!volLd) AliFatal(Form("can't find %s volume",stavnam));
827   //
828   // Loop on all stave nodes, count Chip volumes by checking names
829   Int_t nNodes = volLd->GetNodes()->GetEntries();
830   for (Int_t j=0; j<nNodes; j++) if (strstr(volLd->GetNodes()->At(j)->GetName(),GetITSHalfStavePattern())) nSS++;
831   //
832   return nSS;
833   //
834 }
835
836 //______________________________________________________________________
837 Int_t AliITSUGeomTGeo::ExtractNumberOfModules(Int_t lay) const
838 {
839   // Determines the number of modules in substave in the stave of the layer
840   //
841   // Inputs:
842   //   lay: layer number, starting from 0
843   //
844   // for the setup w/o modules defined the module and the stave or the substave is the same thing
845   if (fgITSModuleName.IsNull()) return 0;
846   char stavnam[30];
847   TGeoVolume* volLd = 0;
848   if (!fgITSHalfStaveName.IsNull()) {
849     snprintf(stavnam, 30, "%s%d", GetITSHalfStavePattern(),lay); 
850     volLd = gGeoManager->GetVolume(stavnam);
851   }
852   if (!volLd) { // no substaves, check staves
853     snprintf(stavnam, 30, "%s%d", GetITSStavePattern(),lay); 
854     volLd = gGeoManager->GetVolume(stavnam);
855   }
856   if (!volLd) return 0;
857   Int_t nMod = 0;
858   //
859   // Loop on all substave nodes, count module volumes by checking names
860   Int_t nNodes = volLd->GetNodes()->GetEntries();
861   for (Int_t j=0; j<nNodes; j++) if (strstr(volLd->GetNodes()->At(j)->GetName(),GetITSModulePattern())) nMod++;
862   //
863   return nMod;
864   //
865 }
866
867 //______________________________________________________________________
868 Int_t AliITSUGeomTGeo::ExtractNChipsPerModule(Int_t lay, int &nrow)  const
869 {
870   // Determines the number of chips per module on the (sub)stave in the Upgrade Geometry
871   // Also extract the layout: span of module centers in Z and X
872   // Inputs:
873   //   lay: layer number from 0
874   // MS
875   Int_t numberOfChips = 0;
876   char stavnam[30];
877   TGeoVolume* volLd = 0;
878   if (!fgITSModuleName.IsNull()) {
879     snprintf(stavnam, 30, "%s%d", GetITSModulePattern(),lay); 
880     volLd = gGeoManager->GetVolume(stavnam);
881   }
882   if (!volLd) { // no modules on this layer, check substaves
883     if (!fgITSHalfStaveName.IsNull()) {
884       snprintf(stavnam, 30, "%s%d", GetITSHalfStavePattern(),lay); 
885       volLd = gGeoManager->GetVolume(stavnam);
886     }
887   }
888   if (!volLd) { // no substaves on this layer, check staves
889     snprintf(stavnam, 30, "%s%d", GetITSStavePattern(),lay);
890     volLd = gGeoManager->GetVolume(stavnam);
891   }
892   if (!volLd) AliFatal(Form("can't find volume containing chips on layer %d",lay));
893   //
894   // Loop on all stave nodes, count Chip volumes by checking names
895   Int_t nNodes = volLd->GetNodes()->GetEntries();
896   //
897   double xmin=1e9,xmax=-1e9, zmin=1e9,zmax=-1e9;
898   double lab[3],loc[3]={0,0,0};
899   double dx=-1,dz=-1;
900   for (Int_t j=0; j<nNodes; j++) {
901     //    AliInfo(Form("L%d %d of %d %s %s -> %d",lay,j,nNodes,volLd->GetNodes()->At(j)->GetName(),GetITSChipPattern(),numberOfChips));
902     TGeoNodeMatrix* node = (TGeoNodeMatrix*)volLd->GetNodes()->At(j);
903     if (!strstr(node->GetName(),GetITSChipPattern())) continue;
904     node->LocalToMaster(loc,lab);
905     if (lab[0]>xmax) xmax=lab[0];
906     if (lab[0]<xmin) xmin=lab[0];    
907     if (lab[2]>zmax) zmax=lab[2];
908     if (lab[2]<zmin) zmin=lab[2];    
909     //
910     numberOfChips++;
911     //
912     if (dx<0) {
913       TGeoShape* chShape = node->GetVolume()->GetShape();
914       TGeoBBox* bbox = dynamic_cast<TGeoBBox*>(chShape);
915       if (!bbox) {
916         AliFatal(Form("Chip %s volume is of unprocessed shape %s",node->GetName(),chShape->IsA()->GetName()));
917       }
918       else {
919         dx = 2*bbox->GetDX();
920         dz = 2*bbox->GetDZ();
921       }
922     }
923   }
924   //
925   double spanX = xmax-xmin;
926   double spanZ = zmax-zmin;  
927   nrow = TMath::Nint(spanX/dx + 1);
928   int ncol = TMath::Nint(spanZ/dz + 1);
929   if (nrow*ncol != numberOfChips) 
930     AliError(Form("Inconsistency between Nchips=%d and Nrow*Ncol=%d*%d->%d\n"
931                   "Extracted chip dimensions (x,z): %.4f %.4f, Module Span: %.4f %.4f",
932                   numberOfChips,nrow,ncol,nrow*ncol,
933                   dx,dz,spanX,spanZ));
934   return numberOfChips;
935   //
936 }
937
938 //______________________________________________________________________
939 Int_t AliITSUGeomTGeo::ExtractLayerChipType(Int_t lay)  const
940 {
941   // Determines the layer detector type the Upgrade Geometry
942   //
943   // Inputs:
944   //   lay: layer number from 0
945   // Outputs:
946   //   none
947   // Return:
948   //   detector type id for the layer
949   // MS
950   char stavnam[30];
951   snprintf(stavnam, 30, "%s%d", GetITSLayerPattern(),lay);
952   TGeoVolume* volLd = gGeoManager->GetVolume(stavnam);
953   if (!volLd) {AliFatal(Form("can't find %s volume",stavnam)); return -1;}
954   //
955   return volLd->GetUniqueID();
956   //
957 }
958
959 //______________________________________________________________________
960 UInt_t AliITSUGeomTGeo::ComposeChipTypeID(UInt_t segmId)
961 {
962   if (segmId>=kMaxSegmPerChipType) AliFatalClass(Form("Id=%d is >= max.allowed %d",segmId,kMaxSegmPerChipType));
963   return segmId + kChipTypePix*kMaxSegmPerChipType;
964 }
965
966 //______________________________________________________________________
967 void AliITSUGeomTGeo::Print(Option_t *) const
968 {
969   // print
970   printf("Geometry version %d, NLayers:%d NChips:%d\n",fVersion,fNLayers,fNChips);
971   if (fVersion==kITSVNA) return;
972   for (int i=0;i<fNLayers;i++) {
973     printf("Lr%2d\tNStav:%2d\tNChips:%2d (%dx%-2d)\tNMod:%d\tNSubSt:%d\tNSt:%3d\tChipType:%3d\tChip#:%5d:%-5d\tWrapVol:%d\n",
974            i,fNStaves[i],fNChipsPerModule[i],fNChipRowsPerModule[i],
975            fNChipRowsPerModule[i] ? fNChipsPerModule[i]/fNChipRowsPerModule[i] : 0,
976            fNModules[i],fNHalfStaves[i],fNStaves[i],
977            fLrChipType[i],GetFirstChipIndex(i),GetLastChipIndex(i),fLr2Wrapper[i]);
978   }
979 }
980
981 //______________________________________________________________________
982 void AliITSUGeomTGeo::FetchMatrices()
983 {
984   // store pointer on often used matrices for faster access
985   if (!gGeoManager) AliFatal("Geometry is not loaded");
986   fMatSens = new TObjArray(fNChips);
987   fMatSens->SetOwner(kTRUE);
988   for (int i=0;i<fNChips;i++) fMatSens->AddAt(new TGeoHMatrix(*ExtractMatrixSens(i)),i);
989   CreateT2LMatrices();
990 }
991
992 //______________________________________________________________________
993 void AliITSUGeomTGeo::CreateT2LMatrices()
994 {
995   // create tracking to local (Sensor!) matrices
996   fMatT2L  = new TObjArray(fNChips);  
997   fMatT2L->SetOwner(kTRUE);
998   TGeoHMatrix matLtoT;
999   double locA[3]={-100,0,0},locB[3]={100,0,0},gloA[3],gloB[3];
1000   for (int isn=0;isn<fNChips;isn++) {
1001     const TGeoHMatrix* matSens = GetMatrixSens(isn);
1002     if (!matSens) {AliFatal(Form("Failed to get matrix for sensor %d",isn)); return;}
1003     matSens->LocalToMaster(locA,gloA);
1004     matSens->LocalToMaster(locB,gloB);
1005     double dx = gloB[0]-gloA[0];
1006     double dy = gloB[1]-gloA[1];
1007     double t = (gloB[0]*dx+gloB[1]*dy)/(dx*dx+dy*dy),x=gloB[0]-dx*t,y=gloB[1]-dy*t;
1008     TGeoHMatrix* t2l = new TGeoHMatrix();
1009     t2l->RotateZ(ATan2(y,x)*RadToDeg()); // rotate in direction of normal to the sensor plane
1010     t2l->SetDx(x);
1011     t2l->SetDy(y);
1012     t2l->MultiplyLeft(&matSens->Inverse());
1013     fMatT2L->AddAt(t2l,isn);
1014     /*
1015     const double *gtrans = matSens->GetTranslation();
1016     memcpy(&rotMatrix[0], matSens->GetRotationMatrix(), 9*sizeof(Double_t));
1017     Double_t al = -ATan2(rotMatrix[1],rotMatrix[0]);
1018     Double_t rSens = Sqrt(gtrans[0]*gtrans[0] + gtrans[1]*gtrans[1]);
1019     Double_t tanAl = ATan2(gtrans[1],gtrans[0]) - Pi()/2; //angle of tangent
1020     Double_t alTr = tanAl - al;
1021     //
1022     // The X axis of tracking frame must always look outward
1023     loc[1] = rSens/2;
1024     matSens->LocalToMaster(loc,glo);
1025     double rPos = Sqrt(glo[0]*glo[0] + glo[1]*glo[1]);
1026     Bool_t rotOutward = rPos>rSens ? kFALSE : kTRUE;
1027     //
1028     // Transformation matrix
1029     matLtoT.Clear();
1030     matLtoT.SetDx(-rSens*Sin(alTr)); // translation
1031     matLtoT.SetDy(0.);
1032     matLtoT.SetDz(gtrans[2]);
1033     // Rotation matrix
1034     rotMatrix[0]= 0;  rotMatrix[1]= 1;  rotMatrix[2]= 0; // + rotation
1035     rotMatrix[3]=-1;  rotMatrix[4]= 0;  rotMatrix[5]= 0;
1036     rotMatrix[6]= 0;  rotMatrix[7]= 0;  rotMatrix[8]= 1;
1037     //
1038     TGeoRotation rot;
1039     rot.SetMatrix(rotMatrix);
1040     matLtoT.MultiplyLeft(&rot);
1041     if (rotOutward) matLtoT.RotateZ(180.);
1042     // Inverse transformation Matrix
1043     fMatT2L->AddAt(new TGeoHMatrix(matLtoT.Inverse()),isn);
1044     */
1045   }
1046   //
1047 }
1048
1049 //______________________________________________________________________
1050 Int_t AliITSUGeomTGeo::ExtractVolumeCopy(const char* name, const char* prefix) const
1051 {
1052   // extract Number following the prefix in the name string
1053   TString nms = name;
1054   if (!nms.BeginsWith(prefix)) return -1;
1055   nms.Remove(0,strlen(prefix));
1056   if (!isdigit(nms.Data()[0])) return -1;
1057   return nms.Atoi();
1058   //
1059 }