]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/UPGRADE/v0/AliITSupgrade.cxx
Merge branch 'master' into dev
[u/mrichter/AliRoot.git] / ITS / UPGRADE / v0 / AliITSupgrade.cxx
CommitLineData
1d9af2d5 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
c93255fe 18#include <TTree.h>
1d9af2d5 19#include <TArrayD.h> //new constructor
20#include <TFile.h>
0ac80088 21#include <TGeoManager.h>
1d9af2d5 22#include <TGeoVolume.h> //CreateGeometry()
0ac80088 23#include <TGeoMatrix.h>
1d9af2d5 24#include <TVirtualMC.h> //->gMC in StepManager
25#include <TPDGCode.h> //StepHistory
26#include <TClonesArray.h>
2accaff0 27#include <TGeoGlobalMagField.h>
1d9af2d5 28#include "AliRun.h" //CreateMaterials()
29#include "AliMC.h" //CreateMaterials()
30#include "AliMagF.h" //CreateMaterials()
1d9af2d5 31#include "AliITSupgrade.h" //class header
32#include "AliITShit.h" // StepManager()
33#include "AliITSDigitUpgrade.h"
34#include "AliTrackReference.h" // StepManager()
a19b4e8c 35#include "AliITSsegmentationUpgrade.h"
36
1d9af2d5 37ClassImp(AliITSupgrade)
38//__________________________________________________________________________________________________
39 AliITSupgrade::AliITSupgrade():
40 AliITS(),
41 fWidths(0),
42 fRadii(0),
0ac80088 43 fNSectors(20),
1d9af2d5 44 fRadiiCu(0),
45 fWidthsCu(0),
46 fCopper(0),
47 fBeampipe(0),
48 fRadiusBP(0),
49 fWidthBP(0),
50 fHalfLengthBP(0),
51 fNlayers(0),
52 fHalfLength(0),
53 fSdigits(0),
a19b4e8c 54 fDigitArray(0),
1d9af2d5 55 fSegmentation(0x0)
56{
57 //
58 // default ctor
387f1a9e 59 //
60}
61
62//__________________________________________________________________________________________________
63AliITSupgrade::AliITSupgrade(const char *name,const char *title, Bool_t isBeamPipe):
64 AliITS(name,title),
65 fWidths(0),
66 fRadii(0),
0ac80088 67 fNSectors(20),
387f1a9e 68 fRadiiCu(0),
69 fWidthsCu(0),
70 fCopper(0),
71 fBeampipe(isBeamPipe),
72 fRadiusBP(0),
73 fWidthBP(0),
74 fHalfLengthBP(0),
75 fNlayers(0),
76 fHalfLength(0),
77 fSdigits(0),
a19b4e8c 78 fDigitArray(0),
387f1a9e 79 fSegmentation(0x0)
80{
a19b4e8c 81 // Default values are used in order to simulate the standard ITS material budget (see The ALICE Collaboration et al 2008 JINST 3 S08002 - Figure 3.2).
82 // The cell segmentation is chosen to achieve the tracking resolution as described in Table 3.2 of The ALICE Collaboration et al 2008 JINST 3 S08002,
83 // apart from SPD0 where the 9 um value is considered : resolution*sqrt(12)
84 // Cilinder lenghts are chosen according to Table 1 of the Alignment paper : http://iopscience.iop.org/1748-0221/5/03/P03003
85
387f1a9e 86 fNlayers = 6;
87 fWidths.Set(fNlayers);
88 fRadii.Set(fNlayers);
89 fRadiiCu.Set(fNlayers);
90 fWidthsCu.Set(fNlayers);
91 fCopper.Set(fNlayers);
92 fHalfLength.Set(fNlayers);
93
387f1a9e 94
a19b4e8c 95 Double_t xsz[6]={31.18*1e-04,41.6*1e-04,121.2e-04,121.4*1e-04,69.3*1e-04,69.3*1e-04};
387f1a9e 96 Double_t zsz[6]={416*1e-04,416*1e-04,97*1e-04,97*1e-04,2875*1e-04,2875*1e-04};
a19b4e8c 97 Double_t halfL[6]={14.1,14.1,22.2,29.7,43.1,48.9};
387f1a9e 98 Double_t r[6]={4.,7.6,14.9,23.8,39.1,43.6};
99 Double_t thick[6]={75.*1e-04,150.*1e-04,150.*1e-04,150.*1e-04,150.*1e-04,150.*1e-04};
100
101 Int_t nlayers =6;
102 TArrayD xsizeSi(nlayers);
103 TArrayD zsizeSi(nlayers);
104
105 // adding beam pipe upgrade
106 if(isBeamPipe){
9b615954 107 fRadiusBP = 2.0;
387f1a9e 108 fWidthBP = 0.08;
9b615954 109 fHalfLengthBP = 400;
387f1a9e 110 }
111
112 // setting the geonetry parameters and the segmentation
113 Int_t npixHalf[6], npixR[6];
a19b4e8c 114 Double_t halfLmod[6], xszInt[6];
387f1a9e 115 Int_t c[6]={1,1,1,1,1,1};
116
117 for(Int_t i=0; i<nlayers; i++){
118 // recalc length
a19b4e8c 119 npixHalf[i]=(Int_t)(halfL[i]/zsz[i]);
120 halfLmod[i]=npixHalf[i]*zsz[i];
387f1a9e 121 // recalc segmentation
122 npixR[i] = (Int_t)(2*TMath::Pi()*r[i]/xsz[i]);
123 xszInt[i]= 2*TMath::Pi()*r[i]/npixR[i];
124 xsizeSi.AddAt(xszInt[i],i);
125 zsizeSi.AddAt(zsz[i],i);
126
a19b4e8c 127 fHalfLength.AddAt(halfLmod[i],i);
387f1a9e 128
129 fRadii.AddAt(r[i],i);
130 fWidths.AddAt(thick[i],i);
131
132 fRadiiCu.AddAt(r[i]+thick[i],i);
133 fWidthsCu.AddAt(0.015,i);//micron
134 fCopper.AddAt(c[i],i);
135 }
136
137 SetFullSegmentation(xsizeSi,zsizeSi);
1d9af2d5 138}
139//__________________________________________________________________________________________________
140AliITSupgrade::AliITSupgrade(const char *name,const char *title, TArrayD widths, TArrayD radii,TArrayD halfLengths, TArrayD radiiCu, TArrayD widthsCu, TArrayS copper,Bool_t bp,Double_t radiusBP, Double_t widthBP, Double_t halfLengthBP):
141 AliITS(name,title),
142 fWidths(0),
143 fRadii(0),
0ac80088 144 fNSectors(20),
1d9af2d5 145 fRadiiCu(radiiCu),
146 fWidthsCu(widthsCu),
147 fCopper(copper),
148 fBeampipe(bp),
387f1a9e 149 fRadiusBP(0),
150 fWidthBP(0),
151 fHalfLengthBP(0),
1d9af2d5 152 fNlayers(widths.GetSize()),
153 fHalfLength(halfLengths),
154 fSdigits(0),
a19b4e8c 155 fDigitArray(0),
1d9af2d5 156 fSegmentation(0x0)
157{
1d9af2d5 158 //
a19b4e8c 159 // constructor to set all the geometrical/segmentation specs from outside
160 //
1d9af2d5 161 for(Int_t i=0;i<fNlayers;i++){
162 fWidths.Set(i+1);
163 fWidths.AddAt(widths.At(i),i);
164 fRadii.Set(i+1);
165 fRadii.AddAt(radii.At(i),i);
166 AliDebug(1,"Creating Volume");
167 }
387f1a9e 168
169 if(bp){
170 fRadiusBP=radiusBP;
171 fWidthBP=widthBP;
172 fHalfLengthBP=halfLengthBP;
173 }
1d9af2d5 174 Init();
175
176}
177
178AliITSupgrade::~AliITSupgrade(){
179
180 if(fSdigits) {fSdigits->Delete(); delete fSdigits;}
a19b4e8c 181 if(fDigitArray) {fDigitArray->Delete(); delete fDigitArray;}
1d9af2d5 182}
183
184
185
186//_________________________________________________________________________________________________
187void AliITSupgrade::AddAlignableVolumes()const
188{
189 //not needed
190 return;
191}
192
193//__________________________________________________________________________________________________
194
195void AliITSupgrade::CreateMaterials()
196{
197 //
198 // Definition of ITS materials
199 //
200
201 AliInfo("Start ITS materials");
202 //data from PDG booklet 2002 density [gr/cm^3] rad len [cm] abs len [cm]
203 Float_t aAir[4]={12,14,16,36} , zAir[4]={6,7,8,18} , wAir[4]={0.000124,0.755267,0.231781,0.012827} , dAir=0.00120479; Int_t nAir=4;//mixture
204 Float_t aBe = 9.012 , zBe = 4 , dBe = 1.848 , radBe = 65.19/dBe , absBe = 77.8/dBe ; // UPGRADE -> beryllium beampipe
205 Float_t aSi = 28.085 , zSi = 14 , dSi = 2.329 , radSi = 21.82/dSi , absSi = 108.4/dSi ; // UPGRADE -> Si tube
206 Float_t aCu = 63.54, zCu = 29 , dCu= 8.96 , radCu = 12.86/dCu, absCu= 137.3/dCu;//Upgrade -> Copper Tube
207
208 Int_t matId=0; //tmp material id number
209 Int_t unsens = 0, sens=1; //sensitive or unsensitive medium
210 Int_t itgfld = 3; //type of field intergration 0 no field -1 user in guswim 1 Runge Kutta 2 helix 3 const field along z
211 Float_t maxfld = 5.; //max field value
212 Float_t tmaxfd = -10.0; //max deflection angle due to magnetic field in one step
213 Float_t deemax = - 0.2; //max fractional energy loss in one step
214 Float_t stemax = - 0.1; //max step allowed [cm]
a346a81b 215 Float_t epsil = 0.001; //abs tracking precision [cm]
1d9af2d5 216 Float_t stmin = - 0.001; //min step size [cm] in continius process transport, negative value: choose it automatically
217 Float_t tmaxfdSi = 0.1; // .10000E+01; // Degree
218 Float_t stemaxSi = 0.0075; // .10000E+01; // cm
219 Float_t deemaxSi = 0.1; // 0.30000E-02; // Fraction of particle's energy 0<deemax<=1
220 Float_t epsilSi = 1.0E-4;// .10000E+01;
221 Float_t stminSi = 0.0; // cm "Default value used"
222
223 Float_t epsilBe = .001; // Tracking precision,
224 Float_t stemaxBe = -0.01; // Maximum displacement for multiple scat
225 Float_t tmaxfdBe = -20.; // Maximum angle due to field deflection
226 Float_t deemaxBe = -.3; // Maximum fractional energy loss, DLS
227 Float_t stminBe = -.8;
2accaff0 228 Int_t isxfld = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Integ();//from CreateMaterials in STRUCT/AliPIPEv3.cxx
229 Float_t sxmgmx = ((AliMagF*)TGeoGlobalMagField::Instance()->GetField())->Max();//from CreateMaterials in STRUCT/AliPIPEv3.cxx
1d9af2d5 230
231
232 AliMixture(++matId,"UpgradeAir" ,aAir ,zAir ,dAir ,nAir ,wAir );
233 AliMedium(kAir ,"UpgradeAir" ,matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
234
235 AliMaterial(++matId,"UpgradeSi" ,aSi ,zSi ,dSi ,radSi ,absSi );
a346a81b 236 AliMedium(kSi ,"UpgradeSi" , matId, sens, isxfld, sxmgmx, tmaxfdSi, stemaxSi, deemaxSi, epsilSi, stminSi);
1d9af2d5 237
238 AliMaterial(++matId,"UpgradeBe" ,aBe ,zBe ,dBe ,radBe ,absBe );
2accaff0 239 AliMedium(kBe ,"UpgradeBe" , matId, unsens, isxfld, sxmgmx, tmaxfdBe, stemaxBe, deemaxBe, epsilBe, stminBe);
1d9af2d5 240
241 AliMaterial(++matId, "UpgradeCu", aCu, zCu, dCu, radCu, absCu);
242 AliMedium(kCu, "UpgradeCu", matId, unsens, itgfld, maxfld, tmaxfd, stemax, deemax, epsil, stmin);
243
244 AliInfo("End ITS materials");
245
246}//void AliITS::CreateMaterials()
247//__________________________________________________________________________________________________
248void AliITSupgrade::CreateGeometry()
249{
250 //
251 //Creates detailed geometry simulation (currently GEANT volumes tree)
252 //
253 AliInfo("Start ITS upgrade preliminary version building");
254 if(!gMC->IsRootGeometrySupported()) return;
255 TGeoVolumeAssembly *vol= CreateVol();
256 gGeoManager->GetVolume("ALIC")->AddNode(vol,0);
257 AliInfo("Stop ITS upgrade preliminary version building");
387f1a9e 258 PrintSummary();
1d9af2d5 259}
260//__________________________________________________________________________________________________
1d9af2d5 261void AliITSupgrade::StepManager()
262{
263 // Full Step Manager.
264 // Arguments: none
265 // Returns: none
a346a81b 266 // StepHistory(); return; //uncomment to print tracks history
1d9af2d5 267 // StepCount(); return; //uncomment to count photons
0ac80088 268
269 if(!fSegmentation) AliFatal("No segmentation available");
270
1d9af2d5 271 if(!(this->IsActive())) return;
272 if(!(gMC->TrackCharge())) return;
273 TString volumeName=gMC->CurrentVolName();
0ac80088 274 if(volumeName.Contains("Be")) return;
275 if(volumeName.Contains("Cu")) return;
276 if(gMC->IsTrackExiting()) {
1d9af2d5 277 AddTrackReference(gAlice->GetMCApp()->GetCurrentTrackNumber(), AliTrackReference::kITS);
278 } // if Outer ITS mother Volume
279
280 static TLorentzVector position, momentum; // Saves on calls to construtors
281 static AliITShit hit;// Saves on calls to constructors
282
283
284 Int_t status = 0;
285
286 //
287 // Track status
288 if(gMC->IsTrackInside()) status += 1;
289 if(gMC->IsTrackEntering()) status += 2;
290 if(gMC->IsTrackExiting()) status += 4;
291 if(gMC->IsTrackOut()) status += 8;
292 if(gMC->IsTrackDisappeared()) status += 16;
293 if(gMC->IsTrackStop()) status += 32;
294 if(gMC->IsTrackAlive()) status += 64;
295
296 //
297 // Fill hit structure.
298 //
0ac80088 299 Int_t copy=-1;
300 gMC->CurrentVolID(copy);
301
302 volumeName.Remove(0,12); // remove letters to get the layer number
303 hit.SetModule(fSegmentation->GetIdIndex(volumeName.Atoi(),copy)); // layer and sector information are together in the IdIndex (if copy=0 the idIndex is the layer));
1d9af2d5 304 hit.SetTrack(gAlice->GetMCApp()->GetCurrentTrackNumber());
305
306 gMC->TrackPosition(position);
307 gMC->TrackMomentum(momentum);
308 hit.SetPosition(position);
309
310 hit.SetTime(gMC->TrackTime());
311 hit.SetMomentum(momentum);
312 hit.SetStatus(status);
313 hit.SetEdep(gMC->Edep());
314 hit.SetShunt(GetIshunt());
315 if(gMC->IsTrackEntering()){
316 hit.SetStartPosition(position);
317 hit.SetStartTime(gMC->TrackTime());
318 hit.SetStartStatus(status);
319 return; // don't save entering hit.
320 }
387f1a9e 321 // Fill hit structure with this new hit.
1d9af2d5 322
323 new((*fHits)[fNhits++]) AliITShit(hit); // Use Copy Construtor.
324 // Save old position... for next hit.
325 hit.SetStartPosition(position);
326 hit.SetStartTime(gMC->TrackTime());
327 hit.SetStartStatus(status);
328 return;
329
330}
331//__________________________________________________________________________________________________
332TGeoVolumeAssembly * AliITSupgrade::CreateVol()
333{
a19b4e8c 334//
335// method to create the Upgrade Geometry (silicon, copper cylinders and beampipe)
336//
337 TGeoVolumeAssembly *vol = new TGeoVolumeAssembly("ITSupgrade");
1d9af2d5 338 TGeoMedium *si =gGeoManager->GetMedium("ITS_UpgradeSi");
339 TGeoMedium *cu =gGeoManager->GetMedium("ITS_UpgradeCu");
340 TGeoMedium *be =gGeoManager->GetMedium("ITS_UpgradeBe");
341 for(Int_t ivol=0;ivol<fNlayers;ivol++){
1d9af2d5 342
0ac80088 343 if(fNSectors<1){
344 TGeoVolume *layer=gGeoManager->MakeTube(Form("LayerSilicon%i",ivol),si , fRadii.At(ivol) , fRadii.At(ivol)+fWidths.At(ivol) , fHalfLength.At(ivol)); //upgraded situation
345
346 TGeoVolume *layerCu=gGeoManager->MakeTube(Form("LayerCu%i",ivol),cu , fRadiiCu.At(ivol) , fRadiiCu.At(ivol)+fWidthsCu.At(ivol) , fHalfLength.At(ivol) ); //upgraded situation
347
348 vol ->AddNode(layer,ivol);
349 if(fCopper.At(ivol)){
350 vol->AddNode(layerCu,ivol);
351 }
0ac80088 352 }else{
353
354 TGeoVolume *layer = gGeoManager->MakeTubs(Form("LayerSilicon%i",ivol),si, fRadii.At(ivol), fRadii.At(ivol)+fWidths.At(ivol) ,fHalfLength.At(ivol),0,(360./fNSectors));
355 TGeoVolume *layerCu = gGeoManager->MakeTubs(Form("LayerCu%i",ivol),cu , fRadiiCu.At(ivol) , fRadiiCu.At(ivol)+fWidthsCu.At(ivol) , fHalfLength.At(ivol) ,0,(360./fNSectors));
356
357
358 for(Int_t i=0;i<fNSectors;i++){
359 TGeoRotation *rot1 = new TGeoRotation(" ",0.0,0.0,360.*i/fNSectors);//sector rotation
360 TGeoRotation *rot2 = new TGeoRotation(" ",0.0,0.0,360.*i/fNSectors);//
361 vol->AddNode(layer,i,rot1);
362 if(fCopper.At(ivol)){
363 vol->AddNode(layerCu,i,rot2);
364 }
365 }
1d9af2d5 366 }
367 }
0ac80088 368
9b615954 369
0ac80088 370 if(fBeampipe) {
371 TGeoVolume *beampipe=gGeoManager->MakeTube("BeamPipe", be , fRadiusBP , fRadiusBP+ fWidthBP , fHalfLengthBP ); //upgraded situation
9b615954 372 vol->AddNode(beampipe,fNlayers);
0ac80088 373 }
374 return vol;
1d9af2d5 375}
376//_________________________________________________________________________________________________
377void AliITSupgrade::SetFullSegmentation(TArrayD xsize,TArrayD zsize){
a19b4e8c 378//
379// Upgrade detector specs specified in used constructor are stored in a file .
380// for further usage (loc<->master<->tracking reference system changes)
0ac80088 381
a19b4e8c 382 Bool_t check=kFALSE;
0ac80088 383 for(Int_t lay = 0; lay< xsize.GetSize(); lay++){
384 Double_t arch = fRadii.At(lay)*(TMath::Pi()*2/fNSectors);
385 Int_t nPixRPhi = (Int_t)(arch/xsize.At(lay));
386 Int_t nPixZed = (Int_t)((2*fHalfLength.At(lay))/zsize.At(lay));
a19b4e8c 387 if(nPixRPhi>9999)check=kTRUE;
388 if(nPixZed>99999)check=kTRUE;
0ac80088 389 }
a19b4e8c 390 if(check) AliFatal(" Segmentation is too small!! ");
0ac80088 391 if(fSegmentation) fSegmentation->SetNSectors(fNSectors);
392 TArrayD nSect(1);
393 nSect.AddAt(fNSectors,0);
1d9af2d5 394 TFile *file= TFile::Open("Segmentation.root","recreate");
395 file->WriteObjectAny(&xsize,"TArrayD","CellSizeX");
0ac80088 396 file->WriteObjectAny(&zsize,"TArrayD","CellSizeZ");
397 file->WriteObjectAny(&nSect,"TArrayD","nSectors");
1d9af2d5 398 file->Close();
399}
400//_________________________________________________________________________________________________
401void AliITSupgrade::StepHistory()
9b615954 402{
1d9af2d5 403 // This methode is invoked from StepManager() in order to print out
a346a81b 404 TString volumeName=gMC->CurrentVolName();
405 if(!volumeName.Contains("Silicon")) return;
1d9af2d5 406 static Int_t iStepN;
407 const char *sParticle;
408 switch(gMC->TrackPid()){
409 case kProton: sParticle="PROTON" ;break;
410 case kNeutron: sParticle="neutron" ;break;
411 case kGamma: sParticle="gamma" ;break;
412 case kPi0: sParticle="Pi0" ;break;
413 case kPiPlus: sParticle="Pi+" ;break;
414 case kPiMinus: sParticle="Pi-" ;break;
415 case kElectron: sParticle="electron" ;break;
416 default: sParticle="not known" ;break;
417 }
418
9b615954 419 TString flag="funny combination";
420 if(gMC->IsTrackAlive()) {
1d9af2d5 421 if(gMC->IsTrackEntering()) flag="enters to";
422 else if(gMC->IsTrackExiting()) flag="exits from";
423 else if(gMC->IsTrackInside()) flag="inside";
424 else if(gMC->IsTrackStop()) flag="stopped in";
9b615954 425 }
426
1d9af2d5 427 Int_t vid=0,copy=0;
428 TString path=gMC->CurrentVolName(); path.Prepend("-");path.Prepend(gMC->CurrentVolOffName(1));//current volume and his mother are always there
429 vid=gMC->CurrentVolOffID(2,copy); if(vid) {path.Prepend("-");path.Prepend(gMC->VolName(vid));}
430 vid=gMC->CurrentVolOffID(3,copy); if(vid) {path.Prepend("-");path.Prepend(gMC->VolName(vid));}
a346a81b 431
1d9af2d5 432
a346a81b 433 AliInfo(Form("\n Step %i: %s (%i) %s %s m=%.6f GeV q=%.1f dEdX=%.4f Etot=%.4f",iStepN,sParticle,gMC->TrackPid(),flag.Data(),path.Data(),gMC->TrackMass(),gMC->TrackCharge(),gMC->Edep()*1e9,gMC->Etot()));
1d9af2d5 434
435 Double_t gMcTrackPos[3]; gMC->TrackPosition(gMcTrackPos[0],gMcTrackPos[1],gMcTrackPos[2]);
436 Double_t gMcTrackPosLoc[3]; gMC->Gmtod(gMcTrackPos,gMcTrackPosLoc,1);
a346a81b 437 TString v(volumeName.Data());
438 v.ReplaceAll("LayerSilicon","");
439 Int_t ilayer = v.Atoi();
440 Double_t rXY = TMath::Sqrt(gMcTrackPos[0]*gMcTrackPos[0]+gMcTrackPos[1]*gMcTrackPos[1]);
441 AliInfo(Form("gMC Track Position (MARS) x: %5.3lf, y: %5.3lf, z: %5.3lf (r: %5.3lf) (deltaR %5.5lf - width %5.5f)",gMcTrackPos[0],gMcTrackPos[1],gMcTrackPos[2], rXY , rXY - (fRadii.At(ilayer)),fWidths.At(ilayer)));
1d9af2d5 442
443
444 AliDebug(10,Form("Step %i: tid=%i flags alive=%i disap=%i enter=%i exit=%i inside=%i out=%i stop=%i new=%i",
445 iStepN, gAlice->GetMCApp()->GetCurrentTrackNumber(),
446 gMC->IsTrackAlive(), gMC->IsTrackDisappeared(),gMC->IsTrackEntering(), gMC->IsTrackExiting(),
447 gMC->IsTrackInside(),gMC->IsTrackOut(), gMC->IsTrackStop(), gMC->IsNewTrack()));
448
449 Float_t a,z,den,rad,abs; a=z=den=rad=abs=-1;
450 Int_t mid=gMC->CurrentMaterial(a,z,den,rad,abs);
451 AliDebug(10, Form("Step %i: mid=%i a=%7.2f z=%7.2f den=%9.4f rad=%9.2f abs=%9.2f\n\n",iStepN,mid,a,z,den,rad,abs));
6cb8ce4c 452/*
1d9af2d5 453 TArrayI proc; gMC->StepProcesses(proc);
6cb8ce4c 454
a346a81b 455 AliInfo("Processes in this step:");
1d9af2d5 456 for ( int i = 0 ; i < proc.GetSize(); i++)
457 {
a346a81b 458 printf(Form(" - %s - ",TMCProcessName[proc.At(i)]));
1d9af2d5 459 }
a346a81b 460 printf("\n");
461 AliInfo("End process list");
6cb8ce4c 462*/
1d9af2d5 463 iStepN++;
464}//StepHistory()
465//______________________________________________________
466void AliITSupgrade::Hits2SDigits(){
467
9b615954 468 // Interface method invoked from AliSimulation to create a list of sdigits corresponding to list of hits.
469 // Every hit generates one or more sdigits.
1d9af2d5 470 // Arguments: none
471 // Returns: none
472 AliDebug(1,"Start Hits2SDigits.");
473
474 if(!fSegmentation)fSegmentation=new AliITSsegmentationUpgrade();
475
476 if(!GetLoader()->TreeH()) {GetLoader()->LoadHits();}
477 if(!GetLoader()->TreeS())
478
479 for(Int_t iEvt=0;iEvt < GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvt++){//events loop
480 GetLoader()->GetRunLoader()->GetEvent(iEvt); //get next event
481 {
482 GetLoader()->MakeTree("S");
483 MakeBranch("S");
484 }
485 SetTreeAddress();
486 Int_t nSdigit[10]; for(Int_t i=0;i<10;i++) nSdigit[i] =0;
9b615954 487 for(Int_t iEnt=0;iEnt<GetLoader()->TreeH()->GetEntries();iEnt++){ // prims loop
488 GetLoader()->TreeH()->GetEntry(iEnt);
1d9af2d5 489 Hit2SumDig(Hits(),SDigitsList(),nSdigit);//convert this hit to list of sdigits
490 }//prims loop
491 GetLoader()->TreeS()->Fill();
492 GetLoader()->WriteSDigits("OVERWRITE");
493 SDigitsReset();
494 }//events loop
495 GetLoader()->UnloadHits();
496 GetLoader()->UnloadSDigits();
497 AliDebug(1,"Stop Hits2SDigits.");
498
499}
500//____________________________
a19b4e8c 501void AliITSupgrade::Hit2SumDig(TClonesArray *hits,const TObjArray *pSDig, Int_t *nSdigit)
1d9af2d5 502{
503 // Adds sdigits of this hit to the list
504 // Returns: none
505
506 AliDebug(1,"Start Hit2SumDig");
9b615954 507
1d9af2d5 508 if(!fSegmentation){
509 AliDebug(1,"Segmentation Not inizialized!!");
510 return ;
511 }
9b615954 512
3303c513 513 TClonesArray *pSdigList[8]; // is is the max number of layers allowed in the tracking
514 for(Int_t il=0; il<8; il++) pSdigList[il]=0x0;
9b615954 515
1d9af2d5 516 for(Int_t i=0;i<fNlayers;i++){
517 pSdigList[i]=(TClonesArray*)(*pSDig)[i];
9b615954 518 if(pSdigList[i]->GetEntries()!=0 && nSdigit[i]==0) {
519 AliDebug(1,Form("Entries of pSdigList %d; layer: %d,",pSdigList[i]->GetEntries(),i));
520 AliErrorClass(" -> Some of sdigits lists is not empty"); //in principle those lists should be empty
521 }
1d9af2d5 522 }
523
1d9af2d5 524 for(Int_t iHit=0;iHit<hits->GetEntries();iHit++){ //hits loop
525 AliITShit *hit = (AliITShit*)hits->At(iHit);
526 Double_t xz[2];
0ac80088 527
14b8503a 528 Int_t module=99;
0ac80088 529 if(!fSegmentation->GlobalToDet(fSegmentation->GetLayerFromIdIndex(hit->GetModule()),hit->GetXG(),hit->GetYG(),hit->GetZG(),xz[0],xz[1],module)) continue;
1d9af2d5 530 AliITSDigitUpgrade digit;
531 digit.SetSignal(hit->GetIonization());
532 digit.SetNelectrons(hit->GetIonization()/(3.62*1e-09));
0ac80088 533 digit.SetLayer(fSegmentation->GetLayerFromIdIndex(hit->GetModule()));
534 digit.SetModule(fSegmentation->GetModuleFromIdIndex(hit->GetModule()));//set the module (=sector) of ITSupgrade
9b615954 535
1d9af2d5 536 digit.SetTrackID(hit->GetTrack());
537
538 Int_t xpix = 999;
1d9af2d5 539 Int_t zpix = 999; // shift at the next line to have zpix always positive. Such numbers are used to build the Pixel Id in the layer (> 0!)
0ac80088 540 fSegmentation->DetToPixID(xz[0], xz[1],fSegmentation->GetLayerFromIdIndex(hit->GetModule()), xpix, zpix);
1d9af2d5 541 digit.SetPixId(xpix,zpix);
0ac80088 542 new((*pSdigList[fSegmentation->GetLayerFromIdIndex(hit->GetModule())])[nSdigit[fSegmentation->GetLayerFromIdIndex(hit->GetModule())]++]) AliITSDigitUpgrade(digit);
9b615954 543
1d9af2d5 544 }
545
546 AliDebug(1,"Stop Hit2SumDig.");
547}
548//_______________________________________________________________________________________________
549void AliITSupgrade::MakeBranch(Option_t *option){
1d9af2d5 550 //Create Tree branches
551 AliDebug(1,Form("Start with option= %s.",option));
552
553 const Int_t kBufSize = 4000;
554
555 const char *cH = strstr(option,"H");
556 const char *cD = strstr(option,"D");
6ea6235b 557 //const char *cR = strstr(option,"R");
1d9af2d5 558 const char *cS = strstr(option,"S");
559
560 if(cH&&fLoader->TreeH()){
561
562 HitCreate();
563 MakeBranchInTree(fLoader->TreeH(),"ITSupgrade",&fHits,kBufSize,0);
564 }
565
566
567 if(cS&&fLoader->TreeS()){
568 SDigitsCreate();
569 for(Int_t i=0;i<fNlayers;i++) MakeBranchInTree(fLoader->TreeS(),Form("Layer%d",i),&((*fSdigits)[i]),kBufSize,0);
570 }
571
572
573 if(cD&&fLoader->TreeD()){
574 DigitsCreate();
a19b4e8c 575 for(Int_t i=0;i<fNlayers;i++) MakeBranchInTree(fLoader->TreeD(),Form("Layer%d",i),&((*fDigitArray)[i]),kBufSize,0);
1d9af2d5 576 }
1d9af2d5 577 AliDebug(1,"Stop.");
578}
579//____________________________________________________________________________________________________
580void AliITSupgrade::SetTreeAddress()
581{
582 //Set branch address for the Hits and Digits Tree.
583 AliDebug(1,"Start.");
584 if(fLoader->TreeH() && fLoader->TreeH()->GetBranch("ITSupgrade" )){
585 HitCreate();
586 fLoader->TreeH()->SetBranchAddress("ITSupgrade",&fHits);
587
588 }
1d9af2d5 589
590 if(fLoader->TreeS() && fLoader->TreeS()->GetBranch("Layer0" )){
591 SDigitsCreate();
592 for(int i=0;i<fNlayers;i++){
593 fLoader->TreeS()->SetBranchAddress(Form("Layer%d",i),&((*fSdigits)[i]));
594 }
595 }
1d9af2d5 596
597 if(fLoader->TreeD() && fLoader->TreeD()->GetBranch("Layer0")){
598 DigitsCreate();
599 for(int i=0;i<fNlayers;i++){
a19b4e8c 600 fLoader->TreeD()->SetBranchAddress(Form("Layer%d",i),&((*fDigitArray)[i]));
1d9af2d5 601 }
602 }
1d9af2d5 603 AliDebug(1,"Stop.");
604}
387f1a9e 605//____________________________________________________________________________________________________
606void AliITSupgrade::PrintSummary()
607{
608
609 // pdg booklet or http://pdg.lbl.gov/2010/reviews/rpp2010-rev-atomic-nuclear-prop.pdf
610
611 // X0 as X0[g/cm^2]/density
612 Double_t beX0= 65.19/1.848;
613 Double_t cuX0= 12.86/8.96;
614 Double_t siX0= 21.82/2.329;
615
616 if(!fSegmentation) fSegmentation= new AliITSsegmentationUpgrade();
617 printf(" %10s %10s %10s %10s %10s %10s %10s %10s\n","Name","R [cm]","x/X0 [%]","Phi res[um]","Z res[um]","X segm[um]","Y segm[um]","widths [um]");
ebaa7604 618
619 if(fBeampipe) printf(" Beampipe %10.3f %10f %10s %10s %10s %10s %10.1f\n", fRadiusBP, 100.*fWidthBP/beX0,"-","-","-","-",fWidthBP*1.e+4);
620 else printf(" * No Beam Pipe Upgrade \n");
387f1a9e 621 for(Int_t i=0; i< fNlayers; i++){
622 printf(" Si Layer %1i %10.3f %10f %10.2f %10.2f %10.2f %10.2f %10.1f\n",i, fRadii.At(i), 100.*fWidths.At(i)/siX0, (fSegmentation->GetCellSizeX(i) *1.e+4)/TMath::Sqrt(12.), (fSegmentation->GetCellSizeZ(i)*1.e+4)/TMath::Sqrt(12.),fSegmentation->GetCellSizeX(i) *1.e+4, fSegmentation->GetCellSizeZ(i)*1.e+4,fWidths.At(i)*1.e+4);
623 printf(" Cu Layer %10.3f %10f %10s %10s %10s %10s %10.1f\n", fRadiiCu.At(i), 100.*fWidthsCu.At(i)/cuX0,"-","-","-","-",fWidthsCu.At(i)*1.e+4);
624 }
625}
626//____________________________________________________________________________________________________
627void AliITSupgrade::SetSegmentationX(Double_t x, Int_t lay){
a19b4e8c 628//
629// Virtual dimension in X direction is written into the file
630// (-> see SetFullSegmentation)
631
387f1a9e 632TFile *f = TFile::Open("Segmentation.root","UPDATE");
633if(!f) {
634 AliError("\n\n\n Segmentation.root file does not exist. The segmentation is 0! \n\n\n");
635 return;
636 }
637else {
638 TArrayD *a = (TArrayD*)f->Get("CellSizeX");
639 if(!a) {
640 AliError("CessSizeX array does not exist!");
641 return;
642 }
643 a->AddAt(x,lay);
644 f->WriteObjectAny(a,"TArrayD","CellSizeX");
645 f->Close();
646 delete f;
647 }
648}
649//____________________________________________________________________________________________________
650void AliITSupgrade::SetSegmentationZ(Double_t z, Int_t lay){
a19b4e8c 651//
652// Virtual dimension in X direction is written into the file
653// (-> see SetFullSegmentation)
654//
655
387f1a9e 656
657TFile *f = TFile::Open("Segmentation.root","UPDATE");
658if(!f) {
659 AliError("\n\n\n Segmentation.root file does not exist. The segmentation is 0! \n\n\n");
660 return;
661 }
662else {
663 TArrayD *a = (TArrayD*)f->Get("CellSizeZ");
664 if(!a) {
665 AliError("CellSizeZ array does not exist!");
666 return;
667 }
668 a->AddAt(z,lay);
669 f->WriteObjectAny(a,"TArrayD","CellSizeZ");
670 f->Close();
671 delete f;
672 }
673
674}