TObjArray replaced with TClonesArray in AliITSpList - Typo corrected
[u/mrichter/AliRoot.git] / ITS / AliITSsimulationSPDdubna.cxx
index afa479f..bf76b68 100644 (file)
-#include <iostream.h>
-#include <TRandom.h>
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+/*
+$Id$
+*/
+
+#include <Riostream.h>
 #include <TH1.h>
 #include <TMath.h>
-#include <TString.h>
 #include <TParticle.h>
-
-
-#include "AliRun.h"
+#include <TRandom.h>
+#include <TString.h>
 #include "AliITS.h"
+#include "AliITSMapA2.h" 
+#include "AliITSdigitSPD.h"
+#include "AliITSgeom.h"
 #include "AliITShit.h"
-#include "AliITSdigit.h"
 #include "AliITSmodule.h"
-#include "AliITSMapA2.h" 
 #include "AliITSpList.h"
+#include "AliITSCalibrationSPD.h"
+#include "AliITSsegmentationSPD.h"
 #include "AliITSsimulationSPDdubna.h"
-#include "AliITSsegmentation.h"
-#include "AliITSresponse.h"
-
-
+#include "AliLog.h"
+#include "AliRun.h"
 
+//#define DEBUG
 
 ClassImp(AliITSsimulationSPDdubna)
 ////////////////////////////////////////////////////////////////////////
+// Version: 1
+// Modified by Bjorn S. Nilsen
 // Version: 0
 // Written by Boris Batyunya
 // December 20 1999
 //
-// AliITSsimulationSPDdubna is the simulation of SPDs
+// AliITSsimulationSPDdubna is to do the simulation of SPDs.
 //______________________________________________________________________
-
-
-AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(){
-    // constructor
-
-    fResponse = 0;
-    fSegmentation = 0;
-    fMapA2 = 0;
-    fpList = 0;
-    fModule = 0;
-    fEvent = 0;
-    fHis = 0;
-    fNoise = 0.;
-    fBaseline = 0.;
-    fNPixelsZ = 0;
-    fNPixelsX = 0;
+AliITSsimulationSPDdubna::AliITSsimulationSPDdubna():
+AliITSsimulation(),
+fHis(0),
+fSPDname(),
+fCoupling(0){
+    // Default constructor.
+    // Inputs:
+    //    none.
+    // Outputs:
+    //    none.
+    // Return:
+    //    A default constructed AliITSsimulationSPDdubna class.
+
+    AliDebug(1,Form("Calling degault constructor"));
 }
 //______________________________________________________________________
-AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(AliITSsegmentation *seg,
-                                                  AliITSresponse *resp){
+AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(AliITSDetTypeSim *dettyp, Int_t cup):
+AliITSsimulation(dettyp),
+fHis(0),
+fSPDname(),
+fCoupling(cup){
     // standard constructor
-
-    fHis = 0;
-    fResponse = resp;
-    fSegmentation = seg;
-    fModule = 0;
-    fEvent = 0;
-
-    fNPixelsZ=fSegmentation->Npz();
-    fNPixelsX=fSegmentation->Npx();
-
-    fResponse->GetNoiseParam(fNoise,fBaseline);
-
-    fMapA2 = new AliITSMapA2(fSegmentation);
-
-    fpList = new AliITSpList(fNPixelsZ+1,fNPixelsX+1);
-
+    // Inputs:
+    //    AliITSsegmentation *seg  A pointer to the segmentation class
+    //                             to be used for this simulation
+    //    AliITSCalibration     *resp A pointer to the responce class to
+    //                             be used for this simulation
+    //    Int_t              cup   The type of coupling to be used
+    //                             =1 uses SetCoupling, =2 uses SetCouplingOld
+    //                             With diffusion tured off
+    //                             =3 uses SetCoupling, =4 uses SetCouplingOld
+    //                             with diffusion on other, no coupling.
+    // Outputs:
+    //    none.
+    // Return:
+    //    A default constructed AliITSsimulationSPDdubna class.
+
+    AliDebug(1,
+            Form("Calling degault constructor cup=%d",cup));
+    if(cup==1||cup==2){ // For the moment, remove defusion if Coupling is
+        // set.
+      AliITSCalibrationSPD* res = (AliITSCalibrationSPD*)GetCalibrationModel(fDetType->GetITSgeom()->GetStartSPD());
+      res->SetTemperature(0.0);
+      res->SetDistanceOverVoltage(0.0);
+    } // end if
+    Init();
+}
+//______________________________________________________________________
+void AliITSsimulationSPDdubna::Init(){
+    // Initilization
+    // Inputs:
+    //    none.
+    // Outputs:
+    //    none.
+    // Return:
+    //    none.
+    const Double_t kmictocm = 1.0e-4; // convert microns to cm.
+
+    SetModuleNumber(0);
+    SetEventNumber(0);
+    SetMap(new AliITSpList(GetNPixelsZ(),GetNPixelsX()));
+    AliITSCalibrationSPD* res = (AliITSCalibrationSPD*)GetCalibrationModel(fDetType->GetITSgeom()->GetStartSPD());
+    AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0);
+    res->SetDistanceOverVoltage(kmictocm*seg->Dy(),50.0);
 }
 //______________________________________________________________________
 AliITSsimulationSPDdubna::~AliITSsimulationSPDdubna(){
     // destructor
-
-    delete fMapA2;
+    // Inputs:
+    //    none.
+    // Outputs:
+    //    none.
+    // Return:
+    //     none.
 
     if (fHis) {
-       fHis->Delete(); 
-       delete fHis;     
+        fHis->Delete(); 
+        delete fHis;     
     } // end if fHis
 }
 //______________________________________________________________________
 AliITSsimulationSPDdubna::AliITSsimulationSPDdubna(const 
                                                   AliITSsimulationSPDdubna 
-                                                  &source){
-    //     Copy Constructor 
-    if(&source == this) return;
-    this->fMapA2 = source.fMapA2;
-    this->fNoise = source.fNoise;
-    this->fBaseline = source.fBaseline;
-    this->fNPixelsX = source.fNPixelsX;
-    this->fNPixelsZ = source.fNPixelsZ;
-    this->fHis = source.fHis;
+                                                  &s) : AliITSsimulation(s){
+    //     Copy Constructor
+    // Inputs:
+    //    AliITSsimulationSPDdubna &s The original class for which
+    //                                this class is a copy of
+    // Outputs:
+    //    none.
+    // Return:
+
+    *this = s;
     return;
 }
 //______________________________________________________________________
 AliITSsimulationSPDdubna&  AliITSsimulationSPDdubna::operator=(const 
-                                           AliITSsimulationSPDdubna &source){
+                                           AliITSsimulationSPDdubna &s){
+    //    Assignment operator
+    // Inputs:
+    //    AliITSsimulationSPDdubna &s The original class for which
+    //                                this class is a copy of
+    // Outputs:
+    //    none.
+    // Return:
+
+    if(&s == this) return *this;
+    this->fHis = s.fHis;
+    fCoupling  = s.fCoupling;
+    fSPDname   = s.fSPDname;
+    return *this;
+}
+//______________________________________________________________________
+AliITSsimulation&  AliITSsimulationSPDdubna::operator=(const 
+                                           AliITSsimulation &s){
     //    Assignment operator
-    if(&source == this) return *this;
-    this->fMapA2 = source.fMapA2;
-    this->fNoise = source.fNoise;
-    this->fBaseline = source.fBaseline;
-    this->fNPixelsX = source.fNPixelsX;
-    this->fNPixelsZ = source.fNPixelsZ;
-    this->fHis = source.fHis;
+    // Inputs:
+    //    AliITSsimulationSPDdubna &s The original class for which
+    //                                this class is a copy of
+    // Outputs:
+    //    none.
+    // Return:
+
+    if(&s == this) return *this;
+    Error("AliITSsimulationSPDdubna","Not allowed to make a = with "
+          "AliITSsimulationSPDdubna","Using default creater instead");
+
     return *this;
 }
 //______________________________________________________________________
 void AliITSsimulationSPDdubna::InitSimulationModule(Int_t module, Int_t event){
     //  This function creates maps to build the list of tracks for each
-    //  summable digit.
-    //
+    //  summable digit. Inputs defined by base class.
     //  Inputs:
     //    Int_t module   // Module number to be simulated
     //    Int_t event    // Event number to be simulated
-    //
     //  Outputs:
     //    none
-    //
     //  Returns:
     //    none
 
-    fModule = module;
-    fEvent  = event;
-    fMapA2->ClearMap();
-    fpList->ClearMap();
+    AliDebug(1,Form("(module=%d,event=%d)",module,event));
+    SetModuleNumber(module);
+    SetEventNumber(event);
+    ClearMap();
 }
 //_____________________________________________________________________
-void AliITSsimulationSPDdubna::SDigitiseModule(AliITSmodule *mod, Int_t mask,
-                                              Int_t event){
-    //  This function begins the work of creating S-Digits
-    //
+void AliITSsimulationSPDdubna::SDigitiseModule(AliITSmodule *mod,Int_t,
+                                               Int_t event){
+    //  This function begins the work of creating S-Digits.  Inputs defined
+    //  by base class.
     //  Inputs:
     //    AliITSmodule *mod  //  module
-    //    Int_t mask         //  mask to be applied to the module
-    //
+    //    Int_t              //  not used
+    //    Int_t event        //  Event number
     //  Outputs:
     //    none
-    //
     //  Return:
     //    test              //  test returns kTRUE if the module contained hits
     //                      //  test returns kFALSE if it did not contain hits
 
-    Int_t module = 0;
-
-    if(!(mod->GetNhits())) return;// if module has no hits don't create Sdigits
-    fModule = mod->GetIndex();
-    HitToSDigit(mod, module, mask, fpList);
-    WriteSDigits(fpList);
-    fMapA2->ClearMap();
-    fpList->ClearMap();
+    AliDebug(1,Form("(mod=%p, ,event=%d)",mod,event));
+    if(!(mod->GetNhits())){
+        AliDebug(1,Form("In event %d module %d there are %d hits returning.",
+                       event, mod->GetIndex(),mod->GetNhits()));
+        return;// if module has no hits don't create Sdigits
+    } // end if
+    SetModuleNumber(mod->GetIndex());
+    SetEventNumber(event);
+    HitToSDigit(mod);
+    WriteSDigits();
+    ClearMap();
 }
 //______________________________________________________________________
-void AliITSsimulationSPDdubna::WriteSDigits(AliITSpList *pList){
+void AliITSsimulationSPDdubna::WriteSDigits(){
     //  This function adds each S-Digit to pList
-    //
     //  Inputs:
-    //    AliITSpList *pList
-    //
+    //    none.
     //  Outputs:
-    //    none
-    //
+    //    none.
     //  Return:
     //    none
     Int_t ix, nix, iz, niz;
     static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
 
-    pList->GetMaxMapIndex(niz, nix);
+    AliDebug(1,Form("Writing SDigits for module %d",GetModuleNumber()));
+    GetMap()->GetMaxMapIndex(niz, nix);
     for(iz=0; iz<niz; iz++)for(ix=0; ix<nix; ix++){
-       if(pList->GetSignalOnly(iz,ix)>0.0){
-           aliITS->AddSumDigit(*(pList->GetpListItem(iz,ix)));
-       } // end if pList
+        if(GetMap()->GetSignalOnly(iz,ix)>0.0){
+            aliITS->AddSumDigit(*(GetMap()->GetpListItem(iz,ix)));
+           if(AliDebugLevel()>0) {
+             AliDebug(1,Form("%d, %d",iz,ix));
+             cout << *(GetMap()->GetpListItem(iz,ix)) << endl;
+            } // end if GetDebug
+        } // end if GetMap()->GetSignalOnly(iz,ix)>0.0
     } // end for iz,ix
     return; 
 }
 //______________________________________________________________________
 void AliITSsimulationSPDdubna::FinishSDigitiseModule(){
     //  This function calls SDigitsToDigits which creates Digits from SDigits
-    //
     //  Inputs:
     //    none
-    //
     //  Outputs:
     //    none
     //  Return
     //    none
 
-    SDigitsToDigits(fModule, fpList);
+    AliDebug(1,"()");
+    pListToDigits(); // Charge To Signal both adds noise and
+    ClearMap();
     return;
 }
 //______________________________________________________________________
-void AliITSsimulationSPDdubna::SDigitsToDigits(Int_t module,
-                                              AliITSpList *pList){
-    //  This function adds electronic noise to the S-Digits and then adds them
-    // to  a new pList
-    //
-    //  Inputs:
-    //    Int_t       module  // module number
-    //    AliITSpList *pList  // pList
-    //
-    //  Outputs:
-    //    pList is passed along to the functions ChargeToSignal and GetList
-    //
-    //  Return:
-    //    none
-
-    fModule = module;
-    ChargeToSignal(pList); // Charge To Signal both adds noise and
-    fMapA2->ClearMap();
-    pList->ClearMap();
-}
-//______________________________________________________________________
-void AliITSsimulationSPDdubna::DigitiseModule(AliITSmodule *mod, Int_t module,
-                                             Int_t dummy){
+void AliITSsimulationSPDdubna::DigitiseModule(AliITSmodule *mod,Int_t,
+                                              Int_t){
     //  This function creates Digits straight from the hits and then adds
     //  electronic noise to the digits before adding them to pList
-    //
-    //  Inputs:
-    //    AliITSmodule *mod    // module
-    //    Int_t        module  // module number  Dummy.
-    //    Int_t        dummy
-    //
-    //  Outputs:
-    //    Each of the input variables is passed along to HitToSDigit
-    //
-    //  Return:
-    //    none
-
-    fModule = mod->GetIndex();  //This calls the module for HitToSDigit
-    HitToSDigit(mod,fModule, dummy, fpList);
-    ChargeToSignal(fpList);
-    fMapA2->ClearMap();
-    fpList->ClearMap();
-}
-//______________________________________________________________________
-void AliITSsimulationSPDdubna::UpdateMapSignal(Int_t iz, Int_t ix, Int_t trk,
-                                              Int_t ht, Int_t module,
-                                              Double_t signal,
-                                              AliITSpList *pList){
-    //  This function adds a signal to the pList from the pList class
-    //
-    //  Inputs:
-    //    Int_t       iz     // row number
-    //    Int_t       ix     // column number
-    //    Int_t       trk    // track number
-    //    Int_t       ht     // hit number
-    //    Double_t    signal // signal strength
-    //    AliITSpList *pList // pList
-    //
-    //  Outputs:
-    //    All of the inputs are passed to AliITSpList::AddSignal
-    //    Int_t    ix  // row number
-    //    Int_t    iz  // column number
-    //    Double_t sig // signal strength
-    //          // These three variables are defined to preserve the
-    //          // assignments used in the function AliITSMapA2::AddSignal
-    //
-    //  Return:
-    //    none
-
-    fMapA2->AddSignal(iz, ix, signal);
-    pList->AddSignal(iz,ix, trk, ht, fModule, signal);
-}
-//______________________________________________________________________
-void AliITSsimulationSPDdubna::UpdateMapNoise(Int_t iz,
-                                             Int_t ix, Int_t fModule,
-                                             Double_t sig, Float_t noise,
-                                             AliITSpList *pList){
-    //  This function adds noise to data in the MapA2 as well as the pList
-    //
+    //  Each of the input variables is passed along to HitToSDigit
     //  Inputs:
-    //    Int_t       iz       // row number
-    //    Int_t       ix       // column number
-    //    Int_t       mod     // module number
-    //    Double_t    sig     // signal strength
-    //    Double_t    noise   // electronic noise generated by ChargeToSignal
-    //    AliITSpList *pList  // pList
-    //
+    //    AliITSmodule *mod     module
+    //    Int_t                 Dummy.
+    //    Int_t                 Dummy
     //  Outputs:
-    //    All of the inputs are passed to AliITSMapA2::AddSignal or
-    //    AliITSpList::AddNoise
-    //
+    //     none.
     //  Return:
-    //    none
+    //    none.
 
-    fMapA2->AddSignal(iz, ix, sig);
-    pList->AddNoise(iz,ix, fModule, noise);
+    AliDebug(1,Form("(mod=%p,,)",mod));
+    HitToSDigit(mod);
+    pListToDigits();
+    ClearMap();
 }
 //______________________________________________________________________
-void AliITSsimulationSPDdubna::HitToDigit(AliITSmodule *mod, Int_t module,
-                                         Int_t dummy){
-    DigitiseModule(mod, module, dummy);
+void AliITSsimulationSPDdubna::HitToSDigit(AliITSmodule *mod){
+    // Does the charge distributions using Gaussian diffusion charge charing.
+    // Inputs:
+    //    AliITSmodule *mod  Pointer to this module
+    // Output:
+    //    none.
+    // Return:
+    //    none.
+    const Double_t kmictocm = 1.0e-4; // convert microns to cm.
+    TObjArray *hits = mod->GetHits();
+    Int_t nhits = hits->GetEntriesFast();
+    Int_t h,ix,iz,i;
+    Int_t idtrack;
+    Double_t x0=0.0,x1=0.0,y0=0.0,y1=0.0,z0=0.0,z1=0.0,de=0.0;
+    Double_t x,y,z,t,tp,st,dt=0.2,el,sig;
+    AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0);
+    AliITSCalibrationSPD* res = (AliITSCalibrationSPD*)GetCalibrationModel(fDetType->GetITSgeom()->GetStartSPD());
+    Double_t thick = kmictocm*seg->Dy();
+
+    AliDebug(1,Form("(mod=%p) fCoupling=%d",mod,fCoupling));
+    if(nhits<=0) return;
+    for(h=0;h<nhits;h++){
+      if(AliDebugLevel()>0) {
+       AliDebug(1,Form("Hits, %d", h));
+       cout << *(mod->GetHit(h)) << endl;
+      } // end if GetDebug
+        if(!mod->LineSegmentL(h,x0,x1,y0,y1,z0,z1,de,idtrack)) continue;
+        st = TMath::Sqrt(x1*x1+y1*y1+z1*z1);
+        if(st>0.0){
+            st = (Double_t)((Int_t)(st/kmictocm)); // number of microns
+            if(st<=1.0) st = 1.0;
+            dt = 1.0/st;
+            for(t=0.0;t<1.0;t+=dt){ // Integrate over t
+                tp  = t+0.5*dt;
+                x   = x0+x1*tp;
+                y   = y0+y1*tp;
+                z   = z0+z1*tp;
+                if(!(seg->LocalToDet(x,z,ix,iz))) continue; // outside
+                el  = res->GeVToCharge((Double_t)(dt*de));
+                if(GetDebug(1)){
+                    if(el<=0.0) cout<<"el="<<el<<" dt="<<dt
+                                    <<" de="<<de<<endl;
+                } // end if GetDebug
+                sig = res->SigmaDiffusion1D(thick + y);
+                SpreadCharge(x,z,ix,iz,el,sig,idtrack,h);
+            } // end for t
+        } else { // st == 0.0 deposit it at this point
+            x   = x0;
+            y   = y0;
+            z   = z0;
+            if(!(seg->LocalToDet(x,z,ix,iz))) continue; // outside
+            el  = res->GeVToCharge((Double_t)de);
+            sig = res->SigmaDiffusion1D(thick + y);
+            SpreadCharge(x,z,ix,iz,el,sig,idtrack,h);
+        } // end if st>0.0
+        // Coupling
+        switch (fCoupling) {
+        default:
+            break;
+        case 1: case 3:
+            // x is column and z is row (see AliITSsegmentationSPD::GetPadIxz)
+            for(i=0;i<GetMap()->GetEntries();i++) 
+                if(GetMap()->GetpListItem(i)==0) continue;
+                else{
+                    GetMap()->GetMapIndex(
+                              GetMap()->GetpListItem(i)->GetIndex(),iz,ix);
+                    SetCoupling(iz,ix,idtrack,h);
+                } // end for i
+            break;
+        case 2: case 4:
+            // x is column and z is row (see AliITSsegmentationSPD::GetPadIxz)
+            for(i=0;i<GetMap()->GetEntries();i++) 
+                if(GetMap()->GetpListItem(i)==0) continue;
+                else{
+                    GetMap()->GetMapIndex(
+                                GetMap()->GetpListItem(i)->GetIndex(),iz,ix);
+                    SetCouplingOld(iz,ix,idtrack,h);
+                } // end for i
+            break;
+        } // end switch
+    } // Loop over all hits h
+    if(GetDebug(2))Info("HitToSDigit","Finished fCoupling=%d",fCoupling);
 }
 //______________________________________________________________________
-void AliITSsimulationSPDdubna::HitToSDigit(AliITSmodule *mod, Int_t module,
-                                           Int_t dummy, AliITSpList *pList){
-    // digitize module 
-    const Float_t kEnToEl = 2.778e+8; // GeV->charge in electrons 
-                                      // for 3.6 eV/pair 
-    const Float_t kconv = 10000.;     // cm -> microns
-
-    Float_t spdLength = fSegmentation->Dz();
-    Float_t spdWidth = fSegmentation->Dx();
-    Float_t spdThickness = fSegmentation->Dy();
-    Float_t difCoef, dum;       
-    fResponse->DiffCoeff(difCoef,dum); 
-    if(spdThickness > 290) difCoef = 0.00613;  
-
-    Float_t zPix0 = 1e+6;
-    Float_t xPix0 = 1e+6;
-    Float_t yPrev = 1e+6;   
-
-    Float_t zPitch = fSegmentation->Dpz(0);
-    Float_t xPitch = fSegmentation->Dpx(0);
-  
-    TObjArray *fHits = mod->GetHits();
-    module = mod->GetIndex();
-    Int_t nhits = fHits->GetEntriesFast();
-    if (!nhits) return;
-
-    cout<<"len,wid,thickness,nx,nz,pitchx,pitchz,difcoef ="<<spdLength<<","
-       <<spdWidth<<","<<spdThickness<<","<<fNPixelsX<<","<<fNPixelsZ<<","
-       <<xPitch<<","<<zPitch<<","<<difCoef<<endl;
-    //  Array of pointers to the label-signal list
-    Int_t indexRange[4] = {0,0,0,0};
-
-    // Fill detector maps with GEANT hits
-    // loop over hits in the module
-    static Bool_t first;
-    Int_t lasttrack=-2;
-    Int_t hit, iZi, jz, jx;
-    Int_t idhit=-1; //!
-    cout<<"SPDdubna: module,nhits ="<<module<<","<<nhits<<endl;
-    for (hit=0;hit<nhits;hit++) {
-        AliITShit *iHit = (AliITShit*) fHits->At(hit);
-       //Int_t layer = iHit->GetLayer();
-        Float_t yPix0 = -spdThickness/2; 
-
-       // work with the idtrack=entry number in the TreeH
-       //Int_t idhit,idtrack; //!
-       //mod->GetHitTrackAndHitIndex(hit,idtrack,idhit);  //!    
-       //Int_t idtrack=mod->GetHitTrackIndex(hit);  
-        // or store straight away the particle position in the array
-       // of particles : 
-       if(iHit->StatusEntering()) idhit=hit;
-        Int_t itrack = iHit->GetTrack();
-        Int_t dray = 0;
-   
-       if (lasttrack != itrack || hit==(nhits-1)) first = kTRUE; 
-
-       //Int_t parent = iHit->GetParticle()->GetFirstMother();
-        Int_t partcode = iHit->GetParticle()->GetPdgCode();
-
-       //  partcode (pdgCode): 11 - e-, 13 - mu-, 22 - gamma, 111 - pi0,
-       // 211 - pi+,  310 - K0s, 321 - K+, 2112 - n, 2212 - p, 3122 - lambda
-
-       Float_t pmod = iHit->GetParticle()->P(); // total momentum at the
-                                                  // vertex
-        pmod *= 1000;
-
-        if(partcode == 11 && pmod < 6) dray = 1; // delta ray is e-
-                                                 // at p < 6 MeV/c
-
-       //  Get hit z and x(r*phi) cordinates for each module (detector)
-       //  in local system.
-
-       Float_t zPix = kconv*iHit->GetZL();
-       Float_t xPix = kconv*iHit->GetXL();
-       Float_t yPix = kconv*iHit->GetYL();
-
-       // Get track status
-       Int_t status = iHit->GetTrackStatus();      
-
-       // Check boundaries
-       if(zPix  > spdLength/2) {
-           //cout<<"!!! SPD: z outside ="<<zPix<<endl;
-           zPix = spdLength/2 - 10;
-       }
-       if(zPix  < 0 && zPix < -spdLength/2) {
-           //cout<<"!!! SPD: z outside ="<<zPix<<endl;
-           zPix = -spdLength/2 + 10;
-       }
-       if(xPix  > spdWidth/2) {
-           //cout<<"!!! SPD: x outside ="<<xPix<<endl;
-           xPix = spdWidth/2 - 10;
-       }
-       if(xPix  < 0 && xPix < -spdWidth/2) {
-           //cout<<"!!! SPD: x outside ="<<xPix<<endl;
-           xPix = -spdWidth/2 + 10;
-       }
-       Int_t trdown = 0;
-
-       // enter Si or after event in Si
-       if (status == 66 ) {  
-           zPix0 = zPix;
-           xPix0 = xPix;
-           yPrev = yPix; 
-       } // end if status == 66
-
-       Float_t depEnergy = iHit->GetIonization();
-       // skip if the input point to Si       
-
-       if(depEnergy <= 0.) continue;        
-
-       // if track returns to the opposite direction:
-       if (yPix < yPrev) {
-            trdown = 1;
-       } // end if yPix < yPrev
-
-       // take into account the holes diffusion inside the Silicon
-       // the straight line between the entrance and exit points in Si is
-       // divided into the several steps; the diffusion is considered 
-       // for each end point of step and charge
-       // is distributed between the pixels through the diffusion.
-
-       //  ---------- the diffusion in Z (beam) direction -------
-       Float_t charge = depEnergy*kEnToEl;         // charge in e-
-       Float_t drPath = 0.;   
-       Float_t tang = 0.;
-       Float_t sigmaDif = 0.; 
-       Float_t zdif = zPix - zPix0;
-       Float_t xdif = xPix - xPix0;
-       Float_t ydif = TMath::Abs(yPix - yPrev);
-       Float_t ydif0 = TMath::Abs(yPrev - yPix0);
-
-       if(ydif < 1) continue; // ydif is not zero
-
-       Float_t projDif = sqrt(xdif*xdif + zdif*zdif);
-
-       Int_t ndZ = (Int_t)TMath::Abs(zdif/zPitch) + 1;
-       Int_t ndX = (Int_t)TMath::Abs(xdif/xPitch) + 1; 
-
-       // number of the steps along the track:
-       Int_t nsteps = ndZ;
-       if(ndX > ndZ) nsteps = ndX;
-       if(nsteps < 20) nsteps = 20;  // minimum number of the steps 
-
-       if (projDif < 5 ) {
-           drPath = (yPix-yPix0)*1.e-4;  
-           drPath = TMath::Abs(drPath);        // drift path in cm
-           sigmaDif = difCoef*sqrt(drPath);    // sigma diffusion in cm
-           sigmaDif = sigmaDif*kconv;         // sigma diffusion in microns
-           nsteps = 1;
-       }  // end if projDif < 5
-
-       if(projDif > 5) tang = ydif/projDif;
-       Float_t dCharge = charge/nsteps;       // charge in e- for one step
-       Float_t dZ = zdif/nsteps;
-       Float_t dX = xdif/nsteps;
-
-       for (iZi = 1; iZi <= nsteps;iZi++) {
-           Float_t dZn = iZi*dZ;
-           Float_t dXn = iZi*dX;
-           Float_t zPixn = zPix0 + dZn;
-           Float_t xPixn = xPix0 + dXn;
-
-           if(projDif >= 5) {
-               Float_t dProjn = sqrt(dZn*dZn+dXn*dXn);
-               drPath = dProjn*tang*1.e-4; // drift path for iZi+1 step in cm 
-               if(trdown == 0) {
-                   drPath = TMath::Abs(drPath) + ydif0*1.e-4;
-               }// end if trdow ==0
-               if(trdown == 1) {
-                   drPath = ydif0*1.e-4 - TMath::Abs(drPath);
-                   drPath = TMath::Abs(drPath);
-               } // end if trdown == 1
-               sigmaDif = difCoef*sqrt(drPath);    
-               sigmaDif = sigmaDif*kconv;       // sigma diffusion in microns
-           } // end if projdif >= 5
-
-           zPixn = (zPixn + spdLength/2.);  
-           xPixn = (xPixn + spdWidth/2.);  
-            Int_t nZpix, nXpix;
-            fSegmentation->GetPadIxz(xPixn,zPixn,nXpix,nZpix);
-           zPitch = fSegmentation->Dpz(nZpix);
-            fSegmentation->GetPadTxz(xPixn,zPixn);
-           // set the window for the integration
-           Int_t jzmin = 1;  
-           Int_t jzmax = 3; 
-           if(nZpix == 1) jzmin =2;
-           if(nZpix == fNPixelsZ) jzmax = 2; 
-
-           Int_t jxmin = 1;  
-           Int_t jxmax = 3; 
-           if(nXpix == 1) jxmin =2;
-           if(nXpix == fNPixelsX) jxmax = 2; 
-
-           Float_t zpix = nZpix; 
-           Float_t dZright = zPitch*(zpix - zPixn);
-           Float_t dZleft = zPitch - dZright;
-
-           Float_t xpix = nXpix; 
-           Float_t dXright = xPitch*(xpix - xPixn);
-           Float_t dXleft = xPitch - dXright;
-
-           Float_t dZprev = 0.;
-           Float_t dZnext = 0.;
-           Float_t dXprev = 0.;
-           Float_t dXnext = 0.;
-
-           for(jz=jzmin; jz <=jzmax; jz++) {
-               if(jz == 1) {
-                   dZprev = -zPitch - dZleft;
-                   dZnext = -dZleft;
-               } else if(jz == 2) {
-                   dZprev = -dZleft;
-                   dZnext = dZright;
-               } else if(jz == 3) {
-                   dZprev = dZright;
-                   dZnext = dZright + zPitch;
-               } // end if jz
-               // kz changes from 1 to the fNofPixels(270)  
-               Int_t kz = nZpix + jz -2; 
-
-               Float_t zArg1 = dZprev/sigmaDif;
-               Float_t zArg2 = dZnext/sigmaDif;
-               Float_t zProb1 = TMath::Erfc(zArg1);
-               Float_t zProb2 = TMath::Erfc(zArg2);
-               Float_t dZCharge =0.5*(zProb1-zProb2)*dCharge; 
-
-
-               // ----------- holes diffusion in X(r*phi) direction  --------
-
-               if(dZCharge > 1.) { 
-                   for(jx=jxmin; jx <=jxmax; jx++) {
-                       if(jx == 1) {
-                           dXprev = -xPitch - dXleft;
-                           dXnext = -dXleft;
-                       } else if(jx == 2) {
-                           dXprev = -dXleft;
-                           dXnext = dXright;
-                       } else if(jx == 3) {
-                           dXprev = dXright;
-                           dXnext = dXright + xPitch;
-                       }  // end if jx
-                       Int_t kx = nXpix + jx -2;
-                       Float_t xArg1 = dXprev/sigmaDif;
-                       Float_t xArg2 = dXnext/sigmaDif;
-                       Float_t xProb1 = TMath::Erfc(xArg1);
-                       Float_t xProb2 = TMath::Erfc(xArg2);
-                       Float_t dXCharge =0.5*(xProb1-xProb2)*dZCharge;
-
-                       if(dXCharge > 1.) {
-                           Int_t index = kz-1;
-                           if (first) {
-                               indexRange[0]=indexRange[1]=index;
-                               indexRange[2]=indexRange[3]=kx-1;
-                               first=kFALSE;
-                           } // end if first
-                           indexRange[0]=TMath::Min(indexRange[0],kz-1);
-                           indexRange[1]=TMath::Max(indexRange[1],kz-1);
-                           indexRange[2]=TMath::Min(indexRange[2],kx-1);
-                           indexRange[3]=TMath::Max(indexRange[3],kx-1);
-
-                           // build the list of digits for this module 
-                           Double_t signal=fMapA2->GetSignal(index,kx-1);
-                           signal+=dXCharge;
-                           fMapA2->SetHit(index,kx-1,(double)signal);
-
-                           // The calling sequence for UpdateMapSignal was 
-                           // moved into the (dx > 1 e-) loop because it 
-                           // needs to call signal which is defined inside 
-                           // this loop
-                           fModule   = module;//Defined because functions 
-                                              // called by UpdateMapSignal 
-                                              // expect module to be an 
-                                              // integer
-                           UpdateMapSignal(kz-1,kx-1,
-                                           mod->GetHitTrackIndex(hit),
-                                           hit,fModule,signal,pList);
-                       }      // dXCharge > 1 e-
-                   }       // jx loop
-               }       // dZCharge > 1 e-
-           }        // jz loop
-       }         // iZi loop
-        if (status == 65) {   // the step is inside of Si
-           zPix0 = zPix;
-           xPix0 = xPix;
-       } // end if status == 65
-       yPrev = yPix;
-    }   // hit loop inside the module
+void AliITSsimulationSPDdubna::SpreadCharge(Double_t x0,Double_t z0,
+                                            Int_t ix0,Int_t iz0,
+                                           Double_t el,Double_t sig,Int_t t,
+                                           Int_t hi){
+    // Spreads the charge over neighboring cells. Assume charge is distributed
+    // as charge(x,z) = (el/2*pi*sig*sig)*exp(-arg)
+    // arg=((x-x0)*(x-x0)/2*sig*sig)+((z-z0*z-z0)/2*sig*sig)
+    // Defined this way, the integral over all x and z is el.
+    // Inputs:
+    //    Double_t x0   x position of point where charge is liberated
+    //    Double_t y0   y position of point where charge is liberated
+    //    Double_t z0   z position of point where charge is liberated
+    //    Int_t    ix0  row of cell corresponding to point x0
+    //    Int_t    iz0  columb of cell corresponding to point z0
+    //    Double_t el   number of electrons liberated in this step
+    //    Double_t sig  Sigma difusion for this step (y0 dependent)
+    //    Int_t    t    track number
+    //    Int_t    ti   hit track index number
+    //    Int_t    hi   hit "hit" index number
+    // Outputs:
+    //     none.
+    // Return:
+    //     none.
+    const Int_t knx = 3,knz = 2;
+    const Double_t kRoot2 = 1.414213562; // Sqrt(2).
+    const Double_t kmictocm = 1.0e-4; // convert microns to cm.
+    Int_t ix,iz,ixs,ixe,izs,ize;
+    Float_t x,z;
+    Double_t x1,x2,z1,z2,s,sp;
+    AliITSsegmentationSPD* seg = (AliITSsegmentationSPD*)GetSegmentationModel(0);
+
+
+    if(GetDebug(4)) Info("SpreadCharge","(x0=%e,z0=%e,ix0=%d,iz0=%d,el=%e,"
+                         "sig=%e,t=%d,i=%d)",x0,z0,ix0,iz0,el,sig,t,hi);
+    if(sig<=0.0) { // if sig<=0 No diffusion to simulate.
+        GetMap()->AddSignal(iz0,ix0,t,hi,GetModuleNumber(),el);
+        if(GetDebug(2)){
+            cout << "sig<=0.0=" << sig << endl;
+        } // end if GetDebug
+        return;
+    } // end if
+    sp = 1.0/(sig*kRoot2);
+    if(GetDebug(2)){
+        cout << "sig=" << sig << " sp=" << sp << endl;
+    } // end if GetDebug
+    ixs = TMath::Max(-knx+ix0,0);
+    ixe = TMath::Min(knx+ix0,seg->Npx()-1);
+    izs = TMath::Max(-knz+iz0,0);
+    ize = TMath::Min(knz+iz0,seg->Npz()-1);
+    for(ix=ixs;ix<=ixe;ix++) for(iz=izs;iz<=ize;iz++){
+        seg->DetToLocal(ix,iz,x,z); // pixel center
+        x1  = x;
+        z1  = z;
+        x2  = x1 + 0.5*kmictocm*seg->Dpx(ix); // Upper
+        x1 -= 0.5*kmictocm*seg->Dpx(ix);  // Lower
+        z2  = z1 + 0.5*kmictocm*seg->Dpz(iz); // Upper
+        z1 -= 0.5*kmictocm*seg->Dpz(iz);  // Lower
+        x1 -= x0; // Distance from where track traveled
+        x2 -= x0; // Distance from where track traveled
+        z1 -= z0; // Distance from where track traveled
+        z2 -= z0; // Distance from where track traveled
+        s   = 0.25; // Correction based on definision of Erfc
+        s  *= TMath::Erfc(sp*x1) - TMath::Erfc(sp*x2);
+        if(GetDebug(3)){
+            cout <<"el="<<el<<" ix0="<<ix0<<" ix="<<ix<<" x0="<<x<<
+                " iz0="<<iz0<<" iz="<<iz<<" z0="<<z<< 
+                " sp*x1="<<sp*x1<<" sp*x2="<<sp*x2<<" s="<<s;
+        } // end if GetDebug
+        s  *= TMath::Erfc(sp*z1) - TMath::Erfc(sp*z2);
+        if(GetDebug(3)){
+            cout<<" sp*z1="<<sp*z1<<" sp*z2="<<sp*z2<<" s="<<s<< endl;
+        } // end if GetDebug
+        GetMap()->AddSignal(iz,ix,t,hi,GetModuleNumber(),s*el);
+    } // end for ix, iz
 }
 //______________________________________________________________________
-void AliITSsimulationSPDdubna::ChargeToSignal(AliITSpList *pList){
+void AliITSsimulationSPDdubna::pListToDigits(){
     // add noise and electronics, perform the zero suppression and add the
     // digit to the list
-
-    AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");  
-
-    Float_t threshold = (float)fResponse->MinVal();
-
-    Int_t    digits[3], tracks[3], hits[3], gi, j1;
-    Float_t  charges[3];
-    Float_t  electronics;
-    Float_t  signal;
-    Float_t  phys; 
+    // Inputs:
+    //    none.
+    // Outputs:
+    //    none.
+    // Return:
+    //    none.
+    static AliITS *aliITS = (AliITS*)gAlice->GetModule("ITS");
+    Int_t j,ix,iz;
+    Double_t  electronics;
     Double_t sig;
-//    Int_t    module = 0;
-    for(Int_t iz=0; iz<fNPixelsZ; iz++){
-       for(Int_t ix=0; ix<fNPixelsX; ix++){
-           electronics = fBaseline + fNoise*gRandom->Gaus();
-           signal = (float)pList->GetSignalOnly(iz,ix);
-           sig = Double_t (signal);  // sig will be passed along to 
-                                     // UpdateMapNoise this is necessary so 
-                                     // that a signal without electronic
-                                     // noise is passed along
-           signal += electronics;
-           gi =iz*fNPixelsX+ix; // global index
-           if (signal > threshold) {
-               digits[0]=iz;
-               digits[1]=ix;
-               digits[2]=1;
-               for(j1=0;j1<3;j1++){
-                   if (pList->GetTrack(iz,ix,gi)) {
-                       //b.b.       tracks[j1]=-3;
-                       tracks[j1] = (Int_t)(pList->GetTrack(iz,ix,j1)+j1);
-                       hits[j1] = (Int_t)(pList->GetHit(iz,ix,j1)+j1+6);
-                   }else {
-                       tracks[j1]=-2; //noise
-                       hits[j1] = -1;
-                   } // end if pList
-                   charges[j1] = 0;
-               } // end for j1
-
-               if(tracks[0] == tracks[1] && tracks[0] == tracks[2]) {
-                   tracks[1] = -3;
-                   hits[1] = -1;
-                   tracks[2] = -3;
-                   hits[2] = -1;
-               } else if(tracks[0] == tracks[1] && tracks[0] != tracks[2]) {
-                   tracks[1] = -3;
-                   hits[1] = -1;   
-               } else if(tracks[0] == tracks[2] && tracks[0] != tracks[1]) {
-                   tracks[2] = -3;
-                   hits[2] = -1;   
-               } else if(tracks[1] == tracks[2] && tracks[0] != tracks[1]) {
-                   tracks[2] = -3;
-                   hits[2] = -1;   
-               } // end if
-
-               phys = 0;
-
-               UpdateMapNoise(iz,ix,fModule,sig,electronics,pList);
-               aliITS->AddSimDigit(0, phys, digits, tracks, hits, charges);
-           } // 
-       } // 
-    } //
+    const Int_t    nmaxtrk=AliITSdigitSPD::GetNTracks();
+    static AliITSdigitSPD dig;
+    AliITSCalibrationSPD* res = (AliITSCalibrationSPD*)GetCalibrationModel(fDetType->GetITSgeom()->GetStartSPD());
+    if(GetDebug(1)) Info("pListToDigits","()");
+    for(iz=0; iz<GetNPixelsZ(); iz++) for(ix=0; ix<GetNPixelsX(); ix++){
+        // Apply Noise/Dead channals and the like
+        if(res->IsPixelDead(GetModuleNumber(),ix,iz)) continue;
+        electronics = res->ApplyBaselineAndNoise();
+        UpdateMapNoise(ix,iz,electronics);
+        //
+        // Apply Threshold and write Digits.
+        sig = GetMap()->GetSignalOnly(iz,ix);
+        FillHistograms(ix,iz,sig+electronics);
+        if(GetDebug(3)){
+            cout<<sig<<"+"<<electronics<<">threshold("<<ix<<","<<iz
+                <<")="<<GetThreshold() <<endl;
+        } // end if GetDebug
+        if (sig+electronics <= GetThreshold()) continue;
+        dig.SetCoord1(iz);
+        dig.SetCoord2(ix);
+        dig.SetSignal(1);
+        dig.SetSignalSPD((Int_t) GetMap()->GetSignal(iz,ix));
+        for(j=0;j<nmaxtrk;j++){
+            if (j<GetMap()->GetNEntries()) {
+                dig.SetTrack(j,GetMap()->GetTrack(iz,ix,j));
+                dig.SetHit(j,GetMap()->GetHit(iz,ix,j));
+            }else { // Default values
+                dig.SetTrack(j,-3);
+                dig.SetHit(j,-1);
+            } // end if GetMap()
+        } // end for j
+        if(GetDebug(3)){
+            cout<<iz<<","<<ix<<","<<*(GetMap()->GetpListItem(iz,ix))<<endl;
+        } // end if GetDebug
+        aliITS->AddSimDigit(0,&dig);
+    } //  for ix/iz
 }
 //______________________________________________________________________
 void AliITSsimulationSPDdubna::CreateHistograms(){
     // create 1D histograms for tests
-
-    printf("SPD - create histograms\n");
-
-    fHis=new TObjArray(fNPixelsZ);
-    TString spdName("spd_");
-    for (Int_t i=0;i<fNPixelsZ;i++) {
-       Char_t pixelz[4];
-       sprintf(pixelz,"%d",i+1);
-       spdName.Append(pixelz);
-       //PH       (*fHis)[i] = new TH1F(spdName.Data(),"SPD maps",
-       //PH                       fNPixelsX,0.,(Float_t) fNPixelsX);
-       fHis->AddAt(new TH1F(spdName.Data(),"SPD maps",
-                            fNPixelsX,0.,(Float_t) fNPixelsX), i);
+    // Inputs:
+    //    none.
+    // Outputs:
+    //    none.
+    // Return:
+    //     none.
+
+    if(GetDebug(1)) Info("CreateHistograms","create histograms");
+
+    fHis = new TObjArray(GetNPixelsZ());
+    TString fSPDname("spd_");
+    for(Int_t i=0;i<GetNPixelsZ();i++) {
+        Char_t pixelz[4];
+        sprintf(pixelz,"%d",i);
+        fSPDname.Append(pixelz);
+        fHis->AddAt(new TH1F(fSPDname.Data(),"SPD maps",
+                             GetNPixelsX(),0.,(Double_t)GetNPixelsX()),i);
     } // end for i
 }
 //______________________________________________________________________
+void AliITSsimulationSPDdubna::FillHistograms(Int_t ix,Int_t iz,Double_t v){
+    // Fill the histogram
+    // Inputs:
+    //    none.
+    // Outputs:
+    //    none.
+    // Return:
+    //     none.
+
+    if(!GetHistArray()) return; // Only fill if setup.
+    if(GetDebug(2)) Info("FillHistograms","fill histograms");
+    GetHistogram(iz)->Fill(ix,v);
+}
+//______________________________________________________________________
 void AliITSsimulationSPDdubna::ResetHistograms(){
-    //
     // Reset histograms for this detector
-    //
-
-    for ( int i=0;i<fNPixelsZ;i++ ) {
-       //PH    if ((*fHis)[i])    ((TH1F*)(*fHis)[i])->Reset();
-       if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
+    // Inputs:
+    //    none.
+    // Outputs:
+    //    none.
+    // Return:
+    //     none.
+
+    if(!GetHistArray()) return; // Only fill if setup.
+    if(GetDebug(2)) Info("FillHistograms","fill histograms");
+    for ( int i=0;i<GetNPixelsZ();i++ ) {
+        if (fHis->At(i))    ((TH1F*)fHis->At(i))->Reset();
     } // end for i
 }
+
+//______________________________________________________________________
+void AliITSsimulationSPDdubna::SetCoupling(Int_t row, Int_t col, Int_t ntrack,
+                                     Int_t idhit) {
+    //  Take into account the coupling between adiacent pixels.
+    //  The parameters probcol and probrow are the probability of the
+    //  signal in one pixel shared in the two adjacent pixels along
+    //  the column and row direction, respectively.
+    //  Note pList is goten via GetMap() and module is not need any more.
+    //  Otherwise it is identical to that coded by Tiziano Virgili (BSN).
+    //Begin_Html
+    /*
+      <img src="picts/ITS/barimodel_3.gif">
+      </pre>
+      <br clear=left>
+      <font size=+2 color=red>
+      <a href="mailto:tiziano.virgili@cern.ch"></a>.
+      </font>
+      <pre>
+    */
+    //End_Html
+    // Inputs:
+    //    Int_t row            z cell index
+    //    Int_t col            x cell index
+    //    Int_t ntrack         track incex number
+    //    Int_t idhit          hit index number
+    // Outputs:
+    //    none.
+    // Return:
+    //     none.
+    Int_t j1,j2,flag=0;
+    Double_t pulse1,pulse2;
+    Double_t couplR=0.0,couplC=0.0;
+    Double_t xr=0.;
+
+    GetCouplings(couplR,couplC);
+    if(GetDebug(3)) Info("SetCoupling","(row=%d,col=%d,ntrack=%d,idhit=%d) "
+                         "Calling SetCoupling couplR=%e couplC=%e",
+                         row,col,ntrack,idhit,couplR,couplC);
+    j1 = row;
+    j2 = col;
+    pulse1 = GetMap()->GetSignalOnly(row,col);
+    pulse2 = pulse1;
+    for (Int_t isign=-1;isign<=1;isign+=2){// loop in row direction
+        do{
+            j1 += isign;
+            //   pulse1 *= couplR; 
+            xr = gRandom->Rndm();
+            //if ((j1<0)||(j1>GetNPixelsZ()-1)||(pulse1<GetThreshold())){
+            if ((j1<0) || (j1>GetNPixelsZ()-1) || (xr>couplR)){
+                j1 = row;
+                flag = 1;
+            }else{
+                UpdateMapSignal(col,j1,ntrack,idhit,pulse1);
+                //  flag = 0;
+                flag = 1; // only first next!!
+            } // end if
+        } while(flag == 0);
+        // loop in column direction
+        do{
+            j2 += isign;
+            // pulse2 *= couplC; 
+            xr = gRandom->Rndm();
+            //if((j2<0)||j2>(GetNPixelsX()-1)||pulse2<GetThreshold()){
+            if ((j2<0) || (j2>GetNPixelsX()-1) || (xr>couplC)){
+                j2 = col;
+                flag = 1;
+            }else{
+                UpdateMapSignal(j2,row,ntrack,idhit,pulse2);
+                //  flag = 0;
+                flag = 1; // only first next!!
+            } // end if
+        } while(flag == 0);
+    } // for isign
+}
+//______________________________________________________________________
+void AliITSsimulationSPDdubna::SetCouplingOld(Int_t row, Int_t col,
+                Int_t ntrack,Int_t idhit) {
+    //  Take into account the coupling between adiacent pixels.
+    //  The parameters probcol and probrow are the fractions of the
+    //  signal in one pixel shared in the two adjacent pixels along
+    //  the column and row direction, respectively.
+    //Begin_Html
+    /*
+      <img src="picts/ITS/barimodel_3.gif">
+      </pre>
+      <br clear=left>
+      <font size=+2 color=red>
+      <a href="mailto:Rocco.Caliandro@ba.infn.it"></a>.
+      </font>
+      <pre>
+    */
+    //End_Html
+    // Inputs:
+    //    Int_t row            z cell index
+    //    Int_t col            x cell index
+    //    Int_t ntrack         track incex number
+    //    Int_t idhit          hit index number
+    //    Int_t module         module number
+    // Outputs:
+    //    none.
+    // Return:
+    //     none.
+    Int_t j1,j2,flag=0;
+    Double_t pulse1,pulse2;
+    Double_t couplR=0.0,couplC=0.0;
+
+    GetCouplings(couplR,couplC);
+    if(GetDebug(3)) Info("SetCouplingOld","(row=%d,col=%d,ntrack=%d,idhit=%d) "
+                         "Calling SetCoupling couplR=%e couplC=%e",
+                         row,col,ntrack,idhit,couplR,couplC);
+    j1 = row;
+    j2 = col;
+    pulse1 = GetMap()->GetSignalOnly(row,col);
+    pulse2 = pulse1;
+    for (Int_t isign=-1;isign<=1;isign+=2){// loop in row direction
+        do{
+            j1 += isign;
+            pulse1 *= couplR;
+            if ((j1<0)||(j1>GetNPixelsZ()-1)||(pulse1<GetThreshold())){
+                pulse1 = GetMap()->GetSignalOnly(row,col);
+                j1 = row;
+                flag = 1;
+            }else{
+                UpdateMapSignal(col,j1,ntrack,idhit,pulse1);
+                flag = 0;
+            } // end if
+        } while(flag == 0);
+        // loop in column direction
+        do{
+            j2 += isign;
+            pulse2 *= couplC;
+            if((j2<0)||(j2>(GetNPixelsX()-1))||(pulse2<GetThreshold())){
+                pulse2 = GetMap()->GetSignalOnly(row,col);
+                j2 = col;
+                flag = 1;
+            }else{
+                UpdateMapSignal(j2,row,ntrack,idhit,pulse2);
+                flag = 0;
+            } // end if
+        } while(flag == 0);
+    } // for isign
+}