Update of library for low mass resonances (Yiota Foka)
[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
a82b8484 216 // to account for unexpected equal entrance and
217 // exit coordinates
218 if (x1l==x2l) x2l=x2l+x2l*0.000001;
219 if (z1l==z2l) z2l=z2l+z2l*0.000001;
220
221 // to account for tracks at the edge of the sensitive area
222 // of SPDs
223 if (x1l<0) return;
224 if (x2l<0) return;
225 if (z1l<0) return;
226 if (z2l<0) return;
0aedb77e 227
a82b8484 228
0aedb77e 229 if ((r1==r2) && (c1==c2))
230 {
231 // no charge sharing
232 npixel = 1;
233 frowpixel[npixel-1] = r1;
234 fcolpixel[npixel-1] = c1;
235 fenepixel[npixel-1] = etot;
236 }
237 else {
238 // charge sharing
239 ChargeSharing(x1l,z1l,x2l,z2l,c1,r1,c2,r2,etot,
240 npixel,frowpixel,fcolpixel,fenepixel);
241
242 }
243
244
245 for (Int_t npix=0;npix<npixel;npix++)
246 {
247 row = frowpixel[npix];
248 col = fcolpixel[npix];
249 ene = fenepixel[npix];
250 UpdateMap(row,col,ene);
251 GetList(ntrack,pList,row,col);
252 // Starting capacitive coupling effect
253 SetCoupling(row,col,ntrack,pList);
254 }
255 x1l=x2l;
256 y1l=y2l;
257 z1l=z2l;
258 }
259}
260
261//_________________________________________________________________________
262
263void AliITSsimulationSPDbari::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
264 Float_t z2l,Int_t c1,Int_t r1,Int_t c2,
265 Int_t r2,Float_t etot,
266 Int_t &npixel,Int_t *frowpixel,
267 Int_t *fcolpixel,Double_t *fenepixel){
268 //
269 // Take into account the geometrical charge sharing when the track
270 // crosses more than one pixel.
271 //
272 //Begin_Html
273 /*
4b365332 274 <img src="picts/ITS/barimodel_2.gif">
0aedb77e 275 </pre>
276 <br clear=left>
277 <font size=+2 color=red>
278 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
279 </font>
280 <pre>
281 */
282 //End_Html
283
284
285 Float_t xa,za,xb,zb,dx,dz,dtot,dm,refr,refm,refc;
286 Float_t refn=0.;
287 Float_t arefm, arefr, arefn, arefc, azb, az2l, axb, ax2l;
288 Int_t dirx,dirz,rb,cb;
289
290
291 Int_t flag,flagrow,flagcol;
292
293 Double_t epar;
294
295
296 npixel = 0;
297 xa = x1l;
298 za = z1l;
299 dx = TMath::Abs(x1l-x2l);
300 dz = TMath::Abs(z1l-z2l);
301 dtot = TMath::Sqrt((dx*dx)+(dz*dz));
302 dm = (x2l - x1l) / (z2l - z1l);
303
304 dirx = (Int_t) ((x2l - x1l) / dx);
305 dirz = (Int_t) ((z2l - z1l) / dz);
306
307
308 // calculate the x coordinate of the pixel in the next column
309 // and the z coordinate of the pixel in the next row
310
311 Float_t xpos, zpos;
312
313 fSegmentation->GetPadCxz(c1, r1-1, xpos, zpos);
314
315 Float_t xsize = fSegmentation->Dpx(0);
316 Float_t zsize = fSegmentation->Dpz(r1);
317
318 if (dirx == 1) refr = xpos+xsize/2.;
319 else refr = xpos-xsize/2.;
320
321 if (dirz == 1) refn = zpos+zsize/2.;
322 else refn = zpos-zsize/2.;
323
324
325 flag = 0;
326 flagrow = 0;
327 flagcol = 0;
328 do
329 {
330
331 // calculate the x coordinate of the intersection with the pixel
332 // in the next cell in row direction
333
334 refm = (refn - z1l)*dm + x1l;
335
336 // calculate the z coordinate of the intersection with the pixel
337 // in the next cell in column direction
338
339 refc = (refr - x1l)/dm + z1l;
340
341
342 arefm = refm * dirx;
343 arefr = refr * dirx;
344 arefn = refn * dirz;
345 arefc = refc * dirz;
346
347
348 if ((arefm < arefr) && (arefn < arefc)){
349
350 // the track goes in the pixel in the next cell in row direction
351 xb = refm;
352 zb = refn;
353 cb = c1;
354 rb = r1 + dirz;
355 azb = zb * dirz;
356 az2l = z2l * dirz;
357 if (rb == r2) flagrow=1;
358 if (azb > az2l) {
359 zb = z2l;
360 xb = x2l;
361 }
362
363 // shift to the pixel in the next cell in row direction
364 Float_t zsizeNext = fSegmentation->Dpz(rb);
365 refn += zsizeNext*dirz;
366
367 }
368 else {
369
370 // the track goes in the pixel in the next cell in column direction
371 xb = refr;
372 zb = refc;
373 cb = c1 + dirx;
374 rb = r1;
375 axb = xb * dirx;
376 ax2l = x2l * dirx;
377 if (cb == c2) flagcol=1;
378 if (axb > ax2l) {
379 zb = z2l;
380 xb = x2l;
381 }
382
383 // shift to the pixel in the next cell in column direction
384 Float_t xsizeNext = fSegmentation->Dpx(cb);
385 refr += xsizeNext*dirx;
386
387 }
388
389 //calculate the energy lost in the crossed pixel
390 epar = TMath::Sqrt((xb-xa)*(xb-xa)+(zb-za)*(zb-za));
391 epar = etot*(epar/dtot);
392
393 //store row, column and energy lost in the crossed pixel
a82b8484 394
395
0aedb77e 396 frowpixel[npixel] = r1;
397 fcolpixel[npixel] = c1;
398 fenepixel[npixel] = epar;
399 npixel++;
400
401 // the exit point of the track is reached
402 if (epar == 0) flag = 1;
403 if ((r1 == r2) && (c1 == c2)) flag = 1;
404 if (flag!=1) {
405 r1 = rb;
406 c1 = cb;
407 xa = xb;
408 za = zb;
409 }
410
411 } while (flag==0);
412
413}
414//___________________________________________________________________________
415void AliITSsimulationSPDbari::SetCoupling(Int_t row, Int_t col, Int_t ntrack,
416 Float_t **pList) {
417 //
418 // Take into account the coupling between adiacent pixels.
419 // The parameters probcol and probrow are the fractions of the
420 // signal in one pixel shared in the two adjacent pixels along
421 // the column and row direction, respectively.
422 //
423 //Begin_Html
424 /*
4b365332 425 <img src="picts/ITS/barimodel_3.gif">
0aedb77e 426 </pre>
427 <br clear=left>
428 <font size=+2 color=red>
429 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
430 </font>
431 <pre>
432 */
433 //End_Html
434
435
436 Int_t j1,j2,flag=0;
437 Double_t pulse1,pulse2;
438
439
440 j1 = row;
441 j2 = col;
442
443 pulse1 = fMapA2->GetSignal(row,col);
444 pulse2 = pulse1;
445
446 for (Int_t isign=-1;isign<=1;isign+=2)
447 {
448
449// loop in row direction
450
451 do
452 {
453 j1 += isign;
454 pulse1 *= fCouplRow;
455
456 if ((j1 < 0) || (j1 > fNPixelsZ-1) || (pulse1 < fThresh))
457 {
458 pulse1 = fMapA2->GetSignal(row,col);
459 j1 = row;
460 flag = 1;
461 }
462 else{
463 UpdateMap(j1,col,pulse1);
464 GetList(ntrack,pList,j1,col);
465 flag = 0;
466 }
467
468 } while(flag == 0);
469
470
471// loop in column direction
472
473 do
474 {
475 j2 += isign;
476 pulse2 *= fCouplCol;
477
478 if ((j2 < 0) || (j2 > (fNPixelsX-1)) || (pulse2 < fThresh))
479 {
480 pulse2 = fMapA2->GetSignal(row,col);
481 j2 = col;
482 flag = 1;
483 }
484 else{
485 UpdateMap(row,j2,pulse2);
486 GetList(ntrack,pList,row,j2);
487 flag = 0;
488 }
489
490 } while(flag == 0);
491
492 }
493
494}
495//___________________________________________________________________________
496void AliITSsimulationSPDbari::CreateDigit(Int_t nhits, Int_t module, Float_t
497**pList) {
498 //
499 // The pixels are fired if the energy deposited inside them is above
500 // the threshold parameter ethr. Fired pixed are interpreted as digits
501 // and stored in the file digitfilename.
502 //
503
504 AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
505
506
507 Int_t digits[3];
508 Int_t tracks[3];
509 Int_t hits[3];
510 Float_t charges[3];
511 Int_t gi,j1;
512
513
514 if (nhits > 0) {
515
516 for (Int_t r=1;r<=fNPixelsZ;r++) {
517 for (Int_t c=1;c<=fNPixelsX;c++) {
518
519 // check if the deposited energy in a pixel is above the threshold
520 Float_t signal = (Float_t) fMapA2->GetSignal(r,c);
521 if ( signal > fThresh) {
522 digits[0] = r-1; // digits starts from 0
523 digits[1] = c-1; // digits starts from 0
524 digits[2] = 1;
525 gi =r*fNPixelsX+c; // global index
526 for(j1=0;j1<3;j1++){
527 tracks[j1] = (Int_t)(*(pList[gi]+j1));
528 charges[j1] = 0;
529 }
530 Float_t phys = 0;
531 aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
532 if(pList[gi]) delete [] pList[gi];
533 }//endif of threshold condition
534 }
535 }// enddo on pixels
536 }
537
538}
539//_____________________________________________________________________________
540
541void AliITSsimulationSPDbari::GetList(Int_t label,Float_t **pList,
542 Int_t row, Int_t col) {
543 // loop over nonzero digits
544
545 Int_t ix = col;
546 Int_t iz = row;
547 Int_t globalIndex;
548 Float_t signal;
549 Float_t highest,middle,lowest;
550
551
552 signal=fMapA2->GetSignal(iz,ix);
553
554 globalIndex = iz*fNPixelsX+ix; // globalIndex starts from 1
555
556
557 if(!pList[globalIndex])
558 {
559 //
560 // Create new list (6 elements - 3 signals and 3 tracks + total sig)
561 //
562
563 pList[globalIndex] = new Float_t [6];
564
565
566 // set list to -2
567 *(pList[globalIndex]) = -2.;
568 *(pList[globalIndex]+1) = -2.;
569 *(pList[globalIndex]+2) = -2.;
570 *(pList[globalIndex]+3) = 0.;
571 *(pList[globalIndex]+4) = 0.;
572 *(pList[globalIndex]+5) = 0.;
573
574 *pList[globalIndex] = (float)label;
575 *(pList[globalIndex]+3) = signal;
576 }
577 else{
578
579
580 // check the signal magnitude
581 highest = *(pList[globalIndex]+3);
582 middle = *(pList[globalIndex]+4);
583 lowest = *(pList[globalIndex]+5);
584
585
586 signal -= (highest+middle+lowest);
587
588
589 //
590 // compare the new signal with already existing list
591 //
592 if(signal<lowest) return; // neglect this track
593
594 if (signal>highest)
595 {
596 *(pList[globalIndex]+5) = middle;
597 *(pList[globalIndex]+4) = highest;
598 *(pList[globalIndex]+3) = signal;
599 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
600 *(pList[globalIndex]+1) = *pList[globalIndex];
601 *(pList[globalIndex]) = label;
602 }
603 else if (signal>middle)
604 {
605 *(pList[globalIndex]+5) = middle;
606 *(pList[globalIndex]+4) = signal;
607 *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
608 *(pList[globalIndex]+1) = label;
609 }
610 else
611 {
612 *(pList[globalIndex]+5) = signal;
613 *(pList[globalIndex]+2) = label;
614 }
615 }
616}
617//_________________________________________________________________________
618void AliITSsimulationSPDbari::SetFluctuations(Float_t **pList) {
619 //
620 // Set the electronic noise and threshold non-uniformities to all the
621 // pixels in a detector.
622 // The parameter fSigma is the squared sum of the sigma due to noise
623 // and the sigma of the threshold distribution among pixels.
624 //
625 //Begin_Html
626 /*
4b365332 627 <img src="picts/ITS/barimodel_1.gif">
0aedb77e 628 </pre>
629 <br clear=left>
630 <font size=+2 color=red>
631 <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
632 </font>
633 <pre>
634 */
635 //End_Html
636
637
8813cd63 638 TRandom random;
0aedb77e 639 Double_t signal;
640
641 Int_t iz,ix;
642 for(iz=1;iz<=fNPixelsZ;iz++){
643 for(ix=1;ix<=fNPixelsX;ix++){
8813cd63 644 signal = fSigma*random.Gaus();
0aedb77e 645 fMapA2->SetHit(iz,ix,signal);
646
647 // insert in the label-signal list the pixels fired only by noise
648 if ( signal > fThresh) {
649 Int_t globalIndex = iz*fNPixelsX+ix;
650 pList[globalIndex] = new Float_t [6];
651 *(pList[globalIndex]) = -2.;
652 *(pList[globalIndex]+1) = -2.;
653 *(pList[globalIndex]+2) = -2.;
654 *(pList[globalIndex]+3) = signal;
655 *(pList[globalIndex]+4) = 0.;
656 *(pList[globalIndex]+5) = 0.;
657 }
658 } // end of loop on pixels
659 } // end of loop on pixels
660
661 }
662//____________________________________________
663
664void AliITSsimulationSPDbari::CreateHistograms() {
665 // CreateHistograms
666
667 Int_t i;
668 for(i=0;i<fNPixelsZ;i++) {
8813cd63 669 TString spdname("spd_");
0aedb77e 670 Char_t candnum[4];
671 sprintf(candnum,"%d",i+1);
8813cd63 672 spdname.Append(candnum);
673 (*fHis)[i] = new TH1F(spdname.Data(),"SPD maps",
0aedb77e 674 fNPixelsX,0.,(Float_t) fNPixelsX);
0aedb77e 675 }
676
677}
678
679//____________________________________________
680
681void AliITSsimulationSPDbari::ResetHistograms() {
682 //
683 // Reset histograms for this detector
684 //
685 Int_t i;
686 for(i=0;i<fNPixelsZ;i++ ) {
687 if ((*fHis)[i]) ((TH1F*)(*fHis)[i])->Reset();
688 }
689
690}