New ITS code for new structure and simulations.
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSPD.cxx
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
11 ClassImp(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
21 AliITSsimulationSPD::AliITSsimulationSPD(){
22   // constructor
23   fResponse     = 0;
24   fSegmentation = 0;
25   fHis          = 0;
26   fNoise        = 0.;
27   fBaseline     = 0.;
28 }
29 //_____________________________________________________________________________
30
31 AliITSsimulationSPD::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
49 AliITSsimulationSPD::~AliITSsimulationSPD() { 
50   // destructor
51
52   delete fMapA2;
53
54   if (fHis) {
55      fHis->Delete(); 
56      delete fHis;     
57   }                
58 }
59
60 //__________________________________________________________________________
61 AliITSsimulationSPD::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 //_________________________________________________________________________
74 AliITSsimulationSPD& 
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
88 void 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 //---------------------------------------------
393 void 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 //---------------------------------------------
472 void 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
514 void 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
532 void 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 }