Divide by zero error fix in function F. Sigma difusion = 0 now set to same
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSSD.cxx
CommitLineData
b0f5e3fc 1
2#include <stdio.h>
3#include <TObjArray.h>
4
b0f5e3fc 5#include "AliITSsimulationSSD.h"
6#include "AliITSdictSSD.h"
7#include "AliITSdcsSSD.h"
8#include "AliITS.h"
9#include "AliRun.h"
10
11
12ClassImp(AliITSsimulationSSD);
13//------------------------------------------------------------
14AliITSsimulationSSD::AliITSsimulationSSD(AliITSsegmentation *seg,
15 AliITSresponse *resp){
16 // Constructor
17
18 fSegmentation = seg;
19 fResponse = resp;
20 fDCS = new AliITSdcsSSD(seg,resp);
21
22 fNstrips = fSegmentation->Npx();
23 fPitch = fSegmentation->Dpx(0);
24
e8189707 25 fP = new TArrayF(fNstrips+1);
26 fN = new TArrayF(fNstrips+1);
b0f5e3fc 27
e8189707 28 fTracksP = new AliITSdictSSD[fNstrips+1];
29 fTracksN = new AliITSdictSSD[fNstrips+1];
30
b0f5e3fc 31
32 fSteps = 10; // still hard-wired - set in SetDetParam and get it via
33 // fDCS together with the others eventually
34
b0f5e3fc 35}
36//___________________________________________________________________________
37AliITSsimulationSSD& AliITSsimulationSSD::operator=(AliITSsimulationSSD
38 &source){
39// Operator =
40 if(this==&source) return *this;
41
42 this->fDCS = new AliITSdcsSSD(*(source.fDCS));
43 this->fN = new TArrayF(*(source.fN));
44 this->fP = new TArrayF(*(source.fP));
45 this->fTracksP = new AliITSdictSSD(*(source.fTracksP));
46 this->fTracksN = new AliITSdictSSD(*(source.fTracksN));
47 this->fNstrips = source.fNstrips;
48 this->fPitch = source.fPitch;
49 this->fSteps = source.fSteps;
50 return *this;
51}
52//_____________________________________________________________
53AliITSsimulationSSD::AliITSsimulationSSD(AliITSsimulationSSD &source){
54 // copy constructor
55 *this = source;
56}
57//____________________________________________________________________________
58AliITSsimulationSSD::~AliITSsimulationSSD() {
59 // anihilator
e8189707 60
b0f5e3fc 61 if(fP) delete fP;
62 if(fN) delete fN;
63
64 if(fTracksP) delete fTracksP;
65 if(fTracksN) delete fTracksN;
66
67 delete fDCS;
68}
69//_______________________________________________________________
70//
71// Hit to digit
72//_______________________________________________________________
73//
74void AliITSsimulationSSD::DigitiseModule(AliITSmodule *mod,Int_t module,
75 Int_t dummy) {
76 // Digitizes one SSD module of hits.
77
78 TObjArray *hits = mod->GetHits();
79 Int_t nhits = hits->GetEntriesFast();
b0f5e3fc 80 if (!nhits) return;
81
82 Int_t i;
83 for(i=0; i<fNstrips; i++) {
84 (*fP)[i] = 0;
85 (*fN)[i] = 0;
86 fTracksP[i].ZeroTracks();
87 fTracksN[i].ZeroTracks();
88 }
89
90 for(i=0; i<nhits; i++) {
91 Int_t idtrack=mod->GetHitTrackIndex(i);
92 HitToDigit(i,idtrack,nhits,hits);
93 }
94
95
e8189707 96
b0f5e3fc 97 ApplyNoise();
98 ApplyCoupling();
99 ApplyThreshold();
100 ApplyDAQ();
101
102
103}
104
105//---------------------------------------------------------------
106
107void AliITSsimulationSSD::HitToDigit(Int_t & hitNo,Int_t idtrack,
108 Int_t nhits,TObjArray *hits) {
109 // Turns one or more hits in an SSD module into one or more digits.
110
111 Int_t stripP, stripN, i;
112 Float_t dsP, dsN;
113 Float_t sP, sN;
114 Float_t eP, eN;
e8189707 115 Float_t arrayEP[10]; // hard-wired number of steps
bba587ca 116 Float_t arrayEN[10];
e8189707 117 Int_t track = -1;
b0f5e3fc 118
119 Float_t ionization = 0;
120 Float_t signal;
121
122 AliITSdictSSD *dict;
123
124
125 // check if this is the right order !!!!!
126
127 AliITShit *hitI = (AliITShit*)hits->At(hitNo++);
128 AliITShit *hitE = (AliITShit*)hits->At(hitNo);
129
130
131 while (!((hitE->StatusExiting()) ||
132 (hitE->StatusDisappeared()) ||
133 (hitE->StatusStop()))) {
134
135 if (++hitNo<nhits) {
136 ionization = hitE->GetIonization();
137 hitE = (AliITShit*)hits->At(hitNo);
138 }
139 }
140
141
142 if (hitI->GetTrack() == hitE->GetTrack())
e8189707 143 //track = idtrack;
b0f5e3fc 144 track = hitI->GetTrack();
145 else
146 printf("!!! Emergency !!!\n");
147
148
149 ionization += hitE->GetIonization();
150
151 const Float_t kconvm=10000.; // cm -> microns
152
153 Float_t xI, yI, zI;
154 hitI->GetPositionL(xI, yI, zI);
155
b0f5e3fc 156 xI *= kconvm;
157 yI *= kconvm;
158 zI *= kconvm;
159
160 Float_t xE, yE, zE;
161 hitE->GetPositionL(xE, yE, zE);
162
b0f5e3fc 163 xE *= kconvm;
164 yE *= kconvm;
165 zE *= kconvm;
166
167 Float_t dx = (xE - xI);
168 Float_t dz = (zE - zI);
169
170
171 // Debuging
172 /*
e8189707 173 fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
b0f5e3fc 174
175 printf("%5d %8.3f %8.3f %8.3f %8.3f %d %d %d\n",
176 hitNo, xI, zI, dx, dz,
177 stripP, stripN, track);
178 printf("%10.5f %10d \n", ionization, hitI->fTrack);
179 */
180
181 // end of debuging
182
183
184 eP=0;
185 eN=0;
b0f5e3fc 186
e8189707 187 for (i=0; i<fSteps; i++) {
b0f5e3fc 188
189 // arrayEP[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
190 // arrayEN[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
191 arrayEP[i] = ionization/fSteps;
192 arrayEN[i] = ionization/fSteps;
193
194 eP += arrayEP[i];
195 eN += arrayEN[i];
196 }
197
198 const Float_t kconv = 1.0e9 / 3.6; // GeV -> e-hole pairs
199
200 for(i=0; i<fSteps; i++) {
201
202 arrayEP[i] = kconv * arrayEP[i] * (ionization / eP);
203 arrayEN[i] = kconv * arrayEN[i] * (ionization / eN);
204 }
205
206 dx /= fSteps;
207 dz /= fSteps;
208
209 Float_t sigmaP, sigmaN;
210 fResponse->SigmaSpread(sigmaP,sigmaN);
211
212 //printf("SigmaP SigmaN %f %f\n",sigmaP, sigmaN);
213
214 Float_t noiseP, noiseN;
215 fResponse->GetNoiseParam(noiseP,noiseN);
216
217 //printf("NoiseP NoiseN %f %f\n",noiseP, noiseN);
218
219 for(i=0; i<fSteps; i++) {
220
221 Int_t j;
222
e8189707 223 fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
224 //printf("hitNo %d i xI zI stripP stripN %d %f %f %d %d\n",hitNo,i,xI, zI, stripP, stripN);
b0f5e3fc 225 dsP = Get2Strip(1,stripP,xI, zI); // Between 0-1
226 dsN = Get2Strip(0,stripN,xI, zI); // Between 0-1
227
228 sP = sigmaP * sqrt(300. * i / (fSteps));
8fdc89c8 229 if(sP<=0.0) sP = sigmaP*sqrt(300.);
b0f5e3fc 230 sN = sigmaN * sqrt(300. * i /(fSteps-i));
8fdc89c8 231 if(sN<=0.0) sN = sigmaN*sqrt(300.);
b0f5e3fc 232
233
234 sP = (i<2 && dsP>0.3 && dsP<0.7)? 20. : sP; // square of (microns)
235 sN = (i>fSteps-2 && dsN>0.3 && dsN<0.7)? 20. : sN; // square of (microns)
236
237 sP = (i==2 && dsP>0.4 && dsP<0.6)? 15. : sP; // square of (microns)
238 sN = (i==8 && dsN>0.4 && dsN<0.6)? 15. : sN; // square of (microns)
239
e8189707 240 for (j=-1; j<2; j++) {
b0f5e3fc 241 if (stripP+j<0 || stripP+j>fNstrips) continue;
b0f5e3fc 242 signal = arrayEP[i] * TMath::Abs( (F(j+0.5-dsP,sP)-F(j-0.5-dsP,sP)) );
243 //printf("SimSSD::HitsToDigits:%d arrayEP[%d]=%e signal=%e\n",j,i,arrayEP[i],signal);
244 if (signal > noiseP/fSteps) {
245 (*fP)[stripP+j] += signal;
246 dict = (fTracksP+stripP+j);
247 (*dict).AddTrack(track);
248 }
249 } // end for j loop over neighboring strips
e8189707 250
251 for (j=-1; j<2; j++) {
b0f5e3fc 252 if (stripN+j<0 || stripN+j>fNstrips) continue;
b0f5e3fc 253 signal = arrayEN[i] * TMath::Abs( (F(j+0.5-dsN,sN)-F(j-0.5-dsN,sN)) );
254 //printf("SimSSD::HitsToDigits:%d arrayEN[%d]=%e signal=%e\n",j,i,arrayEN[i],signal);
255 if (signal > noiseN/fSteps) {
256 (*fN)[stripN+j] += signal;
257 dict = (fTracksN+stripN+j); //co to jest
258 (*dict).AddTrack(track);
259 }
260 } // end for j loop over neighboring strips
261
262 xI += dx;
263 zI += dz;
264 }
265
266
267}
268
269
270//____________________________________________________________________
271//
272// Private Methods for Simulation
273//______________________________________________________________________
274//
275
276void AliITSsimulationSSD::ApplyNoise() {
277 // Apply Noise.
278 Float_t noiseP, noiseN;
279 fResponse->GetNoiseParam(noiseP,noiseN);
280
281 Int_t i;
282 for(i = 0; i<fNstrips; i++) {
283 (*fP)[i] += gRandom->Gaus(0,noiseP);
284 (*fN)[i] += gRandom->Gaus(0,noiseN);
285 }
286}
287
288//_________________________________________________________________________
289
290void AliITSsimulationSSD::ApplyCoupling() {
291 // Apply the effecto of electronic coupling between channels
b0f5e3fc 292 Int_t i;
293 for(i = 1; i<fNstrips-1; i++) {
294 (*fP)[i] += (*fP)[i-1]*fDCS->GetCouplingPL() + (*fP)[i+1]*fDCS->GetCouplingPR();
295 (*fN)[i] += (*fN)[i-1]*fDCS->GetCouplingNL() + (*fN)[i+1]*fDCS->GetCouplingNR();
296 }
297}
298
299//__________________________________________________________________________
300
301void AliITSsimulationSSD::ApplyThreshold() {
302 // Applies the effect of a threshold on the signals for digitization.
303 Float_t noiseP, noiseN;
304 fResponse->GetNoiseParam(noiseP,noiseN);
305
306 // or introduce the SetThresholds in fResponse
307
308 Int_t i;
309 for(i=0; i<fNstrips; i++) {
310 (*fP)[i] = ((*fP)[i] > noiseP*4) ? (*fP)[i] : 0;
311 (*fN)[i] = ((*fN)[i] > noiseN*4) ? (*fN)[i] : 0;
b0f5e3fc 312 }
313
314}
315
316//__________________________________________________________________________
317
318void AliITSsimulationSSD::ApplyDAQ() {
319 // Converts simulated signals to simulated ADC counts
320 AliITS *its=(AliITS*)gAlice->GetModule("ITS");
321
322 Float_t noiseP, noiseN;
323 fResponse->GetNoiseParam(noiseP,noiseN);
324
e8189707 325 char opt[30],dummy[20];
326 fResponse->ParamOptions(opt,dummy);
327
b0f5e3fc 328 Int_t i,j;
e8189707 329 if (strstr(opt,"SetInvalid")) {
330 // Set signal = 0 if invalid strip
331 for(i=0; i<fNstrips; i++) {
332 if (!(fDCS->IsValidP(i))) (*fP)[i] = 0;
333 if (!(fDCS->IsValidN(i))) (*fN)[i] = 0;
334 }
b0f5e3fc 335 }
336
e8189707 337 Int_t digits[3], tracks[3], hits[3];
b0f5e3fc 338 Float_t charges[3];
339 Float_t phys=0;
e8189707 340 for(i=0;i<3;i++) tracks[i]=0;
b0f5e3fc 341 for(i=0; i<fNstrips; i++) {
e8189707 342 if( (strstr(opt,"SetInvalid") && (*fP)[i] < noiseP*4) || !(*fP)[i]) continue;
b0f5e3fc 343 digits[0]=1;
344 digits[1]=i;
345 digits[2]=(int)(*fP)[i];
346 for(j=0; j<(fTracksP+i)->GetNTracks(); j++) {
e8189707 347 if(j>2) continue;
348 if((fTracksP+i)->GetNTracks()) tracks[j]=(fTracksP+i)->GetTrack(j);
349 else tracks[j]=-2;
350 charges[j] = 0;
351 hits[j] = -1;
b0f5e3fc 352 }
e8189707 353 its->AddSimDigit(2,phys,digits,tracks,hits,charges);
b0f5e3fc 354
b0f5e3fc 355 }
356
357
358 for(i=0; i<fNstrips; i++) {
e8189707 359 if( (strstr(opt,"SetInvalid") && (*fN)[i] < noiseN*4)|| !(*fN)[i]) continue;
b0f5e3fc 360 digits[0]=0;
361 digits[1]=i;
362 digits[2]=(int)(*fN)[i];
e8189707 363 for( j=0; j<(fTracksN+i)->GetNTracks(); j++) {
364 if(j>2) continue;
365 if((fTracksN+i)->GetNTracks()) tracks[j]=(fTracksN+i)->GetTrack(j);
366 else tracks[j]=-2;
367 charges[j] = 0;
368 hits[j] = -1;
b0f5e3fc 369 }
e8189707 370 its->AddSimDigit(2,phys,digits,tracks,hits,charges);
b0f5e3fc 371
b0f5e3fc 372 }
373
374}
375
376
377//____________________________________________________________________________
378
379Float_t AliITSsimulationSSD::F(Float_t x, Float_t s) {
380 // Computes the integral of a gaussian at the mean valuse x with sigma s.
e8189707 381
b0f5e3fc 382 //printf("SDD:F(%e,%e)\n",x,s);
383 return 0.5*TMath::Erf(x * fPitch / s) ;
384}
385
386//______________________________________________________________________
387
388Float_t AliITSsimulationSSD::Get2Strip(Int_t flag, Int_t iStrip, Float_t x, Float_t z){
389 // Returns the relative space between two strips.
390
391 // flag==1 for Pside, 0 for Nside
392
393 Float_t stereoP, stereoN;
394 fSegmentation->Angles(stereoP,stereoN);
395
396 Float_t tanP=TMath::Tan(stereoP);
397 Float_t tanN=TMath::Tan(stereoN);
398
399 Float_t dx = fSegmentation->Dx();
400 Float_t dz = fSegmentation->Dz();
401
402
403 x += dx/2;
404 z += dz/2;
405
406 if (flag) return (x - z*tanP) / fPitch - iStrip; // from 0 to 1
407 else return (x - tanN*(dz - z)) / fPitch - iStrip;
408}
409//____________________________________________________________________________