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.20  2002/10/14 14:57:00  hristov
  19 Merging the VirtualMC branch to the main development branch (HEAD)
  20
  21 Revision 1.18.8.1  2002/07/24 09:27:50  alibrary
  22 Updating on VirtualMC
  23
  24 Revision 1.19  2002/05/31 21:07:42  mariana
  25 Fix memory leak
  26
  27 Revision 1.18  2001/08/24 21:06:37  nilsen
  28 Added more documentation, fixed up some coding violations, and some
  29 forward declorations.
  30
  31 Revision 1.17  2001/07/27 08:06:48  hristov
  32 Use global gRandom generator (M.Ivanov)
  33
  34 Revision 1.16  2001/02/08 23:57:00  nilsen
  35 Fixed up some informational printouts.
  36
  37 Revision 1.15  2001/02/07 20:23:21  nilsen
  38 Fixed bug with HP and no unget in iostream.h. Now using putback instead.
  39 Other changes and fixes also included.
  40
  41 Revision 1.14  2001/02/03 00:00:29  nilsen
  42 New version of AliITSgeom and related files. Now uses automatic streamers,
  43 set up for new formatted .det file which includes detector information.
  44 Additional smaller modifications are still to come.
  45
  46 Revision 1.11  2000/10/02 16:32:35  barbera
  47 Forward declaration added
  48
  49 Revision 1.4.4.15  2000/10/02 15:52:05  barbera
  50 Forward declaration added
  51
  52 Revision 1.10  2000/09/05 14:25:50  nilsen
  53 Made fixes for HP compiler. All function parameter default values placed
  54 in .h file. Fixed the usual problem with HP compilers and the "for(Int_t i..."
  55 business. Replaced casting (Double_t [3][3]) to (Double_t (*)[3]) for HP.
  56 Lastly removed all "const" before function parameters which were 2 dim. arrays,
  57 because on HP root generates some strange code (?). Thanks Peter for the
  58 changes.
  59
  60 Revision 1.9  2000/08/29 20:19:03  nilsen
  61 Removed dependency on structure AliITSeomS and replaced it with class
  62 AliITSgeomMatrix. Added many new functions with many new arguments. Most
  63 in the form of in line functions for speed.
  64
  65 Revision 1.4.4.6  2000/06/04 16:33:32  Nilsen
  66 A restructured AliITSgeom class. Now used AliITSgeomMatrix.
  67
  68 Revision 1.4.4.5  2000/03/04 23:42:39  Nilsen
  69 Updated the comments/documentations and improved the maintainability of the
  70 code.
  71
  72 Revision 1.4.4.4  2000/03/02 21:27:07  Nilsen
  73 Added two functions, SetByAngles and SetTrans.
  74
  75 Revision 1.4.4.3  2000/01/23 03:09:10  Nilsen
  76 // fixed compiler warnings for new function LtLErrorMatrix(...)
  77
  78 Revision 1.4.4.2  2000/01/19 23:18:20  Nilsen
  79 Added transformations of Error matrix to AliITSgeom and fixed some typos
  80 in AliITS.h and AliITShitIndex.h
  81
  82 Revision 1.4.4.1  2000/01/12 19:03:32  Nilsen
  83 This is the version of the files after the merging done in December 1999.
  84 See the ReadMe110100.txt file for details
  85
  86 Revision 1.4  1999/10/15 07:03:20  fca
  87 Fixed bug in GetModuleId(Int_t index,Int_t &lay,Int_t &lad, Int_t &det) and
  88 a typo in the creator. aliroot need to be rerun to get a fixed geometry.
  89
  90 Revision 1.3  1999/10/04 15:20:12  fca
  91 Correct syntax accepted by g++ but not standard for static members, remove minor warnings
  92
  93 Revision 1.2  1999/09/29 09:24:20  fca
  94 Introduction of the Copyright and cvs Log
  95
  96 */
  97
  98 ///////////////////////////////////////////////////////////////////////
  99 // ITS geometry manipulation routines.                               //
 100 // Created April 15 1999.                                            //
 101 // version: 0.0.0                                                    //
 102 // By: Bjorn S. Nilsen                                               //
 103 // version: 0.0.1                                                    //
 104 // Updated May 27 1999.                                              //
 105 // Added Cylindrical random and global based changes.               //
 106 // Added  function PrintComparison.                                  //
 107 ///////////////////////////////////////////////////////////////////////
 108
 109
 110 ////////////////////////////////////////////////////////////////////////
 111 //     The local coordinate system by, default, is show in the following
 112 // figures. Also shown are the ladder numbering scheme.
 113 //Begin_Html
 114 /*
 115 <img src="picts/ITS/AliITSgeomMatrix_L1.gif">
 116 </pre>
 117 <br clear=left>
 118 <font size=+2 color=blue>
 119 <p>This shows the relative geometry differences between the ALICE Global
 120 coordinate system and the local detector coordinate system.
 121 </font>
 122 <pre>
 123
 124 <pre>
 125 <img src="picts/ITS/its1+2_convention_front_5.gif">
 126 </pre>
 127 <br clear=left>
 128 <font size=+2 color=blue>
 129 <p>This shows the front view of the SPDs and the orientation of the local
 130 pixel coordinate system. Note that the inner pixel layer has its y coordinate
 131 in the opposite direction from all of the other layers.
 132 </font>
 133 <pre>
 134
 135 <pre>
 136 <img src="picts/ITS/its3+4_convention_front_5.gif">
 137 </pre>
 138 <br clear=left>
 139 <font size=+2 color=blue>
 140 <p>This shows the front view of the SDDs and the orientation of the local
 141 pixel coordinate system.
 142 </font>
 143 <pre>
 144
 145 <pre>
 146 <img src="picts/ITS/its5+6_convention_front_5.gif">
 147 </pre>
 148 <br clear=left>
 149 <font size=+2 color=blue>
 150 <p>This shows the front view of the SSDs and the orientation of the local
 151 pixel coordinate system.
 152 </font>
 153 <pre>
 154 */
 155 //End_Html
 156 //
 157 ////////////////////////////////////////////////////////////////////////
 158
 159 ////////////////////////////////////////////////////////////////////////
 160 //
 161 // version: 0
 162 // Written by Bjorn S. Nilsen
 163 //
 164 // Data Members:
 165 //
 166 // Int_t fNlayers
 167 //     The number of ITS layers for this geometry. By default this
 168 //  is 6, but can be modified by the creator function if there are
 169 // more layers defined.
 170 //
 171 // Int_t *fNlad
 172 //     A pointer to an array fNlayers long containing the number of
 173 // ladders for each layer. This array is typically created and filled
 174 // by the AliITSgeom creator function.
 175 //
 176 // Int_t *fNdet
 177 //     A pointer to an array fNlayers long containing the number of
 178 // active detector volumes for each ladder. This array is typically
 179 // created and filled by the AliITSgeom creator function.
 180 //
 181 // AliITSgeomMatrix *fGm
 182 //     A pointer to an array of AliITSgeomMatrix classes. One element
 183 // per module (detector) in the ITS. AliITSgeomMatrix basicly contains
 184 // all of the necessary information about the detector and it's coordinate
 185 // transformations.
 186 //
 187 // TObjArray *fShape
 188 //     A pointer to an array of TObjects containing the detailed shape
 189 // information for each type of detector used in the ITS. For example
 190 // I have created AliITSgeomSPD, AliITSgeomSDD, and AliITSgeomSSD as
 191 // example structures, derived from TObjects, to hold the detector
 192 // information. I would recommend that one element in each of these
 193 // structures, that which describes the shape of the active volume,
 194 // be one of the ROOT classes derived from TShape. In this way it would
 195 // be easy to have the display program display the correct active
 196 // ITS volumes. See the example classes AliITSgeomSPD, AliITSgeomSDD,
 197 // and AliITSgeomSSD for a more detailed example.
 198 ////////////////////////////////////////////////////////////////////////
 199 #include <Riostream.h>
 200 #include <stdlib.h>
 201 #include <stdio.h>
 202 #include <string.h>
 203 #include <ctype.h>
 204
 205 #include <TSystem.h>
 206 #include <TRandom.h>
 207
 208 #include "AliITSgeom.h"
 209 #include "AliITSgeomSPD.h"
 210 #include "AliITSgeomSDD.h"
 211 #include "AliITSgeomSSD.h"
 212
 213 ClassImp(AliITSgeom)
 214
 215 //______________________________________________________________________
 216 AliITSgeom::AliITSgeom(){
 217     //     The default constructor for the AliITSgeom class. It, by default,
 218     // sets fNlayers to zero and zeros all pointers.
 219     // Do not allocate anything zero everything.
 220
 221     fTrans   = 0; // standard GEANT global/local coordinate system.
 222     fNlayers = 0;
 223     fNlad    = 0;
 224     fNdet    = 0;
 225     fGm      = 0;
 226     fShape   = 0;
 227     strcpy(fVersion,"test");
 228     return;
 229 }
 230 //______________________________________________________________________
 231 AliITSgeom::AliITSgeom(Int_t itype,Int_t nlayers,Int_t *nlads,Int_t *ndets,
 232                        Int_t mods){
 233     //     A simple constructor to set basic geometry class variables
 234     // Inputs:
 235     // Int_t itype   the type of transformation kept.
 236     //               bit 0 => Standard GEANT
 237     //               bit 1 => ITS tracking
 238     //               bit 2 => A change in the coordinate system has been made.
 239     //               others are still to be defined as needed.
 240     // Int_t nlayers The number of ITS layers also set the size of the arrays
 241     // Int_t *nlads  an array of the number of ladders for each layer. This
 242     //               array must be nlayers long.
 243     // Int_t *ndets  an array of the number of detectors per ladder for each
 244     //               layer. This array must be nlayers long.
 245     // Int_t mods    The number of modules. Typicaly the sum of all the
 246     //               detectors on every layer and ladder.
 247     // Outputs:
 248     // none
 249     Int_t i;
 250
 251     fTrans    = itype;
 252     fNlayers  = nlayers;
 253     fNlad     = new Int_t[nlayers];
 254     fNdet     = new Int_t[nlayers];
 255     for(i=0;i<nlayers;i++){fNlad[i] = nlads[i];fNdet[i] = ndets[i];}
 256     fNmodules = mods;
 257     fGm       = new TObjArray(mods,0);
 258     fShape    = new TObjArray(5); // default value
 259     for(i=0;i<5;i++) fShape->AddAt(0,i);
 260     strcpy(fVersion,"test");
 261     return;
 262 }
 263 //______________________________________________________________________
 264 void AliITSgeom::CreatMatrix(Int_t mod,Int_t lay,Int_t lad,Int_t det,
 265                              AliITSDetector idet,const Double_t tran[3],
 266                              const Double_t rot[10]){
 267     // Given the translation vector tran[3] and the rotation matrix rot[1],
 268     // this function creates and adds to the TObject Array fGm the
 269     // AliITSgeomMatrix object.
 270     // Inputs are:
 271     // Int_t           mod     The module number. The location in TObjArray
 272     // Int_t           lay     The layer where this module is
 273     // Int_t           lad     On which ladder this module is
 274     // Int_t           det     Which detector on this ladder this module is
 275     // AliITSDetector idet     The type of detector see AliITSgeom.h
 276     // Double_t       tran[3]  The translation vector
 277     // Double_t       rot[10]  The rotation matrix.
 278     // Outputs are:
 279     //   none
 280     // The rot[10] matrix is set up like:
 281     /*   / rot[0]  rot[1]  rot[2] \
 282     //  |  rot[3]  rot[4]  rot[5]  |
 283     //   \ rot[6]  rot[7]  rot[8] /  if(rot[9]!=0) then the Identity matrix
 284     // is used regardless of the values in rot[0]-rot[8].
 285     */
 286     Int_t id[3];
 287     Double_t r[3][3] = {{1.0,0.0,0.0},{0.0,1.0,0.0},{0.0,0.0,1.0}};
 288
 289     if(fGm->At(mod)!=0) delete fGm->At(mod);
 290     id[0] = lay; id[1] = lad; id[2] = det;
 291     if(rot[9]!=0.0) { // null rotation
 292         r[0][0] = rot[0]; r[0][1] = rot[1]; r[0][2] = rot[2];
 293         r[1][0] = rot[3]; r[1][1] = rot[4]; r[1][2] = rot[5];
 294         r[2][0] = rot[6]; r[2][1] = rot[7]; r[2][2] = rot[8];
 295     } // end if
 296     fGm->AddAt(new AliITSgeomMatrix(idet,id,r,tran),mod);
 297 }
 298 //______________________________________________________________________
 299 AliITSgeom::~AliITSgeom(){
 300     //     The destructor for the AliITSgeom class. If the arrays fNlad,
 301     // fNdet, or fGm have had memory allocated to them, there pointer values
 302     // are non zero, then this memory space is freed and they are set
 303     // to zero. In addition, fNlayers is set to zero. The destruction of
 304     // TObjArray fShape is, by default, handled by the TObjArray destructor.
 305
 306    if(fGm!=0){
 307      //for(Int_t i=0;i<fNlayers;i++) delete fGm->At(i);
 308       fGm->Delete();
 309       delete fGm;
 310    } // end if fGm!=0
 311    if(fNlad!=0) delete[] fNlad;
 312    if(fNdet!=0) delete[] fNdet;
 313    fNlayers = 0;
 314    fNlad    = 0;
 315    fNdet    = 0;
 316    fGm      = 0;
 317    return;
 318 }
 319 //______________________________________________________________________
 320 void AliITSgeom::ReadNewFile(const char *filename){
 321     // It is generaly preferred to define the geometry in AliITSgeom
 322     // directly from the GEANT geometry, see AliITSvPPRasymm.cxx for
 323     // and example. Under some circumstances this may not be possible.
 324     // This function will read in a formatted file for all of the
 325     // information needed to define the geometry in AliITSgeom.
 326     // Unlike the older file format, this file may contain comments
 327     // and the order of the data does not need to be completely
 328     // respected. A file can be created using the function WriteNewFile
 329     // defined below.
 330     // Inputs are:
 331     // const char *filename The file name of the file to be read in.
 332     // Outputs are:
 333     //  none
 334     Int_t ncmd=9;
 335     const char *cmda[]={"Version"        ,"fTrans"  ,"fNmodules",
 336                         "fNlayers"       ,"fNladers","fNdetectors",
 337                         "fNDetectorTypes","fShape"  ,"Matrix"};
 338     Int_t i,j,lNdetTypes,ldet;
 339     char cmd[20],c;
 340     AliITSgeomSPD *spd;
 341     AliITSgeomSDD *sdd;
 342     AliITSgeomSSD *ssd;
 343     AliITSgeomMatrix *m;
 344     ifstream *fp;
 345     char *filtmp;
 346
 347     filtmp = gSystem->ExpandPathName(filename);
 348     cout << "AliITSgeom, Reading New .det file " << filtmp << endl;
 349     fp = new ifstream(filtmp,ios::in);  // open file to write
 350     while(fp->get(c)!=NULL){ // for ever loop
 351         if(c==' ') continue; // remove blanks
 352         if(c=='\n') continue;
 353         if(c=='#' || c=='!'){for(;fp->get(c)!=NULL,c!='\n';); continue;}
 354         if(c=='/'){
 355             fp->get(c);{
 356                 if(c=='/'){for(;fp->get(c)!=NULL,c!='\n';);continue;}
 357                 if(c=='*'){
 358                 NotYet:
 359                     for(;fp->get(c)!=NULL,c!='*';);
 360                     fp->get(c);{
 361                         if(c=='/') continue;
 362                         goto NotYet;
 363                     } //
 364                 } // end if c=='*'
 365             } // end if second /
 366         } // end if first /
 367         fp->putback(c);
 368         *fp >> cmd;
 369         for(i=0;i<ncmd;i++) if(strcmp(cmd,cmda[i])==0) break;
 370         switch (i){
 371         case 0:   // Version
 372             *fp >> fVersion;
 373             break;
 374         case 1:  // fTrans
 375             *fp >> fTrans;
 376             break;
 377         case 2:  // fNModules
 378             *fp >> fNmodules;
 379             if(fGm!=0){
 380                 for(j=0;j<fGm->GetEntriesFast();j++) delete fGm->At(j);
 381                 delete fGm;
 382             } // end if
 383             fGm = new TObjArray(fNmodules,0);
 384             break;
 385         case 3:  // fNlayers
 386             *fp >> fNlayers;
 387             if(fNlad!=0) delete fNlad;
 388             if(fNdet!=0) delete fNdet;
 389             fNlad = new Int_t[fNlayers];
 390             fNdet = new Int_t[fNlayers];
 391             break;
 392         case 4:  // fNladers
 393             for(j=0;j<fNlayers;j++) *fp >> fNlad[j];
 394             break;
 395         case 5:  // fNdetectors
 396             for(j=0;j<fNlayers;j++) *fp >> fNdet[j];
 397             break;
 398         case 6:  // fNDetectorTypes
 399             *fp >> lNdetTypes;
 400             if(fShape!=0){
 401                 for(j=0;j<fShape->GetEntriesFast();j++) delete fShape->At(j);
 402                 delete fShape;
 403             } // end if
 404             fShape = new TObjArray(lNdetTypes,0);
 405             break;
 406         case 7:  // fShape
 407             *fp >> ldet;
 408             if(fShape==0) fShape = new TObjArray(5,0);
 409             switch (ldet){
 410             case kSPD :
 411                 ReSetShape(ldet,(TObject*) new AliITSgeomSPD());
 412                 spd = (AliITSgeomSPD*) (fShape->At(ldet));
 413                 *fp >> *spd;
 414                 spd = 0;
 415                 break;
 416             case kSDD :
 417                 ReSetShape(ldet,(TObject*) new AliITSgeomSDD());
 418                 sdd = (AliITSgeomSDD*) (fShape->At(ldet));
 419                 *fp >> *sdd;
 420                 sdd = 0;
 421                 break;
 422             case kSSD : case kSSDp :
 423                 ReSetShape(ldet,(TObject*) new AliITSgeomSSD());
 424                 ssd = (AliITSgeomSSD*) (fShape->At(ldet));
 425                 *fp >> *ssd;
 426                 ssd = 0;
 427                 break;
 428             default:
 429                 Error("ReadNewFile","Unknown fShape type number=%d c=%c",ldet,c);
 430                 for(;fp->get(c)==NULL,c!='\n';); // skip to end of line.
 431                 break;
 432             } // end switch
 433             break;
 434         case 8:  // Matrix
 435             *fp >> ldet;
 436             if(fGm==0) fGm = new TObjArray(2270,0);
 437             if(fGm->At(ldet)!=0) delete (fGm->At(ldet));
 438             fGm->AddAt((TObject*)new AliITSgeomMatrix(),ldet);
 439             m = (AliITSgeomMatrix*) fGm->At(ldet);
 440             *fp >> *m;
 441             m = 0;
 442             break;
 443         default:
 444             Error("ReadNewFile","Data line i=%d c=%c",i,c);
 445             for(;fp->get(c)==NULL,c!='\n';); // skip this line
 446             break;
 447         } // end switch i
 448     } // end while
 449     delete fp;
 450
 451     return;
 452 }
 453 //______________________________________________________________________
 454 void AliITSgeom::WriteNewFile(const char *filename){
 455     // Writes AliITSgeom, AliITSgeomMatrix, and the defined AliITSgeomS*D
 456     // classes to a file in a format that is more readable and commendable.
 457     // Inputs are:
 458     // const char *filename The file name of the file to be write to.
 459     // Outputs are:
 460     //  none
 461     ofstream *fp;
 462     Int_t i;
 463     char *filtmp;
 464
 465     filtmp = gSystem->ExpandPathName(filename);
 466     cout << "AliITSgeom, Writing New .det file " << filtmp << endl;
 467     fp = new ofstream(filtmp,ios::out);  // open file to write
 468     *fp << "//Comment lines begin with two //, one #, or one !" << endl;
 469     *fp << "#Blank lines are skipped including /* and */ sections." << endl;
 470     *fp << "!and, in principle the order of the lines is not important" <<endl;
 471     *fp << "/* In AliITSgeom.h are defined an enumerated type called" << endl;
 472     *fp << " AliITSDetectors These are kSPD=" << (Int_t) kSPD ;
 473     *fp << ", kSDD=" << (Int_t) kSDD << ", kSSD=" << (Int_t) kSSD;
 474     *fp << ", kSSDp=" << (Int_t) kSSDp << ", and kSDDp=" << (Int_t) kSDDp;
 475     *fp << "*/" << endl;
 476     *fp << "Version " << fVersion << endl;//This should be consistent with the
 477                                            // geometry version.
 478     *fp << "fTrans " << fTrans << endl;
 479     *fp << "fNmodules " << fNmodules << endl;
 480     *fp << "fNlayers " << fNlayers << endl;
 481     *fp << "fNladers ";
 482     for(i=0;i<fNlayers;i++) *fp << fNlad[i] << " ";
 483     *fp << endl;
 484     *fp << "fNdetectors ";
 485     for(i=0;i<fNlayers;i++) *fp << fNdet[i] << " ";
 486     *fp << endl;
 487     *fp << "fNDetectorTypes " << fShape->GetEntriesFast() << endl;
 488     for(i=0;i<fShape->GetEntriesFast();i++){
 489         if(!IsShapeDefined(i)) continue; // only print out used shapes.
 490         switch (i){
 491         case kSPD :
 492             *fp << "fShape " << (Int_t) kSPD << " ";
 493             *fp << *((AliITSgeomSPD*)(fShape->At(i)));
 494             break;
 495         case kSDD :
 496             *fp << "fShape " << (Int_t) kSDD << " ";
 497             *fp << *((AliITSgeomSDD*)(fShape->At(i)));
 498             break;
 499         case kSSD : case kSSDp :
 500             *fp << "fShape " << i << " ";
 501             *fp << *((AliITSgeomSSD*)(fShape->At(i)));
 502             break;
 503         default:
 504             Error("AliITSgeom::WriteNewFile","Unknown Shape value");
 505         } // end switch (i)
 506     } // end for i
 507     for(i=0;i<fNmodules;i++){
 508         *fp << "Matrix " << i << " ";
 509         *fp << *GetGeomMatrix(i);
 510     } // end for i
 511     *fp << "//End of File" << endl;;
 512
 513     delete fp;
 514     return;
 515 }
 516 //______________________________________________________________________
 517 AliITSgeom::AliITSgeom(const char *filename){
 518     //     The constructor for the AliITSgeom class. All of the data to fill
 519     // this structure is read in from the file given my the input filename.
 520     // Inputs are:
 521     // const char *filename The file name of the file to be read in.
 522     // Outputs are:
 523     //  none
 524     FILE     *pf=0;
 525     Int_t    i,lm=0,id[3];
 526     Int_t    l,a,d;
 527     Float_t  x,y,z,o,p,q,r,s,t;
 528     Double_t rot6[6],tran[3];
 529     char     buf[200],*buff=0; // input character buffer;
 530     char *filtmp;
 531
 532     filtmp = gSystem->ExpandPathName(filename);
 533     cout << "AliITSgeom reading old .det file " << filtmp << endl;
 534     fShape = 0;
 535     strcpy(fVersion,"DefauleV5");
 536     pf = fopen(filtmp,"r");
 537
 538     fNlayers = 6; // set default number of ladders
 539  TryAgain:
 540     fNlad    = new Int_t[fNlayers];
 541     fNdet    = new Int_t[fNlayers];
 542     fNmodules = 0;
 543     // find the number of ladders and detectors in this geometry.
 544     for(i=0;i<fNlayers;i++){fNlad[i]=fNdet[i]=0;} // zero out arrays
 545     while(fgets(buf,200,pf)!=NULL){ // for ever loop
 546         for(i=0;i<200;i++)if(buf[i]!=' '){ // remove blank spaces.
 547             buff = &(buf[i]);
 548             break;
 549         } // end for i
 550         // remove blank lines and comments.
 551         if(buff[0]=='\n'||buff[0]=='#'||buff[0]=='!'||
 552            (buff[0]=='/'&&buff[1]=='/')) continue;
 553         if(isalpha(buff[0])) { // must be the new file formated file.
 554             fclose(pf);
 555             delete[] fNlad;delete[] fNdet;
 556             ReadNewFile(filename);
 557             return;
 558         } // end if isalpha(buff[0])
 559         sscanf(buff,"%d %d %d %f %f %f %f %f %f %f %f %f",
 560                &l,&a,&d,&x,&y,&z,&o,&p,&q,&r,&s,&t);
 561         if(l>lm) lm = l;
 562         if(l<1 || l>fNlayers) {
 563             printf("error in file %s layer=%d min. is 1 max is %d\n",
 564                    filename,l,fNlayers);
 565             continue;
 566         }// end if l
 567         fNmodules++;
 568         if(l<=fNlayers&&fNlad[l-1]<a) fNlad[l-1] = a;
 569         if(l<=fNlayers&&fNdet[l-1]<d) fNdet[l-1] = d;
 570     } // end while ever loop
 571     if(lm>fNlayers){
 572         delete[] fNlad;
 573         delete[] fNdet;
 574         fNlayers = lm;
 575         goto TryAgain;
 576     } // end if lm>fNlayers
 577     // counted the number of ladders and detectors now allocate space.
 578     fGm = new TObjArray(fNmodules,0);
 579
 580     // Set up Shapes for a default configuration of 6 layers.
 581     fTrans   = 0; // standard GEANT global/local coordinate system.
 582     // prepare to read in transforms
 583     lm = 0; // reuse lm as counter of modules.
 584     rewind(pf); // start over reading file
 585     while(fgets(buf,200,pf)!=NULL){ // for ever loop
 586         for(i=0;i<200;i++)if(buf[i]!=' '){ // remove blank spaces.
 587             buff = &(buf[i]);
 588             break;
 589         } // end for i
 590         // remove blank lines and comments.
 591         if(buff[0]=='\n'||buff[0]=='#'||buff[0]=='!'||
 592            (buff[0]=='/'&&buff[1]=='/')) continue;
 593         x = y = z = o = p = q = r = s = t = 0.0;
 594         sscanf(buff,"%d %d %d %f %f %f %f %f %f %f %f %f",
 595                &l,&a,&d,&x,&y,&z,&o,&p,&q,&r,&s,&t);
 596         if(l<1 || l>fNlayers) {
 597             printf("error in file %s layer=%d min. is 1 max is %d/n",
 598                    filename,l,fNlayers);
 599             continue;
 600         }// end if l
 601         id[0] = l;id[1] = a;id[2] = d;
 602         tran[0] = tran[1] = tran[2]  = 0.0;
 603         tran[0] = (Double_t)x;tran[1] = (Double_t)y;tran[2] = (Double_t)z;
 604         rot6[0] = rot6[1] = rot6[2] = rot6[3] = rot6[4] = rot6[5] =0.0;
 605         rot6[0] = (Double_t)o;rot6[1] = (Double_t)p;rot6[2] = (Double_t)q;
 606         rot6[3] = (Double_t)r;rot6[4] = (Double_t)s;rot6[5] = (Double_t)t;
 607         switch (l){
 608         case 1: case 2: // layer 1 or2 SPD
 609             fGm->AddAt(new AliITSgeomMatrix(rot6,kSPD,id,tran),lm++);
 610             break;
 611         case 3: case 4: // layer 3 or 4 SDD
 612             fGm->AddAt(new AliITSgeomMatrix(rot6,kSDD,id,tran),lm++);
 613             break;
 614         case 5: case 6: // layer 5 or 6 SSD
 615             fGm->AddAt(new AliITSgeomMatrix(rot6,kSSD,id,tran),lm++);
 616             break;
 617         } // end switch
 618     } // end while ever loop
 619     fclose(pf);
 620 }
 621 //______________________________________________________________________
 622 AliITSgeom::AliITSgeom(AliITSgeom &source){
 623     //     The copy constructor for the AliITSgeom class. It calls the
 624     // = operator function. See the = operator function for more details.
 625     // Inputs are:
 626     // AliITSgeom &source  The AliITSgeom class with which to make this
 627     //                     a copy of.
 628     // Outputs are:
 629     // none.
 630
 631     *this = source;  // Just use the = operator for now.
 632     return;
 633 }
 634 //______________________________________________________________________
 635 AliITSgeom& AliITSgeom::operator=(AliITSgeom &source){
 636     //     The = operator function for the AliITSgeom class. It makes an
 637     // independent copy of the class in such a way that any changes made
 638     // to the copied class will not affect the source class in any way.
 639     // This is required for many ITS alignment studies where the copied
 640     // class is then modified by introducing some misalignment.
 641     // Inputs are:
 642     // AliITSgeom &source  The AliITSgeom class with which to make this
 643     //                     a copy of.
 644     // Outputs are:
 645     // return  *this       The a new copy of source.
 646    Int_t i;
 647
 648    if(this == &source) return *this; // don't assign to ones self.
 649
 650    // if there is an old structure allocated delete it first.
 651    if(this->fGm != 0){
 652       for(i=0;i<this->fNmodules;i++) delete this->fGm->At(i);
 653       delete this->fGm;
 654    } // end if fGm != 0
 655    if(fNlad != 0) delete[] fNlad;
 656    if(fNdet != 0) delete[] fNdet;
 657
 658    this->fTrans    = source.fTrans;
 659    this->fNmodules = source.fNmodules;
 660    this->fNlayers = source.fNlayers;
 661    this->fNlad = new Int_t[fNlayers];
 662    for(i=0;i<this->fNlayers;i++) this->fNlad[i] = source.fNlad[i];
 663    this->fNdet = new Int_t[fNlayers];
 664    for(i=0;i<this->fNlayers;i++) this->fNdet[i] = source.fNdet[i];
 665    this->fShape = new TObjArray(*(source.fShape));//This does not make a proper copy.
 666    this->fGm = new TObjArray(this->fNmodules,0);
 667    for(i=0;i<this->fNmodules;i++){
 668        this->fGm->AddAt(new AliITSgeomMatrix(*(
 669            (AliITSgeomMatrix*)(source.fGm->At(i)))),i);
 670    } // end for i
 671    return *this;
 672 }
 673 //______________________________________________________________________
 674 Int_t AliITSgeom::GetModuleIndex(Int_t lay,Int_t lad,Int_t det){
 675     //      This routine computes the module index number from the layer,
 676     // ladder, and detector numbers. The number of ladders and detectors
 677     // per layer is determined when this geometry package is constructed,
 678     // see AliITSgeom(const char *filename) for specifics.
 679     // Inputs are:
 680     // Int_t lay  The layer number. Starting from 1.
 681     // Int_t lad  The ladder number. Starting from 1.
 682     // Int_t det  The detector number. Starting from 1.
 683     // Outputs are:
 684     // return the module index number, starting from zero.
 685     Int_t i,j,k,id[3];
 686
 687     i = fNdet[lay-1] * (lad-1) + det - 1;
 688     j = 0;
 689     for(k=0;k<lay-1;k++) j += fNdet[k]*fNlad[k];
 690     i = i+j;
 691     GetGeomMatrix(i)->GetIndex(id);
 692     if(id[0]==lay&&id[1]==lad&&id[2]==det) return i;
 693     // Array of modules fGm is not in expected order. Search for this index
 694     for(i=0;i<fNmodules;i++){
 695         GetGeomMatrix(i)->GetIndex(id);
 696         if(id[0]==lay&&id[1]==lad&&id[2]==det) return i;
 697     } // end for i
 698     // This layer ladder and detector combination does not exist return -1.
 699     return -1;
 700 }
 701 //______________________________________________________________________
 702 void AliITSgeom::GetModuleId(Int_t index,Int_t &lay,Int_t &lad,Int_t &det){
 703     //      This routine computes the layer, ladder and detector number
 704     // given the module index number. The number of ladders and detectors
 705     // per layer is determined when this geometry package is constructed,
 706     // see AliITSgeom(const char *filename) for specifics.
 707     // Inputs are:
 708     // Int_t index  The module index number, starting from zero.
 709     // Outputs are:
 710     // Int_t lay    The layer number. Starting from 1.
 711     // Int_t lad    The ladder number. Starting from 1.
 712     // Int_t det    The detector number. Starting from 1.
 713     Int_t id[3];
 714
 715     GetGeomMatrix(index)->GetIndex(id);
 716     lay = id[0]; lad = id[1]; det = id[2];
 717     return;
 718
 719     // The old way kept for posterity.
 720 /*
 721     Int_t i,j,k;
 722     j = 0;
 723     for(k=0;k<fNlayers;k++){
 724         j += fNdet[k]*fNlad[k];
 725         if(j>index)break;
 726     } // end for k
 727     lay = k+1;
 728     i = index -j + fNdet[k]*fNlad[k];
 729     j = 0;
 730     for(k=0;k<fNlad[lay-1];k++){
 731         j += fNdet[lay-1];
 732         if(j>i)break;
 733     } // end for k
 734     lad = k+1;
 735     det = 1+i-fNdet[lay-1]*k;
 736     return;
 737 */
 738 }
 739 //______________________________________________________________________
 740 Int_t AliITSgeom::GetStartDet(Int_t dtype){
 741     // returns the starting module index value for a give type of detector id.
 742     // This assumes that the detector types are different on different layers
 743     // and that they are not mixed up.
 744     // Inputs are:
 745     // Int_t dtype A detector type number. 0 for SPD, 1 for SDD, and 2 for SSD.
 746     // outputs:
 747     // return the module index for the first occurance of that detector type.
 748
 749     switch(dtype){
 750     case 0:
 751         return GetModuleIndex(1,1,1);
 752         break;
 753     case 1:
 754         return GetModuleIndex(3,1,1);
 755         break;
 756     case 2:
 757         return GetModuleIndex(5,1,1);
 758         break;
 759     default:
 760         Warning("GetStartDet","undefined detector type %d",dtype);
 761         return 0;
 762     } // end switch
 763
 764     Warning("GetStartDet","undefined detector type %d",dtype);
 765     return 0;
 766 }
 767 //______________________________________________________________________
 768 Int_t AliITSgeom::GetLastDet(Int_t dtype){
 769     // returns the last module index value for a give type of detector id.
 770     // This assumes that the detector types are different on different layers
 771     // and that they are not mixed up.
 772     // Inputs are:
 773     // Int_t dtype A detector type number. 0 for SPD, 1 for SDD, and 2 for SSD.
 774     // outputs are:
 775     // return the module index for the last occurance of that detector type.
 776
 777     switch(dtype){
 778     case 0:
 779         return GetLastSPD();
 780         break;
 781     case 1:
 782         return GetLastSDD();
 783         break;
 784     case 2:
 785         return GetLastSSD();
 786         break;
 787     default:
 788         Warning("GetLastDet","undefined detector type %d",dtype);
 789         return 0;
 790     } // end switch
 791
 792     Warning("GetLastDet","undefined detector type %d",dtype);
 793     return 0;
 794 }
 795 //______________________________________________________________________
 796 void AliITSgeom::PrintComparison(FILE *fp,AliITSgeom *other){
 797     //     This function was primarily created for diagnostic reasons. It
 798     // print to a file pointed to by the file pointer fp the difference
 799     // between two AliITSgeom classes. The format of the file is basicly,
 800     // define d? to be the difference between the same element of the two
 801     // classes. For example dfrx = this->GetGeomMatrix(i)->frx
 802     // - other->GetGeomMatrix(i)->frx.
 803     // if(at least one of dfx0, dfy0, dfz0,dfrx,dfry,dfrz are non zero) then
 804     // print layer ladder detector dfx0 dfy0 dfz0 dfrx dfry dfrz
 805     // if(at least one of the 9 elements of dfr[] are non zero) then print
 806     // layer ladder detector dfr[0] dfr[1] dfr[2]
 807     //                       dfr[3] dfr[4] dfr[5]
 808     //                       dfr[6] dfr[7] dfr[8]
 809     // Only non zero values are printed to save space. The differences are
 810     // typical written to a file because there are usually a lot of numbers
 811     // printed out and it is usually easier to read them in some nice editor
 812     // rather than zooming quickly past you on a screen. fprintf is used to
 813     // do the printing. The fShapeIndex difference is not printed at this time.
 814     // Inputs are:
 815     // FILE *fp           A file pointer to an opened file for writing in which
 816     //                    the results of the comparison will be written.
 817     // AliITSgeom *other  The other AliITSgeom class to which this one is
 818     //                    being compared.
 819     // outputs are:
 820     // none
 821     Int_t    i,j,idt[3],ido[3];
 822     Double_t tt[3],to[3];  // translation
 823     Double_t rt[3],ro[3];  // phi in radians
 824     Double_t mt[3][3],mo[3][3]; // matrixes
 825     AliITSgeomMatrix *gt,*go;
 826     Bool_t   t;
 827
 828     for(i=0;i<this->fNmodules;i++){
 829         gt  =  this->GetGeomMatrix(i);
 830         go  = other->GetGeomMatrix(i);
 831         gt->GetIndex(idt);
 832         go->GetIndex(ido);
 833         t = kFALSE;
 834         for(i=0;i<3;i++) t = t&&idt[i]!=ido[i];
 835         if(t) fprintf(fp,"%4.4d %1.1d %2.2d %2.2d %1.1d %2.2d %2.2d\n",i,
 836                       idt[0],idt[1],idt[2],ido[0],ido[1],ido[2]);
 837         gt->GetTranslation(tt);
 838         go->GetTranslation(to);
 839         gt->GetAngles(rt);
 840         go->GetAngles(ro);
 841         t = kFALSE;
 842         for(i=0;i<3;i++) t = t&&tt[i]!=to[i];
 843         if(t) fprintf(fp,"%1.1d %2.2d %2.2d dTrans=%f %f %f drot=%f %f %f\n",
 844                       idt[0],idt[1],idt[2],
 845                       tt[0]-to[0],tt[1]-to[1],tt[2]-to[2],
 846                       rt[0]-ro[0],rt[1]-ro[1],rt[2]-ro[2]);
 847         t = kFALSE;
 848         gt->GetMatrix(mt);
 849         go->GetMatrix(mo);
 850         for(i=0;i<3;i++)for(j=0;j<3;j++)  t = mt[i][j] != mo[i][j];
 851         if(t){
 852             fprintf(fp,"%1.1d %2.2d %2.2d dfr= %e %e %e\n",
 853                     idt[0],idt[1],idt[2],
 854                     mt[0][0]-mo[0][0],mt[0][1]-mo[0][1],mt[0][2]-mo[0][2]);
 855             fprintf(fp,"        dfr= %e %e %e\n",
 856                     mt[1][0]-mo[1][0],mt[1][1]-mo[1][1],mt[1][2]-mo[1][2]);
 857             fprintf(fp,"        dfr= %e %e %e\n",
 858                     mt[2][0]-mo[2][0],mt[2][1]-mo[2][1],mt[2][2]-mo[2][2]);
 859         } // end if t
 860     } // end for i
 861     return;
 862 }
 863 //______________________________________________________________________
 864 void AliITSgeom::PrintData(FILE *fp,Int_t lay,Int_t lad,Int_t det){
 865     //     This function prints out the coordinate transformations for
 866     // the particular detector defined by layer, ladder, and detector
 867     // to the file pointed to by the File pointer fp. fprintf statements
 868     // are used to print out the numbers. The format is
 869     // layer ladder detector Trans= fx0 fy0 fz0 rot= frx fry frz
 870     // Shape=fShapeIndex
 871     //                         dfr= fr[0] fr[1] fr[2]
 872     //                         dfr= fr[3] fr[4] fr[5]
 873     //                         dfr= fr[6] fr[7] fr[8]
 874     // By indicating which detector, some control over the information
 875     // is given to the user. The output it written to the file pointed
 876     // to by the file pointer fp. This can be set to stdout if you want.
 877     // Inputs are:
 878     // FILE *fp           A file pointer to an opened file for writing in which
 879     //                    the results of the comparison will be written.
 880     // Int_t lay          The layer number. Starting from 1.
 881     // Int_t lad          The ladder number. Starting from 1.
 882     // Int_t det          The detector number. Starting from 1.
 883     // outputs are:
 884     // none
 885     AliITSgeomMatrix *gt;
 886     Double_t t[3],r[3],m[3][3];
 887
 888     gt = this->GetGeomMatrix(GetModuleIndex(lay,lad,det));
 889     gt->GetTranslation(t);
 890     gt->GetAngles(r);
 891     fprintf(fp,"%1.1d %2.2d %2.2d Trans=%f %f %f rot=%f %f %f Shape=%d\n",
 892             lay,lad,det,t[0],t[1],t[2],r[0],r[1],r[2],
 893             gt->GetDetectorIndex());
 894     gt->GetMatrix(m);
 895     fprintf(fp,"        dfr= %e %e %e\n",m[0][0],m[0][1],m[0][2]);
 896     fprintf(fp,"        dfr= %e %e %e\n",m[1][0],m[1][1],m[1][2]);
 897     fprintf(fp,"        dfr= %e %e %e\n",m[2][0],m[2][1],m[2][2]);
 898     return;
 899 }
 900 //______________________________________________________________________
 901 ofstream & AliITSgeom::PrintGeom(ofstream &rb){
 902     //     Stream out an object of class AliITSgeom to standard output.
 903     // Intputs are:
 904     // ofstream &rb    The output streaming buffer.
 905     // Outputs are:
 906     // ofstream &rb    The output streaming buffer.
 907     Int_t i;
 908
 909     rb.setf(ios::scientific);
 910     rb << fTrans << " ";
 911     rb << fNmodules << " ";
 912     rb << fNlayers << " ";
 913     for(i=0;i<fNlayers;i++) rb << fNlad[i] << " ";
 914     for(i=0;i<fNlayers;i++) rb << fNdet[i] << "\n";
 915     for(i=0;i<fNmodules;i++) {
 916         rb <<setprecision(16) << *(GetGeomMatrix(i)) << "\n";
 917     } // end for i
 918     return rb;
 919 }
 920 //______________________________________________________________________
 921 ifstream & AliITSgeom::ReadGeom(ifstream &rb){
 922     //     Stream in an object of class AliITSgeom from standard input.
 923     // Intputs are:
 924     // ifstream &rb    The input streaming buffer.
 925     // Outputs are:
 926     // ifstream &rb    The input streaming buffer.
 927     Int_t i;
 928
 929     fNlad = new Int_t[fNlayers];
 930     fNdet = new Int_t[fNlayers];
 931     if(fGm!=0){
 932         for(i=0;i<fNmodules;i++) delete GetGeomMatrix(i);
 933         delete fGm;
 934     } // end if fGm!=0
 935
 936     rb >> fTrans >> fNmodules >> fNlayers;
 937     fNlad = new Int_t[fNlayers];
 938     fNdet = new Int_t[fNlayers];
 939     for(i=0;i<fNlayers;i++) rb >> fNlad[i];
 940     for(i=0;i<fNlayers;i++) rb >> fNdet[i];
 941     fGm = new TObjArray(fNmodules,0);
 942     for(i=0;i<fNmodules;i++){
 943         fGm->AddAt(new AliITSgeomMatrix,i);
 944         rb >> *(GetGeomMatrix(i));
 945     } // end for i
 946     return rb;
 947 }
 948 //______________________________________________________________________
 949 //     The following routines modify the transformation of "this"
 950 // geometry transformations in a number of different ways.
 951 //______________________________________________________________________
 952 void AliITSgeom::GlobalChange(const Float_t *tran,const Float_t *rot){
 953     //     This function performs a Cartesian translation and rotation of
 954     // the full ITS from its default position by an amount determined by
 955     // the three element arrays tran and rot. If every element
 956     // of tran and rot are zero then there is no change made
 957     // the geometry. The change is global in that the exact same translation
 958     // and rotation is done to every detector element in the exact same way.
 959     // The units of the translation are those of the Monte Carlo, usually cm,
 960     // and those of the rotation are in radians. The elements of tran
 961     // are tran[0] = x, tran[1] = y, and tran[2] = z.
 962     // The elements of rot are rot[0] = rx, rot[1] = ry, and
 963     // rot[2] = rz. A change in x will move the hole ITS in the ALICE
 964     // global x direction, the same for a change in y. A change in z will
 965     // result in a translation of the ITS as a hole up or down the beam line.
 966     // A change in the angles will result in the inclination of the ITS with
 967     // respect to the beam line, except for an effective rotation about the
 968     // beam axis which will just rotate the ITS as a hole about the beam axis.
 969     // Intputs are:
 970     // Float_t *tran   A 3 element array representing the global translations.
 971     //                 the elements are x,y,z in cm.
 972     // Float_t *rot    A 3 element array representing the global rotation
 973     //                 angles about the three axis x,y,z in radians
 974     // Outputs are:
 975     // none.
 976     Int_t    i,j;
 977     Double_t t[3],r[3];
 978     AliITSgeomMatrix *g;
 979
 980     fTrans = (fTrans && 0xfffd) + 2;  // set bit 1 true.
 981     for(i=0;i<fNmodules;i++){
 982         g = this->GetGeomMatrix(i);
 983         g->GetTranslation(t);
 984         g->GetAngles(r);
 985         for(j=0;j<3;j++){
 986             t[j] += tran[j];
 987             r[j] += rot[j];
 988         } // end for j
 989         g->SetTranslation(t);
 990         g->SetAngles(r);
 991     } // end for i
 992     return;
 993 }
 994 //______________________________________________________________________
 995 void AliITSgeom::GlobalCylindericalChange(const Float_t *tran,
 996                                           const Float_t *rot){
 997     //     This function performs a cylindrical translation and rotation of
 998     // each ITS element by a fixed about in radius, rphi, and z from its
 999     // default position by an amount determined by the three element arrays
1000     // tran and rot. If every element of tran and
1001     // rot are zero then there is no change made the geometry. The
1002     // change is global in that the exact same distance change in translation
1003     // and rotation is done to every detector element in the exact same way.
1004     // The units of the translation are those of the Monte Carlo, usually cm,
1005     // and those of the rotation are in radians. The elements of tran
1006     // are tran[0] = r, tran[1] = rphi, and tran[2] = z.
1007     // The elements of rot are rot[0] = rx, rot[1] = ry, and
1008     // rot[2] = rz. A change in r will results in the increase of the
1009     // radius of each layer by the same about. A change in rphi will results in
1010     // the rotation of each layer by a different angle but by the same
1011     // circumferential distance. A change in z will result in a translation
1012     // of the ITS as a hole up or down the beam line. A change in the angles
1013     // will result in the inclination of the ITS with respect to the beam
1014     // line, except for an effective rotation about the beam axis which will
1015     // just rotate the ITS as a hole about the beam axis.
1016     // Intputs are:
1017     // Float_t *tran   A 3 element array representing the global translations.
1018     //                 the elements are r,theta,z in cm/radians.
1019     // Float_t *rot    A 3 element array representing the global rotation
1020     //                 angles about the three axis x,y,z in radians
1021     // Outputs are:
1022     // none.
1023     Int_t    i,j;
1024     Double_t t[3],ro[3],r,r0,phi,rphi;
1025     AliITSgeomMatrix *g;
1026
1027     fTrans = (fTrans && 0xfffd) + 2;  // set bit 1 true.
1028     for(i=0;i<fNmodules;i++){
1029         g = this->GetGeomMatrix(i);
1030         g->GetTranslation(t);
1031         g->GetAngles(ro);
1032         r = r0= TMath::Hypot(t[1],t[0]);
1033         phi   = TMath::ATan2(t[1],t[0]);
1034         rphi  = r0*phi;
1035         r    += tran[0];
1036         rphi += tran[1];
1037         phi   = rphi/r0;
1038         t[0]  = r*TMath::Cos(phi);
1039         t[1]  = r*TMath::Sin(phi);
1040         t[2] += tran[2];
1041         for(j=0;j<3;j++){
1042             ro[j] += rot[j];
1043         } // end for j
1044         g->SetTranslation(t);
1045         g->SetAngles(ro);
1046     } // end for i
1047     return;
1048 }
1049 //______________________________________________________________________
1050 void AliITSgeom::RandomChange(const Float_t *stran,const Float_t *srot){
1051     //     This function performs a Gaussian random displacement and/or
1052     // rotation about the present global position of each active
1053     // volume/detector of the ITS. The sigma of the random displacement
1054     // is determined by the three element array stran, for the
1055     // x y and z translations, and the three element array srot,
1056     // for the three rotation about the axis x y and z.
1057     // Intputs are:
1058     // Float_t *stran  A 3 element array representing the global translations
1059     //                 variances. The elements are x,y,z in cm.
1060     // Float_t *srot   A 3 element array representing the global rotation
1061     //                 angles variances about the three axis x,y,z in radians.
1062     // Outputs are:
1063     // none.
1064     Int_t    i,j;
1065     Double_t t[3],r[3];
1066     AliITSgeomMatrix *g;
1067
1068     fTrans = (fTrans && 0xfffd) + 2;  // set bit 1 true.
1069     for(i=0;i<fNmodules;i++){
1070         g = this->GetGeomMatrix(i);
1071         g->GetTranslation(t);
1072         g->GetAngles(r);
1073         for(j=0;j<3;j++){
1074             t[j] += gRandom->Gaus(0.0,stran[j]);
1075             r[j] += gRandom->Gaus(0.0, srot[j]);
1076         } // end for j
1077         g->SetTranslation(t);
1078         g->SetAngles(r);
1079     } // end for i
1080     return;
1081 }
1082 //______________________________________________________________________
1083 void AliITSgeom::RandomCylindericalChange(const Float_t *stran,
1084                                           const Float_t *srot){
1085     //     This function performs a Gaussian random displacement and/or
1086     // rotation about the present global position of each active
1087     // volume/detector of the ITS. The sigma of the random displacement
1088     // is determined by the three element array stran, for the
1089     // r rphi and z translations, and the three element array srot,
1090     // for the three rotation about the axis x y and z. This random change
1091     // in detector position allow for the simulation of a random uncertainty
1092     // in the detector positions of the ITS.
1093     // Intputs are:
1094     // Float_t *stran  A 3 element array representing the global translations
1095     //                 variances. The elements are r,theta,z in cm/readians.
1096     // Float_t *srot   A 3 element array representing the global rotation
1097     //                 angles variances about the three axis x,y,z in radians.
1098     // Outputs are:
1099     // none.
1100     Int_t    i,j;
1101     Double_t t[3],ro[3],r,r0,phi,rphi;
1102     TRandom ran;
1103     AliITSgeomMatrix *g;
1104
1105     fTrans = (fTrans && 0xfffd) + 2;  // set bit 1 true.
1106     for(i=0;i<fNmodules;i++){
1107         g = this->GetGeomMatrix(i);
1108         g->GetTranslation(t);
1109         g->GetAngles(ro);
1110         r = r0= TMath::Hypot(t[1],t[0]);
1111         phi   = TMath::ATan2(t[1],t[0]);
1112         rphi  = r0*phi;
1113         r    += ran.Gaus(0.0,stran[0]);
1114         rphi += ran.Gaus(0.0,stran[1]);
1115         phi   = rphi/r0;
1116         t[0]  = r*TMath::Cos(phi);
1117         t[1]  = r*TMath::Sin(phi);
1118         t[2] += ran.Gaus(0.0,stran[2]);
1119         for(j=0;j<3;j++){
1120             ro[j] += ran.Gaus(0.0, srot[j]);
1121         } // end for j
1122         g->SetTranslation(t);
1123         g->SetAngles(ro);
1124     } // end for i
1125     return;
1126 }
1127 //______________________________________________________________________
1128 void AliITSgeom::GeantToTracking(AliITSgeom &source){
1129     //     Copy the geometry data but change it to go between the ALICE
1130     // Global coordinate system to that used by the ITS tracking. A slightly
1131     // different coordinate system is used when tracking. This coordinate
1132     // system is only relevant when the geometry represents the cylindrical
1133     // ALICE ITS geometry. For tracking the Z axis is left alone but X-> -Y
1134     // and Y-> X such that X always points out of the ITS cylinder for every
1135     // layer including layer 1 (where the detectors are mounted upside down).
1136     // Inputs are:
1137     // AliITSgeom &source  The AliITSgeom class with which to make this
1138     //                     a copy of.
1139     // Outputs are:
1140     // return  *this       The a new copy of source.
1141     //Begin_Html
1142     /*
1143       <img src="picts/ITS/AliITSgeomMatrix_T1.gif">
1144     */
1145     //End_Html
1146     Int_t    i,j,k,l,id[3];
1147     Double_t r0[3][3],r1[3][3];
1148     Double_t a0[3][3] = {{0.,+1.,0.},{-1.,0.,0.},{0.,0.,+1.}};
1149     Double_t a1[3][3] = {{0.,-1.,0.},{+1.,0.,0.},{0.,0.,+1.}};
1150
1151     *this = source;  // copy everything
1152     for(i=0;i<GetIndexMax();i++){
1153         GetGeomMatrix(i)->GetIndex(id);
1154         GetGeomMatrix(i)->GetMatrix(r0);
1155         if(id[0]==1){ // Layer 1 is treated different from the others.
1156             for(j=0;j<3;j++) for(k=0;k<3;k++){
1157                 r1[j][k] = 0.;
1158                 for(l=0;l<3;l++) r1[j][k] += a0[j][l]*r0[l][k];
1159             } // end for j,k
1160         }else{
1161             for(j=0;j<3;j++) for(k=0;k<3;k++){
1162                 r1[j][k] = 0.;
1163                 for(l=0;l<3;l++) r1[j][k] += a1[j][l]*r0[l][k];
1164             } // end for j,k
1165         } // end if
1166         GetGeomMatrix(i)->SetMatrix(r1);
1167     } // end for i
1168     this->fTrans = (this->fTrans && 0xfffe) + 1;  // set bit 0 true.
1169     return;
1170 }
1171 //______________________________________________________________________
1172 Int_t AliITSgeom::GetNearest(const Double_t g[3],Int_t lay){
1173     //      Finds the Detector (Module) that is nearest the point g [cm] in
1174     // ALICE Global coordinates. If layer !=0 then the search is restricted
1175     // to Detectors (Modules) in that particular layer.
1176     // Inputs are:
1177     // Double_t g[3]  The ALICE Cartesean global coordinate from which the
1178     //                distance is to be calculated with.
1179     // Int_t lay      The layer to restrict the search to. If layer=0 then
1180     //                all layers are searched. Default is lay=0.
1181     // Outputs are:
1182     // return         The module number representing the nearest module.
1183     Int_t    i,l,a,e,in=0;
1184     Double_t d,dn=1.0e10;
1185     Bool_t   t=lay!=0; // skip if lay = 0 default value check all layers.
1186
1187     for(i=0;i<fNmodules;i++){
1188         if(t){GetModuleId(i,l,a,e);if(l!=lay) continue;}
1189         if((d=GetGeomMatrix(i)->Distance2(g))<dn){
1190             dn = d;
1191             in = i;
1192         } // end if
1193     } // end for i
1194     return in;
1195 }
1196 //______________________________________________________________________
1197 void AliITSgeom::GetNearest27(const Double_t g[3],Int_t n[27],Int_t lay){
1198     //      Finds 27 Detectors (Modules) that are nearest the point g [cm] in
1199     // ALICE Global coordinates. If layer !=0 then the search is restricted
1200     // to Detectors (Modules) in that particular layer. The number 27 comes
1201     // from including the nearest detector and all those around it (up, down,
1202     // left, right, forwards, backwards, and the corners).
1203     // Inputs are:
1204     // Double_t g[3]  The ALICE Cartesean global coordinate from which the
1205     //                distance is to be calculated with.
1206     // Int_t lay      The layer to restrict the search to. If layer=0 then
1207     //                all layers are searched. Default is lay=0.
1208     // Outputs are:
1209     // Int_t n[27]    The module number representing the nearest 27 modules
1210     //                in order.
1211     Int_t    i,l,a,e,in[27]={0,0,0,0,0,0,0,0,0,
1212                              0,0,0,0,0,0,0,0,0,
1213                              0,0,0,0,0,0,0,0,0,};
1214     Double_t d,dn[27]={1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,
1215                        1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,
1216                        1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,
1217                        1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,1.0e10,
1218                        1.0e10,1.0e10,1.0e10};
1219     Bool_t   t=(lay!=0); // skip if lay = 0 default value check all layers.
1220
1221     for(i=0;i<fNmodules;i++){
1222         if(t){GetModuleId(i,l,a,e);if(l!=lay) continue;}
1223         for(a=0;a<27;a++){
1224             d = GetGeomMatrix(i)->Distance2(g);
1225             if(d<dn[a]){
1226                 for(e=26;e>a;e--){dn[e] = dn[e-1];in[e] = in[e-1];}
1227                 dn[a] = d; in[a] = i;
1228             } // end if d<dn[i]
1229         } // end for a
1230     } // end for i
1231     for(i=0;i<27;i++) n[i] = in[i];
1232 }
1233 //----------------------------------------------------------------------