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