ITS/AliITSsimulationSPDdubna.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 $Id$
  18 */
  19 #include <Riostream.h>
  20 #include <TRandom.h>
  21 #include <TH1.h>
  22 #include <TMath.h>
  23 #include <TString.h>
  24 #include <TParticle.h>
  25
  26 #include "AliRun.h"
  27 #include "AliITS.h"
  28 #include "AliITShit.h"
  29 #include "AliITSdigit.h"
  30 #include "AliITSmodule.h"
  31 #include "AliITSMapA2.h"
  32 #include "AliITSpList.h"
  33 #include "AliITSsimulationSPDdubna.h"
  34 #include "AliITSsegmentationSPD.h"
  35 #include "AliITSresponseSPDdubna.h"
  36
  37 //#define DEBUG
  38
  39 ClassImp(AliITSsimulationSPDdubna)
  40 ////////////////////////////////////////////////////////////////////////
  41 // Version: 1
  42 // Modified by Bjorn S. Nilsen
  43 // Version: 0
  44 // Written by Boris Batyunya
  45 // December 20 1999
  46 //
  47 // AliITSsimulationSPDdubna is to do the simulation of SPDs.
  48 //______________________________________________________________________
  49 AliITSsimulationSPDdubna::AliITSsimulationSPDdubna() : AliITSsimulation(){
  50     // Default constructor.
  51     // Inputs:
  52     //    none.
  53     // Outputs:
  54     //    none.
  55     // Return:
  56     //    A default constructed AliITSsimulationSPDdubna class.
  57
  58     fHis          = 0;
  59 }
  60 //______________________________________________________________________
  61 AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(AliITSsegmentation *seg,
  62                                                    AliITSresponse *resp){
  63     // standard constructor
  64     // Inputs:
  65     //    AliITSsegmentation *seg  A pointer to the segmentation class
  66     //                             to be used for this simulation
  67     //    AliITSresponse     *resp A pointer to the responce class to
  68     //                             be used for this simulation
  69     // Outputs:
  70     //    none.
  71     // Return:
  72     //    A default constructed AliITSsimulationSPDdubna class.
  73
  74     fHis = 0;
  75     Init(seg,resp);
  76 }
  77 //______________________________________________________________________
  78 void AliITSsimulationSPDdubna::Init(AliITSsegmentation *seg,
  79                                                    AliITSresponse *resp){
  80     // Initilization
  81     // Inputs:
  82     //    AliITSsegmentation *seg  A pointer to the segmentation class
  83     //                             to be used for this simulation
  84     //    AliITSresponse     *resp A pointer to the responce class to
  85     //                             be used for this simulation
  86     // Outputs:
  87     //    none.
  88     // Return:
  89     //    none.
  90     const Double_t kmictocm = 1.0e-4; // convert microns to cm.
  91
  92     SetResponseModel(resp);
  93     SetSegmentationModel(seg);
  94     SetModuleNumber(0);
  95     SetEventNumber(0);
  96     SetMap(new AliITSpList(GetNPixelsZ(),GetNPixelsX()));
  97
  98     GetResp()->SetDistanceOverVoltage(kmictocm*GetSeg()->Dy(),50.0);
  99 }
 100 //______________________________________________________________________
 101 AliITSsimulationSPDdubna::~AliITSsimulationSPDdubna(){
 102     // destructor
 103     // Inputs:
 104     //    none.
 105     // Outputs:
 106     //    none.
 107     // Return:
 108     //     none.
 109
 110     if (fHis) {
 111         fHis->Delete();
 112         delete fHis;
 113     } // end if fHis
 114 }
 115 //______________________________________________________________________
 116 AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(const
 117                                                    AliITSsimulationSPDdubna
 118                                                    &s) : AliITSsimulation(s){
 119     //     Copy Constructor
 120     // Inputs:
 121     //    AliITSsimulationSPDdubna &s The original class for which
 122     //                                this class is a copy of
 123     // Outputs:
 124     //    none.
 125     // Return:
 126
 127     *this = s;
 128     return;
 129 }
 130 //______________________________________________________________________
 131 AliITSsimulationSPDdubna&  AliITSsimulationSPDdubna::operator=(const
 132                                            AliITSsimulationSPDdubna &s){
 133     //    Assignment operator
 134     // Inputs:
 135     //    AliITSsimulationSPDdubna &s The original class for which
 136     //                                this class is a copy of
 137     // Outputs:
 138     //    none.
 139     // Return:
 140
 141     if(&s == this) return *this;
 142     this->fHis      = s.fHis;
 143     return *this;
 144 }
 145 //______________________________________________________________________
 146 void AliITSsimulationSPDdubna::InitSimulationModule(Int_t module, Int_t event){
 147     //  This function creates maps to build the list of tracks for each
 148     //  summable digit. Inputs defined by base class.
 149     //  Inputs:
 150     //    Int_t module   // Module number to be simulated
 151     //    Int_t event    // Event number to be simulated
 152     //  Outputs:
 153     //    none
 154     //  Returns:
 155     //    none
 156
 157     SetModuleNumber(module);
 158     SetEventNumber(event);
 159     ClearMap();
 160 }
 161 //_____________________________________________________________________
 162 void AliITSsimulationSPDdubna::SDigitiseModule(AliITSmodule *mod, Int_t mask,
 163                                                Int_t event){
 164     //  This function begins the work of creating S-Digits.  Inputs defined
 165     //  by base class.
 166     //  Inputs:
 167     //    AliITSmodule *mod  //  module
 168     //    Int_t mask         //  mask to be applied to the module
 169     //    Int_t event        // Event number
 170     //  Outputs:
 171     //    none
 172     //  Return:
 173     //    test              //  test returns kTRUE if the module contained hits
 174     //                      //  test returns kFALSE if it did not contain hits
 175
 176     mask = mod->GetNhits();
 177     if(!mask) return;// if module has no hits don't create Sdigits
 178     SetModuleNumber(mod->GetIndex());
 179     SetEventNumber(event);
 180     HitToSDigit(mod);
 181     WriteSDigits();
 182     ClearMap();
 183 }
 184 //______________________________________________________________________
 185 void AliITSsimulationSPDdubna::WriteSDigits(){
 186     //  This function adds each S-Digit to pList
 187     //  Inputs:
 188     //    none.
 189     //  Outputs:
 190     //    none.
 191     //  Return:
 192     //    none
 193     Int_t ix, nix, iz, niz;
 194     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
 195
 196     GetMap()->GetMaxMapIndex(niz, nix);
 197     for(iz=0; iz<niz; iz++)for(ix=0; ix<nix; ix++){
 198         if(GetMap()->GetSignalOnly(iz,ix)>0.0){
 199             aliITS->AddSumDigit(*(GetMap()->GetpListItem(iz,ix)));
 200 #ifdef DEBUG
 201             cout <<"SDigits " << iz << "," << ix << "," <<
 202                 *(GetMap()->GetpListItem(iz,ix)) << endl;
 203 #endif
 204         } // end if GetMap()->GetSignalOnly(iz,ix)>0.0
 205     } // end for iz,ix
 206     return;
 207 }
 208 //______________________________________________________________________
 209 void AliITSsimulationSPDdubna::FinishSDigitiseModule(){
 210     //  This function calls SDigitsToDigits which creates Digits from SDigits
 211     //  Inputs:
 212     //    none
 213     //  Outputs:
 214     //    none
 215     //  Return
 216     //    none
 217
 218     SDigitsToDigits();
 219     return;
 220 }
 221 //______________________________________________________________________
 222 void AliITSsimulationSPDdubna::SDigitsToDigits(){
 223     //  This function adds electronic noise to the S-Digits and then adds them
 224     // to  a new pList
 225     //  Inputs:
 226     //    none.
 227     //  Outputs:
 228     //    none.
 229     //  Return:
 230     //    none
 231
 232     ChargeToSignal(); // Charge To Signal both adds noise and
 233     ClearMap();
 234 }
 235 //______________________________________________________________________
 236 void AliITSsimulationSPDdubna::DigitiseModule(AliITSmodule *mod, Int_t module,
 237                                               Int_t dummy){
 238     //  This function creates Digits straight from the hits and then adds
 239     //  electronic noise to the digits before adding them to pList
 240     //  Each of the input variables is passed along to HitToSDigit
 241     //  Inputs:
 242     //    AliITSmodule *mod     module
 243     //    Int_t        module   module number  Dummy.
 244     //    Int_t        dummy
 245     //  Outputs:
 246     //     none.
 247     //  Return:
 248     //    none.
 249
 250     //This calls the module for HitToSDigit
 251     fModule = dummy = module = mod->GetIndex();
 252     HitToSDigit(mod);
 253     ChargeToSignal();
 254     ClearMap();
 255 }
 256 //______________________________________________________________________
 257 void AliITSsimulationSPDdubna::UpdateMapSignal(Int_t iz,Int_t ix,Int_t trk,
 258                                                Int_t ht,Double_t signal){
 259     //  This function adds a signal to the pList from the pList class
 260     //  Inputs:
 261     //    Int_t       iz     // row number
 262     //    Int_t       ix     // column number
 263     //    Int_t       trk    // track number
 264     //    Int_t       ht     // hit number
 265     //    Double_t    signal // signal strength
 266     //  Outputs:
 267     //    none.
 268     //  Return:
 269     //    none
 270
 271     GetMap()->AddSignal(iz,ix,trk,ht,GetModuleNumber(),signal);
 272 }
 273 //______________________________________________________________________
 274 void AliITSsimulationSPDdubna::UpdateMapNoise(Int_t iz,Int_t ix,Float_t noise){
 275     //  This function adds noise to data in the MapA2 as well as the pList
 276     //  Inputs:
 277     //    Int_t       iz       // row number
 278     //    Int_t       ix       // column number
 279     //    Int_t       mod     // module number
 280     //    Double_t    sig     // signal strength
 281     //    Double_t    noise   // electronic noise generated by ChargeToSignal
 282     //  Outputs:
 283     //    All of the inputs are passed to AliITSMapA2::AddSignal or
 284     //    AliITSpList::AddNoise
 285     //  Return:
 286     //    none
 287
 288     GetMap()->AddNoise(iz,ix,GetModuleNumber(),noise);
 289 }
 290 //______________________________________________________________________
 291 void AliITSsimulationSPDdubna::HitToDigit(AliITSmodule *mod){
 292     // Standard interface to DigitiseModule
 293     // Inputs:
 294     //    AliITSmodule *mod  Pointer to this module
 295     // Outputs:
 296     //    none.
 297     // Return:
 298     //    none.
 299
 300     DigitiseModule(mod,GetModuleNumber(),0);
 301 }
 302 //______________________________________________________________________
 303 void AliITSsimulationSPDdubna::HitToSDigit(AliITSmodule *mod){
 304     // Does the charge distributions using Gaussian diffusion charge charing.
 305     // Inputs:
 306     //    AliITSmodule *mod  Pointer to this module
 307     // Output:
 308     //    none.
 309     // Return:
 310     //    none.
 311     const Double_t kmictocm = 1.0e-4; // convert microns to cm.
 312     TObjArray *hits = mod->GetHits();
 313     Int_t nhits = hits->GetEntriesFast();
 314     Int_t h,ix,iz;
 315     Int_t idtrack;
 316     Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0;
 317     Double_t x,y,z,t,tp,st,dt=0.2,el,sig;
 318     Double_t thick = kmictocm*GetSeg()->Dy();
 319
 320     if(nhits<=0) return;
 321     for(h=0;h<nhits;h++){
 322 #ifdef DEBUG
 323         cout << "Hits=" << h << "," << *(mod->GetHit(h)) << endl;
 324 #endif
 325         if(mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)){
 326           st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
 327           if(st>0.0){
 328             st = (Double_t)((Int_t)(1.0E+04*st)); // number of microns
 329             if(st<=0.0) st = 1.0;
 330             dt = 1.0/st;
 331             for(t=0;t<1.0;t+=dt){ // Integrate over t
 332               tp  = t+0.5*dt;
 333               el  = GetResp()->GeVToCharge((Float_t)(dt*de));
 334 #ifdef DEBUG
 335               if(el<=0.0) cout << "el="<<el<<" dt="<<dt<<" de="<<de<<endl;
 336 #endif
 337               x   = x0+x1*tp;
 338               y   = y0+y1*tp;
 339               z   = z0+z1*tp;
 340               GetSeg()->LocalToDet(x,z,ix,iz);
 341               sig = GetResp()->SigmaDiffusion1D(thick + y);
 342               SpreadCharge(x,z,ix,iz,el,sig,idtrack,h);
 343             } // end for t
 344           } else { // st == 0.0 deposit it at this point
 345             el  = GetResp()->GeVToCharge((Float_t)de);
 346             x   = x0;
 347             y   = y0;
 348             z   = z0;
 349             GetSeg()->LocalToDet(x,z,ix,iz);
 350             sig = GetResp()->SigmaDiffusion1D(thick + y);
 351             SpreadCharge(x,z,ix,iz,el,sig,idtrack,h);
 352           } // end if st>0.0
 353         } // end if mod->LineSegmentL...
 354     } // Loop over all hits h
 355 }
 356 //______________________________________________________________________
 357 void AliITSsimulationSPDdubna::SpreadCharge(Double_t x0,Double_t z0,
 358                                             Int_t ix0,Int_t iz0,
 359                                             Double_t el,Double_t sig,Int_t t,
 360                                             Int_t hi){
 361     // Spreads the charge over neighboring cells. Assume charge is distributed
 362     // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg)
 363     // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig)
 364     // Defined this way, the integral over all x and z is el.
 365     // Inputs:
 366     //    Double_t x0   x position of point where charge is liberated
 367     //    Double_t y0   y position of point where charge is liberated
 368     //    Double_t z0   z position of point where charge is liberated
 369     //    Int_t    ix0  row of cell corresponding to point x0
 370     //    Int_t    iz0  columb of cell corresponding to point z0
 371     //    Double_t el   number of electrons liberated in this step
 372     //    Double_t sig  Sigma difusion for this step (y0 dependent)
 373     //    Int_t    t    track number
 374     //    Int_t    ti   hit track index number
 375     //    Int_t    hi   hit "hit" index number
 376     // Outputs:
 377     //     none.
 378     // Return:
 379     //     none.
 380     const Int_t knx = 3,knz = 2;
 381     const Double_t kRoot2 = 1.414213562; // Sqrt(2).
 382     const Double_t kmictocm = 1.0e-4; // convert microns to cm.
 383     Int_t ix,iz,ixs,ixe,izs,ize;
 384     Float_t x,z;
 385     Double_t x1,x2,z1,z2,s,sp;
 386
 387     if(sig<=0.0) {
 388         GetMap()->AddSignal(iz0,ix0,t,hi,GetModuleNumber(),el);
 389         return;
 390     } // end if
 391     sp = 1.0/(sig*kRoot2);
 392 #ifdef DEBUG
 393     cout << "sig=" << sig << " sp=" << sp << endl;
 394 #endif
 395     ixs = TMath::Max(-knx+ix0,0);
 396     ixe = TMath::Min(knx+ix0,GetSeg()->Npx()-1);
 397     izs = TMath::Max(-knz+iz0,0);
 398     ize = TMath::Min(knz+iz0,GetSeg()->Npz()-1);
 399     for(ix=ixs;ix<=ixe;ix++) for(iz=izs;iz<=ize;iz++){
 400         GetSeg()->DetToLocal(ix,iz,x,z); // pixel center
 401         x1  = x;
 402         z1  = z;
 403         x2  = x1 + 0.5*kmictocm*GetSeg()->Dpx(ix); // Upper
 404         x1 -= 0.5*kmictocm*GetSeg()->Dpx(ix);  // Lower
 405         z2  = z1 + 0.5*kmictocm*GetSeg()->Dpz(iz); // Upper
 406         z1 -= 0.5*kmictocm*GetSeg()->Dpz(iz);  // Lower
 407         x1 -= x0; // Distance from where track traveled
 408         x2 -= x0; // Distance from where track traveled
 409         z1 -= z0; // Distance from where track traveled
 410         z2 -= z0; // Distance from where track traveled
 411         s   = 0.25; // Correction based on definision of Erfc
 412         s  *= TMath::Erfc(sp*x1) - TMath::Erfc(sp*x2);
 413 #ifdef DEBUG
 414         cout << "el=" << el << " ix0=" << ix0 << " ix=" << ix << " x0="<< x <<
 415             " iz0=" << iz0 << " iz=" << iz << " z0=" << z  <<
 416             " sp*x1=" << sp*x1 <<" sp*x2=" << sp*x2 << " s=" << s;
 417 #endif
 418         s  *= TMath::Erfc(sp*z1) - TMath::Erfc(sp*z2);
 419 #ifdef DEBUG
 420         cout << " sp*z1=" << sp*z1 <<" sp*z2=" << sp*z2 << " s=" << s << endl;
 421 #endif
 422         GetMap()->AddSignal(iz,ix,t,hi,GetModuleNumber(),s*el);
 423     } // end for ix, iz
 424 }
 425 //______________________________________________________________________
 426 /*
 427 Double_t *AliITSsimulationSPDdubna::CreateFindCellEdges(Double_t x0,
 428                             Double_t x1,Double_t z0,Double_t z1,Int_t &n){
 429     // Note: This function is a potensial source for a memory leak. The memory
 430     // pointed to in its return, must be deleted.
 431     // Inputs:
 432     //    Double_t x0   The starting location of the track step in x
 433     //    Double_t x1   The distance allong x for the track step
 434     //    Double_t z0   The starting location of the track step in z
 435     //    Double_t z1   The distance allong z for the track step
 436     // Output:
 437     //    Int_t &n      The size of the array returned. Minimal n=2.
 438     // Return:
 439     //    The pointer to the array of track steps.
 440     Int_t ix0,ix1,ix,iz0,iz1,iz,i;
 441     Double_t x,z,lx,ux,lz,uz,a,b,c,d;
 442     Double_t *t;
 443
 444     GetSeg()->LocalToDet(x0,z0,ix0,iz0);
 445     GetSeg()->LocalToDet(x1,z1,ix1,iz1);
 446     n      = 2 + TMath::Abs(ix1-ix0) + TMath::Abs(iz1-iz0);
 447     t      = new Double_t[n];
 448     t[0]   = 0.0;
 449     t[n-1] = 1.0;
 450     x      = x0;
 451     z      = z0;
 452     for(i=1;i<n-1;i++){
 453         GetSeg()->LocalToDet(x,z,ix,iz);
 454         GetSeg()->CellBoundries(ix,iz,lx,ux,lz,uz);
 455         a = (lx-x0)/x1;
 456         if(a<=t[i-1]) a = 1.0;
 457         b = (ux-x0)/x1;
 458         if(b<=t[i-1]) b = 1.0;
 459         c = (lz-z0)/z1;
 460         if(c<=t[i-1]) c = 1.0;
 461         d = (uz-z0)/z1;
 462         if(d<=t[i-1]) d = 1.0;
 463         t[i] = TMath::Min(TMath::Min(TMath::Min(a,b),c),d);
 464         x = x0+x1*(t[i]*1.00000001);
 465         z = z0+z1*(t[i]*1.00000001);
 466         i++;
 467     } // end for i
 468     return t;
 469 }
 470 */
 471 //______________________________________________________________________
 472 void AliITSsimulationSPDdubna::HitToSDigitOld(AliITSmodule *mod){
 473     // digitize module Old method
 474     // Inputs:
 475     //    AliITSmodule *mod  Pointer to this module
 476     // Output:
 477     //    none.
 478     // Return:
 479     //    none.
 480     const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons
 481                                       // for 3.6 eV/pair
 482     const Float_t kconv   = 10000.;   // cm -> microns
 483
 484     Float_t spdLength    = GetSeg()->Dz();
 485     Float_t spdWidth     = GetSeg()->Dx();
 486     Float_t spdThickness = GetSeg()->Dy();
 487     Float_t difCoef=GetResp()->SigmaDiffusion1D(1.0);
 488     Float_t zPix0  = 1e+6;
 489     Float_t xPix0  = 1e+6;
 490     Float_t yPrev  = 1e+6;
 491     Float_t zPitch = GetSeg()->Dpz(0);
 492     Float_t xPitch = GetSeg()->Dpx(0);
 493     //  Array of pointers to the label-signal list
 494     Int_t indexRange[4] = {0,0,0,0};
 495
 496     // Fill detector maps with GEANT hits
 497     // loop over hits in the module
 498     static Bool_t first;
 499     Int_t lasttrack=-2;
 500     Int_t hit, iZi, jz, jx;
 501     Int_t idhit=-1; //!
 502     TObjArray *fHits = mod->GetHits();
 503     Int_t nhits      = fHits->GetEntriesFast();
 504     Float_t yPix0 = -spdThickness/2;
 505     if (!nhits) return;
 506 #ifdef DEBUG
 507     cout<<"len,wid,thickness,nx,nz,pitchx,pitchz,difcoef ="<<spdLength<<","
 508         <<spdWidth<<","<<spdThickness<<","<<GetNPixelsX()<<","<<GetNPixelsZ()<<","
 509         <<xPitch<<","<<zPitch<<","<<difCoef<<endl;
 510     cout<<"SPDdubna: module,nhits ="<<GetModuleNumber()<<","<<nhits<<endl;
 511 #endif
 512     for (hit=0;hit<nhits;hit++) {
 513         AliITShit *iHit = (AliITShit*) fHits->At(hit);
 514 #ifdef DEBUG
 515         cout << "Hits=" << hit << "," << *iHit << endl;
 516 #endif
 517         //Int_t layer = iHit->GetLayer();
 518
 519         // work with the idtrack=entry number in the TreeH
 520         //Int_t idhit,idtrack; //!
 521         //mod->GetHitTrackAndHitIndex(hit,idtrack,idhit);  //!
 522         //Int_t idtrack=mod->GetHitTrackIndex(hit);
 523         // or store straight away the particle position in the array
 524         // of particles :
 525         if(iHit->StatusEntering()) idhit=hit;
 526         Int_t itrack = iHit->GetTrack();
 527         Int_t dray = 0;
 528
 529         if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE;
 530
 531         //Int_t parent = iHit->GetParticle()->GetFirstMother();
 532         Int_t partcode = iHit->GetParticle()->GetPdgCode();
 533
 534         //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0,
 535         // 211 - pi+,  310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
 536
 537         Float_t pmod = iHit->GetParticle()->P(); // total momentum at the
 538                                                    // vertex
 539         pmod *= 1000;
 540
 541         if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
 542                                                  // at p < 6 MeV/c
 543
 544         //  Get hit z and x(r*phi) cordinates for each module (detector)
 545         //  in local system.
 546
 547         Float_t zPix = kconv*iHit->GetZL();
 548         Float_t xPix = kconv*iHit->GetXL();
 549         Float_t yPix = kconv*iHit->GetYL();
 550
 551         // Get track status
 552         Int_t status = iHit->GetTrackStatus();
 553
 554         // Check boundaries
 555         if(zPix  > spdLength/2) {
 556 #ifdef DEBUG
 557             cout<<"!!! SPD: z outside ="<<zPix<<endl;
 558 #endif
 559             zPix = spdLength/2 - 10;
 560         }
 561         if(zPix  < 0 && zPix < -spdLength/2) {
 562 #ifdef DEBUG
 563             cout<<"!!! SPD: z outside ="<<zPix<<endl;
 564 #endif
 565             zPix = -spdLength/2 + 10;
 566         }
 567         if(xPix  > spdWidth/2) {
 568 #ifdef DEBUG
 569             cout<<"!!! SPD: x outside ="<<xPix<<endl;
 570 #endif
 571             xPix = spdWidth/2 - 10;
 572         }
 573         if(xPix  < 0 && xPix < -spdWidth/2) {
 574 #ifdef DEBUG
 575             cout<<"!!! SPD: x outside ="<<xPix<<endl;
 576 #endif
 577             xPix = -spdWidth/2 + 10;
 578         }
 579         Int_t trdown = 0;
 580
 581         // enter Si or after event in Si
 582         if (status == 66 ) {
 583             zPix0 = zPix;
 584             xPix0 = xPix;
 585             yPrev = yPix;
 586         } // end if status == 66
 587
 588         Float_t depEnergy = iHit->GetIonization();
 589         // skip if the input point to Si
 590
 591         if(depEnergy <= 0.) continue;
 592
 593         // if track returns to the opposite direction:
 594         if (yPix < yPrev) {
 595             trdown = 1;
 596         } // end if yPix < yPrev
 597
 598         // take into account the holes diffusion inside the Silicon
 599         // the straight line between the entrance and exit points in Si is
 600         // divided into the several steps; the diffusion is considered
 601         // for each end point of step and charge
 602         // is distributed between the pixels through the diffusion.
 603
 604         //  ---------- the diffusion in Z (beam) direction -------
 605         Float_t charge = depEnergy*kEnToEl;         // charge in e-
 606         Float_t drPath = 0.;
 607         Float_t tang = 0.;
 608         Float_t sigmaDif = 0.;
 609         Float_t zdif = zPix - zPix0;
 610         Float_t xdif = xPix - xPix0;
 611         Float_t ydif = TMath::Abs(yPix - yPrev);
 612         Float_t ydif0 = TMath::Abs(yPrev - yPix0);
 613
 614         if(ydif < 1) continue; // ydif is not zero
 615
 616         Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
 617
 618         Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
 619         Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1;
 620
 621         // number of the steps along the track:
 622         Int_t nsteps = ndZ;
 623         if(ndX > ndZ) nsteps = ndX;
 624         if(nsteps < 20) nsteps = 20;  // minimum number of the steps
 625
 626         if (projDif < 5 ) {
 627             drPath = (yPix-yPix0)*1.e-4;
 628             drPath = TMath::Abs(drPath);        // drift path in cm
 629             sigmaDif = difCoef*sqrt(drPath);    // sigma diffusion in cm
 630             sigmaDif = sigmaDif*kconv;         // sigma diffusion in microns
 631             nsteps = 1;
 632         }  // end if projDif < 5
 633
 634         if(projDif > 5) tang = ydif/projDif;
 635         Float_t dCharge = charge/nsteps;       // charge in e- for one step
 636         Float_t dZ = zdif/nsteps;
 637         Float_t dX = xdif/nsteps;
 638
 639         for (iZi = 1; iZi <= nsteps;iZi++) {
 640             Float_t dZn = iZi*dZ;
 641             Float_t dXn = iZi*dX;
 642             Float_t zPixn = zPix0 + dZn;
 643             Float_t xPixn = xPix0 + dXn;
 644
 645             if(projDif >= 5) {
 646                 Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
 647                 drPath = dProjn*tang*1.e-4; // drift path for iZi+1 step in cm
 648                 if(trdown == 0) {
 649                     drPath = TMath::Abs(drPath) + ydif0*1.e-4;
 650                 }// end if trdow ==0
 651                 if(trdown == 1) {
 652                     drPath = ydif0*1.e-4 - TMath::Abs(drPath);
 653                     drPath = TMath::Abs(drPath);
 654                 } // end if trdown == 1
 655                 sigmaDif = difCoef*sqrt(drPath);
 656                 sigmaDif = sigmaDif*kconv;       // sigma diffusion in microns
 657             } // end if projdif >= 5
 658
 659             zPixn = (zPixn + spdLength/2.);
 660             xPixn = (xPixn + spdWidth/2.);
 661             Int_t nZpix, nXpix;
 662             GetSeg()->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
 663             zPitch = GetSeg()->Dpz(nZpix);
 664             GetSeg()->GetPadTxz(xPixn,zPixn);
 665             // set the window for the integration
 666             Int_t jzmin = 1;
 667             Int_t jzmax = 3;
 668             if(nZpix == 1) jzmin =2;
 669             if(nZpix == GetNPixelsZ()) jzmax = 2;
 670
 671             Int_t jxmin = 1;
 672             Int_t jxmax = 3;
 673             if(nXpix == 1) jxmin =2;
 674             if(nXpix == GetNPixelsX()) jxmax = 2;
 675
 676             Float_t zpix = nZpix;
 677             Float_t dZright = zPitch*(zpix - zPixn);
 678             Float_t dZleft = zPitch - dZright;
 679
 680             Float_t xpix = nXpix;
 681             Float_t dXright = xPitch*(xpix - xPixn);
 682             Float_t dXleft = xPitch - dXright;
 683
 684             Float_t dZprev = 0.;
 685             Float_t dZnext = 0.;
 686             Float_t dXprev = 0.;
 687             Float_t dXnext = 0.;
 688
 689             for(jz=jzmin; jz <=jzmax; jz++) {
 690                 if(jz == 1) {
 691                     dZprev = -zPitch - dZleft;
 692                     dZnext = -dZleft;
 693                 } else if(jz == 2) {
 694                     dZprev = -dZleft;
 695                     dZnext = dZright;
 696                 } else if(jz == 3) {
 697                     dZprev = dZright;
 698                     dZnext = dZright + zPitch;
 699                 } // end if jz
 700                 // kz changes from 1 to the fNofPixels(270)
 701                 Int_t kz = nZpix + jz -2;
 702
 703                 Float_t zArg1 = dZprev/sigmaDif;
 704                 Float_t zArg2 = dZnext/sigmaDif;
 705                 Float_t zProb1 = TMath::Erfc(zArg1);
 706                 Float_t zProb2 = TMath::Erfc(zArg2);
 707                 Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge;
 708
 709
 710                 // ----------- holes diffusion in X(r*phi) direction  --------
 711
 712                 if(dZCharge > 1.) {
 713                     for(jx=jxmin; jx <=jxmax; jx++) {
 714                         if(jx == 1) {
 715                             dXprev = -xPitch - dXleft;
 716                             dXnext = -dXleft;
 717                         } else if(jx == 2) {
 718                             dXprev = -dXleft;
 719                             dXnext = dXright;
 720                         } else if(jx == 3) {
 721                             dXprev = dXright;
 722                             dXnext = dXright + xPitch;
 723                         }  // end if jx
 724                         Int_t kx = nXpix + jx -2;
 725                         Float_t xArg1 = dXprev/sigmaDif;
 726                         Float_t xArg2 = dXnext/sigmaDif;
 727                         Float_t xProb1 = TMath::Erfc(xArg1);
 728                         Float_t xProb2 = TMath::Erfc(xArg2);
 729                         Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge;
 730
 731                         if(dXCharge > 1.) {
 732                             if (first) {
 733                                 indexRange[0]=indexRange[1]=kz-1;
 734                                 indexRange[2]=indexRange[3]=kx-1;
 735                                 first=kFALSE;
 736                             } // end if first
 737                             indexRange[0]=TMath::Min(indexRange[0],kz-1);
 738                             indexRange[1]=TMath::Max(indexRange[1],kz-1);
 739                             indexRange[2]=TMath::Min(indexRange[2],kx-1);
 740                             indexRange[3]=TMath::Max(indexRange[3],kx-1);
 741                             // The calling sequence for UpdateMapSignal was
 742                             // moved into the (dx > 1 e-) loop because it
 743                             // needs to call signal which is defined inside
 744                             // this loop
 745                             UpdateMapSignal(kz-1,kx-1,
 746                              ((AliITShit*)(mod->GetHit(hit)))->GetTrack(),
 747                                             hit,dXCharge);
 748                         } // dXCharge > 1 e-
 749                     } // jx loop
 750                 } // dZCharge > 1 e-
 751             } // jz loop
 752         } // iZi loop
 753         if (status == 65) {   // the step is inside of Si
 754             zPix0 = zPix;
 755             xPix0 = xPix;
 756         } // end if status == 65
 757         yPrev = yPix;
 758     }   // hit loop inside the module
 759 }
 760 //______________________________________________________________________
 761 void AliITSsimulationSPDdubna::ChargeToSignal(){
 762     // add noise and electronics, perform the zero suppression and add the
 763     // digit to the list
 764     // Inputs:
 765     //    none.
 766     // Outputs:
 767     //    none.
 768     // Return:
 769     //    none.
 770     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
 771     Float_t threshold = GetResp()->GetThreshold();
 772     Int_t j,ix,iz;
 773     Float_t  electronics;
 774     Double_t sig;
 775     const Int_t    nmaxtrk=AliITSdigitSPD::GetNTracks();
 776     static AliITSdigitSPD dig;
 777
 778     for(iz=0; iz<GetNPixelsZ(); iz++) for(ix=0; ix<GetNPixelsX(); ix++){
 779       if(GetResp()->IsPixelDead(GetModuleNumber(),ix,iz)) continue;
 780       electronics = GetResp()->ApplyBaselineAndNoise();
 781       sig = GetMap()->GetSignalOnly(iz,ix);
 782       UpdateMapNoise(iz,ix,electronics);
 783 #ifdef DEBUG
 784       cout << sig << "+" << electronics <<">threshold=" << threshold  << endl;
 785 #endif
 786       if (sig+electronics <= threshold) continue;
 787       dig.SetCoord1(iz);
 788       dig.SetCoord2(ix);
 789       dig.SetSignal(1);
 790       dig.SetSignalSPD((Int_t) GetMap()->GetSignal(iz,ix));
 791       for(j=0;j<nmaxtrk;j++){
 792         if (j<GetMap()->GetNEnteries()) {
 793           dig.SetTrack(j,GetMap()->GetTrack(iz,ix,j));
 794           dig.SetHit(j,GetMap()->GetHit(iz,ix,j));
 795         }else { // Default values
 796           dig.SetTrack(j,-3);
 797           dig.SetHit(j,-1);
 798         } // end if GetMap()
 799       } // end for j
 800 #ifdef DEBUG
 801       cout<<iz<<","<<ix<<","<<*(GetMap()->GetpListItem(iz,ix)) << endl;
 802 #endif
 803       aliITS->AddSimDigit(0,&dig);
 804     } //  for ix/iz
 805 }
 806 //______________________________________________________________________
 807 void AliITSsimulationSPDdubna::CreateHistograms(){
 808     // create 1D histograms for tests
 809     // Inputs:
 810     //    none.
 811     // Outputs:
 812     //    none.
 813     // Return:
 814     //     none.
 815
 816     printf("SPD - create histograms\n");
 817
 818     fHis=new TObjArray(GetNPixelsZ());
 819     TString spdName("spd_");
 820     for (Int_t i=0;i<GetNPixelsZ();i++) {
 821         Char_t pixelz[4];
 822         sprintf(pixelz,"%d",i);
 823         spdName.Append(pixelz);
 824         fHis->AddAt(new TH1F(spdName.Data(),"SPD maps",
 825                              GetNPixelsX(),0.,(Float_t) GetNPixelsX()), i);
 826     } // end for i
 827 }
 828 //______________________________________________________________________
 829 void AliITSsimulationSPDdubna::ResetHistograms(){
 830     // Reset histograms for this detector
 831     // Inputs:
 832     //    none.
 833     // Outputs:
 834     //    none.
 835     // Return:
 836     //     none.
 837
 838     for ( int i=0;i<GetNPixelsZ();i++ ) {
 839         if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
 840     } // end for i
 841 }
 842
 843 //______________________________________________________________________
 844 void AliITSsimulationSPDdubna::SetCoupling(Int_t row, Int_t col, Int_t ntrack,
 845                                       Int_t idhit) {
 846     //  Take into account the coupling between adiacent pixels.
 847     //  The parameters probcol and probrow are the probability of the
 848     //  signal in one pixel shared in the two adjacent pixels along
 849     //  the column and row direction, respectively.
 850     //Begin_Html
 851     /*
 852       <img src="picts/ITS/barimodel_3.gif">
 853       </pre>
 854       <br clear=left>
 855       <font size=+2 color=red>
 856       <a href="mailto:tiziano.virgili@cern.ch"></a>.
 857       </font>
 858       <pre>
 859     */
 860     //End_Html
 861     // Inputs:
 862     //    Int_t row            z cell index
 863     //    Int_t col            x cell index
 864     //    Int_t ntrack         track incex number
 865     //    Int_t idhit          hit index number
 866     //    Int_t module         module number
 867     // Outputs:
 868     //    none.
 869     // Return:
 870     //     none.
 871     Int_t j1,j2,flag=0;
 872     Double_t pulse1,pulse2;
 873     Float_t couplR=0.0,couplC=0.0;
 874     Double_t xr=0.;
 875     AliITSpList *pList = GetMap();
 876
 877     GetCouplings(couplR,couplC);
 878     j1 = row;
 879     j2 = col;
 880     pulse1 = pList->GetSignalOnly(row,col);
 881     pulse2 = pulse1;
 882     for (Int_t isign=-1;isign<=1;isign+=2){// loop in row direction
 883       do{
 884         j1 += isign;
 885         //   pulse1 *= couplR;
 886         xr = gRandom->Rndm();
 887         //   if ((j1<0) || (j1>GetNPixelsZ()-1) || (pulse1<GetThreshold())){
 888         if ((j1<0) || (j1>GetNPixelsZ()-1) || (xr>couplR)){
 889           j1 = row;
 890           flag = 1;
 891         }else{
 892           UpdateMapSignal(j1,col,ntrack,idhit,pulse1);
 893           //  flag = 0;
 894           flag = 1; // only first next!!
 895         } // end if
 896       } while(flag == 0);
 897       // loop in column direction
 898       do{
 899         j2 += isign;
 900         // pulse2 *= couplC;
 901         xr = gRandom->Rndm();
 902         //  if ((j2<0) || (j2>(GetNPixelsX()-1)) || (pulse2<GetThreshold())){
 903         if ((j2<0) || (j2>GetNPixelsX()-1) || (xr>couplC)){
 904           j2 = col;
 905           flag = 1;
 906         }else{
 907           UpdateMapSignal(row,j2,ntrack,idhit,pulse2);
 908           //  flag = 0;
 909           flag = 1; // only first next!!
 910         } // end if
 911       } while(flag == 0);
 912     } // for isign
 913 }
 914 //______________________________________________________________________
 915 void AliITSsimulationSPDdubna::SetCouplingOld(Int_t row, Int_t col,
 916                 Int_t ntrack,Int_t idhit) {
 917     //  Take into account the coupling between adiacent pixels.
 918     //  The parameters probcol and probrow are the fractions of the
 919     //  signal in one pixel shared in the two adjacent pixels along
 920     //  the column and row direction, respectively.
 921     //Begin_Html
 922     /*
 923       <img src="picts/ITS/barimodel_3.gif">
 924       </pre>
 925       <br clear=left>
 926       <font size=+2 color=red>
 927       <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
 928       </font>
 929       <pre>
 930     */
 931     //End_Html
 932     // Inputs:
 933     //    Int_t row            z cell index
 934     //    Int_t col            x cell index
 935     //    Int_t ntrack         track incex number
 936     //    Int_t idhit          hit index number
 937     //    Int_t module         module number
 938     // Outputs:
 939     //    none.
 940     // Return:
 941     //     none.
 942     Int_t j1,j2,flag=0;
 943     Double_t pulse1,pulse2;
 944     Float_t couplR=0.0,couplC=0.0;
 945     AliITSpList *pList = GetMap();
 946
 947     GetCouplings(couplR,couplC);
 948     j1 = row;
 949     j2 = col;
 950     pulse1 = pList->GetSignalOnly(row,col);
 951     pulse2 = pulse1;
 952     for (Int_t isign=-1;isign<=1;isign+=2){// loop in row direction
 953       do{
 954         j1 += isign;
 955         pulse1 *= couplR;
 956         if ((j1<0) || (j1>GetNPixelsZ()-1) || (pulse1<GetThreshold())){
 957           pulse1 = pList->GetSignalOnly(row,col);
 958           j1 = row;
 959           flag = 1;
 960         }else{
 961           UpdateMapSignal(j1,col,ntrack,idhit,pulse1);
 962           flag = 0;
 963         } // end if
 964       } while(flag == 0);
 965       // loop in column direction
 966       do{
 967         j2 += isign;
 968         pulse2 *= couplC;
 969         if ((j2<0) || (j2>(GetNPixelsX()-1)) || (pulse2<GetThreshold())){
 970           pulse2 = pList->GetSignalOnly(row,col);
 971           j2 = col;
 972           flag = 1;
 973         }else{
 974           UpdateMapSignal(row,j2,ntrack,idhit,pulse2);
 975           flag = 0;
 976         } // end if
 977       } while(flag == 0);
 978     } // for isign
 979 }