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