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