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