Release version of ITS code
[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;
126
127
128
129 // Array of pointers to the label-signal list
130
e8189707 131 Int_t maxNDigits = fNPixelsX*fNPixelsZ + fNPixelsX ;;
132 Float_t **pList = new Float_t* [maxNDigits];
133 memset(pList,0,sizeof(Float_t*)*maxNDigits);
b0f5e3fc 134 Int_t indexRange[4] = {0,0,0,0};
135
136 // Fill detector maps with GEANT hits
137 // loop over hits in the module
b0f5e3fc 138 static Bool_t first=kTRUE;
e8189707 139 Int_t lasttrack=-2;
140 Int_t hit, iZi, jz, jx;
141 for (hit=0;hit<nhits;hit++) {
b0f5e3fc 142 AliITShit *iHit = (AliITShit*) fHits->At(hit);
e8189707 143 Int_t layer = iHit->GetLayer();
144
b0f5e3fc 145
146 // work with the idtrack=entry number in the TreeH
e8189707 147 Int_t idhit,idtrack;
148 mod->GetHitTrackAndHitIndex(hit,idtrack,idhit);
149 //Int_t idtrack=mod->GetHitTrackIndex(hit);
b0f5e3fc 150 // or store straight away the particle position in the array
151 // of particles :
e8189707 152 Int_t itrack = iHit->GetTrack();
b0f5e3fc 153
154 // Get hit z and x(r*phi) cordinates for each module (detector)
155 // in local system.
156
e8189707 157 Float_t zPix = kconv*iHit->GetZL();
158 Float_t xPix = kconv*iHit->GetXL();
159 Float_t yPix = kconv*iHit->GetYL();
b0f5e3fc 160
161 // Get track status
e8189707 162 Int_t status = iHit->GetTrackStatus();
163
164 // Check boundaries
165 if(TMath::Abs(zPix) > spdLength/2.) {
166 printf("!! Zpix outside = %f\n",zPix);
167 if(status == 66) zP0=100;
168 continue;
169 }
170
171
172 if (TMath::Abs(xPix) > spdWidth/2.) {
173 printf("!! Xpix outside = %f\n",xPix);
174 if (status == 66) xP0=100;
175 continue;
176 }
177
178 Float_t zP = (zPix + spdLength/2.)/1000.;
179 Float_t xP = (xPix + spdWidth/2.)/1000.;
b0f5e3fc 180
e8189707 181 Int_t trdown = 0;
b0f5e3fc 182
183 // enter Si or after event in Si
184 if (status == 66 ) {
185 zPix0 = zPix;
186 xPix0 = xPix;
187 yPrev = yPix;
188 }
189 // enter Si only
190 if (layer == 1 && status == 66 && yPix > 71.) {
191 yPix0 = yPix;
e8189707 192 zP0 = zP;
193 xP0 = xP;
b0f5e3fc 194 }
195 // enter Si only
196 if (layer == 2 && status == 66 && yPix < -71.) {
197 yPix0 = yPix;
e8189707 198 zP0 = zP;
199 xP0 = xP;
200 }
201 Float_t depEnergy = iHit->GetIonization();
b0f5e3fc 202 // skip if the input point to Si
e8189707 203 if(depEnergy <= 0.) continue;
204 // skip if the input point is outside of Si, but the next
205 // point is inside of Si
206 if(zP0 > 90 || xP0 > 90) continue;
b0f5e3fc 207 // if track returns to the opposite direction:
208 if (layer == 1 && yPix > yPrev) {
209 yPix0 = yPrev;
210 trdown = 1;
211 }
212 if (layer == 2 && yPix < yPrev) {
213 yPix0 = yPrev;
214 trdown = 1;
e8189707 215 }
b0f5e3fc 216
217 // take into account the holes diffusion inside the Silicon
218 // the straight line between the entrance and exit points in Si is
219 // divided into the several steps; the diffusion is considered
220 // for each end point of step and charge
221 // is distributed between the pixels through the diffusion.
222
223
224 // ---------- the diffusion in Z (beam) direction -------
225
e8189707 226 Float_t charge = depEnergy*kEnToEl; // charge in e-
227 Float_t drPath = 0.;
228 Float_t tang = 0.;
229 Float_t sigmaDif = 0.;
230 Float_t zdif = zPix - zPix0;
231 Float_t xdif = xPix - xPix0;
232 Float_t ydif = yPix - yPix0;
b0f5e3fc 233
e8189707 234 if(TMath::Abs(ydif) < 0.1) continue; // Ydif is not zero
b0f5e3fc 235
e8189707 236 Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
237 Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
238 Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1;
b0f5e3fc 239
240 // number of the steps along the track:
e8189707 241 Int_t nsteps = ndZ;
b0f5e3fc 242 if(ndX > ndZ) nsteps = ndX;
243 if(nsteps < 6) nsteps = 6; // minimum number of the steps
244
e8189707 245 if(TMath::Abs(projDif) > 5.0) tang = ydif/projDif;
246 Float_t dCharge = charge/nsteps; // charge in e- for one step
247 Float_t dZ = zdif/nsteps;
248 Float_t dX = xdif/nsteps;
b0f5e3fc 249
250 if (TMath::Abs(projDif) < 5.0 ) {
e8189707 251 drPath = ydif*1.e-4;
b0f5e3fc 252 drPath = TMath::Abs(drPath); // drift path in cm
253 sigmaDif = difCoef*sqrt(drPath); // sigma diffusion in cm
254 }
255
e8189707 256 for (iZi = 1;iZi <= nsteps;iZi++) {
257 Float_t dZn = iZi*dZ;
258 Float_t dXn = iZi*dX;
259 Float_t zPixn = zPix0 + dZn;
260 Float_t xPixn = xPix0 + dXn;
b0f5e3fc 261
262 if(TMath::Abs(projDif) >= 5.) {
e8189707 263 Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
b0f5e3fc 264 if(trdown == 0) {
e8189707 265 drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm
b0f5e3fc 266 drPath = TMath::Abs(drPath);
267 }
268 if(trdown == 1) {
e8189707 269 Float_t dProjn = projDif/nsteps;
270 drPath = (projDif-(iZi-1)*dProjn)*tang*1.e-4;
b0f5e3fc 271 drPath = TMath::Abs(drPath);
272 }
273 sigmaDif = difCoef*sqrt(drPath);
274 sigmaDif = sigmaDif*kconv; // sigma diffusion in microns
275 }
e8189707 276 zPixn = (zPixn + spdLength/2.);
277 xPixn = (xPixn + spdWidth/2.);
278 Int_t nZpix, nXpix;
279 fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
b0f5e3fc 280 zPitch = fSegmentation->Dpz(nZpix);
e8189707 281 fSegmentation->GetPadTxz(xPixn,zPixn);
b0f5e3fc 282 // set the window for the integration
e8189707 283 Int_t jzmin = 1;
284 Int_t jzmax = 3;
b0f5e3fc 285 if(nZpix == 1) jzmin =2;
286 if(nZpix == fNPixelsZ) jzmax = 2;
287
e8189707 288 Int_t jxmin = 1;
289 Int_t jxmax = 3;
b0f5e3fc 290 if(nXpix == 1) jxmin =2;
291 if(nXpix == fNPixelsX) jxmax = 2;
292
e8189707 293 Float_t zpix = nZpix;
294 Float_t dZright = zPitch*(zpix - zPixn);
295 Float_t dZleft = zPitch - dZright;
296
297 Float_t xpix = nXpix;
298 Float_t dXright = xPitch*(xpix - xPixn);
299 Float_t dXleft = xPitch - dXright;
300
301 Float_t dZprev = 0.;
302 Float_t dZnext = 0.;
303 Float_t dXprev = 0.;
304 Float_t dXnext = 0.;
b0f5e3fc 305
306 for(jz=jzmin; jz <=jzmax; jz++) {
307 if(jz == 1) {
308 dZprev = -zPitch - dZleft;
309 dZnext = -dZleft;
310 }
311 if(jz == 2) {
312 dZprev = -dZleft;
313 dZnext = dZright;
314 }
315 if(jz == 3) {
316 dZprev = dZright;
317 dZnext = dZright + zPitch;
318 }
e8189707 319 // kz changes from 1 to the fNofPixels(270)
320 Int_t kz = nZpix + jz -2;
b0f5e3fc 321
e8189707 322 Float_t zArg1 = dZprev/sigmaDif;
323 Float_t zArg2 = dZnext/sigmaDif;
324 Float_t zProb1 = TMath::Erfc(zArg1);
325 Float_t zProb2 = TMath::Erfc(zArg2);
326 Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge;
b0f5e3fc 327
b0f5e3fc 328
329 // ----------- holes diffusion in X(r*phi) direction --------
330
331 if(dZCharge > 1.) {
332 for(jx=jxmin; jx <=jxmax; jx++) {
333 if(jx == 1) {
334 dXprev = -xPitch - dXleft;
335 dXnext = -dXleft;
336 }
337 if(jx == 2) {
338 dXprev = -dXleft;
339 dXnext = dXright;
340 }
341 if(jx == 3) {
342 dXprev = dXright;
343 dXnext = dXright + xPitch;
344 }
e8189707 345 Int_t kx = nXpix + jx -2;
b0f5e3fc 346
e8189707 347 Float_t xArg1 = dXprev/sigmaDif;
348 Float_t xArg2 = dXnext/sigmaDif;
349 Float_t xProb1 = TMath::Erfc(xArg1);
350 Float_t xProb2 = TMath::Erfc(xArg2);
351 Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge;
b0f5e3fc 352
353 if(dXCharge > 1.) {
e8189707 354 Int_t index = kz-1;
355
b0f5e3fc 356 if (first) {
357 indexRange[0]=indexRange[1]=index;
e8189707 358 indexRange[2]=indexRange[3]=kx-1;
b0f5e3fc 359 first=kFALSE;
360 }
361
e8189707 362 indexRange[0]=TMath::Min(indexRange[0],kz-1);
363 indexRange[1]=TMath::Max(indexRange[1],kz-1);
364 indexRange[2]=TMath::Min(indexRange[2],kx-1);
365 indexRange[3]=TMath::Max(indexRange[3],kx-1);
366
b0f5e3fc 367 // build the list of digits for this module
e8189707 368 Double_t signal=fMapA2->GetSignal(index,kx-1);
b0f5e3fc 369 signal+=dXCharge;
e8189707 370 fMapA2->SetHit(index,kx-1,(double)signal);
b0f5e3fc 371 } // dXCharge > 1 e-
372 } // jx loop
373 } // dZCharge > 1 e-
374 } // jz loop
375 } // iZi loop
376
377 if (status == 65) { // the step is inside of Si
378 zPix0 = zPix;
379 xPix0 = xPix;
380 }
381 yPrev = yPix; //ch
382
e8189707 383 if (lasttrack != itrack || hit==(nhits-1)) {
384 GetList(itrack,idhit,pList,indexRange);
b0f5e3fc 385 first=kTRUE;
386 }
e8189707 387
388 lasttrack=itrack;
b0f5e3fc 389 } // hit loop inside the module
390
391
392 // introduce the electronics effects and do zero-suppression
393 ChargeToSignal(pList);
394
395 // clean memory
396
397 fMapA2->ClearMap();
398
399
400}
401
402//---------------------------------------------
e8189707 403void AliITSsimulationSPD::GetList(Int_t label,Int_t idhit,Float_t **pList,Int_t *indexRange)
404{
405 // lop over nonzero digits
b0f5e3fc 406
e8189707 407
408 //set protection
409 for(int k=0;k<4;k++) {
410 if (indexRange[k] < 0) indexRange[k]=0;
411 }
b0f5e3fc 412
e8189707 413 for(Int_t iz=indexRange[0];iz<indexRange[1]+1;iz++){
414 for(Int_t ix=indexRange[2];ix<indexRange[3]+1;ix++){
b0f5e3fc 415
e8189707 416 Float_t signal=fMapA2->GetSignal(iz,ix);
417 if (!signal) continue;
418
419 Int_t globalIndex = iz*fNPixelsX+ix; // GlobalIndex starts from 0!
b0f5e3fc 420 if(!pList[globalIndex]){
421
422 //
e8189707 423 // Create new list (9 elements - 3 signals and 3 tracks + 3 hits)
b0f5e3fc 424 //
425
e8189707 426 pList[globalIndex] = new Float_t [9];
b0f5e3fc 427
e8189707 428 // set list to -3
b0f5e3fc 429
e8189707 430 *pList[globalIndex] = -3.;
431 *(pList[globalIndex]+1) = -3.;
432 *(pList[globalIndex]+2) = -3.;
b0f5e3fc 433 *(pList[globalIndex]+3) = 0.;
434 *(pList[globalIndex]+4) = 0.;
435 *(pList[globalIndex]+5) = 0.;
e8189707 436 *(pList[globalIndex]+6) = -1.;
437 *(pList[globalIndex]+7) = -1.;
438 *(pList[globalIndex]+8) = -1.;
b0f5e3fc 439
440
441 *pList[globalIndex] = (float)label;
442 *(pList[globalIndex]+3) = signal;
e8189707 443 *(pList[globalIndex]+6) = (float)idhit;
b0f5e3fc 444 }
445 else{
446
447 // check the signal magnitude
448
e8189707 449 Float_t highest = *(pList[globalIndex]+3);
450 Float_t middle = *(pList[globalIndex]+4);
451 Float_t lowest = *(pList[globalIndex]+5);
b0f5e3fc 452
453 signal -= (highest+middle+lowest);
454
455 //
456 // compare the new signal with already existing list
457 //
458
459 if(signal<lowest) continue; // neglect this track
460
461 if (signal>highest){
462 *(pList[globalIndex]+5) = middle;
463 *(pList[globalIndex]+4) = highest;
464 *(pList[globalIndex]+3) = signal;
465
466 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
467 *(pList[globalIndex]+1) = *pList[globalIndex];
468 *pList[globalIndex] = label;
e8189707 469
470 *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
471 *(pList[globalIndex]+7) = *(pList[globalIndex]+6);
472 *(pList[globalIndex]+6) = idhit;
b0f5e3fc 473 }
474 else if (signal>middle){
475 *(pList[globalIndex]+5) = middle;
476 *(pList[globalIndex]+4) = signal;
477
478 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
479 *(pList[globalIndex]+1) = label;
e8189707 480
481 *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
482 *(pList[globalIndex]+7) = idhit;
b0f5e3fc 483 }
484 else{
485 *(pList[globalIndex]+5) = signal;
486 *(pList[globalIndex]+2) = label;
e8189707 487 *(pList[globalIndex]+8) = idhit;
b0f5e3fc 488 }
489 }
490 } // end of loop pixels in x
491 } // end of loop over pixels in z
492
493
494}
495
496
497//---------------------------------------------
e8189707 498void AliITSsimulationSPD::ChargeToSignal(Float_t **pList)
499{
500 // add noise and electronics, perform the zero suppression and add the
501 // digit to the list
b0f5e3fc 502
503 AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
504
505
506 TRandom *random = new TRandom();
507 Float_t threshold = (float)fResponse->MinVal();
508
e8189707 509 Int_t digits[3], tracks[3], hits[3],gi,j1;
b0f5e3fc 510 Float_t charges[3];
511 Float_t electronics;
512 Float_t signal,phys;
e8189707 513 for(Int_t iz=0;iz<fNPixelsZ;iz++){
514 for(Int_t ix=0;ix<fNPixelsX;ix++){
b0f5e3fc 515 electronics = fBaseline + fNoise*random->Gaus();
516 signal = (float)fMapA2->GetSignal(iz,ix);
517 signal += electronics;
e8189707 518 gi =iz*fNPixelsX+ix; // global index
b0f5e3fc 519 if (signal > threshold) {
520 digits[0]=iz;
521 digits[1]=ix;
522 digits[2]=1;
b0f5e3fc 523 for(j1=0;j1<3;j1++){
e8189707 524 if (pList[gi]) {
525 tracks[j1] = (Int_t)(*(pList[gi]+j1));
526 hits[j1] = (Int_t)(*(pList[gi]+j1+6));
527 }else {
528 tracks[j1]=-2; //noise
529 hits[j1] = -1;
530 }
b0f5e3fc 531 charges[j1] = 0;
532 }
533 phys=0;
e8189707 534 aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
b0f5e3fc 535 }
e8189707 536 if(pList[gi]) delete [] pList[gi];
b0f5e3fc 537 }
538 }
539 delete [] pList;
540
b0f5e3fc 541}
542
543
544//____________________________________________
545
e8189707 546void AliITSsimulationSPD::CreateHistograms()
547{
548 // create 1D histograms for tests
549
550 printf("SPD - create histograms\n");
b0f5e3fc 551
e8189707 552 for (Int_t i=0;i<fNPixelsZ;i++) {
b0f5e3fc 553 TString *spdname = new TString("spd_");
e8189707 554 Char_t pixelz[4];
555 sprintf(pixelz,"%d",i+1);
556 spdname->Append(pixelz);
b0f5e3fc 557 (*fHis)[i] = new TH1F(spdname->Data(),"SPD maps",
558 fNPixelsX,0.,(Float_t) fNPixelsX);
559 delete spdname;
560 }
561
562}
563
564//____________________________________________
565
e8189707 566void AliITSsimulationSPD::ResetHistograms()
567{
b0f5e3fc 568 //
569 // Reset histograms for this detector
570 //
e8189707 571 for ( int i=0;i<fNPixelsZ;i++ ) {
b0f5e3fc 572 if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();
573 }
574
575}
e8189707 576
577
578
579
580
581
582
583
584