ITS/AliITSgeom.cxx

   1 ///////////////////////////////////////////////////////////////////////
   2 // ITS geometry manimulaiton routines.                               //
   3 // Created April 15 1999.                                            //
   4 // version: 0.0.0                                                    //
   5 // By: Bjorn S. Nilsen                                               //
   6 // version: 0.0.1                                                    //
   7 // Updated May 27 1999.                                              //
   8 // Added Cylinderical random and global based changes.               //
   9 // Added  function PrintComparison.                                  //
  10 ///////////////////////////////////////////////////////////////////////
  11 #include <iostream.h>
  12 #include <fstream.h>
  13 #include <iomanip.h>
  14 #include <stdio.h>
  15 #include "AliITSgeom.h"
  16 #include "TRandom.h"
  17
  18 ClassImp(AliITSgeom)
  19
  20 //_____________________________________________________________________
  21 AliITSgeom::AliITSgeom(){
  22 ////////////////////////////////////////////////////////////////////////
  23 //     The default constructor for the AliITSgeom class. It, by default,
  24 // sets fNlayers to zero and zeros all pointers.
  25 ////////////////////////////////////////////////////////////////////////
  26   // Default constructor.
  27   // Do not allocate anything zero everything
  28    fNlayers = 0;
  29    fNlad    = 0;
  30    fNdet    = 0;
  31    fg       = 0;
  32    fShape   = 0;
  33    return;
  34 }
  35
  36 //_____________________________________________________________________
  37 AliITSgeom::~AliITSgeom(){
  38 ////////////////////////////////////////////////////////////////////////
  39 //     The destructor for the AliITSgeom class. If the arrays fNlad,
  40 // fNdet, or fg have had memory allocated to them, there pointer values
  41 // are non zero, then this memory space is freed and they are set
  42 // to zero. In addition, fNlayers is set to zero. The destruction of
  43 // TObjArray fShape is, by default, handled by the TObjArray destructor.
  44 ////////////////////////////////////////////////////////////////////////
  45   // Default destructor.
  46   // if arrays exist delet them. Then set everything to zero.
  47    if(fg!=0){
  48       for(Int_t i=0;i<fNlayers;i++) delete[] fg[i];
  49       delete[] fg;
  50    } // end if fg!=0
  51    if(fNlad!=0) delete[] fNlad;
  52    if(fNdet!=0) delete[] fNdet;
  53    fNlayers = 0;
  54    fNlad    = 0;
  55    fNdet    = 0;
  56    fg       = 0;
  57    return;
  58 }
  59
  60 //_____________________________________________________________________
  61 AliITSgeom::AliITSgeom(const char *filename){
  62 ////////////////////////////////////////////////////////////////////////
  63 //     The constructor for the AliITSgeom class. All of the data to fill
  64 // this structure is read in from the file given my the input filename.
  65 ////////////////////////////////////////////////////////////////////////
  66    FILE     *pf;
  67    Int_t    i;
  68    ITS_geom *g;
  69    Int_t    l,a,d;
  70    Float_t  x,y,z,o,p,q,r,s,t;
  71    Double_t or,pr,qr,rr,sr,tr; // Radians
  72    Double_t lr[9];
  73    Double_t si; // sin(angle)
  74    Double_t PI = TMath::Pi(), byPI = PI/180.;
  75
  76    pf = fopen(filename,"r");
  77
  78    fNlayers = 6; // set default number of ladders
  79    fNlad    = new Int_t[fNlayers];
  80    fNdet    = new Int_t[fNlayers];
  81    // find the number of laders and detectors in this geometry.
  82    for(i=0;i<fNlayers;i++){fNlad[i]=fNdet[i]=0;} // zero out arrays
  83    for(;;){ // for ever loop
  84       i = fscanf(pf,"%d %d %d %f %f %f %f %f %f %f %f %f",
  85                      &l,&a,&d,&x,&y,&z,&o,&p,&q,&r,&s,&t);
  86       if(i==EOF) break;
  87       if(l<1 || l>fNlayers) {
  88          printf("error in file %s layer=%d min is 1 max is %d/n",
  89                  filename,l,fNlayers);
  90          continue;
  91       }// end if l
  92       if(fNlad[l-1]<a) fNlad[l-1] = a;
  93       if(fNdet[l-1]<d) fNdet[l-1] = d;
  94    } // end for ever loop
  95    // counted the number of laders and detectors now allocate space.
  96    fg = new ITS_geom* [fNlayers];
  97    for(i=0;i<fNlayers;i++){
  98       fg[i] = 0;
  99       l = fNlad[i]*fNdet[i];
 100       fg[i] = new ITS_geom[l]; // allocate space for transforms
 101    } // end for i
 102
 103    // Set up Shapes for a default configuration of 6 layers.
 104    fShape = new TObjArray;
 105    AddShape((TObject *) new AliITSgeomSPD());  // shape 0
 106    AddShape((TObject *) new AliITSgeomSDD());  // shape 1
 107    AddShape((TObject *) new AliITSgeomSPD());  // shape 2
 108
 109    // prepair to read in transforms
 110    rewind(pf); // start over reading file
 111    for(;;){ // for ever loop
 112       i = fscanf(pf,"%d %d %d %f %f %f %f %f %f %f %f %f",
 113                      &l,&a,&d,&x,&y,&z,&o,&p,&q,&r,&s,&t);
 114       if(i==EOF) break;
 115       if(l<1 || l>fNlayers) {
 116          printf("error in file %s layer=%d min is 1 max is %d/n",
 117                  filename,l,fNlayers);
 118          continue;
 119       }// end if l
 120       l--; a--; d--; // shift layer, lader, and detector counters to zero base
 121       i = d + a*fNdet[l]; // position of this detector
 122       g = &(fg[l][i]);
 123
 124       or = byPI*o;
 125       pr = byPI*p;
 126       qr = byPI*q;
 127       rr = byPI*r;
 128       sr = byPI*s;
 129       tr = byPI*t;
 130
 131       g->fx0   = x;
 132       g->fy0   = y;
 133       g->fz0   = z;
 134       si    = sin(or);if(o== 90.0) si = +1.0;
 135                       if(o==270.0) si = -1.0;
 136                       if(o==  0.0||o==180.) si = 0.0;
 137       lr[0] = si * cos(pr);
 138       lr[1] = si * sin(pr);
 139       lr[2] = cos(or);if(o== 90.0||o==270.) lr[2] = 0.0;
 140                       if(o== 0.0)           lr[2] = +1.0;
 141                       if(o==180.0)          lr[2] = -1.0;
 142       si    =  sin(qr);if(q== 90.0) si = +1.0;
 143                        if(q==270.0) si = -1.0;
 144                        if(q==  0.0||q==180.) si = 0.0;
 145       lr[3] = si * cos(rr);
 146       lr[4] = si * sin(rr);
 147       lr[5] = cos(qr);if(q== 90.0||q==270.) lr[5] = 0.0;
 148                       if(q==  0.0)          lr[5] = +1.0;
 149                       if(q==180.0)          lr[5] = -1.0;
 150       si    = sin(sr);if(r== 90.0) si = +1.0;
 151                       if(r==270.0) si = -1.0;
 152                       if(r==  0.0||r==180.) si = 0.0;
 153       lr[6] = si * cos(tr);
 154       lr[7] = si * sin(tr);
 155       lr[8] = cos(sr);if(r== 90.0||r==270.0) lr[8] =  0.0;
 156                       if(r==  0.0)           lr[8] = +1.0;
 157                       if(r==180.0)           lr[8] = -1.0;
 158       // Normalize these elements
 159       for(a=0;a<3;a++){// reuse float si and integers a and d.
 160          si = 0.0;
 161          for(d=0;d<3;d++) si += lr[3*a+d]*lr[3*a+d];
 162          si = TMath::Sqrt(1./si);
 163          for(d=0;d<3;d++) g->fr[3*a+d] = lr[3*a+d] = si*lr[3*a+d];
 164       } // end for a
 165       // get angles from matrix up to a phase of 180 degrees.
 166       or     = atan2(lr[7],lr[8]);if(or<0.0) or += 2.0*PI;
 167       pr     = asin(lr[2]);       if(pr<0.0) pr += 2.0*PI;
 168       qr     = atan2(lr[3],lr[0]);if(qr<0.0) qr += 2.0*PI;
 169       g->frx = or;
 170       g->fry = pr;
 171       g->frz = qr;
 172       // l = layer-1 at this point.
 173            if(l==0||l==1) g->fShapeIndex = 0; // SPD's
 174       else if(l==2||l==3) g->fShapeIndex = 1; // SDD's
 175       else if(l==4||l==5) g->fShapeIndex = 2; // SSD's
 176    } // end for ever loop
 177    fclose(pf);
 178 }
 179
 180 //________________________________________________________________________
 181 AliITSgeom::AliITSgeom(AliITSgeom &source){
 182 ////////////////////////////////////////////////////////////////////////
 183 //     The copy constructor for the AliITSgeom class. It calls the
 184 // = operator function. See the = operator function for more details.
 185 ////////////////////////////////////////////////////////////////////////
 186    source = *this;  // Just use the = operator for now.
 187    return;
 188 }
 189
 190 //________________________________________________________________________
 191 void AliITSgeom::operator=(AliITSgeom &source){
 192 ////////////////////////////////////////////////////////////////////////
 193 //     The = operator function for the AliITSgeom class. It makes an
 194 // independent copy of the class in such a way that any changes made
 195 // to the copied class will not affect the source class in any way.
 196 // This is required for many ITS alignment studies where the copied
 197 // class is then modified by introducing some misalignment.
 198 ////////////////////////////////////////////////////////////////////////
 199    Int_t i,j,k;
 200
 201    if(this == &source) return; // don't assign to ones self.
 202
 203    // if there is an old structure allocated delete it first.
 204    if(fg != 0){
 205       for(i=0;i<fNlayers;i++) delete[] fg[i];
 206       delete[] fg;
 207    } // end if fg != 0
 208    if(fNlad != 0) delete[] fNlad;
 209    if(fNdet != 0) delete[] fNdet;
 210
 211    fNlayers = source.fNlayers;
 212    fNlad = new Int_t[fNlayers];
 213    for(i=0;i<fNlayers;i++) fNlad[i] = source.fNlad[i];
 214    fNdet = new Int_t[fNlayers];
 215    for(i=0;i<fNlayers;i++) fNdet[i] = source.fNdet[i];
 216    fShape = new TObjArray(*(source.fShape));//This does not make a proper copy.
 217    fg = new ITS_geom* [fNlayers];
 218    for(i=0;i<fNlayers;i++){
 219       fg[i] = new ITS_geom[fNlad[i]*fNdet[i]];
 220       for(j=0;j<(fNlad[i]*fNdet[i]);j++){
 221           fg[i][j].fShapeIndex = source.fg[i][j].fShapeIndex;
 222           fg[i][j].fx0 = source.fg[i][j].fx0;
 223           fg[i][j].fy0 = source.fg[i][j].fy0;
 224           fg[i][j].fz0 = source.fg[i][j].fz0;
 225           fg[i][j].frx = source.fg[i][j].frx;
 226           fg[i][j].fry = source.fg[i][j].fry;
 227           fg[i][j].frz = source.fg[i][j].frz;
 228           for(k=0;k<9;k++) fg[i][j].fr[k] = source.fg[i][j].fr[k];
 229       } // end for j
 230    } // end for i
 231    return;
 232 }
 233
 234
 235 //________________________________________________________________________
 236 void AliITSgeom::GtoL(Int_t lay,Int_t lad,Int_t det,
 237                        const Float_t *g,Float_t *l){
 238 ////////////////////////////////////////////////////////////////////////
 239 //     The function that does the global ALICE Cartesian coordinate
 240 // to local active volume detector Cartesian coordinate transformation.
 241 // The local detector coordinate system is determined by the layer,
 242 // ladder, and detector numbers. The global coordinates are entered by
 243 // the three element Float_t array g and the local coordinate values
 244 // are returned by the three element Float_t array l. The order of the
 245 // three elements are g[0]=x, g[1]=y, and g[2]=z, similarly for l.
 246 ////////////////////////////////////////////////////////////////////////
 247    Double_t x,y,z;
 248    ITS_geom *gl;
 249
 250    lay--; lad--; det--;
 251    gl = &(fg[lay][fNdet[lay]*lad+det]);
 252
 253    x    = g[0] - gl->fx0;
 254    y    = g[1] - gl->fy0;
 255    z    = g[2] - gl->fz0;
 256    l[0] = gl->fr[0]*x + gl->fr[1]*y + gl->fr[2]*z;
 257    l[1] = gl->fr[3]*x + gl->fr[4]*y + gl->fr[5]*z;
 258    l[2] = gl->fr[6]*x + gl->fr[7]*y + gl->fr[8]*z;
 259    return;
 260 }
 261
 262 //________________________________________________________________________
 263 void AliITSgeom::GtoL(const Int_t *id,const Float_t *g,Float_t *l){
 264 ////////////////////////////////////////////////////////////////////////
 265 //     The function that does the local active volume detector Cartesian
 266 // coordinate to global ALICE Cartesian coordinate transformation.
 267 // The local detector coordinate system is determined by the layer,
 268 // ladder, and detector numbers. The local coordinates are entered by
 269 // the three element Float_t array l and the global coordinate values
 270 // are returned by the three element Float_t array g. The order of the
 271 // three elements are l[0]=x, l[1]=y, and l[2]=z, similarly for g.
 272 ////////////////////////////////////////////////////////////////////////
 273    Int_t    lay,lad,det;
 274    Double_t x,y,z;
 275    ITS_geom *gl;
 276
 277    lay = id[0]; lad = id[1]; det = id[2];
 278    lay--; lad--; det--;
 279    gl = &(fg[lay][fNdet[lay]*lad+det]);
 280
 281    x    = g[0] - gl->fx0;
 282    y    = g[1] - gl->fy0;
 283    z    = g[2] - gl->fz0;
 284    l[0] = gl->fr[0]*x + gl->fr[1]*y + gl->fr[2]*z;
 285    l[1] = gl->fr[3]*x + gl->fr[4]*y + gl->fr[5]*z;
 286    l[2] = gl->fr[6]*x + gl->fr[7]*y + gl->fr[8]*z;
 287    return;
 288 }
 289 //________________________________________________________________________
 290 void AliITSgeom::GtoL(Int_t index,const Float_t *g,Float_t *l){
 291 ////////////////////////////////////////////////////////////////////////
 292 //     The function that does the local active volume detector Cartesian
 293 // coordinate to global ALICE Cartesian coordinate transformation.
 294 // The local detector coordinate system is determined by the detector
 295 // index numbers (see GetModuleIndex and GetModuleID). The local
 296 // coordinates are entered by the three element Float_t array l and the
 297 // global coordinate values are returned by the three element Float_t array g.
 298 // The order of the three elements are l[0]=x, l[1]=y, and l[2]=z, similarly
 299 // for g.
 300 ////////////////////////////////////////////////////////////////////////
 301    Int_t    lay,lad,det;
 302    Double_t x,y,z;
 303    ITS_geom *gl;
 304
 305    this->GetModuleId(index,lay,lad,det);
 306    lay--; lad--; det--;
 307    gl = &(fg[lay][fNdet[lay]*lad+det]);
 308
 309    x    = g[0] - gl->fx0;
 310    y    = g[1] - gl->fy0;
 311    z    = g[2] - gl->fz0;
 312    l[0] = gl->fr[0]*x + gl->fr[1]*y + gl->fr[2]*z;
 313    l[1] = gl->fr[3]*x + gl->fr[4]*y + gl->fr[5]*z;
 314    l[2] = gl->fr[6]*x + gl->fr[7]*y + gl->fr[8]*z;
 315    return;
 316 }
 317
 318 //________________________________________________________________________
 319 void AliITSgeom::LtoG(Int_t lay,Int_t lad,Int_t det,
 320                       const Float_t *l,Float_t *g){
 321 ////////////////////////////////////////////////////////////////////////
 322 //     The function that does the local active volume detector Cartesian
 323 // coordinate to global ALICE Cartesian coordinate transformation.
 324 // The local detector coordinate system is determined by the layer,
 325 // ladder, and detector numbers. The local coordinates are entered by
 326 // the three element Float_t array l and the global coordinate values
 327 // are returned by the three element Float_t array g. The order of the
 328 // three elements are l[0]=x, l[1]=y, and l[2]=z, similarly for g.
 329 ////////////////////////////////////////////////////////////////////////
 330    Double_t x,y,z;
 331    ITS_geom *gl;
 332
 333    lay--; lad--; det--;
 334    gl   = &(fg[lay][fNdet[lay]*lad+det]);
 335
 336    x    = gl->fr[0]*l[0] + gl->fr[3]*l[1] + gl->fr[6]*l[2];
 337    y    = gl->fr[1]*l[0] + gl->fr[4]*l[1] + gl->fr[7]*l[2];
 338    z    = gl->fr[2]*l[0] + gl->fr[5]*l[1] + gl->fr[8]*l[2];
 339    g[0] = x + gl->fx0;
 340    g[1] = y + gl->fy0;
 341    g[2] = z + gl->fz0;
 342    return;
 343 }
 344
 345 //________________________________________________________________________
 346 void AliITSgeom::LtoG(const Int_t *id,const Float_t *l,Float_t *g){
 347 ////////////////////////////////////////////////////////////////////////
 348 //     The function that does the local active volume detector Cartesian
 349 // coordinate to global ALICE Cartesian coordinate transformation.
 350 // The local detector coordinate system is determined by the three
 351 // element array Id containing as it's three elements Id[0]=layer,
 352 // Id[1]=ladder, and Id[2]=detector numbers. The local coordinates
 353 // are entered by the three element Float_t array l and the global
 354 // coordinate values are returned by the three element Float_t array g.
 355 // The order of the three elements are l[0]=x, l[1]=y, and l[2]=z,
 356 // similarly for g.
 357 ////////////////////////////////////////////////////////////////////////
 358    Int_t    lay,lad,det;
 359    Double_t x,y,z;
 360    ITS_geom *gl;
 361
 362    lay = id[0]; lad = id[1]; det = id[2];
 363    lay--; lad--; det--;
 364    gl   = &(fg[lay][fNdet[lay]*lad+det]);
 365
 366    x    = gl->fr[0]*l[0] + gl->fr[3]*l[1] + gl->fr[6]*l[2];
 367    y    = gl->fr[1]*l[0] + gl->fr[4]*l[1] + gl->fr[7]*l[2];
 368    z    = gl->fr[2]*l[0] + gl->fr[5]*l[1] + gl->fr[8]*l[2];
 369    g[0] = x + gl->fx0;
 370    g[1] = y + gl->fy0;
 371    g[2] = z + gl->fz0;
 372    return;
 373 }
 374 //________________________________________________________________________
 375 void AliITSgeom::LtoG(Int_t index,const Float_t *l,Float_t *g){
 376 ////////////////////////////////////////////////////////////////////////
 377 //     The function that does the local active volume detector Cartesian
 378 // coordinate to global ALICE Cartesian coordinate transformation.
 379 // The local detector coordinate system is determined by the detector
 380 // index number (see GetModuleIndex and GetModuleId). The local coordinates
 381 // are entered by the three element Float_t array l and the global
 382 // coordinate values are returned by the three element Float_t array g.
 383 // The order of the three elements are l[0]=x, l[1]=y, and l[2]=z,
 384 // similarly for g.
 385 ////////////////////////////////////////////////////////////////////////
 386    Int_t    lay,lad,det;
 387    Double_t x,y,z;
 388    ITS_geom *gl;
 389
 390    this->GetModuleId(index,lay,lad,det);
 391    lay--; lad--; det--;
 392    gl   = &(fg[lay][fNdet[lay]*lad+det]);
 393
 394    x    = gl->fr[0]*l[0] + gl->fr[3]*l[1] + gl->fr[6]*l[2];
 395    y    = gl->fr[1]*l[0] + gl->fr[4]*l[1] + gl->fr[7]*l[2];
 396    z    = gl->fr[2]*l[0] + gl->fr[5]*l[1] + gl->fr[8]*l[2];
 397    g[0] = x + gl->fx0;
 398    g[1] = y + gl->fy0;
 399    g[2] = z + gl->fz0;
 400    return;
 401 }
 402 //________________________________________________________________________
 403 void AliITSgeom::GtoLMomentum(Int_t lay,Int_t lad,Int_t det,
 404                               const Float_t *g,Float_t *l){
 405 ////////////////////////////////////////////////////////////////////////
 406 //     The function that does the global ALICE Cartesian momentum
 407 // to local active volume detector Cartesian momentum transformation.
 408 // The local detector coordinate system is determined by the layer,
 409 // ladder, and detector numbers. The global momentums are entered by
 410 // the three element Float_t array g and the local momentums values
 411 // are returned by the three element Float_t array l. The order of the
 412 // three elements are g[0]=x, g[1]=y, and g[2]=z, similarly for l.
 413 ////////////////////////////////////////////////////////////////////////
 414    Double_t px,py,pz;
 415    ITS_geom *gl;
 416
 417    lay--; lad--; det--;
 418    gl = &(fg[lay][fNdet[lay]*lad+det]);
 419
 420    px   = g[0];
 421    py   = g[1];
 422    pz   = g[2];
 423    l[0] = gl->fr[0]*px + gl->fr[1]*py + gl->fr[2]*pz;
 424    l[1] = gl->fr[3]*px + gl->fr[4]*py + gl->fr[5]*pz;
 425    l[2] = gl->fr[6]*px + gl->fr[7]*py + gl->fr[8]*pz;
 426    return;
 427 }
 428 //________________________________________________________________________
 429 void AliITSgeom::LtoGMomentum(Int_t lay,Int_t lad,Int_t det,
 430                               const Float_t *l,Float_t *g){
 431 ////////////////////////////////////////////////////////////////////////
 432 //     The function that does the local active volume detector Cartesian
 433 // momentum to global ALICE Cartesian momentum transformation.
 434 // The local detector momentum system is determined by the layer,
 435 // ladder, and detector numbers. The locall momentums are entered by
 436 // the three element Float_t array l and the global momentum values
 437 // are returned by the three element Float_t array g. The order of the
 438 // three elements are l[0]=x, l[1]=y, and l[2]=z, similarly for g.
 439 ////////////////////////////////////////////////////////////////////////
 440    Double_t px,py,pz;
 441    ITS_geom *gl;
 442
 443    lay--; lad--; det--;
 444    gl   = &(fg[lay][fNdet[lay]*lad+det]);
 445
 446    px   = gl->fr[0]*l[0] + gl->fr[3]*l[1] + gl->fr[6]*l[2];
 447    py   = gl->fr[1]*l[0] + gl->fr[4]*l[1] + gl->fr[7]*l[2];
 448    pz   = gl->fr[2]*l[0] + gl->fr[5]*l[1] + gl->fr[8]*l[2];
 449    g[0] = px;
 450    g[1] = py;
 451    g[2] = pz;
 452    return;
 453 }
 454 //___________________________________________________________________________
 455 Int_t AliITSgeom::GetModuleIndex(Int_t lay,Int_t lad,Int_t det){
 456     Int_t i,j,k;
 457
 458     i = fNdet[lay-1] * (lad-1) + det - 1;
 459     j = 0;
 460     for(k=0;k<lay-1;k++) j += fNdet[k]*fNlad[k];
 461     return (i+j);
 462 }
 463 //___________________________________________________________________________
 464 void AliITSgeom::GetModuleId(Int_t index,Int_t &lay,Int_t &lad,Int_t &det){
 465     Int_t i,j,k;
 466
 467     j = 0;
 468     for(k=0;k<fNlayers;k++){
 469         j += fNdet[k]*fNlad[k];
 470         if(index>j)break;
 471     } // end for k
 472     lay = k+1;
 473     i = index -j + fNdet[k]*fNlad[k];
 474     j = 0;
 475     for(k=0;k<fNlad[lay-1];k++){
 476         j += fNdet[k];
 477         if(i>fNdet[k])break;
 478     } // end for k
 479     lad = k+1;
 480     det = 1+i-fNdet[lay-1]*k;
 481     return;
 482 }
 483 //___________________________________________________________________________
 484 void AliITSgeom::GlobalChange(Float_t *tran,Float_t *rot){
 485 ////////////////////////////////////////////////////////////////////////
 486 //     This function performs a Cartesian translation and rotation of
 487 // the full ITS from its default position by an amount determined by
 488 // the three element arrays dtranslation and drotation. If every element
 489 // of dtranslation and drotation are zero then there is no change made
 490 // the geometry. The change is global in that the exact same translation
 491 // and rotation is done to every detector element in the exact same way.
 492 // The units of the translation are those of the Monte Carlo, usually cm,
 493 // and those of the rotation are in radians. The elements of dtranslation
 494 // are dtranslation[0] = x, dtranslation[1] = y, and dtranslation[2] = z.
 495 // The elements of drotation are drotation[0] = rx, drotation[1] = ry, and
 496 // drotation[2] = rz. A change in x will move the hole ITS in the ALICE
 497 // global x direction, the same for a change in y. A change in z will
 498 // result in a translation of the ITS as a hole up or down the beam line.
 499 // A change in the angles will result in the inclination of the ITS with
 500 // respect to the beam line, except for an effective rotation about the
 501 // beam axis which will just rotate the ITS as a hole about the beam axis.
 502 ////////////////////////////////////////////////////////////////////////
 503    Int_t    i,j,k,l;
 504    Double_t rx,ry,rz;
 505    Double_t sx,cx,sy,cy,sz,cz;
 506    ITS_geom *gl;
 507
 508    for(i=0;i<fNlayers;i++){
 509       for(j=0;j<fNlad[i];j++) for(k=0;k<fNdet[i];k++){
 510          l = fNdet[i]*j+k; // resolved index
 511          gl = &(fg[i][l]);
 512          gl->fx0 += tran[0];
 513          gl->fy0 += tran[1];
 514          gl->fz0 += tran[2];
 515          gl->frx +=  rot[0];
 516          gl->fry +=  rot[1];
 517          gl->frz +=  rot[2];
 518          rx = gl->frx; ry = gl->fry; rz = gl->frz;
 519          sx = sin(rx); cx = cos(rx);
 520          sy = sin(ry); cy = cos(ry);
 521          sz = sin(rz); cz = cos(rz);
 522          gl->fr[0] =  cz*cy;
 523          gl->fr[1] = -cz*sy*sx - sz*cx;
 524          gl->fr[2] = -cz*sy*cx + sz*sx;
 525          gl->fr[3] =  sz*cy;
 526          gl->fr[4] = -sz*sy*sx + cz*cx;
 527          gl->fr[5] = -sz*sy*cx - cz*sx;
 528          gl->fr[6] =  sy;
 529          gl->fr[7] =  cy*sx;
 530          gl->fr[8] =  cy*cx;
 531       } // end for j,k
 532    } // end for i
 533    return;
 534 }
 535
 536 //___________________________________________________________________________
 537 void AliITSgeom::GlobalCylindericalChange(Float_t *tran,Float_t *rot){
 538 ////////////////////////////////////////////////////////////////////////
 539 //     This function performs a cylindrical translation and rotation of
 540 // each ITS element by a fixed about in radius, rphi, and z from its
 541 // default position by an amount determined by the three element arrays
 542 // dtranslation and drotation. If every element of dtranslation and
 543 // drotation are zero then there is no change made the geometry. The
 544 // change is global in that the exact same distance change in translation
 545 // and rotation is done to every detector element in the exact same way.
 546 // The units of the translation are those of the Monte Carlo, usually cm,
 547 // and those of the rotation are in radians. The elements of dtranslation
 548 // are dtranslation[0] = r, dtranslation[1] = rphi, and dtranslation[2] = z.
 549 // The elements of drotation are drotation[0] = rx, drotation[1] = ry, and
 550 // drotation[2] = rz. A change in r will results in the increase of the
 551 // radius of each layer by the same about. A change in rphi will results in
 552 // the rotation of each layer by a different angle but by the same
 553 // circumferential distance. A change in z will result in a translation
 554 // of the ITS as a hole up or down the beam line. A change in the angles
 555 // will result in the inclination of the ITS with respect to the beam
 556 // line, except for an effective rotation about the beam axis which will
 557 // just rotate the ITS as a hole about the beam axis.
 558 ////////////////////////////////////////////////////////////////////////
 559    Int_t    i,j,k,l;
 560    Double_t rx,ry,rz,r,phi,rphi; // phi in radians
 561    Double_t sx,cx,sy,cy,sz,cz,r0;
 562    ITS_geom *gl;
 563
 564 //   printf("trans=%f %f %f rot=%f %f %f\n",tran[0],tran[1],tran[2],
 565 //        rot[0],rot[1],rot[2]);
 566    for(i=0;i<fNlayers;i++){
 567       for(j=0;j<fNlad[i];j++) for(k=0;k<fNdet[i];k++){
 568          l     = fNdet[i]*j+k; // resolved index
 569          gl    = &(fg[i][l]);
 570          r = r0= TMath::Hypot(gl->fy0,gl->fx0);
 571          phi   = atan2(gl->fy0,gl->fx0);
 572          rphi  = r0*phi;
 573          r    += tran[0];
 574          rphi += tran[1];
 575          phi   = rphi/r0;
 576          gl->fx0  = r*TMath::Cos(phi);
 577          gl->fy0  = r*TMath::Sin(phi);
 578          gl->fz0 += tran[2];
 579          gl->frx +=  rot[0];
 580          gl->fry +=  rot[1];
 581          gl->frz +=  rot[2];
 582          rx = gl->frx; ry = gl->fry; rz = gl->frz;
 583          sx = sin(rx); cx = cos(rx);
 584          sy = sin(ry); cy = cos(ry);
 585          sz = sin(rz); cz = cos(rz);
 586          gl->fr[0] =  cz*cy;
 587          gl->fr[1] = -cz*sy*sx - sz*cx;
 588          gl->fr[2] = -cz*sy*cx + sz*sx;
 589          gl->fr[3] =  sz*cy;
 590          gl->fr[4] = -sz*sy*sx + cz*cx;
 591          gl->fr[5] = -sz*sy*cx - cz*sx;
 592          gl->fr[6] =  sy;
 593          gl->fr[7] =  cy*sx;
 594          gl->fr[8] =  cy*cx;
 595       } // end for j,k
 596    } // end for i
 597    return;
 598 }
 599
 600 //___________________________________________________________________________
 601 void AliITSgeom::RandomChange(Float_t *stran,Float_t *srot){
 602 ////////////////////////////////////////////////////////////////////////
 603 //     This function performs a Gaussian random displacement and/or
 604 // rotation about the present global position of each active
 605 // volume/detector of the ITS. The sigma of the random displacement
 606 // is determined by the three element array stranslation, for the
 607 // x y and z translations, and the three element array srotation,
 608 // for the three rotation about the axis x y and z.
 609 ////////////////////////////////////////////////////////////////////////
 610    Int_t    i,j,k,l;
 611    Double_t rx,ry,rz;
 612    Double_t sx,cx,sy,cy,sz,cz;
 613    TRandom  ran;
 614    ITS_geom *gl;
 615
 616    for(i=0;i<fNlayers;i++){
 617       for(j=0;j<fNlad[i];j++) for(k=0;k<fNdet[i];k++){
 618          l = fNdet[i]*j+k; // resolved index
 619          gl = &(fg[i][l]);
 620          gl->fx0 += ran.Gaus(0.0,stran[0]);
 621          gl->fy0 += ran.Gaus(0.0,stran[1]);
 622          gl->fz0 += ran.Gaus(0.0,stran[2]);
 623          gl->frx += ran.Gaus(0.0, srot[0]);
 624          gl->fry += ran.Gaus(0.0, srot[1]);
 625          gl->frz += ran.Gaus(0.0, srot[2]);
 626          rx = gl->frx; ry = gl->fry; rz = gl->frz;
 627          sx = sin(rx); cx = cos(rx);
 628          sy = sin(ry); cy = cos(ry);
 629          sz = sin(rz); cz = cos(rz);
 630          gl->fr[0] =  cz*cy;
 631          gl->fr[1] = -cz*sy*sx - sz*cx;
 632          gl->fr[2] = -cz*sy*cx + sz*sx;
 633          gl->fr[3] =  sz*cy;
 634          gl->fr[4] = -sz*sy*sx + cz*cx;
 635          gl->fr[5] = -sz*sy*cx - cz*sx;
 636          gl->fr[6] =  sy;
 637          gl->fr[7] =  cy*sx;
 638          gl->fr[8] =  cy*cx;
 639       } // end for j,k
 640    } // end for i
 641    return;
 642 }
 643
 644 //___________________________________________________________________________
 645 void AliITSgeom::RandomCylindericalChange(Float_t *stran,Float_t *srot){
 646 ////////////////////////////////////////////////////////////////////////
 647 //     This function performs a Gaussian random displacement and/or
 648 // rotation about the present global position of each active
 649 // volume/detector of the ITS. The sigma of the random displacement
 650 // is determined by the three element array stranslation, for the
 651 // r rphi and z translations, and the three element array srotation,
 652 // for the three rotation about the axis x y and z. This random change
 653 // in detector position allow for the simulation of a random uncertainty
 654 // in the detector positions of the ITS.
 655 ////////////////////////////////////////////////////////////////////////
 656    Int_t     i,j,k,l;
 657    Double_t  rx,ry,rz,r,phi,x,y;  // phi in radians
 658    Double_t  sx,cx,sy,cy,sz,cz,r0;
 659    TRandom   ran;
 660    ITS_geom  *gl;
 661
 662 //   printf("trans=%f %f %f rot=%f %f %f\n",stran[0],stran[1],stran[2],
 663 //        srot[0],srot[1],srot[2]);
 664    for(i=0;i<fNlayers;i++){
 665       for(j=0;j<fNlad[i];j++) for(k=0;k<fNdet[i];k++){
 666          l     = fNdet[i]*j+k; // resolved index
 667          gl    = &(fg[i][l]);
 668          x     = gl->fx0;
 669          y     = gl->fy0;
 670          r = r0= TMath::Hypot(y,x);
 671          phi   = TMath::ATan2(y,x);
 672 //       if(phi<0.0) phi += 2.0*TMath::Pi();
 673          r    += ran.Gaus(0.0,stran[0]);
 674          phi  += ran.Gaus(0.0,stran[1])/r0;
 675 //       printf("fx0=%f fy0=%f rcos(phi)=%f rsin(phi)=%f\n",gl->fx0,gl->fy0,
 676 //              r*TMath::Cos(phi),r*TMath::Sin(phi));
 677          gl->fx0  = r*TMath::Cos(phi);
 678          gl->fy0  = r*TMath::Sin(phi);
 679 //       printf("r0=%f r=%f hypot=%f phi0=%f phi=%f ATan2=%f\n",
 680 //              r0,r,TMath::Hypot(gl->fy0,gl->fx0),
 681 //              phi0,phi,TMath::ATan2(gl->fy0,gl->fx0));
 682          gl->fz0 += ran.Gaus(0.0,stran[2]);
 683          gl->frx += ran.Gaus(0.0, srot[0]);
 684          gl->fry += ran.Gaus(0.0, srot[1]);
 685          gl->frz += ran.Gaus(0.0, srot[2]);
 686          rx = gl->frx; ry = gl->fry; rz = gl->frz;
 687          sx = sin(rx); cx = cos(rx);
 688          sy = sin(ry); cy = cos(ry);
 689          sz = sin(rz); cz = cos(rz);
 690          gl->fr[0] =  cz*cy;
 691          gl->fr[1] = -cz*sy*sx - sz*cx;
 692          gl->fr[2] = -cz*sy*cx + sz*sx;
 693          gl->fr[3] =  sz*cy;
 694          gl->fr[4] = -sz*sy*sx + cz*cx;
 695          gl->fr[5] = -sz*sy*cx - cz*sx;
 696          gl->fr[6] =  sy;
 697          gl->fr[7] =  cy*sx;
 698          gl->fr[8] =  cy*cx;
 699       } // end for j,k
 700    } // end for i
 701    return;
 702 }
 703
 704 //___________________________________________________________________________
 705 void AliITSgeom::SetByAngles(Int_t lay,Int_t lad,Int_t det,
 706                              Float_t rx,Float_t ry,Float_t rz){
 707 ////////////////////////////////////////////////////////////////////////
 708 //     This function computes a new rotation matrix based on the angles
 709 // rx, ry, and rz (in radians) for a give detector on the give ladder
 710 // in the give layer. A new
 711 // fg[layer-1][(fNlad[layer-1]*(ladder-1)+detector-1)].fr[] array is
 712 // computed.
 713 ////////////////////////////////////////////////////////////////////////
 714    ITS_geom *g;
 715    Double_t  sx,cx,sy,cy,sz,cz;
 716
 717    lay--; lad--; det--; // set to zero base now.
 718    g = &(fg[lay][fNdet[lay]*lad+det]);
 719
 720    sx = sin(rx); cx = cos(rx);
 721    sy = sin(ry); cy = cos(ry);
 722    sz = sin(rz); cz = cos(rz);
 723    g->frx   = rx;
 724    g->fry   = ry;
 725    g->frz   = rz;
 726    g->fr[0] =  cz*cy;
 727    g->fr[1] = -cz*sy*sx - sz*cx;
 728    g->fr[2] = -cz*sy*cx + sz*sx;
 729    g->fr[3] =  sz*cy;
 730    g->fr[4] = -sz*sy*sx + cz*cx;
 731    g->fr[5] = -sz*sy*cx - cz*sx;
 732    g->fr[6] =  sy;
 733    g->fr[7] =  cy*sx;
 734    g->fr[8] =  cy*cx;
 735    return;
 736 }
 737
 738 //___________________________________________________________________________
 739 void AliITSgeom::GetRotMatrix(Int_t lay,Int_t lad,Int_t det,Float_t *mat){
 740 ////////////////////////////////////////////////////////////////////////
 741 //     Returns, in the Float_t array pointed to by mat, the full rotation
 742 // matrix for the give detector defined by layer, ladder, and detector.
 743 // It returns all nine elements of fr in the ITS_geom structure. See the
 744 // description of the ITS_geom structure for further details of this
 745 // rotation matrix.
 746 ////////////////////////////////////////////////////////////////////////
 747    Int_t    i;
 748    ITS_geom *g;
 749
 750    lay--; lad--; det--; // shift to base 0
 751    g = &(fg[lay][fNdet[lay]*lad+det]);
 752    for(i=0;i<9;i++) mat[i] = g->fr[i];
 753    return;
 754 }
 755
 756 //___________________________________________________________________________
 757 void AliITSgeom::PrintComparison(FILE *fp,AliITSgeom *other){
 758 ////////////////////////////////////////////////////////////////////////
 759 //     This function was primarily created for diagnostic reasons. It
 760 // print to a file pointed to by the file pointer fp the difference
 761 // between two AliITSgeom classes. The format of the file is basicly,
 762 // define d? to be the difference between the same element of the two
 763 // classes. For example dfrx = this->fg[i][j].frx - other->fg[i][j].frx.
 764 // if(at least one of dfx0, dfy0, dfz0,dfrx,dfry,dfrz are non zero) then print
 765 // layer ladder detector dfx0 dfy0 dfz0 dfrx dfry dfrz
 766 // if(at least one of the 9 elements of dfr[] are non zero) then print
 767 // layer ladder detector dfr[0] dfr[1] dfr[2]
 768 //                       dfr[3] dfr[4] dfr[5]
 769 //                       dfr[6] dfr[7] dfr[8]
 770 // Only non zero values are printed to save space. The differences are
 771 // typical written to a file because there are usually a lot of numbers
 772 // printed out and it is usually easier to read them in some nice editor
 773 // rather than zooming quickly past you on a screen. fprintf is used to
 774 // do the printing. The fShapeIndex difference is not printed at this time.
 775 ////////////////////////////////////////////////////////////////////////
 776    Int_t    i,j,k,l;
 777    Double_t xt,yt,zt,xo,yo,zo;
 778    Double_t rxt,ryt,rzt,rxo,ryo,rzo;  // phi in radians
 779    ITS_geom *gt,*go;
 780    Bool_t   t;
 781
 782    for(i=0;i<this->fNlayers;i++){
 783       for(j=0;j<this->fNlad[i];j++) for(k=0;k<this->fNdet[i];k++){
 784          l   = this->fNdet[i]*j+k; // resolved index
 785          gt  = &(this->fg[i][l]);
 786          go  = &(other->fg[i][l]);
 787          xt  = gt->fx0; yt  = gt->fy0; zt  = gt->fz0;
 788          xo  = go->fx0; yo  = go->fy0; zo  = go->fz0;
 789          rxt = gt->frx; ryt = gt->fry; rzt = gt->frz;
 790          rxo = go->frx; ryo = go->fry; rzo = go->frz;
 791          if(!(xt==xo&&yt==yo&&zt==zo&&rxt==rxo&&ryt==ryo&&rzt==rzo))
 792          fprintf(fp,"%1.1d %2.2d %2.2d dTrans=%f %f %f drot=%f %f %f\n",
 793                  i+1,j+1,k+1,xt-xo,yt-yo,zt-zo,rxt-rxo,ryt-ryo,rzt-rzo);
 794          t = kFALSE;
 795          for(i=0;i<9;i++) t = gt->fr[i] != go->fr[i];
 796          if(t){
 797              fprintf(fp,"%1.1d %2.2d %2.2d dfr= %e %e %e\n",i+1,j+1,k+1,
 798                  gt->fr[0]-go->fr[0],gt->fr[1]-go->fr[1],gt->fr[2]-go->fr[2]);
 799              fprintf(fp,"        dfr= %e %e %e\n",
 800                  gt->fr[3]-go->fr[3],gt->fr[4]-go->fr[4],gt->fr[5]-go->fr[5]);
 801              fprintf(fp,"        dfr= %e %e %e\n",
 802                  gt->fr[6]-go->fr[6],gt->fr[7]-go->fr[7],gt->fr[8]-go->fr[8]);
 803          }
 804       } // end for j,k
 805    } // end for i
 806    return;
 807 }
 808
 809 //___________________________________________________________________________
 810 void AliITSgeom::PrintData(FILE *fp,Int_t lay,Int_t lad,Int_t det){
 811 ////////////////////////////////////////////////////////////////////////
 812 //     This function prints out the coordinate transformations for
 813 // the particular detector defined by layer, ladder, and detector
 814 // to the file pointed to by the File pointer fp. fprinf statements
 815 // are used to print out the numbers. The format is
 816 // layer ladder detector Trans= fx0 fy0 fz0 rot= frx fry frz Shape=fShapeIndex
 817 //                         dfr= fr[0] fr[1] fr[2]
 818 //                         dfr= fr[3] fr[4] fr[5]
 819 //                         dfr= fr[6] fr[7] fr[8]
 820 // By indicating which detector, some control over the information
 821 // is given to the user. The output it written to the file pointed
 822 // to by the file pointer fp. This can be set to stdout if you want.
 823 ////////////////////////////////////////////////////////////////////////
 824    Int_t    i,j,k,l;
 825    ITS_geom *gt;
 826
 827    i  = lay-1;
 828    j  = lad-1;
 829    k  = det-1;
 830    l  = this->fNdet[i]*j+k; // resolved index
 831    gt = &(this->fg[i][l]);
 832    fprintf(fp,"%1.1d %2.2d %2.2d Trans=%f %f %f rot=%f %f %f Shape=%d\n",
 833            i+1,j+1,k+1,gt->fx0,gt->fy0,gt->fz0,gt->frx,gt->fry,gt->frz,
 834            gt->fShapeIndex);
 835    fprintf(fp,"        dfr= %e %e %e\n",gt->fr[0],gt->fr[1],gt->fr[2]);
 836    fprintf(fp,"        dfr= %e %e %e\n",gt->fr[3],gt->fr[4],gt->fr[5]);
 837    fprintf(fp,"        dfr= %e %e %e\n",gt->fr[6],gt->fr[7],gt->fr[8]);
 838    return;
 839 }
 840 //___________________________________________________________________________
 841 void AliITSgeom::Streamer(TBuffer &R__b){
 842 ////////////////////////////////////////////////////////////////////////
 843 //     The default Streamer function "written by ROOT" doesn't write out
 844 // the arrays referenced by pointers. Therefore, a specific Streamer function
 845 // has to be written. This function should not be modified but instead added
 846 // on to so that older versions can still be read. The proper handling of
 847 // the version dependent streamer function hasn't been written do to the lack
 848 // of finding an example at the time of writting.
 849 ////////////////////////////////////////////////////////////////////////
 850    // Stream an object of class AliITSgeom.
 851     Int_t i,j,k;
 852
 853    if (R__b.IsReading()) {
 854       Version_t R__v = R__b.ReadVersion(); if (R__v) { }
 855       TObject::Streamer(R__b);
 856       R__b >> fNlayers;
 857       if(fNlad!=0) delete[] fNlad;
 858       if(fNdet!=0) delete[] fNdet;
 859       fNlad = new Int_t[fNlayers];
 860       fNdet = new Int_t[fNlayers];
 861       for(i=0;i<fNlayers;i++) R__b >> fNlad[i];
 862       for(i=0;i<fNlayers;i++) R__b >> fNdet[i];
 863       if(fg!=0){
 864           for(i=0;i<fNlayers;i++) delete[] fg[i];
 865           delete[] fg;
 866       } // end if fg!=0
 867       fg = new ITS_geom*[fNlayers];
 868       for(i=0;i<fNlayers;i++){
 869           fg[i] = new ITS_geom[fNlad[i]*fNdet[i]];
 870           for(j=0;j<fNlad[i]*fNdet[i];j++){
 871               R__b >> fg[i][j].fShapeIndex;
 872               R__b >> fg[i][j].fx0;
 873               R__b >> fg[i][j].fy0;
 874               R__b >> fg[i][j].fz0;
 875               R__b >> fg[i][j].frx;
 876               R__b >> fg[i][j].fry;
 877               R__b >> fg[i][j].frz;
 878               for(k=0;k<9;k++) R__b >> fg[i][j].fr[k];
 879           } // end for j
 880       } // end for i
 881       R__b >> fShape;
 882    } else {
 883       R__b.WriteVersion(AliITSgeom::IsA());
 884       TObject::Streamer(R__b);
 885       R__b << fNlayers;
 886       for(i=0;i<fNlayers;i++) R__b << fNlad[i];
 887       for(i=0;i<fNlayers;i++) R__b << fNdet[i];
 888       for(i=0;i<fNlayers;i++) for(j=0;j<fNlad[i]*fNdet[i];j++){
 889           R__b << fg[i][j].fShapeIndex;
 890           R__b << fg[i][j].fx0;
 891           R__b << fg[i][j].fy0;
 892           R__b << fg[i][j].fz0;
 893           R__b << fg[i][j].frx;
 894           R__b << fg[i][j].fry;
 895           R__b << fg[i][j].frz;
 896           for(k=0;k<9;k++) R__b << fg[i][j].fr[k];
 897       } // end for i,j
 898       R__b << fShape;
 899    }
 900 }
 901
 902 //___________________________________________________________________________
 903 ofstream & AliITSgeom::PrintGeom(ofstream &R__b){
 904 ////////////////////////////////////////////////////////////////////////
 905 //     The default Streamer function "written by ROOT" doesn't write out
 906 // the arrays referenced by pointers. Therefore, a specific Streamer function
 907 // has to be written. This function should not be modified but instead added
 908 // on to so that older versions can still be read. The proper handling of
 909 // the version dependent streamer function hasn't been written do to the lack
 910 // of finding an example at the time of writting.
 911 ////////////////////////////////////////////////////////////////////////
 912    // Stream an object of class AliITSgeom.
 913     Int_t i,j,k;
 914
 915     R__b.setf(ios::scientific);
 916     R__b << fNlayers << " ";
 917     for(i=0;i<fNlayers;i++) R__b << fNlad[i] << " ";
 918     for(i=0;i<fNlayers;i++) R__b << fNdet[i] << "\n";
 919     for(i=0;i<fNlayers;i++) for(j=0;j<fNlad[i]*fNdet[i];j++){
 920         R__b <<setprecision(16) << fg[i][j].fShapeIndex << " ";
 921         R__b <<setprecision(16) << fg[i][j].fx0 << " ";
 922         R__b <<setprecision(16) << fg[i][j].fy0 << " ";
 923         R__b <<setprecision(16) << fg[i][j].fz0 << " ";
 924         R__b <<setprecision(16) << fg[i][j].frx << " ";
 925         R__b <<setprecision(16) << fg[i][j].fry << " ";
 926         R__b <<setprecision(16) << fg[i][j].frz << "\n";
 927         for(k=0;k<9;k++) R__b <<setprecision(16) << fg[i][j].fr[k] << " ";
 928         R__b << "\n";
 929       } // end for i,j
 930 //      R__b << fShape;
 931       return R__b;
 932 }
 933
 934 //___________________________________________________________________________
 935 ifstream & AliITSgeom::ReadGeom(ifstream &R__b){
 936 ////////////////////////////////////////////////////////////////////////
 937 //     The default Streamer function "written by ROOT" doesn't write out
 938 // the arrays referenced by pointers. Therefore, a specific Streamer function
 939 // has to be written. This function should not be modified but instead added
 940 // on to so that older versions can still be read. The proper handling of
 941 // the version dependent streamer function hasn't been written do to the lack
 942 // of finding an example at the time of writting.
 943 ////////////////////////////////////////////////////////////////////////
 944    // Stream an object of class AliITSgeom.
 945     Int_t i,j,k;
 946
 947       R__b >> fNlayers;
 948       if(fNlad!=0) delete[] fNlad;
 949       if(fNdet!=0) delete[] fNdet;
 950       fNlad = new Int_t[fNlayers];
 951       fNdet = new Int_t[fNlayers];
 952       for(i=0;i<fNlayers;i++) R__b >> fNlad[i];
 953       for(i=0;i<fNlayers;i++) R__b >> fNdet[i];
 954       if(fg!=0){
 955           for(i=0;i<fNlayers;i++) delete[] fg[i];
 956           delete[] fg;
 957       } // end if fg!=0
 958       fg = new ITS_geom*[fNlayers];
 959       for(i=0;i<fNlayers;i++){
 960           fg[i] = new ITS_geom[fNlad[i]*fNdet[i]];
 961           for(j=0;j<fNlad[i]*fNdet[i];j++){
 962               R__b >> fg[i][j].fShapeIndex;
 963               R__b >> fg[i][j].fx0;
 964               R__b >> fg[i][j].fy0;
 965               R__b >> fg[i][j].fz0;
 966               R__b >> fg[i][j].frx;
 967               R__b >> fg[i][j].fry;
 968               R__b >> fg[i][j].frz;
 969               for(k=0;k<9;k++) R__b >> fg[i][j].fr[k];
 970           } // end for j
 971       } // end for i
 972 //      R__b >> fShape;
 973       return R__b;
 974 }