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