Release version of ITS code
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSPD.cxx
1 #include <iostream.h>
2 #include <TRandom.h>
3 #include <TH1.h>
4 #include <TMath.h>
5
6
7
8 #include "AliRun.h"
9 #include "AliITS.h"
10 #include "AliITSMapA2.h" 
11 #include "AliITSsimulationSPD.h"
12
13
14
15 ClassImp(AliITSsimulationSPD)
16 ////////////////////////////////////////////////////////////////////////
17 // Version: 0
18 // Written by Boris Batyunya
19 // December 20 1999
20 //
21 // AliITSsimulationSPD is the simulation of SPDs
22 //________________________________________________________________________
23
24
25 AliITSsimulationSPD::AliITSsimulationSPD()
26 {
27   // constructor
28   fResponse = 0;
29   fSegmentation = 0;
30   fHis = 0;
31   fNoise=0.;
32   fBaseline=0.;
33 }
34
35
36 //_____________________________________________________________________________
37
38 AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg, AliITSresponse *resp) {
39   // standard constructor
40
41       fResponse = resp;
42       fSegmentation = seg;
43
44       fResponse->GetNoiseParam(fNoise,fBaseline);
45
46       fMapA2 = new AliITSMapA2(fSegmentation);
47
48       //
49       fNPixelsZ=fSegmentation->Npz();
50       fNPixelsX=fSegmentation->Npx();
51
52 }
53
54 //_____________________________________________________________________________
55
56 AliITSsimulationSPD::~AliITSsimulationSPD() { 
57   // destructor
58
59   delete fMapA2;
60
61   if (fHis) {
62      fHis->Delete(); 
63      delete fHis;     
64   }                
65 }
66
67
68 //__________________________________________________________________________
69 AliITSsimulationSPD::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 //_________________________________________________________________________
82 AliITSsimulationSPD& 
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
96 void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t module, Int_t dummy)
97 {
98   // digitize module
99
100     const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons 
101                                       // for 3.6 eV/pair 
102     const Float_t kconv = 10000.;     // cm -> microns
103
104     Float_t spdLength = fSegmentation->Dz();
105     Float_t spdWidth = fSegmentation->Dx();
106
107     Float_t difCoef, dum; 
108     fResponse->DiffCoeff(difCoef,dum); 
109
110     Float_t zPix0 = 1e+6;
111     Float_t xPix0 = 1e+6;
112     Float_t yPix0 = 1e+6;
113     Float_t yPrev = 1e+6;   
114     Float_t zP0 = 100.;
115     Float_t xP0 = 100.;
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
131     Int_t maxNDigits = fNPixelsX*fNPixelsZ + fNPixelsX ;; 
132     Float_t  **pList = new Float_t* [maxNDigits]; 
133     memset(pList,0,sizeof(Float_t*)*maxNDigits);
134     Int_t indexRange[4] = {0,0,0,0};
135
136     // Fill detector maps with GEANT hits
137     // loop over hits in the module
138     static Bool_t first=kTRUE;
139     Int_t lasttrack=-2;
140     Int_t hit, iZi, jz, jx;
141     for (hit=0;hit<nhits;hit++) {
142         AliITShit *iHit = (AliITShit*) fHits->At(hit);
143         Int_t layer = iHit->GetLayer();
144
145
146         // work with the idtrack=entry number in the TreeH
147         Int_t idhit,idtrack;
148         mod->GetHitTrackAndHitIndex(hit,idtrack,idhit);    
149         //Int_t idtrack=mod->GetHitTrackIndex(hit);  
150         // or store straight away the particle position in the array
151         // of particles : 
152         Int_t itrack = iHit->GetTrack();
153    
154         //  Get hit z and x(r*phi) cordinates for each module (detector)
155         //  in local system.
156
157         Float_t zPix = kconv*iHit->GetZL();
158         Float_t xPix = kconv*iHit->GetXL();
159         Float_t yPix = kconv*iHit->GetYL();
160
161         // Get track status
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.;  
180
181         Int_t 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              zP0 = zP;
193              xP0 = xP;
194         }
195         // enter Si only
196         if (layer == 2 && status == 66 && yPix < -71.) {    
197              yPix0 = yPix;
198              zP0 = zP;
199              xP0 = xP;
200         }       
201         Float_t depEnergy = iHit->GetIonization();
202         // skip if the input point to Si       
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; 
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;
215         } 
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
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;
233
234         if(TMath::Abs(ydif) < 0.1) continue; // Ydif is not zero
235
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; 
239
240         // number of the steps along the track:
241         Int_t nsteps = ndZ;
242         if(ndX > ndZ) nsteps = ndX;
243         if(nsteps < 6) nsteps = 6;  // minimum number of the steps 
244
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;
249
250         if (TMath::Abs(projDif) < 5.0 ) {
251            drPath = ydif*1.e-4;  
252            drPath = TMath::Abs(drPath);        // drift path in cm
253            sigmaDif = difCoef*sqrt(drPath);    // sigma diffusion in cm        
254         }  
255
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;
261
262             if(TMath::Abs(projDif) >= 5.) {
263               Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
264               if(trdown == 0) {
265                 drPath = dProjn*tang*1.e-4; // drift path for iZi step in cm 
266                 drPath = TMath::Abs(drPath);
267               }
268               if(trdown == 1) {
269                 Float_t dProjn = projDif/nsteps; 
270                 drPath = (projDif-(iZi-1)*dProjn)*tang*1.e-4;
271                 drPath = TMath::Abs(drPath);
272               }
273               sigmaDif = difCoef*sqrt(drPath);    
274               sigmaDif = sigmaDif*kconv;         // sigma diffusion in microns
275             }
276             zPixn = (zPixn + spdLength/2.);  
277             xPixn = (xPixn + spdWidth/2.);  
278             Int_t nZpix, nXpix;
279             fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
280             zPitch = fSegmentation->Dpz(nZpix);
281             fSegmentation->GetPadTxz(xPixn,zPixn);
282             // set the window for the integration
283             Int_t jzmin = 1;  
284             Int_t jzmax = 3; 
285             if(nZpix == 1) jzmin =2;
286             if(nZpix == fNPixelsZ) jzmax = 2; 
287
288             Int_t jxmin = 1;  
289             Int_t jxmax = 3; 
290             if(nXpix == 1) jxmin =2;
291             if(nXpix == fNPixelsX) jxmax = 2; 
292
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.;
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                 } 
319                 // kz changes from 1 to the fNofPixels(270)  
320                 Int_t kz = nZpix + jz -2; 
321
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; 
327
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                      } 
345                      Int_t kx = nXpix + jx -2;  
346
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; 
352
353                      if(dXCharge > 1.) {
354                        Int_t index = kz-1;
355
356                        if (first) {
357                           indexRange[0]=indexRange[1]=index;
358                           indexRange[2]=indexRange[3]=kx-1;
359                           first=kFALSE;
360                        }
361
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
367                        // build the list of digits for this module      
368                        Double_t signal=fMapA2->GetSignal(index,kx-1);
369                        signal+=dXCharge;
370                        fMapA2->SetHit(index,kx-1,(double)signal);
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
383         if (lasttrack != itrack || hit==(nhits-1)) {
384             GetList(itrack,idhit,pList,indexRange);
385             first=kTRUE;
386         }
387
388         lasttrack=itrack;
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 //---------------------------------------------
403 void AliITSsimulationSPD::GetList(Int_t label,Int_t idhit,Float_t **pList,Int_t *indexRange)
404 {
405   // lop over nonzero digits
406
407    
408   //set protection
409   for(int k=0;k<4;k++) {
410      if (indexRange[k] < 0) indexRange[k]=0;
411   }
412
413   for(Int_t iz=indexRange[0];iz<indexRange[1]+1;iz++){
414     for(Int_t ix=indexRange[2];ix<indexRange[3]+1;ix++){
415
416         Float_t signal=fMapA2->GetSignal(iz,ix);
417         if (!signal) continue;
418
419         Int_t globalIndex = iz*fNPixelsX+ix; // GlobalIndex starts from 0!
420         if(!pList[globalIndex]){
421         
422            // 
423            // Create new list (9 elements - 3 signals and 3 tracks + 3 hits)
424            //
425
426            pList[globalIndex] = new Float_t [9];
427
428            // set list to -3 
429
430            *pList[globalIndex] = -3.;
431            *(pList[globalIndex]+1) = -3.;
432            *(pList[globalIndex]+2) = -3.;
433            *(pList[globalIndex]+3) =  0.;
434            *(pList[globalIndex]+4) =  0.;
435            *(pList[globalIndex]+5) =  0.;
436            *(pList[globalIndex]+6) = -1.;
437            *(pList[globalIndex]+7) = -1.;
438            *(pList[globalIndex]+8) = -1.;
439
440
441            *pList[globalIndex] = (float)label;
442            *(pList[globalIndex]+3) = signal;
443            *(pList[globalIndex]+6) = (float)idhit;
444         }
445         else{
446
447           // check the signal magnitude
448
449           Float_t highest = *(pList[globalIndex]+3);
450           Float_t middle = *(pList[globalIndex]+4);
451           Float_t lowest = *(pList[globalIndex]+5);
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;
469
470             *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
471             *(pList[globalIndex]+7) = *(pList[globalIndex]+6);
472             *(pList[globalIndex]+6) = idhit;
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;
480
481             *(pList[globalIndex]+8) = *(pList[globalIndex]+7);
482             *(pList[globalIndex]+7) = idhit;
483           }
484           else{
485             *(pList[globalIndex]+5) = signal;
486             *(pList[globalIndex]+2) = label;
487             *(pList[globalIndex]+8) = idhit;
488           }
489         }
490     } // end of loop pixels in x
491   } // end of loop over pixels in z
492
493
494 }
495
496
497 //---------------------------------------------
498 void AliITSsimulationSPD::ChargeToSignal(Float_t **pList)
499 {
500   // add noise and electronics, perform the zero suppression and add the
501   // digit to the list
502
503   AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
504   
505
506   TRandom *random = new TRandom(); 
507   Float_t threshold = (float)fResponse->MinVal();
508
509   Int_t digits[3], tracks[3], hits[3],gi,j1;
510   Float_t charges[3];
511   Float_t electronics;
512   Float_t signal,phys;
513   for(Int_t iz=0;iz<fNPixelsZ;iz++){
514     for(Int_t ix=0;ix<fNPixelsX;ix++){
515       electronics = fBaseline + fNoise*random->Gaus();
516       signal = (float)fMapA2->GetSignal(iz,ix);
517       signal += electronics;
518       gi =iz*fNPixelsX+ix; // global index
519       if (signal > threshold) {
520          digits[0]=iz;
521          digits[1]=ix;
522          digits[2]=1;
523          for(j1=0;j1<3;j1++){
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            }
531            charges[j1] = 0;
532          }
533          phys=0;
534          aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
535       }
536       if(pList[gi]) delete [] pList[gi];
537     }
538   }
539   delete [] pList;
540
541 }
542
543
544 //____________________________________________
545
546 void AliITSsimulationSPD::CreateHistograms()
547 {
548   // create 1D histograms for tests
549
550       printf("SPD - create histograms\n");
551
552       for (Int_t i=0;i<fNPixelsZ;i++) {
553            TString *spdname = new TString("spd_");
554            Char_t pixelz[4];
555            sprintf(pixelz,"%d",i+1);
556            spdname->Append(pixelz);
557            (*fHis)[i] = new TH1F(spdname->Data(),"SPD maps",
558                               fNPixelsX,0.,(Float_t) fNPixelsX);
559            delete spdname;
560       }
561
562 }
563
564 //____________________________________________
565
566 void AliITSsimulationSPD::ResetHistograms()
567 {
568     //
569     // Reset histograms for this detector
570     //
571     for ( int i=0;i<fNPixelsZ;i++ ) {
572         if ((*fHis)[i])    ((TH1F*)(*fHis)[i])->Reset();
573     }
574
575 }
576
577
578
579
580
581
582
583
584