Coding rule violations corrected.
[u/mrichter/AliRoot.git] / ITS / AliITSvtest.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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //  Inner Traking System version Test                                        //
21 //  This class contains the base procedures for the Inner Tracking System    //
22 //                                                                           //
23 // Authors: R. Barbera, B. S. Nilsen.                                        //
24 // version  Test                                                             //
25 // Created October 16 1999.                                                  //
26 //                                                                           //
27 ///////////////////////////////////////////////////////////////////////////////
28
29 #include <Riostream.h>
30 #include <stdio.h>
31 #include <stdlib.h>
32 #include <TMath.h>
33 #include <TSystem.h>
34 #include <TVirtualMC.h>
35
36 #include "AliRun.h"
37 #include "AliITSGeant3Geometry.h"
38 #include "AliITSgeom.h"
39 #include "AliITSgeomSDD.h"
40 #include "AliITSgeomSPD.h"
41 #include "AliITSgeomSSD.h"
42 #include "AliITShit.h"
43 #include "AliITSvtest.h"
44 #include "AliRun.h"
45 #include "AliMC.h"
46
47 ClassImp(AliITSvtest)
48  
49 //_____________________________________________________________________________
50 AliITSvtest::AliITSvtest() {
51     // Standard constructor for the ITS
52     Int_t i;
53
54     fIdN    = 0;
55     fIdName = 0;
56     fIdSens = 0;
57     SetEUCLID(kFALSE);
58     fGeomDetOut   = kFALSE; // Don't write .det file
59     fGeomDetIn    = kTRUE; // Read .det file
60     fMajorVersion = IsVersion();
61     fMinorVersion = -1;
62     for(i=0;i<60;i++) fRead[i] = '\0';
63     for(i=0;i<60;i++) fWrite[i] = '\0';
64     for(i=0;i<60;i++) fEuclidGeomDet[i] = '\0';
65 }
66 //____________________________________________________________________________
67 AliITSvtest::AliITSvtest(const AliITSvtest &source) : AliITS(source){
68 ////////////////////////////////////////////////////////////////////////
69 //     Copy Constructor for ITS test version.
70 ////////////////////////////////////////////////////////////////////////
71     if(&source == this) return;
72     Warning("Copy Constructor","Not allowed to copy AliITSvtest");
73     return;
74 }
75 //_____________________________________________________________________________
76 AliITSvtest& AliITSvtest::operator=(const AliITSvtest &source){
77 ////////////////////////////////////////////////////////////////////////
78 //    Assignment operator for the ITS version 1.
79 ////////////////////////////////////////////////////////////////////////
80         if(&source == this) return *this;
81         Warning("= operator","Not allowed to copy AliITSvtest");
82         return *this;
83 }
84 //_____________________________________________________________________________
85 AliITSvtest::~AliITSvtest() {
86     // Standard destructor for the ITS
87 }
88 //_____________________________________________________________________________
89 AliITSvtest::AliITSvtest(const char *fileeuc,const char *filetme,
90                          const char *name, const char *title) 
91     : AliITS(name, title){
92     //
93     // Standard constructor for the ITS
94     //
95     fIdN    = 6;
96     fIdName    = new TString[fIdN];
97     fIdName[0] = "ITS1";
98     fIdName[1] = "ITS2";
99     fIdName[2] = "ITS3";
100     fIdName[3] = "ITS4";
101     fIdName[4] = "ITS5";
102     fIdName[5] = "ITS6";
103     fIdSens    = new Int_t[fIdN];
104     for (Int_t i=0;i<fIdN;i++) fIdSens[i] = 0;
105     fMajorVersion = IsVersion();
106     fMinorVersion = 1;
107     SetEUCLID(kFALSE);
108     fGeomDetOut   = kFALSE; // Don't write .det file
109     fGeomDetIn    = kTRUE; // Read .det file
110
111     fEuclidMaterial = filetme;
112     fEuclidGeometry = fileeuc;
113     strncpy(fEuclidGeomDet,"$ALICE_ROOT/ITS/ITSgeometry_PPR.det",60);
114     strncpy(fRead,fEuclidGeomDet,60);
115     strncpy(fWrite,fEuclidGeomDet,60);
116 //  The .det file for the geometry must have the same name as fileeuc with
117 //  .euc replaced by .det.
118 }
119
120  
121 //_____________________________________________________________________________
122 void AliITSvtest::CreateMaterials(){
123   //
124   // Read materials for the ITS
125   //
126     char *filtmp;
127 //
128   filtmp = gSystem->ExpandPathName(fEuclidMaterial.Data());
129 //  FILE *file = fopen(fEuclidMaterial.Data(),"r");
130   FILE *file = fopen(filtmp,"r");
131   if(file) {
132     fclose(file);
133 //    ReadEuclidMedia(fEuclidMaterial.Data(),this);
134     ReadEuclidMedia(filtmp);
135   } else {
136     Error("CreateMaterials"," THE MEDIA FILE %s DOES NOT EXIST !",
137 //        fEuclidMaterial.Data());
138           filtmp);
139     exit(1);
140   } // end if(file)
141 }
142
143 //_____________________________________________________________________________
144 void AliITSvtest::CreateGeometry(){
145 //////////////////////////////////////////////////////////////////////
146 ////////////////////////////////////////////////////////////////////////
147 // Read geometry for the ITS
148 //
149     char topvol[5];
150     char *filtmp;
151 //
152   filtmp = gSystem->ExpandPathName(fEuclidGeometry.Data());
153   FILE *file = fopen(filtmp,"r");
154   delete [] filtmp;
155   if(file) {
156     fclose(file);
157     printf("Ready to read Euclid geometry file\n");
158     ReadEuclid(fEuclidGeometry.Data(),topvol);
159     printf("Read in euclid geometries\n");
160   } else {
161     Error("CreateGeometry"," THE GEOM FILE %s DOES NOT EXIST !",
162           fEuclidGeometry.Data());
163     exit(1);
164   } // end if(file)
165   //
166   //---Place the ITS ghost volume ITSV in its mother volume (ALIC) and make it
167   //     invisible
168   //
169   gMC->Gspos("ITSV",1,"ALIC",0,0,0,0,"ONLY");
170   //
171   //---Outputs the geometry tree in the EUCLID/CAD format
172   
173     if (fEuclidOut) {
174       gMC->WriteEuclid("ITSgeometry", "ITSV", 1, 5);
175     } // end if (fEuclidOut)
176     cout <<"finished with euclid geometrys"<< endl;
177 }
178 //______________________________________________________________________
179 void AliITSvtest::InitAliITSgeom(){
180 //     Based on the geometry tree defined in Geant 3.21, this
181 // routine initilizes the Class AliITSgeom from the Geant 3.21 ITS geometry
182 // sturture.
183 //    if(gMC->IsA()!=TGeant3::Class()) {
184   if(strcmp(gMC->GetName(),"TGeant3")) {
185         Error("InitAliITSgeom",
186                 "Wrong Monte Carlo. InitAliITSgeom uses TGeant3 calls");
187         return;
188     } // end if
189     cout << "Reading Geometry transformation directly from Geant 3." << endl;
190     const Int_t kNlayers = 6;
191     const Int_t kndeep = 9;
192     Int_t itsGeomTreeNames[kNlayers][kndeep],lnam[20],lnum[20];
193     Int_t nlad[kNlayers],ndet[kNlayers];
194     Double_t t[3],r[10];
195     Float_t  par[20],att[20];
196     Int_t    npar,natt,idshape,imat,imed;
197     AliITSGeant3Geometry *ig = new AliITSGeant3Geometry();
198     Int_t mod,lay,lad,det,i,j,k;
199     const char *names[kNlayers][kndeep] = {
200      {"ALIC","ITSV","ITSD","IT12","I12B","I10B","I107","I101","ITS1"}, // lay=1
201      {"ALIC","ITSV","ITSD","IT12","I12B","I20B","I1D7","I1D1","ITS2"}, // lay=2
202      {"ALIC","ITSV","ITSD","IT34","I004","I302","ITS3","    ","    "}, // lay=3
203      {"ALIC","ITSV","ITSD","IT34","I005","I402","ITS4","    ","    "}, // lay=4
204      {"ALIC","ITSV","ITSD","IT56","I565","I562","ITS5","    ","    "}, // lay=5
205      {"ALIC","ITSV","ITSD","IT56","I569","I566","ITS6","    ","    "}};// lay=6
206     Int_t itsGeomTreeCopys[kNlayers][kndeep] = {{1,1,1,1,10, 2, 4, 1, 1},// lay=1
207                                               {1,1,1,1,10, 4, 4, 1, 1},// lay=2
208                                               {1,1,1,1,14, 6, 1, 0, 0},// lay=3
209                                               {1,1,1,1,22, 8, 1, 0, 0},// lay=4
210                                               {1,1,1,1,34,22, 1, 0, 0},// lay=5
211                                               {1,1,1,1,38,25, 1, 0, 0}};//lay=6
212
213     // Sorry, but this is not very pritty code. It should be replaced
214     // at some point with a version that can search through the geometry
215     // tree its self.
216     cout << "Reading Geometry informaton from Geant3 common blocks" << endl;
217     for(i=0;i<20;i++) lnam[i] = lnum[i] = 0;
218     for(i=0;i<kNlayers;i++)for(j=0;j<kndeep;j++)
219         strncpy((char*) &itsGeomTreeNames[i][j],names[i][j],4); 
220     //  itsGeomTreeNames[i][j] = ig->StringToInt(names[i][j]);
221     mod = 0;
222     for(i=0;i<kNlayers;i++){
223         k = 1;
224         for(j=0;j<kndeep;j++) if(itsGeomTreeCopys[i][j]!=0)
225             k *= TMath::Abs(itsGeomTreeCopys[i][j]);
226         mod += k;
227     } // end for i
228
229     if(fITSgeom!=0) delete fITSgeom;
230     nlad[0]=20;nlad[1]=40;nlad[2]=14;nlad[3]=22;nlad[4]=34;nlad[5]=38;
231     ndet[0]=4;ndet[1]=4;ndet[2]=6;ndet[3]=8;ndet[4]=22;ndet[5]=25;
232     fITSgeom = new AliITSgeom(0,6,nlad,ndet,mod);
233     mod = -1;
234     for(lay=1;lay<=kNlayers;lay++){
235         for(j=0;j<kndeep;j++) lnam[j] = itsGeomTreeNames[lay-1][j];
236         for(j=0;j<kndeep;j++) lnum[j] = itsGeomTreeCopys[lay-1][j];
237         switch (lay){
238         case 1: case 2: // layers 1 and 2 are a bit special
239             lad = 0;
240             for(j=1;j<=itsGeomTreeCopys[lay-1][4];j++){
241                 lnum[4] = j;
242                 for(k=1;k<=itsGeomTreeCopys[lay-1][5];k++){
243                     lad++;
244                     lnum[5] = k;
245                     for(det=1;det<=itsGeomTreeCopys[lay-1][6];det++){
246                         lnum[6] = det;
247                         mod++;
248                         ig->GetGeometry(kndeep,lnam,lnum,t,r,idshape,npar,natt,
249                                         par,att,imat,imed);
250                         fITSgeom->CreatMatrix(mod,lay,lad,det,kSPD,t,r);
251                         if(!(fITSgeom->IsShapeDefined((Int_t)kSPD)))
252                             if(fMinorVersion==1){
253                              fITSgeom->ReSetShape(kSPD,
254                                                   new AliITSgeomSPD425Short());
255                             } else if(fMinorVersion==2)
256                              fITSgeom->ReSetShape(kSPD,
257                                                   new AliITSgeomSPD425Short());
258                     } // end for det
259                 } // end for k
260             } // end for j
261             break;
262         case 3: case 4: case 5: case 6: // layers 3-6
263             lnum[6] = 1;
264             for(lad=1;lad<=itsGeomTreeCopys[lay-1][4];lad++){
265                 lnum[4] = lad;
266                 for(det=1;det<=itsGeomTreeCopys[lay-1][5];det++){
267                     lnum[5] = det;
268                     mod++;
269                     ig->GetGeometry(7,lnam,lnum,t,r,idshape,npar,natt,
270                                     par,att,imat,imed);
271                     switch (lay){
272                     case 3: case 4:
273                         fITSgeom->CreatMatrix(mod,lay,lad,det,kSDD,t,r);
274                         if(!(fITSgeom->IsShapeDefined(kSDD))) 
275                             fITSgeom->ReSetShape(kSDD,new AliITSgeomSDD256());
276                             break;
277                         case 5:
278                             fITSgeom->CreatMatrix(mod,lay,lad,det,kSSD,t,r);
279                             if(!(fITSgeom->IsShapeDefined(kSSD))) 
280                                 fITSgeom->ReSetShape(kSSD,new AliITSgeomSSD275and75());
281                             break;
282                         case 6:
283                             fITSgeom->CreatMatrix(mod,lay,lad,det,kSSD,t,r);
284                             if(!(fITSgeom->IsShapeDefined(kSSD))) 
285                                 fITSgeom->ReSetShape(kSSD,new AliITSgeomSSD75and275());
286                             break;
287                         } // end switch
288                 } // end for det
289             } // end for lad
290             break;
291         } // end switch
292     } // end for lay
293     return;
294 }
295 //_____________________________________________________________________________
296 void AliITSvtest::Init(){
297 ////////////////////////////////////////////////////////////////////////
298 //     Initialise the ITS after it has been created.
299 ////////////////////////////////////////////////////////////////////////
300     Int_t i;
301
302     cout << endl;
303     for(i=0;i<29;i++) cout << "*";cout << " ITSvtest_Init ";
304     for(i=0;i<28;i++) cout << "*";cout << endl;
305 //
306     if(fRead[0]=='\0') strncpy(fRead,fEuclidGeomDet,60);
307     if(fWrite[0]=='\0') strncpy(fWrite,fEuclidGeomDet,60);
308     if(fITSgeom!=0) delete fITSgeom;
309     fITSgeom = new AliITSgeom();
310     if(fGeomDetIn) fITSgeom->ReadNewFile(fRead);
311     if(!fGeomDetIn) this->InitAliITSgeom();
312     if(fGeomDetOut) fITSgeom->WriteNewFile(fWrite);
313     AliITS::Init();
314 //
315     for(i=0;i<72;i++) cout << "*";
316     cout << endl;
317 }
318 //_____________________________________________________________________________
319 void AliITSvtest::StepManager(){
320   //
321   // Called for every step in the ITS
322   //
323   Int_t          copy, id;
324   Int_t          copy1,copy2;
325   Float_t        hits[8];
326   Int_t          vol[4];
327   TLorentzVector position, momentum;
328   TClonesArray   &lhits = *fHits;
329   //
330   // Track status
331   vol[3] = 0;
332   if(gMC->IsTrackInside())      vol[3] +=  1;
333   if(gMC->IsTrackEntering())    vol[3] +=  2;
334   if(gMC->IsTrackExiting())     vol[3] +=  4;
335   if(gMC->IsTrackOut())         vol[3] +=  8;
336   if(gMC->IsTrackDisappeared()) vol[3] += 16;
337   if(gMC->IsTrackStop())        vol[3] += 32;
338   if(gMC->IsTrackAlive())       vol[3] += 64;
339   //
340   // Fill hit structure.
341   if(!(gMC->TrackCharge())) return;
342   //
343   // Only entering charged tracks
344   if((id = gMC->CurrentVolID(copy)) == fIdSens[0]) {
345       vol[0] = 1;
346       id = gMC->CurrentVolOffID(0,copy);
347       //detector copy in the ladder = 1<->4  (ITS1)
348       vol[1] = copy;
349       gMC->CurrentVolOffID(1,copy1);
350       //ladder copy in the module   = 1<->2  (I186)
351       gMC->CurrentVolOffID(2,copy2);
352       //module copy in the layer    = 1<->10 (I132)
353       vol[2] = copy1+(copy2-1)*2;//# of ladders in one module  = 2
354   } else if(id == fIdSens[1]){
355       vol[0] = 2;
356       id = gMC->CurrentVolOffID(0,copy);
357       //detector copy in the ladder = 1<->4  (ITS2)
358       vol[1] = copy;
359       gMC->CurrentVolOffID(1,copy1);
360       //ladder copy in the module   = 1<->4  (I131)
361       gMC->CurrentVolOffID(2,copy2);
362       //module copy in the layer    = 1<->10 (I132)
363       vol[2] = copy1+(copy2-1)*4;//# of ladders in one module  = 4
364   } else if(id == fIdSens[2]){
365       vol[0] = 3;
366       id = gMC->CurrentVolOffID(1,copy);
367       //detector copy in the ladder = 1<->5  (ITS3 is inside I314)
368       vol[1] = copy;
369       id = gMC->CurrentVolOffID(2,copy);
370       //ladder copy in the layer    = 1<->12 (I316)
371       vol[2] = copy;
372   } else if(id == fIdSens[3]){
373       vol[0] = 4;
374       id = gMC->CurrentVolOffID(1,copy);
375       //detector copy in the ladder = 1<->8  (ITS4 is inside I414)
376       vol[1] = copy;
377       id = gMC->CurrentVolOffID(2,copy);
378       //ladder copy in the layer    = 1<->22 (I417)
379       vol[2] = copy;
380   }else if(id == fIdSens[4]){
381       vol[0] = 5;
382       id = gMC->CurrentVolOffID(1,copy);
383       //detector copy in the ladder = 1<->23  (ITS5 is inside I562)
384       vol[1] = copy;
385       id = gMC->CurrentVolOffID(2,copy);
386      //ladder copy in the layer    = 1<->34 (I565)
387       vol[2] = copy;
388   }else if(id == fIdSens[5]){
389       vol[0] = 6;
390       id = gMC->CurrentVolOffID(1,copy);
391       //detector copy in the ladder = 1<->26  (ITS6 is inside I566)
392       vol[1] = copy;
393       id = gMC->CurrentVolOffID(2,copy);
394       //ladder copy in the layer = 1<->38 (I569)
395       vol[2] = copy;
396   } else {
397       return; // not an ITS volume?
398   } // end if/else if (gMC->CurentVolID(copy) == fIdSens[i])
399 //
400   gMC->TrackPosition(position);
401   gMC->TrackMomentum(momentum);
402   hits[0]=position[0];
403   hits[1]=position[1];
404   hits[2]=position[2];
405   hits[3]=momentum[0];
406   hits[4]=momentum[1];
407   hits[5]=momentum[2];
408   hits[6]=gMC->Edep();
409   hits[7]=gMC->TrackTime();
410   // Fill hit structure with this new hit.
411   new(lhits[fNhits++]) AliITShit(fIshunt,gAlice->GetMCApp()->GetCurrentTrackNumber(),vol,hits);
412   return;
413 }
414