dee880f3789564bbedae6757e8d0576221176a98
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSSD.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 /* $Id$ */
16
17 #include <stdio.h>
18 #include <stdlib.h>
19 #include <iostream.h>
20 #include <iomanip.h>
21 #include <TObjArray.h>
22 #include <TRandom.h>
23 #include <TMath.h>
24
25 #include "AliITSmodule.h"
26 #include "AliITSMapA2.h"
27 #include "AliITSsegmentationSSD.h"
28 #include "AliITSresponseSSD.h"
29 #include "AliITSdcsSSD.h"
30 #include "AliITS.h"
31 #include "AliRun.h"
32 #include "AliITSgeom.h"
33 #include "AliITSsimulationSSD.h"
34
35 ClassImp(AliITSsimulationSSD);
36 ////////////////////////////////////////////////////////////////////////
37 // Version: 0
38 // Written by Enrico Fragiacomo
39 // July 2000
40 //
41 // AliITSsimulationSSD is the simulation of SSDs.
42
43 //----------------------------------------------------------------------
44 AliITSsimulationSSD::AliITSsimulationSSD(){
45     //default Constructor
46
47     fDCS     = 0;
48     fNstrips = GetSegmentation()->Npx();
49     fPitch   = GetSegmentation()->Dpx(0);
50     fDifConst[0] = fDifConst[1] = 0.0;
51     fDriftVel[0] = fDriftVel[1] = 0.0;
52     fMapA2   = 0;
53 }
54 //----------------------------------------------------------------------
55 AliITSsimulationSSD::AliITSsimulationSSD(AliITSsegmentation *seg,
56                                          AliITSresponse *resp){
57     // Constructor
58
59     fSegmentation    = seg;
60     fResponse        = resp;
61     Float_t noise[2] = {0.,0.};
62     fResponse->GetNoiseParam(noise[0],noise[1]); // retrieves noise parameters
63     fDCS             = new AliITSdcsSSD(seg,resp); 
64
65     fNstrips         = GetSegmentation()->Npx();
66     fPitch           = GetSegmentation()->Dpx(0);
67     SetDriftVelocity(); // use default values in .h file
68     SetIonizeE();       // use default values in .h file
69     SetDiffConst();     // use default values in .h file
70     fMapA2           = new AliITSMapA2(fSegmentation);
71 }
72 //______________________________________________________________________
73 AliITSsimulationSSD& AliITSsimulationSSD::operator=(
74                                                 const AliITSsimulationSSD &s){
75     // Operator =
76
77     if(this==&s) return *this;
78
79     this->fDCS         = new AliITSdcsSSD(*(s.fDCS));
80     this->fMapA2       = s.fMapA2;
81     this->fNstrips     = s.fNstrips;
82     this->fPitch       = s.fPitch;
83     this->fIonE        = s.fIonE;
84     this->fDifConst[0] = s.fDifConst[0];
85     this->fDifConst[1] = s.fDifConst[1];
86     this->fDriftVel[0] = s.fDriftVel[0];
87     this->fDriftVel[1] = s.fDriftVel[1];
88     return *this;
89 }
90 //______________________________________________________________________
91 AliITSsimulationSSD::AliITSsimulationSSD(const AliITSsimulationSSD &source){
92     // copy constructor
93
94     *this = source;
95 }
96 //______________________________________________________________________
97 AliITSsimulationSSD::~AliITSsimulationSSD() {
98     // destructor
99     delete fMapA2;
100     delete fDCS;
101 }
102 //______________________________________________________________________
103 void AliITSsimulationSSD::DigitiseModule(AliITSmodule *mod,Int_t module,
104                                          Int_t dummy) {
105     // Digitizes hits for one SSD module 
106
107     Int_t lay, lad, detect;
108     AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
109     AliITSgeom *geom = aliITS->GetITSgeom();
110     geom->GetModuleId(module,lay, lad, detect);
111     if ( lay == 6 ) GetSegmentation()->SetLayer(6);
112     if ( lay == 5 ) GetSegmentation()->SetLayer(5);
113
114     TObjArray *hits = mod->GetHits();
115     Int_t nhits     = hits->GetEntriesFast();
116     if (!nhits) return;
117   
118     Double_t x0=0.0, y0=0.0, z0=0.0;
119     Double_t x1=0.0, y1=0.0, z1=0.0;
120     Double_t de=0.0;
121     Int_t maxNdigits = 2*fNstrips;
122     Float_t  **pList = new Float_t* [maxNdigits]; 
123     memset(pList,0,sizeof(Float_t*)*maxNdigits);
124     Int_t indexRange[4] = {0,0,0,0};
125     static Bool_t first = kTRUE;
126     Int_t lasttrack     = -2;
127     Int_t idtrack       = -2;
128     
129     for(Int_t i=0; i<nhits; i++) {    
130         // LineSegmentL returns 0 if the hit is entering
131         // If hits is exiting returns positions of entering and exiting hits
132         // Returns also energy loss
133
134         if (mod->LineSegmentL(i, x0, x1, y0, y1, z0, z1, de, idtrack)) {
135             HitToDigit(module, x0, y0, z0, x1, y1, z1, de, indexRange, first);
136        
137             if (lasttrack != idtrack || i==(nhits-1)) {
138                 GetList(idtrack,pList,indexRange);
139                 first=kTRUE;
140             } // end if
141             lasttrack=idtrack;
142         } // end if
143     }  // end loop over hits
144   
145     ApplyNoise();
146     ApplyCoupling();
147
148     ChargeToSignal(pList);
149
150     fMapA2->ClearMap();
151 }
152 //----------------------------------------------------------------------
153 void AliITSsimulationSSD::HitToDigit(Int_t module, Double_t x0, Double_t y0, 
154                                      Double_t z0, Double_t x1, Double_t y1, 
155                                      Double_t z1, Double_t de,
156                                      Int_t *indexRange, Bool_t first) {
157     // Turns hits in SSD module into one or more digits.
158
159     Float_t tang[2] = {0.0,0.0};
160     GetSegmentation()->Angles(tang[0], tang[1]);//stereo<<->tan(stereo)~=stereo
161     Double_t x, y, z;
162     Double_t dex=0.0, dey=0.0, dez=0.0;
163     Double_t pairs; // pair generation energy per step.
164     Double_t sigma[2] = {0.,0.};// standard deviation of the diffusion gaussian
165     Double_t tdrift[2] = {0.,0.}; // time of drift
166     Double_t w;
167     Double_t inf[2], sup[2], par0[2];                 
168
169     // Steps in the module are determined "manually" (i.e. No Geant)
170     // NumOfSteps divide path between entering and exiting hits in steps 
171     Int_t numOfSteps = NumOfSteps(x1, y1, z1, dex, dey, dez);
172   
173     // Enery loss is equally distributed among steps
174     de    = de/numOfSteps;
175     pairs = de/GetIonizeE(); // e-h pairs generated
176
177     for(Int_t j=0; j<numOfSteps; j++) {     // stepping
178         x = x0 + (j+0.5)*dex;
179         y = y0 + (j+0.5)*dey;
180         if ( y > (GetSegmentation()->Dy()/2+10)*1.0E-4 ) {
181             // check if particle is within the detector
182             cout<<"AliITSsimulationSSD::HitToDigit: Warning: hit "
183                 "out of detector y0,y,dey,j ="
184                 <<y0<<","<<y<<","<<dey<<","<<j<<endl;
185             return;
186         } // end if
187         z = z0 + (j+0.5)*dez;
188
189         // calculate drift time
190         // y is the minimum path
191         tdrift[0] = (y+(GetSegmentation()->Dy()*1.0E-4)/2)/GetDriftVelocity(0);
192         tdrift[1] = ((GetSegmentation()->Dy()*1.0E-4)/2-y)/GetDriftVelocity(1);
193
194         for(Int_t k=0; k<2; k++) {   // both sides    remember: 0=Pside 1=Nside
195
196             tang[k]=TMath::Tan(tang[k]);
197
198             // w is the coord. perpendicular to the strips
199             if(k==0) {
200                 w = (x+(GetSegmentation()->Dx()*1.0E-4)/2) -
201                     (z+(GetSegmentation()->Dz()*1.0E-4)/2)*tang[k]; 
202             }else{
203                 w = (x+(GetSegmentation()->Dx()*1.0E-4)/2) + 
204                     (z-(GetSegmentation()->Dz()*1.0E-4)/2)*tang[k];
205             } // end if
206             w = w / (fPitch*1.0E-4); // w is converted in units of pitch
207
208             if((w<(-0.5)) || (w>(fNstrips-0.5))) {
209                 // this check rejects hits in regions not covered by strips
210                 // 0.5 takes into account boundaries 
211                 if(k==0) cout<<"AliITSsimulationSSD::HitToDigit: "
212                              "Warning: no strip in this region of P side"
213                              <<endl;
214                 else cout<<"AliITSsimulationSSD::HitToDigit: "
215                          "Warning: no strip in this region of N side"<<endl;
216                 return;
217             } // end if
218
219             // sigma is the standard deviation of the diffusion gaussian
220             if(tdrift[k]<0) return;
221             sigma[k] = TMath::Sqrt(2*GetDiffConst(k)*tdrift[k]);
222             sigma[k] = sigma[k] /(fPitch*1.0E-4);  //units of Pitch
223             if(sigma[k]==0.0) {         
224                 cout<<"AliITSsimulationSSD::DigitiseModule: Error: sigma=0"
225                     <<endl; 
226                 exit(0);
227             } // end if
228
229             par0[k] = pairs;
230             // we integrate the diffusion gaussian from -3sigma to 3sigma 
231             inf[k] = w - 3*sigma[k]; // 3 sigma from the gaussian average  
232             sup[k] = w + 3*sigma[k]; // 3 sigma from the gaussian average
233                                      // IntegrateGaussian does the actual
234                                      // integration of diffusion gaussian
235             IntegrateGaussian(k, par0[k], w, sigma[k], inf[k], sup[k], 
236                               indexRange, first);
237         }  // end for loop over side (0=Pside, 1=Nside)      
238     } // end stepping
239     //delete seg;
240 }
241 //______________________________________________________________________
242 void AliITSsimulationSSD::ApplyNoise() {
243     // Apply Noise.
244   
245     Double_t signal;
246     Double_t noise[2] = {0.,0.};
247     Float_t a,b;
248     fResponse->GetNoiseParam(a,b); // retrieves noise parameters
249     noise[0] = (Double_t) a; noise[1] = (Double_t) b;
250     for(Int_t k=0;k<2;k++){                    // both sides (0=Pside, 1=Nside)
251         for(Int_t ix=0;ix<fNstrips;ix++){      // loop over strips
252             signal = fMapA2->GetSignal(k,ix); // retrieves signal from map
253             signal += gRandom->Gaus(0,noise[k]);// add noise to signal
254             if(signal<0.) signal=0.0;           // in case noise is negative...
255             fMapA2->SetHit(k,ix,signal); // give back signal to map
256         } // loop over strip 
257     } // loop over k (P or N side)
258 }
259 //______________________________________________________________________
260 void AliITSsimulationSSD::ApplyCoupling() {
261     // Apply the effect of electronic coupling between channels
262     Double_t signalLeft=0, signalRight=0;
263
264     for(Int_t ix=0;ix<fNstrips;ix++){
265         if(ix>0.)signalLeft = fMapA2->GetSignal(0,ix-1)*fDCS->GetCouplingPL();
266         else signalLeft = 0.0;
267         if(ix<(fNstrips-1)) signalRight = fMapA2->GetSignal(0,ix+1)*
268                                                         fDCS->GetCouplingPR();
269         else signalRight = 0.0;
270         fMapA2->AddSignal(0,ix,signalLeft + signalRight);
271     
272         if(ix>0.) signalLeft = fMapA2->GetSignal(1,ix-1)*fDCS->GetCouplingNL();
273         else signalLeft = 0.0;
274         if(ix<(fNstrips-1)) signalRight = fMapA2->GetSignal(1,ix+1)*
275                                                          fDCS->GetCouplingNR();
276         else signalRight = 0.0;
277         fMapA2->AddSignal(1,ix,signalLeft + signalRight);
278     } // loop over strips 
279 }
280 //______________________________________________________________________
281 Float_t AliITSsimulationSSD::F(Float_t av, Float_t x, Float_t s) {
282     // Computes the integral of a gaussian using Error Function
283     Float_t sqrt2 = TMath::Sqrt(2.0);
284     Float_t sigm2 = sqrt2*s;
285     Float_t integral;
286
287     integral = 0.5 * TMath::Erf( (x - av) / sigm2);
288     return integral;
289 }
290 //______________________________________________________________________
291 void AliITSsimulationSSD::IntegrateGaussian(Int_t k,Double_t par, Double_t w,
292                                             Double_t sigma, 
293                                             Double_t inf, Double_t sup,
294                                             Int_t *indexRange, Bool_t first) {
295     // integrate the diffusion gaussian
296     // remind: inf and sup are w-3sigma and w+3sigma
297     //         we could define them here instead of passing them
298     //         this way we are free to introduce asimmetry
299
300     Double_t a=0.0, b=0.0;
301     Double_t dXCharge1 = 0.0, dXCharge2 = 0.0;
302     // dXCharge1 and 2 are the charge to two neighbouring strips
303     // Watch that we only involve at least two strips
304     // Numbers greater than 2 of strips in a cluster depend on
305     //  geometry of the track and delta rays, not charge diffusion!   
306   
307     Double_t strip = TMath::Floor(w);         // clostest strip on the left
308
309     if ( TMath::Abs((strip - w)) < 0.5) { 
310         // gaussian mean is closer to strip on the left
311         a = inf;                         // integration starting point
312         if((strip+0.5)<=sup) {
313             // this means that the tail of the gaussian goes beyond
314             // the middle point between strips ---> part of the signal
315             // is given to the strip on the right
316             b = strip + 0.5;               // integration stopping point
317             dXCharge1 = F( w, b, sigma) - F(w, a, sigma);
318             dXCharge2 = F( w, sup, sigma) - F(w ,b, sigma); 
319         }else { 
320             // this means that all the charge is given to the strip on the left
321             b = sup;
322             dXCharge1 = 0.9973;   // gaussian integral at 3 sigmas
323             dXCharge2 = 0.0;
324         } // end if
325
326         dXCharge1 = par * dXCharge1;// normalize by mean of number of carriers
327         dXCharge2 = par * dXCharge2;
328
329         // for the time being, signal is the charge
330         // in ChargeToSignal signal is converted in ADC channel
331         fMapA2->AddSignal(k,strip,dXCharge1);
332         if(((Int_t) strip) < (fNstrips-1)) {
333             // strip doesn't have to be the last (remind: last=fNstrips-1)
334             // otherwise part of the charge is lost
335             fMapA2->AddSignal(k,(strip+1),dXCharge2);
336         } // end if
337     
338         if(dXCharge1 > 1.) {
339             if (first) {
340                 indexRange[k*2+0] = indexRange[k*2+1]=(Int_t) strip;
341                 first=kFALSE;
342             } // end if first
343
344             indexRange[k*2+0]=TMath::Min(indexRange[k*2+0],(Int_t) strip);
345             indexRange[k*2+1]=TMath::Max(indexRange[k*2+1],(Int_t) strip);
346         }      // dXCharge > 1 e-
347
348     }else{
349         // gaussian mean is closer to strip on the right
350         strip++;     // move to strip on the rigth
351         b = sup;     // now you know where to stop integrating
352         if((strip-0.5)>=inf) { 
353             // tail of diffusion gaussian on the left goes left of
354             // middle point between strips
355             a = strip - 0.5;        // integration starting point
356             dXCharge1 = F(w, b, sigma) - F(w, a, sigma);
357             dXCharge2 = F(w, a, sigma) - F(w, inf, sigma);
358         }else {
359             a = inf;
360             dXCharge1 = 0.9973;   // gaussian integral at 3 sigmas
361             dXCharge2 = 0.0;
362         } // end if
363     
364         dXCharge1 = par * dXCharge1;    // normalize by means of carriers
365         dXCharge2 = par * dXCharge2;
366
367         // for the time being, signal is the charge
368         // in ChargeToSignal signal is converted in ADC channel
369         fMapA2->AddSignal(k,strip,dXCharge1);
370         if(((Int_t) strip) > 0) {
371             // strip doesn't have to be the first
372             // otherwise part of the charge is lost
373             fMapA2->AddSignal(k,(strip-1),dXCharge2);
374         } // end if
375     
376         if(dXCharge1 > 1.) {
377             if (first) {
378                 indexRange[k*2+0]=indexRange[k*2+1]=(Int_t) strip;
379                 first=kFALSE;
380             } // end if first
381
382             indexRange[k*2+0]=TMath::Min(indexRange[k*2+0],(Int_t) strip);
383             indexRange[k*2+1]=TMath::Max(indexRange[k*2+1],(Int_t) strip);
384         }      // dXCharge > 1 e-
385     } // end if
386 }
387 //______________________________________________________________________
388 Int_t AliITSsimulationSSD::NumOfSteps(Double_t x, Double_t y, Double_t z,
389                                  Double_t & dex,Double_t & dey,Double_t & dez){
390     // number of steps
391     // it also returns steps for each coord
392     //AliITSsegmentationSSD *seg = new AliITSsegmentationSSD();
393
394     Double_t step = 25E-4;
395     //step = (Double_t) seg->GetStepSize();  // step size (cm)
396     Int_t numOfSteps = (Int_t) (TMath::Sqrt(x*x+y*y+z*z)/step); 
397
398     if (numOfSteps < 1) numOfSteps = 1;       // one step, at least
399
400     // we could condition the stepping depending on the incident angle
401     // of the track
402     dex = x/numOfSteps;
403     dey = y/numOfSteps;
404     dez = z/numOfSteps;
405
406     return numOfSteps;
407 }
408 //----------------------------------------------------------------------
409 void AliITSsimulationSSD::GetList(Int_t label,Float_t **pList,
410                                   Int_t *indexRange) {
411     // loop over nonzero digits
412     Int_t ix,globalIndex;
413     Float_t signal=0.;
414     Float_t highest,middle,lowest;
415
416     for(Int_t k=0; k<2; k++) {
417         for(ix=indexRange[k*2+0];ix<indexRange[k*2+1]+1;ix++){
418             if(indexRange[k*2+0]<indexRange[k*2+1]) 
419                 signal = fMapA2->GetSignal(k,ix);
420
421             globalIndex = k*fNstrips+ix; // globalIndex starts from 0!
422             if(!pList[globalIndex]){
423                 // 
424                 //Create new list (6 elements-3 signals and 3 tracks+total sig)
425                 //
426                 pList[globalIndex] = new Float_t [6];
427                 // set list to -1
428                 *pList[globalIndex] = -2.;
429                 *(pList[globalIndex]+1) = -2.;
430                 *(pList[globalIndex]+2) = -2.;
431                 *(pList[globalIndex]+3) =  0.;
432                 *(pList[globalIndex]+4) =  0.;
433                 *(pList[globalIndex]+5) =  0.;
434                 *pList[globalIndex] = (float)label;
435                 *(pList[globalIndex]+3) = signal;
436             }else{
437                 // check the signal magnitude
438                 highest = *(pList[globalIndex]+3);
439                 middle = *(pList[globalIndex]+4);
440                 lowest = *(pList[globalIndex]+5);
441                 signal -= (highest+middle+lowest);
442                 //
443                 //  compare the new signal with already existing list
444                 //
445                 if(signal<lowest) continue; // neglect this track
446                 if (signal>highest){
447                     *(pList[globalIndex]+5) = middle;
448                     *(pList[globalIndex]+4) = highest;
449                     *(pList[globalIndex]+3) = signal;
450                     *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
451                     *(pList[globalIndex]+1) = *pList[globalIndex];
452                     *pList[globalIndex] = label;
453                 }else if (signal>middle){
454                     *(pList[globalIndex]+5) = middle;
455                     *(pList[globalIndex]+4) = signal;
456                     *(pList[globalIndex]+2) = *(pList[globalIndex]+1);
457                     *(pList[globalIndex]+1) = label;
458                 }else{
459                     *(pList[globalIndex]+5) = signal;
460                     *(pList[globalIndex]+2) = label;
461                 } // end if
462             } // end if
463         } // end of loop pixels in x
464     } // end of loop over pixels in z
465 }
466 //----------------------------------------------------------------------
467 void AliITSsimulationSSD::ChargeToSignal(Float_t **pList) {
468     // charge to signal
469     AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
470     Float_t threshold = 0.;
471     Int_t digits[3], tracks[3],hits[3],gi,j1;
472     Float_t charges[3];
473     Float_t signal,phys;
474     Float_t noise[2] = {0.,0.};
475
476     fResponse->GetNoiseParam(noise[0],noise[1]);
477   
478     for(Int_t k=0;k<2;k++){         // both sides (0=Pside, 1=Nside)
479
480         // Threshold for zero-suppression
481         // It can be defined in AliITSresponseSSD
482         //             threshold = (Float_t)fResponse->MinVal(k);
483         // I prefer to think adjusting the threshold "manually", looking
484         // at the scope, and considering noise standard deviation
485         threshold = 4.0*noise[k];      // 4 times noise is a choice
486         for(Int_t ix=0;ix<fNstrips;ix++){         // loop over strips
487
488             signal = (Float_t) fMapA2->GetSignal(k,ix);
489             gi =k*fNstrips+ix; // global index
490             if (signal > threshold) {
491                 digits[0]=k;
492                 digits[1]=ix;
493
494                 // convert to ADC signal
495                 // conversion factor are rather arbitrary (need tuning)
496                 // minimum ionizing particle--> ~30000 pairs--> ADC channel 50
497                 signal = signal*50.0/30000.0;        
498                 if(signal>1000.) signal = 1000.0;//if exceeding, accumulate
499                                                  // last one
500                 digits[2]=(Int_t) signal;
501                 for(j1=0;j1<3;j1++){
502                     if (pList[gi]) {
503                         tracks[j1] = (Int_t)(*(pList[gi]+j1));
504                     } else {
505                         tracks[j1]=-2; //noise
506                     } // end if pList
507                     charges[j1] = 0;
508                 } // end for j1
509
510                 phys=0;
511
512                 hits[0]=0;
513                 hits[1]=0;
514                 hits[2]=0;
515                 // finally add digit
516                 aliITS->AddSimDigit(2,phys,digits,tracks,hits,charges);
517             } // end if signal > threshold
518             if(pList[gi]) delete [] pList[gi];
519         } // end for ix
520     } // end for k
521     delete [] pList;
522 }
523 //______________________________________________________________________
524 void AliITSsimulationSSD::Print(ostream *os){
525     //Standard output format for this class
526
527     //AliITSsimulation::Print(os);
528     *os << fNstrips <<","<< fPitch <<","<< fIonE <<",";
529     *os << fDifConst[0] <<","<< fDifConst[1] <<",";
530     *os << fDriftVel[0] <<","<< fDriftVel[1];
531     //*os <<","; fDCS->Print(os);
532     //*os <<","; fMapA2->Print(os);
533 }
534 //______________________________________________________________________
535 void AliITSsimulationSSD::Read(istream *is){
536     // Standard output streaming function.
537
538     //AliITSsimulation::Read(is);
539     *is >> fNstrips >> fPitch >> fIonE;
540     *is >> fDifConst[0] >> fDifConst[1];
541     *is >> fDriftVel[0] >> fDriftVel[1];
542     //fDCS->Read(is);
543     //fMapA2->Read(is);
544 }
545 //______________________________________________________________________
546 ostream &operator<<(ostream &os,AliITSsimulationSSD &source){
547     // Standard output streaming function.
548
549     source.Print(&os);
550     return os;
551 }
552 //______________________________________________________________________
553 istream &operator>>(istream &os,AliITSsimulationSSD &source){
554     // Standard output streaming function.
555
556     source.Read(&os);
557     return os;
558 }