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