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