don't sort clusters after local reco, do this in AliITSUTrackerGlo
[u/mrichter/AliRoot.git] / ITS / UPGRADE / AliITSUv1.cxx
1 /**************************************************************************
2  * Copyright(c) 2007-2009, 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 /* $Id: AliITSUv1.cxx */
18
19
20 //========================================================================
21 //
22 //        Geometry for the Upgrade of the Inner Tracking System
23 //
24 // Mario Sitta (sitta@to.infn.it)
25 // Chinorat Kobdaj (kobdaj@g.sut.ac.th)
26 //
27 //========================================================================
28
29
30
31 // $Log: AliITSUv1.cxx,v $
32
33 #include <TClonesArray.h>
34 #include <TGeoGlobalMagField.h>
35 #include <TGeoManager.h>
36 #include <TGeoMatrix.h>
37 #include <TGeoPhysicalNode.h>
38 #include <TGeoVolume.h>
39 #include <TGeoTube.h>
40 #include <TGeoXtru.h>
41 #include <TLorentzVector.h>
42 #include <TString.h>
43 #include <TVirtualMC.h>
44
45 #include "AliITSU.h"
46 #include "AliITSUHit.h"
47 #include "AliLog.h"
48 #include "AliMC.h"
49 #include "AliMagF.h"
50 #include "AliRun.h"
51 #include "AliTrackReference.h"
52 #include "AliITSv11Geometry.h"
53 #include "AliITSUv1Layer.h"
54 #include "AliITSUv1.h"
55 #include "AliITSUGeomTGeo.h"
56 #include "AliGeomManager.h"
57 using namespace TMath;
58
59
60 ClassImp(AliITSUv1)
61
62 //______________________________________________________________________
63 AliITSUv1::AliITSUv1()
64 :  fNWrapVol(0)
65   ,fWrapRMin(0)
66   ,fWrapRMax(0)
67   ,fWrapZSpan(0)
68   ,fLay2WrapV(0)
69   ,fLayTurbo(0)
70   ,fLayPhi0(0)
71   ,fLayRadii(0)
72   ,fLayZLength(0)
73   ,fStavPerLay(0)
74   ,fUnitPerStave(0)
75   ,fChipThick(0)
76   ,fStaveWidth(0)
77   ,fStaveTilt(0)
78   ,fSensThick(0)
79   ,fChipTypeID(0)
80   ,fBuildLevel(0)
81   ,fUpGeom(0)
82   ,fStaveModelIB(kIBModel0)
83   ,fStaveModelOB(kOBModel0)
84 {
85   //    Standard default constructor
86   // Inputs:
87   //   none.
88   // Outputs:
89   //   none.
90   // Return:
91   //   none.
92 }
93
94 //______________________________________________________________________
95 AliITSUv1::AliITSUv1(const char *title, Int_t nlay)
96   :AliITSU(title,nlay)
97   ,fNWrapVol(0)
98   ,fWrapRMin(0)
99   ,fWrapRMax(0)
100   ,fWrapZSpan(0)
101   ,fLay2WrapV(0)
102   ,fLayTurbo(0)
103   ,fLayPhi0(0)
104   ,fLayRadii(0)
105   ,fLayZLength(0)
106   ,fStavPerLay(0)
107   ,fUnitPerStave(0)
108   ,fChipThick(0)
109   ,fStaveWidth(0)
110   ,fStaveTilt(0)
111   ,fSensThick(0)
112   ,fChipTypeID(0)
113   ,fBuildLevel(0)
114   ,fUpGeom(0)
115   ,fStaveModelIB(kIBModel0)
116   ,fStaveModelOB(kOBModel0)
117 {
118   //    Standard constructor for the Upgrade geometry.
119   // Inputs:
120   //   const char * name   Ignored, set to "ITS"
121   //   const char * title  Arbitrary title
122   //   const Int_t nlay    Number of layers
123   //
124   fLayerName = new TString[fNLayers];
125   //
126   for (Int_t j=0; j<fNLayers; j++)
127     fLayerName[j].Form("%s%d",AliITSUGeomTGeo::GetITSSensorPattern(),j); // See AliITSUv1Layer
128   //
129   fLayTurbo     = new Bool_t[fNLayers];
130   fLayPhi0      = new Double_t[fNLayers];
131   fLayRadii     = new Double_t[fNLayers];
132   fLayZLength   = new Double_t[fNLayers];
133   fStavPerLay   = new Int_t[fNLayers];
134   fUnitPerStave = new Int_t[fNLayers];
135   fChipThick    = new Double_t[fNLayers];
136   fStaveWidth   = new Double_t[fNLayers];
137   fStaveTilt    = new Double_t[fNLayers];
138   fSensThick    = new Double_t[fNLayers];
139   fChipTypeID   = new UInt_t[fNLayers];
140   fBuildLevel   = new Int_t[fNLayers];
141
142
143   fUpGeom = new AliITSUv1Layer*[fNLayers];
144   
145   if (fNLayers > 0) { // if not, we'll Fatal-ize in CreateGeometry
146     for (Int_t j=0; j<fNLayers; j++) {
147       fLayPhi0[j]      = 0;
148       fLayRadii[j]     = 0.;
149       fLayZLength[j]   = 0.;
150       fStavPerLay[j]   = 0;
151       fUnitPerStave[j] = 0;
152       fStaveWidth[j]   = 0.;
153       fSensThick[j]    = 0.;
154       fChipTypeID[j]   = 0;
155       fBuildLevel[j]   = 0;
156       fUpGeom[j]       = 0;
157     }
158   }
159 }
160
161 //______________________________________________________________________
162 AliITSUv1::~AliITSUv1() {
163   //    Standard destructor
164   // Inputs:
165   //   none.
166   // Outputs:
167   //   none.
168   // Return:
169   //   none.
170   delete [] fLayTurbo;
171   delete [] fLayPhi0;
172   delete [] fLayRadii;
173   delete [] fLayZLength;
174   delete [] fStavPerLay;
175   delete [] fUnitPerStave;
176   delete [] fChipThick;
177   delete [] fStaveWidth;
178   delete [] fStaveTilt;
179   delete [] fSensThick;
180   delete [] fChipTypeID;
181   delete [] fBuildLevel;
182   delete [] fUpGeom;
183   delete [] fWrapRMin;
184   delete [] fWrapRMax;
185   delete [] fWrapZSpan;
186   delete [] fLay2WrapV;
187 }
188
189 //______________________________________________________________________
190 void AliITSUv1::AddAlignableVolumes() const
191 {
192   // Creates entries for alignable volumes associating the symbolic volume
193   // name with the corresponding volume path.
194   // 
195   // Records in the alignable entries the transformation matrices converting
196   // TGeo local coordinates (in the RS of alignable volumes) to the tracking
197   // system
198   // For this, this function has to run before the misalignment because we
199   // are using the ideal positions in the AliITSgeom object.
200   AliInfo("Add ITS alignable volumes");
201
202   if (!gGeoManager) { AliFatal("TGeoManager doesn't exist !"); return;  }
203   //
204   TString path = Form("ALIC_1/%s_2",AliITSUGeomTGeo::GetITSVolPattern());
205   TString sname = AliITSUGeomTGeo::ComposeSymNameITS();
206   //
207   AliDebug(1,Form("%s <-> %s",sname.Data(),path.Data()));
208   if( !gGeoManager->SetAlignableEntry(sname.Data(),path.Data()) )
209     AliFatal(Form("Unable to set alignable entry ! %s %s",sname.Data(),path.Data()));    
210   //
211   int lastUID = 0;
212   for (int lr=0; lr<fNLayers; lr++) AddAlignableVolumesLayer(lr,path,lastUID);
213   //
214 }
215
216 //______________________________________________________________________
217 void AliITSUv1::AddAlignableVolumesLayer(int lr, TString& parent,Int_t &lastUID) const
218 {
219   // add alignable volumes for layer and its daughters
220   //
221   TString wrpV = fLay2WrapV[lr]!=-1 ? Form("%s%d_1/",AliITSUGeomTGeo::GetITSWrapVolPattern(),fLay2WrapV[lr]) : "";
222   TString path = Form("%s/%s%s%d_1",parent.Data(),wrpV.Data(),AliITSUGeomTGeo::GetITSLayerPattern(),lr);
223   TString sname = AliITSUGeomTGeo::ComposeSymNameLayer(lr);
224   AliDebug(1,Form("Add %s <-> %s", sname.Data(),path.Data()));
225   if ( !gGeoManager->SetAlignableEntry(sname.Data(),path.Data()) ) 
226     AliFatal(Form("Unable to set alignable entry ! %s : %s",sname.Data(),path.Data()));
227   //
228   const AliITSUv1Layer* lrobj = fUpGeom[lr];
229   int nstaves = lrobj->GetNStavesPerParent();
230   for (int st=0; st<nstaves; st++) AddAlignableVolumesStave(lr,st,path,lastUID);
231   //
232 }
233     
234 //______________________________________________________________________
235 void AliITSUv1::AddAlignableVolumesStave(int lr, int st, TString& parent, Int_t &lastUID) const
236 {
237   // add alignable volumes for stave and its daughters
238   //
239   TString path = Form("%s/%s%d_%d",parent.Data(),AliITSUGeomTGeo::GetITSStavePattern(),lr,st);
240   TString sname = AliITSUGeomTGeo::ComposeSymNameStave(lr,st);
241   AliDebug(1,Form("Add %s <-> %s", sname.Data(),path.Data()));
242   if ( !gGeoManager->SetAlignableEntry(sname.Data(),path.Data()) ) 
243     AliFatal(Form("Unable to set alignable entry ! %s : %s",sname.Data(),path.Data()));
244   //
245   const AliITSUv1Layer* lrobj = fUpGeom[lr];
246   int nsstave = lrobj->GetNHalfStavesPerParent();
247   int start = nsstave>0 ? 0:-1;
248   //
249   for (int sst=start; sst<nsstave; sst++) AddAlignableVolumesHalfStave(lr,st,sst,path,lastUID);
250 }
251
252 //______________________________________________________________________
253 void AliITSUv1::AddAlignableVolumesHalfStave(int lr, int st, int sst, TString& parent, Int_t &lastUID) const
254 {
255   // add alignable volumes for halfstave (if any) and its daughters
256   //
257   TString path = parent;
258   if (sst>=0) {
259     path = Form("%s/%s%d_%d",parent.Data(),AliITSUGeomTGeo::GetITSHalfStavePattern(),lr,sst);
260     TString sname = AliITSUGeomTGeo::ComposeSymNameHalfStave(lr,st,sst);
261     AliDebug(1,Form("Add %s <-> %s", sname.Data(),path.Data()));
262     if ( !gGeoManager->SetAlignableEntry(sname.Data(),path.Data()) ) 
263       AliFatal(Form("Unable to set alignable entry ! %s : %s",sname.Data(),path.Data()));
264     //
265   }
266   const AliITSUv1Layer* lrobj = fUpGeom[lr];
267   int nmodules = lrobj->GetNModulesPerParent();
268   int start = nmodules>0 ? 0:-1;
269   //
270   for (int md=start; md<nmodules; md++) AddAlignableVolumesModule(lr,st,sst,md,path,lastUID);
271 }
272
273 //______________________________________________________________________
274 void AliITSUv1::AddAlignableVolumesModule(int lr, int st, int sst, int md, TString& parent, Int_t &lastUID) const
275 {
276   // add alignable volumes for module (if any) and its daughters
277   //
278   TString path = parent;
279   if (md>=0) {
280     path = Form("%s/%s%d_%d",parent.Data(),AliITSUGeomTGeo::GetITSModulePattern(),lr,md);
281     TString sname = AliITSUGeomTGeo::ComposeSymNameModule(lr,st,sst,md);
282     AliDebug(1,Form("Add %s <-> %s", sname.Data(),path.Data()));
283     if ( !gGeoManager->SetAlignableEntry(sname.Data(),path.Data()) ) 
284       AliFatal(Form("Unable to set alignable entry ! %s : %s",sname.Data(),path.Data()));
285     //
286   }
287   //
288   const AliITSUv1Layer* lrobj = fUpGeom[lr];
289   int nchips = lrobj->GetNChipsPerParent();
290   //
291   for (int ic=0; ic<nchips; ic++) AddAlignableVolumesChip(lr,st,sst,md,ic,path,lastUID);
292 }  
293
294 //______________________________________________________________________
295 void AliITSUv1::AddAlignableVolumesChip(int lr, int st, int sst, int md, int ch, TString& parent, Int_t &lastUID) const
296 {
297   // add alignable volumes for chip
298   //
299   TString path = Form("%s/%s%d_%d",parent.Data(),AliITSUGeomTGeo::GetITSChipPattern(),lr,ch);
300   TString sname = AliITSUGeomTGeo::ComposeSymNameChip(lr,st,sst,md,ch);
301   int modUID = AliITSUGeomTGeo::ChipVolUID( lastUID++ );
302   //
303   AliDebug(1,Form("Add %s <-> %s : ID=%d", sname.Data(),path.Data(),modUID));
304   if ( !gGeoManager->SetAlignableEntry(sname,path.Data(),modUID) )
305     AliFatal(Form("Unable to set alignable entry ! %s : %s | %d",sname.Data(),path.Data(),modUID));
306   //
307 }
308
309 //______________________________________________________________________
310 void AliITSUv1::SetNWrapVolumes(Int_t n)
311 {
312   // book arrays for wrapper volumes
313   if (fNWrapVol) AliFatal(Form("%d wrapper volumes already defined",fNWrapVol));
314   if (n<1) return;
315   fNWrapVol = n;
316   fWrapRMin = new Double_t[fNWrapVol];
317   fWrapRMax = new Double_t[fNWrapVol];
318   fWrapZSpan= new Double_t[fNWrapVol];
319   for (int i=fNWrapVol;i--;) fWrapRMin[i]=fWrapRMax[i]=fWrapZSpan[i]=-1;
320   //
321 }
322
323 //______________________________________________________________________
324 void AliITSUv1::DefineWrapVolume(Int_t id, Double_t rmin,Double_t rmax, Double_t zspan)
325 {
326   // set parameters of id-th wrapper volume
327   if (id>=fNWrapVol||id<0) AliFatal(Form("id=%d of wrapper volume is not in 0-%d range",id,fNWrapVol-1));
328   fWrapRMin[id] = rmin;
329   fWrapRMax[id] = rmax;
330   fWrapZSpan[id] = zspan;
331 }
332
333 //______________________________________________________________________
334 void AliITSUv1::CreateGeometry() {
335
336   // Create the geometry and insert it in the mother volume ITSV
337   TGeoManager *geoManager = gGeoManager;
338
339   TGeoVolume *vALIC = geoManager->GetVolume("ALIC");
340
341   new TGeoVolumeAssembly(AliITSUGeomTGeo::GetITSVolPattern());
342   TGeoVolume *vITSV = geoManager->GetVolume(AliITSUGeomTGeo::GetITSVolPattern());
343   vITSV->SetUniqueID(AliITSUGeomTGeo::GetUIDShift()); // store modID -> midUUID bitshift
344   vALIC->AddNode(vITSV, 2, 0);  // Copy number is 2 to cheat AliGeoManager::CheckSymNamesLUT
345   //
346   const Int_t kLength=100;
347   Char_t vstrng[kLength] = "xxxRS"; //?
348   vITSV->SetTitle(vstrng);
349   //
350   // Check that we have all needed parameters
351   if (fNLayers <= 0) AliFatal(Form("Wrong number of layers (%d)",fNLayers));
352   //
353   for (Int_t j=0; j<fNLayers; j++) {
354     if (fLayRadii[j] <= 0)                   AliFatal(Form("Wrong layer radius for layer %d (%f)",j,fLayRadii[j]));
355     if (fLayZLength[j] <= 0)                 AliFatal(Form("Wrong layer length for layer %d (%f)",j,fLayZLength[j]));
356     if (fStavPerLay[j] <= 0)                 AliFatal(Form("Wrong number of staves for layer %d (%d)",j,fStavPerLay[j]));
357     if (fUnitPerStave[j] <= 0)               AliFatal(Form("Wrong number of chips for layer %d (%d)",j,fUnitPerStave[j]));
358     if (fChipThick[j] < 0)                   AliFatal(Form("Wrong chip thickness for layer %d (%f)",j,fChipThick[j]));
359     if (fLayTurbo[j] && fStaveWidth[j] <= 0) AliFatal(Form("Wrong stave width for layer %d (%f)",j,fStaveWidth[j]));
360     if (fSensThick[j] < 0)                   AliFatal(Form("Wrong sensor thickness for layer %d (%f)",j,fSensThick[j]));
361     //
362     if (j > 0) {
363       if (fLayRadii[j]<=fLayRadii[j-1])      AliFatal(Form("Layer %d radius (%f) is smaller than layer %d radius (%f)",
364                            j,fLayRadii[j],j-1,fLayRadii[j-1]));
365     } // if (j > 0)
366
367     if (fChipThick[j] == 0) AliInfo(Form("Chip thickness for layer %d not set, using default",j));
368     if (fSensThick[j] == 0) AliInfo(Form("Sensor thickness for layer %d not set, using default",j));
369
370   } // for (Int_t j=0; j<fNLayers; j++)
371
372   // Create the wrapper volumes
373   TGeoVolume **wrapVols = 0;
374   if (fNWrapVol) {
375     wrapVols = new TGeoVolume*[fNWrapVol];
376     for (int id=0;id<fNWrapVol;id++) {
377       wrapVols[id] = CreateWrapperVolume(id);
378       vITSV->AddNode(wrapVols[id], 1, 0);
379     }
380   }
381   //
382   fLay2WrapV = new Int_t[fNLayers];
383
384   // Now create the actual geometry
385   for (Int_t j=0; j<fNLayers; j++) {
386     TGeoVolume* dest = vITSV;
387     fLay2WrapV[j] = -1;
388     //
389     if (fLayTurbo[j]) {
390       fUpGeom[j] = new AliITSUv1Layer(j,kTRUE,kFALSE);
391       fUpGeom[j]->SetStaveWidth(fStaveWidth[j]);
392       fUpGeom[j]->SetStaveTilt(fStaveTilt[j]);
393     }
394     else fUpGeom[j] = new AliITSUv1Layer(j,kFALSE);
395     //
396     fUpGeom[j]->SetPhi0(fLayPhi0[j]);
397     fUpGeom[j]->SetRadius(fLayRadii[j]);
398     fUpGeom[j]->SetZLength(fLayZLength[j]);
399     fUpGeom[j]->SetNStaves(fStavPerLay[j]);
400     fUpGeom[j]->SetNUnits(fUnitPerStave[j]);
401     fUpGeom[j]->SetChipType(fChipTypeID[j]);
402     fUpGeom[j]->SetBuildLevel(fBuildLevel[j]);
403     if (j < 3)
404       fUpGeom[j]->SetStaveModel(fStaveModelIB);
405     else
406       fUpGeom[j]->SetStaveModel(fStaveModelOB);
407     AliDebug(1,Form("fBuildLevel: %d\n",fBuildLevel[j]));
408     //
409     if (fChipThick[j] != 0) fUpGeom[j]->SetChipThick(fChipThick[j]);
410     if (fSensThick[j] != 0) fUpGeom[j]->SetSensorThick(fSensThick[j]);
411     //
412     for (int iw=0;iw<fNWrapVol;iw++) {
413       if (fLayRadii[j]>fWrapRMin[iw] && fLayRadii[j]<fWrapRMax[iw]) {
414         AliInfo(Form("Will embed layer %d in wrapper volume %d",j,iw));
415         if (fLayZLength[j]>=fWrapZSpan[iw]) AliFatal(Form("ZSpan %.3f of wrapper volume %d is less than ZSpan %.3f of layer %d",
416                                                           fWrapZSpan[iw],iw,fLayZLength[j],j));
417         dest = wrapVols[iw];
418         fLay2WrapV[j] = iw;
419         break; 
420       }
421     }
422     fUpGeom[j]->CreateLayer(dest);
423   }
424   CreateSuppCyl(kTRUE,wrapVols[0]);
425   CreateSuppCyl(kFALSE,wrapVols[2]);
426
427   delete[] wrapVols; // delete pointer only, not the volumes
428   //
429 }
430
431 //____________________________________________________________
432 //Service Barrel
433 void AliITSUv1::CreateSuppCyl(const Bool_t innerBarrel,TGeoVolume *dest,const TGeoManager *mgr){
434   // Creates the Service Barrel (as a simple cylinder) for IB and OB
435   // Inputs:
436   //         innerBarrel : if true, build IB service barrel, otherwise for OB
437   //         dest        : the mother volume holding the service barrel
438   //         mgr         : the gGeoManager pointer (used to get the material)
439   //
440
441   Double_t rminIB =  4.7;
442   Double_t rminOB = 43.4;
443   Double_t zLenOB ;
444   Double_t cInt = 0.22; //dimensioni cilindro di supporto interno
445   Double_t cExt = 1.00; //dimensioni cilindro di supporto esterno
446 //  Double_t phi1   =  180;
447 //  Double_t phi2   =  360;
448
449
450   TGeoMedium *medCarbonFleece = mgr->GetMedium("ITS_CarbonFleece$");
451
452   if (innerBarrel){
453     zLenOB=((TGeoTube*)(dest->GetShape()))->GetDz();
454 //    TGeoTube*ibSuppSh = new TGeoTubeSeg(rminIB,rminIB+cInt,zLenOB,phi1,phi2);
455     TGeoTube*ibSuppSh = new TGeoTube(rminIB,rminIB+cInt,zLenOB);
456     TGeoVolume *ibSupp = new TGeoVolume("ibSuppCyl",ibSuppSh,medCarbonFleece);
457     dest->AddNode(ibSupp,1);
458   }
459   else {
460     zLenOB=((TGeoTube*)(dest->GetShape()))->GetDz();
461     TGeoTube*obSuppSh=new TGeoTube(rminOB,rminOB+cExt,zLenOB);
462     TGeoVolume *obSupp=new TGeoVolume("obSuppCyl",obSuppSh,medCarbonFleece);
463     dest->AddNode(obSupp,1);
464   }
465
466   return;
467 }
468
469 //______________________________________________________________________
470 void AliITSUv1::CreateMaterials() {
471   // Create ITS materials
472   //     This function defines the default materials used in the Geant
473   // Monte Carlo simulations for the geometries AliITSv1, AliITSv3,
474   // AliITSv11Hybrid.
475   // In general it is automatically replaced by
476   // the CreateMaterials routine defined in AliITSv?. Should the function
477   // CreateMaterials not exist for the geometry version you are using this
478   // one is used. See the definition found in AliITSv5 or the other routine
479   // for a complete definition.
480   // Inputs:
481   //   none.
482   // Outputs:
483   //   none.
484   // Return:
485   //   none.
486
487   Int_t   ifield = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ();
488   Float_t fieldm = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();
489
490   Float_t tmaxfd = 0.1; // 1.0; // Degree
491   Float_t stemax = 1.0; // cm
492   Float_t deemax = 0.1; // 30.0; // Fraction of particle's energy 0<deemax<=1
493   Float_t epsil  = 1.0E-4; // 1.0; // cm
494   Float_t stmin  = 0.0; // cm "Default value used"
495
496   Float_t tmaxfdSi = 0.1; // .10000E+01; // Degree
497   Float_t stemaxSi = 0.0075; //  .10000E+01; // cm
498   Float_t deemaxSi = 0.1; // 0.30000E-02; // Fraction of particle's energy 0<deemax<=1
499   Float_t epsilSi  = 1.0E-4;// .10000E+01;
500   Float_t stminSi  = 0.0; // cm "Default value used"
501
502   Float_t tmaxfdAir = 0.1; // .10000E+01; // Degree
503   Float_t stemaxAir = .10000E+01; // cm
504   Float_t deemaxAir = 0.1; // 0.30000E-02; // Fraction of particle's energy 0<deemax<=1
505   Float_t epsilAir  = 1.0E-4;// .10000E+01;
506   Float_t stminAir  = 0.0; // cm "Default value used"
507
508   // AIR
509   Float_t aAir[4]={12.0107,14.0067,15.9994,39.948};
510   Float_t zAir[4]={6.,7.,8.,18.};
511   Float_t wAir[4]={0.000124,0.755267,0.231781,0.012827};
512   Float_t dAir = 1.20479E-3;
513
514   // Water
515   Float_t aWater[2]={1.00794,15.9994};
516   Float_t zWater[2]={1.,8.};
517   Float_t wWater[2]={0.111894,0.888106};
518   Float_t dWater   = 1.0;
519
520
521   // Kapton
522   Float_t aKapton[4]={1.00794,12.0107, 14.010,15.9994};
523   Float_t zKapton[4]={1.,6.,7.,8.};
524   Float_t wKapton[4]={0.026362,0.69113,0.07327,0.209235};
525   Float_t dKapton   = 1.42;
526  
527   AliMixture(1,"AIR$",aAir,zAir,dAir,4,wAir);
528   AliMedium(1, "AIR$",1,0,ifield,fieldm,tmaxfdAir,stemaxAir,deemaxAir,epsilAir,stminAir);
529
530   AliMixture(2,"WATER$",aWater,zWater,dWater,2,wWater);
531   AliMedium(2, "WATER$",2,0,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin);
532
533   AliMaterial(3,"SI$",0.28086E+02,0.14000E+02,0.23300E+01,0.93600E+01,0.99900E+03);
534   AliMedium(3,  "SI$",3,0,ifield,fieldm,tmaxfdSi,stemaxSi,deemaxSi,epsilSi,stminSi);
535
536   AliMaterial(4,"BERILLIUM$",9.01, 4., 1.848, 35.3, 36.7);// From AliPIPEv3
537   AliMedium(4,  "BERILLIUM$",4,0,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin);
538
539   AliMaterial(5,"COPPER$",0.63546E+02,0.29000E+02,0.89600E+01,0.14300E+01,0.99900E+03);
540   AliMedium(5,  "COPPER$",5,0,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin);
541
542     
543   // needed for STAVE , Carbon, kapton, Epoxy, flexcable
544
545   //AliMaterial(6,"CARBON$",12.0107,6,2.210,999,999);
546   AliMaterial(6,"CARBON$",12.0107,6,2.210/1.3,999,999);
547   AliMedium(6,  "CARBON$",6,0,ifield,fieldm,tmaxfdSi,stemaxSi,deemaxSi,epsilSi,stminSi);
548
549   AliMixture(7,"KAPTON(POLYCH2)$", aKapton, zKapton, dKapton, 4, wKapton);
550   AliMedium(7, "KAPTON(POLYCH2)$",7,0,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin);
551
552
553  
554   // values below modified as compared to source AliITSv11 !
555
556   //AliMaterial(7,"GLUE$",0.12011E+02,0.60000E+01,0.1930E+01/2.015,999,999); // original
557   AliMaterial(15,"GLUE$",12.011,6,1.93/2.015,999,999);  // conform with ATLAS, Corrado, Stefan
558   AliMedium(15,  "GLUE$",15,0,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin);
559
560   // All types of carbon
561   // Unidirectional prepreg
562  AliMaterial(8,"K13D2U2k$",12.0107,6,1.643,999,999);
563  AliMedium(8,  "K13D2U2k$",8,0,ifield,fieldm,tmaxfdSi,stemaxSi,deemaxSi,epsilSi,stminSi);
564  //Impregnated thread
565  AliMaterial(9,"M60J3K$",12.0107,6,2.21,999,999);
566  AliMedium(9,  "M60J3K$",9,0,ifield,fieldm,tmaxfdSi,stemaxSi,deemaxSi,epsilSi,stminSi);
567  //Impregnated thread
568  AliMaterial(10,"M55J6K$",12.0107,6,1.63,999,999);
569  AliMedium(10,  "M55J6K$",10,0,ifield,fieldm,tmaxfdSi,stemaxSi,deemaxSi,epsilSi,stminSi);
570  // Fabric(0/90)
571  AliMaterial(11,"T300$",12.0107,6,1.725,999,999);
572  AliMedium(11,  "T300$",11,0,ifield,fieldm,tmaxfdSi,stemaxSi,deemaxSi,epsilSi,stminSi);
573  //AMEC Thermasol
574  AliMaterial(12,"FGS003$",12.0107,6,1.6,999,999);
575  AliMedium(12,  "FGS003$",12,0,ifield,fieldm,tmaxfdSi,stemaxSi,deemaxSi,epsilSi,stminSi);
576  // Carbon fleece
577  AliMaterial(13,"CarbonFleece$",12.0107,6,0.4,999,999);
578  AliMedium(13,  "CarbonFleece$",13,0,ifield,fieldm,tmaxfdSi,stemaxSi,deemaxSi,epsilSi,stminSi);
579
580   // Flex cable
581   Float_t aFCm[5]={12.0107,1.00794,14.0067,15.9994,26.981538};
582   Float_t zFCm[5]={6.,1.,7.,8.,13.};
583   Float_t wFCm[5]={0.520088819984,0.01983871336,0.0551367996,0.157399667056, 0.247536};
584   //Float_t dFCm = 1.6087;  // original
585   //Float_t dFCm = 2.55;   // conform with STAR
586    Float_t dFCm = 2.595;   // conform with Corrado
587
588   AliMixture(14,"FLEXCABLE$",aFCm,zFCm,dFCm,5,wFCm);
589   AliMedium(14, "FLEXCABLE$",14,0,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin);
590
591   AliMaterial(16,"ALUMINUM$",0.26982E+02,0.13000E+02,0.26989E+01,0.89000E+01,0.99900E+03);
592   AliMedium(16,"ALUMINUM$",16,0,ifield,fieldm,tmaxfd,stemax,deemax,epsil,stmin);
593
594 }
595
596 //______________________________________________________________________
597 void AliITSUv1::DefineLayer(Int_t nlay, double phi0, Double_t r,
598                             Double_t zlen, Int_t nstav,
599                             Int_t nunit, Double_t lthick,
600                             Double_t dthick, UInt_t dettypeID,
601                             Int_t buildLevel)
602 {
603   //     Sets the layer parameters
604   // Inputs:
605   //          nlay    layer number
606   //          phi0    layer phi0
607   //          r       layer radius
608   //          zlen    layer length
609   //          nstav   number of staves
610   //          nunit   IB: number of chips per stave
611   //                  OB: number of modules per half stave
612   //          lthick  chip thickness (if omitted, defaults to 0)
613   //          dthick  sensor thickness (if omitted, defaults to 0)
614   //          dettypeID  ??
615   //          buildLevel (if 0, all geometry is build, used for material budget studies)
616   // Outputs:
617   //   none.
618   // Return:
619   //   none.
620
621   AliInfo(Form("L# %d Phi:%+.3f R:%+7.3f DZ:%7.2f Nst:%2d Nunit:%2d Lthick:%.4f Dthick:%.4f DetID:%d B:%d",
622                nlay,phi0,r,zlen,nstav,nunit,lthick,dthick,dettypeID,buildLevel));
623
624   if (nlay >= fNLayers || nlay < 0) {
625     AliError(Form("Wrong layer number (%d)",nlay));
626     return;
627   }
628   
629   fLayTurbo[nlay] = kFALSE;
630   fLayPhi0[nlay] = phi0;
631   fLayRadii[nlay] = r;
632   fLayZLength[nlay] = zlen;
633   fStavPerLay[nlay] = nstav;
634   fUnitPerStave[nlay] = nunit;
635   fChipThick[nlay] = lthick;
636   fSensThick[nlay] = dthick;
637   fChipTypeID[nlay] = dettypeID;
638   fBuildLevel[nlay] = buildLevel;
639     
640 }
641
642 //______________________________________________________________________
643 void AliITSUv1::DefineLayerTurbo(Int_t nlay, Double_t phi0, Double_t r, Double_t zlen, Int_t nstav,
644                                  Int_t nunit, Double_t width, Double_t tilt,
645                                  Double_t lthick,Double_t dthick,
646                                  UInt_t dettypeID, Int_t buildLevel)
647 {
648   //     Sets the layer parameters for a "turbo" layer
649   //     (i.e. a layer whose staves overlap in phi)
650   // Inputs:
651   //          nlay    layer number
652   //          phi0    phi of 1st stave
653   //          r       layer radius
654   //          zlen    layer length
655   //          nstav   number of staves
656   //          nunit   IB: number of chips per stave
657   //                  OB: number of modules per half stave
658   //          width   stave width
659   //          tilt    layer tilt angle (degrees)
660   //          lthick  chip thickness (if omitted, defaults to 0)
661   //          dthick  sensor thickness (if omitted, defaults to 0)
662   //          dettypeID  ??
663   //          buildLevel (if 0, all geometry is build, used for material budget studies)
664   // Outputs:
665   //   none.
666   // Return:
667   //   none.
668
669   AliInfo(Form("L# %d Phi:%+.3f R:%+7.3f DZ:%7.2f Nst:%2d Nunit:%2d W:%7.4f Tilt:%+.3f Lthick:%.4f Dthick:%.4f DetID:%d B:%d",
670                nlay,phi0,r,zlen,nstav,nunit,width,tilt,lthick,dthick,dettypeID,buildLevel));
671
672   if (nlay >= fNLayers || nlay < 0) {
673     AliError(Form("Wrong layer number (%d)",nlay));
674     return;
675   }
676
677   fLayTurbo[nlay] = kTRUE;
678   fLayPhi0[nlay] = phi0;
679   fLayRadii[nlay] = r;
680   fLayZLength[nlay] = zlen;
681   fStavPerLay[nlay] = nstav;
682   fUnitPerStave[nlay] = nunit;
683   fChipThick[nlay] = lthick;
684   fStaveWidth[nlay] = width;
685   fStaveTilt[nlay] = tilt;
686   fSensThick[nlay] = dthick;
687   fChipTypeID[nlay] = dettypeID;
688   fBuildLevel[nlay] = buildLevel;
689
690 }
691
692 //______________________________________________________________________
693 void AliITSUv1::GetLayerParameters(Int_t nlay, Double_t &phi0,
694                                    Double_t &r, Double_t &zlen,
695                                    Int_t &nstav, Int_t &nmod,
696                                    Double_t &width, Double_t &tilt,
697                                    Double_t &lthick, Double_t &dthick,
698                                    UInt_t &dettype) const
699 {
700   //     Gets the layer parameters
701   // Inputs:
702   //          nlay    layer number
703   // Outputs:
704   //          phi0    phi of 1st stave
705   //          r       layer radius
706   //          zlen    layer length
707   //          nstav   number of staves
708   //          nmod    IB: number of chips per stave
709   //                  OB: number of modules per half stave
710   //          width   stave width
711   //          tilt    stave tilt angle
712   //          lthick  chip thickness
713   //          dthick  sensor thickness
714   //          dettype detector type
715   // Return:
716   //   none.
717
718   if (nlay >= fNLayers || nlay < 0) {
719     AliError(Form("Wrong layer number (%d)",nlay));
720     return;
721   }
722   
723   phi0   = fLayPhi0[nlay];
724   r      = fLayRadii[nlay];
725   zlen   = fLayZLength[nlay];
726   nstav  = fStavPerLay[nlay];
727   nmod   = fUnitPerStave[nlay];
728   width  = fStaveWidth[nlay];
729   tilt   = fStaveTilt[nlay];
730   lthick = fChipThick[nlay];
731   dthick = fSensThick[nlay];
732   dettype= fChipTypeID[nlay];
733 }
734
735 //______________________________________________________________________
736 TGeoVolume* AliITSUv1::CreateWrapperVolume(Int_t id)
737 {
738   //     Creates an air-filled wrapper cylindrical volume 
739   // Inputs:
740   //          volume id
741   // Outputs:
742   //          the wrapper volume
743
744   if (fWrapRMin[id]<0 || fWrapRMax[id]<0 || fWrapZSpan[id]<0) AliFatal(Form("Wrapper volume %d was requested but not defined",id));
745   // Now create the actual shape and volume
746   //
747   TGeoTube *tube = new TGeoTube(fWrapRMin[id], fWrapRMax[id], fWrapZSpan[id]/2.);
748
749   TGeoMedium *medAir = gGeoManager->GetMedium("ITS_AIR$");
750
751   char volnam[30];
752   snprintf(volnam, 29, "%s%d", AliITSUGeomTGeo::GetITSWrapVolPattern(),id);
753
754   TGeoVolume *wrapper = new TGeoVolume(volnam, tube, medAir);
755
756   return wrapper;
757 }
758
759 //______________________________________________________________________
760 void AliITSUv1::Init()
761 {
762   //     Initialise the ITS after it has been created.
763   UpdateInternalGeometry();
764   AliITSU::Init();
765   //  
766 }
767
768 //______________________________________________________________________
769 Bool_t AliITSUv1::IsLayerTurbo(Int_t nlay)
770 {
771   //     Returns true if the layer is a "turbo" layer
772   if ( nlay < 0 || nlay > fNLayers ) {
773     AliError(Form("Wrong layer number %d",nlay));
774     return kFALSE;
775   } 
776   else return fUpGeom[nlay]->IsTurbo();
777 }
778
779 //______________________________________________________________________
780 void AliITSUv1::SetDefaults()
781 {
782   // sets the default segmentation, response, digit and raw cluster classes
783 }
784
785 //______________________________________________________________________
786 void AliITSUv1::StepManager()
787 {
788   //    Called for every step in the ITS, then calles the AliITSUHit class
789   // creator with the information to be recoreded about that hit.
790   //     The value of the macro ALIITSPRINTGEOM if set to 1 will allow the
791   // printing of information to a file which can be used to create a .det
792   // file read in by the routine CreateGeometry(). If set to 0 or any other
793   // value except 1, the default behavior, then no such file is created nor
794   // it the extra variables and the like used in the printing allocated.
795   // Inputs:
796   //   none.
797   // Outputs:
798   //   none.
799   // Return:
800   //   none.
801   //
802   if(!(this->IsActive())) return;
803   if(!(TVirtualMC::GetMC()->TrackCharge())) return;
804   //
805   Int_t copy, lay = 0;
806   Int_t id = TVirtualMC::GetMC()->CurrentVolID(copy);
807
808   Bool_t notSens = kFALSE;
809   while ((lay<fNLayers)  && (notSens = (id!=fIdSens[lay]))) ++lay;
810   //printf("R: %.1f | Lay: %d  NotSens: %d\n",positionRS.Pt(), lay, notSens);
811            
812   if (notSens) return;
813   //
814   if (lay < 0 || lay >= fNLayers) {
815     AliError(Form("Invalid value: lay=%d. Not an ITS sensitive volume",lay));
816     return; // not an ITS sensitive volume.
817   } 
818   //
819   static TLorentzVector position, momentum; // Saves on calls to construtors
820   static AliITSUHit hit;// Saves on calls to constructors
821   //
822   TClonesArray &lhits = *(Hits());
823   Int_t chipID, status = 0;
824   //
825   // Track status
826   if(TVirtualMC::GetMC()->IsTrackInside())      status +=  1;
827   if(TVirtualMC::GetMC()->IsTrackEntering())    status +=  2;
828   if(TVirtualMC::GetMC()->IsTrackExiting()) {
829     AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kITS);
830     status +=  4;
831   } // if Outer ITS mother Volume
832   if(TVirtualMC::GetMC()->IsTrackOut())         status +=  8;
833   if(TVirtualMC::GetMC()->IsTrackDisappeared()) status += 16;
834   if(TVirtualMC::GetMC()->IsTrackStop())        status += 32;
835   if(TVirtualMC::GetMC()->IsTrackAlive())       status += 64;
836   //
837   // retrieve the indices with the volume path
838   //
839   TVirtualMC::GetMC()->TrackPosition(position);
840   int chip=-1,module=-1,sstave=-1,stave=-1,level=0; // volume copies on different levels
841   TVirtualMC::GetMC()->CurrentVolOffID(++level,chip);
842   if (fGeomTGeo->GetNModules(lay)>0)    TVirtualMC::GetMC()->CurrentVolOffID(++level,module);
843   if (fGeomTGeo->GetNHalfStaves(lay)>0) TVirtualMC::GetMC()->CurrentVolOffID(++level,sstave);
844   TVirtualMC::GetMC()->CurrentVolOffID(++level,stave);
845   //
846   chipID = fGeomTGeo->GetChipIndex(lay,stave,sstave,module,chip);
847   // Fill hit structure.
848   //
849   hit.SetChip(chipID);
850   hit.SetTrack(gAlice->GetMCApp()->GetCurrentTrackNumber());
851   TVirtualMC::GetMC()->TrackPosition(position);
852   TVirtualMC::GetMC()->TrackMomentum(momentum);
853   hit.SetPosition(position);
854   hit.SetTime(TVirtualMC::GetMC()->TrackTime());
855   hit.SetMomentum(momentum);
856   hit.SetStatus(status);
857   hit.SetEdep(TVirtualMC::GetMC()->Edep());
858   hit.SetShunt(GetIshunt());
859   if(TVirtualMC::GetMC()->IsTrackEntering()){
860     hit.SetStartPosition(position);
861     hit.SetStartTime(TVirtualMC::GetMC()->TrackTime());
862     hit.SetStartStatus(status);
863     return; // don't save entering hit.
864   } // end if IsEntering
865     // Fill hit structure with this new hit.
866     //Info("StepManager","Calling Copy Constructor");
867   new(lhits[fNhits++]) AliITSUHit(hit); // Use Copy Construtor.
868   // Save old position... for next hit.
869   hit.SetStartPosition(position);
870   hit.SetStartTime(TVirtualMC::GetMC()->TrackTime());
871   hit.SetStartStatus(status);
872
873   return;
874 }
875
876 //______________________________________________________________________
877 void AliITSUv1::SetLayerChipTypeID(Int_t lr, UInt_t id)
878 {
879   // set det type
880   if (!fChipTypeID || fNLayers<=lr) AliFatal(Form("Number of layers %d, %d is manipulated",fNLayers,lr));
881   fChipTypeID[lr] = id;
882 }
883
884 //______________________________________________________________________
885 Int_t AliITSUv1::GetLayerChipTypeID(Int_t lr)
886 {
887   // set det type
888   if (!fChipTypeID || fNLayers<=lr) AliFatal(Form("Number of layers %d, %d is manipulated",fNLayers,lr));
889   return fChipTypeID[lr];
890 }