Divide by zero error fix in function F. Sigma difusion = 0 now set to same
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSSD.cxx
1
2 #include <stdio.h>
3 #include <TObjArray.h>
4
5 #include "AliITSsimulationSSD.h"
6 #include "AliITSdictSSD.h"
7 #include "AliITSdcsSSD.h"
8 #include "AliITS.h"
9 #include "AliRun.h"
10
11
12 ClassImp(AliITSsimulationSSD);
13 //------------------------------------------------------------
14 AliITSsimulationSSD::AliITSsimulationSSD(AliITSsegmentation *seg,
15                                          AliITSresponse *resp){
16   // Constructor
17
18     fSegmentation = seg;
19     fResponse = resp;
20     fDCS = new AliITSdcsSSD(seg,resp); 
21
22     fNstrips = fSegmentation->Npx();
23     fPitch = fSegmentation->Dpx(0);
24     
25     fP = new TArrayF(fNstrips+1); 
26     fN = new TArrayF(fNstrips+1);
27      
28     fTracksP = new AliITSdictSSD[fNstrips+1];
29     fTracksN = new AliITSdictSSD[fNstrips+1];
30
31     
32     fSteps  = 10;   // still hard-wired - set in SetDetParam and get it via  
33                      // fDCS together with the others eventually    
34
35 }
36 //___________________________________________________________________________
37 AliITSsimulationSSD& AliITSsimulationSSD::operator=(AliITSsimulationSSD 
38                                                                       &source){
39 // Operator =
40     if(this==&source) return *this;
41
42     this->fDCS = new AliITSdcsSSD(*(source.fDCS));
43     this->fN   = new TArrayF(*(source.fN));
44     this->fP   = new TArrayF(*(source.fP));
45     this->fTracksP = new AliITSdictSSD(*(source.fTracksP));
46     this->fTracksN = new AliITSdictSSD(*(source.fTracksN));
47     this->fNstrips = source.fNstrips;
48     this->fPitch   = source.fPitch;
49     this->fSteps   = source.fSteps;
50     return *this;
51 }
52 //_____________________________________________________________
53 AliITSsimulationSSD::AliITSsimulationSSD(AliITSsimulationSSD &source){
54   // copy constructor
55    *this = source;
56 }
57 //____________________________________________________________________________
58 AliITSsimulationSSD::~AliITSsimulationSSD() {
59   // anihilator    
60
61     if(fP) delete fP;
62     if(fN) delete fN;
63     
64     if(fTracksP) delete fTracksP;
65     if(fTracksN) delete fTracksN;
66
67     delete fDCS;
68
69 //_______________________________________________________________
70 //
71 // Hit to digit
72 //_______________________________________________________________
73 //
74 void AliITSsimulationSSD::DigitiseModule(AliITSmodule *mod,Int_t module,
75                                          Int_t dummy) {
76   // Digitizes one SSD module of hits.
77     
78     TObjArray *hits = mod->GetHits();
79     Int_t nhits = hits->GetEntriesFast();
80     if (!nhits) return;
81
82     Int_t i;
83     for(i=0; i<fNstrips; i++) {
84        (*fP)[i] = 0;
85        (*fN)[i] = 0;
86        fTracksP[i].ZeroTracks();
87        fTracksN[i].ZeroTracks();
88     }
89     
90     for(i=0; i<nhits; i++) {
91        Int_t idtrack=mod->GetHitTrackIndex(i);  
92        HitToDigit(i,idtrack,nhits,hits);
93     }
94  
95    
96
97     ApplyNoise();
98     ApplyCoupling();    
99     ApplyThreshold();
100     ApplyDAQ();
101       
102
103 }
104
105 //---------------------------------------------------------------
106
107 void AliITSsimulationSSD::HitToDigit(Int_t & hitNo,Int_t idtrack,
108                                      Int_t nhits,TObjArray *hits) {
109   // Turns one or more hits in an SSD module into one or more digits.
110      
111     Int_t      stripP, stripN, i;
112     Float_t    dsP, dsN;
113     Float_t    sP, sN;
114     Float_t    eP, eN;
115     Float_t    arrayEP[10];         // hard-wired number of steps
116     Float_t    arrayEN[10];
117     Int_t      track = -1;
118        
119     Float_t    ionization = 0;
120     Float_t    signal;
121     
122     AliITSdictSSD *dict;  
123     
124  
125     // check if this is the right order !!!!!
126
127     AliITShit *hitI = (AliITShit*)hits->At(hitNo++);
128     AliITShit *hitE = (AliITShit*)hits->At(hitNo);
129
130
131     while (!((hitE->StatusExiting()) || 
132              (hitE->StatusDisappeared()) ||
133              (hitE->StatusStop()))) {
134        
135         if (++hitNo<nhits) {
136            ionization = hitE->GetIonization(); 
137            hitE = (AliITShit*)hits->At(hitNo);
138         }
139     }   
140     
141         
142     if (hitI->GetTrack() == hitE->GetTrack()) 
143        //track = idtrack;
144        track = hitI->GetTrack();
145     else 
146        printf("!!! Emergency !!!\n");
147      
148             
149     ionization += hitE->GetIonization();
150        
151     const Float_t kconvm=10000.;  // cm -> microns
152
153     Float_t xI, yI, zI;
154     hitI->GetPositionL(xI, yI, zI);
155     
156     xI *= kconvm;
157     yI *= kconvm;
158     zI *= kconvm;
159     
160     Float_t xE, yE, zE;
161     hitE->GetPositionL(xE, yE, zE);
162     
163     xE *= kconvm;
164     yE *= kconvm;
165     zE *= kconvm;
166
167     Float_t dx = (xE - xI);
168     Float_t dz = (zE - zI);
169               
170     
171     // Debuging
172     /*
173     fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
174    
175        printf("%5d %8.3f %8.3f %8.3f %8.3f %d %d  %d\n", 
176              hitNo, xI, zI, dx, dz, 
177              stripP, stripN, track);
178      printf("%10.5f %10d \n", ionization, hitI->fTrack); 
179     */ 
180     
181     // end of debuging   
182     
183     
184     eP=0;
185     eN=0;
186     
187     for (i=0; i<fSteps; i++) {
188         
189       //        arrayEP[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
190       //        arrayEN[i] = gRandom->Landau(ionization/fSteps, ionization/(4*fSteps));
191         arrayEP[i] = ionization/fSteps;
192         arrayEN[i] = ionization/fSteps;
193        
194         eP += arrayEP[i];
195         eN += arrayEN[i];
196     } 
197        
198     const Float_t kconv = 1.0e9 / 3.6;  // GeV -> e-hole pairs
199        
200     for(i=0; i<fSteps; i++) {
201     
202         arrayEP[i] = kconv * arrayEP[i] * (ionization / eP); 
203         arrayEN[i] = kconv * arrayEN[i] * (ionization / eN);        
204     }    
205         
206     dx /= fSteps;
207     dz /= fSteps;  
208
209     Float_t sigmaP, sigmaN; 
210     fResponse->SigmaSpread(sigmaP,sigmaN); 
211
212     //printf("SigmaP SigmaN %f %f\n",sigmaP, sigmaN);
213
214     Float_t noiseP, noiseN;
215     fResponse->GetNoiseParam(noiseP,noiseN);
216
217     //printf("NoiseP NoiseN %f %f\n",noiseP, noiseN);
218
219     for(i=0; i<fSteps; i++) {
220         
221         Int_t j;
222        
223         fSegmentation->GetPadIxz(xI,zI,stripP,stripN);
224         //printf("hitNo %d i xI zI stripP stripN %d %f %f %d %d\n",hitNo,i,xI, zI, stripP, stripN);
225         dsP    = Get2Strip(1,stripP,xI, zI); // Between 0-1
226         dsN    = Get2Strip(0,stripN,xI, zI); // Between 0-1
227
228         sP = sigmaP * sqrt(300. * i / (fSteps));
229         if(sP<=0.0) sP = sigmaP*sqrt(300.);
230         sN = sigmaN * sqrt(300. * i /(fSteps-i));
231         if(sN<=0.0) sN = sigmaN*sqrt(300.);
232
233
234         sP = (i<2        && dsP>0.3 && dsP<0.7)? 20. : sP;  // square of (microns) 
235         sN = (i>fSteps-2 && dsN>0.3 && dsN<0.7)? 20. : sN;  // square of (microns) 
236
237         sP = (i==2 && dsP>0.4 && dsP<0.6)? 15. : sP;  // square of (microns) 
238         sN = (i==8 && dsN>0.4 && dsN<0.6)? 15. : sN;  // square of (microns)        
239         
240         for (j=-1; j<2; j++) {
241             if (stripP+j<0 || stripP+j>fNstrips) continue;
242             signal = arrayEP[i] * TMath::Abs( (F(j+0.5-dsP,sP)-F(j-0.5-dsP,sP)) );
243             //printf("SimSSD::HitsToDigits:%d arrayEP[%d]=%e signal=%e\n",j,i,arrayEP[i],signal);
244             if (signal > noiseP/fSteps) {
245                (*fP)[stripP+j] += signal;
246                dict = (fTracksP+stripP+j);   
247                (*dict).AddTrack(track);
248             } 
249         }  // end for j loop over neighboring strips
250
251         for (j=-1; j<2; j++) {
252             if (stripN+j<0 || stripN+j>fNstrips) continue;
253             signal = arrayEN[i] * TMath::Abs( (F(j+0.5-dsN,sN)-F(j-0.5-dsN,sN)) );
254             //printf("SimSSD::HitsToDigits:%d arrayEN[%d]=%e signal=%e\n",j,i,arrayEN[i],signal);
255             if (signal > noiseN/fSteps) {
256                (*fN)[stripN+j] += signal;
257                dict = (fTracksN+stripN+j);    //co to jest
258                (*dict).AddTrack(track);
259             } 
260         }  // end for j loop over neighboring strips
261                 
262         xI += dx; 
263         zI += dz; 
264     }
265     
266     
267 }
268
269
270 //____________________________________________________________________
271 //
272 //  Private Methods for Simulation
273 //______________________________________________________________________
274 //
275
276 void AliITSsimulationSSD::ApplyNoise() {
277   // Apply Noise.
278    Float_t noiseP, noiseN;
279    fResponse->GetNoiseParam(noiseP,noiseN);
280        
281     Int_t i;
282     for(i = 0; i<fNstrips; i++) {
283        (*fP)[i] += gRandom->Gaus(0,noiseP);
284        (*fN)[i] += gRandom->Gaus(0,noiseN);
285     }
286 }
287
288 //_________________________________________________________________________
289
290 void AliITSsimulationSSD::ApplyCoupling() {
291   // Apply the effecto of electronic coupling between channels    
292     Int_t i;
293     for(i = 1; i<fNstrips-1; i++) {
294       (*fP)[i] += (*fP)[i-1]*fDCS->GetCouplingPL() + (*fP)[i+1]*fDCS->GetCouplingPR();
295       (*fN)[i] += (*fN)[i-1]*fDCS->GetCouplingNL() + (*fN)[i+1]*fDCS->GetCouplingNR();
296     }
297 }
298
299 //__________________________________________________________________________
300
301 void AliITSsimulationSSD::ApplyThreshold() {
302   // Applies the effect of a threshold on the signals for digitization.
303    Float_t noiseP, noiseN;
304    fResponse->GetNoiseParam(noiseP,noiseN);
305
306    // or introduce the SetThresholds in fResponse  
307
308     Int_t i;
309     for(i=0; i<fNstrips; i++) {
310        (*fP)[i] = ((*fP)[i] > noiseP*4) ? (*fP)[i] : 0;
311        (*fN)[i] = ((*fN)[i] > noiseN*4) ? (*fN)[i] : 0; 
312     }
313         
314 }
315
316 //__________________________________________________________________________
317
318 void AliITSsimulationSSD::ApplyDAQ() {
319   // Converts simulated signals to simulated ADC counts
320     AliITS *its=(AliITS*)gAlice->GetModule("ITS");
321
322     Float_t noiseP, noiseN;
323     fResponse->GetNoiseParam(noiseP,noiseN);
324
325     char opt[30],dummy[20];
326     fResponse->ParamOptions(opt,dummy);
327
328     Int_t i,j;
329     if (strstr(opt,"SetInvalid")) {
330       // Set signal = 0 if invalid strip
331       for(i=0; i<fNstrips; i++) {
332          if (!(fDCS->IsValidP(i))) (*fP)[i] = 0;
333          if (!(fDCS->IsValidN(i))) (*fN)[i] = 0;
334       }
335     }
336     
337     Int_t digits[3], tracks[3], hits[3];
338     Float_t charges[3];
339     Float_t phys=0;
340     for(i=0;i<3;i++) tracks[i]=0;
341     for(i=0; i<fNstrips; i++) { 
342        if( (strstr(opt,"SetInvalid") && (*fP)[i] < noiseP*4) || !(*fP)[i]) continue;
343           digits[0]=1;
344           digits[1]=i;
345           digits[2]=(int)(*fP)[i];
346           for(j=0; j<(fTracksP+i)->GetNTracks(); j++) {
347             if(j>2) continue;
348             if((fTracksP+i)->GetNTracks()) tracks[j]=(fTracksP+i)->GetTrack(j);
349             else tracks[j]=-2;
350             charges[j] = 0;
351             hits[j] = -1;
352           }
353           its->AddSimDigit(2,phys,digits,tracks,hits,charges);
354           
355     }
356     
357     
358     for(i=0; i<fNstrips; i++) {
359        if( (strstr(opt,"SetInvalid") && (*fN)[i] < noiseN*4)|| !(*fN)[i]) continue;
360           digits[0]=0;
361           digits[1]=i;
362           digits[2]=(int)(*fN)[i];
363           for( j=0; j<(fTracksN+i)->GetNTracks(); j++) {
364             if(j>2) continue;
365             if((fTracksN+i)->GetNTracks()) tracks[j]=(fTracksN+i)->GetTrack(j);
366             else tracks[j]=-2;
367             charges[j] = 0;
368             hits[j] = -1;
369           }
370           its->AddSimDigit(2,phys,digits,tracks,hits,charges);
371           
372     }
373     
374 }
375
376
377 //____________________________________________________________________________
378
379 Float_t AliITSsimulationSSD::F(Float_t x, Float_t s) {
380   // Computes the integral of a gaussian at the mean valuse x with sigma s.
381
382     //printf("SDD:F(%e,%e)\n",x,s);
383     return 0.5*TMath::Erf(x * fPitch / s) ;
384
385
386 //______________________________________________________________________
387
388 Float_t AliITSsimulationSSD::Get2Strip(Int_t flag, Int_t iStrip, Float_t x, Float_t z){
389   // Returns the relative space between two strips.
390
391     // flag==1 for Pside, 0 for Nside
392
393     Float_t stereoP, stereoN;
394     fSegmentation->Angles(stereoP,stereoN);
395     
396     Float_t tanP=TMath::Tan(stereoP);
397     Float_t tanN=TMath::Tan(stereoN);
398  
399     Float_t dx = fSegmentation->Dx();
400     Float_t dz = fSegmentation->Dz();
401
402
403      x += dx/2;
404      z += dz/2; 
405     
406      if (flag) return (x - z*tanP) / fPitch - iStrip;       // from 0 to 1
407      else  return (x - tanN*(dz - z)) / fPitch - iStrip;
408 }
409 //____________________________________________________________________________