new data member: fDebug=kFALSE by default. Setter and Getter defined in the base...
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSPD.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 /* $Id$ */
17
18 #include <Riostream.h>
19 #include <TRandom.h>
20 #include <TH1.h>
21 #include <TMath.h>
22 #include <TString.h>
23 #include <TParticle.h>
24
25 #include "AliRun.h"
26 #include "AliITS.h"
27 #include "AliITShit.h"
28 #include "AliITSdigit.h"
29 #include "AliITSmodule.h"
30 #include "AliITSMapA2.h"
31 #include "AliITSpList.h"
32 #include "AliITSsimulationSPD.h"
33 #include "AliITSsegmentation.h"
34 #include "AliITSresponse.h"
35 #include "AliITSsegmentationSPD.h"
36 #include "AliITSresponseSPD.h"
37
38 //#define DEBUG
39
40 ClassImp(AliITSsimulationSPD)
41 ////////////////////////////////////////////////////////////////////////
42 // Version: 0
43 // Written by Rocco Caliandro
44 // from a model developed with T. Virgili and R.A. Fini
45 // June 15 2000
46 //
47 // AliITSsimulationSPD is the simulation of SPDs
48 //
49 //______________________________________________________________________
50 AliITSsimulationSPD::AliITSsimulationSPD(){
51     // Default constructor
52
53     fResponse     = 0;
54     fSegmentation = 0;
55     fHis          = 0;
56     fMapA2        = 0;
57
58 /*
59     fThresh       = 0.;
60     fSigma        = 0.;
61     fCouplCol     = 0.;
62     fCouplRow     = 0.; */
63 }
64 //______________________________________________________________________
65 AliITSsimulationSPD::AliITSsimulationSPD(AliITSsegmentation *seg,
66                                          AliITSresponse *resp) {
67     // Standard constructor
68
69     fResponse     = 0;
70     fSegmentation = 0;
71     fHis          = 0;
72     fMapA2        = 0;
73     SetDebug(kFALSE);
74
75 /*
76     fThresh       = 0.;
77     fSigma        = 0.;
78     fCouplCol     = 0.;
79     fCouplRow     = 0.*/
80     Init((AliITSsegmentationSPD*)seg,(AliITSresponseSPD*)resp);
81 }
82 //______________________________________________________________________
83 void AliITSsimulationSPD::Init(AliITSsegmentationSPD *seg,
84                                AliITSresponseSPD *resp) {
85     // Initilizes the variables of AliITSsimulation SPD.
86
87     fHis = 0;
88     fResponse     = resp;
89     fSegmentation = seg;
90     fMapA2  = new AliITSMapA2(fSegmentation);
91     fpList  = new AliITSpList(GetNPixelsZ()+1,GetNPixelsX()+1);
92 /*
93     fResponse->Thresholds(fThresh,fSigma);
94     fResponse->GetNoiseParam(fCouplCol,fCouplRow);
95     fNPixelsZ = fSegmentation->Npz();
96     fNPixelsX = fSegmentation->Npx();
97 */
98 }
99 //______________________________________________________________________
100 AliITSsimulationSPD::~AliITSsimulationSPD() { 
101     // destructor
102
103     delete fMapA2;
104 //    delete fpList;
105
106     if (fHis) {
107         fHis->Delete(); 
108         delete fHis;     
109     } // end if
110 }
111 //______________________________________________________________________
112 AliITSsimulationSPD::AliITSsimulationSPD(const AliITSsimulationSPD &source) :
113     AliITSsimulation(source){
114     // Copy Constructor
115
116     if(&source == this) return;
117
118     this->fMapA2    = source.fMapA2;
119     this->fHis      = source.fHis;
120 /*
121     this->fThresh   = source.fThresh;
122     this->fSigma    = source.fSigma;
123     this->fCouplCol = source.fCouplCol;
124     this->fCouplRow = source.fCouplRow;
125     this->fNPixelsX = source.fNPixelsX;
126     this->fNPixelsZ = source.fNPixelsZ;
127 */
128     return;
129 }
130 //______________________________________________________________________
131 AliITSsimulationSPD& AliITSsimulationSPD::operator=(const AliITSsimulationSPD 
132                                                     &source) {
133     //    Assignment operator
134
135     if(&source == this) return *this;
136
137     this->fMapA2    = source.fMapA2;
138     this->fHis      = source.fHis;
139 /*
140     this->fThresh   = source.fThresh;
141     this->fSigma    = source.fSigma;
142     this->fCouplCol = source.fCouplCol;
143     this->fCouplRow = source.fCouplRow;
144     this->fNPixelsX = source.fNPixelsX;
145     this->fNPixelsZ = source.fNPixelsZ;
146 */
147     return *this;
148
149 //______________________________________________________________________
150 void AliITSsimulationSPD::InitSimulationModule(Int_t module,Int_t event){
151     // Creates maps to build the list of tracks for each sumable digit
152     // Inputs:
153     //   Int_t module    // Module number to be simulated
154     //   Int_t event     // Event number to be simulated
155     // Outputs:
156     //   none.
157     // Return
158     //    none.
159  
160     fModule = module;
161     fEvent  = event;
162     fMapA2->ClearMap();
163     fpList->ClearMap();
164 }
165 //______________________________________________________________________
166 void AliITSsimulationSPD::FinishSDigitiseModule(){
167     // Does the Sdigits to Digits work
168     // Inputs:
169     //   none.
170     // Outputs:
171     //   none.
172     // Return:
173     //   none.
174
175     SDigitsToDigits(fModule,fpList);
176 }
177 //______________________________________________________________________
178 void AliITSsimulationSPD::SDigitiseModule(AliITSmodule *mod, Int_t dummy0,
179                                              Int_t dummy1) {
180     // Sum digitize module
181     if (!(mod->GetNhits())) return; // if module has no hits then no Sdigits.
182     Int_t    number     = 10000;
183     Int_t    *frowpixel = new Int_t[number];
184     Int_t    *fcolpixel = new Int_t[number];
185     Double_t *fenepixel = new Double_t[number];
186
187     dummy0 = dummy1; // remove unsued variable warning.
188     fModule = mod->GetIndex();
189
190     // Array of pointers to store the track index of the digits
191     // leave +1, otherwise pList crashes when col=256, row=192
192
193     HitsToAnalogDigits(mod,frowpixel,fcolpixel,fenepixel,fpList);
194
195     WriteSDigits(fpList);
196
197     // clean memory
198     delete[] frowpixel;
199     delete[] fcolpixel;
200     delete[] fenepixel;
201     fMapA2->ClearMap();
202     fpList->ClearMap();
203 }
204 //______________________________________________________________________
205 void AliITSsimulationSPD::DigitiseModule(AliITSmodule *mod, Int_t dummy0,
206                                              Int_t dummy1) {
207     // digitize module. Also need to digitize modules with only noise.
208
209     Int_t    number     = 10000;
210     Int_t    *frowpixel = new Int_t[number];
211     Int_t    *fcolpixel = new Int_t[number];
212     Double_t *fenepixel = new Double_t[number];
213
214     dummy0 = dummy1; // remove unsued variable warning.
215     // Array of pointers to store the track index of the digits
216     // leave +1, otherwise pList crashes when col=256, row=192
217     fModule = mod->GetIndex();
218     // noise setting
219     SetFluctuations(fpList,fModule);
220
221     HitsToAnalogDigits(mod,frowpixel,fcolpixel,fenepixel,fpList);
222
223     // apply mask to SPD module
224     SetMask();
225
226     CreateDigit(fModule,fpList);
227
228     // clean memory
229     delete[] frowpixel;
230     delete[] fcolpixel;
231     delete[] fenepixel;
232     fMapA2->ClearMap();
233     fpList->ClearMap();
234 }
235 //______________________________________________________________________
236 void AliITSsimulationSPD::SDigitsToDigits(Int_t module,AliITSpList *pList) {
237     // sum digits to Digits.
238
239 #ifdef DEBUG
240     cout << "Entering AliITSsimulatinSPD::SDigitsToDigits for module=";
241     cout << module << endl;
242 #endif
243     fModule = module;
244
245     // noise setting
246     SetFluctuations(pList,module);
247
248     fMapA2->ClearMap(); // since noise is in pList aready. Zero Map so that
249     // noise is not doubled when calling FillMapFrompList.
250
251     FillMapFrompList(pList);
252
253     // apply mask to SPD module
254     SetMask();
255
256     CreateDigit(module,pList);
257
258     fMapA2->ClearMap();
259     pList->ClearMap();
260 }
261 //______________________________________________________________________
262 void AliITSsimulationSPD::UpdateMapSignal(Int_t row,Int_t col,Int_t trk,
263                                           Int_t hit,Int_t mod,Double_t ene,
264                                           AliITSpList *pList) {
265     // updates the Map of signal, adding the energy  (ene) released by
266     // the current track
267
268     fMapA2->AddSignal(row,col,ene);
269     pList->AddSignal(row,col,trk,hit,mod,ene);
270 }
271 //______________________________________________________________________
272 void AliITSsimulationSPD::UpdateMapNoise(Int_t row,Int_t col,Int_t mod,
273                                          Double_t ene,AliITSpList *pList) {
274     // updates the Map of noise, adding the energy  (ene) give my noise
275
276     fMapA2->AddSignal(row,col,ene);
277     pList->AddNoise(row,col,mod,ene);
278 }
279 //______________________________________________________________________
280 void AliITSsimulationSPD::HitsToAnalogDigits(AliITSmodule *mod,
281                                              Int_t *frowpixel,Int_t *fcolpixel,
282                                              Double_t *fenepixel,
283                                              AliITSpList *pList) {
284     // Loops over all hits to produce Analog/floting point digits. This
285     // is also the first task in producing standard digits.
286     
287     // loop over hits in the module
288     Int_t hitpos,nhits = mod->GetNhits();
289     for (hitpos=0;hitpos<nhits;hitpos++) {
290         HitToDigit(mod,hitpos,frowpixel,fcolpixel,fenepixel,pList);
291     }// end loop over digits
292 }
293 //______________________________________________________________________
294 void AliITSsimulationSPD::HitToDigit(AliITSmodule *mod,Int_t hitpos,
295                                      Int_t *frowpixel,Int_t *fcolpixel,
296                                      Double_t *fenepixel,AliITSpList *pList) {
297     //  Steering function to determine the digits associated to a given
298     // hit (hitpos)
299     // The digits are created by charge sharing (ChargeSharing) and by
300     // capacitive coupling (SetCoupling). At all the created digits is
301     // associated the track number of the hit (ntrack)
302     Double_t x1l=0.0,y1l=0.0,z1l=0.0,x2l=0.0,y2l=0.0,z2l=0.0;
303     Int_t r1,r2,c1,c2,row,col,npixel = 0;
304     Int_t ntrack;
305     Double_t ene=0.0,etot=0.0;
306     const Float_t kconv = 10000.;     // cm -> microns
307     const Float_t kconv1= 0.277e9;    // GeV -> electrons equivalent
308
309     if(!(mod->LineSegmentL(hitpos,x1l,x2l,y1l,y2l,z1l,z2l,etot,ntrack)))return;
310
311     x2l += x1l; y2l += y1l; z2l += z1l; // Convert to ending coordinate.
312     // positions shifted and converted in microns
313     x1l   = x1l*kconv + fSegmentation->Dx()/2.;
314     z1l   = z1l*kconv + fSegmentation->Dz()/2.;
315     // positions  shifted and converted in microns
316     x2l   = x2l*kconv + fSegmentation->Dx()/2.;
317     z2l   = z2l*kconv + fSegmentation->Dz()/2.;
318     etot *= kconv1; // convert from GeV to electrons equivalent.
319     Int_t module = mod->GetIndex();
320
321     // to account for the effective sensitive area
322     // introduced in geometry 
323     if (z1l<0 || z1l>fSegmentation->Dz()) return;
324     if (z2l<0 || z2l>fSegmentation->Dz()) return;
325     if (x1l<0 || x1l>fSegmentation->Dx()) return;
326     if (x2l<0 || x2l>fSegmentation->Dx()) return;
327
328     //Get the col and row number starting from 1
329     // the x direction is not inverted for the second layer!!!
330     fSegmentation->GetPadIxz(x1l, z1l, c1, r1); 
331     fSegmentation->GetPadIxz(x2l, z2l, c2, r2);
332
333     // to account for unexpected equal entrance and 
334     // exit coordinates
335     if (x1l==x2l) x2l=x2l+x2l*0.1;
336     if (z1l==z2l) z2l=z2l+z2l*0.1;
337
338     if ((r1==r2) && (c1==c2)){
339         // no charge sharing
340         npixel = 1;              
341         frowpixel[npixel-1] = r1;
342         fcolpixel[npixel-1] = c1;
343         fenepixel[npixel-1] = etot;
344     } else {
345         // charge sharing
346         ChargeSharing(x1l,z1l,x2l,z2l,c1,r1,c2,r2,etot,
347                       npixel,frowpixel,fcolpixel,fenepixel);
348     } // end if r1==r2 && c1==c2.
349
350     for (Int_t npix=0;npix<npixel;npix++){
351         row = frowpixel[npix];
352         col = fcolpixel[npix];
353         ene = fenepixel[npix];
354         UpdateMapSignal(row,col,ntrack,hitpos,module,ene,pList); 
355         // Starting capacitive coupling effect
356         SetCoupling(row,col,ntrack,hitpos,module,pList); 
357     } // end for npix
358 }
359 //______________________________________________________________________
360 void AliITSsimulationSPD::ChargeSharing(Float_t x1l,Float_t z1l,Float_t x2l,
361                                         Float_t z2l,Int_t c1,Int_t r1,Int_t c2,
362                                         Int_t r2,Float_t etot,
363                                         Int_t &npixel,Int_t *frowpixel,
364                                         Int_t *fcolpixel,Double_t *fenepixel){
365     //  Take into account the geometrical charge sharing when the track
366     //  crosses more than one pixel.
367     //
368     //Begin_Html
369     /*
370       <img src="picts/ITS/barimodel_2.gif">
371       </pre>
372       <br clear=left>
373       <font size=+2 color=red>
374       <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
375       </font>
376       <pre>
377     */
378     //End_Html
379     //Float_t dm;
380     Float_t xa,za,xb,zb,dx,dz,dtot,refr,refm,refc;
381     Float_t refn=0.;
382     Float_t arefm, arefr, arefn, arefc, azb, az2l, axb, ax2l;
383     Int_t   dirx,dirz,rb,cb;
384     Int_t flag,flagrow,flagcol;
385     Double_t epar;
386
387     npixel = 0;
388     xa     = x1l;
389     za     = z1l;
390 //    dx     = x1l-x2l;
391 //    dz     = z1l-z2l;
392     dx     = x2l-x1l;
393     dz     = z2l-z1l;
394     dtot   = TMath::Sqrt((dx*dx)+(dz*dz));   
395     if (dtot==0.0) dtot = 0.01;
396     dirx   = (Int_t) TMath::Sign((Float_t)1,dx);
397     dirz   = (Int_t) TMath::Sign((Float_t)1,dz);
398
399     // calculate the x coordinate of  the pixel in the next column    
400     // and the z coordinate of  the pixel in the next row
401     Float_t xpos, zpos;
402
403     fSegmentation->GetPadCxz(c1, r1-1, xpos, zpos); 
404
405     Float_t xsize = fSegmentation->Dpx(0);
406     Float_t zsize = fSegmentation->Dpz(r1-1);
407     
408     if (dirx == 1) refr = xpos+xsize/2.;
409     else refr = xpos-xsize/2.;
410
411     if (dirz == 1) refn = zpos+zsize/2.;
412     else refn = zpos-zsize/2.;
413
414     flag = 0;
415     flagrow = 0;
416     flagcol = 0;
417     do{
418         // calculate the x coordinate of the intersection with the pixel
419         // in the next cell in row  direction
420       if(dz!=0)
421         refm = dx*((refn - z1l)/dz) + x1l;
422       else
423         refm = refr+dirx*xsize;
424    
425         // calculate the z coordinate of the intersection with the pixel
426         // in the next cell in column direction
427       if (dx!=0)
428         refc = dz*((refr - x1l)/dx) + z1l;
429       else
430         refc = refn+dirz*zsize;
431
432         arefm = refm * dirx;
433         arefr = refr * dirx;
434         arefn = refn * dirz;
435         arefc = refc * dirz;
436
437         if ((arefm < arefr) && (arefn < arefc)){
438             // the track goes in the pixel in the next cell in row direction
439             xb = refm;
440             zb = refn;
441             cb = c1;
442             rb = r1 + dirz;
443             azb = zb * dirz;
444             az2l = z2l * dirz;
445             if (rb == r2) flagrow=1;
446             if (azb > az2l) {
447                 zb = z2l;
448                 xb = x2l;
449             } // end if
450             // shift to the pixel in the next cell in row direction
451             Float_t zsizeNext = fSegmentation->Dpz(rb-1);
452             //to account for cell at the borders of the detector
453             if(zsizeNext==0) zsizeNext = zsize;
454             refn += zsizeNext*dirz;
455         }else {
456             // the track goes in the pixel in the next cell in column direction
457             xb = refr;
458             zb = refc;
459             cb = c1 + dirx;
460             rb = r1;
461             axb = xb * dirx;
462             ax2l = x2l * dirx;
463             if (cb == c2) flagcol=1;
464             if (axb > ax2l) {
465                 zb = z2l;
466                 xb = x2l;
467             } // end ifaxb > ax2l
468
469             // shift to the pixel in the next cell in column direction
470             Float_t xsizeNext = fSegmentation->Dpx(cb-1);
471             //to account for cell at the borders of the detector
472             if(xsizeNext==0) xsizeNext = xsize;
473             refr += xsizeNext*dirx;
474         } // end if (arefm < arefr) && (arefn < arefc)
475
476         //calculate the energy lost in the crossed pixel      
477         epar = TMath::Sqrt((xb-xa)*(xb-xa)+(zb-za)*(zb-za)); 
478         epar = etot*(epar/dtot);
479
480         //store row, column and energy lost in the crossed pixel
481         frowpixel[npixel] = r1;
482         fcolpixel[npixel] = c1;
483         fenepixel[npixel] = epar;
484         npixel++;
485  
486         // the exit point of the track is reached
487         if (epar == 0) flag = 1;
488         if ((r1 == r2) && (c1 == c2)) flag = 1;
489         if (flag!=1) {
490             r1 = rb;
491             c1 = cb;
492             xa = xb;
493             za = zb;
494         } // end if flag!=1
495     } while (flag==0);
496 }
497 //______________________________________________________________________
498 void AliITSsimulationSPD::SetCoupling(Int_t row, Int_t col, Int_t ntrack,
499                                       Int_t idhit,Int_t module,
500                                       AliITSpList *pList) {
501     //  Take into account the coupling between adiacent pixels.
502     //  The parameters probcol and probrow are the probability of the
503     //  signal in one pixel shared in the two adjacent pixels along
504     //  the column and row direction, respectively.
505     //
506     //Begin_Html
507     /*
508       <img src="picts/ITS/barimodel_3.gif">
509       </pre>
510       <br clear=left>
511       <font size=+2 color=red>
512       <a href="mailto:tiziano.virgili@cern.ch"></a>.
513       </font>
514       <pre>
515     */
516     //End_Html
517     Int_t j1,j2,flag=0;
518     Double_t pulse1,pulse2;
519     Float_t couplR=0.0,couplC=0.0;
520     Double_t xr=0.;
521
522     GetCouplings(couplR,couplC);
523     j1 = row;
524     j2 = col;
525     pulse1 = fMapA2->GetSignal(row,col);
526     pulse2 = pulse1;
527     for (Int_t isign=-1;isign<=1;isign+=2){// loop in row direction
528         do{
529             j1 += isign;
530             //   pulse1 *= couplR; 
531           xr = gRandom->Rndm();
532           //   if ((j1<0) || (j1>GetNPixelsZ()-1) || (pulse1<GetThreshold())){
533             if ((j1<0) || (j1>GetNPixelsZ()-1) || (xr>couplR)){
534                 j1 = row;
535                 flag = 1;
536             }else{
537                 UpdateMapSignal(j1,col,ntrack,idhit,module,pulse1,pList);
538               //  flag = 0;
539               flag = 1; // only first next!!
540             } // end if
541         } while(flag == 0);
542         // loop in column direction
543       do{
544           j2 += isign;
545           // pulse2 *= couplC; 
546           xr = gRandom->Rndm();
547           //  if ((j2<0) || (j2>(GetNPixelsX()-1)) || (pulse2<GetThreshold())){
548             if ((j2<0) || (j2>GetNPixelsX()-1) || (xr>couplC)){
549               j2 = col;
550               flag = 1;
551           }else{
552               UpdateMapSignal(row,j2,ntrack,idhit,module,pulse2,pList);
553               //  flag = 0;
554               flag = 1; // only first next!!
555           } // end if
556       } while(flag == 0);
557     } // for isign
558 }
559 //______________________________________________________________________
560 void AliITSsimulationSPD::SetCouplingOld(Int_t row, Int_t col, Int_t ntrack,
561                                          Int_t idhit,Int_t module,
562                                          AliITSpList *pList) {
563     //  Take into account the coupling between adiacent pixels.
564     //  The parameters probcol and probrow are the fractions of the
565     //  signal in one pixel shared in the two adjacent pixels along
566     //  the column and row direction, respectively.
567     //
568     //Begin_Html
569     /*
570       <img src="picts/ITS/barimodel_3.gif">
571       </pre>
572       <br clear=left>
573       <font size=+2 color=red>
574       <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
575       </font>
576       <pre>
577     */
578     //End_Html
579     Int_t j1,j2,flag=0;
580     Double_t pulse1,pulse2;
581     Float_t couplR=0.0,couplC=0.0;
582
583     GetCouplings(couplR,couplC);
584     j1 = row;
585     j2 = col;
586     pulse1 = fMapA2->GetSignal(row,col);
587     pulse2 = pulse1;
588     for (Int_t isign=-1;isign<=1;isign+=2){// loop in row direction
589         do{
590             j1 += isign;
591             pulse1 *= couplR;
592             if ((j1<0) || (j1>GetNPixelsZ()-1) || (pulse1<GetThreshold())){
593                 pulse1 = fMapA2->GetSignal(row,col);
594                 j1 = row;
595                 flag = 1;
596             }else{
597                 UpdateMapSignal(j1,col,ntrack,idhit,module,pulse1,pList);
598                 flag = 0;
599             } // end if
600         } while(flag == 0);
601         // loop in column direction
602       do{
603           j2 += isign;
604           pulse2 *= couplC;
605           if ((j2<0) || (j2>(GetNPixelsX()-1)) || (pulse2<GetThreshold())){
606               pulse2 = fMapA2->GetSignal(row,col);
607               j2 = col;
608               flag = 1;
609           }else{
610               UpdateMapSignal(row,j2,ntrack,idhit,module,pulse2,pList);
611               flag = 0;
612           } // end if
613       } while(flag == 0);
614     } // for isign
615 }
616 //______________________________________________________________________
617 void AliITSsimulationSPD::CreateDigit(Int_t module,AliITSpList *pList) {
618     // The pixels are fired if the energy deposited inside them is above
619     // the threshold parameter ethr. Fired pixed are interpreted as digits
620     // and stored in the file digitfilename. One also needs to write out
621     // cases when there is only noise (nhits==0).
622
623     static AliITS *aliITS  = (AliITS*)gAlice->GetModule("ITS");
624
625     Int_t size = AliITSdigitSPD::GetNTracks();
626     Int_t * digits = new Int_t[size];
627     Int_t * tracks = new Int_t[size];
628     Int_t * hits = new Int_t[size];
629     Float_t * charges = new Float_t[size]; 
630     Int_t j1;
631
632     module=0; // remove unused variable warning.
633     for(j1=0;j1<size;j1++){tracks[j1]=-3;hits[j1]=-1;charges[j1]=0.0;}
634     for (Int_t r=1;r<=GetNPixelsZ();r++) {
635         for (Int_t c=1;c<=GetNPixelsX();c++) {
636             // check if the deposited energy in a pixel is above the
637             // threshold 
638             Float_t signal = (Float_t) fMapA2->GetSignal(r,c);
639             if ( signal > GetThreshold()) {
640                 digits[0] = r-1;  // digits starts from 0
641                 digits[1] = c-1;  // digits starts from 0
642                 //digits[2] = 1;  
643                 digits[2] =  (Int_t) signal;  // the signal is stored in
644                                               //  electrons
645                 for(j1=0;j1<size;j1++){
646                     if(j1<pList->GetNEnteries()){
647                         tracks[j1] = pList->GetTrack(r,c,j1);
648                         hits[j1]   = pList->GetHit(r,c,j1);
649                         //}else{
650                         //tracks[j1] = -3;
651                         //hits[j1]   = -1;
652                     } // end if
653                     //charges[j1] = 0;
654                 } // end for j1
655                 Float_t phys = 0;
656                 aliITS->AddSimDigit(0,phys,digits,tracks,hits,charges);
657 #ifdef DEBUG
658                 cout << " CreateSPDDigit mod=" << fModule << " r,c=" << r;
659                 cout <<","<<c<< " sig=" << fpList->GetSignalOnly(r,c);
660                 cout << " noise=" << fpList->GetNoise(r,c);
661                 cout << " Msig="<< signal << " Thres=" << GetThreshold()<<endl;
662 #endif
663             } // end if of threshold condition
664         } // for c
665     }// end do on pixels
666     delete [] digits;
667     delete [] tracks;
668     delete [] hits;
669     delete [] charges;
670
671 }
672 //______________________________________________________________________
673 void AliITSsimulationSPD::SetFluctuations(AliITSpList *pList,Int_t module) {
674     //  Set the electronic noise and threshold non-uniformities to all the
675     //  pixels in a detector.
676     //  The parameter fSigma is the squared sum of the sigma due to noise
677     //  and the sigma of the threshold distribution among pixels.
678     //
679     //Begin_Html
680     /*
681       <img src="picts/ITS/barimodel_1.gif">
682       </pre>
683       <br clear=left>
684       <font size=+2 color=red>
685       <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
686       </font>
687       <pre>
688     */
689     //End_Html
690     Float_t  thr=0.0,sigm=0.0;
691     Double_t signal,sigma;
692     Int_t iz,ix;
693
694     GetThresholds(thr,sigm);
695     sigma = (Double_t) sigm;
696     for(iz=1;iz<=GetNPixelsZ();iz++){
697         for(ix=1;ix<=GetNPixelsX();ix++){
698             signal = sigma*gRandom->Gaus(); 
699             fMapA2->SetHit(iz,ix,signal);
700             // insert in the label-signal-hit list the pixels fired
701             // only by noise
702             pList->AddNoise(iz,ix,module,signal);
703         } // end of loop on pixels
704     } // end of loop on pixels
705 }
706 //______________________________________________________________________
707 void AliITSsimulationSPD::SetMask() {
708     //  Apply a mask to the SPD module. 1% of the pixel channels are
709     //  masked. When the database will be ready, the masked pixels
710     //  should be read from it.
711     Double_t signal;
712     Int_t iz,ix,im;
713     Float_t totMask;
714     Float_t perc = ((AliITSresponseSPD*)fResponse)->GetFractionDead();
715     // in this way we get the same set of random numbers for all runs.
716     // This is a cluge for now.
717     static TRandom *rnd = new TRandom();
718
719     totMask= perc*GetNPixelsZ()*GetNPixelsX();
720     for(im=1;im<totMask;im++){
721         do{
722             ix=(Int_t)(rnd->Rndm()*(GetNPixelsX()-1.) + 1.);
723         } while(ix<=0 || ix>GetNPixelsX());
724         do{
725             iz=(Int_t)(rnd->Rndm()*(GetNPixelsZ()-1.) + 1.);
726         } while(iz<=0 || iz>GetNPixelsZ());
727         signal = -1.;
728         fMapA2->SetHit(iz,ix,signal);
729     } // end loop on masked pixels
730 }
731 //______________________________________________________________________
732 void AliITSsimulationSPD::CreateHistograms() {
733     // Create Histograms
734     Int_t i;
735
736     fHis=new TObjArray(GetNPixelsZ());
737     for(i=0;i<GetNPixelsZ();i++) {
738         TString spdname("spd_");
739         Char_t candnum[4];
740         sprintf(candnum,"%d",i+1);
741         spdname.Append(candnum);
742         (*fHis)[i] = new TH1F(spdname.Data(),"SPD maps",
743                               GetNPixelsX(),0.,(Float_t) GetNPixelsX());
744     } // end for i
745 }
746 //______________________________________________________________________
747 void AliITSsimulationSPD::ResetHistograms() {
748     // Reset histograms for this detector
749     Int_t i;
750
751     for(i=0;i<GetNPixelsZ();i++ ) {
752         if ((*fHis)[i])    ((TH1F*)(*fHis)[i])->Reset();
753     } // end for i
754 }
755 //______________________________________________________________________
756 void AliITSsimulationSPD::WriteSDigits(AliITSpList *pList){
757     // Fills the Summable digits Tree
758     Int_t i,ni,j,nj;
759     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
760
761     pList->GetMaxMapIndex(ni,nj);
762     for(i=0;i<ni;i++)for(j=0;j<nj;j++){
763         if(pList->GetSignalOnly(i,j)>0.0){
764             aliITS->AddSumDigit(*(pList->GetpListItem(i,j)));
765 #ifdef DEBUG
766             cout << "pListSPD: " << *(pList->GetpListItem(i,j)) << endl;
767             cout << " CreateSPDSDigit mod=" << fModule << " r,c=";
768             cout << i  <<","<< j << " sig=" << fpList->GetSignalOnly(i,j);
769             cout << " noise=" << fpList->GetNoise(i,j) <<endl;
770 #endif
771         } // end if
772     } // end for i,j
773     return;
774 }
775 //______________________________________________________________________
776 void AliITSsimulationSPD::FillMapFrompList(AliITSpList *pList){
777     // Fills fMap2A from the pList of Summable digits
778     Int_t ix,iz;
779
780     for(iz=0;iz<GetNPixelsZ();iz++)for(ix=0;ix<GetNPixelsX();ix++) 
781         fMapA2->AddSignal(iz,ix,pList->GetSignal(iz,ix));
782     return;
783 }