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