ITS/AliITSgeom.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 $Log$
  18 Revision 1.4.4.15  2000/10/02 15:52:05  barbera
  19 Forward declaration added
  20
  21 Revision 1.10  2000/09/05 14:25:50  nilsen
  22 Made fixes for HP compiler. All function parameter default values placed
  23 in .h file. Fixed the usual problem with HP comilers and the "for(Int_t i..."
  24 business. Replaced casting (Double_t [3][3]) to (Double_t (*)[3]) for HP.
  25 Lastly removed all "const" before function parameters which were 2 dim. arrays,
  26 because on HP root generates some strange code (?). Thanks Peter for the
  27 changes.
  28
  29 Revision 1.9  2000/08/29 20:19:03  nilsen
  30 Removed dependancy on structure AliITSeomS and replaced it with class
  31 AliITSgeomMatrix. Added many new functions with many new arguments. Most
  32 in the form of in line functions for speed.
  33
  34 Revision 1.4.4.6  2000/06/04 16:33:32  Nilsen
  35 A restructured AliITSgeom class. Now used AliITSgeomMatrix.
  36
  37 Revision 1.4.4.5  2000/03/04 23:42:39  Nilsen
  38 Updated the comments/documentations and improved the maintainability of the
  39 code.
  40
  41 Revision 1.4.4.4  2000/03/02 21:27:07  Nilsen
  42 Added two functions, SetByAngles and SetTrans.
  43
  44 Revision 1.4.4.3  2000/01/23 03:09:10  Nilsen
  45 // fixed compiler warnings for new function LtLErrorMatrix(...)
  46
  47 Revision 1.4.4.2  2000/01/19 23:18:20  Nilsen
  48 Added transformations of Error matrix to AliITSgeom and fixed some typos
  49 in AliITS.h and AliITShitIndex.h
  50
  51 Revision 1.4.4.1  2000/01/12 19:03:32  Nilsen
  52 This is the version of the files after the merging done in December 1999.
  53 See the ReadMe110100.txt file for details
  54
  55 Revision 1.4  1999/10/15 07:03:20  fca
  56 Fixed bug in GetModuleId(Int_t index,Int_t &lay,Int_t &lad, Int_t &det) and
  57 a typo in the creator. aliroot need to be rerun to get a fixed geometry.
  58
  59 Revision 1.3  1999/10/04 15:20:12  fca
  60 Correct syntax accepted by g++ but not standard for static members, remove minor warnings
  61
  62 Revision 1.2  1999/09/29 09:24:20  fca
  63 Introduction of the Copyright and cvs Log
  64
  65 */
  66
  67 ///////////////////////////////////////////////////////////////////////
  68 // ITS geometry manipulation routines.                               //
  69 // Created April 15 1999.                                            //
  70 // version: 0.0.0                                                    //
  71 // By: Bjorn S. Nilsen                                               //
  72 // version: 0.0.1                                                    //
  73 // Updated May 27 1999.                                              //
  74 // Added Cylindrical random and global based changes.               //
  75 // Added  function PrintComparison.                                  //
  76 ///////////////////////////////////////////////////////////////////////
  77
  78
  79 ////////////////////////////////////////////////////////////////////////
  80 //     The local coordinate system by, default, is show in the following
  81 // figures. Also shown are the ladder numbering scheme.
  82 //Begin_Html
  83 /*
  84 <img src="picts/ITS/AliITSgeomMatrix_L1.gif">
  85 </pre>
  86 <br clear=left>
  87 <font size=+2 color=blue>
  88 <p>This shows the relative geometry differences between the ALICE Global
  89 coordinate system and the local detector coordinate system.
  90 </font>
  91 <pre>
  92
  93 <pre>
  94 <img src="picts/ITS/its1+2_convention_front_5.gif">
  95 </pre>
  96 <br clear=left>
  97 <font size=+2 color=blue>
  98 <p>This shows the front view of the SPDs and the orientation of the local
  99 pixel coordinate system. Note that the inner pixel layer has its y coordinate
 100 in the opposite direction from all of the other layers.
 101 </font>
 102 <pre>
 103
 104 <pre>
 105 <img src="picts/ITS/its3+4_convention_front_5.gif">
 106 </pre>
 107 <br clear=left>
 108 <font size=+2 color=blue>
 109 <p>This shows the front view of the SDDs and the orientation of the local
 110 pixel coordinate system.
 111 </font>
 112 <pre>
 113
 114 <pre>
 115 <img src="picts/ITS/its5+6_convention_front_5.gif">
 116 </pre>
 117 <br clear=left>
 118 <font size=+2 color=blue>
 119 <p>This shows the front view of the SSDs and the orientation of the local
 120 pixel coordinate system.
 121 </font>
 122 <pre>
 123 */
 124 //End_Html
 125 //
 126 ////////////////////////////////////////////////////////////////////////
 127
 128 ////////////////////////////////////////////////////////////////////////
 129 //
 130 // version: 0
 131 // Written by Bjorn S. Nilsen
 132 //
 133 // Data Members:
 134 //
 135 // Int_t fNlayers
 136 //     The number of ITS layers for this geometry. By default this
 137 //  is 6, but can be modified by the creator function if there are
 138 // more layers defined.
 139 //
 140 // Int_t *fNlad
 141 //     A pointer to an array fNlayers long containing the number of
 142 // ladders for each layer. This array is typically created and filled
 143 // by the AliITSgeom creator function.
 144 //
 145 // Int_t *fNdet
 146 //     A pointer to an array fNlayers long containing the number of
 147 // active detector volumes for each ladder. This array is typically
 148 // created and filled by the AliITSgeom creator function.
 149 //
 150 // AliITSgeomMatrix *fGm
 151 //     A pointer to an array of AliITSgeomMatrix classes. One element
 152 // per module (detector) in the ITS. AliITSgeomMatrix basicly contains
 153 // all of the necessary information about the detector and it's coordinate
 154 // transformations.
 155 //
 156 // TObjArray *fShape
 157 //     A pointer to an array of TObjects containing the detailed shape
 158 // information for each type of detector used in the ITS. For example
 159 // I have created AliITSgeomSPD, AliITSgeomSDD, and AliITSgeomSSD as
 160 // example structures, derived from TObjects, to hold the detector
 161 // information. I would recommend that one element in each of these
 162 // structures, that which describes the shape of the active volume,
 163 // be one of the ROOT classes derived from TShape. In this way it would
 164 // be easy to have the display program display the correct active
 165 // ITS volumes. See the example classes AliITSgeomSPD, AliITSgeomSDD,
 166 // and AliITSgeomSSD for a more detailed example.
 167 ////////////////////////////////////////////////////////////////////////
 168 #include <iostream.h>
 169 #include <fstream.h>
 170 #include <iomanip.h>
 171 #include <stdio.h>
 172 #include <string.h>
 173 #include <ctype.h>
 174 #include <TObject.h>
 175 #include <TRandom.h>
 176
 177 #include "AliITSgeom.h"
 178 #include "AliITSgeomMatrix.h"
 179 #include "AliITSgeomSPD.h"
 180 #include "AliITSgeomSDD.h"
 181 #include "AliITSgeomSSD.h"
 182
 183 ClassImp(AliITSgeom)
 184
 185 //_____________________________________________________________________
 186 AliITSgeom::AliITSgeom(){
 187 ////////////////////////////////////////////////////////////////////////
 188 //     The default constructor for the AliITSgeom class. It, by default,
 189 // sets fNlayers to zero and zeros all pointers.
 190 ////////////////////////////////////////////////////////////////////////
 191   // Default constructor.
 192   // Do not allocate anything zero everything
 193    fTrans   = 0; // standard GEANT global/local coordinate system.
 194    fNlayers = 0;
 195    fNlad    = 0;
 196    fNdet    = 0;
 197    fGm      = 0;
 198    fShape   = 0;
 199    return;
 200 }
 201
 202 //_____________________________________________________________________
 203 AliITSgeom::~AliITSgeom(){
 204 ////////////////////////////////////////////////////////////////////////
 205 //     The destructor for the AliITSgeom class. If the arrays fNlad,
 206 // fNdet, or fGm have had memory allocated to them, there pointer values
 207 // are non zero, then this memory space is freed and they are set
 208 // to zero. In addition, fNlayers is set to zero. The destruction of
 209 // TObjArray fShape is, by default, handled by the TObjArray destructor.
 210 ////////////////////////////////////////////////////////////////////////
 211   // Default destructor.
 212   // if arrays exist delete them. Then set everything to zero.
 213    if(fGm!=0){
 214       for(Int_t i=0;i<fNlayers;i++) delete fGm[i];
 215       delete fGm;
 216    } // end if fGm!=0
 217    if(fNlad!=0) delete[] fNlad;
 218    if(fNdet!=0) delete[] fNdet;
 219    fNlayers = 0;
 220    fNlad    = 0;
 221    fNdet    = 0;
 222    fGm      = 0;
 223    return;
 224 }
 225 //______________________________________________________________________
 226 void AliITSgeom::ReadNewFile(const char *filename){
 227     printf("New file format not defined yet\n");
 228     return;
 229 }
 230 //_____________________________________________________________________
 231 AliITSgeom::AliITSgeom(const char *filename){
 232 ////////////////////////////////////////////////////////////////////////
 233 //     The constructor for the AliITSgeom class. All of the data to fill
 234 // this structure is read in from the file given my the input filename.
 235 ////////////////////////////////////////////////////////////////////////
 236    FILE     *pf;
 237    Int_t    i,lm=0,id[3];
 238    Int_t    l,a,d;
 239    Float_t  x,y,z,o,p,q,r,s,t;
 240    Double_t rot6[6],tran[3];
 241    char     buf[200],*buff=0; // input character buffer;
 242
 243    pf = fopen(filename,"r");
 244
 245    fNlayers = 6; // set default number of ladders
 246 TryAgain:
 247    fNlad    = new Int_t[fNlayers];
 248    fNdet    = new Int_t[fNlayers];
 249    fNmodules = 0;
 250    // find the number of ladders and detectors in this geometry.
 251    for(i=0;i<fNlayers;i++){fNlad[i]=fNdet[i]=0;} // zero out arrays
 252    while(fgets(buf,200,pf)!=NULL){ // for ever loop
 253       for(i=0;i<200;i++)if(buf[i]!=' '){ // remove blank spaces.
 254            buff = &(buf[i]);
 255            break;
 256       } // end for i
 257       // remove blank lines and comments.
 258       if(buff[0]=='\n'||buff[0]=='#'||buff[0]=='!'||
 259          (buff[0]=='/'&&buff[1]=='/')) continue;
 260       if(isalpha(buff[0])) { // must be the new file formated file.
 261             fclose(pf);
 262             delete[] fNlad;delete[] fNdet;
 263             ReadNewFile(filename);
 264             return;
 265       } // end if isalpha(buff[0])
 266       sscanf(buff,"%d %d %d %f %f %f %f %f %f %f %f %f",
 267                   &l,&a,&d,&x,&y,&z,&o,&p,&q,&r,&s,&t);
 268       if(l>lm) lm = l;
 269       if(l<1 || l>fNlayers) {
 270          printf("error in file %s layer=%d min. is 1 max is %d/n",
 271                  filename,l,fNlayers);
 272          continue;
 273       }// end if l
 274       fNmodules++;
 275       if(l<=fNlayers&&fNlad[l-1]<a) fNlad[l-1] = a;
 276       if(l<=fNlayers&&fNdet[l-1]<d) fNdet[l-1] = d;
 277    } // end while ever loop
 278    if(lm>fNlayers){
 279         delete[] fNlad;
 280         delete[] fNdet;
 281         fNlayers = lm;
 282         goto TryAgain;
 283    } // end if lm>fNlayers
 284    // counted the number of ladders and detectors now allocate space.
 285    fGm = new AliITSgeomMatrix*[fNmodules];
 286
 287    // Set up Shapes for a default configuration of 6 layers.
 288    fTrans   = 0; // standard GEANT global/local coordinate system.
 289    fShape = new TObjArray(3);
 290    AddShape((TObject *) new AliITSgeomSPD());  // shape 0
 291    AddShape((TObject *) new AliITSgeomSDD());  // shape 1
 292    AddShape((TObject *) new AliITSgeomSPD());  // shape 2
 293
 294    // prepare to read in transforms
 295    lm = 0; // reuse lm as counter of modules.
 296    rewind(pf); // start over reading file
 297    while(fgets(buf,200,pf)!=NULL){ // for ever loop
 298       for(i=0;i<200;i++)if(buf[i]!=' '){ // remove blank spaces.
 299            buff = &(buf[i]);
 300            break;
 301       } // end for i
 302       // remove blank lines and comments.
 303       if(buff[0]=='\n'||buff[0]=='#'||buff[0]=='!'||
 304         (buff[0]=='/'&&buff[1]=='/')) continue;
 305       x = y = z = o = p = q = r = s = t = 0.0;
 306       sscanf(buff,"%d %d %d %f %f %f %f %f %f %f %f %f",
 307                   &l,&a,&d,&x,&y,&z,&o,&p,&q,&r,&s,&t);
 308       if(l<1 || l>fNlayers) {
 309          printf("error in file %s layer=%d min. is 1 max is %d/n",
 310                  filename,l,fNlayers);
 311          continue;
 312       }// end if l
 313       id[0] = l;id[1] = a;id[2] = d;
 314       tran[0] = tran[1] = tran[2]  = 0.0;
 315       tran[0] = (Double_t)x;tran[1] = (Double_t)y;tran[2] = (Double_t)z;
 316       rot6[0] = rot6[1] = rot6[2] = rot6[3] = rot6[4] = rot6[5] =0.0;
 317       rot6[0] = (Double_t)o;rot6[1] = (Double_t)p;rot6[2] = (Double_t)q;
 318       rot6[3] = (Double_t)r;rot6[4] = (Double_t)s;rot6[5] = (Double_t)t;
 319       switch (l){
 320       case 1: case 2: // layer 1 or2 SPD
 321           fGm[lm++] = new AliITSgeomMatrix(rot6,0,id,tran);
 322           break;
 323       case 3: case 4: // layer 3 or 4 SDD
 324           fGm[lm++] = new AliITSgeomMatrix(rot6,1,id,tran);
 325           break;
 326       case 5: case 6: // layer 5 or 6 SSD
 327           fGm[lm++] = new AliITSgeomMatrix(rot6,2,id,tran);
 328           break;
 329       } // end switch
 330    } // end while ever loop
 331    fclose(pf);
 332 }
 333
 334 //________________________________________________________________________
 335 AliITSgeom::AliITSgeom(AliITSgeom &source){
 336 ////////////////////////////////////////////////////////////////////////
 337 //     The copy constructor for the AliITSgeom class. It calls the
 338 // = operator function. See the = operator function for more details.
 339 ////////////////////////////////////////////////////////////////////////
 340
 341     *this = source;  // Just use the = operator for now.
 342
 343     return;
 344 }
 345
 346 //________________________________________________________________________
 347 void AliITSgeom::operator=(AliITSgeom &source){
 348 ////////////////////////////////////////////////////////////////////////
 349 //     The = operator function for the AliITSgeom class. It makes an
 350 // independent copy of the class in such a way that any changes made
 351 // to the copied class will not affect the source class in any way.
 352 // This is required for many ITS alignment studies where the copied
 353 // class is then modified by introducing some misalignment.
 354 ////////////////////////////////////////////////////////////////////////
 355    Int_t i;
 356
 357    if(this == &source) return; // don't assign to ones self.
 358
 359    // if there is an old structure allocated delete it first.
 360    if(this->fGm != 0){
 361       for(i=0;i<this->fNmodules;i++) delete this->fGm[i];
 362       delete this->fGm;
 363    } // end if fGm != 0
 364    if(fNlad != 0) delete[] fNlad;
 365    if(fNdet != 0) delete[] fNdet;
 366
 367    this->fTrans    = source.fTrans;
 368    this->fNmodules = source.fNmodules;
 369    this->fNlayers = source.fNlayers;
 370    this->fNlad = new Int_t[fNlayers];
 371    for(i=0;i<this->fNlayers;i++) this->fNlad[i] = source.fNlad[i];
 372    this->fNdet = new Int_t[fNlayers];
 373    for(i=0;i<this->fNlayers;i++) this->fNdet[i] = source.fNdet[i];
 374    this->fShape = new TObjArray(*(source.fShape));//This does not make a proper copy.
 375    this->fGm = new AliITSgeomMatrix*[this->fNmodules];
 376    for(i=0;i<this->fNmodules;i++){
 377        this->fGm[i] = new AliITSgeomMatrix(*(source.fGm[i]));
 378    } // end for i
 379    return;
 380 }//_____________________________________________________________________
 381 Int_t AliITSgeom::GetModuleIndex(const Int_t lay,const Int_t lad,
 382                                  const Int_t det){
 383 ////////////////////////////////////////////////////////////////////////
 384 //      This routine computes the module index number from the layer,
 385 // ladder, and detector numbers. The number of ladders and detectors
 386 // per layer is determined when this geometry package is constructed,
 387 // see AliITSgeom(const char *filename) for specifics.
 388 ////////////////////////////////////////////////////////////////////////
 389     Int_t i,j,k,id[3];
 390
 391     i = fNdet[lay-1] * (lad-1) + det - 1;
 392     j = 0;
 393     for(k=0;k<lay-1;k++) j += fNdet[k]*fNlad[k];
 394     i = i+j;
 395     fGm[i]->GetIndex(id);
 396     if(id[0]==lay&&id[1]==lad&&id[2]==det) return i;
 397     // Array of modules fGm is not in expected order. Search for this index
 398     for(i=0;i<fNmodules;i++){
 399         fGm[i]->GetIndex(id);
 400         if(id[0]==lay&&id[1]==lad&&id[2]==det) return i;
 401     } // end for i
 402     // This layer ladder and detector combination does not exist return -1.
 403     return -1;
 404 }
 405 //______________________________________________________________________
 406 void AliITSgeom::GetModuleId(const Int_t index,
 407                              Int_t &lay,Int_t &lad,Int_t &det){
 408 ////////////////////////////////////////////////////////////////////////
 409 //      This routine computes the layer, ladder and detector number
 410 // given the module index number. The number of ladders and detectors
 411 // per layer is determined when this geometry package is constructed,
 412 // see AliITSgeom(const char *filename) for specifics.
 413 ////////////////////////////////////////////////////////////////////////
 414     Int_t id[3];
 415
 416     fGm[index]->GetIndex(id);
 417     lay = id[0]; lad = id[1]; det = id[2];
 418     return;
 419
 420     // The old way kept for posterity.
 421 /*
 422     Int_t i,j,k;
 423     j = 0;
 424     for(k=0;k<fNlayers;k++){
 425         j += fNdet[k]*fNlad[k];
 426         if(j>index)break;
 427     } // end for k
 428     lay = k+1;
 429     i = index -j + fNdet[k]*fNlad[k];
 430     j = 0;
 431     for(k=0;k<fNlad[lay-1];k++){
 432         j += fNdet[lay-1];
 433         if(j>i)break;
 434     } // end for k
 435     lad = k+1;
 436     det = 1+i-fNdet[lay-1]*k;
 437     return;
 438 */
 439 }
 440 //___________________________________________________________________________
 441 Int_t AliITSgeom::GetStartDet(const Int_t dtype){
 442   /////////////////////////////////////////////////////////////////////////
 443   // returns the starting module index value for a give type of detector id
 444   /////////////////////////////////////////////////////////////////////////
 445
 446   switch(dtype){
 447   case 0:
 448      return GetModuleIndex(1,1,1);
 449      break;
 450   case 1:
 451      return GetModuleIndex(3,1,1);
 452      break;
 453   case 2:
 454      return GetModuleIndex(5,1,1);
 455      break;
 456   default:
 457      printf("<AliITSgeom::GetFirstDet> undefined detector type\n");
 458      return 0;
 459   } // end switch
 460
 461   printf("<AliITSgeom::GetFirstDet> undefined detector type\n");
 462   return 0;
 463 }
 464
 465 //___________________________________________________________________________
 466 Int_t AliITSgeom::GetLastDet(const Int_t dtype){
 467   /////////////////////////////////////////////////////////////////////////
 468   // returns the last module index value for a give type of detector id
 469   /////////////////////////////////////////////////////////////////////////
 470
 471   switch(dtype){
 472   case 0:
 473      return GetLastSPD();
 474      break;
 475    case 1:
 476      return GetLastSDD();
 477      break;
 478    case 2:
 479      return GetLastSSD();
 480      break;
 481    default:
 482      printf("<AliITSgeom::GetLastDet> undefined detector type\n");
 483      return 0;
 484   } // end switch
 485
 486   printf("<AliITSgeom::GetLastDet> undefined detector type\n");
 487   return 0;
 488 }
 489
 490 //___________________________________________________________________________
 491 void AliITSgeom::PrintComparison(FILE *fp,AliITSgeom *other){
 492 ////////////////////////////////////////////////////////////////////////
 493 //     This function was primarily created for diagnostic reasons. It
 494 // print to a file pointed to by the file pointer fp the difference
 495 // between two AliITSgeom classes. The format of the file is basicly,
 496 // define d? to be the difference between the same element of the two
 497 // classes. For example dfrx = this->fGm[i][j].frx - other->fGm[i][j].frx.
 498 // if(at least one of dfx0, dfy0, dfz0,dfrx,dfry,dfrz are non zero) then print
 499 // layer ladder detector dfx0 dfy0 dfz0 dfrx dfry dfrz
 500 // if(at least one of the 9 elements of dfr[] are non zero) then print
 501 // layer ladder detector dfr[0] dfr[1] dfr[2]
 502 //                       dfr[3] dfr[4] dfr[5]
 503 //                       dfr[6] dfr[7] dfr[8]
 504 // Only non zero values are printed to save space. The differences are
 505 // typical written to a file because there are usually a lot of numbers
 506 // printed out and it is usually easier to read them in some nice editor
 507 // rather than zooming quickly past you on a screen. fprintf is used to
 508 // do the printing. The fShapeIndex difference is not printed at this time.
 509 ////////////////////////////////////////////////////////////////////////
 510    Int_t    i,j,idt[3],ido[3];
 511    Double_t tt[3],to[3];  // translation
 512    Double_t rt[3],ro[3];  // phi in radians
 513    Double_t mt[3][3],mo[3][3]; // matrixes
 514    AliITSgeomMatrix *gt,*go;
 515    Bool_t   t;
 516
 517    for(i=0;i<this->fNmodules;i++){
 518          gt  =  this->GetGeomMatrix(i);
 519          go  = other->GetGeomMatrix(i);
 520          gt->GetIndex(idt);
 521          go->GetIndex(ido);
 522          t = kFALSE;
 523          for(i=0;i<3;i++) t = t&&idt[i]!=ido[i];
 524          if(t) fprintf(fp,"%4.4d %1.1d %2.2d %2.2d %1.1d %2.2d %2.2d\n",i,
 525                        idt[0],idt[1],idt[2],ido[0],ido[1],ido[2]);
 526          gt->GetTranslation(tt);
 527          go->GetTranslation(to);
 528          gt->GetAngles(rt);
 529          go->GetAngles(ro);
 530          t = kFALSE;
 531          for(i=0;i<3;i++) t = t&&tt[i]!=to[i];
 532          if(t) fprintf(fp,"%1.1d %2.2d %2.2d dTrans=%f %f %f drot=%f %f %f\n",
 533                        idt[0],idt[1],idt[2],
 534                        tt[0]-to[0],tt[1]-to[1],tt[2]-to[2],
 535                        rt[0]-ro[0],rt[1]-ro[1],rt[2]-ro[2]);
 536          t = kFALSE;
 537          gt->GetMatrix(mt);
 538          go->GetMatrix(mo);
 539          for(i=0;i<3;i++)for(j=0;j<3;j++)  t = mt[i][j] != mo[i][j];
 540          if(t){
 541              fprintf(fp,"%1.1d %2.2d %2.2d dfr= %e %e %e\n",
 542                      idt[0],idt[1],idt[2],
 543                  mt[0][0]-mo[0][0],mt[0][1]-mo[0][1],mt[0][2]-mo[0][2]);
 544              fprintf(fp,"        dfr= %e %e %e\n",
 545                  mt[1][0]-mo[1][0],mt[1][1]-mo[1][1],mt[1][2]-mo[1][2]);
 546              fprintf(fp,"        dfr= %e %e %e\n",
 547                  mt[2][0]-mo[2][0],mt[2][1]-mo[2][1],mt[2][2]-mo[2][2]);
 548          } // end if t
 549    } // end for i
 550    return;
 551 }
 552
 553 //___________________________________________________________________________
 554 void AliITSgeom::PrintData(FILE *fp,
 555                            const Int_t lay,const Int_t lad,const Int_t det){
 556 ////////////////////////////////////////////////////////////////////////
 557 //     This function prints out the coordinate transformations for
 558 // the particular detector defined by layer, ladder, and detector
 559 // to the file pointed to by the File pointer fp. fprintf statements
 560 // are used to print out the numbers. The format is
 561 // layer ladder detector Trans= fx0 fy0 fz0 rot= frx fry frz Shape=fShapeIndex
 562 //                         dfr= fr[0] fr[1] fr[2]
 563 //                         dfr= fr[3] fr[4] fr[5]
 564 //                         dfr= fr[6] fr[7] fr[8]
 565 // By indicating which detector, some control over the information
 566 // is given to the user. The output it written to the file pointed
 567 // to by the file pointer fp. This can be set to stdout if you want.
 568 ////////////////////////////////////////////////////////////////////////
 569    AliITSgeomMatrix *gt;
 570    Double_t t[3],r[3],m[3][3];
 571
 572    gt = this->GetGeomMatrix(GetModuleIndex(lay,lad,det));
 573    gt->GetTranslation(t);
 574    gt->GetAngles(r);
 575    fprintf(fp,"%1.1d %2.2d %2.2d Trans=%f %f %f rot=%f %f %f Shape=%d\n",
 576            lay,lad,det,t[0],t[1],t[2],r[0],r[1],r[2],
 577            gt->GetDetectorIndex());
 578    gt->GetMatrix(m);
 579    fprintf(fp,"        dfr= %e %e %e\n",m[0][0],m[0][1],m[0][2]);
 580    fprintf(fp,"        dfr= %e %e %e\n",m[1][0],m[1][1],m[1][2]);
 581    fprintf(fp,"        dfr= %e %e %e\n",m[2][0],m[2][1],m[2][2]);
 582    return;
 583 }
 584 //___________________________________________________________________________
 585 ofstream & AliITSgeom::PrintGeom(ofstream &R__b){
 586 ////////////////////////////////////////////////////////////////////////
 587 //     The default Streamer function "written by ROOT" doesn't write out
 588 // the arrays referenced by pointers. Therefore, a specific Streamer function
 589 // has to be written. This function should not be modified but instead added
 590 // on to so that older versions can still be read. The proper handling of
 591 // the version dependent streamer function hasn't been written do to the lack
 592 // of finding an example at the time of writing.
 593 ////////////////////////////////////////////////////////////////////////
 594    // Stream an object of class AliITSgeom.
 595     Int_t i;
 596
 597     R__b.setf(ios::scientific);
 598     R__b << fTrans << " ";
 599     R__b << fNmodules << " ";
 600     R__b << fNlayers << " ";
 601     for(i=0;i<fNlayers;i++) R__b << fNlad[i] << " ";
 602     for(i=0;i<fNlayers;i++) R__b << fNdet[i] << "\n";
 603     for(i=0;i<fNmodules;i++) {
 604         R__b <<setprecision(16) << *(fGm[i]) << "\n";
 605     } // end for i
 606     return R__b;
 607 }
 608 //___________________________________________________________________________
 609 ifstream & AliITSgeom::ReadGeom(ifstream &R__b){
 610 ////////////////////////////////////////////////////////////////////////
 611 //     The default Streamer function "written by ROOT" doesn't write out
 612 // the arrays referenced by pointers. Therefore, a specific Streamer function
 613 // has to be written. This function should not be modified but instead added
 614 // on to so that older versions can still be read. The proper handling of
 615 // the version dependent streamer function hasn't been written do to the lack
 616 // of finding an example at the time of writing.
 617 ////////////////////////////////////////////////////////////////////////
 618    // Stream an object of class AliITSgeom.
 619       Int_t i;
 620
 621       fNlad = new Int_t[fNlayers];
 622       fNdet = new Int_t[fNlayers];
 623       if(fGm!=0){
 624           for(i=0;i<fNmodules;i++) delete fGm[i];
 625           delete fGm;
 626       } // end if fGm!=0
 627
 628       R__b >> fTrans >> fNmodules >> fNlayers;
 629       fNlad = new Int_t[fNlayers];
 630       fNdet = new Int_t[fNlayers];
 631       for(i=0;i<fNlayers;i++) R__b >> fNlad[i];
 632       for(i=0;i<fNlayers;i++) R__b >> fNdet[i];
 633       fGm = new AliITSgeomMatrix*[fNmodules];
 634       for(i=0;i<fNmodules;i++){
 635           fGm[i] = new AliITSgeomMatrix;
 636           R__b >> *(fGm[i]);
 637       } // end for i
 638       return R__b;
 639 }
 640 //___________________________________________________________________________
 641 void AliITSgeom::Streamer(TBuffer &R__b){
 642 ////////////////////////////////////////////////////////////////////////
 643 //     The default Streamer function "written by ROOT" doesn't write out
 644 // the arrays referenced by pointers. Therefore, a specific Streamer function
 645 // has to be written. This function should not be modified but instead added
 646 // on to so that older versions can still be read. The proper handling of
 647 // the version dependent streamer function hasn't been written do to the lack
 648 // of finding an example at the time of writing.
 649 ////////////////////////////////////////////////////////////////////////
 650     // Stream an object of class AliITSgeom.
 651     Int_t  i,j,k,l;
 652     UInt_t R__s=0, R__c=0;
 653
 654     if (R__b.IsReading()) {
 655         Version_t R__v = R__b.ReadVersion();
 656         if (R__v==1) {
 657             if(fNlad!=0) delete[] fNlad;
 658             if(fNdet!=0) delete[] fNdet;
 659             if(fGm!=0){
 660                 for(i=0;i<fNlayers;i++) delete[] fGm[i];
 661                 delete[] fGm;
 662             } // end if fGm!=0
 663             Int_t    idt,id[3],inmax;
 664             Double_t t[3],r[3],m[9],s[3][3];
 665             TObject::Streamer(R__b);
 666             fTrans = 0;
 667             R__b >> fNlayers;
 668             fNlad = new Int_t[fNlayers];
 669             fNdet = new Int_t[fNlayers];
 670             for(i=0;i<fNlayers;i++) R__b >> fNlad[i];
 671             for(i=0;i<fNlayers;i++) R__b >> fNdet[i];
 672             fNmodules = GetModuleIndex(fNlayers,fNlad[fNlayers-1],
 673                                        fNdet[fNlayers-1]);
 674             fGm = new AliITSgeomMatrix*[fNmodules];
 675             inmax = 0;
 676             for(i=0;i<fNlayers;i++){
 677                 for(j=0;j<fNlad[i]*fNdet[i];j++){
 678                     R__b >> idt;
 679                     R__b >> t[0];
 680                     R__b >> t[1];
 681                     R__b >> t[2];
 682                     R__b >> r[0];
 683                     R__b >> r[1];
 684                     R__b >> r[2];
 685                     for(k=0;k<9;k++) R__b >> m[k];
 686                     for(k=0;k<3;k++)for(l=0;l<3;l++) s[k][l] = m[3*k+l];
 687                     GetModuleId(inmax,id[0],id[1],id[2]);
 688                     fGm[inmax++] = new AliITSgeomMatrix(idt,id,s,t);
 689                 } // end for j
 690             } // end for i
 691             R__b >> fShape;
 692         } else if(R__v==2){
 693             if(fNlad!=0) delete[] fNlad;
 694             if(fNdet!=0) delete[] fNdet;
 695             if(fGm!=0){for(i=0;i<fNmodules;i++) delete fGm[i];delete[] fGm;}
 696             TObject::Streamer(R__b);
 697             R__b >> fTrans;
 698             R__b >> fNlayers;
 699             R__b.ReadArray(fNlad);
 700             R__b.ReadArray(fNdet);
 701             R__b >> fShape;
 702             R__b >> fNmodules;
 703             fGm = new AliITSgeomMatrix*[fNmodules];
 704             for(i=0;i<fNmodules;i++){
 705                 fGm[i] = new AliITSgeomMatrix;
 706                 fGm[i]->Streamer(R__b);
 707             } // end for i
 708             //R__b.ReadArray(fGm);
 709             R__b.CheckByteCount(R__s, R__c, AliITSgeom::IsA());
 710         } // end if R__v==?
 711     } else { // Writing.
 712         R__c = R__b.WriteVersion(AliITSgeom::IsA(), kTRUE);
 713         TObject::Streamer(R__b);
 714         R__b << fTrans;
 715         R__b << fNlayers;
 716         R__b.WriteArray(fNlad, fNlayers);
 717         R__b.WriteArray(fNdet, fNlayers);
 718         R__b << fShape;
 719         R__b << fNmodules;
 720         //R__b.WriteArray(fGm, __COUNTER__);
 721         for(i=0;i<fNmodules;i++){
 722             fGm[i]->Streamer(R__b);
 723         } // end for i
 724         R__b.SetByteCount(R__c, kTRUE);
 725     } // end if reading/writing.
 726 }
 727 //______________________________________________________________________
 728 //     The following routines modify the transformation of "this"
 729 // geometry transformations in a number of different ways.
 730 //______________________________________________________________________
 731 void AliITSgeom::GlobalChange(const Float_t *tran,const Float_t *rot){
 732 ////////////////////////////////////////////////////////////////////////
 733 //     This function performs a Cartesian translation and rotation of
 734 // the full ITS from its default position by an amount determined by
 735 // the three element arrays dtranslation and drotation. If every element
 736 // of dtranslation and drotation are zero then there is no change made
 737 // the geometry. The change is global in that the exact same translation
 738 // and rotation is done to every detector element in the exact same way.
 739 // The units of the translation are those of the Monte Carlo, usually cm,
 740 // and those of the rotation are in radians. The elements of dtranslation
 741 // are dtranslation[0] = x, dtranslation[1] = y, and dtranslation[2] = z.
 742 // The elements of drotation are drotation[0] = rx, drotation[1] = ry, and
 743 // drotation[2] = rz. A change in x will move the hole ITS in the ALICE
 744 // global x direction, the same for a change in y. A change in z will
 745 // result in a translation of the ITS as a hole up or down the beam line.
 746 // A change in the angles will result in the inclination of the ITS with
 747 // respect to the beam line, except for an effective rotation about the
 748 // beam axis which will just rotate the ITS as a hole about the beam axis.
 749 ////////////////////////////////////////////////////////////////////////
 750    Int_t    i,j;
 751    Double_t t[3],r[3];
 752    AliITSgeomMatrix *g;
 753
 754    fTrans = (fTrans && 0xfffd) + 2;  // set bit 1 true.
 755    for(i=0;i<fNmodules;i++){
 756          g = this->GetGeomMatrix(i);
 757          g->GetTranslation(t);
 758          g->GetAngles(r);
 759          for(j=0;j<3;j++){
 760               t[j] += tran[j];
 761               r[j] += rot[j];
 762          } // end for j
 763          g->SetTranslation(t);
 764          g->SetAngles(r);
 765    } // end for i
 766    return;
 767 }
 768 //___________________________________________________________________________
 769 void AliITSgeom::GlobalCylindericalChange(const Float_t *tran,const Float_t *rot){
 770 ////////////////////////////////////////////////////////////////////////
 771 //     This function performs a cylindrical translation and rotation of
 772 // each ITS element by a fixed about in radius, rphi, and z from its
 773 // default position by an amount determined by the three element arrays
 774 // dtranslation and drotation. If every element of dtranslation and
 775 // drotation are zero then there is no change made the geometry. The
 776 // change is global in that the exact same distance change in translation
 777 // and rotation is done to every detector element in the exact same way.
 778 // The units of the translation are those of the Monte Carlo, usually cm,
 779 // and those of the rotation are in radians. The elements of dtranslation
 780 // are dtranslation[0] = r, dtranslation[1] = rphi, and dtranslation[2] = z.
 781 // The elements of drotation are drotation[0] = rx, drotation[1] = ry, and
 782 // drotation[2] = rz. A change in r will results in the increase of the
 783 // radius of each layer by the same about. A change in rphi will results in
 784 // the rotation of each layer by a different angle but by the same
 785 // circumferential distance. A change in z will result in a translation
 786 // of the ITS as a hole up or down the beam line. A change in the angles
 787 // will result in the inclination of the ITS with respect to the beam
 788 // line, except for an effective rotation about the beam axis which will
 789 // just rotate the ITS as a hole about the beam axis.
 790 ////////////////////////////////////////////////////////////////////////
 791    Int_t    i,j;
 792    Double_t t[3],ro[3],r,r0,phi,rphi;
 793    AliITSgeomMatrix *g;
 794
 795    fTrans = (fTrans && 0xfffd) + 2;  // set bit 1 true.
 796    for(i=0;i<fNmodules;i++){
 797          g = this->GetGeomMatrix(i);
 798          g->GetTranslation(t);
 799          g->GetAngles(ro);
 800          r = r0= TMath::Hypot(t[1],t[0]);
 801          phi   = TMath::ATan2(t[1],t[0]);
 802          rphi  = r0*phi;
 803          r    += tran[0];
 804          rphi += tran[1];
 805          phi   = rphi/r0;
 806          t[0]  = r*TMath::Cos(phi);
 807          t[1]  = r*TMath::Sin(phi);
 808          t[2] += tran[2];
 809          for(j=0;j<3;j++){
 810               ro[j] += rot[j];
 811          } // end for j
 812          g->SetTranslation(t);
 813          g->SetAngles(ro);
 814    } // end for i
 815    return;
 816 }
 817 //___________________________________________________________________________
 818 void AliITSgeom::RandomChange(const Float_t *stran,const Float_t *srot){
 819 ////////////////////////////////////////////////////////////////////////
 820 //     This function performs a Gaussian random displacement and/or
 821 // rotation about the present global position of each active
 822 // volume/detector of the ITS. The sigma of the random displacement
 823 // is determined by the three element array stran, for the
 824 // x y and z translations, and the three element array srot,
 825 // for the three rotation about the axis x y and z.
 826 ////////////////////////////////////////////////////////////////////////
 827    Int_t    i,j;
 828    Double_t t[3],r[3];
 829    TRandom ran;
 830    AliITSgeomMatrix *g;
 831
 832    fTrans = (fTrans && 0xfffd) + 2;  // set bit 1 true.
 833    for(i=0;i<fNmodules;i++){
 834          g = this->GetGeomMatrix(i);
 835          g->GetTranslation(t);
 836          g->GetAngles(r);
 837          for(j=0;j<3;j++){
 838               t[j] += ran.Gaus(0.0,stran[j]);
 839               r[j] += ran.Gaus(0.0, srot[j]);
 840          } // end for j
 841          g->SetTranslation(t);
 842          g->SetAngles(r);
 843    } // end for i
 844    return;
 845 }
 846 //___________________________________________________________________________
 847 void AliITSgeom::RandomCylindericalChange(const Float_t *stran,
 848                                           const Float_t *srot){
 849 ////////////////////////////////////////////////////////////////////////
 850 //     This function performs a Gaussian random displacement and/or
 851 // rotation about the present global position of each active
 852 // volume/detector of the ITS. The sigma of the random displacement
 853 // is determined by the three element array stran, for the
 854 // r rphi and z translations, and the three element array srot,
 855 // for the three rotation about the axis x y and z. This random change
 856 // in detector position allow for the simulation of a random uncertainty
 857 // in the detector positions of the ITS.
 858 ////////////////////////////////////////////////////////////////////////
 859    Int_t    i,j;
 860    Double_t t[3],ro[3],r,r0,phi,rphi;
 861    TRandom ran;
 862    AliITSgeomMatrix *g;
 863
 864    fTrans = (fTrans && 0xfffd) + 2;  // set bit 1 true.
 865    for(i=0;i<fNmodules;i++){
 866          g = this->GetGeomMatrix(i);
 867          g->GetTranslation(t);
 868          g->GetAngles(ro);
 869          r = r0= TMath::Hypot(t[1],t[0]);
 870          phi   = TMath::ATan2(t[1],t[0]);
 871          rphi  = r0*phi;
 872          r    += ran.Gaus(0.0,stran[0]);
 873          rphi += ran.Gaus(0.0,stran[1]);
 874          phi   = rphi/r0;
 875          t[0]  = r*TMath::Cos(phi);
 876          t[1]  = r*TMath::Sin(phi);
 877          t[2] += ran.Gaus(0.0,stran[2]);
 878          for(j=0;j<3;j++){
 879               ro[j] += ran.Gaus(0.0, srot[j]);
 880          } // end for j
 881          g->SetTranslation(t);
 882          g->SetAngles(ro);
 883    } // end for i
 884    return;
 885 }
 886 //______________________________________________________________________
 887 void AliITSgeom::GeantToTracking(AliITSgeom &source){
 888 /////////////////////////////////////////////////////////////////////////
 889 //     Copy the geometry data but change it to go between the ALICE
 890 // Global coordinate system to that used by the ITS tracking. A slightly
 891 // different coordinate system is used when tracking. This coordinate
 892 // system is only relevant when the geometry represents the cylindrical
 893 // ALICE ITS geometry. For tracking the Z axis is left alone but X-> -Y
 894 // and Y-> X such that X always points out of the ITS cylinder for every
 895 // layer including layer 1 (where the detectors are mounted upside down).
 896 //Begin_Html
 897 /*
 898 <img src="picts/ITS/AliITSgeomMatrix_T1.gif">
 899 */
 900 //End_Html
 901 ////////////////////////////////////////////////////////////////////////
 902    Int_t    i,j,k,l,id[3];
 903    Double_t R0[3][3],R1[3][3];
 904    Double_t A0[3][3] = {{0.,+1.,0.},{-1.,0.,0.},{0.,0.,+1.}};
 905    Double_t A1[3][3] = {{0.,-1.,0.},{+1.,0.,0.},{0.,0.,+1.}};
 906
 907    *this = source;  // copy everything
 908    for(i=0;i<GetIndexMax();i++){
 909        fGm[i]->GetIndex(id);
 910        fGm[i]->GetMatrix(R0);
 911        if(id[0]==1){ // Layer 1 is treated different from the others.
 912            for(j=0;j<3;j++) for(k=0;k<3;k++){
 913                R1[j][k] = 0.;
 914                for(l=0;l<3;l++) R1[j][k] += A0[j][l]*R0[l][k];
 915            } // end for j,k
 916        }else{
 917            for(j=0;j<3;j++) for(k=0;k<3;k++){
 918                R1[j][k] = 0.;
 919                for(l=0;l<3;l++) R1[j][k] += A1[j][l]*R0[l][k];
 920            } // end for j,k
 921        } // end if
 922        fGm[i]->SetMatrix(R1);
 923    } // end for i
 924    this->fTrans = (this->fTrans && 0xfffe) + 1;  // set bit 0 true.
 925    return;
 926 }
 927 //______________________________________________________________________
 928 Int_t AliITSgeom::GetNearest(const Double_t g[3],const Int_t lay){
 929 ////////////////////////////////////////////////////////////////////////
 930 //      Finds the Detector (Module) that is nearest the point g [cm] in
 931 // ALICE Global coordinates. If layer !=0 then the search is restricted
 932 // to Detectors (Modules) in that particular layer.
 933 ////////////////////////////////////////////////////////////////////////
 934      Int_t    i,l,a,e,in=0;
 935      Double_t d,dn=1.0e10;
 936      Bool_t   t=lay!=0; // skip if lay = 0 default value check all layers.
 937
 938      for(i=0;i<fNmodules;i++){
 939           if(t){GetModuleId(i,l,a,e);if(l!=lay) continue;}
 940           if((d=fGm[i]->Distance2(g))<dn){
 941                dn = d;
 942                in = i;
 943           } // end if
 944      } // end for i
 945      return in;
 946 }
 947 //______________________________________________________________________
 948 void AliITSgeom::GetNearest27(const Double_t g[3],Int_t n[27],const Int_t lay){
 949 ////////////////////////////////////////////////////////////////////////
 950 //      Finds 27 Detectors (Modules) that are nearest the point g [cm] in
 951 // ALICE Global coordinates. If layer !=0 then the search is restricted
 952 // to Detectors (Modules) in that particular layer. The number 27 comes
 953 // from including the nearest detector and all those around it (up, down,
 954 // left, right, forwards, backwards, and the corners).
 955 ////////////////////////////////////////////////////////////////////////
 956      Int_t    i,l,a,e,in[27]={0,0,0,0,0,0,0,0,0,
 957                               0,0,0,0,0,0,0,0,0,
 958                               0,0,0,0,0,0,0,0,0,};
 959      Double_t d,dn[27]={1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,
 960                         1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,
 961                         1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,
 962                         1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,
 963                         1.0e10,1.0e10,1.0e10};
 964      Bool_t   t=(lay!=0); // skip if lay = 0 default value check all layers.
 965
 966      for(i=0;i<fNmodules;i++){
 967           if(t){GetModuleId(i,l,a,e);if(l!=lay) continue;}
 968           for(a=0;a<27;a++){
 969                d = fGm[i]->Distance2(g);
 970                if(d<dn[a]){
 971                    for(e=26;e>a;e--){dn[e] = dn[e-1];in[e] = in[e-1];}
 972                    dn[a] = d; in[a] = i;
 973                } // end if d<dn[i]
 974           } // end for a
 975      } // end for i
 976      for(i=0;i<27;i++) n[i] = in[i];
 977 }
 978 //----------------------------------------------------------------------