SPD Bari model added
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSPDbari.cxx
CommitLineData
0aedb77e 1#include <iostream.h>
2#include <TRandom.h>
3#include <TH1.h>
4#include <TMath.h>
5#include <TString.h>
6#include <TParticle.h>
7
8
9#include "AliRun.h"
10#include "AliITS.h"
11#include "AliITShit.h"
12#include "AliITSdigit.h"
13#include "AliITSmodule.h"
14#include "AliITSMapA2.h"
15#include "AliITSsimulationSPDbari.h"
16#include "AliITSsegmentation.h"
17#include "AliITSresponse.h"
18
19
20ClassImp(AliITSsimulationSPDbari)
21////////////////////////////////////////////////////////////////////////
22// Version: 0
23// Written by Rocco Caliandro
24// from a model developed with T. Virgili and R.A. Fini
25// June 15 2000
26//
27// AliITSsimulationSPD is the simulation of SPDs
28//
29//________________________________________________________________________
30
31AliITSsimulationSPDbari::AliITSsimulationSPDbari(){
32 // constructor
33 fResponse = 0;
34 fSegmentation = 0;
35 fHis = 0;
36 fThresh = 0.;
37 fSigma = 0.;
38 fCouplCol = 0.;
39 fCouplRow = 0.;
40}
41//_____________________________________________________________________________
42
43AliITSsimulationSPDbari::AliITSsimulationSPDbari(AliITSsegmentation *seg, AliITSresponse *resp) {
44 // constructor
45 fResponse = resp;
46 fSegmentation = seg;
47
48 fResponse->Thresholds(fThresh,fSigma);
49 fResponse->GetNoiseParam(fCouplCol,fCouplRow);
50
51 fMapA2 = new AliITSMapA2(fSegmentation);
52
53 //
54 fNPixelsZ=fSegmentation->Npz();
55 fNPixelsX=fSegmentation->Npx();
56
57}
58
59//_____________________________________________________________________________
60
61AliITSsimulationSPDbari::~AliITSsimulationSPDbari() {
62 // destructor
63
64 delete fMapA2;
65
66 if (fHis) {
67 fHis->Delete();
68 delete fHis;
69 }
70}
71
72//__________________________________________________________________________
73AliITSsimulationSPDbari::AliITSsimulationSPDbari(const AliITSsimulationSPDbari &source){
74 // Copy Constructor
75 if(&source == this) return;
76 this->fMapA2 = source.fMapA2;
77 this->fThresh = source.fThresh;
78 this->fSigma = source.fSigma;
79 this->fCouplCol = source.fCouplCol;
80 this->fCouplRow = source.fCouplRow;
81 this->fNPixelsX = source.fNPixelsX;
82 this->fNPixelsZ = source.fNPixelsZ;
83 this->fHis = source.fHis;
84 return;
85}
86
87//_________________________________________________________________________
88AliITSsimulationSPDbari&
89 AliITSsimulationSPDbari::operator=(const AliITSsimulationSPDbari &source) {
90 // Assignment operator
91 if(&source == this) return *this;
92 this->fMapA2 = source.fMapA2;
93 this->fThresh = source.fThresh;
94 this->fSigma = source.fSigma;
95 this->fCouplCol = source.fCouplCol;
96 this->fCouplRow = source.fCouplRow;
97 this->fNPixelsX = source.fNPixelsX;
98 this->fNPixelsZ = source.fNPixelsZ;
99 this->fHis = source.fHis;
100 return *this;
101 }
102//_____________________________________________________________________________
103
104void AliITSsimulationSPDbari::DigitiseModule(AliITSmodule *mod, Int_t module,
105 Int_t dummy) {
106 // digitize module
107
108
109 TObjArray *fHits = mod->GetHits();
110 Int_t nhits = fHits->GetEntriesFast();
111 if (!nhits) return;
112
113
114 printf("Module %d (%d hits) \n",module+1,nhits);
115
116
117 Int_t number=10000;
118 Int_t *frowpixel = new Int_t[number];
119 Int_t *fcolpixel = new Int_t[number];
120 Double_t *fenepixel = new Double_t[number];
121
122 // Array of pointers to store the track index of the digits
123 // leave +1, otherwise pList crashes when col=256, row=192
124 Int_t maxNdigits = fNPixelsX*fNPixelsZ+fNPixelsX+1;
125 Float_t **pList = new Float_t* [maxNdigits];
126 memset(pList,0,sizeof(Float_t*)*maxNdigits);
127
128
129 // noise setting
130 SetFluctuations(pList);
131
132
133 // loop over hits in the module
134 Int_t hitpos;
135 for (hitpos=0;hitpos<nhits;hitpos++) {
136
137 HitToDigit(mod,hitpos,module,frowpixel,fcolpixel,fenepixel,pList);
138
139
140 }// end loop over digits
141
142 CreateDigit(nhits,module,pList);
143
144 // clean memory
145 delete[] frowpixel;
146 delete[] fcolpixel;
147 delete[] fenepixel;
148 fMapA2->ClearMap();
149 delete [] pList;
150
151}
152//_____________________________________________________________________________
153
154void AliITSsimulationSPDbari::UpdateMap( Int_t row, Int_t col, Double_t ene) {
155//
156// updates the Map of signal, adding the energy (ene) released by the current track
157//
158 Double_t signal;
159 signal = fMapA2->GetSignal(row,col);
160 signal += ene;
161 fMapA2->SetHit(row,col,signal);
162
163 }
164//_____________________________________________________________________________
165void AliITSsimulationSPDbari::HitToDigit(AliITSmodule *mod, Int_t hitpos, Int_t module,
166 Int_t *frowpixel, Int_t *fcolpixel,
167 Double_t *fenepixel, Float_t **pList) {
168//
169// Steering function to determine the digits associated to a given hit (hitpos)
170// The digits are created by charge sharing (ChargeSharing) and by
171// capacitive coupling (SetCoupling). At all the created digits is associated
172// the track number of the hit (ntrack)
173//
174
175
176 static Float_t x1l,y1l,z1l;
177 Float_t x2l,y2l,z2l,etot;
178 Int_t layer,r1,r2,c1,c2,row,col,npixel = 0;
179 Int_t ntrack;
180 Double_t ene;
181 const Float_t kconv = 10000.; // cm -> microns
182
183 TObjArray *fHits = mod->GetHits();
184 AliITShit *hit = (AliITShit*) fHits->At(hitpos);
185 layer = hit->GetLayer();
186 etot=hit->GetIonization();
187 ntrack=hit->GetTrack();
188
189
190 if (hit->GetTrackStatus()==66) {
191 hit->GetPositionL(x1l,y1l,z1l);
192 // positions shifted and converted in microns
193 x1l = x1l*kconv + fSegmentation->Dx()/2.;
194 z1l = z1l*kconv + fSegmentation->Dz()/2.;
195 }
196 else {
197 hit->GetPositionL(x2l,y2l,z2l);
198 // positions shifted and converted in microns
199 x2l = x2l*kconv + fSegmentation->Dx()/2.;
200 z2l = z2l*kconv + fSegmentation->Dz()/2.;
201
202
203 // to account for 83750 effective sensitive area
204 // introduced in geometry (Dz=83600)
205 if (z1l>fSegmentation->Dz()) return;
206 if (z2l>fSegmentation->Dz()) return;
207 // to preserve for hit having x>12800
208 if (x1l>fSegmentation->Dx()) return;
209 if (x2l>fSegmentation->Dx()) return;
210
211 //Get the col and row number starting from 1
212 // the x direction is not inverted for the second layer!!!
213 fSegmentation->GetPadIxz(x1l, z1l, c1, r1);
214 fSegmentation->GetPadIxz(x2l, z2l, c2, r2);
215
216
217
218 if ((r1==r2) && (c1==c2))
219 {
220 // no charge sharing
221 npixel = 1;
222 frowpixel[npixel-1] = r1;
223 fcolpixel[npixel-1] = c1;
224 fenepixel[npixel-1] = etot;
225 }
226 else {
227 // charge sharing
228 ChargeSharing(x1l,z1l,x2l,z2l,c1,r1,c2,r2,etot,
229 npixel,frowpixel,fcolpixel,fenepixel);
230
231 }
232
233
234 for (Int_t npix=0;npix<npixel;npix++)
235 {
236 row = frowpixel[npix];
237 col = fcolpixel[npix];
238 ene = fenepixel[npix];
239 UpdateMap(row,col,ene);
240 GetList(ntrack,pList,row,col);
241 // Starting capacitive coupling effect
242 SetCoupling(row,col,ntrack,pList);
243 }
244 x1l=x2l;
245 y1l=y2l;
246 z1l=z2l;
247 }
248}
249
250//_________________________________________________________________________
251
252void AliITSsimulationSPDbari::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
253 Float_t z2l,Int_t c1,Int_t r1,Int_t c2,
254 Int_t r2,Float_t etot,
255 Int_t &npixel,Int_t *frowpixel,
256 Int_t *fcolpixel,Double_t *fenepixel){
257 //
258 // Take into account the geometrical charge sharing when the track
259 // crosses more than one pixel.
260 //
261 //Begin_Html
262 /*
263 <img src="model_2.gif">
264 </pre>
265 <br clear=left>
266 <font size=+2 color=red>
267 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
268 </font>
269 <pre>
270 */
271 //End_Html
272
273
274 Float_t xa,za,xb,zb,dx,dz,dtot,dm,refr,refm,refc;
275 Float_t refn=0.;
276 Float_t arefm, arefr, arefn, arefc, azb, az2l, axb, ax2l;
277 Int_t dirx,dirz,rb,cb;
278
279
280 Int_t flag,flagrow,flagcol;
281
282 Double_t epar;
283
284
285 npixel = 0;
286 xa = x1l;
287 za = z1l;
288 dx = TMath::Abs(x1l-x2l);
289 dz = TMath::Abs(z1l-z2l);
290 dtot = TMath::Sqrt((dx*dx)+(dz*dz));
291 dm = (x2l - x1l) / (z2l - z1l);
292
293 dirx = (Int_t) ((x2l - x1l) / dx);
294 dirz = (Int_t) ((z2l - z1l) / dz);
295
296
297 // calculate the x coordinate of the pixel in the next column
298 // and the z coordinate of the pixel in the next row
299
300 Float_t xpos, zpos;
301
302 fSegmentation->GetPadCxz(c1, r1-1, xpos, zpos);
303
304 Float_t xsize = fSegmentation->Dpx(0);
305 Float_t zsize = fSegmentation->Dpz(r1);
306
307 if (dirx == 1) refr = xpos+xsize/2.;
308 else refr = xpos-xsize/2.;
309
310 if (dirz == 1) refn = zpos+zsize/2.;
311 else refn = zpos-zsize/2.;
312
313
314 flag = 0;
315 flagrow = 0;
316 flagcol = 0;
317 do
318 {
319
320 // calculate the x coordinate of the intersection with the pixel
321 // in the next cell in row direction
322
323 refm = (refn - z1l)*dm + x1l;
324
325 // calculate the z coordinate of the intersection with the pixel
326 // in the next cell in column direction
327
328 refc = (refr - x1l)/dm + z1l;
329
330
331 arefm = refm * dirx;
332 arefr = refr * dirx;
333 arefn = refn * dirz;
334 arefc = refc * dirz;
335
336
337 if ((arefm < arefr) && (arefn < arefc)){
338
339 // the track goes in the pixel in the next cell in row direction
340 xb = refm;
341 zb = refn;
342 cb = c1;
343 rb = r1 + dirz;
344 azb = zb * dirz;
345 az2l = z2l * dirz;
346 if (rb == r2) flagrow=1;
347 if (azb > az2l) {
348 zb = z2l;
349 xb = x2l;
350 }
351
352 // shift to the pixel in the next cell in row direction
353 Float_t zsizeNext = fSegmentation->Dpz(rb);
354 refn += zsizeNext*dirz;
355
356 }
357 else {
358
359 // the track goes in the pixel in the next cell in column direction
360 xb = refr;
361 zb = refc;
362 cb = c1 + dirx;
363 rb = r1;
364 axb = xb * dirx;
365 ax2l = x2l * dirx;
366 if (cb == c2) flagcol=1;
367 if (axb > ax2l) {
368 zb = z2l;
369 xb = x2l;
370 }
371
372 // shift to the pixel in the next cell in column direction
373 Float_t xsizeNext = fSegmentation->Dpx(cb);
374 refr += xsizeNext*dirx;
375
376 }
377
378 //calculate the energy lost in the crossed pixel
379 epar = TMath::Sqrt((xb-xa)*(xb-xa)+(zb-za)*(zb-za));
380 epar = etot*(epar/dtot);
381
382 //store row, column and energy lost in the crossed pixel
383 frowpixel[npixel] = r1;
384 fcolpixel[npixel] = c1;
385 fenepixel[npixel] = epar;
386 npixel++;
387
388 // the exit point of the track is reached
389 if (epar == 0) flag = 1;
390 if ((r1 == r2) && (c1 == c2)) flag = 1;
391 if (flag!=1) {
392 r1 = rb;
393 c1 = cb;
394 xa = xb;
395 za = zb;
396 }
397
398 } while (flag==0);
399
400}
401//___________________________________________________________________________
402void AliITSsimulationSPDbari::SetCoupling(Int_t row, Int_t col, Int_t ntrack,
403 Float_t **pList) {
404 //
405 // Take into account the coupling between adiacent pixels.
406 // The parameters probcol and probrow are the fractions of the
407 // signal in one pixel shared in the two adjacent pixels along
408 // the column and row direction, respectively.
409 //
410 //Begin_Html
411 /*
412 <img src="model_3.gif">
413 </pre>
414 <br clear=left>
415 <font size=+2 color=red>
416 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
417 </font>
418 <pre>
419 */
420 //End_Html
421
422
423 Int_t j1,j2,flag=0;
424 Double_t pulse1,pulse2;
425
426
427 j1 = row;
428 j2 = col;
429
430 pulse1 = fMapA2->GetSignal(row,col);
431 pulse2 = pulse1;
432
433 for (Int_t isign=-1;isign<=1;isign+=2)
434 {
435
436// loop in row direction
437
438 do
439 {
440 j1 += isign;
441 pulse1 *= fCouplRow;
442
443 if ((j1 < 0) || (j1 > fNPixelsZ-1) || (pulse1 < fThresh))
444 {
445 pulse1 = fMapA2->GetSignal(row,col);
446 j1 = row;
447 flag = 1;
448 }
449 else{
450 UpdateMap(j1,col,pulse1);
451 GetList(ntrack,pList,j1,col);
452 flag = 0;
453 }
454
455 } while(flag == 0);
456
457
458// loop in column direction
459
460 do
461 {
462 j2 += isign;
463 pulse2 *= fCouplCol;
464
465 if ((j2 < 0) || (j2 > (fNPixelsX-1)) || (pulse2 < fThresh))
466 {
467 pulse2 = fMapA2->GetSignal(row,col);
468 j2 = col;
469 flag = 1;
470 }
471 else{
472 UpdateMap(row,j2,pulse2);
473 GetList(ntrack,pList,row,j2);
474 flag = 0;
475 }
476
477 } while(flag == 0);
478
479 }
480
481}
482//___________________________________________________________________________
483void AliITSsimulationSPDbari::CreateDigit(Int_t nhits, Int_t module, Float_t
484**pList) {
485 //
486 // The pixels are fired if the energy deposited inside them is above
487 // the threshold parameter ethr. Fired pixed are interpreted as digits
488 // and stored in the file digitfilename.
489 //
490
491 AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
492
493
494 Int_t digits[3];
495 Int_t tracks[3];
496 Int_t hits[3];
497 Float_t charges[3];
498 Int_t gi,j1;
499
500
501 if (nhits > 0) {
502
503 for (Int_t r=1;r<=fNPixelsZ;r++) {
504 for (Int_t c=1;c<=fNPixelsX;c++) {
505
506 // check if the deposited energy in a pixel is above the threshold
507 Float_t signal = (Float_t) fMapA2->GetSignal(r,c);
508 if ( signal > fThresh) {
509 digits[0] = r-1; // digits starts from 0
510 digits[1] = c-1; // digits starts from 0
511 digits[2] = 1;
512 gi =r*fNPixelsX+c; // global index
513 for(j1=0;j1<3;j1++){
514 tracks[j1] = (Int_t)(*(pList[gi]+j1));
515 charges[j1] = 0;
516 }
517 Float_t phys = 0;
518 aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
519 if(pList[gi]) delete [] pList[gi];
520 }//endif of threshold condition
521 }
522 }// enddo on pixels
523 }
524
525}
526//_____________________________________________________________________________
527
528void AliITSsimulationSPDbari::GetList(Int_t label,Float_t **pList,
529 Int_t row, Int_t col) {
530 // loop over nonzero digits
531
532 Int_t ix = col;
533 Int_t iz = row;
534 Int_t globalIndex;
535 Float_t signal;
536 Float_t highest,middle,lowest;
537
538
539 signal=fMapA2->GetSignal(iz,ix);
540
541 globalIndex = iz*fNPixelsX+ix; // globalIndex starts from 1
542
543
544 if(!pList[globalIndex])
545 {
546 //
547 // Create new list (6 elements - 3 signals and 3 tracks + total sig)
548 //
549
550 pList[globalIndex] = new Float_t [6];
551
552
553 // set list to -2
554 *(pList[globalIndex]) = -2.;
555 *(pList[globalIndex]+1) = -2.;
556 *(pList[globalIndex]+2) = -2.;
557 *(pList[globalIndex]+3) = 0.;
558 *(pList[globalIndex]+4) = 0.;
559 *(pList[globalIndex]+5) = 0.;
560
561 *pList[globalIndex] = (float)label;
562 *(pList[globalIndex]+3) = signal;
563 }
564 else{
565
566
567 // check the signal magnitude
568 highest = *(pList[globalIndex]+3);
569 middle = *(pList[globalIndex]+4);
570 lowest = *(pList[globalIndex]+5);
571
572
573 signal -= (highest+middle+lowest);
574
575
576 //
577 // compare the new signal with already existing list
578 //
579 if(signal<lowest) return; // neglect this track
580
581 if (signal>highest)
582 {
583 *(pList[globalIndex]+5) = middle;
584 *(pList[globalIndex]+4) = highest;
585 *(pList[globalIndex]+3) = signal;
586 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
587 *(pList[globalIndex]+1) = *pList[globalIndex];
588 *(pList[globalIndex]) = label;
589 }
590 else if (signal>middle)
591 {
592 *(pList[globalIndex]+5) = middle;
593 *(pList[globalIndex]+4) = signal;
594 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
595 *(pList[globalIndex]+1) = label;
596 }
597 else
598 {
599 *(pList[globalIndex]+5) = signal;
600 *(pList[globalIndex]+2) = label;
601 }
602 }
603}
604//_________________________________________________________________________
605void AliITSsimulationSPDbari::SetFluctuations(Float_t **pList) {
606 //
607 // Set the electronic noise and threshold non-uniformities to all the
608 // pixels in a detector.
609 // The parameter fSigma is the squared sum of the sigma due to noise
610 // and the sigma of the threshold distribution among pixels.
611 //
612 //Begin_Html
613 /*
614 <img src="model_1.gif">
615 </pre>
616 <br clear=left>
617 <font size=+2 color=red>
618 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
619 </font>
620 <pre>
621 */
622 //End_Html
623
624
625 TRandom *random = new TRandom();
626 Double_t signal;
627
628 Int_t iz,ix;
629 for(iz=1;iz<=fNPixelsZ;iz++){
630 for(ix=1;ix<=fNPixelsX;ix++){
631 signal = fSigma*random->Gaus();
632 fMapA2->SetHit(iz,ix,signal);
633
634 // insert in the label-signal list the pixels fired only by noise
635 if ( signal > fThresh) {
636 Int_t globalIndex = iz*fNPixelsX+ix;
637 pList[globalIndex] = new Float_t [6];
638 *(pList[globalIndex]) = -2.;
639 *(pList[globalIndex]+1) = -2.;
640 *(pList[globalIndex]+2) = -2.;
641 *(pList[globalIndex]+3) = signal;
642 *(pList[globalIndex]+4) = 0.;
643 *(pList[globalIndex]+5) = 0.;
644 }
645 } // end of loop on pixels
646 } // end of loop on pixels
647
648 }
649//____________________________________________
650
651void AliITSsimulationSPDbari::CreateHistograms() {
652 // CreateHistograms
653
654 Int_t i;
655 for(i=0;i<fNPixelsZ;i++) {
656 TString *spdname = new TString("spd_");
657 Char_t candnum[4];
658 sprintf(candnum,"%d",i+1);
659 spdname->Append(candnum);
660 (*fHis)[i] = new TH1F(spdname->Data(),"SPD maps",
661 fNPixelsX,0.,(Float_t) fNPixelsX);
662 delete spdname;
663 }
664
665}
666
667//____________________________________________
668
669void AliITSsimulationSPDbari::ResetHistograms() {
670 //
671 // Reset histograms for this detector
672 //
673 Int_t i;
674 for(i=0;i<fNPixelsZ;i++ ) {
675 if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();
676 }
677
678}