Implement fixed due to Boris Batyunya. Including a fix for a bug with the
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSPD.cxx
CommitLineData
b0f5e3fc 1#include <iostream.h>
2#include <TRandom.h>
3#include <TH1.h>
e8189707 4#include <TMath.h>
5
6
b0f5e3fc 7
8#include "AliRun.h"
e8189707 9#include "AliITS.h"
10#include "AliITSMapA2.h"
b0f5e3fc 11#include "AliITSsimulationSPD.h"
e8189707 12
13
b0f5e3fc 14
15ClassImp(AliITSsimulationSPD)
16////////////////////////////////////////////////////////////////////////
17// Version: 0
18// Written by Boris Batyunya
19// December 20 1999
20//
21// AliITSsimulationSPD is the simulation of SPDs
b0f5e3fc 22//________________________________________________________________________
23
e8189707 24
25AliITSsimulationSPD::AliITSsimulationSPD()
26{
b0f5e3fc 27 // constructor
e8189707 28 fResponse = 0;
b0f5e3fc 29 fSegmentation = 0;
e8189707 30 fHis = 0;
31 fNoise=0.;
32 fBaseline=0.;
b0f5e3fc 33}
e8189707 34
35
b0f5e3fc 36//_____________________________________________________________________________
37
38AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg, AliITSresponse *resp) {
e8189707 39 // standard constructor
40
b0f5e3fc 41 fResponse = resp;
42 fSegmentation = seg;
43
44 fResponse->GetNoiseParam(fNoise,fBaseline);
45
b0f5e3fc 46 fMapA2 = new AliITSMapA2(fSegmentation);
e8189707 47
48 //
b0f5e3fc 49 fNPixelsZ=fSegmentation->Npz();
50 fNPixelsX=fSegmentation->Npx();
51
52}
53
54//_____________________________________________________________________________
55
56AliITSsimulationSPD::~AliITSsimulationSPD() {
57 // destructor
58
59 delete fMapA2;
60
61 if (fHis) {
62 fHis->Delete();
63 delete fHis;
64 }
65}
66
e8189707 67
b0f5e3fc 68//__________________________________________________________________________
69AliITSsimulationSPD::AliITSsimulationSPD(const AliITSsimulationSPD &source){
70 // Copy Constructor
71 if(&source == this) return;
72 this->fMapA2 = source.fMapA2;
73 this->fNoise = source.fNoise;
74 this->fBaseline = source.fBaseline;
75 this->fNPixelsX = source.fNPixelsX;
76 this->fNPixelsZ = source.fNPixelsZ;
77 this->fHis = source.fHis;
78 return;
79}
80
81//_________________________________________________________________________
82AliITSsimulationSPD&
83 AliITSsimulationSPD::operator=(const AliITSsimulationSPD &source) {
84 // Assignment operator
85 if(&source == this) return *this;
86 this->fMapA2 = source.fMapA2;
87 this->fNoise = source.fNoise;
88 this->fBaseline = source.fBaseline;
89 this->fNPixelsX = source.fNPixelsX;
90 this->fNPixelsZ = source.fNPixelsZ;
91 this->fHis = source.fHis;
92 return *this;
93 }
94//_____________________________________________________________________________
95
e8189707 96void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t module, Int_t dummy)
97{
b0f5e3fc 98 // digitize module
e8189707 99
100 const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons
b0f5e3fc 101 // for 3.6 eV/pair
102 const Float_t kconv = 10000.; // cm -> microns
103
e8189707 104 Float_t spdLength = fSegmentation->Dz();
b0f5e3fc 105 Float_t spdWidth = fSegmentation->Dx();
106
e8189707 107 Float_t difCoef, dum;
108 fResponse->DiffCoeff(difCoef,dum);
b0f5e3fc 109
110 Float_t zPix0 = 1e+6;
111 Float_t xPix0 = 1e+6;
112 Float_t yPix0 = 1e+6;
e8189707 113 Float_t yPrev = 1e+6;
114 Float_t zP0 = 100.;
115 Float_t xP0 = 100.;
b0f5e3fc 116
117 Float_t zPitch = fSegmentation->Dpz(0);
118 Float_t xPitch = fSegmentation->Dpx(0);
119
120 //cout << "pitch per z: " << zPitch << endl;
121 //cout << "pitch per r*phi: " << xPitch << endl;
122
123 TObjArray *fHits = mod->GetHits();
124 Int_t nhits = fHits->GetEntriesFast();
125 if (!nhits) return;
a3e16987 126 // cout << "module, nhits ="<<module<<","<<nhits<< endl;
b0f5e3fc 127
128 // Array of pointers to the label-signal list
129
e8189707 130 Int_t maxNDigits = fNPixelsX*fNPixelsZ + fNPixelsX ;;
131 Float_t **pList = new Float_t* [maxNDigits];
132 memset(pList,0,sizeof(Float_t*)*maxNDigits);
b0f5e3fc 133 Int_t indexRange[4] = {0,0,0,0};
134
135 // Fill detector maps with GEANT hits
136 // loop over hits in the module
a3e16987 137 static Bool_t first;
e8189707 138 Int_t lasttrack=-2;
139 Int_t hit, iZi, jz, jx;
140 for (hit=0;hit<nhits;hit++) {
b0f5e3fc 141 AliITShit *iHit = (AliITShit*) fHits->At(hit);
e8189707 142 Int_t layer = iHit->GetLayer();
143
b0f5e3fc 144 // work with the idtrack=entry number in the TreeH
e8189707 145 Int_t idhit,idtrack;
146 mod->GetHitTrackAndHitIndex(hit,idtrack,idhit);
147 //Int_t idtrack=mod->GetHitTrackIndex(hit);
b0f5e3fc 148 // or store straight away the particle position in the array
149 // of particles :
e8189707 150 Int_t itrack = iHit->GetTrack();
a3e16987 151 Int_t dray = 0;
b0f5e3fc 152
a3e16987 153 if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE;
154
155 // Int_t parent = iHit->GetParticle()->GetFirstMother();
156 Int_t partcode = iHit->GetParticle()->GetPdgCode();
157
158// partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0, 211 - pi+
159// 310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
160
161 Float_t px = iHit->GetPXL();
162 Float_t py = iHit->GetPYL();
163 Float_t pz = iHit->GetPZL();
164 Float_t pmod = 1000*sqrt(px*px+py*py+pz*pz);
165
166
167 if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
168 // at p < 6 MeV/c
169
170
b0f5e3fc 171 // Get hit z and x(r*phi) cordinates for each module (detector)
172 // in local system.
173
e8189707 174 Float_t zPix = kconv*iHit->GetZL();
175 Float_t xPix = kconv*iHit->GetXL();
176 Float_t yPix = kconv*iHit->GetYL();
b0f5e3fc 177
178 // Get track status
e8189707 179 Int_t status = iHit->GetTrackStatus();
180
181 // Check boundaries
182 if(TMath::Abs(zPix) > spdLength/2.) {
183 printf("!! Zpix outside = %f\n",zPix);
184 if(status == 66) zP0=100;
185 continue;
186 }
187
188
189 if (TMath::Abs(xPix) > spdWidth/2.) {
190 printf("!! Xpix outside = %f\n",xPix);
191 if (status == 66) xP0=100;
192 continue;
193 }
194
195 Float_t zP = (zPix + spdLength/2.)/1000.;
196 Float_t xP = (xPix + spdWidth/2.)/1000.;
b0f5e3fc 197
e8189707 198 Int_t trdown = 0;
b0f5e3fc 199
200 // enter Si or after event in Si
201 if (status == 66 ) {
202 zPix0 = zPix;
203 xPix0 = xPix;
204 yPrev = yPix;
205 }
206 // enter Si only
207 if (layer == 1 && status == 66 && yPix > 71.) {
208 yPix0 = yPix;
e8189707 209 zP0 = zP;
210 xP0 = xP;
b0f5e3fc 211 }
212 // enter Si only
213 if (layer == 2 && status == 66 && yPix < -71.) {
214 yPix0 = yPix;
e8189707 215 zP0 = zP;
216 xP0 = xP;
217 }
218 Float_t depEnergy = iHit->GetIonization();
b0f5e3fc 219 // skip if the input point to Si
e8189707 220 if(depEnergy <= 0.) continue;
221 // skip if the input point is outside of Si, but the next
222 // point is inside of Si
223 if(zP0 > 90 || xP0 > 90) continue;
b0f5e3fc 224 // if track returns to the opposite direction:
225 if (layer == 1 && yPix > yPrev) {
226 yPix0 = yPrev;
227 trdown = 1;
228 }
229 if (layer == 2 && yPix < yPrev) {
230 yPix0 = yPrev;
231 trdown = 1;
e8189707 232 }
b0f5e3fc 233
234 // take into account the holes diffusion inside the Silicon
235 // the straight line between the entrance and exit points in Si is
236 // divided into the several steps; the diffusion is considered
237 // for each end point of step and charge
238 // is distributed between the pixels through the diffusion.
239
240
241 // ---------- the diffusion in Z (beam) direction -------
242
e8189707 243 Float_t charge = depEnergy*kEnToEl; // charge in e-
244 Float_t drPath = 0.;
245 Float_t tang = 0.;
246 Float_t sigmaDif = 0.;
247 Float_t zdif = zPix - zPix0;
248 Float_t xdif = xPix - xPix0;
249 Float_t ydif = yPix - yPix0;
b0f5e3fc 250
e8189707 251 if(TMath::Abs(ydif) < 0.1) continue; // Ydif is not zero
b0f5e3fc 252
e8189707 253 Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
254 Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
255 Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1;
b0f5e3fc 256
257 // number of the steps along the track:
e8189707 258 Int_t nsteps = ndZ;
b0f5e3fc 259 if(ndX > ndZ) nsteps = ndX;
260 if(nsteps < 6) nsteps = 6; // minimum number of the steps
261
e8189707 262 if(TMath::Abs(projDif) > 5.0) tang = ydif/projDif;
263 Float_t dCharge = charge/nsteps; // charge in e- for one step
264 Float_t dZ = zdif/nsteps;
265 Float_t dX = xdif/nsteps;
b0f5e3fc 266
267 if (TMath::Abs(projDif) < 5.0 ) {
e8189707 268 drPath = ydif*1.e-4;
b0f5e3fc 269 drPath = TMath::Abs(drPath); // drift path in cm
270 sigmaDif = difCoef*sqrt(drPath); // sigma diffusion in cm
271 }
272
e8189707 273 for (iZi = 1;iZi <= nsteps;iZi++) {
274 Float_t dZn = iZi*dZ;
275 Float_t dXn = iZi*dX;
276 Float_t zPixn = zPix0 + dZn;
277 Float_t xPixn = xPix0 + dXn;
b0f5e3fc 278
279 if(TMath::Abs(projDif) >= 5.) {
e8189707 280 Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
b0f5e3fc 281 if(trdown == 0) {
e8189707 282 drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm
b0f5e3fc 283 drPath = TMath::Abs(drPath);
284 }
285 if(trdown == 1) {
e8189707 286 Float_t dProjn = projDif/nsteps;
287 drPath = (projDif-(iZi-1)*dProjn)*tang*1.e-4;
b0f5e3fc 288 drPath = TMath::Abs(drPath);
289 }
290 sigmaDif = difCoef*sqrt(drPath);
291 sigmaDif = sigmaDif*kconv; // sigma diffusion in microns
292 }
e8189707 293 zPixn = (zPixn + spdLength/2.);
294 xPixn = (xPixn + spdWidth/2.);
295 Int_t nZpix, nXpix;
296 fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
b0f5e3fc 297 zPitch = fSegmentation->Dpz(nZpix);
e8189707 298 fSegmentation->GetPadTxz(xPixn,zPixn);
b0f5e3fc 299 // set the window for the integration
e8189707 300 Int_t jzmin = 1;
301 Int_t jzmax = 3;
b0f5e3fc 302 if(nZpix == 1) jzmin =2;
303 if(nZpix == fNPixelsZ) jzmax = 2;
304
e8189707 305 Int_t jxmin = 1;
306 Int_t jxmax = 3;
b0f5e3fc 307 if(nXpix == 1) jxmin =2;
308 if(nXpix == fNPixelsX) jxmax = 2;
309
e8189707 310 Float_t zpix = nZpix;
311 Float_t dZright = zPitch*(zpix - zPixn);
312 Float_t dZleft = zPitch - dZright;
313
314 Float_t xpix = nXpix;
315 Float_t dXright = xPitch*(xpix - xPixn);
316 Float_t dXleft = xPitch - dXright;
317
318 Float_t dZprev = 0.;
319 Float_t dZnext = 0.;
320 Float_t dXprev = 0.;
321 Float_t dXnext = 0.;
b0f5e3fc 322
323 for(jz=jzmin; jz <=jzmax; jz++) {
324 if(jz == 1) {
325 dZprev = -zPitch - dZleft;
326 dZnext = -dZleft;
327 }
328 if(jz == 2) {
329 dZprev = -dZleft;
330 dZnext = dZright;
331 }
332 if(jz == 3) {
333 dZprev = dZright;
334 dZnext = dZright + zPitch;
335 }
e8189707 336 // kz changes from 1 to the fNofPixels(270)
337 Int_t kz = nZpix + jz -2;
b0f5e3fc 338
e8189707 339 Float_t zArg1 = dZprev/sigmaDif;
340 Float_t zArg2 = dZnext/sigmaDif;
341 Float_t zProb1 = TMath::Erfc(zArg1);
342 Float_t zProb2 = TMath::Erfc(zArg2);
343 Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge;
b0f5e3fc 344
b0f5e3fc 345
346 // ----------- holes diffusion in X(r*phi) direction --------
347
348 if(dZCharge > 1.) {
349 for(jx=jxmin; jx <=jxmax; jx++) {
350 if(jx == 1) {
351 dXprev = -xPitch - dXleft;
352 dXnext = -dXleft;
353 }
354 if(jx == 2) {
355 dXprev = -dXleft;
356 dXnext = dXright;
357 }
358 if(jx == 3) {
359 dXprev = dXright;
360 dXnext = dXright + xPitch;
361 }
e8189707 362 Int_t kx = nXpix + jx -2;
b0f5e3fc 363
e8189707 364 Float_t xArg1 = dXprev/sigmaDif;
365 Float_t xArg2 = dXnext/sigmaDif;
366 Float_t xProb1 = TMath::Erfc(xArg1);
367 Float_t xProb2 = TMath::Erfc(xArg2);
368 Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge;
b0f5e3fc 369
370 if(dXCharge > 1.) {
e8189707 371 Int_t index = kz-1;
372
b0f5e3fc 373 if (first) {
374 indexRange[0]=indexRange[1]=index;
e8189707 375 indexRange[2]=indexRange[3]=kx-1;
b0f5e3fc 376 first=kFALSE;
377 }
378
e8189707 379 indexRange[0]=TMath::Min(indexRange[0],kz-1);
380 indexRange[1]=TMath::Max(indexRange[1],kz-1);
381 indexRange[2]=TMath::Min(indexRange[2],kx-1);
382 indexRange[3]=TMath::Max(indexRange[3],kx-1);
383
b0f5e3fc 384 // build the list of digits for this module
e8189707 385 Double_t signal=fMapA2->GetSignal(index,kx-1);
b0f5e3fc 386 signal+=dXCharge;
e8189707 387 fMapA2->SetHit(index,kx-1,(double)signal);
b0f5e3fc 388 } // dXCharge > 1 e-
389 } // jx loop
390 } // dZCharge > 1 e-
391 } // jz loop
392 } // iZi loop
393
394 if (status == 65) { // the step is inside of Si
395 zPix0 = zPix;
396 xPix0 = xPix;
397 }
a3e16987 398 yPrev = yPix;
b0f5e3fc 399
a3e16987 400 if(dray == 0) {
e8189707 401 GetList(itrack,idhit,pList,indexRange);
b0f5e3fc 402 }
e8189707 403
404 lasttrack=itrack;
b0f5e3fc 405 } // hit loop inside the module
406
407
408 // introduce the electronics effects and do zero-suppression
409 ChargeToSignal(pList);
410
411 // clean memory
412
413 fMapA2->ClearMap();
414
415
416}
417
418//---------------------------------------------
e8189707 419void AliITSsimulationSPD::GetList(Int_t label,Int_t idhit,Float_t **pList,Int_t *indexRange)
420{
421 // lop over nonzero digits
b0f5e3fc 422
e8189707 423
424 //set protection
425 for(int k=0;k<4;k++) {
426 if (indexRange[k] < 0) indexRange[k]=0;
427 }
b0f5e3fc 428
e8189707 429 for(Int_t iz=indexRange[0];iz<indexRange[1]+1;iz++){
430 for(Int_t ix=indexRange[2];ix<indexRange[3]+1;ix++){
b0f5e3fc 431
e8189707 432 Float_t signal=fMapA2->GetSignal(iz,ix);
a3e16987 433
e8189707 434 if (!signal) continue;
435
436 Int_t globalIndex = iz*fNPixelsX+ix; // GlobalIndex starts from 0!
b0f5e3fc 437 if(!pList[globalIndex]){
a3e16987 438
b0f5e3fc 439 //
e8189707 440 // Create new list (9 elements - 3 signals and 3 tracks + 3 hits)
b0f5e3fc 441 //
442
e8189707 443 pList[globalIndex] = new Float_t [9];
b0f5e3fc 444
e8189707 445 // set list to -3
b0f5e3fc 446
e8189707 447 *pList[globalIndex] = -3.;
448 *(pList[globalIndex]+1) = -3.;
449 *(pList[globalIndex]+2) = -3.;
b0f5e3fc 450 *(pList[globalIndex]+3) = 0.;
451 *(pList[globalIndex]+4) = 0.;
452 *(pList[globalIndex]+5) = 0.;
e8189707 453 *(pList[globalIndex]+6) = -1.;
454 *(pList[globalIndex]+7) = -1.;
455 *(pList[globalIndex]+8) = -1.;
b0f5e3fc 456
457
458 *pList[globalIndex] = (float)label;
459 *(pList[globalIndex]+3) = signal;
e8189707 460 *(pList[globalIndex]+6) = (float)idhit;
b0f5e3fc 461 }
462 else{
463
464 // check the signal magnitude
465
e8189707 466 Float_t highest = *(pList[globalIndex]+3);
467 Float_t middle = *(pList[globalIndex]+4);
468 Float_t lowest = *(pList[globalIndex]+5);
b0f5e3fc 469
470 signal -= (highest+middle+lowest);
471
472 //
473 // compare the new signal with already existing list
474 //
475
476 if(signal<lowest) continue; // neglect this track
477
478 if (signal>highest){
479 *(pList[globalIndex]+5) = middle;
480 *(pList[globalIndex]+4) = highest;
481 *(pList[globalIndex]+3) = signal;
482
483 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
484 *(pList[globalIndex]+1) = *pList[globalIndex];
485 *pList[globalIndex] = label;
e8189707 486
487 *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
488 *(pList[globalIndex]+7) = *(pList[globalIndex]+6);
489 *(pList[globalIndex]+6) = idhit;
b0f5e3fc 490 }
491 else if (signal>middle){
492 *(pList[globalIndex]+5) = middle;
493 *(pList[globalIndex]+4) = signal;
494
495 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
496 *(pList[globalIndex]+1) = label;
e8189707 497
498 *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
499 *(pList[globalIndex]+7) = idhit;
b0f5e3fc 500 }
501 else{
502 *(pList[globalIndex]+5) = signal;
503 *(pList[globalIndex]+2) = label;
e8189707 504 *(pList[globalIndex]+8) = idhit;
b0f5e3fc 505 }
506 }
507 } // end of loop pixels in x
508 } // end of loop over pixels in z
509
510
511}
512
513
514//---------------------------------------------
e8189707 515void AliITSsimulationSPD::ChargeToSignal(Float_t **pList)
516{
517 // add noise and electronics, perform the zero suppression and add the
518 // digit to the list
b0f5e3fc 519
520 AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
521
522
523 TRandom *random = new TRandom();
524 Float_t threshold = (float)fResponse->MinVal();
525
e8189707 526 Int_t digits[3], tracks[3], hits[3],gi,j1;
b0f5e3fc 527 Float_t charges[3];
528 Float_t electronics;
529 Float_t signal,phys;
e8189707 530 for(Int_t iz=0;iz<fNPixelsZ;iz++){
531 for(Int_t ix=0;ix<fNPixelsX;ix++){
b0f5e3fc 532 electronics = fBaseline + fNoise*random->Gaus();
533 signal = (float)fMapA2->GetSignal(iz,ix);
534 signal += electronics;
e8189707 535 gi =iz*fNPixelsX+ix; // global index
b0f5e3fc 536 if (signal > threshold) {
537 digits[0]=iz;
538 digits[1]=ix;
539 digits[2]=1;
b0f5e3fc 540 for(j1=0;j1<3;j1++){
e8189707 541 if (pList[gi]) {
542 tracks[j1] = (Int_t)(*(pList[gi]+j1));
543 hits[j1] = (Int_t)(*(pList[gi]+j1+6));
544 }else {
545 tracks[j1]=-2; //noise
546 hits[j1] = -1;
547 }
b0f5e3fc 548 charges[j1] = 0;
549 }
a3e16987 550
551 if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) {
552 tracks[1] = -3;
553 hits[1] = -1;
554 tracks[2] = -3;
555 hits[2] = -1;
556 }
557 if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) {
558 tracks[1] = -3;
559 hits[1] = -1;
560 }
561 if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) {
562 tracks[2] = -3;
563 hits[2] = -1;
564 }
565 if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) {
566 tracks[2] = -3;
567 hits[2] = -1;
568 }
569
b0f5e3fc 570 phys=0;
e8189707 571 aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
b0f5e3fc 572 }
e8189707 573 if(pList[gi]) delete [] pList[gi];
b0f5e3fc 574 }
575 }
576 delete [] pList;
577
b0f5e3fc 578}
579
580
581//____________________________________________
582
e8189707 583void AliITSsimulationSPD::CreateHistograms()
584{
585 // create 1D histograms for tests
586
587 printf("SPD - create histograms\n");
b0f5e3fc 588
e8189707 589 for (Int_t i=0;i<fNPixelsZ;i++) {
b0f5e3fc 590 TString *spdname = new TString("spd_");
e8189707 591 Char_t pixelz[4];
592 sprintf(pixelz,"%d",i+1);
593 spdname->Append(pixelz);
b0f5e3fc 594 (*fHis)[i] = new TH1F(spdname->Data(),"SPD maps",
595 fNPixelsX,0.,(Float_t) fNPixelsX);
596 delete spdname;
597 }
598
599}
600
601//____________________________________________
602
e8189707 603void AliITSsimulationSPD::ResetHistograms()
604{
b0f5e3fc 605 //
606 // Reset histograms for this detector
607 //
e8189707 608 for ( int i=0;i<fNPixelsZ;i++ ) {
b0f5e3fc 609 if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();
610 }
611
612}
e8189707 613
614
615
616
617
618
619
620
621